diff --git a/ICsolvers/N803_shared_drug_q.m b/ICsolvers/N803_shared_drug_q.m index a826c37..6cb99c4 100644 --- a/ICsolvers/N803_shared_drug_q.m +++ b/ICsolvers/N803_shared_drug_q.m @@ -1,7 +1,7 @@ %% N803_shared_drug_q.m - solves model for 2 cohorts with shared drug and q % % /--------------------------------------------------------------\ -% | Date: 06/30/2022 | +% | Date: 07/01/2022 | % | Author: Jonathan Cody | % | Affiliation: Purdue University | % | Weldon School of Biomedical Engineering | @@ -54,8 +54,8 @@ % Y_OUT(:,3) = V at points in 'SoluTimes' [log fold change] cohort 2 % Y_OUT(:,4) = T8 at points in 'SoluTimes' [fold change] cohort 2 % -% PARS(1,:) = parameters for cohort 1 (see code) -% PARS(2,:) = parameters for cohort 2 (see code) +% PARS(1,:) = [ parameters , initial conditions ] for cohort 1 (see code) +% PARS(2,:) = [ parameters , initial conditions ] for cohort 2 (see code) % %% ======================================================================== % FUNCTION @@ -223,7 +223,7 @@ % V E0-8 B0-2 X C R initial values Yic = [ Vi E0 E B0 B 0 0 Ri Ri ] ; -if any( [ Yic aE0 aB0 p0 gE0 gB0 sE sB ] < 0 ) +if any( [ Pars Yic ] < 0 ) error('Negative parameters or initial values.') end @@ -240,7 +240,7 @@ Y_HI = ones(sum(idHi),1)*Y_COH(end,:) ;% constant Y for later times Y_COH = [ Y_LO ; Y_COH ; Y_HI ] ;%#ok % total 'solution' matrix Y_Cohort{c} = Y_COH ; - P_Cohort{c} = Pars ; + P_Cohort{c} = [ Pars Yic ] ; end end diff --git a/ICsolvers/N803_shared_reg_fA.m b/ICsolvers/N803_shared_reg_fA.m index f7a0e64..55cf0eb 100644 --- a/ICsolvers/N803_shared_reg_fA.m +++ b/ICsolvers/N803_shared_reg_fA.m @@ -1,7 +1,7 @@ %% N803_shared_reg_fA.m - solves model for 2 cohorts with shared reg and fA % % /--------------------------------------------------------------\ -% | Date: 06/30/2022 | +% | Date: 07/01/2022 | % | Author: Jonathan Cody | % | Affiliation: Purdue University | % | Weldon School of Biomedical Engineering | @@ -54,8 +54,8 @@ % Y_OUT(:,3) = V at points in 'SoluTimes' [log fold change] cohort 2 % Y_OUT(:,4) = T8 at points in 'SoluTimes' [fold change] cohort 2 % -% PARS(1,:) = parameters for cohort 1 (see code) -% PARS(2,:) = parameters for cohort 2 (see code) +% PARS(1,:) = [ parameters , initial conditions ] for cohort 1 (see code) +% PARS(2,:) = [ parameters , initial conditions ] for cohort 2 (see code) % %% ======================================================================== % FUNCTION @@ -233,7 +233,7 @@ % V E0-8 B0-2 X C R initial values Yic = [ Vi E0 E B0 B 0 0 Ri Ri ] ; -if any( [ Yic aE0 aB0 p0 gE0 gB0 sE sB ] < 0 ) +if any( [ Pars Yic ] < 0 ) error('Negative parameters or initial values.') end @@ -250,7 +250,7 @@ Y_HI = ones(sum(idHi),1)*Y_COH(end,:) ;% constant Y for later times Y_COH = [ Y_LO ; Y_COH ; Y_HI ] ;%#ok % total 'solution' matrix Y_Cohort{c} = Y_COH ; - P_Cohort{c} = Pars ; + P_Cohort{c} = [ Pars Yic ] ; end end diff --git a/ICsolvers/N803_single.m b/ICsolvers/N803_single.m index e443bf2..6132764 100644 --- a/ICsolvers/N803_single.m +++ b/ICsolvers/N803_single.m @@ -1,7 +1,7 @@ %% N803_single.m - solves model for 1 cohort % % /--------------------------------------------------------------\ -% | Date: 06/30/2022 | +% | Date: 07/01/2022 | % | Author: Jonathan Cody | % | Affiliation: Purdue University | % | Weldon School of Biomedical Engineering | @@ -49,7 +49,7 @@ % Y_OUT(:,1) = V at points in 'SoluTimes' [log fold change] % Y_OUT(:,2) = T8 at points in 'SoluTimes' [fold change] % -% Pars = parameters for (including calculated) +% Pars = [ parameters , initial conditions ] % %% ======================================================================== % FUNCTION @@ -194,7 +194,7 @@ % V E0-8 B0-2 X C R initial values Yic = [ Vi E0 E B0 B 0 0 1 1 ] ; -if any( [ Yic aE0 aB0 p0 gE0 gB0 sE sB ] < 0 ) +if any( [ Pars Yic ] < 0 ) error('Negative parameters or initial values.') end @@ -212,4 +212,6 @@ Y_HI = ones(sum(idHi),1)*Y_OUT(end,:) ;% constant Y for later times Y_OUT = [ Y_LO ; Y_OUT ; Y_HI ] ;% total 'solution' matrix +Pars = [ Pars Yic ] ;% adding initial conditions to parameter output + end \ No newline at end of file diff --git a/ICsolvers/N803_single_no_B.m b/ICsolvers/N803_single_no_B.m new file mode 100644 index 0000000..bb00ac4 --- /dev/null +++ b/ICsolvers/N803_single_no_B.m @@ -0,0 +1,254 @@ +%% N803_single_no_B.m - solves model for 2 cohorts without bystander T8 +% +% /--------------------------------------------------------------\ +% | Date: 07/01/2022 | +% | Author: Jonathan Cody | +% | Affiliation: Purdue University | +% | Weldon School of Biomedical Engineering | +% | Pienaar Computational Systems Pharmacology Lab | +% \--------------------------------------------------------------/ +% +% Nomenclature: V = SIV virions [#/無] +% T8 = total CD8+ T cells [#/無] +% E0 = resting SIV-specific CD8+ T cells [#/無] +% Ea = active SIV-specific CD8+ T cells [#/無] +% B0 = resting bystander CD8+ T cells [#/無] +% Ba = active bystander CD8+ T cells [#/無] +% X = N803 at absorption site [pmol/kg] +% C = N803 plasma concentration [pM] +% R = regulation [] (dimensionless quantity) +% +%% ======================================================================== +% INPUTS +% ======================================================================== +% +% SoluTimes = ascending vector of days at which to evaluate solution +% +% DoseTimes{c} = ascending vector of days at which to administer doses +% (elements of 'DoseTimes' must also be in 'SoluTimes') +% +% AllPars = vector of parameters (see list in function) +% +%% ======================================================================== +% OPTIONS +% ======================================================================== +% +% SkipTimes{c} = [min max] time point beyond which to skip model solving +% (outputs before 'min' will be made equal to output at 'min') +% (leave as [] to ignore and solve for all 'SoluTimes') +% +% oneCohort = scalar to run model for just one cohort ('1' or '2') +% (leave as [] to ignore and solve for both cohorts) +% +% All additional inputs will be passed as a cell vector to 'N803_model_2' +% and used to define options (see function for list) +% EX: N803_single(SoluTimes,DoseTimes,AllPars,'AbsTol',1e-2} +% will set ode solver absolute tolerance to 1e-2 +% +%% ======================================================================== +% OUTPUTS +% ======================================================================== +% +% Y_OUT(:,1) = V at points in 'SoluTimes' [log fold change] cohort 1 +% Y_OUT(:,2) = T8 at points in 'SoluTimes' [fold change] cohort 1 +% Y_OUT(:,3) = V at points in 'SoluTimes' [log fold change] cohort 2 +% Y_OUT(:,4) = T8 at points in 'SoluTimes' [fold change] cohort 2 +% +% PARS(1,:) = [ parameters , initial conditions ] for cohort 1 (see code) +% PARS(2,:) = [ parameters , initial conditions ] for cohort 2 (see code) +% +%% ======================================================================== +% FUNCTION +% ======================================================================== +function [Y_OUT,PARS] = ... + N803_single_no_B(SoluTimes,DoseTimes,AllPars,SkipTimes,oneCohort,... + varargin) + +if isempty(oneCohort) ; RunCohort = [ 1 1 ] ; +elseif oneCohort == 1 ; RunCohort = [ 1 0 ] ; +else ; RunCohort = [ 0 1 ] ; +end + +Y_Cohort = cell(1,2) ;% cell for storing outputs +P_Cohort = cell(1,2) ;% cell for storing parmeters + +nShared = 4 ;% number of shared parameters +nVaried = 19 ;% number of varied parameters + +% Rename inputed parameters ----------------------------------------------- +Xi = AllPars(01) ;% X initial condition [pmol/kg] +ka = AllPars(02) ;% N803 absorption rate constant [/d] +ke = AllPars(03) ;% N803 elimination rate constant [/d] +vd = AllPars(04) ;% N803 'volume of distribution'/'bioavailability' [L/kg] + +fE = 1 ;% initial frequency: E/(E+B) +psi = 0 ;% Ba/Ea killing rate ratio [gB0/gE0] +V50B = 0 ;% 50% viral stimulation saturation for B [#/mL] +mBn = 5 ;% normalized Ba reversion rate constant [] +pB = 1 ;% Ba proliferation rate constant [/d] +gB2 = 0 ;% initial B killing regulation [] +aB2 = 0 ;% initial B activation regulation [] + +%% Do for each cohort (NOT indenting loop) ================================ +for c = 1:2 + +n = nShared + nVaried*(c-1) ; +% Rename inputed parameters ----------------------------------------------- +Vi = AllPars(01+n) ;% V initial value [log(#/mL)] +EBi = AllPars(02+n) ;% E+B initial value [#/無] +fEA = AllPars(03+n) ;% initial frequency: Ea/E + +q = AllPars(04+n) ;% V growth rate (if E+B were absent) [/d] +V50E = AllPars(05+n) ;% 50% viral stimulation saturation for E [#/mL] +mEn = AllPars(06+n) ;% normalized Ea reversion rate constant [] + +EB50 = AllPars(07+n) ;% 50% E+B proliferation saturation [#/無] +pE = AllPars(08+n) ;% Ea proliferation rate constant [/d] +d = AllPars(09+n) ;% E0/B0 death rate constant [/d] +dA = AllPars(10+n) ;% Ea/Ba death rate constant [/d] + +C50 = AllPars(11+n) ;% 50% N803 stimulation concentration [pM] (Cohort 1) +pm = AllPars(12+n) ;% E0/B0 maximum proliferation rate [] +aE1 = AllPars(13+n) ;% E activation stimulation factor [] +aB1 = AllPars(14+n) ;% B activation stimulation factor [] + +sig = AllPars(15+n) ;% sB/sE regulation generation rate ratio +dR = AllPars(16+n) ;% R decay rate constant [/d] +gE2 = AllPars(17+n) ;% initial E killing regulation [] +p2 = AllPars(18+n) ;% initial E0/B0 proliferation regulation [] +aE2 = AllPars(19+n) ;% initial E activation regulation [] + +%% ------------------------------------------------------------------------ +% Calculate some initial conditions & parameters -------------------------- + +Vi = 10^(Vi-3) ;% V initial value [#/無] +V50E = V50E/1000 ;% 50% viral stimulation saturation for E [#/無] +V50B = V50B/1000 ;% 50% viral stimulation saturation for B [#/無] + +% restrict mE and mB such that initial activation aE and aB are positive +UE = (2*pE/(pE+dA))^7 ; +UB = 2*pB/(pB+dA) ; +mE = mEn*dA/(UE-1) ;% Ea reversion rate constant [/d] +mB = mBn*dA/(UB-1) ;% Ba reversion rate constant [/d] + +% solve for initial ratios below (based on active steady-state) +ZE = UE/(mE+dA) ; +for i = 1:7 + ZE = ZE + (2*pE)^(i-1)/(pE+dA)^i ;% EAi/aEi/E0i +end +ZB = 1/(pB+dA) + UB/(mB+dA) ;% BAi/aBi/B0i + +WE = 1 ; +for i = 1:7 + WE = WE + (mE+dA)*(pE+dA)^(i-1)/(2*pE)^i ;% EAi/E8i +end +WB = 1 + (mB+dA)/(2*pB) ;% BAi/B2i + +QE = mE/WE - 1/ZE ;% collection +QB = mB/WB - 1/ZB ;% collection + +fBA = 1/( 1 + QB/QE*(1-fEA)/fEA ) ;% initial frequency: Ba/B + +% solve for E and B initial conditions +Ei = EBi*fE ;% initial E +Bi = EBi*(1-fE) ;% initial B +EA = Ei*fEA ;% initial Ea +E0 = Ei*(1-fEA) ;% initial E0 +BA = Bi*fBA ;% initial Ea +B0 = Bi*(1-fBA) ;% initial E0 + +E = zeros(1,8) ;% initial E1-E8 +E(8) = EA/WE ;% E8 +E(7) = E(8) * (mE+dA)/(2*pE) ;% E7 +for i = 6:-1:1 + E(i) = E(i+1) * (pE+dA) / (2*pE) ;% E6 to E1 +end +B = BA/WB * [ (mB+dA) / (2*pB) , 1 ] ;% initial B1-B2 + +% solve for rate constants +aEi = E(1) / E0 * (pE+dA) ;% initial activation rate for E0 +aBi = aEi * (UE*mE/(mE+dA)-1) / (UB*mB/(mB+dA)-1) ;% for B0 +aE0 = aEi * (V50E+Vi)/Vi * (1+aE2) ;% E0 activation rate constant [/d] +aB0 = aBi * (V50B+Vi)/Vi * (1+aB2) ;% B0 activation rate constant [/d] + +pi = d - QE*EA/E0 ;% initial proliferation rate for E0/B0 +p0 = pi * (EB50+EBi)/EB50 * (1+p2) ;% E0/B0 proliferation rate con [/d] +p1 = pm/p0 ;% E0/B0 proliferation stimulation factor + +gE0 = q / ( EA/(1+gE2) + psi*BA/(1+gB2) ) ;% Ea killing rate con [無/#-d] +gB0 = psi*gE0 ;% Ba killing rate constant [無/#-d] + +sE = dR / ( Vi/(V50E+Vi) + sig*Vi/(V50B+Vi) ) ;% reg gen by E act +sB = sig*sE ;% regulation generation by B activation +Ri = 1 ;% initial regulation + +%% ------------------------------------------------------------------------ +% Prepare parameter and initial value vectors and call 'N803_model_2' ----- + +Pars(01) = q ;% V growth rate (if E+B were absent) [/d] +Pars(02) = gE0 ;% Ea killing rate constant [無/#-d] +Pars(03) = gB0 ;% Ba killing rate constant [無/#-d] + +Pars(04) = V50E ;% 50% viral stimulation saturation for E [#/無] +Pars(05) = V50B ;% 50% viral stimulation saturation for B [#/無] +Pars(06) = aE0 ;% E0 activation rate constant [/d] +Pars(07) = aB0 ;% B0 activation rate constant [/d] +Pars(08) = mE ;% Ea reversion rate constant [/d] +Pars(09) = mB ;% Ba reversion rate constant [/d] + +Pars(10) = EB50 ;% 50% E+B proliferation saturation [#/無] +Pars(11) = p0 ;% E0/B0 proliferation rate constant [/d] +Pars(12) = pE ;% Ea proliferation rate constant [/d] +Pars(13) = pB ;% Ba proliferation rate constant [/d] +Pars(14) = d ;% E0/B0 death rate constant [/d] +Pars(15) = dA ;% Ea/Ba death rate constant [/d] + +Pars(16) = Xi ;% X initial condition [pmol/kg] +Pars(17) = ka ;% N803 absorption rate constant [/d] +Pars(18) = ke ;% N803 elimination rate constant [/d] +Pars(19) = vd ;% N803 'volume of distribution'/'bioavailability' [L/kg] +Pars(20) = C50 ;% 50% N803 stimulation concentration [pM] +Pars(21) = p1 ;% E0/B0 proliferation stimulation factor [] +Pars(22) = aE1 ;% E activation stimulation factor [] +Pars(23) = aB1 ;% B activation stimulation factor [] + +Pars(24) = sE ;% R generation due to E0 activation [/d] +Pars(25) = sB ;% R generation due to B0 activation [/d] +Pars(26) = dR ;% R decay rate constant [/d] +Pars(27) = gE2 ;% E killing regulation factor [] +Pars(28) = gB2 ;% B killing regulation factor [] +Pars(29) = p2 ;% E0/B0 proliferation regulation factor [] +Pars(30) = aE2 ;% E activation regulation factor [] +Pars(31) = aB2 ;% B activation regulation factor [] + +% V E0-8 B0-2 X C R initial values +Yic = [ Vi E0 E B0 B 0 0 Ri Ri ] ; + +if any( [ Pars Yic ] < 0 ) + error('Negative parameters or initial values.') +end + +% If 'SkipTimes' is empty, do not skip any times +if isempty(SkipTimes{c}) ; SkipTimes{c} = [-inf inf] ; end + +idLo = SoluTimes < SkipTimes{c}(1) ;% index of early times to skip soln +idHi = SoluTimes > SkipTimes{c}(2) ;% index of later times to skip soln +idSol = ~ ( idLo | idHi ) ;% index of times in 'SoluTimes' to solve + +if RunCohort(c) == 1 + Y_COH = N803_model_2(SoluTimes(idSol),DoseTimes{c},Pars,Yic,varargin) ; + Y_LO = ones(sum(idLo),1)*Y_COH(1 ,:) ;% constant Y for early times + Y_HI = ones(sum(idHi),1)*Y_COH(end,:) ;% constant Y for later times + Y_COH = [ Y_LO ; Y_COH ; Y_HI ] ;%#ok % total 'solution' matrix + Y_Cohort{c} = Y_COH ; + P_Cohort{c} = [ Pars Yic ] ; +end + +end + +%% Prepare main outputs 'Y_OUT' and 'PARS' ================================ + +Y_OUT = [ Y_Cohort{1} , Y_Cohort{2} ] ; +PARS = [ P_Cohort{1} ; P_Cohort{2} ] ; + +end \ No newline at end of file diff --git a/ICsolvers/N803_varied_reg_fA.m b/ICsolvers/N803_varied_reg_fA.m new file mode 100644 index 0000000..cc63047 --- /dev/null +++ b/ICsolvers/N803_varied_reg_fA.m @@ -0,0 +1,283 @@ +%% N803_varied_reg_fA.m - solves model for 2 cohorts with shared reg and fA +% +% /--------------------------------------------------------------\ +% | Date: 07/01/2022 | +% | Author: Jonathan Cody | +% | Affiliation: Purdue University | +% | Weldon School of Biomedical Engineering | +% | Pienaar Computational Systems Pharmacology Lab | +% \--------------------------------------------------------------/ +% +% Nomenclature: V = SIV virions [#/無] +% T8 = total CD8+ T cells [#/無] +% E0 = resting SIV-specific CD8+ T cells [#/無] +% Ea = active SIV-specific CD8+ T cells [#/無] +% B0 = resting bystander CD8+ T cells [#/無] +% Ba = active bystander CD8+ T cells [#/無] +% X = N803 at absorption site [pmol/kg] +% C = N803 plasma concentration [pM] +% R = regulation [] (dimensionless quantity) +% +%% ======================================================================== +% INPUTS +% ======================================================================== +% +% SoluTimes = ascending vector of days at which to evaluate solution +% +% DoseTimes{c} = ascending vector of days at which to administer doses +% (elements of 'DoseTimes' must also be in 'SoluTimes') +% +% AllPars = vector of parameters (see list in function) +% +%% ======================================================================== +% OPTIONS +% ======================================================================== +% +% SkipTimes{c} = [min max] time point beyond which to skip model solving +% (outputs before 'min' will be made equal to output at 'min') +% (leave as [] to ignore and solve for all 'SoluTimes') +% +% oneCohort = scalar to run model for just one cohort ('1' or '2') +% (leave as [] to ignore and solve for both cohorts) +% +% All additional inputs will be passed as a cell vector to 'N803_model_2' +% and used to define options (see function for list) +% EX: N803_single(SoluTimes,DoseTimes,AllPars,'AbsTol',1e-2} +% will set ode solver absolute tolerance to 1e-2 +% +%% ======================================================================== +% OUTPUTS +% ======================================================================== +% +% Y_OUT(:,1) = V at points in 'SoluTimes' [log fold change] cohort 1 +% Y_OUT(:,2) = T8 at points in 'SoluTimes' [fold change] cohort 1 +% Y_OUT(:,3) = V at points in 'SoluTimes' [log fold change] cohort 2 +% Y_OUT(:,4) = T8 at points in 'SoluTimes' [fold change] cohort 2 +% +% PARS(1,:) = [ parameters , initial conditions ] for cohort 1 (see code) +% PARS(2,:) = [ parameters , initial conditions ] for cohort 2 (see code) +% +%% ======================================================================== +% FUNCTION +% ======================================================================== +function [Y_OUT,PARS] = ... + N803_varied_reg_fA(SoluTimes,DoseTimes,AllPars,SkipTimes,oneCohort,... + varargin) + +if isempty(oneCohort) ; RunCohort = [ 1 1 ] ; +elseif oneCohort == 1 ; RunCohort = [ 1 0 ] ; +else ; RunCohort = [ 0 1 ] ; +end + +Y_Cohort = cell(1,2) ;% cell for storing outputs +P_Cohort = cell(1,2) ;% cell for storing parmeters + +% Rename inputed parameters ----------------------------------------------- +Vi(1) = AllPars(01) ;% V initial value [log(#/mL)] (cohort 1) +Vi(2) = AllPars(02) ;% V initial value [log(#/mL)] (cohort 2) +EBi(1) = AllPars(03) ;% E+B initial value [#/無] (cohort 1) +EBi(2) = AllPars(04) ;% E+B initial value [#/無] (cohort 2) +fE(1) = AllPars(05) ;% initial frequency: E/(E+B) (cohort 1) +fE(2) = AllPars(06) ;% initial frequency: E/(E+B) (cohort 2) + +q = AllPars(07) ;% V growth rate (if E+B were absent) [/d] +psi = AllPars(08) ;% ratio of base killing rates gB0/gE0 +V50E = AllPars(09) ;% 50% viral stimulation saturation for E [#/mL] +V50B = AllPars(10) ;% 50% viral stimulation saturation for B [#/mL] +mEn = AllPars(11) ;% normalized Ea reversion rate constant [] +mBn = AllPars(12) ;% normalized Ba reversion rate constant [] + +EB50 = AllPars(13) ;% 50% E+B proliferation saturation [#/無] +p0n = AllPars(14) ;% nomalized E0/B0 proliferation rate constant [/d] +pE = AllPars(15) ;% Ea proliferation rate constant [/d] +pB = AllPars(16) ;% Ba proliferation rate constant [/d] +d = AllPars(17) ;% E0/B0 death rate constant [/d] +dA = AllPars(18) ;% Ea/Ba death rate constant [/d] + +Xi = AllPars(19) ;% X initial condition [pmol/kg] +ka = AllPars(20) ;% N803 absorption rate constant [/d] +ke = AllPars(21) ;% N803 elimination rate constant [/d] +vd = AllPars(22) ;% N803 'volume of distribution'/'bioavailability' [L/kg] +C50 = AllPars(23) ;% 50% N803 stimulation concentration [pM] (Cohort 1) +pm = AllPars(24) ;% E0/B0 maximum proliferation rate [] +aE1 = AllPars(25) ;% E activation stimulation factor [] +aB1 = AllPars(26) ;% B activation stimulation factor [] + +dR(1) = AllPars(27) ;% R decay rate constant [/d] (cohort 1) +dR(2) = AllPars(28) ;% R decay rate constant [/d] (cohort 2) +sig(1) = AllPars(29) ;% ratio of initial regulation due to B/E (cohort 1) +sig(2) = AllPars(30) ;% ratio of initial regulation due to B/E (cohort 2) +p2(1) = AllPars(31) ;% E0/B0 proliferation regulation factor [] (cohort 1) +p2(2) = AllPars(32) ;% E0/B0 proliferation regulation factor [] (cohort 2) +aE2s = AllPars(33) ;% E activation regulation factor [] (sampled) +aB2s = AllPars(34) ;% B activation regulation factor [] (sampled) +gE2s = AllPars(35) ;% E killing regulation factor [] (sampled) +gB2(1) = AllPars(36) ;% B killing regulation factor [] (cohort 1) +gB2(2) = AllPars(37) ;% B killing regulation factor [] (cohort 2) + +%% ------------------------------------------------------------------------ +% Calculate some initial conditions & parameters -------------------------- + +Vi = 10.^(Vi - 3) ;% converting V initial value [#/無] +V50E = V50E/1000 ;% 50% viral stimulation saturation for E [#/無] +V50B = V50B/1000 ;% 50% viral stimulation saturation for B [#/無] +YE = Vi ./ (Vi+ V50E) ;% initial V/(V50E+V) [cohort 1,2] +YB = Vi ./ (Vi+ V50B) ;% initial V/(V50B+V) [cohort 1,2] + +% restrict mE and mB such that initial activation aE and aB are positive +UE = (2*pE/(pE+dA))^7 ; +UB = 2*pB/(pB+dA) ; +mE = mEn*dA/(UE-1) ;% Ea reversion rate constant [/d] +mB = mBn*dA/(UB-1) ;% Ba reversion rate constant [/d] + +% solve for initial ratios below (based on active steady-state) +ZE = UE/(mE+dA) ; +for i = 1:7 + ZE = ZE + (2*pE)^(i-1)/(pE+dA)^i ;% EAi/aEi/E0i +end +ZB = 1/(pB+dA) + UB/(mB+dA) ;% BAi/aBi/B0i + +WE = 1 ; +for i = 1:7 + WE = WE + (mE+dA)*(pE+dA)^(i-1)/(2*pE)^i ;% EAi/E8i +end +WB = 1 + (mB+dA)/(2*pB) ;% BAi/B2i + +QE = mE/WE - 1/ZE ;% collection +QB = mB/WB - 1/ZB ;% collection + +% restrict p0 to ensure same sign for E0/EA and for B0/BA +p0_max = min( d*(1+p2).*(EB50+EBi)/EB50 ) ;% maximum possible value for p0 +p0 = p0n * p0_max ;% E0/B0 prolif rate constant [/d] +p1 = pm/p0 ;% E0/B0 proliferation stimulation factor +pi = p0*EB50./(EB50+EBi)./(1+p2) ;% initial E0/B0 proliferation rate [/d] + +% calculate aE0,aB0 (same for each cohort) and aE2,aB2 (different) +aEi = (d-pi)/QE/ZE ;% initial E0 activation rate [/d] +aBi = (d-pi)/QB/ZB ;% initial E0 activation rate [/d] + +aE_YE_ratio = aEi ./ YE ; +if aE_YE_ratio(1) > aE_YE_ratio(2) ; s = 1 ; c = 2 ; +else ; s = 2 ; c = 1 ; +end +aE2(s) = aE2s ;% regulation of E activation (sampled) +aE2(c) = aE_YE_ratio(s) / aE_YE_ratio(c) * (1+aE2s) - 1 ;% (calculated) + +aB_YB_ratio = aBi ./ YB ; +if aB_YB_ratio(1) > aB_YB_ratio(2) ; s = 1 ; c = 2 ; +else ; s = 2 ; c = 1 ; +end +aB2(s) = aB2s ;% regulation of B activation (sampled) +aB2(c) = aB_YB_ratio(s) / aB_YB_ratio(c) * (1+aB2s) - 1 ;% (calculated) + +aE0 = aE_YE_ratio .* (1+aE2) ;% base activation constant for E +aB0 = aB_YB_ratio .* (1+aB2) ;% base activation constant for B + +% solve for initial values of E0,Ea,B0,Ba each cohort +Ei = EBi.*fE ;% initial E +Bi = EBi.*(1-fE) ;% initial B +E0 = Ei./(1+ZE*aEi) ;% initial E0 +B0 = Bi./(1+ZB*aBi) ;% initial B0 +EA = E0.*(ZE*aEi) ;% initial EA +BA = B0.*(ZB*aBi) ;% initial BA +Ri = 1 ;% regulation normalized to initial value + +% calculate gE0,gB0 (same for each cohort) and gE2 (different) +if EA(1) < EA(2) ; s = 1 ; c = 2 ; +else ; s = 2 ; c = 1 ; +end +gE2(s) = gE2s ;% regulation of E killing (sampled) +beta = psi .* BA ./ (1+gB2) ; +eps(s) = EA(s) / (1+gE2s) ; +eps(c) = eps(s) + beta(s) - beta(c) ; +gE2(c) = EA(c) / eps(c) - 1 ;% regulation of E killing (calculated) + +gE0 = q ./ ( eps + beta ) ;% base killing constant for E +gB0 = psi * gE0 ;% base killing constant for B + +% calculate sE,sB (different for each cohort) +sE = dR ./ ( YE + sig.*YB ) ;% reg gen by E act +sB = sig.*sE ;% regulation generation by B activation + +%% Do for each cohort (NOT indenting loop) ================================ +for c = 1:2 + +% solve for initial E1-8 and B1-2 +E = zeros(1,8) ;% initial E1-E8 +E(8) = EA(c)/WE ;% E8 +E(7) = E(8) * (mE+dA)/(2*pE) ;% E7 +for i = 6:-1:1 + E(i) = E(i+1) * (pE+dA) / (2*pE) ;% E6 to E1 +end +B = BA(c)/WB * [ (mB+dA) / (2*pB) , 1 ] ;% initial B1-B2 + +%% ------------------------------------------------------------------------ +% Prepare parameter and initial value vectors and call 'N803_model_2' ----- + +Pars(01) = q ;% V growth rate (if E+B were absent) [/d] +Pars(02) = gE0(1) ;% Ea killing rate constant [無/#-d] +Pars(03) = gB0(1) ;% Ba killing rate constant [無/#-d] + +Pars(04) = V50E ;% 50% viral stimulation saturation for E [#/無] +Pars(05) = V50B ;% 50% viral stimulation saturation for B [#/無] +Pars(06) = aE0(1) ;% E0 activation rate constant [/d] +Pars(07) = aB0(1) ;% B0 activation rate constant [/d] +Pars(08) = mE ;% Ea reversion rate constant [/d] +Pars(09) = mB ;% Ba reversion rate constant [/d] + +Pars(10) = EB50 ;% 50% E+B proliferation saturation [#/無] +Pars(11) = p0 ;% E0/B0 proliferation rate constant [/d] +Pars(12) = pE ;% Ea proliferation rate constant [/d] +Pars(13) = pB ;% Ba proliferation rate constant [/d] +Pars(14) = d ;% E0/B0 death rate constant [/d] +Pars(15) = dA ;% Ea/Ba death rate constant [/d] + +Pars(16) = Xi ;% X initial condition [pmol/kg] +Pars(17) = ka ;% N803 absorption rate constant [/d] +Pars(18) = ke ;% N803 elimination rate constant [/d] +Pars(19) = vd ;% N803 'volume of distribution'/'bioavailability' [L/kg] +Pars(20) = C50 ;% 50% N803 stimulation concentration [pM] +Pars(21) = p1 ;% E0/B0 proliferation stimulation factor [] +Pars(22) = aE1 ;% E activation stimulation factor [] +Pars(23) = aB1 ;% B activation stimulation factor [] + +Pars(24) = sE(c) ;% R generation due to E0 activation [/d] +Pars(25) = sB(c) ;% R generation due to B0 activation [/d] +Pars(26) = dR(c) ;% R decay rate constant [/d] +Pars(27) = gE2(c) ;% E killing regulation factor [] +Pars(28) = gB2(c) ;% B killing regulation factor [] +Pars(29) = p2(c) ;% E0/B0 proliferation regulation factor [] +Pars(30) = aE2(c) ;% E activation regulation factor [] +Pars(31) = aB2(c) ;% B activation regulation factor [] + +% V E0-8 B0-2 X C R initial values +Yic = [ Vi(c) E0(c) E B0(c) B 0 0 Ri Ri ] ; + +if any( [ Pars Yic ] < 0 ) + error('Negative parameters or initial values.') +end + +% If 'SkipTimes' is empty, do not skip any times +if isempty(SkipTimes{c}) ; SkipTimes{c} = [-inf inf] ; end + +idLo = SoluTimes < SkipTimes{c}(1) ;% index of early times to skip soln +idHi = SoluTimes > SkipTimes{c}(2) ;% index of later times to skip soln +idSol = ~ ( idLo | idHi ) ;% index of times in 'SoluTimes' to solve + +if RunCohort(c) == 1 + Y_COH = N803_model_2(SoluTimes(idSol),DoseTimes{c},Pars,Yic,varargin) ; + Y_LO = ones(sum(idLo),1)*Y_COH(1 ,:) ;% constant Y for early times + Y_HI = ones(sum(idHi),1)*Y_COH(end,:) ;% constant Y for later times + Y_COH = [ Y_LO ; Y_COH ; Y_HI ] ;%#ok % total 'solution' matrix + Y_Cohort{c} = Y_COH ; + P_Cohort{c} = [ Pars Yic ] ; +end + +end + +%% Prepare main outputs 'Y_OUT' and 'PARS' ================================ + +Y_OUT = [ Y_Cohort{1} , Y_Cohort{2} ] ; +PARS = [ P_Cohort{1} ; P_Cohort{2} ] ; + +end \ No newline at end of file