Add files via upload

ElsjePienaarGroup · Oct 29, 2021 · 53272a9 · 53272a9
1 parent 3f763ad
commit 53272a9
Showing 1 changed file with 76 additions and 60 deletions.
diff --git 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) ;