Skip to content

Commit

Permalink
Add files via upload
Browse files Browse the repository at this point in the history
  • Loading branch information
@cody2
cody2 authored Dec 24, 2022
1 parent 5791e28 commit efe6e5f
Showing 1 changed file with 39 additions and 25 deletions.
64 changes: 39 additions & 25 deletions N803_model_2.m
View file
Original file line number Diff line number Diff line change
@@ -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 |
Expand All @@ -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
%
%% ========================================================================
Expand All @@ -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'
Expand All @@ -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
Expand All @@ -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 -----------------------------------------------
Expand Down Expand Up @@ -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 ------------------------------------------

Expand Down Expand Up @@ -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

%% ------------------------------------------------------------------------
Expand All @@ -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) ] ;

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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 ;
Expand All @@ -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 ;
Expand Down

0 comments on commit efe6e5f

Please sign in to comment.