diff --git a/ICsolvers/N803_shared_drug_q.m b/ICsolvers/N803_shared_drug_q.m deleted file mode 100644 index 6cb99c4..0000000 --- a/ICsolvers/N803_shared_drug_q.m +++ /dev/null @@ -1,253 +0,0 @@ -%% N803_shared_drug_q.m - solves model for 2 cohorts with shared drug and q -% -% /--------------------------------------------------------------\ -% | 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_shared_drug_q(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 = 09 ;% number of shared parameters -nVaried = 21 ;% 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] - -q = AllPars(05) ;% V growth rate (if E+B were absent) [/d] -C50 = AllPars(06) ;% 50% N803 stimulation concentration [pM] (Cohort 1) -pm = AllPars(07) ;% E0/B0 maximum proliferation rate [] -aE1 = AllPars(08) ;% E activation stimulation factor [] -aB1 = AllPars(09) ;% B activation stimulation factor [] - -%% 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 [#/無] -fE = AllPars(03+n) ;% initial frequency: E/(E+B) -fEA = AllPars(04+n) ;% initial frequency: Ea/E - -psi = AllPars(05+n) ;% Ba/Ea killing rate ratio [gB0/gE0] -V50E = AllPars(06+n) ;% 50% viral stimulation saturation for E [#/mL] -V50B = AllPars(07+n) ;% 50% viral stimulation saturation for B [#/mL] -mEn = AllPars(08+n) ;% normalized Ea reversion rate constant [] -mBn = AllPars(09+n) ;% normalized Ba reversion rate constant [] - -EB50 = AllPars(10+n) ;% 50% E+B proliferation saturation [#/無] -pE = AllPars(11+n) ;% Ea proliferation rate constant [/d] -pB = AllPars(12+n) ;% Ba proliferation rate constant [/d] -d = AllPars(13+n) ;% E0/B0 death rate constant [/d] -dA = AllPars(14+n) ;% Ea/Ba death rate constant [/d] - -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 [] -gB2 = AllPars(18+n) ;% initial B killing regulation [] -p2 = AllPars(19+n) ;% initial E0/B0 proliferation regulation [] -aE2 = AllPars(20+n) ;% initial E activation regulation [] -aB2 = AllPars(21+n) ;% initial B 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