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 dc6a94f commit f2e491b
Showing 1 changed file with 284 additions and 0 deletions.
284 changes: 284 additions & 0 deletions N803_model_2.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,284 @@
% N803_model_2.m - solves host-level model for SIV virions, CD8+ T cells,
% and NK cells during a subcutaneously administered regimen of N803

% /--------------------------------------------------------------\
% | Date: 07/12/2021 |
% | Author: Jonathan Cody |
% | Affiliation: Pienaar Computational Systems Pharmacology Lab |
% | Weldon School of Biomedical Engineering |
% | Purdue University |
% \--------------------------------------------------------------/
%
% Version 2A
%
% Nomenclature: V = SIV virions [#/無]
% K = NK cells [#/無]
% E = CD8+ T cells [#/無] (E0 resting, EA active)
% X = N803 at absorption site [pmol/kg]
% C = N803 plasma concentration [pM]
% R = regulation [] (dimensionless quantity)
% T = tolerance [] (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')

% P = vector of parameters (see list in function)

% EquatedPars = vector of indices in 'P' used for equating parameter values
% (for example, EquatedPars(18) = 17 will set P(18) = P(17)
% (if empty, EquatedPars is ignored)

% timeOut = scalar time limit [seconds] for ODE solver 'ode15s'
% (if exceeded during 'ode15s' execution, an error is thrown)
% (if empty, timeOut = 10)

% fullOut = switch to determine output variables (see end of code)
% (if empty, fullOut is ignored)

%% ========================================================================
% OUTPUTS
% ========================================================================

% Y_OUT(:,1) = V at points in 'SoluTimes' [log fold change]
% Y_OUT(:,2) = E at points in 'SoluTimes' [fold change]
% Y_OUT(:,3) = K at points in 'SoluTimes' [fold change]

% If 'fullOut' is non-empty, more columns are included (see end of code)

%% ========================================================================
% FUNCTION
% ========================================================================

function Y_OUT = N803_model_2(SoluTimes,DoseTimes,Pars,EquatedPars,fullOut)

% convert inputs to column vectors
if size(SoluTimes ,1)==1 ; SoluTimes = SoluTimes' ; end
if size(DoseTimes ,1)==1 ; DoseTimes = DoseTimes' ; end
if size(Pars ,1)==1 ; Pars = Pars' ; end

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

% equate parameters in 'P' based on indices in 'EquatedPars'
if ~isempty(EquatedPars) ; Pars = Pars(EquatedPars) ; end

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

Vi = Pars(01) ;% V initial condition [log(#/mL)]
Ki = Pars(02) ;% K initial condition [#/無]
Ei = Pars(03) ;% E initial condition [#/無]
fEA = Pars(04) ;% active E initial frequency

gK_Ki = Pars(05)*Pars(06) ;% K initial killing rate [/d]
gE_EAi = Pars(06) ;% E initial killing rate [/d]
Km = Pars(07) ;% K max proliferating concentration [#/無]
Em = Pars(08) ;% E max proliferating concentration [#/無]
V50 = Pars(09) ;% V half-sat concentration for E activation [#/uL]
p = Pars(10) ;% EA proliferation rate [/d]
fm = Pars(11) ;% fraction of EA that revert to memory
dK = Pars(12) ;% K death rate constant [/d]
d0 = Pars(13) ;% E0 death rate constant [/d]
dA = Pars(14) ;% EA death rate constant [/d]

Xi = Pars(15) ;% X initial condition [pmol/kg]
ka = Pars(16) ;% N803 absorption rate constant [/d]
ke = Pars(17) ;% N803 elimination rate constant [/d]
vd = Pars(18) ;% N803 'volume of distribution'/'bioavailability' [L/kg]

C50 = Pars(19) ;% N803 half-saturation concentration [pM]
rK1_dK = Pars(20) ;% K maximum growth rate [/d]
rE1_d0 = Pars(21) ;% E0 maximum growth rate [/d]
a1 = Pars(22) ;% E activation stimulation factor []

dR = Pars(23) ;% regulation decay rate constant [/d]
rK2 = Pars(24) ;% K proliferation regulation factor []
rE2 = Pars(25) ;% E0 proliferation regulation factor []
a2 = Pars(26) ;% E activation regulation factor []

N = Pars(27) ;% # of tolerance variables (delays + acting)
del = Pars(28) ;% tolerance rate constant [/d]
tau = Pars(29) ;% tolerance recovery []
eta = Pars(30) ;% tolerance effect factor []

%% ------------------------------------------------------------------------
% Declare initial conditions & calculated parameters ----------------------

% some parameters are solved for by assuming an initially steady state

Vi = 10^(Vi-3) ;% V initial condition [#/無]
E0i = Ei*(1-fEA) ;% E0 initial condition [#/無]
EAi = Ei*( fEA) ;% EA initial condition [#/無] (sum of E1 to E8)

F5i = Vi/(Vi+V50) / (1+a2) ;% initial E activation rate modifiers []

% solve for activation rate a [/d] and E1 to E8 initial conditions [#/無]
MAT = diag( [0 -p*ones(1,7) -dA] ) ;
MAT(1,2:9) = ones(1,8) ;
MAT(2,1) = F5i * E0i ;
MAT(3:9,2:8) = MAT(3:9,2:8) + diag( 2*p*ones(1,7) ) ;
vec = [ EAi ; zeros(8,1) ] ;
sol = MAT\vec ;

a = sol(1) ;% E activation rate constant [/d]
pV = gE_EAi + gK_Ki ;% V proliferation rate constant [/d]
gE = gE_EAi/EAi ;% E killing rate constant []
gK = gK_Ki/Ki ;% K killing rate constant []
rK = (1+rK2) * (Km+Ki)/Km * dK ;% K proliferation rate constant [/d]
rE = (1+rE2) * (Em+Ei)/Em * ( d0 + a*F5i - fm*dA*sol(9)/E0i ) ;% E0 prolif
rK1 = rK1_dK/dK ;% K proliferation stimulation factor []
rE1 = rE1_d0/d0 ;% E0 proliferation stimulation factor []
psi = dR/F5i ;% R generation rate constant [/d]

% declare initial conditions for Y (see 'system' function)
% V K E0 E1-8 X C R tolerance
Yic = [ Vi Ki E0i sol(2:9)' 0 0 1 zeros(1,N) ] ;

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

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

for n = 1:length(Pieces)-1

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

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

[~,Y_TEMP] = ode15s(@system,SoluTimes(Piece),Yic) ;% solution piece

if length(Piece) == 2 % exclude timepoints added by ode15s
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(12) = Xi ;
end

%% ------------------------------------------------------------------------
% 'system' function -------------------------------------------------------

function dY = system(~,Y)

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

V = Y(01) ;% SIV virions [#/無]
K = Y(02) ;% NK cells [#/無]
E0 = Y(03) ;% resting CD8+ T cells [#/無]
EA = Y(04:11) ;% active CD8+ T cells [#/無]
X = Y(12) ;% N803 at absorption site [pmol/kg]
C = Y(13) ;% N803 plasma concentration [pM]
R = Y(14) ;% immune regulation []
F = terms(Y') ;% modifiers for 'rK','rE','a'

dY = Y ;% dY/dt

% calculate rates of change for V,K,E,X,C,R
dY(01) = pV*V - gK*K*V - gE*sum(EA)*V ;% dV/dt
dY(02) = rK*K *F(13) - dK*K ;% dK/dt
dY(03) = rE*E0*F(14) - d0*E0 + fm*dA*EA(8) - a*E0*F(15) ;% dE0/dt
dY(04) = - p*EA(1) + a*E0*F(15) ;% dE1/dt
dY(05:10) = 2*p*EA(1:6) - p*EA(2:7) ;% dE2/dt to dE7/dt
dY(11) = 2*p*EA(7) - dA*EA(8) ;% dE8/dt
dY(12) = - ka*X ;% dX/dt
dY(13) = ka*X/vd - ke*C ;% dC/dt
dY(14) = psi*F(15) - dR*R ;% dR/dt

% calculate rates of change for tolerance
if N ~= 0
dY(15) = del * ( F(1) - Y(15) ) ;
dY(16:end-2) = del * ( Y(15:end-3) - Y(16:end-2) ) ;
dY(end-1) = del * ( Y( end-2) - (1+tau)*Y( end-1) ) ;
dY(end) = del * ( tau*Y( end-1) - tau *Y( end ) ) ;
end
end

%% ------------------------------------------------------------------------
% 'terms' function --------------------------------------------------------

function F = terms(Y)

V = Y(:,01) ;% SIV virions [#/無]
K = Y(:,02) ;% NK cells [#/無]
E = sum(Y(:,03:11),2) ;% CD8+ T cells [#/無]
C = Y(:,13) ;% N803 plasma concentration [pM]
R = Y(:,14) ;% immune regulation []
T = Y(:,13+N)+Y(:,14+N) ;% drug tolerance []

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

F(:,01) = C ./ (C50+C) ;% N803 effect []
F(:,02) = 1 ./ (1+eta*T) ;% tolerance effect []
F(:,03) = F(:,1) .* F(:,2) ;% N803 effect (with tolerance) []
F(:,04) = 1 + rK1*F(:,3) ;% K proliferation stimulation []
F(:,05) = 1 + rE1*F(:,3) ;% E0 proliferation stimulation []
F(:,06) = 1 + a1*F(:,3) ;% activation stimulation []
F(:,07) = 1 ./ (1+rK2*R) ;% K proliferation regulation []
F(:,08) = 1 ./ (1+rE2*R) ;% E0 proliferation regulation []
F(:,09) = 1 ./ (1+ a2*R) ;% activation regulation []
F(:,10) = Km./ (Km +K) ;% K proliferation density dependence []
F(:,11) = Em./ (Em +E) ;% E0 proliferation density dependence []
F(:,12) = V ./ (V50+V) ;% activation antigen dependence []
F(:,13) = F(:,4) .* F(:,7) .* F(:,10) ;% K proliferation modifier []
F(:,14) = F(:,5) .* F(:,8) .* F(:,11) ;% E0 proliferation modifier []
F(:,15) = F(:,6) .* F(:,9) .* F(:,12) ;% activation modifier []

end

%% ------------------------------------------------------------------------
% Prepare main outputs 'Y_OUT' -------------------------------------------

Y = abs(Y) ;% removing negatives resulting from numerical error

% V,E,K [fold change]
Y_OUT = [ log10(Y(:,1)/Vi) ... V [log fold change]
, sum(Y(:,3:11),2)/Ei ... E [fold change]
, Y(:,2)/Ki ];% K [fold change]

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

if isempty(fullOut) ; return ; end

Y_OUT = [ Y_OUT zeros(Pieces(end),20-3) ] ;

Y_OUT(:,4) = Y(:,03)/E0i ;% E0 [fold change]
Y_OUT(:,5) = sum(Y(:,04:11),2)/EAi ;% EA [fold change]
Y_OUT(:,6) = Y(:,12) ;% X [pmol/kg]
Y_OUT(:,7) = Y(:,13) ;% C [pM]
Y_OUT(:,8) = Y(:,14) ;% R []
Y_OUT(:,9) = (Y(:,13+N)+Y(:,14+N)) * (N~=0) ;% T []

F = terms(Y) ;% modifiers for 'rK','rE','a'

Y_OUT(:,10) = F(:,01) ;% N803 effect []
Y_OUT(:,11) = F(:,02) ;% tolerance effect []
Y_OUT(:,12) = F(:,03) ;% N803 effect (with tolerance) []
Y_OUT(:,13) = F(:,04) ;% K proliferation stimulation []
Y_OUT(:,14) = F(:,05) ;% E0 proliferation stimulation []
Y_OUT(:,15) = F(:,06) ;% activation stimulation []
Y_OUT(:,16) = F(:,07)/F(1,07) ;% K proliferation regulation [fold change]
Y_OUT(:,17) = F(:,08)/F(1,08) ;% E0 proliferation regulation [fold change]
Y_OUT(:,18) = F(:,09)/F(1,09) ;% activation regulation [fold change]
Y_OUT(:,19) = F(:,10)/F(1,10) ;% K prolif density dependence [fold change]
Y_OUT(:,20) = F(:,11)/F(1,11) ;% E0 prolif density dependence [fold change]
Y_OUT(:,21) = F(:,12)/F(1,12) ;% activ antigen dependence [fold change]
Y_OUT(:,22) = F(:,13)/F(1,13) ;% K proliferation modifier [fold change]
Y_OUT(:,23) = F(:,14)/F(1,14) ;% E0 proliferation modifier [fold change]
Y_OUT(:,24) = F(:,15)/F(1,15) ;% activation modifier [fold change]

% NK killing proportion
Y_OUT(:,25) = gK*Y(:,2) ./ ( gK*Y(:,2) + gE*sum(Y(:,04:11),2) ) ;

end

0 comments on commit f2e491b

Please sign in to comment.