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

% /--------------------------------------------------------------\
% | Date: 06/29/2019 |
% | Author: Jonathan Cody |
% | Affiliation: Pienaar Computational Systems Pharmacology Lab |
% | Weldon School of Biomedical Engineering |
% | Purdue University |
% \--------------------------------------------------------------/

% Nomenclature: V = SIV virions (dominant strain) [#/無]
% W = SIV virions (escape strain) [#/無]
% E = CD8+ T cells [#/無]
% K = NK cells [#/無]
% X = N803 at absorption site [pmol/kg]
% C = N803 plasma concentration [pM]
% TOL = tolerance [] (dimensionless quantity)
% REG = regulation [] (dimensionless quantity)

%% ========================================================================
% REQUIRED 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 corresponding to list below
% (parameter pairs are for E and K, respectively)

% P(01) = V+W initial condition [log(#/mL)] (will convert to [#/無])
% P(02) = W initial frequency
% P(03) = E initial condition [#/無]
% P(04) = K initial condition [#/無]

% P(05) = X initial condition [pmol/kg]
% P(06) = N803 absorption rate constant [/d]
% P(07) = N803 elimination rate constant [/d]
% P(08) = N803 'volume of distribution'/'bioavailability' [L/kg]

% P(09) = N803 half-saturation concentration [pM]
% P(10) = # of tolerance variables (delays + acting)
% P(11) = # of regulation variables (delays + acting)
% P(12) = tolerance rate constant [/d]
% P(13) = regulation rate constant [/d]
% P(14) = W targeting (by E) factor []

% P(15) = killing rate constant [無/#-d]
% P(16) [*P(15)] (will convert to [無/#-d])
% P(17) = death rate constant [/d]
% P(18)
% P(19) = max proliferating concentration [#/無]
% P(20)

% P(21) = killing stimulation factor []
% P(22)
% P(23) = proliferation stimulation factor []
% P(24)
% P(25) = tolerance effect factor []
% P(26)

% P(27) = killing regulation factor []
% P(28)
% P(29) = proliferation regulation factor []
% P(30)

%% ========================================================================
% OPTIONAL INPUTS
% ========================================================================

% 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)

% normOut = scalar used to change units of Y_MAIN (see OUTPUTS)
% (if empty, normOut is ignored)

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

% if normOut = 1
% Y_MAIN(:,1) = V+W at points in 'SoluTimes' [log(#/mL)] [log fold change]
% Y_MAIN(:,2) = E at points in 'SoluTimes' [#/無] [fold change]
% Y_MAIN(:,3) = K at points in 'SoluTimes' [#/無] [fold change]
% (see INPUTS)

% extra = structure variable of extra outputs (see end of file)

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

function [Y_MAIN,extra] = N803_model_1(SoluTimes,DoseTimes,P,...
EquatedPars,timeOut,normOut)

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

% declare reference index in 'SoluTimes' containing start time, all dose
% times, and end time (used for piecewise solution of ODEs)
PieceIn = [ 1 ; find(ismember(SoluTimes,DoseTimes)) ; length(SoluTimes) ] ;

% declare index in Y (see 'system' function) for C, last tolerance
% variable, and last regulation variable
TolRegIn = [ 6 ; 6+P(10) ; 6+P(10)+P(11) ] ;

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

% set default 'timeOut' to 10 seconds
if isempty(timeOut) ; timeOut = 10 ; end

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

% convert parameters
P(01) = 10^(P(01)-3) ;% [log(#/mL)] to [#/無]
P(16) = P(16)*P(15) ;% [*P(15)] to [#/無-d]

% declare initial conditions for Y (see 'system' function)
Yic = zeros(1,TolRegIn(3)) ;% zero except for...
Yic(1:4) = [ P(01)*[1-P(02) P(02)] P(03:04)' ] ;% V,W,E,K

% declare proliferation rate constants 'q_V','q_W','r_E','r_K' [/d]
qVW = [ 1 1 ; P(14) 1 ] * ( P(15:16).* P(03:04) ) ;% q for [V;W]
rEK = P(17:18).* ( P(19:20) + P(03:04) )./ P(19:20) ;% r for [E;K]

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

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

for n = 1:length(PieceIn)-1

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

% solve ODEs (see MATLAB documentation for 'ode15s')
[~,Y_TEMP] = ode15s(@system,SoluTimes(PieceIn(n):PieceIn(n+1)),Yic) ;

% for solutions of 2 timepoints, exclude the transition timepoints that
% ode15s adds by default
if length(SoluTimes(PieceIn(n):PieceIn(n+1))) == 2
Y_TEMP = [ Y_TEMP(1,:) ; Y_TEMP(end,:) ] ;
end

% add 'Y_TEMP' to full solution 'Y'
Y(PieceIn(n):PieceIn(n+1),:) = Y_TEMP ;

% update initial condition for next solution window
Yic = Y(PieceIn(n+1),:) + [ 0 0 0 0 P(05) Yic(6:end)*0 ] ;
end

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

function dY = system(~,Y)

% Y(1) = V = SIV virions (dominant strain) [#/無]
% Y(2) = W = SIV virions (escape strain) [#/無]
% Y(3) = E = CD8+ T cells [#/無]
% Y(4) = K = NK cells [#/無]
% Y(5) = X = N803 at absorption site [pmol/kg]
% Y(6) = C = N803 plasma concentration [pM]
% Y(7 :5+P(10) ) = tolerance delay variables
% Y( 6+P(10) ) = tolerance acting variable
% Y(7+P(10):5+P(10)+P(11)) = regulation delay variables
% Y( 6+P(10)+P(11)) = regulation acting variable

% throw error if 'solvertime' exceeds 'timeOut' limit
if toc(solvertime) >= timeOut
error('ODE solver stalled.')
end

% percolate dY/dt
dY = Y ;

%% --------------------------------------------------------------------
% Update rates of change for V,W,E,K,X,C ------------------------------

% call 'terms' function to calculate updated terms
Y_terms = terms(Y) ;

% Proliferation/Transfer - Killing/Decay/Transfer
dY(1:2) = qVW.*Y(1:2) - ([ 1 1 ; ...
P(14) 1 ] * Y_terms(8:9)).*Y(1:2) ;% dV:W/dt
dY(3:4) = Y_terms(10:11).*Y(3:4) - P(17:18).*Y(3:4) ;% dE:K/dt
dY(5) = - P(06)*Y(5) ;% dX/dt
dY(6) = P(06)/P(08)*Y(5) - P(07)*Y(6) ;% dC/dt

%% --------------------------------------------------------------------
% Calculate rates of change for tolerance & regulation variables ------

% for tolerance and regulation
for m = 1:2
% if number of variables is not zero
if P(09+m) ~= 0
% update variable rates of change
dY(TolRegIn(m)+1:TolRegIn(m+1)) = P(11+m) * ...
( [ Y_terms(1) ; Y(TolRegIn(m)+1:TolRegIn(m+1)-1) ] ...
- Y(TolRegIn(m)+1:TolRegIn(m+1)) ) ;
end
end

end

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

function Y_terms = terms(Y)

% percolate vector for storing terms
Y_terms = zeros(12,1) ;

% Y_terms(1) = N803 effect []
Y_terms(1) = Y(6) / ( P(09) + Y(6) ) ;

% Y_terms(2) = N803 effect for E (with tolerance) []
% Y_terms(3) = N803 effect for K (with tolerance) []
Y_terms(2:3) = [Y_terms(1);Y_terms(1)]./ ...
( 1 + P(25:26) * Y(TolRegIn(2)) ) ;

% Y_terms(4) = killing modifier for E (with effect & regulation) []
% Y_terms(5) = killing modifier for K (with effect & regulation) []
% Y_terms(6) = proliferation modifier for E (w/ effect & regulation) []
% Y_terms(7) = proliferation modifier for K (w/ effect & regulation) []
Y_terms(4:7) = ( 1 + P(21:24).* [Y_terms(2:3);Y_terms(2:3)] )./ ...
( 1 + P(27:30) * Y(TolRegIn(3)) ) ;

% Y_terms(8) = total killing rate applied to V by E [/d]
% Y_terms(9) = total killing rate applied to V by K [/d]
Y_terms(8:9) = P(15:16).* Y_terms(4:5).* Y(3:4) ;

% Y_terms(10) = proliferation rate of E (with density-dependence) [/d]
% Y_terms(11) = proliferation rate of K (with density-dependence) [/d]
Y_terms(10:11) = rEK.* Y_terms(6:7).* P(19:20)./ (P(19:20) + Y(3:4)) ;

end

%% ------------------------------------------------------------------------
% Prepare main outputs 'Y_MAIN' -------------------------------------------

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

if normOut == 1 % V+W,E,K [fold change]
Y_MAIN = [ (Y(:,1)+Y(:,2)) / P(01) , Y(:,3:4)./ P(03:04)' ] ;

else % V+W [#/mL] E,K [#/無]
Y_MAIN = [ (Y(:,1)+Y(:,2)) * 1000 , Y(:,3:4) ] ;
end

Y_MAIN(:,1) = log10( Y_MAIN(:,1) ) ;% convert to log-scale

%% ------------------------------------------------------------------------
% Prepare 'extra' outputs data structure ----------------------------------

% if only one output is requested ('Y_MAIN'), exit function
if nargout == 1 ; return ; end

% percolate matrix for storing terms
Y_TERMS = zeros(PieceIn(end),11) ;

% populate 'Y_TERMS'
for n = 1:PieceIn(end)
Y_TERMS(n,:) = terms(Y(n,:)')' ;
end

% store extra variables in 'extra' structure
extra = struct('Y',{Y}, ... all model variables (see 'system')
'Y_TERMS',{Y_TERMS}, ... extra model terms (see 'terms')
'P',{P}, ... model parameters (with conversions)
'qVW',{qVW}, ... rate constants [ 'q_V' ; 'q_W' ]
'rEK',{rEK} );% rate constants [ 'r_E' ; 'r_K' ]

end

0 comments on commit d343a05

Please sign in to comment.