diff --git a/N803_model_1.m b/N803_model_1.m new file mode 100644 index 0000000..2ba0e1b --- /dev/null +++ b/N803_model_1.m @@ -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 \ No newline at end of file