Add files via upload

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 [#/無] (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,43 +126,34 @@
 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]
 p1 = pm/p0         ;% E0/B0 proliferation stimulation factor 
 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 []
 z = aBi./YB  ;
 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

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 [#/無]
 
 % 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]