From fa59fece5bffd9b9a1f2e59c97fb7cf7f02f6752 Mon Sep 17 00:00:00 2001 From: "Cody, Jonathan William" Date: Sun, 8 Jan 2023 11:54:04 -0600 Subject: [PATCH] Delete N803_model_2S.m --- N803_model_2S.m | 353 ------------------------------------------------ 1 file changed, 353 deletions(-) delete mode 100644 N803_model_2S.m diff --git a/N803_model_2S.m b/N803_model_2S.m deleted file mode 100644 index 4c0e754..0000000 --- a/N803_model_2S.m +++ /dev/null @@ -1,353 +0,0 @@ -%% 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 \ No newline at end of file