diff --git a/N803_model_2.m b/N803_model_2.m index 07a2e7f..d88276d 100644 --- a/N803_model_2.m +++ b/N803_model_2.m @@ -1,17 +1,13 @@ %% N803_model_2.m - solves ODE model for N-803 treatment of SIV % -% Extra outputs are absolute instead of fold-change (otherwise identical) -% % /--------------------------------------------------------------\ -% | Date: 02/14/2022 | +% | Date: 06/26/2022 | % | Author: Jonathan Cody | % | Affiliation: Purdue University | % | Weldon School of Biomedical Engineering | % | Pienaar Computational Systems Pharmacology Lab | % \--------------------------------------------------------------/ % -% Version 2G with expanded outputs -% % Nomenclature: V = SIV virions [#/無] % T8 = total CD8+ T cells [#/無] % E0 = resting SIV-specific CD8+ T cells [#/無] @@ -21,7 +17,6 @@ % X = N803 at absorption site [pmol/kg] % C = N803 plasma concentration [pM] % R = regulation [] (dimensionless quantity) -% T = tolerance [] (dimensionless quantity) % %% ======================================================================== % INPUTS @@ -32,24 +27,16 @@ % 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) +% Pars = vector of parameters (see list on line 92) +% +% Yic = vector of initial conditions (see list on line 241) % %% ======================================================================== % 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) % -% 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 = 60) -% -% warnOff = switch to suppress display of errors +% varargin = cell vector used to add optional inputs (see line 54) +% EX: N803_model_2(SoluTimes,DoseTimes,Pars,Yic,'InputName',[InputValue]) % %% ======================================================================== % OUTPUTS @@ -57,151 +44,124 @@ % % Y_OUT(:,1) = V at points in 'SoluTimes' [log fold change] % Y_OUT(:,2) = T8 at points in 'SoluTimes' [fold change] -% Y_OUT(:,3) = R at points in 'SoluTimes' [fold change] -% -% If 'fullOut' is non-empty, more columns are included (see code) -% -% CalcPars = vector of calculated parameters % %% ======================================================================== % FUNCTION % ======================================================================== -function [Y_OUT,CalcPars] = N803_model_2(SoluTimes,DoseTimes,Pars,... - EquatedPars,fullOut,timeOut,warnOff) +function Y_OUT = N803_model_2(SoluTimes,DoseTimes,Pars,Yic,varargin) + +% 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(varargin) +if strcmp('EqualPars', varargin(n)) ; EqualPars = varargin{n+1} ; end +if strcmp('fullOut', varargin(n)) ; fullOut = varargin{n+1} ; end +if strcmp('timeOut', varargin(n)) ; timeOut = varargin{n+1} ; end +if strcmp('warnOff', varargin(n)) ; warnOff = varargin{n+1} ; end +if strcmp('AbsTol', varargin(n)) ; AbsTol = varargin{n+1} ; end +if strcmp('RelTol', varargin(n)) ; RelTol = varargin{n+1} ; end +if strcmp('odeSolver', varargin(n)) ; odeSolver = varargin{n+1} ; end +end -% convert inputs to column vectors +% check inputs vector orientation if size(SoluTimes ,1)==1 ; SoluTimes = SoluTimes' ; end if size(DoseTimes ,1)==1 ; DoseTimes = DoseTimes' ; end -if size(Pars ,1)==1 ; Pars = Pars' ; 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(EquatedPars) ; Pars = Pars(EquatedPars) ; end - -% set default evaluation time out to 60 seconds -if isempty(timeOut) ; timeOut = 60 ; end +if ~isempty(EqualPars) ; Pars = Pars(EqualPars) ; end % suppress warnings (if requested) -if ~isempty(warnOff) - warning off -end +if ~isempty(warnOff) ; warning off ; end %% ------------------------------------------------------------------------ % Rename inputed parameters ----------------------------------------------- -Vi = Pars(01) ;% V initial value [log(#/mL)] -EBi = Pars(02) ;% E+B initial value [#/無] -fE = Pars(03) ;% initial frequency: E/(E+B) -fA = Pars(04) ;% initial frequency: (Ea+Ba)/(E+B) -q = Pars(05) ;% V growth rate (if E+B were absent) [/d] -psi = Pars(06) ;% Ba/Ea killing rate ratio [gB0/gE0] -EB50 = Pars(07) ;% 50% E+B proliferation saturation [#/無] -V50E = Pars(08) ;% 50% viral stimulation saturation for E [#/mL] -V50B = Pars(09) ;% 50% viral stimulation saturation for B [#/mL] -pE = Pars(10) ;% Ea proliferation rate constant [/d] -pB = Pars(11) ;% Ba proliferation rate constant [/d] -d = Pars(12) ;% E0/B0 death rate constant [/d] -dA = Pars(13) ;% Ea/Ba death rate constant [/d] -mEn = Pars(14) ;% normalized Ea reversion rate constant [] -mBn = Pars(15) ;% normalized Ba reversion rate constant [] +q = Pars(01) ;% V growth rate (if E+B were absent) [/d] +gE0 = Pars(02) ;% Ea killing rate constant [/d] +gB0 = Pars(03) ;% Ba killing rate constant [/d] + +V50E = Pars(04) ;% 50% viral stimulation saturation for E [#/無] +V50B = Pars(05) ;% 50% viral stimulation saturation for B [#/無] +aE0 = Pars(06) ;% E activation rate constant [/d] +aB0 = Pars(07) ;% B activation rate constant [/d] +mE = Pars(08) ;% Ea reversion rate constant [/d] +mB = Pars(09) ;% Ba reversion rate constant [/d] + +EB50 = Pars(10) ;% 50% E+B proliferation saturation [#/無] +p0 = Pars(11) ;% E0/P0 proliferation rate constant [/d] +pE = Pars(12) ;% Ea proliferation rate constant [/d] +pB = Pars(13) ;% Ba proliferation rate constant [/d] +d = Pars(14) ;% E0/B0 death rate constant [/d] +dA = Pars(15) ;% Ea/Ba death rate constant [/d] + Xi = Pars(16) ;% X initial condition [pmol/kg] ka = Pars(17) ;% N803 absorption rate constant [/d] ke = Pars(18) ;% N803 elimination rate constant [/d] vd = Pars(19) ;% N803 'volume of distribution'/'bioavailability' [L/kg] C50 = Pars(20) ;% 50% N803 stimulation concentration [pM] -pm = Pars(21) ;% E0/B0 maximum proliferation rate [] +p1 = Pars(21) ;% E0/B0 proliferation stimulation factor aE1 = Pars(22) ;% E activation stimulation factor [] 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] -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 ---------------------- - -Vi = 10^(Vi-3) ;% V initial value [#/無] -V50E = V50E/1000 ;% 50% viral stimulation saturation for E [#/無] -V50B = V50B/1000 ;% 50% viral stimulation saturation for B [#/無] - -% restrict mE and mB such that phi (see below) is positive -kE = (2*pE/(dA+pE))^7 ; -kB = 2*pB/(dA+pB) ; -mE = mEn*dA/(kE-1) ;% Ea reversion rate constant [/d] -mB = mBn*dA/(kB-1) ;% Ba reversion rate constant [/d] - -% ratio between initial activation rate of E and of B -phi = (kB*mB/(mB+dA) - 1) / (kE*mE/(mE+dA) - 1) ; - -% solve for aE,aB initial values and E and B initial values -x = fzero(@(x)icFun(x),[-15 5]) ; -[~,aEi,aBi,Ei,Bi] = icFun(x) ; - -% regulation initial value -Ri = ( sE*aEi*(1+aE2R) + sB*aBi*(1+aB2R) ) / dR ; - -% Ea/Ba killing rate constants gE0/gB0 -gE0 = q / ( sum(Ei(2:9))/(1+gE2R) + psi*sum(Bi(2:3))/(1+gB2R) ) ; -gB0 = psi*gE0 ; -p0 = (d+aEi-mE*Ei(9)/Ei(1)) * (1+p2R) * (EB50+EBi)/EB50 ;% E0/B0 prolif -aE0 = aEi * (V50E+Vi)/Vi * (1+aE2R) ;% E activation rate constant [/d] -aB0 = aBi * (V50B+Vi)/Vi * (1+aB2R) ;% B activation rate constant [/d] -p1 = pm/p0 ;% E0/B0 proliferation stimulation factor -p2 = p2R/Ri ;% E0/B0 proliferation regulation factor -aE2 = aE2R/Ri ;% E activation regulation factor [] -aB2 = aB2R/Ri ;% B activation regulation factor [] -gE2 = gE2R/Ri ;% E killing regulation factor [] -gB2 = gB2R/Ri ;% B killing regulation factor [] - -% declare initial conditions for Y (see 'system' function) -% V E0-8 B0-2 X C R T -Yic = [ Vi Ei' Bi' 0 0 Ri*ones(1,nR) zeros(1,nT) ] ; +sE = Pars(24) ;% R generation due to E activation [/d] +sB = Pars(25) ;% R generation due to B activation [/d] +dR = Pars(26) ;% R decay rate constant [/d] +gE2 = Pars(27) ;% E killing regulation factor [] +gB2 = Pars(28) ;% B killing regulation factor [] +p2 = Pars(29) ;% E0/B0 proliferation regulation factor [] +aE2 = Pars(30) ;% E activation regulation factor [] +aB2 = Pars(31) ;% B 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 - [~,Y_TEMP] = ode15s(@system,SoluTimes(Piece),Yic) ;% 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 ode15s + 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(14) = Xi ; + Yic(14) = Yic(14) + Xi ;% add drug dose as initial condition end %% ------------------------------------------------------------------------ -% Prepare main outputs 'Y_OUT' ------------------------------------------- +% Prepare main outputs 'Y_OUT' -------------------------------------------- Y = abs(Y) ;% removing negatives resulting from numerical error -Y_OUT = [ log10(Y(:,1)/Vi) ... V [log fold change] - , sum(Y(:,2:13),2)/EBi ... E [fold change] - , Y(:,15+nR) ];% R - -% Calculated parameters -CalcPars = [mE mB gE0 gB0 p0 aE0 aB0 p1 p2 aE2 aB2 gE2 gB2 Ri] ; +Y_OUT = [ log10( Y(:,1) / Y(1,1) )... V [log fold change] + , sum(Y(:,2:13),2) / sum(Y(1,2:13)) ];% E+B [fold change] %% ------------------------------------------------------------------------ % Add columns to 'Y_OUT' (if requested) ----------------------------------- @@ -209,95 +169,67 @@ if isempty(fullOut) ; return ; end -Y_OUT = [ Y_OUT zeros(Pieces(end),32) ] ; +Y_OUT = [ Y_OUT, Y(:,2:end), zeros(Pieces(end),33) ] ; + +F = terms(Y) ;% see 'terms' function +gE = gE0*F(:,08) ;% specific killing rate [#/無-d] +gB = gB0*F(:,09) ;% bystander killing rate [#/無-d] +p = p0*F(:,13) ;% E0/B0 proliferation rate [/d] +aE = aE0*F(:,14) ;% specific activation rate [/d] +aB = aB0*F(:,15) ;% bystander activation rate [/d] + +Y_OUT(:,19) = F(:,01) ;% N803 effect [] E0 = Y(:,02) ;% resting specific CD8+ T cells [#/無] 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 [#/無] -Eap = sum( Y(:,03:09) ,2) ;% Ea that is proliferating -Bap = Y(:,12) ;% Ba that is proliferating +Ea = sum( Y(:,03:10) ,2) ;% active specific CD8+ T cells [#/無] +Ba = sum( Y(:,12:13) ,2) ;% active bystander T cells [#/無] +Y_OUT(:,20) = Ea ;% Ea +Y_OUT(:,21) = Ba ;% Ba +Y_OUT(:,22) = (E0+Ea) ;% E +Y_OUT(:,23) = (B0+Ba) ;% B +Y_OUT(:,24) = (E0+Ea)./(E0+B0+Ea+Ba) ;% E frequency in (E+B) +Y_OUT(:,25) = (Ea)./(E0+Ea) ;% activated frequency in E +Y_OUT(:,26) = (Ba)./(B0+Ba) ;% activated frequency in B + Eaf = Y(:,10) ;% last generation of Ea Baf = Y(:,13) ;% last generation of Ba -R = Y(:,15+nR) ;% immune regulation [] -[F,T] = terms(Y) ;% see 'terms' function -gE = gE0*F(:,10) ;% specific killing rate [#/無-d] -gB = gB0*F(:,11) ;% bystander killing rate [#/無-d] -p = p0*F(:,15) ;% E0/B0 proliferation rate [/d] -aE = aE0*F(:,16) ;% specific activation rate [/d] -aB = aB0*F(:,17) ;% bystander activation rate [/d] -RgenE = sE*aE0* F(:,13).* F(:,05) ;% R gen by aE -RgenB = sB*aB0*( F(:,14) + F(:,06) ) ;% R gen by aB - -Y_OUT(:,04) = Y(:,14) ;% X [pmol/kg] -Y_OUT(:,05) = Y(:,15) ;% C [pM] -Y_OUT(:,06) = T ;% T [] -Y_OUT(:,07) = F(:,01) ;% N803 effect [] -Y_OUT(:,08) = F(:,02) ;% tolerance effect [] -Y_OUT(:,09) = F(:,03) ;% N803 effect (with tolerance) [] - -Y_OUT(:,10) = (E0+Ea) ;% E -Y_OUT(:,11) = (B0+Ba) ;% B -Y_OUT(:,12) = (E0+Ea)./(E0+B0+Ea+Ba) ;% E frequency in (E+B) -Y_OUT(:,13) = (Ea+Ba)./(E0+B0+Ea+Ba) ;% activated frequency in (E+B) -Y_OUT(:,14) = (Ea)./(E0+Ea) ;% activated frequency in E -Y_OUT(:,15) = (Ba)./(B0+Ba) ;% activated frequency in B - -Y_OUT(:,16) = E0 ;% E0 -Y_OUT(:,17) = B0 ;% B0 -Y_OUT(:,18) = p - d ;% E0/B0 expansion ERC1 -Y_OUT(:,19) = mE*Eaf./E0 - aE ;% E0 activation ERC1 -Y_OUT(:,20) = mB*Baf./B0 - aB ;% B0 activation ERC1 -Y_OUT(:,21) = log10( F(:,04) ) ;% C multiplier for p* (log) -Y_OUT(:,22) = log10( F(:,07) ) ;% R multiplier for p* (log) -Y_OUT(:,23) = log10( F(:,12) ) ;% density-dep multiplier for p* (log) - -Y_OUT(:,24) = Ea ;% Ea -Y_OUT(:,25) = Ba ;% Ba -Y_OUT(:,26) = (pE*Eap - dA*Ea) ./ Ea ;% EA expansion ERC1 -Y_OUT(:,27) = (pB*Bap - dA*Ba) ./ Ba ;% BA expansion ERC1 -Y_OUT(:,28) = (aE.*E0 - mE*Eaf) ./ Ea ;% EA activation ERC1 -Y_OUT(:,29) = (aB.*B0 - mB*Baf) ./ Ba ;% BA activation ERC1 -Y_OUT(:,30) = log10( F(:,05) ) ;% C multiplier for aE* (log) -Y_OUT(:,31) = log10( ( F(:,06) + F(:,14) )./F(:,14) ) ;% same for aB* (log) -Y_OUT(:,32) = log10( F(:,08) ) ;% R multiplier for aE* (log) -Y_OUT(:,33) = log10( F(:,09) ) ;% R multiplier for aB* (log) -Y_OUT(:,34) = log10( F(:,13) ) ;% V multiplier for aE* (log) -Y_OUT(:,35) = log10( F(:,14) ) ;% V multiplier for aB* (log) - -Y_OUT(:,36) = q - gE.*Ea - gB.*Ba ;% V ERC1 -Y_OUT(:,37) = log10( gE.*Ea ) ;% E killing ERC1 (log) -Y_OUT(:,38) = log10( gB.*Ba ) ;% B killing ERC1 (log) -Y_OUT(:,39) = log10( F(:,10) ) ;% R multiplier for gE* (log) -Y_OUT(:,40) = log10( F(:,11) ) ;% R multiplier for gB* (log) -Y_OUT(:,41) = gB.*Ba ./ (gE.*Ea + gB.*Ba) ;% B killing fraction - -Y_OUT(:,42) = log10( RgenE ) ;% R gen by aE (log) -Y_OUT(:,43) = log10( RgenB ) ;% R gen by aB (log) -Y_OUT(:,44) = RgenB ./ (RgenE + RgenB) ;% R gen by aB fraction +Y_OUT(:,27) = p - d ;% E0/B0 expansion ERC1 +Y_OUT(:,28) = mE*Eaf./E0 - aE ;% E0 activation ERC1 +Y_OUT(:,29) = mB*Baf./B0 - aB ;% B0 activation ERC1 +Y_OUT(:,30) = log10( F(:,02) ) ;% C multiplier for p* (log) +Y_OUT(:,31) = log10( F(:,05) ) ;% R multiplier for p* (log) +Y_OUT(:,32) = log10( F(:,10) ) ;% density-dep multiplier for p* (log) + +Eap = sum( Y(:,03:09) ,2) ;% Ea that is proliferating +Bap = Y(:,12) ;% Ba that is proliferating +Y_OUT(:,33) = (pE*Eap - dA*Ea) ./ Ea ;% EA expansion ERC1 +Y_OUT(:,34) = (pB*Bap - dA*Ba) ./ Ba ;% BA expansion ERC1 +Y_OUT(:,35) = (aE.*E0 - mE*Eaf) ./ Ea ;% EA activation ERC1 +Y_OUT(:,36) = (aB.*B0 - mB*Baf) ./ Ba ;% BA activation ERC1 +Y_OUT(:,37) = log10( F(:,03) ) ;% C multiplier for aE* (log) +Y_OUT(:,38) = log10( ( F(:,04) + F(:,12) )./F(:,12) ) ;% same for aB* (log) +Y_OUT(:,39) = log10( F(:,06) ) ;% R multiplier for aE* (log) +Y_OUT(:,40) = log10( F(:,07) ) ;% R multiplier for aB* (log) +Y_OUT(:,41) = log10( F(:,11) ) ;% V multiplier for aE* (log) +Y_OUT(:,42) = log10( F(:,12) ) ;% V multiplier for aB* (log) + +Y_OUT(:,43) = q - gE.*Ea - gB.*Ba ;% V ERC1 +Y_OUT(:,44) = log10( gE.*Ea ) ;% E killing ERC1 (log) +Y_OUT(:,45) = log10( gB.*Ba ) ;% B killing ERC1 (log) +Y_OUT(:,46) = log10( F(:,08) ) ;% R multiplier for gE* (log) +Y_OUT(:,47) = log10( F(:,09) ) ;% R multiplier for gB* (log) +Y_OUT(:,48) = gB.*Ba ./ (gE.*Ea + gB.*Ba) ;% B killing fraction + +RgenE = sE*aE0* F(:,11).* F(:,03) ;% R gen by aE +RgenB = sB*aB0*( F(:,12) + F(:,04) ) ;% R gen by aB +Y_OUT(:,49) = log10( RgenE ) ;% R gen by aE (log) +Y_OUT(:,50) = log10( RgenB ) ;% R gen by aB (log) +Y_OUT(:,51) = RgenB ./ (RgenE + RgenB) ;% R gen by aB fraction %% ======================================================================== % NESTED FUNCTIONS % ======================================================================== -function [targ,aEi,aBi,Ei,Bi] = icFun(x) - - % declare initial activation rates aE and aB - aEi = 10^x ; - aBi = aEi/phi ; - - % solve for E0-E8 and B0-B2 initial conditions - MAT = zeros(12) ; - MAT( 1:13: 92) = [ aEi,2*pE*ones(1,7) ] ; - MAT( 13:13:130) = [ -(dA+pE)*ones(1,7),-(dA+mE),aBi,2*pB ] ; - MAT(129:13:142) = [ -(dA+pB),-(dA+mB) ] ; - MAT(11:12,:) = [ ones(1,12) ; [ ones(1,9) 0 0 0 ] ] ; - Sol = MAT \ [ zeros(10,1) ; EBi ; fE*EBi ] ; - Ei = Sol(1:9) ; - Bi = Sol(10:12) ; - - % ensure initial conditions satisfy fA - targ = Ei(1) + Bi(1) - (1-fA)*EBi ; -end function dY = system(~,Y) @@ -307,85 +239,163 @@ V = Y(01) ;% SIV virions [#/無] E0 = Y(02) ;% resting specific CD8+ T cells [#/無] - Ea = Y(03:10) ;% active specific CD8+ T cells [#/無] + Ea = Y(03:10) ;% active specific CD8+ T cells [#/無] (E1 to E8) B0 = Y(11) ;% resting bystander CD8+ T cells [#/無] - Ba = Y(12:13) ;% active bystander T cells [#/無] + Ba = Y(12:13) ;% active bystander T cells [#/無] (B1 to B2) X = Y(14) ;% N803 at absorption site [pmol/kg] C = Y(15) ;% N803 plasma concentration [pM] + R = Y(16:17) ;% immune regulation [] F = terms(Y') ;% modifiers for 'gE0','gB0','p0','aE0','aB0' - gE = gE0*F(10) ;% specific killing rate [#/無-d] - gB = gB0*F(11) ;% bystander killing rate [#/無-d] - p = p0*F(15) ;% E0/B0 proliferation rate [/d] - aE = aE0*F(16) ;% specific activation rate [/d] - aB = aB0*F(17) ;% bystander activation rate [/d] + gE = gE0*F(08) ;% specific killing rate [#/無-d] + gB = gB0*F(09) ;% bystander killing rate [#/無-d] + p = p0*F(13) ;% E0/B0 proliferation rate [/d] + aE = aE0*F(14) ;% specific activation rate [/d] + aB = aB0*F(15) ;% bystander activation rate [/d] dY = Y ;% dY/dt % 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(14) = - ka*X ;% dX/dt - dY(15) = ka*X/vd - ke*C ;% dC/dt + 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(14) = - ka*X ;% dX/dt + dY(15) = ka*X/vd - ke*C ;% dC/dt + dY(16) = F(16) - dR*R(1) ;% dR1/dt + dY(17) = dR*R(1) - dR*R(2) ;% dR2/dt + +end + +function J = jacob(~,Y) - % 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 + V = Y(01) ;% SIV virions [#/無] + E0 = Y(02) ;% resting specific CD8+ T cells [#/無] + Ea = sum(Y(03:10)) ;% active specific CD8+ T cells [#/無] + B0 = Y(11) ;% resting bystander CD8+ T cells [#/無] + Ba = Y(12)+Y(13) ;% active bystander T cells [#/無] + C = Y(15) ;% N803 plasma concentration [pM] + F = terms(Y') ;% modifiers for 'gE0','gB0','p0','aE0','aB0' - % 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 + gE = gE0*F(08) ;% specific killing rate [#/無-d] + gB = gB0*F(09) ;% bystander killing rate [#/無-d] + p = p0*F(13) ;% E0/B0 proliferation rate [/d] + aE = aE0*F(14) ;% specific activation rate [/d] + aB = aB0*F(15) ;% bystander activation rate [/d] + + J = zeros(17) ; + + % Partial derivatives of gE,gB,p,aE,aB w.r.t. V,E,B,C,R + dpdR = - p * p2 * F(05) ; + daEdR = - aE * aE2 * F(06) ; + daBdR = - aB * aB2 * F(07) ; + dgEdR = - gE * gE2 * F(08) ; + dgBdR = - gB * gB2 * F(09) ; + dpdEB = - p * F(10) / EB50 ; + daEdV = aE0 * V50E / (V50E+V)^2 * F(03) * F(06) ; + daBdV = aB0 * V50B / (V50B+V)^2 * F(07) ; + dOdC = C50 / (C50+C)^2 ; + dpdC = p0 * F(10) * p1 * dOdC * F(05) ; + daEdC = aE0 * F(11) * aE1 * dOdC * F(06) ; + daBdC = aB0 * aB1 * dOdC * F(07) ; + + % Partial derivatives of f1 = dV/dt + J(01,01 ) = q - gE*Ea - gB*Ba ; + J(01,02:10) = - gE * V ; + J(01,11:13) = - gB * V ; + J(01,17 ) = - ( dgEdR*Ea + dgBdR*Ba ) * V ; + + % Partial derivatives of f2 = dE0/dt + J(02,01 ) = - daEdV * E0 ; + J(02,02 ) = p + dpdEB*E0 - d - aE ; + J(02,03:13) = dpdEB*E0 ; + J(02,10 ) = dpdEB*E0 + mE ; + J(02,15 ) = ( dpdC - daEdC ) * E0 ; + 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 ; + + % Partial derivatives of f4-10 = dE2-8/dt + J(39:19:153) = - dA - pE ; + J(40:19:154) = 2*pE ; + J(10,10) = - dA - mE ; + + % Partial derivatives of f11 = dB0/dt + J(11,01 ) = - daBdV*B0 ; + J(11,02:12) = dpdEB*B0 ; + J(11,11 ) = p + dpdEB*B0 - d - aB ; + J(11,13 ) = dpdEB*B0 + mB ; + J(11,15 ) = ( dpdC - daBdC ) * B0 ; + J(11,17 ) = ( dpdR - daBdR ) * B0 ; + + % Partial derivatives of f12 = dB1/dt + J(12,01) = daBdV * B0 ; + J(12,11) = aB ; + J(12,12) = - dA - pB ; + J(12,15) = daBdC * B0 ; + J(12,17) = daBdR * B0 ; + + % Partial derivatives of f12-13 = dB1-2/dt + J(13,12) = 2*pB ; + J(13,13) = - dA - mB ; + + % Partial derivatives of f14 = dX/dt + J(14,14) = -ka ; + + % Partial derivatives of f15 = dC/dt + J(15,14) = ka/vd ; + J(15,15) = -ke ; + + % Partial derivatives of f16 = dR1/dt + J(16,01) = sE*daEdV/F(06) + sB*daBdV/F(07) ; + J(16,15) = sE*daEdC/F(06) + sB*daBdC/F(07) ; + J(16,16) = -dR ; + + % Partial derivatives of f17 = dR2/dt + J(17,16) = dR ; + J(17,17) = -dR ; end -function [F,T] = terms(Y) +function F = 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(:,15+nR) ;% immune regulation [] - T = (Y(:,14+nR+nT)+Y(:,15+nR+nT)) ;% drug tolerance [] + V = Y(:,01) ;% SIV virions [#/無] + EB = sum(Y(:,02:13),2) ;% CD8+ T cells [#/無] + C = Y(:,15) ;% N803 plasma concentration [pM] + R = Y(:,17) ;% immune regulation [] % percolate vector for storing terms - F = zeros(size(Y,1),18) ; + F = zeros(size(Y,1),16) ; F(:,01) = C./ (C50+C) ;% N803 effect [] - F(:,02) = 1./ (1+eta*T) ;% tolerance effect [] - F(:,03) = F(:,1).* F(:,2) ;% N803 effect (with tolerance) [] - F(:,04) = 1 + p1*F(:,3) ;% E0/B0 proliferation stimulation [] - F(:,05) = 1 + aE1*F(:,3) ;% E activation stimulation [] - F(:,06) = aB1*F(:,3) ;% N803-induced B activation [] + F(:,02) = 1 + p1*F(:,1) ;% E0/B0 proliferation stimulation [] + F(:,03) = 1 + aE1*F(:,1) ;% E activation stimulation [] + F(:,04) = aB1*F(:,1) ;% N803-induced B activation [] - F(:,07) = 1./ (1+ p2*R) ;% E0/B0 proliferation regulation [] - F(:,08) = 1./ (1+aE2*R) ;% Ea activation regulation [] - F(:,09) = 1./ (1+aB2*R) ;% Ba activation regulation [] - F(:,10) = 1./ (1+gE2*R) ;% Ea killing regulation [] - F(:,11) = 1./ (1+gB2*R) ;% Ba killing regulation [] + F(:,05) = 1./ (1+ p2*R) ;% E0/B0 proliferation regulation [] + F(:,06) = 1./ (1+aE2*R) ;% Ea activation regulation [] + F(:,07) = 1./ (1+aB2*R) ;% Ba activation regulation [] + F(:,08) = 1./ (1+gE2*R) ;% Ea killing regulation [] + F(:,09) = 1./ (1+gB2*R) ;% Ba killing regulation [] - F(:,12) = EB50./(EB50+EB) ;% E0/B0 proliferation density dependence [] - F(:,13) = V./ (V50E+V) ;% Ea activation antigen dependence [] - F(:,14) = V./ (V50B+V) ;% Ba activation antigen dependence [] + F(:,10) = EB50./(EB50+EB) ;% E0/B0 proliferation density dependence [] + F(:,11) = V./ (V50E+V) ;% Ea activation antigen dependence [] + F(:,12) = V./ (V50B+V) ;% Ba activation antigen dependence [] - F(:,15) = F(:,12).* F(:,04).* F(:,07) ;% E0/B0 prolif modifier [] - F(:,16) = F(:,13).* F(:,05).* F(:,08) ;% Ea activation modifier [] - F(:,17) = (F(:,14)+F(:,06)).* F(:,09) ;% Ba activation modifier [] + F(:,13) = F(:,10).* F(:,02).* F(:,05) ;% E0/B0 prolif modifier [] + F(:,14) = F(:,11).* F(:,03).* F(:,06) ;% Ea activation modifier [] + F(:,15) = (F(:,12)+F(:,04)).* F(:,07) ;% Ba activation modifier [] - F(:,18) = sE*aE0* F(:,13).* F(:,05) ... - + sB*aB0*( F(:,14) + F(:,06) ) ;% regulation gen [/d] + F(:,16) = sE*aE0* F(:,11).* F(:,03) ... + + sB*aB0*( F(:,12) + F(:,04) ) ;% regulation gen [/d] end end \ No newline at end of file