diff --git a/N803_model_2.m b/N803_model_2.m index d88276d..43ccc5d 100644 --- a/N803_model_2.m +++ b/N803_model_2.m @@ -1,7 +1,7 @@ %% N803_model_2.m - solves ODE model for N-803 treatment of SIV % % /--------------------------------------------------------------\ -% | Date: 06/26/2022 | +% | Date: 06/27/2022 | % | Author: Jonathan Cody | % | Affiliation: Purdue University | % | Weldon School of Biomedical Engineering | @@ -90,18 +90,18 @@ % Rename inputed parameters ----------------------------------------------- q = Pars(01) ;% V growth rate (if E+B were absent) [/d] -gE0 = Pars(02) ;% Ea killing rate constant [/d] -gB0 = Pars(03) ;% Ba killing rate constant [/d] +gE0 = Pars(02) ;% Ea killing rate constant [無/#-d] +gB0 = Pars(03) ;% Ba killing rate constant [無/#-d] V50E = Pars(04) ;% 50% viral stimulation saturation for E [#/無] V50B = Pars(05) ;% 50% viral stimulation saturation for B [#/無] -aE0 = Pars(06) ;% E activation rate constant [/d] -aB0 = Pars(07) ;% B activation rate constant [/d] +aE0 = Pars(06) ;% E0 activation rate constant [/d] +aB0 = Pars(07) ;% B0 activation rate constant [/d] mE = Pars(08) ;% Ea reversion rate constant [/d] mB = Pars(09) ;% Ba reversion rate constant [/d] EB50 = Pars(10) ;% 50% E+B proliferation saturation [#/無] -p0 = Pars(11) ;% E0/P0 proliferation rate constant [/d] +p0 = Pars(11) ;% E0/B0 proliferation rate constant [/d] pE = Pars(12) ;% Ea proliferation rate constant [/d] pB = Pars(13) ;% Ba proliferation rate constant [/d] d = Pars(14) ;% E0/B0 death rate constant [/d] @@ -112,12 +112,12 @@ ke = Pars(18) ;% N803 elimination rate constant [/d] vd = Pars(19) ;% N803 'volume of distribution'/'bioavailability' [L/kg] C50 = Pars(20) ;% 50% N803 stimulation concentration [pM] -p1 = Pars(21) ;% E0/B0 proliferation stimulation factor +p1 = Pars(21) ;% E0/B0 proliferation stimulation factor [] aE1 = Pars(22) ;% E activation stimulation factor [] aB1 = Pars(23) ;% B activation stimulation factor [] -sE = Pars(24) ;% R generation due to E activation [/d] -sB = Pars(25) ;% R generation due to B activation [/d] +sE = Pars(24) ;% R generation due to E0 activation [/d] +sB = Pars(25) ;% R generation due to B0 activation [/d] dR = Pars(26) ;% R decay rate constant [/d] gE2 = Pars(27) ;% E killing regulation factor [] gB2 = Pars(28) ;% B killing regulation factor [] @@ -221,8 +221,8 @@ Y_OUT(:,47) = log10( F(:,09) ) ;% R multiplier for gB* (log) Y_OUT(:,48) = gB.*Ba ./ (gE.*Ea + gB.*Ba) ;% B killing fraction -RgenE = sE*aE0* F(:,11).* F(:,03) ;% R gen by aE -RgenB = sB*aB0*( F(:,12) + F(:,04) ) ;% R gen by aB +RgenE = sE* F(:,11).* F(:,03) ;% R gen by aE +RgenB = sB*( F(:,12) + F(:,04) ) ;% R gen by aB Y_OUT(:,49) = log10( RgenE ) ;% R gen by aE (log) Y_OUT(:,50) = log10( RgenB ) ;% R gen by aB (log) Y_OUT(:,51) = RgenB ./ (RgenE + RgenB) ;% R gen by aB fraction @@ -355,8 +355,8 @@ J(15,15) = -ke ; % Partial derivatives of f16 = dR1/dt - J(16,01) = sE*daEdV/F(06) + sB*daBdV/F(07) ; - J(16,15) = sE*daEdC/F(06) + sB*daBdC/F(07) ; + J(16,01) = sE/aE0*daEdV/F(06) + sB/aB0*daBdV/F(07) ; + J(16,15) = sE/aE0*daEdC/F(06) + sB/aB0*daBdC/F(07) ; J(16,16) = -dR ; % Partial derivatives of f17 = dR2/dt @@ -386,16 +386,16 @@ F(:,08) = 1./ (1+gE2*R) ;% Ea killing regulation [] F(:,09) = 1./ (1+gB2*R) ;% Ba killing regulation [] - F(:,10) = EB50./(EB50+EB) ;% E0/B0 proliferation density dependence [] - F(:,11) = V./ (V50E+V) ;% Ea activation antigen dependence [] - F(:,12) = V./ (V50B+V) ;% Ba activation antigen dependence [] + F(:,10) = EB50./(EB50+EB) ;% E0/B0 proliferation density dependence [] + F(:,11) = V./ (V50E+V) ;% Ea activation antigen dependence [] + F(:,12) = V./ (V50B+V) ;% Ba activation antigen dependence [] - F(:,13) = F(:,10).* F(:,02).* F(:,05) ;% E0/B0 prolif modifier [] - F(:,14) = F(:,11).* F(:,03).* F(:,06) ;% Ea activation modifier [] - F(:,15) = (F(:,12)+F(:,04)).* F(:,07) ;% Ba activation modifier [] + F(:,13) = F(:,10).* F(:,02).* F(:,05) ;% E0/B0 prolif modifier [] + F(:,14) = F(:,11).* F(:,03).* F(:,06) ;% Ea activation modifier [] + F(:,15) = (F(:,12)+F(:,04)).* F(:,07) ;% Ba activation modifier [] - F(:,16) = sE*aE0* F(:,11).* F(:,03) ... - + sB*aB0*( F(:,12) + F(:,04) ) ;% regulation gen [/d] + F(:,16) = sE* F(:,11).* F(:,03) ... + + sB*( F(:,12) + F(:,04) ) ;% regulation gen [/d] end end \ No newline at end of file diff --git a/N803_single.m b/N803_single.m new file mode 100644 index 0000000..ad44d84 --- /dev/null +++ b/N803_single.m @@ -0,0 +1,196 @@ +%% N803_shared_drug.m - solves model for 2 cohorts with shared drug pars +% +% /--------------------------------------------------------------\ +% | Date: 06/27/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 = 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 +% ======================================================================== +% +% varargin = cell vector used to add optional inputs (see 'N803_model_2') +% +%% ======================================================================== +% OUTPUTS +% ======================================================================== +% +% 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) +% +%% ======================================================================== +% FUNCTION +% ======================================================================== +function [Y_OUT,Pars] = ... + N803_single(SoluTimes,DoseTimes,AllPars,varargin) + +% Rename inputed parameters ----------------------------------------------- +Vi = AllPars(01) ;% V initial value [log(#/mL)] +EBi = AllPars(02) ;% E+B initial value [#/無] +fE = AllPars(03) ;% initial frequency: E/(E+B) +fEA = AllPars(04) ;% initial frequency: Ea/E +q = AllPars(05) ;% V growth rate (if E+B were absent) [/d] +psi = AllPars(06) ;% Ba/Ea killing rate ratio [gB0/gE0] +V50E = AllPars(07) ;% 50% viral stimulation saturation for E [#/mL] +V50B = AllPars(08) ;% 50% viral stimulation saturation for B [#/mL] +mEn = AllPars(09) ;% normalized Ea reversion rate constant [] +mBn = AllPars(10) ;% normalized Ba reversion rate constant [] +EB50 = AllPars(11) ;% 50% E+B proliferation saturation [#/無] +pE = AllPars(12) ;% Ea proliferation rate constant [/d] +pB = AllPars(13) ;% Ba proliferation rate constant [/d] +d = AllPars(14) ;% E0/B0 death rate constant [/d] +dA = AllPars(15) ;% Ea/Ba death rate constant [/d] +Xi = AllPars(16) ;% X initial condition [pmol/kg] +ka = AllPars(17) ;% N803 absorption rate constant [/d] +ke = AllPars(18) ;% N803 elimination rate constant [/d] +vd = AllPars(19) ;% N803 'volume of distribution'/'bioavailability' [L/kg] +C50 = AllPars(20) ;% 50% N803 stimulation concentration [pM] (Cohort 1) +pm = AllPars(21) ;% E0/B0 maximum proliferation rate [] +aE1 = AllPars(22) ;% E activation stimulation factor [] +aB1 = AllPars(23) ;% B activation stimulation factor [] +sig = AllPars(24) ;% sB/sE regulation generation rate ratio [/d] +dR = AllPars(25) ;% R decay rate constant [/d] +gE2 = AllPars(26) ;% initial E killing regulation [] +gB2 = AllPars(27) ;% initial B killing regulation [] +p2 = AllPars(28) ;% initial E0/B0 proliferation regulation [] +aE2 = AllPars(29) ;% initial E activation regulation [] +aB2 = AllPars(30) ;% 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 = B(1) / B0 * (pB+dA) ;% initial activation rate 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) ) ;% regulation gen by E act +sB = sig*sE ;% regulation generation by B activation + +%% ------------------------------------------------------------------------ +% 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 1 1 ] ; + +if any( [ Yic aE0 aB0 p0 gE0 gB0 sE sB ] < 0 ) + error('Negative parameters or initial values.') +end + +Y_OUT = N803_model_2(SoluTimes,DoseTimes,Pars,Yic,varargin) ; + +end \ No newline at end of file