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 Oct 29, 2021
1 parent d9d2f5a commit 3f763ad
Showing 1 changed file with 35 additions and 32 deletions.
67 changes: 35 additions & 32 deletions N803_model_2.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,19 +2,17 @@
% during a subcutaneously administered regimen of N803

% /--------------------------------------------------------------\
% | Date: 10/26/2021 |
% | Date: 10/27/2021 |
% | Author: Jonathan Cody |
% | Affiliation: Pienaar Computational Systems Pharmacology Lab |
% | Weldon School of Biomedical Engineering |
% | Purdue University |
% \--------------------------------------------------------------/

% Version 2C
% Version 2D
%
% TOODOO
% -> reconsider p1 calculation
% -> remove regulation feedback
% -> add death to all activated cells

% Nomenclature: V = SIV virions [#/無]
% E0 = resting SIV-specific CD8+ T cells [#/無]
Expand Down Expand Up @@ -90,8 +88,8 @@
pB = Pars(11) ;% Ba proliferation rate constant [/d]
d = Pars(12) ;% E0/B0 death rate constant [/d]
dA = Pars(13) ;% Ea/Ba death rate constant [/d]
mE = Pars(14) ;% Ea reversion to memory rate constant [/d]
mB = Pars(15) ;% Ba reversion to memory rate constant [/d]
mEn = Pars(14) ;% normalized Ea reversion rate constant []
mBn = Pars(15) ;% normalized Ba reversion rate constant []
Xi = Pars(16) ;% X initial condition [pmol/kg]
ka = Pars(17) ;% N803 absorption rate constant [/d]
ke = Pars(18) ;% N803 elimination rate constant [/d]
Expand Down Expand Up @@ -120,24 +118,29 @@
V50E = V50E/1000 ;% 50% viral stimulation saturation for E [#/無]
V50B = V50B/1000 ;% 50% viral stimulation saturation for B [#/無]

% restrict mE and mB such that phi (see below) is positive
kE = (2*pE/(dA+pE))^7 ;
kB = 2*pB/(dA+pB) ;
mE = mEn*dA/(kE-1) ;% Ea reversion rate constant [/d]
mB = mBn*dA/(kB-1) ;% Ba reversion rate constant [/d]

% ratio between initial activation rate of E and of B
phi = (2*mB/(mB+dA) - 1) / (2^7*mE/(mE+dA) - 1) ;
phi = (kB*mB/(mB+dA) - 1) / (kE*mE/(mE+dA) - 1) ;

% ratio between initial activation rate of E and initial regulation
ome = dR / (sE + sB/phi) ;
% solve for aE,aB initial values and E and B initial values
x = fzero(@(x)icFun(x),[-8 2]) ;
[~,aEi,aBi,Ei,Bi] = icFun(x) ;

% solve for aE,aB initial values and initial values (HVL cohort)
% options = optimset('TolX',1e-30) ;
x = fzero(@(x)icFun(x),[-6 6]) ;%,options) ;
[~,aEi,aBi,Ei,Bi,Ri] = icFun(x) ;
% regulation initial value
Ri = ( sE*aEi*(1+aE2R) + sB*aBi*(1+aB2R) ) / dR ;

% Ea/Ba killing rate constants gE0/gB0
gE0 = q / ( sum(Ei(2:9))/(1+gE2R) + psi*sum(Bi(2:3))/(1+gB2R) ) ;
gB0 = psi*gE0 ;
p0 = (d+aEi-mE*Ei(9)/Ei(1)) * (1+p2R) * (EB50+EBi)/EB50 ;% E0/B0 prolif
aE0 = aEi * (V50E+Vi)/Vi * (1+aE2R) ;% E activation rate constant [/d]
aB0 = aBi * (V50B+Vi)/Vi * (1+aB2R) ;% B activation rate constant [/d]
p1 = pm/p0 - 1 ;% E0/B0 proliferation stimulation factor
p1 = pm/p0 ;% E0/B0 proliferation stimulation factor
p2 = p2R/Ri ;% E0/B0 proliferation regulation factor
aE2 = aE2R/Ri ;% E activation regulation factor []
aB2 = aB2R/Ri ;% B activation regulation factor []
Expand Down Expand Up @@ -229,7 +232,7 @@
Y_OUT(:,34) = log10( F(:,18)/F(1,18) ) ;% lfc in regulation generation

% Ba regulation stimulation proportion
Y_OUT(:,35) = sB*aB0*F(:,17) ./ F(:,18) ;
Y_OUT(:,35) = sB*aB0*( F(:,14) + F(:,06) ) ./ F(:,18) ;

% lfc in killing
Kill = gB0*Ba./(1+gB2*R) + gE0*Ea./(1+gE2*R) ;
Expand All @@ -241,18 +244,17 @@
%% ========================================================================
% NESTED FUNCTIONS
% ========================================================================
function [targ,aEi,aBi,Ei,Bi,Ri] = icFun(x)
function [targ,aEi,aBi,Ei,Bi] = icFun(x)

% declare initial activation rates aE and aB
Ri = 10^x ;
aEi = ome*Ri ;
aEi = 10^x ;
aBi = aEi/phi ;

% solve for E0-E8 and B0-B2 initial conditions
MAT = zeros(12) ;
MAT( 1:13: 92) = [ aEi,2*pE*ones(1,7) ] ;
MAT( 13:13:130) = [ -pE*ones(1,7),-(dA+mE),aBi,2*pB ] ;
MAT(129:13:142) = [ -pB,-(dA+mB) ] ;
MAT( 13:13:130) = [ -(dA+pE)*ones(1,7),-(dA+mE),aBi,2*pB ] ;
MAT(129:13:142) = [ -(dA+pB),-(dA+mB) ] ;
MAT(11:12,:) = [ ones(1,12) ; [ ones(1,9) 0 0 0 ] ] ;
Sol = MAT \ [ zeros(10,1) ; EBi ; fE*EBi ] ;
Ei = Sol(1:9) ;
Expand Down Expand Up @@ -288,17 +290,17 @@
dY = Y ;% dY/dt

% calculate rates of change for V,E,B,X,C,R
dY(01) = q*V - gE*sum(Ea)*V - gB*sum(Ba)*V ;% dV/dt
dY(02) = p*E0 - aE*E0 - d*E0 + mE*Ea(8) ;% dE0/dt
dY(03) = + aE*E0 - pE*Ea(1) ;% dE1/dt
dY(4:9)= 2*pE*Ea(1:6) - pE*Ea(2:7) ;% dE2-7/dt
dY(10) = 2*pE*Ea(7) - dA*Ea(8) - mE*Ea(8) ;% dE8/dt
dY(11) = p*B0 - aB*B0 - d*B0 + mB*Ba(2) ;% dB0/dt
dY(12) = + aB*B0 - pB*Ba(1) ;% dB1/dt
dY(13) = 2*pB*Ba(1) - dA*Ba(2) - mB*Ba(2) ;% dB2/dt
dY(14) = - ka*X ;% dX/dt
dY(15) = ka*X/vd - ke*C ;% dC/dt
dY(16) = F(18) - dR*R ;% dR/dt
dY(01) = q*V - gE*sum(Ea)*V - gB*sum(Ba)*V ;% dV/dt
dY(02) = p*E0 - aE*E0 - d*E0 + mE*Ea(8) ;% dE0/dt
dY(03) = + aE*E0 - dA*Ea(1) - pE*Ea(1) ;% dE1/dt
dY(4:9)= 2*pE*Ea(1:6) - dA*Ea(2:7) - pE*Ea(2:7) ;% dE2-7/dt
dY(10) = 2*pE*Ea(7) - dA*Ea(8) - mE*Ea(8) ;% dE8/dt
dY(11) = p*B0 - aB*B0 - d*B0 + mB*Ba(2) ;% dB0/dt
dY(12) = + aB*B0 - dA*Ba(2) - pB*Ba(1) ;% dB1/dt
dY(13) = 2*pB*Ba(1) - dA*Ba(2) - mB*Ba(2) ;% dB2/dt
dY(14) = - ka*X ;% dX/dt
dY(15) = ka*X/vd - ke*C ;% dC/dt
dY(16) = F(18) - dR*R ;% dR/dt

% calculate rates of change for tolerance
if N ~= 0
Expand Down Expand Up @@ -342,7 +344,8 @@
F(:,16) = F(:,13).* F(:,05).* F(:,08) ;% Ea activation modifier []
F(:,17) = (F(:,14)+F(:,06)).* F(:,09) ;% Ba activation modifier []

F(:,18) = sE*aE0*F(:,16) + sB*aB0*F(:,17) ;% regulation gen [/d]
F(:,18) = sE*aE0* F(:,13).* F(:,05) ...
+ sB*aB0*( F(:,14) + F(:,06) ) ;% regulation gen [/d]
end

end

0 comments on commit 3f763ad

Please sign in to comment.