Skip to content

Commit

Permalink
Add files via upload
Browse files Browse the repository at this point in the history
  • Loading branch information
@cody2
cody2 authored Dec 25, 2022
1 parent 56d8f63 commit 954e2e8
Show file tree
Hide file tree
Showing 2 changed files with 208 additions and 208 deletions.
226 changes: 113 additions & 113 deletions ICsolvers/N803_shared.m
View file
Original file line number Diff line number Diff line change
@@ -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 |
Expand All @@ -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)
Expand Down Expand Up @@ -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.')
Expand Down
Loading

0 comments on commit 954e2e8

Please sign in to comment.