From 53272a93f8e25eb84dd7da9b1f4e8db5884cd690 Mon Sep 17 00:00:00 2001 From: "Cody, Jonathan William" Date: Thu, 28 Oct 2021 20:16:35 -0400 Subject: [PATCH] Add files via upload --- N803_model_2.m | 136 +++++++++++++++++++++++++++---------------------- 1 file changed, 76 insertions(+), 60 deletions(-) diff --git a/N803_model_2.m b/N803_model_2.m index 29d9357..36a425c 100644 --- a/N803_model_2.m +++ b/N803_model_2.m @@ -1,19 +1,15 @@ -% N803_model_2.m - solves host-level model for SIV virions and CD8+ T cells -% during a subcutaneously administered regimen of N803 - +%% N803_model_2.m - solves ODE model for N-803 treatment of SIV +% % /--------------------------------------------------------------\ -% | Date: 10/27/2021 | +% | Date: 10/28/2021 | % | Author: Jonathan Cody | -% | Affiliation: Pienaar Computational Systems Pharmacology Lab | +% | Affiliation: Purdue University | % | Weldon School of Biomedical Engineering | -% | Purdue University | +% | Pienaar Computational Systems Pharmacology Lab | % \--------------------------------------------------------------/ - -% Version 2D % -% TOODOO -% -> reconsider p1 calculation - +% Version 2D+ +% % Nomenclature: V = SIV virions [#/無] % E0 = resting SIV-specific CD8+ T cells [#/無] % Ea = active SIV-specific CD8+ T cells [#/無] @@ -23,43 +19,49 @@ % 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) - +% +%% ======================================================================== +% OPTIONS +% ======================================================================== +% NOTE: Set these inputs to [] to ignore or default them. +% EXAMPLE: Y_OUT = N803_model_2(SoluTimes,DoseTimes,Pars,[],[],[],[]) +% % 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) - +% +% fullOut = switch to determine output variables (see code) +% % 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) - +% (if empty, timeOut = 60) +% +% warnOff = switch to suppress display of errors +% %% ======================================================================== % OUTPUTS % ======================================================================== - +% % Y_OUT(:,1) = V at points in 'SoluTimes' [log fold change] % Y_OUT(:,2) = E at points in 'SoluTimes' [fold change] - -% If 'fullOut' is non-empty, more columns are included (see end of code) - +% +% If 'fullOut' is non-empty, more columns are included (see code) +% %% ======================================================================== % FUNCTION % ======================================================================== - -function Y_OUT = N803_model_2(SoluTimes,DoseTimes,Pars,EquatedPars,fullOut) +function Y_OUT = N803_model_2(SoluTimes,DoseTimes,Pars,... + EquatedPars,fullOut,timeOut,warnOff) % convert inputs to column vectors if size(SoluTimes ,1)==1 ; SoluTimes = SoluTimes' ; end @@ -72,6 +74,14 @@ % equate parameters in 'P' based on indices in 'EquatedPars' if ~isempty(EquatedPars) ; Pars = Pars(EquatedPars) ; end +% set default evaluation time out to 60 seconds +if isempty(timeOut) ; timeOut = 60 ; end + +% suppress warnings (if requested) +if ~isempty(warnOff) + warning off +end + %% ------------------------------------------------------------------------ % Rename inputed parameters ----------------------------------------------- @@ -100,16 +110,17 @@ aB1 = Pars(23) ;% B activation stimulation factor [] sE = 1e3 ;% R generation due to E activation [/d] sB = Pars(24) ;% R generation due to B activation [/d] -dR = Pars(25) ;% R decay rate constant [/d] -p2R = Pars(26) ;% initial E0/B0 proliferation regulation [] -aE2R = Pars(27) ;% initial E activation regulation [] -aB2R = Pars(28) ;% initial B activation regulation [] -gE2R = Pars(29) ;% initial E killing regulation [] -gB2R = Pars(30) ;% initial B killing regulation [] -N = Pars(31) ;% # of tolerance variables (delays + acting) -del = Pars(32) ;% tolerance rate constant [/d] -tau = Pars(33) ;% tolerance recovery [] -eta = Pars(34) ;% tolerance effect factor [] +nR = Pars(25) ;% # of regulation variables (delays + acting) (MUST BE >0) +dR = Pars(26) ;% R decay rate constant [/d] +p2R = Pars(27) ;% initial E0/B0 proliferation regulation [] +aE2R = Pars(28) ;% initial E activation regulation [] +aB2R = Pars(29) ;% initial B activation regulation [] +gE2R = Pars(30) ;% initial E killing regulation [] +gB2R = Pars(31) ;% initial B killing regulation [] +nT = Pars(32) ;% # of tolerance variables (delays + acting) (MUST BE >2) +dT = Pars(33) ;% tolerance rate constant [/d] +tau = Pars(34) ;% tolerance recovery [] +eta = Pars(35) ;% tolerance effect factor [] %% ------------------------------------------------------------------------ % Declare initial conditions & calculated parameters ---------------------- @@ -148,8 +159,8 @@ gB2 = gB2R/Ri ;% B killing regulation factor [] % declare initial conditions for Y (see 'system' function) -% V E0-8 B0-2 X C R tolerance -Yic = [ Vi Ei' Bi' 0 0 Ri zeros(1,N) ] ; +% V E0-8 B0-2 X C R T +Yic = [ Vi Ei' Bi' 0 0 Ri*ones(1,nR) zeros(1,nT) ] ; %% ------------------------------------------------------------------------ % Solve for each solution period ------------------------------------------ @@ -194,8 +205,8 @@ B0 = Y(:,11) ;% resting bystander CD8+ T cells [#/無] Ea = sum( Y(:,03:10) ,2) ;% active specific CD8+ T cells [#/無] Ba = sum( Y(:,12:13) ,2) ;% active bystander T cells [#/無] -R = Y(:,16) ;% immune regulation [] -F = terms(Y) ;% modifiers for 'gE0','gB0','p0','aE0','aB0' +R = Y(:,15+nR) ;% immune regulation [] +[F,T] = terms(Y) ;% modifiers for 'gE0','gB0','p0','aE0','aB0' Y_OUT(:,03) = (E0+Ea) / (E0(1)+Ea(1)) ;% E [fold change] Y_OUT(:,04) = (B0+Ba) / (B0(1)+Ba(1)) ;% B [fold change] @@ -210,7 +221,7 @@ Y_OUT(:,13) = Y(:,14) ;% X [pmol/kg] Y_OUT(:,14) = Y(:,15) ;% C [pM] Y_OUT(:,15) = R / R(1) ;% R [fold change] -Y_OUT(:,16) = (Y(:,15+N)+Y(:,16+N)) * (N~=0) ;% T [] +Y_OUT(:,16) = T ;% T [] Y_OUT(:,17) = F(:,01) ;% N803 effect [] Y_OUT(:,18) = F(:,02) ;% tolerance effect [] @@ -232,14 +243,14 @@ Y_OUT(:,34) = log10( F(:,18)/F(1,18) ) ;% lfc in regulation generation % Ba regulation stimulation proportion -Y_OUT(:,35) = sB*aB0*( F(:,14) + F(:,06) ) ./ F(:,18) ; +Y_OUT(:,35) = sB*aB0*(F(:,14)+F(:,06)) ./ F(:,18) ; % lfc in killing Kill = gB0*Ba./(1+gB2*R) + gE0*Ea./(1+gE2*R) ; Y_OUT(:,36) = log10( Kill / Kill(1) ) ; % Ba killing proportion -Y_OUT(:,37) = gB0*Ba./(1+gB2*R) ./ Kill ; +Y_OUT(:,37) = gB0*Ba./(1+gB2*R) ./ Kill ; %% ======================================================================== % NESTED FUNCTIONS @@ -266,7 +277,7 @@ function dY = system(~,Y) - if toc(solvertime) >= 10 + if toc(solvertime) >= timeOut error('ODE solver stalled.') end @@ -277,8 +288,6 @@ Ba = Y(12:13) ;% active bystander T cells [#/無] X = Y(14) ;% N803 at absorption site [pmol/kg] C = Y(15) ;% N803 plasma concentration [pM] - R = Y(16) ;% immune regulation [] - T = Y(17:16+N) ;% drug tolerance [] F = terms(Y') ;% modifiers for 'gE0','gB0','p0','aE0','aB0' gE = gE0*F(10) ;% specific killing rate [/d] @@ -300,24 +309,31 @@ dY(13) = 2*pB*Ba(1) - dA*Ba(2) - mB*Ba(2) ;% dB2/dt dY(14) = - ka*X ;% dX/dt dY(15) = ka*X/vd - ke*C ;% dC/dt - dY(16) = F(18) - dR*R ;% dR/dt - - % calculate rates of change for tolerance - if N ~= 0 - dY(17) = del*( F(1) - T(1) ) ;% T1/dt - dY(18:end-2) = del*( T(1:N-3) - T(2:N-2) ) ; - dY(end-1) = del*( T( N-2) - (1+tau)*T( N-1) ) ;% TN-1/dt - dY(end) = del*( tau*T( N-1) - tau *T( N ) ) ;% TN/dt + + % calculate rate of change for regulation variables + R = Y(16:15+nR) ;% immune regulation [] + dY(16) = F(18) - dR*R(1) ;% dR1/dt + if nR > 1 + dY(17:15+nR) = dR*( R(1:nR-1) - R(2:nR) ) ; + end + + % calculate rates of change for tolerance variables + T = Y(16+nR:15+nR+nT) ;% drug tolerance [] + if nT ~= 0 + dY(16+nR) = dT*( F(1) - T(1) ) ; + dY(17+nR:end-2) = dT*( T(1:nT-3) - T(2:nT-2) ) ; + dY(end-1) = dT*( T( nT-2) - (1+tau)*T( nT-1) ) ; + dY(end) = dT*( tau*T( nT-1) - tau *T( nT ) ) ; end end -function F = terms(Y) +function [F,T] = terms(Y) V = Y(:,01) ;% SIV virions [#/無] EB = sum(Y(:,02:13),2) ;% CD8+ T cells [#/無] C = Y(:,15) ;% N803 plasma concentration [pM] - R = Y(:,16) ;% immune regulation [] - T = Y(:,15+N)+Y(:,16+N) ;% drug tolerance [] + R = Y(:,15+nR) ;% immune regulation [] + T = (Y(:,14+nR+nT)+Y(:,15+nR+nT)) ;% drug tolerance [] % percolate vector for storing terms F = zeros(size(Y,1),18) ;