diff --git a/N803_model_2.m b/N803_model_2.m index a8e5255..6509329 100644 --- a/N803_model_2.m +++ b/N803_model_2.m @@ -1,7 +1,7 @@ %% N803_model_2.m - solves ODE model for N-803 treatment of SIV % % /--------------------------------------------------------------\ -% | Date: 07/08/2022 | +% | Date: 08/14/2022 | % | Author: Jonathan Cody | % | Affiliation: Purdue University | % | Weldon School of Biomedical Engineering | @@ -27,15 +27,15 @@ % 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 92) +% Pars = vector of parameters (see list on line 96) % -% Yic = vector of initial conditions (see list on line 241) +% Yic = vector of initial conditions (see list on line 273) % %% ======================================================================== % OPTIONS % ======================================================================== % -% Opts = cell vector used to add optional inputs (see list on line 54) +% Opts = cell vector used to add optional inputs (see list on line 53) % EX: Opts = {'AbsTol',1e-2} will set ode solver absolute tolerance to 1e-2 % %% ======================================================================== @@ -53,10 +53,12 @@ % 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) +fullOut = false ;% if true, more outputs are added (see line 183) +ivDosing = [] ;% vector of intravenous dose strengths +isNaive = false ;% if true, all IC are zero (except for B0) timeOut = 60 ;% scalar time limit [seconds] for ODE solver % if exceeded, an error is thrown -warnOff = [] ;% if non-empty, warning messages are suppressed +warnOff = false ;% if true, 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' @@ -65,6 +67,8 @@ 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('ivDosing', Opts(n)) ; ivDosing = Opts{n+1} ; end +if strcmp('isNaive', Opts(n)) ; isNaive = 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 @@ -84,7 +88,7 @@ if ~isempty(EqualPars) ; Pars = Pars(EqualPars) ; end % suppress warnings (if requested) -if ~isempty(warnOff) ; warning off ; end +if warnOff ; warning off ; end %% ------------------------------------------------------------------------ % Rename inputed parameters ----------------------------------------------- @@ -125,6 +129,11 @@ aE2 = Pars(30) ;% E activation regulation factor [] aB2 = Pars(31) ;% B activation regulation factor [] +if isNaive % If SIV-naive, calculate initial B0 + Yic = Yic*0 ; + Yic(11) = EB50*( p0/d - 1 ) ; +end + %% ------------------------------------------------------------------------ % Solve for each solution period ------------------------------------------ @@ -152,7 +161,12 @@ Y(Piece,:) = Y_TEMP ; % add 'Y_TEMP' to full solution 'Y' Yic = Y_TEMP(end,:) ;% update ICs for next solution piece - Yic(14) = Yic(14) + Xi ;% add drug dose as initial condition + + if isempty(ivDosing) + Yic(14) = Yic(14) + Xi ;% add SC drug dose as initial condition + else + Yic(15) = Yic(15) + ivDosing(n)/vd ;% add IV drug dose + end end %% ------------------------------------------------------------------------ @@ -166,7 +180,7 @@ %% ------------------------------------------------------------------------ % Add columns to 'Y_OUT' (if requested) ----------------------------------- -if isempty(fullOut) ; return ; end +if ~fullOut ; return ; end Y_OUT = [ Y_OUT, Y, zeros(Pieces(end),46) ] ; @@ -194,7 +208,7 @@ Y_OUT(:,23) = Ba ;% Ba Y_OUT(:,24) = (E0+Ea) ;% E Y_OUT(:,25) = (B0+Ba) ;% B -Y_OUT(:,26) = (E0+Ea)./(E0+B0+Ea+Ba) ;% E frequency in (E+B) +Y_OUT(:,26) = (B0+Ba)./(E0+B0+Ea+Ba) ;% B frequency in (E+B) Y_OUT(:,27) = (Ea)./(E0+Ea) ;% activated frequency in E Y_OUT(:,28) = (Ba)./(B0+Ba) ;% activated frequency in B Y_OUT(:,29) = log10( Ea ) ;% Ea (log) @@ -277,13 +291,13 @@ % calculate rates of change for V,E,B,X,C,R dY(01) = q*V - gE*sum(Ea)*V - gB*sum(Ba)*V ;% dV/dt - dY(02) = p*E0 - aE*E0 - d*E0 + mE*Ea(8) ;% dE0/dt - dY(03) = + aE*E0 - dA*Ea(1) - pE*Ea(1) ;% dE1/dt - dY(4:9)= 2*pE*Ea(1:6) - dA*Ea(2:7) - pE*Ea(2:7) ;% dE2-7/dt - dY(10) = 2*pE*Ea(7) - dA*Ea(8) - mE*Ea(8) ;% dE8/dt - dY(11) = p*B0 - aB*B0 - d*B0 + mB*Ba(2) ;% dB0/dt - dY(12) = + aB*B0 - dA*Ba(1) - pB*Ba(1) ;% dB1/dt - dY(13) = 2*pB*Ba(1) - dA*Ba(2) - mB*Ba(2) ;% dB2/dt + dY(02) = p*E0 - aE*E0 - d*E0 + mE*Ea(8) ;% dE0/dt + dY(03) = + 2*aE*E0 - dA*Ea(1) - pE*Ea(1) ;% dE1/dt + dY(4:9)= 2*pE*Ea(1:6) - dA*Ea(2:7) - pE*Ea(2:7) ;% dE2-7/dt + dY(10) = 2*pE*Ea(7) - dA*Ea(8) - mE*Ea(8) ;% dE8/dt + dY(11) = p*B0 - aB*B0 - d*B0 + mB*Ba(2) ;% dB0/dt + dY(12) = + 2*aB*B0 - dA*Ba(1) - pB*Ba(1) ;% dB1/dt + 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) = s - dR*R(1) ;% dR1/dt @@ -338,10 +352,10 @@ J(02,17 ) = ( dpdR - daEdR ) * E0 ; % Partial derivatives of f3 = dE1/dt (save f3/dE1) - J(03,01 ) = daEdV * E0 ; - J(03,02 ) = aE ; - J(03,15 ) = daEdC * E0 ; - J(03,17 ) = daEdR * E0 ; + J(03,01 ) = 2*daEdV * E0 ; + J(03,02 ) = 2*aE ; + J(03,15 ) = 2*daEdC * E0 ; + J(03,17 ) = 2*daEdR * E0 ; % Partial derivatives of f4-10 = dE2-8/dt J(37:18:145) = - dA - pE ; @@ -357,11 +371,11 @@ J(11,17 ) = ( dpdR - daBdR ) * B0 ; % Partial derivatives of f12 = dB1/dt - J(12,01) = daBdV * B0 ; - J(12,11) = aB ; + J(12,01) = 2*daBdV * B0 ; + J(12,11) = 2*aB ; J(12,12) = - dA - pB ; - J(12,15) = daBdC * B0 ; - J(12,17) = daBdR * B0 ; + J(12,15) = 2*daBdC * B0 ; + J(12,17) = 2*daBdR * B0 ; % Partial derivatives of f12-13 = dB1-2/dt J(13,12) = 2*pB ;