From 954e2e8e005ff7a219d3ec71ca8105b22380fb2e Mon Sep 17 00:00:00 2001 From: "Cody, Jonathan William" Date: Sun, 25 Dec 2022 08:39:31 -0600 Subject: [PATCH] Add files via upload --- ICsolvers/N803_shared.m | 226 ++++++++++++++++++++-------------------- ICsolvers/N803_single.m | 190 ++++++++++++++++----------------- 2 files changed, 208 insertions(+), 208 deletions(-) diff --git a/ICsolvers/N803_shared.m b/ICsolvers/N803_shared.m index fae86b7..0510cb1 100644 --- a/ICsolvers/N803_shared.m +++ b/ICsolvers/N803_shared.m @@ -1,7 +1,7 @@ %% N803_shared.m - solves model for 2 cohorts with shared parameters % % /--------------------------------------------------------------\ -% | Date: 08/18/2022 | +% | Date: 12/25/2022 | % | Author: Jonathan Cody | % | Affiliation: Purdue University | % | Weldon School of Biomedical Engineering | @@ -10,10 +10,10 @@ % % 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 [#/無] +% S0 = resting SIV-specific CD8+ T cells [#/無] +% Sa = active SIV-specific CD8+ T cells [#/無] +% N0 = resting non-SIV-specific CD8+ T cells [#/無] +% Na = active non-SIV-specific CD8+ T cells [#/無] % X = N803 at absorption site [pmol/kg] % C = N803 plasma concentration [pM] % R = regulation [] (dimensionless quantity) @@ -75,155 +75,155 @@ % 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) -fEn = AllPars(06) ;% normalized E/(E+B) (co 2) -fE(2) = fE(1)+(0.3-fE(1))*fEn ;% 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] +SNi(1) = AllPars(03) ;% S+N initial value [#/無] (cohort 1) +SNi(2) = AllPars(04) ;% S+N initial value [#/無] (cohort 2) +fS(1) = AllPars(05) ;% initial frequency: S/(S+N) (cohort 1) +fSn = AllPars(06) ;% normalized S/(S+N) (co 2) +fS(2) = fS(1)+(0.3-fS(1))*fSn ;% initial frequency: S/(S+N) (cohort 2) + +q = AllPars(07) ;% V growth rate (if S+N were absent) [/d] +psi = AllPars(08) ;% ratio of base killing rates gN0/gS0 +V50S = AllPars(09) ;% 50% viral stimulation saturation for S [#/mL] +V50N = AllPars(10) ;% 50% viral stimulation saturation for N [#/mL] +mSn = AllPars(11) ;% normalized Sa reversion rate constant [] +mNn = AllPars(12) ;% normalized Na reversion rate constant [] + +SN50 = AllPars(13) ;% 50% S+N proliferation saturation [#/無] +p0n = AllPars(14) ;% nomalized S0/N0 proliferation rate constant [/d] +pS = AllPars(15) ;% Sa proliferation rate constant [/d] +pN = AllPars(16) ;% Na proliferation rate constant [/d] +d = AllPars(17) ;% S0/N0 death rate constant [/d] +dA = AllPars(18) ;% Sa/Na 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 [] +pm = AllPars(24) ;% S0/N0 maximum proliferation rate [] +aS1 = AllPars(25) ;% S activation stimulation factor [] +aN1 = AllPars(26) ;% N activation stimulation factor [] dR = AllPars(27) ;% R decay rate constant [/d] -sig = AllPars(28) ;% ratio of initial regulation due to B/E -p2 = AllPars(29) ;% E0/B0 proliferation regulation factor [] -gB2 = AllPars(30) ;% B killing regulation factor [] (cohort 1) +sig = AllPars(28) ;% ratio of initial regulation due to N/S +p2 = AllPars(29) ;% S0/N0 proliferation regulation factor [] +gN2 = AllPars(30) ;% N killing regulation factor [] (cohort 1) %% ------------------------------------------------------------------------ % 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] +V50S = V50S/1000 ;% 50% viral stimulation saturation for S [#/無] +V50N = V50N/1000 ;% 50% viral stimulation saturation for N [#/無] +YS = Vi ./ (Vi+ V50S) ;% initial V/(V50S+V) [cohort 1,2] +YN = Vi ./ (Vi+ V50N) ;% initial V/(V50N+V) [cohort 1,2] -% calculate initial R, and sE,sB -z = YE + sig*YB ; +% calculate initial R, and sS,sN +z = YS + sig*YN ; Ri = [ 1 , (z(2)/z(1)) ] ;% initial regulation [cohort 1,2] R = Ri(2) ;% initial regulation (cohort 2) -sE = dR / ( YE(1) + sig*YB(1) ) ;% R generation due to E0 activation [/d] -sB = sig*sE ;% R generation due to B0 activation [/d] - -% restrict mE and mB such that initial activation aE and aB are positive -UE = 2*(2*pE/(pE+dA))^7 ; -UB = 2* 2*pB/(pB+dA) ; -mE = mEn*dA/(UE-1) ;% Ea reversion rate constant [/d] -mB = mBn*dA/(UB-1) ;% Ba reversion rate constant [/d] - -% restrict p0 to ensure same sign for E0/EA and for B0/BA -p0_max = min( d*(1+p2*Ri).*(EB50+EBi)/EB50 ) ;% maximum 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*Ri) ;% initial E0/B0 prolif rate [/d] - -% calculate aE0,aB0 and aE2,aB2 -aEi = (d-pi)/( mE*UE/(dA+mE) - 1 ) ;% initial E0 activation rate [/d] -aBi = (d-pi)/( mB*UB/(dA+mB) - 1 ) ;% initial B0 activation rate [/d] -z = aEi./YE ; -aE2 = ( z(1)-z(2) ) / ( R*z(2)-z(1) ) ;% E0 activation rate constant [/d] -aE0 = aEi(1) / YE(1) * (1+aE2) ;% E activation regulation factor [] -z = aBi./YB ; -aB2 = ( z(1)-z(2) ) / ( R*z(2)-z(1) ) ;% B0 activation rate constant [/d] -aB0 = aBi(1) / YB(1) * (1+aB2) ;% B activation regulation factor [] +sS = dR / ( YS(1) + sig*YN(1) ) ;% R generation due to S0 activation [/d] +sN = sig*sS ;% R generation due to N0 activation [/d] + +% restrict mS and mN such that initial activation aS and aN are positive +US = 2*(2*pS/(pS+dA))^7 ; +UN = 2* 2*pN/(pN+dA) ; +mS = mSn*dA/(US-1) ;% Sa reversion rate constant [/d] +mN = mNn*dA/(UN-1) ;% Na reversion rate constant [/d] + +% restrict p0 to ensure same sign for S0/SA and for N0/NA +p0_max = min( d*(1+p2*Ri).*(SN50+SNi)/SN50 ) ;% maximum value for p0 +p0 = p0n * p0_max ;% S0/N0 prolif rate constant [/d] +p1 = pm/p0 ;% S0/N0 proliferation stimulation factor +pi = p0*SN50./(SN50+SNi)./(1+p2*Ri) ;% initial S0/N0 prolif rate [/d] + +% calculate aS0,aN0 and aS2,aN2 +aSi = (d-pi)/( mS*US/(dA+mS) - 1 ) ;% initial S0 activation rate [/d] +aNi = (d-pi)/( mN*UN/(dA+mN) - 1 ) ;% initial N0 activation rate [/d] +z = aSi./YS ; +aS2 = ( z(1)-z(2) ) / ( R*z(2)-z(1) ) ;% S0 activation rate constant [/d] +aS0 = aSi(1) / YS(1) * (1+aS2) ;% S activation regulation factor [] +z = aNi./YN ; +aN2 = ( z(1)-z(2) ) / ( R*z(2)-z(1) ) ;% N0 activation rate constant [/d] +aN0 = aNi(1) / YN(1) * (1+aN2) ;% N activation regulation factor [] % solve for initial ratios below (based on active steady-state) -ZE = UE/(mE+dA) ; +ZS = US/(mS+dA) ; for i = 1:7 - ZE = ZE + 2*(2*pE)^(i-1)/(pE+dA)^i ;% EAi/aEi/E0i + ZS = ZS + 2*(2*pS)^(i-1)/(pS+dA)^i ;% SAi/aSi/S0i end -ZB = UB/(mB+dA) + 2/(pB+dA) ;% BAi/aBi/B0i - -% 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 - -% calculate gE0,gB0 and gE2 -beta = psi * BA ./ (1+gB2*Ri) ; z = beta(1) - beta(2) ; -a = z*R ; b = EA(1)*R - EA(2) + z*(1+R) ; c = EA(1) - EA(2) + z ; -gE2 = ( -b + sqrt(b^2 - 4*a*c) ) / (2*a) ;% E killing regulation factor [] -gE0 = q / ( EA(1)/(1+gE2) + beta(1) ) ;% Ea killing rate constant [無/#-d] -gB0 = psi * gE0 ;% Ba killing rate constant [無/#-d] +ZN = UN/(mN+dA) + 2/(pN+dA) ;% NAi/aNi/N0i + +% solve for initial values of S0,Sa,N0,Na each cohort +Si = SNi.*fS ;% initial S +Ni = SNi.*(1-fS) ;% initial N +S0 = Si./(1+ZS*aSi) ;% initial S0 +N0 = Ni./(1+ZN*aNi) ;% initial N0 +SA = S0.*(ZS*aSi) ;% initial SA +NA = N0.*(ZN*aNi) ;% initial NA + +% calculate gS0,gN0 and gS2 +beta = psi * NA ./ (1+gN2*Ri) ; z = beta(1) - beta(2) ; +a = z*R ; b = SA(1)*R - SA(2) + z*(1+R) ; c = SA(1) - SA(2) + z ; +gS2 = ( -b + sqrt(b^2 - 4*a*c) ) / (2*a) ;% S killing regulation factor [] +gS0 = q / ( SA(1)/(1+gS2) + beta(1) ) ;% Sa killing rate constant [無/#-d] +gN0 = psi * gS0 ;% Na killing rate constant [無/#-d] %% 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(1) = 2*aEi(c)*E0(c) / (dA+pE) ;% E1 +% solve for initial S1-8 and N1-2 +S = zeros(1,8) ;% initial S1-S8 +S(1) = 2*aSi(c)*S0(c) / (dA+pS) ;% S1 for i = 2:7 - E(i) = 2*pE*E(i-1) / (dA+pE) ;% E2 to E7 + S(i) = 2*pS*S(i-1) / (dA+pS) ;% S2 to S7 end -E(8) = 2*pE*E(7) / (dA+mE) ;% E8 +S(8) = 2*pS*S(7) / (dA+mS) ;% S8 -B(1) = 2*aBi(c)*B0(c) / (dA+pB) ;% B1 -B(2) = 2*pB*B(1) / (dA+mB) ;% B2 +N(1) = 2*aNi(c)*N0(c) / (dA+pN) ;% N1 +N(2) = 2*pN*N(1) / (dA+mN) ;% N2 %% ------------------------------------------------------------------------ % 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(01) = q ;% V growth rate (if S+N were absent) [/d] +Pars(02) = gS0 ;% Sa killing rate constant [無/#-d] +Pars(03) = gN0 ;% Na 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(04) = V50S ;% 50% viral stimulation saturation for S [#/無] +Pars(05) = V50N ;% 50% viral stimulation saturation for N [#/無] +Pars(06) = aS0 ;% S0 activation rate constant [/d] +Pars(07) = aN0 ;% N0 activation rate constant [/d] +Pars(08) = mS ;% Sa reversion rate constant [/d] +Pars(09) = mN ;% Na 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(10) = SN50 ;% 50% S+N proliferation saturation [#/無] +Pars(11) = p0 ;% S0/N0 proliferation rate constant [/d] +Pars(12) = pS ;% Sa proliferation rate constant [/d] +Pars(13) = pN ;% Na proliferation rate constant [/d] +Pars(14) = d ;% S0/N0 death rate constant [/d] +Pars(15) = dA ;% Sa/Na 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(21) = p1 ;% S0/N0 proliferation stimulation factor [] +Pars(22) = aS1 ;% S activation stimulation factor [] +Pars(23) = aN1 ;% N activation stimulation factor [] -Pars(24) = sE ;% R generation due to E0 activation [/d] -Pars(25) = sB ;% R generation due to B0 activation [/d] +Pars(24) = sS ;% R generation due to S0 activation [/d] +Pars(25) = sN ;% R generation due to N0 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(c) E0(c) E B0(c) B 0 0 Ri(c) Ri(c) ] ; +Pars(27) = gS2 ;% S killing regulation factor [] +Pars(28) = gN2 ;% N killing regulation factor [] +Pars(29) = p2 ;% S0/N0 proliferation regulation factor [] +Pars(30) = aS2 ;% S activation regulation factor [] +Pars(31) = aN2 ;% N activation regulation factor [] + +% V S0-8 N0-2 X C R initial values +Yic = [ Vi(c) S0(c) S N0(c) N 0 0 Ri(c) Ri(c) ] ; % if any( [ Pars Yic ] < 0 ) % error('Negative parameters or initial values.') diff --git a/ICsolvers/N803_single.m b/ICsolvers/N803_single.m index 117338b..88e5d1b 100644 --- a/ICsolvers/N803_single.m +++ b/ICsolvers/N803_single.m @@ -1,7 +1,7 @@ %% N803_single.m - solves model for 1 cohort % % /--------------------------------------------------------------\ -% | Date: 09/25/2022 | +% | Date: 12/25/2022 | % | Author: Jonathan Cody | % | Affiliation: Purdue University | % | Weldon School of Biomedical Engineering | @@ -10,10 +10,10 @@ % % 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 [#/無] +% S0 = resting SIV-specific CD8+ T cells [#/無] +% Sa = active SIV-specific CD8+ T cells [#/無] +% N0 = resting non-SIV-specific CD8+ T cells [#/無] +% Na = active non-SIV-specific CD8+ T cells [#/無] % X = N803 at absorption site [pmol/kg] % C = N803 plasma concentration [pM] % R = regulation [] (dimensionless quantity) @@ -59,140 +59,140 @@ % 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] +SNi = AllPars(02) ;% S+N initial value [#/無] +fS = AllPars(03) ;% initial frequency: S/(S+N) +fSA = AllPars(04) ;% initial frequency: Sa/S +q = AllPars(05) ;% V growth rate (if S+N were absent) [/d] +psi = AllPars(06) ;% Na/Sa killing rate ratio [gN0/gS0] +V50S = AllPars(07) ;% 50% viral stimulation saturation for S [#/mL] +V50N = AllPars(08) ;% 50% viral stimulation saturation for N [#/mL] +mSn = AllPars(09) ;% normalized Sa reversion rate constant [] +mNn = AllPars(10) ;% normalized Na reversion rate constant [] +SN50 = AllPars(11) ;% 50% S+N proliferation saturation [#/無] +pS = AllPars(12) ;% Sa proliferation rate constant [/d] +pN = AllPars(13) ;% Na proliferation rate constant [/d] +d = AllPars(14) ;% S0/N0 death rate constant [/d] +dA = AllPars(15) ;% Sa/Na 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] -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 +pm = AllPars(21) ;% S0/N0 maximum proliferation rate [] +aS1 = AllPars(22) ;% S activation stimulation factor [] +aN1 = AllPars(23) ;% N activation stimulation factor [] +sig = AllPars(24) ;% sN/sS regulation generation rate ratio 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 [] +gS2 = AllPars(26) ;% initial S killing regulation [] +gN2 = AllPars(27) ;% initial N killing regulation [] +p2 = AllPars(28) ;% initial S0/N0 proliferation regulation [] +aS2 = AllPars(29) ;% initial S activation regulation [] +aN2 = AllPars(30) ;% initial N 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 [#/無] +V50S = V50S/1000 ;% 50% viral stimulation saturation for S [#/無] +V50N = V50N/1000 ;% 50% viral stimulation saturation for N [#/無] -% restrict mE and mB such that initial activation aE and aB are positive -UE = 2*(2*pE/(pE+dA))^7 ; -UB = 2* 2*pB/(pB+dA) ; -mE = mEn*dA/(UE-1) ;% Ea reversion rate constant [/d] -mB = mBn*dA/(UB-1) ;% Ba reversion rate constant [/d] +% restrict mS and mN such that initial activation aS and aN are positive +US = 2*(2*pS/(pS+dA))^7 ; +UN = 2* 2*pN/(pN+dA) ; +mS = mSn*dA/(US-1) ;% Sa reversion rate constant [/d] +mN = mNn*dA/(UN-1) ;% Na reversion rate constant [/d] % solve for initial ratios below (based on active steady-state) -ZE = UE/(mE+dA) ; +ZS = US/(mS+dA) ; for i = 1:7 - ZE = ZE + 2*(2*pE)^(i-1)/(pE+dA)^i ;% EAi/aEi/E0i + ZS = ZS + 2*(2*pS)^(i-1)/(pS+dA)^i ;% SAi/aSi/S0i end -ZB = UB/(mB+dA) + 2/(pB+dA) ;% BAi/aBi/B0i +ZN = UN/(mN+dA) + 2/(pN+dA) ;% NAi/aNi/N0i -WE = 1 ; +WS = 1 ; for i = 1:7 - WE = WE + (mE+dA)*(pE+dA)^(i-1)/(2*pE)^i ;% EAi/E8i + WS = WS + (mS+dA)*(pS+dA)^(i-1)/(2*pS)^i ;% SAi/S8i end -WB = 1 + (mB+dA)/(2*pB) ;% BAi/B2i +WN = 1 + (mN+dA)/(2*pN) ;% NAi/N2i -QE = mE/WE - 1/ZE ;% collection -QB = mB/WB - 1/ZB ;% collection +QS = mS/WS - 1/ZS ;% collection +QN = mN/WN - 1/ZN ;% collection -fBA = 1/( 1 + QB/QE*(1-fEA)/fEA ) ;% initial frequency: Ba/B +fNA = 1/( 1 + QN/QS*(1-fSA)/fSA ) ;% initial frequency: Na/N -% 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 +% solve for S and N initial conditions +Si = SNi*fS ;% initial S +Ni = SNi*(1-fS) ;% initial N +SA = Si*fSA ;% initial Sa +S0 = Si*(1-fSA) ;% initial S0 +NA = Ni*fNA ;% initial Sa +N0 = Ni*(1-fNA) ;% initial S0 -E = zeros(1,8) ;% initial E1-E8 -E(8) = EA/WE ;% E8 -E(7) = E(8) * (mE+dA)/(2*pE) ;% E7 +S = zeros(1,8) ;% initial S1-S8 +S(8) = SA/WS ;% S8 +S(7) = S(8) * (mS+dA)/(2*pS) ;% S7 for i = 6:-1:1 - E(i) = E(i+1) * (pE+dA) / (2*pE) ;% E6 to E1 + S(i) = S(i+1) * (pS+dA) / (2*pS) ;% S6 to S1 end -B = BA/WB * [ (mB+dA) / (2*pB) , 1 ] ;% initial B1-B2 +N = NA/WN * [ (mN+dA) / (2*pN) , 1 ] ;% initial N1-N2 % solve for rate constants -aEi = E(1) / (2*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] +aSi = S(1) / (2*S0) * (pS+dA) ;% initial activation rate for S0 +aNi = aSi * (US*mS/(mS+dA)-1) / (UN*mN/(mN+dA)-1) ;% for N0 +aS0 = aSi * (V50S+Vi)/Vi * (1+aS2) ;% S0 activation rate constant [/d] +aN0 = aNi * (V50N+Vi)/Vi * (1+aN2) ;% N0 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 +pi = d - QS*SA/S0 ;% initial proliferation rate for S0/N0 +p0 = pi * (SN50+SNi)/SN50 * (1+p2) ;% S0/N0 proliferation rate con [/d] +p1 = pm/p0 ;% S0/N0 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] +gS0 = q / ( SA/(1+gS2) + psi*NA/(1+gN2) ) ;% Sa killing rate con [無/#-d] +gN0 = psi*gS0 ;% Na 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 +sS = dR / ( Vi/(V50S+Vi) + sig*Vi/(V50N+Vi) ) ;% regulation gen by S act +sN = sig*sS ;% regulation generation by N 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(01) = q ;% V growth rate (if S+N were absent) [/d] +Pars(02) = gS0 ;% Sa killing rate constant [無/#-d] +Pars(03) = gN0 ;% Na 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(04) = V50S ;% 50% viral stimulation saturation for S [#/無] +Pars(05) = V50N ;% 50% viral stimulation saturation for N [#/無] +Pars(06) = aS0 ;% S0 activation rate constant [/d] +Pars(07) = aN0 ;% N0 activation rate constant [/d] +Pars(08) = mS ;% Sa reversion rate constant [/d] +Pars(09) = mN ;% Na 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(10) = SN50 ;% 50% S+N proliferation saturation [#/無] +Pars(11) = p0 ;% S0/N0 proliferation rate constant [/d] +Pars(12) = pS ;% Sa proliferation rate constant [/d] +Pars(13) = pN ;% Na proliferation rate constant [/d] +Pars(14) = d ;% S0/N0 death rate constant [/d] +Pars(15) = dA ;% Sa/Na 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(21) = p1 ;% S0/N0 proliferation stimulation factor [] +Pars(22) = aS1 ;% S activation stimulation factor [] +Pars(23) = aN1 ;% N activation stimulation factor [] -Pars(24) = sE ;% R generation due to E0 activation [/d] -Pars(25) = sB ;% R generation due to B0 activation [/d] +Pars(24) = sS ;% R generation due to S0 activation [/d] +Pars(25) = sN ;% R generation due to N0 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 ] ; +Pars(27) = gS2 ;% S killing regulation factor [] +Pars(28) = gN2 ;% N killing regulation factor [] +Pars(29) = p2 ;% S0/N0 proliferation regulation factor [] +Pars(30) = aS2 ;% S activation regulation factor [] +Pars(31) = aN2 ;% N activation regulation factor [] + +% V S0-8 N0-2 X C R initial values +Yic = [ Vi S0 S N0 N 0 0 1 1 ] ; if any( [ Pars Yic ] < 0 ) error('Negative parameters or initial values.')