Add files via upload

N803_model_2S.m

@@ -0,0 +1,353 @@
+%% N803_model_2S.m - solves ODE model for N-803 treatment of SIV
+%
+%     /--------------------------------------------------------------\
+%     | Date:         12/27/2022                                     |
+%     | Author:       Jonathan Cody                                  |
+%     | Affiliation:  Purdue University                              |
+%     |               Weldon School of Biomedical Engineering        |
+%     |               Pienaar Computational Systems Pharmacology Lab |
+%     \--------------------------------------------------------------/
+%
+% Model 2S is a simplified version of Model 2 that combines SIV-specific
+% and non-SIV-specific CD8+ T cells
+%
+% Nomenclature: V  = SIV virions [#/無]
+%               T8 = total CD8+ T cells [#/無]
+%               E0 = resting CD8+ T cells [#/無]
+%               EA = active  CD8+ T cells [#/無]
+%               X  = N803 at absorption site [pmol/kg]
+%               C  = N803 plasma concentration [pM]
+%               R  = regulation [] (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')
+%
+% Pars = vector of parameters (see list on line 93)
+%
+% Yic = vector of initial conditions (see list on line 230)
+%
+%% ========================================================================
+%	OPTIONS
+%  ========================================================================
+%
+% Opts = cell vector used to add optional inputs (see list on line 54)
+% EX: Opts = {'AbsTol',1e-2} will set ode solver absolute tolerance to 1e-2
+%
+%% ========================================================================
+%	OUTPUTS
+%  ========================================================================
+%
+% Y_OUT(:,1) = V at points in 'SoluTimes' [log fold change]
+% Y_OUT(:,2) = T8 at points in 'SoluTimes' [fold change]
+%
+%% ========================================================================
+%	FUNCTION
+%  ========================================================================
+function Y_OUT = N803_model_2S(SoluTimes,DoseTimes,Pars,Yic,Opts)
+
+% Declare optional input default values
+EqualPars = []     ;% vector of indices in 'Pars' for equating parameters
+                    % EX: EqualPars(28) = 27 will set P(28) = P(27)
+fullOut   = []     ;% if non-empty, more outputs are added (see line 170)
+timeOut   = 60     ;% scalar time limit [seconds] for ODE solver
+                    % if exceeded, an error is thrown
+warnOff   = []     ;% if non-empty, warning messages are suppressed
+AbsTol    = 1e-6   ;% absolute error tolerance (see 'odeset')
+RelTol    = 1e-3   ;% relative error tolerance (see 'odeset')
+odeSolver = '15s'  ;% ode solver function: 'ode15s' or 'ode23s'
+
+% Overwrite optional inputs, if specified
+for n = 1:length(Opts)
+if strcmp('EqualPars', Opts(n))  ; EqualPars = Opts{n+1}  ; end
+if strcmp('fullOut',   Opts(n))  ; fullOut   = Opts{n+1}  ; end
+if strcmp('timeOut',   Opts(n))  ; timeOut   = Opts{n+1}  ; end
+if strcmp('warnOff',   Opts(n))  ; warnOff   = Opts{n+1}  ; end
+if strcmp('AbsTol',    Opts(n))  ; AbsTol    = Opts{n+1}  ; end
+if strcmp('RelTol',    Opts(n))  ; RelTol    = Opts{n+1}  ; end
+if strcmp('odeSolver', Opts(n))  ; odeSolver = Opts{n+1}  ; end
+end
+
+% check inputs vector orientation
+if size(SoluTimes ,1)==1 ; SoluTimes = SoluTimes' ; end
+if size(DoseTimes ,1)==1 ; DoseTimes = DoseTimes' ; end
+if size(Yic       ,2)==1 ; Yic = Yic'             ; end
+
+% declare index in 'SoluTimes' for start, dose, and end times
+Pieces = [ 1 ; find(ismember(SoluTimes,DoseTimes)) ; length(SoluTimes) ] ;
+
+% equate parameters in 'P' based on indices in 'EquatedPars'
+if ~isempty(EqualPars) ;  Pars = Pars(EqualPars)  ; end 
+
+% suppress warnings (if requested)
+if ~isempty(warnOff) ;  warning off  ; end
+
+%% ------------------------------------------------------------------------
+% Rename inputed parameters -----------------------------------------------
+
+q    = Pars(01) ;% V growth rate (if E+B were absent) [/d]
+g0   = Pars(02) ;% EA killing rate constant [無/#-d]
+V50E = Pars(03) ;% 50% viral stimulation saturation for E [#/無]
+V50R = Pars(04) ;% 50% viral stimulation saturation for R [#/無]
+a0   = Pars(05) ;% E0 activation rate constant [/d]
+m    = Pars(06) ;% EA reversion rate constant [/d]
+
+E50  = Pars(07) ;% 50% E0 proliferation saturation [#/無]
+p0   = Pars(08) ;% E0 proliferation rate constant [/d]
+pA   = Pars(09) ;% EA proliferation rate constant [/d]
+d    = Pars(10) ;% E0 death rate constant [/d]
+dA   = Pars(11) ;% EA death rate constant [/d]
+
+Xi   = Pars(12) ;% X initial condition [pmol/kg]
+ka   = Pars(13) ;% N803 absorption rate constant [/d]
+ke   = Pars(14) ;% N803 elimination rate constant [/d]
+vd   = Pars(15) ;% N803 'volume of distribution'/'bioavailability' [L/kg]
+C50  = Pars(16) ;% 50% N803 stimulation concentration [pM]
+p1   = Pars(17) ;% E0 proliferation stimulation factor []
+a1   = Pars(18) ;% E0 activation stimulation factor []
+s1   = Pars(19) ;% R generation stimulation factor []
+
+s0   = Pars(20) ;% R generation maximum base rate [/d]
+dR   = Pars(21) ;% R decay rate constant [/d]
+g2   = Pars(22) ;% EA killing regulation factor []
+p2   = Pars(23) ;% E0 proliferation regulation factor []
+a2   = Pars(24) ;% E0 activation regulation factor []
+
+%% ------------------------------------------------------------------------
+% Solve for each solution period ------------------------------------------
+
+Y = zeros(Pieces(end),length(Yic)) ;% percolate output matrix
+
+odeOptions = odeset('AbsTol',AbsTol,'RelTol',RelTol,'Jacobian',@jacob) ;
+
+for n = 1:length(Pieces)-1
+
+    solvertime = tic ;% timer (for catching a stalled execution)
+
+    Piece = Pieces(n):Pieces(n+1) ;% index for solution piece
+
+    switch odeSolver  % solve ODE for time window indexed by Piece
+        case '15s'
+            [~,Y_TEMP] = ode15s(@system,SoluTimes(Piece),Yic,odeOptions) ;
+        case '23s'
+            [~,Y_TEMP] = ode23s(@system,SoluTimes(Piece),Yic,odeOptions) ;
+    end
+
+    if length(Piece) == 2 % exclude timepoints added by ODE solver
+        Y_TEMP = [ Y_TEMP(1,:) ; Y_TEMP(end,:) ] ;
+    end
+
+    Y(Piece,:) = Y_TEMP ; % add 'Y_TEMP' to full solution 'Y'
+
+    Yic     = Y_TEMP(end,:) ;% update ICs for next solution piece
+    Yic(11) = Yic(11) + Xi  ;% add drug dose as initial condition
+end
+
+%% ------------------------------------------------------------------------
+% Prepare main outputs 'Y_OUT' --------------------------------------------
+
+Y = abs(Y) ;% removing negatives resulting from numerical error
+
+Y_OUT = [ log10(     Y(:,1)       /     Y(1,1)     )... V [log fold change]
+        ,        sum(Y(:,2:10),2) / sum(Y(1,2:10))  ];% E+B [fold change]
+
+%% ------------------------------------------------------------------------
+% Add columns to 'Y_OUT' (if requested) -----------------------------------
+
+if isempty(fullOut) ; return ; end
+
+Y_OUT = [ Y_OUT,  Y,  zeros(Pieces(end),27) ] ;
+
+V  = Y(:,01)  ;% SIV virions [#/無]
+E0 = Y(:,02)  ;% resting specific CD8+ T cells [#/無]
+EA = sum( Y(:,03:10) ,2)  ;% active specific CD8+ T cells [#/無]
+Eaf = Y(:,10)   ;% last generation of EA
+Eap = EA - Eaf  ;% EA that is proliferating
+F  = terms(Y')  ;% modifiers for 'p0','a0','s0','g0'
+
+p = F(11)  ;% E0 proliferation rate [/d]
+a = F(12)  ;% E0 activation rate [/d]
+s = F(13)  ;% regulation generation rate [/d]
+g = F(14)  ;% EA killing rate [#/無-d]
+
+Y_OUT(:,17) = F(:,01)        ;% N803 effect []
+
+Y_OUT(:,18) = (E0+EA)  ;% CD8+ T cells [#/無]
+Y_OUT(:,19) = EA       ;% EA
+Y_OUT(:,20) = (EA)./(E0+EA)  ;% activated frequency in E
+Y_OUT(:,21) = log10( EA )    ;% EA (log)
+
+Y_OUT(:,22) = p.*E0    ;% E0 proliferation term [#/無/d]
+Y_OUT(:,23) = d.*E0    ;% E0 death term [#/無/d]
+Y_OUT(:,24) = p        ;% E0 proliferation rate [/d]
+Y_OUT(:,25) = F(:,02)  ;% C multiplier for p* (log)
+Y_OUT(:,26) = F(:,05)  ;% R multiplier for p* (log)
+Y_OUT(:,27) = F(:,08)  ;% density-dep multiplier for p* (log)
+
+Y_OUT(:,28) = a.*E0    ;% E0 activation term [#/無/d]
+Y_OUT(:,29) = pA*Eap   ;% EA proliferation term [#/無/d]
+Y_OUT(:,30) = dA*EA    ;% EA death term [#/無/d]
+Y_OUT(:,31) = m*Eaf    ;% EA reversion term [#/無/d]
+Y_OUT(:,32) = a        ;% E0 activation rate [/d]
+Y_OUT(:,33) = F(:,03)  ;% C multiplier for a*
+Y_OUT(:,34) = F(:,06)  ;% R multiplier for a*
+Y_OUT(:,35) = F(:,09)  ;% V multiplier for a*
+
+Y_OUT(:,36) = q*V        ;% V growth term [#/無/d]
+Y_OUT(:,37) = g.*EA.*V   ;% E killing term [#/無/d]
+Y_OUT(:,38) = g.*EA      ;% E killing rate [/d]
+Y_OUT(:,39) = F(:,07)    ;% R multiplier for g*
+
+Y_OUT(:,40) = s                     ;% R generation
+Y_OUT(:,41) = 1 + F(:,04)./F(:,10)  ;% C multiplier for s*
+Y_OUT(:,42) = F(:,10)               ;% V multiplier for s*
+
+% Converting most outputs to log-value
+logID = [ 3:12 , 15:16 , 22:42 ]  ;
+Y_OUT(:,logID) = log10( Y_OUT(:,logID) )  ;
+
+% Converting virus to #/mL
+Y_OUT(:,3) = Y_OUT(:,3) + 3 ;
+
+%% ========================================================================
+%	NESTED FUNCTIONS
+%  ========================================================================
+
+function dY = system(~,Y)
+
+    if toc(solvertime) >= timeOut
+        error('ODE solver stalled.')
+    end
+
+    V  = Y(01)      ;% SIV virions [#/無]
+    E0 = Y(02)      ;% resting specific CD8+ T cells [#/無]
+    EA = Y(03:10)   ;% active specific CD8+ T cells [#/無] (E1 to E8)
+    X  = Y(11)      ;% N803 at absorption site [pmol/kg]
+    C  = Y(12)      ;% N803 plasma concentration [pM]
+    R  = Y(13:14)   ;% immune regulation []
+    F  = terms(Y')  ;% modifiers for 'p0','a0','s0','g0'
+
+    p = F(11)  ;% E0 proliferation rate [/d]
+    a = F(12)  ;% E0 activation rate [/d]
+    s = F(13)  ;% regulation generation rate [/d]
+    g = F(14)  ;% EA killing rate [#/無-d]
+
+    dY = Y          ;% dY/dt
+
+    % calculate rates of change for V,E,B,X,C,R
+    dY(01) = q*V - g*sum(EA)*V  ;% dV/dt
+    dY(02) = p*E0  - a*E0  - d*E0        + m*EA(8)     ;% dE0/dt
+    dY(03) =     + 2*a*E0  - dA*EA(1)    - pA*EA(1)    ;% dE1/dt
+    dY(4:9)= 2*pA*EA(1:6)  - dA*EA(2:7)  - pA*EA(2:7)  ;% dE2-7/dt
+    dY(10) = 2*pA*EA(7)    - dA*EA(8)    - m*EA(8)     ;% dE8/dt
+    dY(11) =         - ka*X     ;% dX/dt
+    dY(12) = ka*X/vd - ke*C     ;% dC/dt
+    dY(13) = s       - dR*R(1)  ;% dR1/dt
+    dY(14) = dR*R(1) - dR*R(2)  ;% dR2/dt
+
+end
+
+function J = jacob(~,Y)
+
+    V  = Y(01)      ;% SIV virions [#/無]
+    E0 = Y(02)      ;% resting specific CD8+ T cells [#/無]
+    EA = sum(Y(03:10)) ;% active specific CD8+ T cells [#/無] (E1 to E8)
+    C  = Y(12)      ;% N803 plasma concentration [pM]
+    F  = terms(Y')  ;% modifiers for 'p0','a0','s0','g0'
+
+    p = F(11)  ;% E0 proliferation rate [/d]
+    a = F(12)  ;% E0 activation rate [/d]
+    g = F(14)  ;% EA killing rate [#/無-d]
+
+    J = zeros(14) ;
+
+    % Partial derivatives of g,p,a,s w.r.t. V,E,C,R
+    dpdR = - p * p2 * F(05) ;
+    dadR = - a * a2 * F(06) ;
+    dgdR = - g * g2 * F(07) ;
+    dpdE = - p * F(08) / E50 ;
+    dadV = a0 * V50E / (V50E+V)^2 * F(03) * F(06) ;
+    dsdV = s0 * V50R / (V50R+V)^2 ;
+    dOdC = C50 / (C50+C)^2 ;
+    dpdC = p0 * F(08) * p1 * dOdC * F(05) ;
+    dadC = a0 * F(09) * a1 * dOdC * F(06) ;
+    dsdC = s0         * s1 * dOdC ;
+
+    % Partial derivatives of f1 = dV/dt
+    J(01,01   ) = q - g*EA ;
+    J(01,02:10) = - g * V ;
+    J(01,14   ) = - dgdR*EA * V ;
+
+    % Partial derivatives of f2 = dE0/dt
+    J(02,01   ) = - dadV * E0 ;
+    J(02,02   ) = p + dpdE*E0 - d - a ;
+    J(02,03:09) = dpdE*E0 ;
+    J(02,10   ) = dpdE*E0 + m ;
+    J(02,12   ) = ( dpdC - dadC ) * E0 ;
+    J(02,14   ) = ( dpdR - dadR ) * E0 ;
+
+    % Partial derivatives of f3 = dE1/dt (save f3/dE1)
+    J(03,01   ) = 2*dadV * E0 ;
+    J(03,02   ) = 2*a ;
+    J(03,12   ) = 2*dadC * E0 ;
+    J(03,14   ) = 2*dadR * E0 ;
+
+    % Partial derivatives of f4-10 = dE2-8/dt
+    J(31:15:121) = - dA - pA ;
+    J(32:15:122) = 2*pA ;
+    J(10,10) = - dA - m ;
+
+    % Partial derivatives of f11 = dX/dt
+    J(11,11) = -ka ; 
+
+    % Partial derivatives of f12 = dC/dt
+    J(12,11) = ka/vd ; 
+    J(12,12) = -ke ; 
+
+    % Partial derivatives of f13 = dR1/dt
+    J(13,01) = dsdV ;
+    J(13,12) = dsdC ;
+    J(13,13) = -dR ;
+
+    % Partial derivatives of f14 = dR2/dt
+    J(14,13) = dR ;
+    J(14,14) = -dR ;
+end
+
+function F = terms(Y)
+
+    V  = Y(:,01)            ;% SIV virions [#/無]
+    E  = sum(Y(:,02:10),2)  ;% CD8+ T cells [#/無]
+    C  = Y(:,12)            ;% N803 plasma concentration [pM]
+    R  = Y(:,14)            ;% immune regulation []
+
+    % percolate vector for storing terms
+    F = zeros(size(Y,1),14) ;
+
+    F(:,01) = C./ (C50+C)     ;% N803 effect []
+
+    F(:,02) = 1 + p1*F(:,1)  ;% E0 proliferation stimulation []
+    F(:,03) = 1 + a1*F(:,1)  ;% E0 activation stimulation []
+    F(:,04) =     s1*F(:,1)  ;% N803-induced R generation []
+
+    F(:,05) = 1./ (1+p2*R)   ;% E0 proliferation regulation []
+    F(:,06) = 1./ (1+a2*R)   ;% E0 activation regulation []
+    F(:,07) = 1./ (1+g2*R)   ;% EA killing regulation []
+
+    F(:,08) = E50./(E50+E)   ;% E0 proliferation density dependence []
+    F(:,09) = V./ (V50E+V)   ;% E0 activation antigen dependence []
+    F(:,10) = V./ (V50R+V)   ;% R  generation antigen dependence []
+
+    F(:,11) = p0* F(:,08).* F(:,02).* F(:,05)  ;% E0 proliferation [/d]
+    F(:,12) = a0* F(:,09).* F(:,03).* F(:,06)  ;% E0 activation rate [/d]
+    F(:,13) = s0* (F(:,10)+F(:,04))            ;% R generation rate [/d]
+    F(:,14) = g0* F(:,07)  ;% EA killing rate [#/無-d]
+
+end
+
+end