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 May 15, 2023
1 parent df6636f commit 4de9462
Show file tree
Hide file tree
Showing 21 changed files with 4,928 additions and 600 deletions.
403 changes: 222 additions & 181 deletions N803_collector_2.m
View file

Large diffs are not rendered by default.

643 changes: 308 additions & 335 deletions N803_model_2.m
View file

Large diffs are not rendered by default.

412 changes: 412 additions & 0 deletions N803_model_2A.m
View file

Large diffs are not rendered by default.

347 changes: 347 additions & 0 deletions N803_model_2B.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,347 @@
%% N803_model_2B.m - solves ODE model for N-803 treatment of SIV
%
% /--------------------------------------------------------------\
% | Date: 05/07/2023 |
% | Author: Jonathan Cody |
% | Affiliation: Purdue University |
% | Weldon School of Biomedical Engineering |
% | Pienaar Computational Systems Pharmacology Lab |
% \--------------------------------------------------------------/
%
% Model 2B is a version of model 2 that models all CD8+ T cells as the
% SIV-specific cells from Model 2, which now only divide 4 times after
% activation. Also, regulation is generated with independent parameters.
%
% Nomenclature: V = SIV virions [#/mL]
% E = CD8+ T cells [#/uL]
% X = N803 at absorption site [pmol/kg]
% C = N803 plasma concentration [pM]
% R = immune regulation [] (dimensionless quantity)
%
%% ========================================================================
% INPUTS
% ========================================================================
%
% SoluTimes = ascending vector of days at which to evaluate solution
%
% DoseTimes = ascending vector of days at which to administer doses
% ( members of 'DoseTimes' must be members of 'SoluTimes' )
%
% Pars = vector of parameters (see list on line 91)
%
% Yic = vector of initial conditions (see list on line 172)
%
%% ========================================================================
% OPTIONS
% ========================================================================
%
% Opts = cell vector used to add optional inputs (see list on line 54)
% EX: Opts = {'AbsTol',1e-2} will set ode solver absolute tolerance to 1e-2
%
%% ========================================================================
% OUTPUTS
% ========================================================================
%
% Y_OUT(:,1) = SIV virions [log(#/mL)] at points in 'SoluTimes'
% Y_OUT(:,2) = total CD8+ T cells [#/uL] at points in 'SoluTimes'
%
%% ========================================================================
% FUNCTION
% ========================================================================
function Y_OUT = N803_model_2B(SoluTimes,DoseTimes,Pars,Yic,Opts)

% declare optional input default values
fullOut = false ;% if true, more outputs are added to Y_OUT (line 303)
timeOut = 30 ;% scalar time limit [seconds] for ODE solver
% if exceeded, an error is thrown
warnOff = true ;% if true, warning messages are suppressed
AbsTol = 1e-6 ;% absolute error tolerance (see 'odeset')
RelTol = 1e-3 ;% relative error tolerance (see 'odeset')
odeSolver = '15s' ;% ode solver function: 'ode15s' or 'ode23s'

% overwrite optional inputs, if specified
for n = 1:length(Opts)
if strcmp('fullOut', Opts(n)) ; fullOut = Opts{n+1} ; end
if strcmp('timeOut', Opts(n)) ; timeOut = Opts{n+1} ; end
if strcmp('warnOff', Opts(n)) ; warnOff = Opts{n+1} ; end
if strcmp('AbsTol', Opts(n)) ; AbsTol = Opts{n+1} ; end
if strcmp('RelTol', Opts(n)) ; RelTol = Opts{n+1} ; end
if strcmp('odeSolver', Opts(n)) ; odeSolver = Opts{n+1} ; end
end

% check input vector orientation
if size(SoluTimes ,1)==1 ; SoluTimes = SoluTimes' ; end
if size(DoseTimes ,1)==1 ; DoseTimes = DoseTimes' ; end
if size(Yic ,2)==1 ; Yic = Yic' ; end

% check that members of 'DoseTimes' are members of 'SoluTimes'
if length( unique( [ SoluTimes ; DoseTimes ] ) ) > length( SoluTimes )
error('All members of DoseTimes must be members of SoluTimes!')
end

% declare index in 'SoluTimes' for start, dose, and end times
Pieces = [ 1 ; find(ismember(SoluTimes,DoseTimes)) ; length(SoluTimes) ] ;

% suppress warnings (if requested)
if warnOff ; warning off ; end

%% ------------------------------------------------------------------------
% Rename inputed parameters -----------------------------------------------

q = Pars(01) ;% V growth rate constant [/d]
g = Pars(02) ;% E killing rate constant [uL/#-d]
V50E = Pars(03) ;% 50% V saturation for E activation [#/mL]
V50R = Pars(04) ;% 50% V saturation for R generation [#/mL]
a = Pars(05) ;% E activation rate constant [/d]
m = Pars(06) ;% E reversion rate constant [/d]
E50 = Pars(07) ;% 50% E saturation for proliferation [#/uL]
p = Pars(08) ;% E resting proliferation rate constant [/d]
pA = Pars(09) ;% E active proliferation rate constant [/d]
d = Pars(10) ;% E resting death rate constant [/d]
dA = Pars(11) ;% E active death rate constant [/d]
Xi = Pars(12) ;% X initial condition (at dose) [pmol/kg]
ka = Pars(13) ;% N803 absorption rate constant [/d]
ke = Pars(14) ;% N803 elimination rate constant [/d]
vd = Pars(15) ;% N803 'volume of distribution'/'bioavailability' [L/kg]
C50 = Pars(16) ;% 50% N803 saturation (drug effects) [pM]
ro = Pars(17) ;% E proliferation stimulation factor []
al = Pars(18) ;% E activation stimulation factor []
si = Pars(19) ;% R generation stimulation factor []
s = Pars(20) ;% R generation rate constant [/d]
dR = Pars(21) ;% R decay rate constant [/d]
la = Pars(22) ;% E killing regulation factor []
ph = Pars(23) ;% E proliferation regulation factor []
ze = Pars(24) ;% E activation regulation factor []

%% ------------------------------------------------------------------------
% Solve for each solution period ------------------------------------------

Y = zeros(Pieces(end),length(Yic)) ;% percolate output matrix

odeOptions = odeset('AbsTol',AbsTol,'RelTol',RelTol,'Jacobian',@jacob) ;

for n = 1:length(Pieces)-1

solvertime = tic ;% timer (for catching a stalled execution)

Piece = Pieces(n):Pieces(n+1) ;% index for solution piece

switch odeSolver % solve ODE for time window indexed by Piece
case '15s'
[~,Y_TEMP] = ode15s(@system,SoluTimes(Piece),Yic,odeOptions) ;
case '23s'
[~,Y_TEMP] = ode23s(@system,SoluTimes(Piece),Yic,odeOptions) ;
end

if length(Piece) == 2 % exclude timepoints added by ODE solver
Y_TEMP = [ Y_TEMP(1,:) ; Y_TEMP(end,:) ] ;
end

Y(Piece,:) = Y_TEMP ; % add 'Y_TEMP' to full solution 'Y'

Yic = Y_TEMP(end,:) ;% update ICs for next solution piece
Yic(07) = Yic(07) + Xi ;% add drug dose as initial condition
end

% Extremely small values can create numerical artifacts where the value
% becomes negative (mathematically impossible). This occurs for 'bad'
% parameter sets during calibration, and it occurs for the drug (where
% extremely small values are irrelevant). Taking the absolute value of Y
% prevents future errors with the 'log10' function.
Y = abs(Y) ;

Y_OUT = [ log10( Y(:,1) ) ... SIV virions [log(#/mL)]
, sum( Y(:,2:6), 2) ];% total CD8+ T cells [#/uL]

%% ------------------------------------------------------------------------
% Add columns to 'Y_OUT' (if requested) -----------------------------------

if fullOut
Y_OUT = outputs(Y_OUT,Y) ;
end

%% ========================================================================
% 'system' nested function
% ========================================================================
function dY = system(~,Y)

if toc(solvertime) >= timeOut
error('ODE solver stalled.')
end

V = Y(01) ;% SIV virions [#/mL]
E0 = Y(02) ;% resting CD8+ T cells [#/uL]
EA = Y(03:06) ;% active CD8+ T cells [#/uL] (E1 to E4)
X = Y(07) ;% N803 at absorption site [pmol/kg]
C = Y(08) ;% N803 plasma concentration [pM]
R = Y(09:10) ;% immune regulation []
F = terms(Y') ;% see 'terms' function

P = F(11) ;% E resting proliferation rate [/d]
A = F(12) ;% E activation rate [/d]
G = F(13) ;% E killing rate [#/uL-d]
S = F(14) ;% R generation rate [/d]

dY = Y ;% dY/dt

% calculate rates of change for V,E,X,C,R
dY(01) = q*V - G*sum(EA)*V ;% dV/dt
dY(02) = P*E0 - d *E0 + m *EA(4) - A *E0 ;% dE0/dt
dY(03) = - dA*EA(1) + 2*A *E0 - pA*EA(1) ;% dE1/dt
dY(4:5)= - dA*EA(2:3) + 2*pA*EA(1:2) - pA*EA(2:3) ;% dE2-3/dt
dY(06) = - dA*EA(4) + 2*pA*EA(3) - m *EA(4) ;% dE4/dt
dY(07) = - ka*X ;% dX/dt
dY(08) = ka*X/vd - ke*C ;% dC/dt
dY(09) = S - dR*R(1) ;% dR1/dt
dY(10) = dR*R(1) - dR*R(2) ;% dR2/dt
end
%% ========================================================================
% 'jacob' nested function
% ========================================================================
function J = jacob(~,Y)

V = Y(01) ;% SIV virions [#/mL]
E0 = Y(02) ;% resting CD8+ T cells [#/uL]
EA = sum(Y(03:06)) ;% active CD8+ T cells [#/uL]
C = Y(08) ;% N803 plasma concentration [pM]
F = terms(Y') ;% see 'terms' function

P = F(11) ;% E proliferation rate [/d]
A = F(12) ;% E activation rate [/d]
G = F(13) ;% E killing rate [#/uL-d]

J = zeros(10) ;

% Partial derivatives of P,A,G,S w.r.t. V,E,C,R
dPdR = - P * ph * F(05) ;
dAdR = - A * ze * F(06) ;
dGdR = - G * la * F(07) ;
dPdE = - P * F(08) / E50 ;
dAdV = A * V50E / V / (V50E+V) ;
dSdV = s * V50R / (V50R+V)^2 ;
dOdC = C50 / (C50+C)^2 ;
dPdC = p * F(08) * ro * dOdC * F(05) ;
dAdC = a * F(09) * al * dOdC * F(06) ;
dSdC = s * si * dOdC ;

% Partial derivatives of dV/dt
J(01,01 ) = q - G*EA ;% wrt V
J(01,03:06) = - G *V ;% wrt E
J(01,10 ) = - dGdR*EA*V ;% wrt R

% Partial derivatives of dE0/dt
J(02,01 ) = - dAdV*E0 ;% wrt V
J(02,02 ) = dPdE*E0 + P - d - A ;% wrt E0
J(02,03:05) = dPdE*E0 ;% wrt E1-E3
J(02,06 ) = dPdE*E0 + m ;% wrt E4
J(02,08 ) = dPdC*E0 - dAdC*E0 ;% wrt C
J(02,10 ) = dPdR*E0 - dAdR*E0 ;% wrt R

% Partial derivatives of dEA/dt
J(03,01 ) = 2*dAdV * E0 ;% dE1/dt wrt V
J(03,02 ) = 2*A ;% dE1/dt wrt E1
J(03,08 ) = 2*dAdC * E0 ;% dE1/dt wrt C
J(03,10 ) = 2*dAdR * E0 ;% dE1/dt wrt R
J(23:11:45) = - dA - pA ;% dEi/dt wrt Ei for i=1-3
J(24:11:46) = 2*pA ;% dEi/dt wrt Ei-1 for i=2-4
J(06,06 ) = - dA - m ;% dE8/dt wrt E4

% Partial derivatives of dX/dt
J(07,07) = -ka ;% wrt X

% Partial derivatives of dC/dt
J(08,07) = ka/vd ;% wrt X
J(08,08) = -ke ;% wrt C

% Partial derivatives of dR/dt
J(09,01) = dSdV ;% dR1/dt wrt V
J(09,08) = dSdC ;% dR1/dt wrt C
J(09,09) = -dR ;% dR1/dt wrt R1
J(10,09) = dR ;% dR2/dt wrt R1
J(10,10) = -dR ;% dR2/dt wrt R2
end
%% ========================================================================
% 'terms' nested function
% ========================================================================
function F = terms(Y)

V = Y(:,01) ;% SIV virions [#/mL]
E = sum(Y(:,02:06),2) ;% CD8+ T cells [#/uL]
C = Y(:,08) ;% N803 plasma concentration [pM]
R = Y(:,10) ;% immune regulation []

% percolate vector for storing terms
F = zeros(size(Y,1),14) ;

F(:,01) = C./ (C50+C) ;% N803 effect []

F(:,02) = 1 + ro*F(:,1) ;% E proliferation stimulation []
F(:,03) = 1 + al*F(:,1) ;% E activation stimulation []
F(:,04) = si*F(:,1) ;% N803-induced R generation []

F(:,05) = 1./ (1+ph*R) ;% E proliferation regulation []
F(:,06) = 1./ (1+ze*R) ;% E activation regulation []
F(:,07) = 1./ (1+la*R) ;% E killing regulation []

F(:,08) = E50./(E50+E) ;% E proliferation density dependence []
F(:,09) = V./ (V50E+V) ;% E activation antigen dependence []
F(:,10) = V./ (V50R+V) ;% R generation antigen dependence []

F(:,11) = p* F(:,08).* F(:,02).* F(:,05) ;% E proliferation [/d]
F(:,12) = a* F(:,09).* F(:,03).* F(:,06) ;% E activation rate [/d]
F(:,13) = g* F(:,07) ;% E killing rate [#/uL-d]
F(:,14) = s* (F(:,10)+F(:,04)) ;% R generation rate [/d]
end
%% ========================================================================
% 'outputs' nested function
% ========================================================================
function Y_OUT = outputs(Y_OUT,Y)

% NOTE: There are some zero elements of 'Y_OUT' in order to keep
% the indices matched with 'N803_model_2.m'.

Y_OUT = [ Y_OUT, Y, zeros( size(Y,1) , 47) ] ;

V = Y(:,01) ;% SIV virions [#/mL]
E0 = Y(:,02) ;% resting CD8+ T cells [#/uL]
EA = sum( Y(:,3:6) ,2) ;% active CD8+ T cells [#/uL]
Eaf = Y(:,6) ;% last generation of EA
Eap = EA - Eaf ;% EA that is proliferating

F = terms(Y) ;% see 'terms' function
P = F(:,11) ;% E resting proliferation rate [/d]
A = F(:,12) ;% E activation rate [/d]
G = F(:,13) ;% E killing rate [#/uL-d]
S = F(:,14) ;% R generation rate [/d]

Y_OUT(:,19) = F(:,01) ;% N803 effect []
Y_OUT(:,20) = E0+EA ;% total CD8+ T cells [#/uL]
Y_OUT(:,21) = E0+EA ;% total E
Y_OUT(:,23) = EA ;% active E
Y_OUT(:,26) = EA./(E0+EA) ;% activated frequency in E

Y_OUT(:,28) = P.*E0 ;% E0 proliferation term [#/uL/d]
Y_OUT(:,30) = d.*E0 ;% E0 death term [#/uL/d]
Y_OUT(:,32) = P ;% E0 proliferation rate [/d]
Y_OUT(:,34) = F(:,02) ;% C multiplier for P
Y_OUT(:,35) = F(:,05) ;% R multiplier for P
Y_OUT(:,36) = F(:,08) ;% density-dep multiplier for P

Y_OUT(:,38) = A.*E0 ;% E0 activation term [#/uL/d]
Y_OUT(:,40) = pA*Eap ;% EA proliferation term [#/uL/d]
Y_OUT(:,41) = dA*EA ;% EA death term [#/uL/d]
Y_OUT(:,43) = m*Eaf ;% EA reversion term [#/uL/d]
Y_OUT(:,45) = A ;% E0 activation rate [/d]
Y_OUT(:,47) = F(:,03) ;% C multiplier for A
Y_OUT(:,49) = F(:,06) ;% R multiplier for A
Y_OUT(:,51) = F(:,09) ;% V multiplier for A

Y_OUT(:,53) = q*V ;% V growth term [#/mL/d]
Y_OUT(:,54) = G.*EA.*V ;% E killing term [#/mL/d]
Y_OUT(:,55) = G.*EA ;% E killing rate [/d]
Y_OUT(:,56) = F(:,07) ;% R multiplier for G

Y_OUT(:,57) = S ;% R generation rate [/d]
Y_OUT(:,59) = s*si*F(:,01)./S ;% R gen by C fraction
end
end
Loading

0 comments on commit 4de9462

Please sign in to comment.