From 3f763ad2057d8235385dd808efec6412eee25bbb Mon Sep 17 00:00:00 2001 From: "Cody, Jonathan William" Date: Thu, 28 Oct 2021 20:16:22 -0400 Subject: [PATCH] Add files via upload --- N803_model_2.m | 67 ++++++++++++++++++++++++++------------------------ 1 file changed, 35 insertions(+), 32 deletions(-) diff --git a/N803_model_2.m b/N803_model_2.m index c38d1e6..29d9357 100644 --- a/N803_model_2.m +++ b/N803_model_2.m @@ -2,19 +2,17 @@ % during a subcutaneously administered regimen of N803 % /--------------------------------------------------------------\ -% | Date: 10/26/2021 | +% | Date: 10/27/2021 | % | Author: Jonathan Cody | % | Affiliation: Pienaar Computational Systems Pharmacology Lab | % | Weldon School of Biomedical Engineering | % | Purdue University | % \--------------------------------------------------------------/ -% Version 2C +% Version 2D % % TOODOO % -> reconsider p1 calculation -% -> remove regulation feedback -% -> add death to all activated cells % Nomenclature: V = SIV virions [#/µL] % E0 = resting SIV-specific CD8+ T cells [#/µL] @@ -90,8 +88,8 @@ 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] -mE = Pars(14) ;% Ea reversion to memory rate constant [/d] -mB = Pars(15) ;% Ba reversion to memory rate constant [/d] +mEn = Pars(14) ;% normalized Ea reversion rate constant [] +mBn = Pars(15) ;% normalized Ba reversion rate constant [] Xi = Pars(16) ;% X initial condition [pmol/kg] ka = Pars(17) ;% N803 absorption rate constant [/d] ke = Pars(18) ;% N803 elimination rate constant [/d] @@ -120,16 +118,21 @@ V50E = V50E/1000 ;% 50% viral stimulation saturation for E [#/µL] V50B = V50B/1000 ;% 50% viral stimulation saturation for B [#/µL] +% 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 = (2*mB/(mB+dA) - 1) / (2^7*mE/(mE+dA) - 1) ; +phi = (kB*mB/(mB+dA) - 1) / (kE*mE/(mE+dA) - 1) ; -% ratio between initial activation rate of E and initial regulation -ome = dR / (sE + sB/phi) ; +% solve for aE,aB initial values and E and B initial values +x = fzero(@(x)icFun(x),[-8 2]) ; +[~,aEi,aBi,Ei,Bi] = icFun(x) ; -% solve for aE,aB initial values and initial values (HVL cohort) -% options = optimset('TolX',1e-30) ; -x = fzero(@(x)icFun(x),[-6 6]) ;%,options) ; -[~,aEi,aBi,Ei,Bi,Ri] = 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) ) ; @@ -137,7 +140,7 @@ 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 - 1 ;% E0/B0 proliferation stimulation factor +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 [] @@ -229,7 +232,7 @@ Y_OUT(:,34) = log10( F(:,18)/F(1,18) ) ;% lfc in regulation generation % Ba regulation stimulation proportion -Y_OUT(:,35) = sB*aB0*F(:,17) ./ F(:,18) ; +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) ; @@ -241,18 +244,17 @@ %% ======================================================================== % NESTED FUNCTIONS % ======================================================================== -function [targ,aEi,aBi,Ei,Bi,Ri] = icFun(x) +function [targ,aEi,aBi,Ei,Bi] = icFun(x) % declare initial activation rates aE and aB - Ri = 10^x ; - aEi = ome*Ri ; + 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) = [ -pE*ones(1,7),-(dA+mE),aBi,2*pB ] ; - MAT(129:13:142) = [ -pB,-(dA+mB) ] ; + 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) ; @@ -288,17 +290,17 @@ 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 - pE*Ea(1) ;% dE1/dt - dY(4:9)= 2*pE*Ea(1:6) - 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 - 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(18) - dR*R ;% dR/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(2) - 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(18) - dR*R ;% dR/dt % calculate rates of change for tolerance if N ~= 0 @@ -342,7 +344,8 @@ F(:,16) = F(:,13).* F(:,05).* F(:,08) ;% Ea activation modifier [] F(:,17) = (F(:,14)+F(:,06)).* F(:,09) ;% Ba activation modifier [] - F(:,18) = sE*aE0*F(:,16) + sB*aB0*F(:,17) ;% regulation gen [/d] + F(:,18) = sE*aE0* F(:,13).* F(:,05) ... + + sB*aB0*( F(:,14) + F(:,06) ) ;% regulation gen [/d] end end \ No newline at end of file