From e1d5fa75d29ac71299ebf90d0b567983886fe8a6 Mon Sep 17 00:00:00 2001 From: "Cody, Jonathan William" Date: Sun, 25 Dec 2022 08:38:38 -0600 Subject: [PATCH] Add files via upload --- N803_model_2.m | 372 ++++++++++++++++++++++++------------------------- N803_shared.m | 255 +++++++++++++++++++++++++++++++++ N803_single.m | 217 +++++++++++++++++++++++++++++ 3 files changed, 658 insertions(+), 186 deletions(-) create mode 100644 N803_shared.m create mode 100644 N803_single.m diff --git a/N803_model_2.m b/N803_model_2.m index 6509329..28d06b2 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: 08/14/2022 | +% | Date: 12/25/2022 | % | Author: Jonathan Cody | % | Affiliation: Purdue University | % | Weldon School of Biomedical Engineering | @@ -10,10 +10,10 @@ % % 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 [#/無] -% Ba = active bystander CD8+ T cells [#/無] +% S0 = resting SIV-specific CD8+ T cells [#/無] +% Sa = active SIV-specific CD8+ T cells [#/無] +% N0 = resting non-SIV-specific CD8+ T cells [#/無] +% Na = active non-SIV-specific CD8+ T cells [#/無] % X = N803 at absorption site [pmol/kg] % C = N803 plasma concentration [pM] % R = regulation [] (dimensionless quantity) @@ -55,7 +55,7 @@ % EX: EqualPars(28) = 27 will set P(28) = P(27) 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) +isNaive = false ;% if true, all IC are zero (except for N0) timeOut = 60 ;% scalar time limit [seconds] for ODE solver % if exceeded, an error is thrown warnOff = false ;% if true, warning messages are suppressed @@ -93,45 +93,45 @@ %% ------------------------------------------------------------------------ % Rename inputed parameters ----------------------------------------------- -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] +q = Pars(01) ;% V growth rate (if S+N were absent) [/d] +gS0 = Pars(02) ;% Sa killing rate constant [無/#-d] +gN0 = Pars(03) ;% Na 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) ;% E0 activation rate constant [/d] -aB0 = Pars(07) ;% B0 activation rate constant [/d] -mE = Pars(08) ;% Ea reversion rate constant [/d] -mB = Pars(09) ;% Ba reversion rate constant [/d] +V50S = Pars(04) ;% 50% viral stimulation saturation for S [#/無] +V50N = Pars(05) ;% 50% viral stimulation saturation for N [#/無] +aS0 = Pars(06) ;% S0 activation rate constant [/d] +aN0 = Pars(07) ;% N0 activation rate constant [/d] +mS = Pars(08) ;% Sa reversion rate constant [/d] +mN = Pars(09) ;% Na reversion rate constant [/d] -EB50 = Pars(10) ;% 50% E+B proliferation saturation [#/無] -p0 = Pars(11) ;% E0/B0 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] +SN50 = Pars(10) ;% 50% S+N proliferation saturation [#/無] +p0 = Pars(11) ;% S0/N0 proliferation rate constant [/d] +pS = Pars(12) ;% Sa proliferation rate constant [/d] +pN = Pars(13) ;% Na proliferation rate constant [/d] +d = Pars(14) ;% S0/N0 death rate constant [/d] +dA = Pars(15) ;% Sa/Na 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] -p1 = Pars(21) ;% E0/B0 proliferation stimulation factor [] -aE1 = Pars(22) ;% E activation stimulation factor [] -aB1 = Pars(23) ;% B activation stimulation factor [] +p1 = Pars(21) ;% S0/N0 proliferation stimulation factor [] +aS1 = Pars(22) ;% S activation stimulation factor [] +aN1 = Pars(23) ;% N activation stimulation factor [] -sE = Pars(24) ;% R generation due to E0 activation [/d] -sB = Pars(25) ;% R generation due to B0 activation [/d] +sS = Pars(24) ;% R generation due to S0 activation [/d] +sN = Pars(25) ;% R generation due to N0 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 [] +gS2 = Pars(27) ;% S killing regulation factor [] +gN2 = Pars(28) ;% N killing regulation factor [] +p2 = Pars(29) ;% S0/N0 proliferation regulation factor [] +aS2 = Pars(30) ;% S activation regulation factor [] +aN2 = Pars(31) ;% N activation regulation factor [] -if isNaive % If SIV-naive, calculate initial B0 +if isNaive % If SIV-naive, calculate initial N0 Yic = Yic*0 ; - Yic(11) = EB50*( p0/d - 1 ) ; + Yic(11) = SN50*( p0/d - 1 ) ; end %% ------------------------------------------------------------------------ @@ -175,7 +175,7 @@ Y = abs(Y) ;% removing negatives resulting from numerical error 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] + , sum(Y(:,2:13),2) / sum(Y(1,2:13)) ];% S+N [fold change] %% ------------------------------------------------------------------------ % Add columns to 'Y_OUT' (if requested) ----------------------------------- @@ -185,73 +185,73 @@ Y_OUT = [ Y_OUT, Y, zeros(Pieces(end),46) ] ; V = Y(:,01) ;% SIV virions [#/無] -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 [#/無] -Eaf = Y(:,10) ;% last generation of Ea -Baf = Y(:,13) ;% last generation of Ba -Eap = Ea - Eaf ;% Ea that is proliferating -Bap = Y(:,12) ;% Ba that is proliferating +S0 = Y(:,02) ;% resting specific CD8+ T cells [#/無] +N0 = Y(:,11) ;% resting non-SIV-specific CD8+ T cells [#/無] +Sa = sum( Y(:,03:10) ,2) ;% active specific CD8+ T cells [#/無] +Na = sum( Y(:,12:13) ,2) ;% active non-SIV-specific T cells [#/無] +Saf = Y(:,10) ;% last generation of Sa +Naf = Y(:,13) ;% last generation of Na +Sap = Sa - Saf ;% Sa that is proliferating +Nap = Y(:,12) ;% Na that is proliferating F = terms(Y) ;% see 'terms' function -p = F(:,13) ;% E0/B0 proliferation rate [/d] -aE = F(:,14) ;% E0 activation rate [/d] -aB = F(:,15) ;% B0 activation rate [/d] -gE = F(:,16) ;% Ea killing rate [#/無-d] -gB = F(:,17) ;% Ba killing rate [#/無-d] +p = F(:,13) ;% S0/N0 proliferation rate [/d] +aS = F(:,14) ;% S0 activation rate [/d] +aN = F(:,15) ;% N0 activation rate [/d] +gS = F(:,16) ;% Sa killing rate [#/無-d] +gN = F(:,17) ;% Na killing rate [#/無-d] Y_OUT(:,20) = F(:,01) ;% N803 effect [] -Y_OUT(:,21) = (E0+B0+Ea+Ba) ;% CD8+ T cells [#/無] -Y_OUT(:,22) = Ea ;% Ea -Y_OUT(:,23) = Ba ;% Ba -Y_OUT(:,24) = (E0+Ea) ;% E -Y_OUT(:,25) = (B0+Ba) ;% 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) -Y_OUT(:,30) = log10( Ba ) ;% Ba (log) - -Y_OUT(:,31) = p.*E0 ;% E0 proliferation term [#/無/d] -Y_OUT(:,32) = p.*B0 ;% B0 proliferation term [#/無/d] -Y_OUT(:,33) = d.*E0 ;% E0 death term [#/無/d] -Y_OUT(:,34) = d.*B0 ;% B0 death term [#/無/d] -Y_OUT(:,35) = p ;% E0/B0 proliferation rate [/d] +Y_OUT(:,21) = (S0+N0+Sa+Na) ;% CD8+ T cells [#/無] +Y_OUT(:,22) = Sa ;% Sa +Y_OUT(:,23) = Na ;% Na +Y_OUT(:,24) = (S0+Sa) ;% S +Y_OUT(:,25) = (N0+Na) ;% N +Y_OUT(:,26) = (N0+Na)./(S0+N0+Sa+Na) ;% N frequency in (S+N) +Y_OUT(:,27) = (Sa)./(S0+Sa) ;% activated frequency in S +Y_OUT(:,28) = (Na)./(N0+Na) ;% activated frequency in N +Y_OUT(:,29) = log10( Sa ) ;% Sa (log) +Y_OUT(:,30) = log10( Na ) ;% Na (log) + +Y_OUT(:,31) = p.*S0 ;% S0 proliferation term [#/無/d] +Y_OUT(:,32) = p.*N0 ;% N0 proliferation term [#/無/d] +Y_OUT(:,33) = d.*S0 ;% S0 death term [#/無/d] +Y_OUT(:,34) = d.*N0 ;% N0 death term [#/無/d] +Y_OUT(:,35) = p ;% S0/N0 proliferation rate [/d] Y_OUT(:,36) = F(:,02) ;% C multiplier for p* (log) Y_OUT(:,37) = F(:,05) ;% R multiplier for p* (log) Y_OUT(:,38) = F(:,10) ;% density-dep multiplier for p* (log) -Y_OUT(:,39) = aE.*E0 ;% E0 activation term [#/無/d] -Y_OUT(:,40) = aB.*B0 ;% B0 activation term [#/無/d] -Y_OUT(:,41) = pE*Eap ;% Ea proliferation term [#/無/d] -Y_OUT(:,42) = pB*Bap ;% Ba proliferation term [#/無/d] -Y_OUT(:,43) = dA*Ea ;% Ea death term [#/無/d] -Y_OUT(:,44) = dA*Ba ;% Ba death term [#/無/d] -Y_OUT(:,45) = mE*Eaf ;% Ea reversion term [#/無/d] -Y_OUT(:,46) = mB*Baf ;% Ba reversion term [#/無/d] -Y_OUT(:,47) = aE ;% E0 activation rate [/d] -Y_OUT(:,48) = aB ;% B0 activation rate [/d] -Y_OUT(:,49) = F(:,03) ;% C multiplier for aE* -Y_OUT(:,50) = ( F(:,04) + F(:,12) )./F(:,12) ;% C multiplier for aE* -Y_OUT(:,51) = F(:,06) ;% R multiplier for aE* -Y_OUT(:,52) = F(:,07) ;% R multiplier for aB* -Y_OUT(:,53) = F(:,11) ;% V multiplier for aE* -Y_OUT(:,54) = F(:,12) ;% V multiplier for aB* +Y_OUT(:,39) = aS.*S0 ;% S0 activation term [#/無/d] +Y_OUT(:,40) = aN.*N0 ;% N0 activation term [#/無/d] +Y_OUT(:,41) = pS*Sap ;% Sa proliferation term [#/無/d] +Y_OUT(:,42) = pN*Nap ;% Na proliferation term [#/無/d] +Y_OUT(:,43) = dA*Sa ;% Sa death term [#/無/d] +Y_OUT(:,44) = dA*Na ;% Na death term [#/無/d] +Y_OUT(:,45) = mS*Saf ;% Sa reversion term [#/無/d] +Y_OUT(:,46) = mN*Naf ;% Na reversion term [#/無/d] +Y_OUT(:,47) = aS ;% S0 activation rate [/d] +Y_OUT(:,48) = aN ;% N0 activation rate [/d] +Y_OUT(:,49) = F(:,03) ;% C multiplier for aS* +Y_OUT(:,50) = ( F(:,04) + F(:,12) )./F(:,12) ;% C multiplier for aS* +Y_OUT(:,51) = F(:,06) ;% R multiplier for aS* +Y_OUT(:,52) = F(:,07) ;% R multiplier for aN* +Y_OUT(:,53) = F(:,11) ;% V multiplier for aS* +Y_OUT(:,54) = F(:,12) ;% V multiplier for aN* Y_OUT(:,55) = q*V ;% V growth term [#/無/d] -Y_OUT(:,56) = gE.*Ea.*V ;% E killing term [#/無/d] -Y_OUT(:,57) = gB.*Ba.*V ;% B killing term [#/無/d] -Y_OUT(:,58) = gE.*Ea ;% E killing rate [/d] -Y_OUT(:,59) = gB.*Ba ;% B killing rate [/d] -Y_OUT(:,60) = F(:,08) ;% R multiplier for gE* -Y_OUT(:,61) = F(:,09) ;% R multiplier for gB* -Y_OUT(:,62) = gB.*Ba ./ (gE.*Ea + gB.*Ba) ;% B killing fraction - -Y_OUT(:,63) = F(:,18) ;% R gen by aE -Y_OUT(:,64) = F(:,19) ;% R gen by aB -Y_OUT(:,65) = F(:,19) ./ (F(:,18) + F(:,19)) ;% R gen by aB fraction +Y_OUT(:,56) = gS.*Sa.*V ;% S killing term [#/無/d] +Y_OUT(:,57) = gN.*Na.*V ;% N killing term [#/無/d] +Y_OUT(:,58) = gS.*Sa ;% S killing rate [/d] +Y_OUT(:,59) = gN.*Na ;% N killing rate [/d] +Y_OUT(:,60) = F(:,08) ;% R multiplier for gS* +Y_OUT(:,61) = F(:,09) ;% R multiplier for gN* +Y_OUT(:,62) = gN.*Na ./ (gS.*Sa + gN.*Na) ;% N killing fraction + +Y_OUT(:,63) = F(:,18) ;% R gen by aS +Y_OUT(:,64) = F(:,19) ;% R gen by aN +Y_OUT(:,65) = F(:,19) ./ (F(:,18) + F(:,19)) ;% R gen by aN fraction % Converting most outputs to log-value logID = [ 3:15 , 18:19 , 31:61 , 63:64 ] ; @@ -271,33 +271,33 @@ 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) - B0 = Y(11) ;% resting bystander CD8+ T cells [#/無] - Ba = Y(12:13) ;% active bystander T cells [#/無] (B1 to B2) + S0 = Y(02) ;% resting specific CD8+ T cells [#/無] + Sa = Y(03:10) ;% active specific CD8+ T cells [#/無] (S1 to S8) + N0 = Y(11) ;% resting non-SIV-specific CD8+ T cells [#/無] + Na = Y(12:13) ;% active non-SIV-specific T cells [#/無] (N1 to N2) 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' + F = terms(Y') ;% modifiers for 'gS0','gN0','p0','aS0','aN0' - p = F(13) ;% E0/B0 proliferation rate [/d] - aE = F(14) ;% E0 activation rate [/d] - aB = F(15) ;% B0 activation rate [/d] - gE = F(16) ;% Ea killing rate [#/無-d] - gB = F(17) ;% Ba killing rate [#/無-d] + p = F(13) ;% S0/N0 proliferation rate [/d] + aS = F(14) ;% S0 activation rate [/d] + aN = F(15) ;% N0 activation rate [/d] + gS = F(16) ;% Sa killing rate [#/無-d] + gN = F(17) ;% Na killing rate [#/無-d] s = F(18)+F(19) ;% regulation generation 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) = + 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 + % calculate rates of change for V,S,N,X,C,R + dY(01) = q*V - gS*sum(Sa)*V - gN*sum(Na)*V ;% dV/dt + dY(02) = p*S0 - aS*S0 - d*S0 + mS*Sa(8) ;% dS0/dt + dY(03) = + 2*aS*S0 - dA*Sa(1) - pS*Sa(1) ;% dS1/dt + dY(4:9)= 2*pS*Sa(1:6) - dA*Sa(2:7) - pS*Sa(2:7) ;% dS2-7/dt + dY(10) = 2*pS*Sa(7) - dA*Sa(8) - mS*Sa(8) ;% dS8/dt + dY(11) = p*N0 - aN*N0 - d*N0 + mN*Na(2) ;% dN0/dt + dY(12) = + 2*aN*N0 - dA*Na(1) - pN*Na(1) ;% dN1/dt + dY(13) = 2*pN*Na(1) - dA*Na(2) - mN*Na(2) ;% dN2/dt dY(14) = - ka*X ;% dX/dt dY(15) = ka*X/vd - ke*C ;% dC/dt dY(16) = s - dR*R(1) ;% dR1/dt @@ -308,78 +308,78 @@ 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 [#/無] - B0 = Y(11) ;% resting bystander CD8+ T cells [#/無] - Ba = Y(12)+Y(13) ;% active bystander T cells [#/無] + S0 = Y(02) ;% resting specific CD8+ T cells [#/無] + Sa = sum(Y(03:10)) ;% active specific CD8+ T cells [#/無] + N0 = Y(11) ;% resting non-SIV-specific CD8+ T cells [#/無] + Na = Y(12)+Y(13) ;% active non-SIV-specific T cells [#/無] C = Y(15) ;% N803 plasma concentration [pM] - F = terms(Y') ;% modifiers for 'gE0','gB0','p0','aE0','aB0' + F = terms(Y') ;% modifiers for 'gS0','gN0','p0','aS0','aN0' - p = F(13) ;% E0/B0 proliferation rate [/d] - aE = F(14) ;% E0 activation rate [/d] - aB = F(15) ;% B0 activation rate [/d] - gE = F(16) ;% Ea killing rate [#/無-d] - gB = F(17) ;% Ba killing rate [#/無-d] + p = F(13) ;% S0/N0 proliferation rate [/d] + aS = F(14) ;% S0 activation rate [/d] + aN = F(15) ;% N0 activation rate [/d] + gS = F(16) ;% Sa killing rate [#/無-d] + gN = F(17) ;% Na killing rate [#/無-d] J = zeros(17) ; - % Partial derivatives of gE,gB,p,aE,aB w.r.t. V,E,B,C,R + % Partial derivatives of gS,gN,p,aS,aN w.r.t. V,S,N,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) ; + daSdR = - aS * aS2 * F(06) ; + daNdR = - aN * aN2 * F(07) ; + dgSdR = - gS * gS2 * F(08) ; + dgNdR = - gN * gN2 * F(09) ; + dpdSN = - p * F(10) / SN50 ; + daSdV = aS0 * V50S / (V50S+V)^2 * F(03) * F(06) ; + daNdV = aN0 * V50N / (V50N+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) ; + daSdC = aS0 * F(11) * aS1 * dOdC * F(06) ; + daNdC = aN0 * aN1 * 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 ; + J(01,01 ) = q - gS*Sa - gN*Na ; + J(01,02:10) = - gS * V ; + J(01,11:13) = - gN * V ; + J(01,17 ) = - ( dgSdR*Sa + dgNdR*Na ) * 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 f2 = dS0/dt + J(02,01 ) = - daSdV * S0 ; + J(02,02 ) = p + dpdSN*S0 - d - aS ; + J(02,03:13) = dpdSN*S0 ; + J(02,10 ) = dpdSN*S0 + mS ; + J(02,15 ) = ( dpdC - daSdC ) * S0 ; + J(02,17 ) = ( dpdR - daSdR ) * S0 ; - % Partial derivatives of f3 = dE1/dt (save f3/dE1) - 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 f3 = dS1/dt (save f3/dS1) + J(03,01 ) = 2*daSdV * S0 ; + J(03,02 ) = 2*aS ; + J(03,15 ) = 2*daSdC * S0 ; + J(03,17 ) = 2*daSdR * S0 ; - % Partial derivatives of f4-10 = dE2-8/dt - J(37:18:145) = - dA - pE ; - J(38:18:146) = 2*pE ; - J(10,10) = - dA - mE ; + % Partial derivatives of f4-10 = dS2-8/dt + J(37:18:145) = - dA - pS ; + J(38:18:146) = 2*pS ; + J(10,10) = - dA - mS ; - % 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 f11 = dN0/dt + J(11,01 ) = - daNdV*N0 ; + J(11,02:12) = dpdSN*N0 ; + J(11,11 ) = p + dpdSN*N0 - d - aN ; + J(11,13 ) = dpdSN*N0 + mN ; + J(11,15 ) = ( dpdC - daNdC ) * N0 ; + J(11,17 ) = ( dpdR - daNdR ) * N0 ; - % Partial derivatives of f12 = dB1/dt - J(12,01) = 2*daBdV * B0 ; - J(12,11) = 2*aB ; - J(12,12) = - dA - pB ; - J(12,15) = 2*daBdC * B0 ; - J(12,17) = 2*daBdR * B0 ; + % Partial derivatives of f12 = dN1/dt + J(12,01) = 2*daNdV * N0 ; + J(12,11) = 2*aN ; + J(12,12) = - dA - pN ; + J(12,15) = 2*daNdC * N0 ; + J(12,17) = 2*daNdR * N0 ; - % Partial derivatives of f12-13 = dB1-2/dt - J(13,12) = 2*pB ; - J(13,13) = - dA - mB ; + % Partial derivatives of f12-13 = dN1-2/dt + J(13,12) = 2*pN ; + J(13,13) = - dA - mN ; % Partial derivatives of f14 = dX/dt J(14,14) = -ka ; @@ -389,8 +389,8 @@ J(15,15) = -ke ; % Partial derivatives of f16 = dR1/dt - J(16,01) = sE/aE0*daEdV/F(06) + sB/aB0*daBdV/F(07) ; - J(16,15) = sE/aE0*daEdC/F(06) + sB/aB0*daBdC/F(07) ; + J(16,01) = sS/aS0*daSdV/F(06) + sN/aN0*daNdV/F(07) ; + J(16,15) = sS/aS0*daSdC/F(06) + sN/aN0*daNdC/F(07) ; J(16,16) = -dR ; % Partial derivatives of f17 = dR2/dt @@ -401,7 +401,7 @@ function F = terms(Y) V = Y(:,01) ;% SIV virions [#/無] - EB = sum(Y(:,02:13),2) ;% CD8+ T cells [#/無] + SN = sum(Y(:,02:13),2) ;% CD8+ T cells [#/無] C = Y(:,15) ;% N803 plasma concentration [pM] R = Y(:,17) ;% immune regulation [] @@ -410,28 +410,28 @@ F(:,01) = C./ (C50+C) ;% N803 effect [] - 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(:,02) = 1 + p1*F(:,1) ;% S0/N0 proliferation stimulation [] + F(:,03) = 1 + aS1*F(:,1) ;% S activation stimulation [] + F(:,04) = aN1*F(:,1) ;% N803-induced N activation [] - 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(:,05) = 1./ (1+ p2*R) ;% S0/N0 proliferation regulation [] + F(:,06) = 1./ (1+aS2*R) ;% Sa activation regulation [] + F(:,07) = 1./ (1+aN2*R) ;% Na activation regulation [] + F(:,08) = 1./ (1+gS2*R) ;% Sa killing regulation [] + F(:,09) = 1./ (1+gN2*R) ;% Na killing regulation [] - 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(:,10) = SN50./(SN50+SN) ;% S0/N0 proliferation density dependence [] + F(:,11) = V./ (V50S+V) ;% Sa activation antigen dependence [] + F(:,12) = V./ (V50N+V) ;% Na activation antigen dependence [] - F(:,13) = p0 * F(:,10).* F(:,02).* F(:,05) ;% E0/B0 proliferation [/d] - F(:,14) = aE0* F(:,11).* F(:,03).* F(:,06) ;% E0 activation rate [/d] - F(:,15) = aB0* (F(:,12)+F(:,04)).* F(:,07) ;% B0 activation rate [/d] - F(:,16) = gE0* F(:,08) ;% Ea killing rate [#/無-d] - F(:,17) = gB0* F(:,09) ;% Ba killing rate [#/無-d] + F(:,13) = p0 * F(:,10).* F(:,02).* F(:,05) ;% S0/N0 proliferation [/d] + F(:,14) = aS0* F(:,11).* F(:,03).* F(:,06) ;% S0 activation rate [/d] + F(:,15) = aN0* (F(:,12)+F(:,04)).* F(:,07) ;% N0 activation rate [/d] + F(:,16) = gS0* F(:,08) ;% Sa killing rate [#/無-d] + F(:,17) = gN0* F(:,09) ;% Na killing rate [#/無-d] - F(:,18) = sE* F(:,11).* F(:,03) ;% regulation generation by aE [/d] - F(:,19) = sB* (F(:,12) + F(:,04)) ;% regulation generation by aB [/d] + F(:,18) = sS* F(:,11).* F(:,03) ;% regulation generation by aS [/d] + F(:,19) = sN* (F(:,12) + F(:,04)) ;% regulation generation by aN [/d] end end \ No newline at end of file diff --git a/N803_shared.m b/N803_shared.m new file mode 100644 index 0000000..0510cb1 --- /dev/null +++ b/N803_shared.m @@ -0,0 +1,255 @@ +%% N803_shared.m - solves model for 2 cohorts with shared parameters +% +% /--------------------------------------------------------------\ +% | Date: 12/25/2022 | +% | Author: Jonathan Cody | +% | Affiliation: Purdue University | +% | Weldon School of Biomedical Engineering | +% | Pienaar Computational Systems Pharmacology Lab | +% \--------------------------------------------------------------/ +% +% Nomenclature: V = SIV virions [#/無] +% T8 = total CD8+ T cells [#/無] +% S0 = resting SIV-specific CD8+ T cells [#/無] +% Sa = active SIV-specific CD8+ T cells [#/無] +% N0 = resting non-SIV-specific CD8+ T cells [#/無] +% Na = active non-SIV-specific 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{c} = ascending vector of days at which to administer doses +% (elements of 'DoseTimes' must also be in 'SoluTimes') +% +% AllPars = vector of parameters (see list in function) +% +%% ======================================================================== +% OPTIONS +% ======================================================================== +% +% SkipTimes{c} = [min max] time point beyond which to skip model solving +% (outputs before 'min' will be made equal to output at 'min') +% (leave as [] to ignore and solve for all 'SoluTimes') +% +% oneCohort = scalar to run model for just one cohort ('1' or '2') +% (leave as [] to ignore and solve for both cohorts) +% +% All additional inputs will be passed as a cell vector to 'N803_model_2' +% and used to define options (see function for list) +% EX: N803_single(SoluTimes,DoseTimes,AllPars,'AbsTol',1e-2} +% will set ode solver absolute tolerance to 1e-2 +% +%% ======================================================================== +% OUTPUTS +% ======================================================================== +% +% Y_OUT(:,1) = V at points in 'SoluTimes' [log fold change] cohort 1 +% Y_OUT(:,2) = T8 at points in 'SoluTimes' [fold change] cohort 1 +% Y_OUT(:,3) = V at points in 'SoluTimes' [log fold change] cohort 2 +% Y_OUT(:,4) = T8 at points in 'SoluTimes' [fold change] cohort 2 +% +% PARS(1,:) = parameters for cohort 1 (see code) +% PARS(2,:) = parameters for cohort 2 (see code) +% +%% ======================================================================== +% FUNCTION +% ======================================================================== +function [Y_OUT,PARS] = ... + N803_shared(SoluTimes,DoseTimes,AllPars,SkipTimes,oneCohort,... + varargin) + +if isempty(oneCohort) ; RunCohort = [ 1 1 ] ; +elseif oneCohort == 1 ; RunCohort = [ 1 0 ] ; +else ; RunCohort = [ 0 1 ] ; +end + +Y_Cohort = cell(1,2) ;% cell for storing outputs +P_Cohort = cell(1,2) ;% cell for storing parmeters + +% Rename inputed parameters ----------------------------------------------- +Vi(1) = AllPars(01) ;% V initial value [log(#/mL)] (cohort 1) +Vi(2) = AllPars(02) ;% V initial value [log(#/mL)] (cohort 2) +SNi(1) = AllPars(03) ;% S+N initial value [#/無] (cohort 1) +SNi(2) = AllPars(04) ;% S+N initial value [#/無] (cohort 2) +fS(1) = AllPars(05) ;% initial frequency: S/(S+N) (cohort 1) +fSn = AllPars(06) ;% normalized S/(S+N) (co 2) +fS(2) = fS(1)+(0.3-fS(1))*fSn ;% initial frequency: S/(S+N) (cohort 2) + +q = AllPars(07) ;% V growth rate (if S+N were absent) [/d] +psi = AllPars(08) ;% ratio of base killing rates gN0/gS0 +V50S = AllPars(09) ;% 50% viral stimulation saturation for S [#/mL] +V50N = AllPars(10) ;% 50% viral stimulation saturation for N [#/mL] +mSn = AllPars(11) ;% normalized Sa reversion rate constant [] +mNn = AllPars(12) ;% normalized Na reversion rate constant [] + +SN50 = AllPars(13) ;% 50% S+N proliferation saturation [#/無] +p0n = AllPars(14) ;% nomalized S0/N0 proliferation rate constant [/d] +pS = AllPars(15) ;% Sa proliferation rate constant [/d] +pN = AllPars(16) ;% Na proliferation rate constant [/d] +d = AllPars(17) ;% S0/N0 death rate constant [/d] +dA = AllPars(18) ;% Sa/Na death rate constant [/d] + +Xi = AllPars(19) ;% X initial condition [pmol/kg] +ka = AllPars(20) ;% N803 absorption rate constant [/d] +ke = AllPars(21) ;% N803 elimination rate constant [/d] +vd = AllPars(22) ;% N803 'volume of distribution'/'bioavailability' [L/kg] +C50 = AllPars(23) ;% 50% N803 stimulation concentration [pM] (Cohort 1) +pm = AllPars(24) ;% S0/N0 maximum proliferation rate [] +aS1 = AllPars(25) ;% S activation stimulation factor [] +aN1 = AllPars(26) ;% N activation stimulation factor [] + +dR = AllPars(27) ;% R decay rate constant [/d] +sig = AllPars(28) ;% ratio of initial regulation due to N/S +p2 = AllPars(29) ;% S0/N0 proliferation regulation factor [] +gN2 = AllPars(30) ;% N killing regulation factor [] (cohort 1) + +%% ------------------------------------------------------------------------ +% Calculate some initial conditions & parameters -------------------------- + +Vi = 10.^(Vi - 3) ;% converting V initial value [#/無] +V50S = V50S/1000 ;% 50% viral stimulation saturation for S [#/無] +V50N = V50N/1000 ;% 50% viral stimulation saturation for N [#/無] +YS = Vi ./ (Vi+ V50S) ;% initial V/(V50S+V) [cohort 1,2] +YN = Vi ./ (Vi+ V50N) ;% initial V/(V50N+V) [cohort 1,2] + +% calculate initial R, and sS,sN +z = YS + sig*YN ; +Ri = [ 1 , (z(2)/z(1)) ] ;% initial regulation [cohort 1,2] +R = Ri(2) ;% initial regulation (cohort 2) +sS = dR / ( YS(1) + sig*YN(1) ) ;% R generation due to S0 activation [/d] +sN = sig*sS ;% R generation due to N0 activation [/d] + +% restrict mS and mN such that initial activation aS and aN are positive +US = 2*(2*pS/(pS+dA))^7 ; +UN = 2* 2*pN/(pN+dA) ; +mS = mSn*dA/(US-1) ;% Sa reversion rate constant [/d] +mN = mNn*dA/(UN-1) ;% Na reversion rate constant [/d] + +% restrict p0 to ensure same sign for S0/SA and for N0/NA +p0_max = min( d*(1+p2*Ri).*(SN50+SNi)/SN50 ) ;% maximum value for p0 +p0 = p0n * p0_max ;% S0/N0 prolif rate constant [/d] +p1 = pm/p0 ;% S0/N0 proliferation stimulation factor +pi = p0*SN50./(SN50+SNi)./(1+p2*Ri) ;% initial S0/N0 prolif rate [/d] + +% calculate aS0,aN0 and aS2,aN2 +aSi = (d-pi)/( mS*US/(dA+mS) - 1 ) ;% initial S0 activation rate [/d] +aNi = (d-pi)/( mN*UN/(dA+mN) - 1 ) ;% initial N0 activation rate [/d] +z = aSi./YS ; +aS2 = ( z(1)-z(2) ) / ( R*z(2)-z(1) ) ;% S0 activation rate constant [/d] +aS0 = aSi(1) / YS(1) * (1+aS2) ;% S activation regulation factor [] +z = aNi./YN ; +aN2 = ( z(1)-z(2) ) / ( R*z(2)-z(1) ) ;% N0 activation rate constant [/d] +aN0 = aNi(1) / YN(1) * (1+aN2) ;% N activation regulation factor [] + +% solve for initial ratios below (based on active steady-state) +ZS = US/(mS+dA) ; +for i = 1:7 + ZS = ZS + 2*(2*pS)^(i-1)/(pS+dA)^i ;% SAi/aSi/S0i +end +ZN = UN/(mN+dA) + 2/(pN+dA) ;% NAi/aNi/N0i + +% solve for initial values of S0,Sa,N0,Na each cohort +Si = SNi.*fS ;% initial S +Ni = SNi.*(1-fS) ;% initial N +S0 = Si./(1+ZS*aSi) ;% initial S0 +N0 = Ni./(1+ZN*aNi) ;% initial N0 +SA = S0.*(ZS*aSi) ;% initial SA +NA = N0.*(ZN*aNi) ;% initial NA + +% calculate gS0,gN0 and gS2 +beta = psi * NA ./ (1+gN2*Ri) ; z = beta(1) - beta(2) ; +a = z*R ; b = SA(1)*R - SA(2) + z*(1+R) ; c = SA(1) - SA(2) + z ; +gS2 = ( -b + sqrt(b^2 - 4*a*c) ) / (2*a) ;% S killing regulation factor [] +gS0 = q / ( SA(1)/(1+gS2) + beta(1) ) ;% Sa killing rate constant [無/#-d] +gN0 = psi * gS0 ;% Na killing rate constant [無/#-d] + +%% Do for each cohort (NOT indenting loop) ================================ +for c = 1:2 + +% solve for initial S1-8 and N1-2 +S = zeros(1,8) ;% initial S1-S8 +S(1) = 2*aSi(c)*S0(c) / (dA+pS) ;% S1 +for i = 2:7 + S(i) = 2*pS*S(i-1) / (dA+pS) ;% S2 to S7 +end +S(8) = 2*pS*S(7) / (dA+mS) ;% S8 + +N(1) = 2*aNi(c)*N0(c) / (dA+pN) ;% N1 +N(2) = 2*pN*N(1) / (dA+mN) ;% N2 + +%% ------------------------------------------------------------------------ +% Prepare parameter and initial value vectors and call 'N803_model_2' ----- + +Pars(01) = q ;% V growth rate (if S+N were absent) [/d] +Pars(02) = gS0 ;% Sa killing rate constant [無/#-d] +Pars(03) = gN0 ;% Na killing rate constant [無/#-d] + +Pars(04) = V50S ;% 50% viral stimulation saturation for S [#/無] +Pars(05) = V50N ;% 50% viral stimulation saturation for N [#/無] +Pars(06) = aS0 ;% S0 activation rate constant [/d] +Pars(07) = aN0 ;% N0 activation rate constant [/d] +Pars(08) = mS ;% Sa reversion rate constant [/d] +Pars(09) = mN ;% Na reversion rate constant [/d] + +Pars(10) = SN50 ;% 50% S+N proliferation saturation [#/無] +Pars(11) = p0 ;% S0/N0 proliferation rate constant [/d] +Pars(12) = pS ;% Sa proliferation rate constant [/d] +Pars(13) = pN ;% Na proliferation rate constant [/d] +Pars(14) = d ;% S0/N0 death rate constant [/d] +Pars(15) = dA ;% Sa/Na death rate constant [/d] + +Pars(16) = Xi ;% X initial condition [pmol/kg] +Pars(17) = ka ;% N803 absorption rate constant [/d] +Pars(18) = ke ;% N803 elimination rate constant [/d] +Pars(19) = vd ;% N803 'volume of distribution'/'bioavailability' [L/kg] +Pars(20) = C50 ;% 50% N803 stimulation concentration [pM] +Pars(21) = p1 ;% S0/N0 proliferation stimulation factor [] +Pars(22) = aS1 ;% S activation stimulation factor [] +Pars(23) = aN1 ;% N activation stimulation factor [] + +Pars(24) = sS ;% R generation due to S0 activation [/d] +Pars(25) = sN ;% R generation due to N0 activation [/d] +Pars(26) = dR ;% R decay rate constant [/d] +Pars(27) = gS2 ;% S killing regulation factor [] +Pars(28) = gN2 ;% N killing regulation factor [] +Pars(29) = p2 ;% S0/N0 proliferation regulation factor [] +Pars(30) = aS2 ;% S activation regulation factor [] +Pars(31) = aN2 ;% N activation regulation factor [] + +% V S0-8 N0-2 X C R initial values +Yic = [ Vi(c) S0(c) S N0(c) N 0 0 Ri(c) Ri(c) ] ; + +% if any( [ Pars Yic ] < 0 ) +% error('Negative parameters or initial values.') +% end + +% If 'SkipTimes' is empty, do not skip any times +if isempty(SkipTimes{c}) ; SkipTimes{c} = [-inf inf] ; end + +idLo = SoluTimes < SkipTimes{c}(1) ;% index of early times to skip soln +idHi = SoluTimes > SkipTimes{c}(2) ;% index of later times to skip soln +idSol = ~ ( idLo | idHi ) ;% index of times in 'SoluTimes' to solve + +if RunCohort(c) == 1 + Y_COH = N803_model_2(SoluTimes(idSol),DoseTimes{c},Pars,Yic,varargin) ; + Y_LO = ones(sum(idLo),1)*Y_COH(1 ,:) ;% constant Y for early times + Y_HI = ones(sum(idHi),1)*Y_COH(end,:) ;% constant Y for later times + Y_COH = [ Y_LO ; Y_COH ; Y_HI ] ;%#ok % total 'solution' matrix + Y_Cohort{c} = Y_COH ; + P_Cohort{c} = [ Pars Yic ] ; +end + +end + +%% Prepare main outputs 'Y_OUT' and 'PARS' ================================ + +Y_OUT = [ Y_Cohort{1} , Y_Cohort{2} ] ; +PARS = [ P_Cohort{1} ; P_Cohort{2} ] ; + +end \ No newline at end of file diff --git a/N803_single.m b/N803_single.m new file mode 100644 index 0000000..88e5d1b --- /dev/null +++ b/N803_single.m @@ -0,0 +1,217 @@ +%% N803_single.m - solves model for 1 cohort +% +% /--------------------------------------------------------------\ +% | Date: 12/25/2022 | +% | Author: Jonathan Cody | +% | Affiliation: Purdue University | +% | Weldon School of Biomedical Engineering | +% | Pienaar Computational Systems Pharmacology Lab | +% \--------------------------------------------------------------/ +% +% Nomenclature: V = SIV virions [#/無] +% T8 = total CD8+ T cells [#/無] +% S0 = resting SIV-specific CD8+ T cells [#/無] +% Sa = active SIV-specific CD8+ T cells [#/無] +% N0 = resting non-SIV-specific CD8+ T cells [#/無] +% Na = active non-SIV-specific 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') +% +% AllPars = vector of parameters (see list in function) +% +%% ======================================================================== +% OPTIONS +% ======================================================================== +% +% SkipTimes = [min max] time point beyond which to skip model solving +% (outputs before 'min' will be made equal to output at 'min') +% (leave as [] to ignore) +% +% All additional inputs will be passed as a cell vector to 'N803_model_2' +% and used to define options (see function for list) +% EX: N803_single(SoluTimes,DoseTimes,AllPars,'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] +% +% Pars = [ parameters , initial conditions ] +% +%% ======================================================================== +% FUNCTION +% ======================================================================== +function [Y_OUT,Pars] = ... + N803_single(SoluTimes,DoseTimes,AllPars,SkipTimes,varargin) + +% Rename inputed parameters ----------------------------------------------- +Vi = AllPars(01) ;% V initial value [log(#/mL)] +SNi = AllPars(02) ;% S+N initial value [#/無] +fS = AllPars(03) ;% initial frequency: S/(S+N) +fSA = AllPars(04) ;% initial frequency: Sa/S +q = AllPars(05) ;% V growth rate (if S+N were absent) [/d] +psi = AllPars(06) ;% Na/Sa killing rate ratio [gN0/gS0] +V50S = AllPars(07) ;% 50% viral stimulation saturation for S [#/mL] +V50N = AllPars(08) ;% 50% viral stimulation saturation for N [#/mL] +mSn = AllPars(09) ;% normalized Sa reversion rate constant [] +mNn = AllPars(10) ;% normalized Na reversion rate constant [] +SN50 = AllPars(11) ;% 50% S+N proliferation saturation [#/無] +pS = AllPars(12) ;% Sa proliferation rate constant [/d] +pN = AllPars(13) ;% Na proliferation rate constant [/d] +d = AllPars(14) ;% S0/N0 death rate constant [/d] +dA = AllPars(15) ;% Sa/Na death rate constant [/d] +Xi = AllPars(16) ;% X initial condition [pmol/kg] +ka = AllPars(17) ;% N803 absorption rate constant [/d] +ke = AllPars(18) ;% N803 elimination rate constant [/d] +vd = AllPars(19) ;% N803 'volume of distribution'/'bioavailability' [L/kg] +C50 = AllPars(20) ;% 50% N803 stimulation concentration [pM] +pm = AllPars(21) ;% S0/N0 maximum proliferation rate [] +aS1 = AllPars(22) ;% S activation stimulation factor [] +aN1 = AllPars(23) ;% N activation stimulation factor [] +sig = AllPars(24) ;% sN/sS regulation generation rate ratio +dR = AllPars(25) ;% R decay rate constant [/d] +gS2 = AllPars(26) ;% initial S killing regulation [] +gN2 = AllPars(27) ;% initial N killing regulation [] +p2 = AllPars(28) ;% initial S0/N0 proliferation regulation [] +aS2 = AllPars(29) ;% initial S activation regulation [] +aN2 = AllPars(30) ;% initial N activation regulation [] + +%% ------------------------------------------------------------------------ +% Calculate some initial conditions & parameters -------------------------- + +Vi = 10^(Vi-3) ;% V initial value [#/無] +V50S = V50S/1000 ;% 50% viral stimulation saturation for S [#/無] +V50N = V50N/1000 ;% 50% viral stimulation saturation for N [#/無] + +% restrict mS and mN such that initial activation aS and aN are positive +US = 2*(2*pS/(pS+dA))^7 ; +UN = 2* 2*pN/(pN+dA) ; +mS = mSn*dA/(US-1) ;% Sa reversion rate constant [/d] +mN = mNn*dA/(UN-1) ;% Na reversion rate constant [/d] + +% solve for initial ratios below (based on active steady-state) +ZS = US/(mS+dA) ; +for i = 1:7 + ZS = ZS + 2*(2*pS)^(i-1)/(pS+dA)^i ;% SAi/aSi/S0i +end +ZN = UN/(mN+dA) + 2/(pN+dA) ;% NAi/aNi/N0i + +WS = 1 ; +for i = 1:7 + WS = WS + (mS+dA)*(pS+dA)^(i-1)/(2*pS)^i ;% SAi/S8i +end +WN = 1 + (mN+dA)/(2*pN) ;% NAi/N2i + +QS = mS/WS - 1/ZS ;% collection +QN = mN/WN - 1/ZN ;% collection + +fNA = 1/( 1 + QN/QS*(1-fSA)/fSA ) ;% initial frequency: Na/N + +% solve for S and N initial conditions +Si = SNi*fS ;% initial S +Ni = SNi*(1-fS) ;% initial N +SA = Si*fSA ;% initial Sa +S0 = Si*(1-fSA) ;% initial S0 +NA = Ni*fNA ;% initial Sa +N0 = Ni*(1-fNA) ;% initial S0 + +S = zeros(1,8) ;% initial S1-S8 +S(8) = SA/WS ;% S8 +S(7) = S(8) * (mS+dA)/(2*pS) ;% S7 +for i = 6:-1:1 + S(i) = S(i+1) * (pS+dA) / (2*pS) ;% S6 to S1 +end +N = NA/WN * [ (mN+dA) / (2*pN) , 1 ] ;% initial N1-N2 + +% solve for rate constants +aSi = S(1) / (2*S0) * (pS+dA) ;% initial activation rate for S0 +aNi = aSi * (US*mS/(mS+dA)-1) / (UN*mN/(mN+dA)-1) ;% for N0 +aS0 = aSi * (V50S+Vi)/Vi * (1+aS2) ;% S0 activation rate constant [/d] +aN0 = aNi * (V50N+Vi)/Vi * (1+aN2) ;% N0 activation rate constant [/d] + +pi = d - QS*SA/S0 ;% initial proliferation rate for S0/N0 +p0 = pi * (SN50+SNi)/SN50 * (1+p2) ;% S0/N0 proliferation rate con [/d] +p1 = pm/p0 ;% S0/N0 proliferation stimulation factor + +gS0 = q / ( SA/(1+gS2) + psi*NA/(1+gN2) ) ;% Sa killing rate con [無/#-d] +gN0 = psi*gS0 ;% Na killing rate constant [無/#-d] + +sS = dR / ( Vi/(V50S+Vi) + sig*Vi/(V50N+Vi) ) ;% regulation gen by S act +sN = sig*sS ;% regulation generation by N activation + +%% ------------------------------------------------------------------------ +% Prepare parameter and initial value vectors and call 'N803_model_2' ----- + +Pars(01) = q ;% V growth rate (if S+N were absent) [/d] +Pars(02) = gS0 ;% Sa killing rate constant [無/#-d] +Pars(03) = gN0 ;% Na killing rate constant [無/#-d] + +Pars(04) = V50S ;% 50% viral stimulation saturation for S [#/無] +Pars(05) = V50N ;% 50% viral stimulation saturation for N [#/無] +Pars(06) = aS0 ;% S0 activation rate constant [/d] +Pars(07) = aN0 ;% N0 activation rate constant [/d] +Pars(08) = mS ;% Sa reversion rate constant [/d] +Pars(09) = mN ;% Na reversion rate constant [/d] + +Pars(10) = SN50 ;% 50% S+N proliferation saturation [#/無] +Pars(11) = p0 ;% S0/N0 proliferation rate constant [/d] +Pars(12) = pS ;% Sa proliferation rate constant [/d] +Pars(13) = pN ;% Na proliferation rate constant [/d] +Pars(14) = d ;% S0/N0 death rate constant [/d] +Pars(15) = dA ;% Sa/Na death rate constant [/d] + +Pars(16) = Xi ;% X initial condition [pmol/kg] +Pars(17) = ka ;% N803 absorption rate constant [/d] +Pars(18) = ke ;% N803 elimination rate constant [/d] +Pars(19) = vd ;% N803 'volume of distribution'/'bioavailability' [L/kg] +Pars(20) = C50 ;% 50% N803 stimulation concentration [pM] +Pars(21) = p1 ;% S0/N0 proliferation stimulation factor [] +Pars(22) = aS1 ;% S activation stimulation factor [] +Pars(23) = aN1 ;% N activation stimulation factor [] + +Pars(24) = sS ;% R generation due to S0 activation [/d] +Pars(25) = sN ;% R generation due to N0 activation [/d] +Pars(26) = dR ;% R decay rate constant [/d] +Pars(27) = gS2 ;% S killing regulation factor [] +Pars(28) = gN2 ;% N killing regulation factor [] +Pars(29) = p2 ;% S0/N0 proliferation regulation factor [] +Pars(30) = aS2 ;% S activation regulation factor [] +Pars(31) = aN2 ;% N activation regulation factor [] + +% V S0-8 N0-2 X C R initial values +Yic = [ Vi S0 S N0 N 0 0 1 1 ] ; + +if any( [ Pars Yic ] < 0 ) + error('Negative parameters or initial values.') +end + +% If 'SkipTimes' is empty, do not skip any times +if isempty(SkipTimes) ; SkipTimes = [-inf inf] ; end + +idLo = SoluTimes < SkipTimes(1) ;% index of early times to skip solution +idHi = SoluTimes > SkipTimes(2) ;% index of later times to skip solution +idSol = ~ ( idLo | idHi ) ;% index of times in 'SoluTimes' to solve + +% solve model +Y_OUT = N803_model_2(SoluTimes(idSol),DoseTimes,Pars,Yic,varargin) ; + +Y_LO = ones(sum(idLo),1)*Y_OUT(1 ,:) ;% constant Y for early times +Y_HI = ones(sum(idHi),1)*Y_OUT(end,:) ;% constant Y for later times +Y_OUT = [ Y_LO ; Y_OUT ; Y_HI ] ;% total 'solution' matrix + +Pars = [ Pars Yic ] ;% adding initial conditions to parameter output + +end \ No newline at end of file