Skip to content

Commit

Permalink
Add files via upload
Browse files Browse the repository at this point in the history
  • Loading branch information
cody2 authored Jun 27, 2022
1 parent acc0c66 commit c203863
Show file tree
Hide file tree
Showing 2 changed files with 217 additions and 21 deletions.
42 changes: 21 additions & 21 deletions N803_model_2.m
Original file line number Diff line number Diff line change
@@ -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 |
Expand Down Expand Up @@ -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]
Expand All @@ -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 []
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
196 changes: 196 additions & 0 deletions N803_single.m
Original file line number Diff line number Diff line change
@@ -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

0 comments on commit c203863

Please sign in to comment.