From ca3afa2859d2e008527a872518c4641af6ab7f5f Mon Sep 17 00:00:00 2001 From: "Cody, Jonathan William" Date: Mon, 9 Jan 2023 10:38:48 -0600 Subject: [PATCH] Delete N803_model_1.m --- N803_model_1.m | 297 ------------------------------------------------- 1 file changed, 297 deletions(-) delete mode 100644 N803_model_1.m diff --git a/N803_model_1.m b/N803_model_1.m deleted file mode 100644 index 39659b5..0000000 --- a/N803_model_1.m +++ /dev/null @@ -1,297 +0,0 @@ -% 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: 12/05/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) = maximum expansion rate [/d] (will convert to rho) -% P(24) -% P(25) = tolerance effect factor [] -% P(26) - -% P(27) = killing regulation factor [] -% P(28) -% P(29) = proliferation regulation factor [] -% P(30) -% P(31) = tolerance recovery [] - -%% ======================================================================== -% 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) - -warning off - -% 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 tolerance variables -TolIn = 7:(6+P(10)) ; - -% declare index in Y (see 'system' function) for regulation variables -RegIn = (7+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] - -% convert maximum expansion rate [/d] to prolif. stimulation factor [] -P(23:24) = P(23:24)./P(17:18) ; - -% declare initial conditions for Y (see 'system' function) -Yic = zeros(1,6+P(10)+P(11)) ;% 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 : (6+P(10)) ) = tolerance variables - % Y( (7+P(10)) : (6+P(10)+P(11)) ) = regulation variables - - % 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 ------ - - % if number of tolerance variables is not zero - if P(10) ~= 0 - % update variable rates of change - dY(TolIn(1:end-2)) = P(12) * [ Y_terms(1) ; Y(TolIn(1:end-3)) ] ... - - P(12) * Y(TolIn(1:end-2)) ; - dY(TolIn(end-1)) = P(12) * Y(TolIn(end-2)) ... - - P(12) * (1+P(31)) * Y(TolIn(end-1)) ; - dY(TolIn(end)) = P(12) * P(31) * Y(TolIn(end-1)) ... - - P(12) * P(31) * Y(TolIn(end)) ; - end - - % if number of regulation variables is not zero - if P(11) ~= 0 - % update variable rates of change - dY(RegIn) = P(13) * [ Y_terms(1) ; Y(RegIn(1:end-1)) ] ... - - P(13) * Y(RegIn) ; - end - -end - -%% ------------------------------------------------------------------------ -% 'terms' function -------------------------------------------------------- - -function Y_terms = terms(Y) - - % percolate vector for storing terms - Y_terms = zeros(11,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) * sum(Y((5:6)+P(10))) ) ; - - % 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(6+P(10)+P(11)) ) ; - - % 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