diff --git a/ICsolvers/N803_shared.m b/ICsolvers/N803_shared.m index 605247d..fae86b7 100644 --- a/ICsolvers/N803_shared.m +++ b/ICsolvers/N803_shared.m @@ -1,7 +1,7 @@ %% N803_shared.m - solves model for 2 cohorts with shared parameters % % /--------------------------------------------------------------\ -% | Date: 07/02/2022 | +% | Date: 08/18/2022 | % | Author: Jonathan Cody | % | Affiliation: Purdue University | % | Weldon School of Biomedical Engineering | @@ -79,7 +79,7 @@ EBi(2) = AllPars(04) ;% E+B initial value [#/µL] (cohort 2) fE(1) = AllPars(05) ;% initial frequency: E/(E+B) (cohort 1) fEn = AllPars(06) ;% normalized E/(E+B) (co 2) -fE(2) = fE(1)+(0.5-fE(1))*fEn ;% initial frequency: E/(E+B) (cohort 2) +fE(2) = fE(1)+(0.3-fE(1))*fEn ;% initial frequency: E/(E+B) (cohort 2) q = AllPars(07) ;% V growth rate (if E+B were absent) [/d] psi = AllPars(08) ;% ratio of base killing rates gB0/gE0 @@ -126,27 +126,11 @@ sB = sig*sE ;% R generation due to B0 activation [/d] % restrict mE and mB such that initial activation aE and aB are positive -UE = (2*pE/(pE+dA))^7 ; -UB = 2*pB/(pB+dA) ; +UE = 2*(2*pE/(pE+dA))^7 ; +UB = 2* 2*pB/(pB+dA) ; mE = mEn*dA/(UE-1) ;% Ea reversion rate constant [/d] mB = mBn*dA/(UB-1) ;% Ba reversion rate constant [/d] -% solve for initial ratios below (based on active steady-state) -ZE = UE/(mE+dA) ; -for i = 1:7 - ZE = ZE + (2*pE)^(i-1)/(pE+dA)^i ;% EAi/aEi/E0i -end -ZB = 1/(pB+dA) + UB/(mB+dA) ;% BAi/aBi/B0i - -WE = 1 ; -for i = 1:7 - WE = WE + (mE+dA)*(pE+dA)^(i-1)/(2*pE)^i ;% EAi/E8i -end -WB = 1 + (mB+dA)/(2*pB) ;% BAi/B2i - -QE = mE/WE - 1/ZE ;% collection -QB = mB/WB - 1/ZB ;% collection - % restrict p0 to ensure same sign for E0/EA and for B0/BA p0_max = min( d*(1+p2*Ri).*(EB50+EBi)/EB50 ) ;% maximum value for p0 p0 = p0n * p0_max ;% E0/B0 prolif rate constant [/d] @@ -154,8 +138,8 @@ pi = p0*EB50./(EB50+EBi)./(1+p2*Ri) ;% initial E0/B0 prolif rate [/d] % calculate aE0,aB0 and aE2,aB2 -aEi = (d-pi)/QE/ZE ;% initial E0 activation rate [/d] -aBi = (d-pi)/QB/ZB ;% initial B0 activation rate [/d] +aEi = (d-pi)/( mE*UE/(dA+mE) - 1 ) ;% initial E0 activation rate [/d] +aBi = (d-pi)/( mB*UB/(dA+mB) - 1 ) ;% initial B0 activation rate [/d] z = aEi./YE ; aE2 = ( z(1)-z(2) ) / ( R*z(2)-z(1) ) ;% E0 activation rate constant [/d] aE0 = aEi(1) / YE(1) * (1+aE2) ;% E activation regulation factor [] @@ -163,6 +147,13 @@ aB2 = ( z(1)-z(2) ) / ( R*z(2)-z(1) ) ;% B0 activation rate constant [/d] aB0 = aBi(1) / YB(1) * (1+aB2) ;% B activation regulation factor [] +% solve for initial ratios below (based on active steady-state) +ZE = UE/(mE+dA) ; +for i = 1:7 + ZE = ZE + 2*(2*pE)^(i-1)/(pE+dA)^i ;% EAi/aEi/E0i +end +ZB = UB/(mB+dA) + 2/(pB+dA) ;% BAi/aBi/B0i + % solve for initial values of E0,Ea,B0,Ba each cohort Ei = EBi.*fE ;% initial E Bi = EBi.*(1-fE) ;% initial B @@ -183,12 +174,14 @@ % solve for initial E1-8 and B1-2 E = zeros(1,8) ;% initial E1-E8 -E(8) = EA(c)/WE ;% E8 -E(7) = E(8) * (mE+dA)/(2*pE) ;% E7 -for i = 6:-1:1 - E(i) = E(i+1) * (pE+dA) / (2*pE) ;% E6 to E1 +E(1) = 2*aEi(c)*E0(c) / (dA+pE) ;% E1 +for i = 2:7 + E(i) = 2*pE*E(i-1) / (dA+pE) ;% E2 to E7 end -B = BA(c)/WB * [ (mB+dA) / (2*pB) , 1 ] ;% initial B1-B2 +E(8) = 2*pE*E(7) / (dA+mE) ;% E8 + +B(1) = 2*aBi(c)*B0(c) / (dA+pB) ;% B1 +B(2) = 2*pB*B(1) / (dA+mB) ;% B2 %% ------------------------------------------------------------------------ % Prepare parameter and initial value vectors and call 'N803_model_2' ----- @@ -232,9 +225,9 @@ % V E0-8 B0-2 X C R initial values Yic = [ Vi(c) E0(c) E B0(c) B 0 0 Ri(c) Ri(c) ] ; -if any( [ Pars Yic ] < 0 ) - error('Negative parameters or initial values.') -end +% 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 diff --git a/ICsolvers/N803_single.m b/ICsolvers/N803_single.m index 6132764..117338b 100644 --- a/ICsolvers/N803_single.m +++ b/ICsolvers/N803_single.m @@ -1,7 +1,7 @@ %% N803_single.m - solves model for 1 cohort % % /--------------------------------------------------------------\ -% | Date: 07/01/2022 | +% | Date: 09/25/2022 | % | Author: Jonathan Cody | % | Affiliation: Purdue University | % | Weldon School of Biomedical Engineering | @@ -97,17 +97,17 @@ V50B = V50B/1000 ;% 50% viral stimulation saturation for B [#/µL] % restrict mE and mB such that initial activation aE and aB are positive -UE = (2*pE/(pE+dA))^7 ; -UB = 2*pB/(pB+dA) ; +UE = 2*(2*pE/(pE+dA))^7 ; +UB = 2* 2*pB/(pB+dA) ; mE = mEn*dA/(UE-1) ;% Ea reversion rate constant [/d] mB = mBn*dA/(UB-1) ;% Ba reversion rate constant [/d] % solve for initial ratios below (based on active steady-state) ZE = UE/(mE+dA) ; for i = 1:7 - ZE = ZE + (2*pE)^(i-1)/(pE+dA)^i ;% EAi/aEi/E0i + ZE = ZE + 2*(2*pE)^(i-1)/(pE+dA)^i ;% EAi/aEi/E0i end -ZB = 1/(pB+dA) + UB/(mB+dA) ;% BAi/aBi/B0i +ZB = UB/(mB+dA) + 2/(pB+dA) ;% BAi/aBi/B0i WE = 1 ; for i = 1:7 @@ -137,7 +137,7 @@ B = BA/WB * [ (mB+dA) / (2*pB) , 1 ] ;% initial B1-B2 % solve for rate constants -aEi = E(1) / E0 * (pE+dA) ;% initial activation rate for E0 +aEi = E(1) / (2*E0) * (pE+dA) ;% initial activation rate for E0 aBi = aEi * (UE*mE/(mE+dA)-1) / (UB*mB/(mB+dA)-1) ;% for B0 aE0 = aEi * (V50E+Vi)/Vi * (1+aE2) ;% E0 activation rate constant [/d] aB0 = aBi * (V50B+Vi)/Vi * (1+aB2) ;% B0 activation rate constant [/d]