diff --git a/N803_model_2.m b/N803_model_2.m index 08819e3..07a2e7f 100644 --- a/N803_model_2.m +++ b/N803_model_2.m @@ -1,19 +1,19 @@ %% N803_model_2.m - solves ODE model for N-803 treatment of SIV % +% Extra outputs are absolute instead of fold-change (otherwise identical) +% % /--------------------------------------------------------------\ -% | Date: 11/05/2021 | +% | Date: 02/14/2022 | % | Author: Jonathan Cody | % | Affiliation: Purdue University | % | Weldon School of Biomedical Engineering | % | Pienaar Computational Systems Pharmacology Lab | % \--------------------------------------------------------------/ % -% Version 2D+ (also allowing regulation calibration) -% -% 12/03/2021 - added some extra terms to the full output matrix -% 11/27/2021 - replaced mu with fA -> initial frequency: (Ea+Ba)/(E+B) +% Version 2G with expanded outputs % % Nomenclature: V = SIV virions [#/無] +% T8 = total CD8+ T cells [#/無] % E0 = resting SIV-specific CD8+ T cells [#/無] % Ea = active SIV-specific CD8+ T cells [#/無] % B0 = resting bystander CD8+ T cells [#/無] @@ -56,16 +56,18 @@ % ======================================================================== % % Y_OUT(:,1) = V at points in 'SoluTimes' [log fold change] -% Y_OUT(:,2) = E at points in 'SoluTimes' [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 = N803_model_2(SoluTimes,DoseTimes,Pars,... - EquatedPars,fullOut,timeOut,warnOff) +function [Y_OUT,CalcPars] = N803_model_2(SoluTimes,DoseTimes,Pars,... + EquatedPars,fullOut,timeOut,warnOff) % convert inputs to column vectors if size(SoluTimes ,1)==1 ; SoluTimes = SoluTimes' ; end @@ -196,76 +198,83 @@ Y_OUT = [ log10(Y(:,1)/Vi) ... V [log fold change] , sum(Y(:,2:13),2)/EBi ... E [fold change] - , Y(:,15+nR)/Y(1,15+nR) ];% R [fold change] + , Y(:,15+nR) ];% R + +% Calculated parameters +CalcPars = [mE mB gE0 gB0 p0 aE0 aB0 p1 p2 aE2 aB2 gE2 gB2 Ri] ; %% ------------------------------------------------------------------------ % Add columns to 'Y_OUT' (if requested) ----------------------------------- +% NOTE: ERC1 = effective first-order growth rate constant if isempty(fullOut) ; return ; end Y_OUT = [ Y_OUT zeros(Pieces(end),32) ] ; -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 [#/無] -R = Y(:,15+nR) ;% immune regulation [] -[F,T] = terms(Y) ;% modifiers for 'gE0','gB0','p0','aE0','aB0' - -Y_OUT(:,04) = (E0+Ea) / (E0(1)+Ea(1)) ;% E [fold change] -Y_OUT(:,05) = (B0+Ba) / (B0(1)+Ba(1)) ;% B [fold change] -Y_OUT(:,06) = E0 / E0(1) ;% E0 [fold change] -Y_OUT(:,07) = B0 / B0(1) ;% B0 [fold change] -Y_OUT(:,08) = Ea / Ea(1) ;% Ea [fold change] -Y_OUT(:,09) = Ba / Ba(1) ;% Ba [fold change] -Y_OUT(:,10) = (E0+Ea)./(E0+B0+Ea+Ba) ;% E frequency in (E+B) -Y_OUT(:,11) = (Ea+Ba)./(E0+B0+Ea+Ba) ;% activated frequency in (E+B) -Y_OUT(:,12) = (Ea)./(E0+Ea) ;% activated frequency in E -Y_OUT(:,13) = (Ba)./(B0+Ba) ;% activated frequency in B -Y_OUT(:,14) = Y(:,14) ;% X [pmol/kg] -Y_OUT(:,15) = Y(:,15) ;% C [pM] -Y_OUT(:,16) = T ;% T [] - -Y_OUT(:,17) = F(:,01) ;% N803 effect [] -Y_OUT(:,18) = F(:,02) ;% tolerance effect [] -Y_OUT(:,19) = F(:,03) ;% N803 effect (with tolerance) [] -Y_OUT(:,20) = log10( F(:,04) ) ;% C-lfc in p -Y_OUT(:,21) = log10( F(:,05) ) ;% C-lfc in aE -Y_OUT(:,22) = log10( ( F(:,06) + F(:,14) )./F(:,14) ) ;% C-lfc in aB -Y_OUT(:,23) = log10( F(:,07)/F(1,07) ) ;% R-lfc in p -Y_OUT(:,24) = log10( F(:,08)/F(1,08) ) ;% R-lfc in aE -Y_OUT(:,25) = log10( F(:,09)/F(1,09) ) ;% R-lfc in aB -Y_OUT(:,26) = log10( F(:,10)/F(1,10) ) ;% R-lfc in gE -Y_OUT(:,27) = log10( F(:,11)/F(1,11) ) ;% R-lfc in gB -Y_OUT(:,28) = log10( F(:,12)/F(1,12) ) ;% EB-lfc in p -Y_OUT(:,29) = log10( F(:,13)/F(1,13) ) ;% V-lfc in aE -Y_OUT(:,30) = log10( F(:,14)/F(1,14) ) ;% V-lfc in aB -Y_OUT(:,31) = log10( F(:,15)/F(1,15) ) ;% lfc in p -Y_OUT(:,32) = log10( F(:,16)/F(1,16) ) ;% lfc in aE -Y_OUT(:,33) = log10( F(:,17)/F(1,17) ) ;% lfc in aB -Y_OUT(:,34) = log10( F(:,18)/F(1,18) ) ;% lfc in regulation generation - -% Ba regulation stimulation proportion -Y_OUT(:,35) = sB*aB0*(F(:,14)+F(:,06)) ./ F(:,18) ; - -% lfc in killing -Kill = gB0*Ba./(1+gB2*R) + gE0*Ea./(1+gE2*R) ; -Y_OUT(:,36) = log10( Kill / Kill(1) ) ; - -% Ba killing proportion -Y_OUT(:,37) = gB0*Ba./(1+gB2*R) ./ Kill ; - -% Average E proliferation rate -Y_OUT(:,38) = ( p0*F(:,15).*E0 + pE*sum( Y(:,03:9) ,2) ) ./ (E0+Ea) ; - -% Average B proliferation rate -Y_OUT(:,39) = ( p0*F(:,15).*B0 + pB*Y(:,12) ) ./ (B0+Ba) ; - -% Average E killing rate -Y_OUT(:,40) = gE0*Ea./(1+gE2*R) ./ (E0+Ea) ; - -% Average B killing rate -Y_OUT(:,41) = gB0*Ba./(1+gB2*R) ./ (B0+Ba) ; +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 +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 %% ======================================================================== % NESTED FUNCTIONS @@ -305,8 +314,8 @@ C = Y(15) ;% N803 plasma concentration [pM] 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] + 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]