stochastic corrections

cblb/models.py

@@ -29,7 +29,7 @@ def not_cell_stochastic(state, params, Omega):
     # presume that the molecules are degraded in the same strain as they are produced
     N_Y = N_X
 
-    Omega *= N_X # reaction space volume is proportional to the number of cells
+    #Omega *= N_X # reaction space volume is proportional to the number of cells
 
     gamma_L_X *= Omega
     eta_x *= Omega
@@ -39,17 +39,17 @@ def not_cell_stochastic(state, params, Omega):
 
     p = [0]*5
 
-    #p[0] = N_X * gamma_L_X * (y ** n_y)/(1 + (theta_L_X*y)**n_y ) / Omega
-    p[0] = gamma_L_X * (y ** n_y)/(1 + (theta_L_X*y)**n_y ) / Omega # N_x already included in reaction space volume (Omega)
+    p[0] = N_X * gamma_L_X * (y ** n_y)/(1 + (theta_L_X*y)**n_y ) / Omega
+    #p[0] = gamma_L_X * (y ** n_y)/(1 + (theta_L_X*y)**n_y ) / Omega # N_x already included in reaction space volume (Omega)
 
-    #p[1] = N_X * delta_L * L_X 
-    p[1] = delta_L * L_X # N_x already included in reaction space volume (Omega)
+    p[1] = N_X * delta_L * L_X 
+    #p[1] = delta_L * L_X # N_x already included in reaction space volume (Omega)
 
-    #p[2] = N_X * (eta_x * (1/(1+ (omega_x*L_X)**m_x)))
-    p[2] = (eta_x * (1/(1+ (omega_x*L_X)**m_x))) # N_x already included in reaction space volume (Omega)
+    p[2] = N_X * (eta_x * (1/(1+ (omega_x*L_X)**m_x)))
+    #p[2] = (eta_x * (1/(1+ (omega_x*L_X)**m_x))) # N_x already included in reaction space volume (Omega)
 
-    #p[3] = N_Y * (delta_x * x)
-    p[3] = (delta_x * x) # N_y already included in reaction space volume (Omega)
+    p[3] = N_Y * (delta_x * x)
+    #p[3] = (delta_x * x) # N_y already included in reaction space volume (Omega)
 
     p[4] = rho_x * x
 
@@ -75,19 +75,19 @@ def yes_cell_stochastic(state, params, Omega):
     # presume that the molecules are degraded in the same strain as they are produced
     N_Y = N_X
 
-    Omega *= N_X # reaction space volume is proportional to the number of cells
+    #Omega *= N_X # reaction space volume is proportional to the number of cells
 
     gamma_x *= Omega
     theta_x /= Omega
 
 
     p = [0]*3
 
-    #p[0] = N_X * gamma_x * (y ** n_y)/(1 + (theta_x*y)**n_y )
-    p[0] = gamma_x * (y ** n_y)/(1 + (theta_x*y)**n_y ) # N_x already included in reaction space volume (Omega)
+    p[0] = N_X * gamma_x * (y ** n_y)/(1 + (theta_x*y)**n_y )
+    #p[0] = gamma_x * (y ** n_y)/(1 + (theta_x*y)**n_y ) # N_x already included in reaction space volume (Omega)
 
-    #p[1] = N_Y * (delta_x * x) 
-    p[1] = delta_x * x # N_y already included in reaction space volume (Omega)
+    p[1] = N_Y * (delta_x * x) 
+    #p[1] = delta_x * x # N_y already included in reaction space volume (Omega)
 
     p[2] = rho_x * x
 

cblb/run_ssa_model_clb.py

@@ -93,8 +93,11 @@ def simulate_stochastic_clb(params, Y0, Omega, T_end, dt = 1):
 Y0[42:58] = 1 # number of cells
 
 
-# reaction space volume for a single cell
-Omega = 10
+# reaction space volume for the whole cell population
+# N_cells should be set to 1
+Omega = 5
+
+
 
 t_end = 500
 
@@ -137,8 +140,8 @@ def simulate_stochastic_clb(params, Y0, Omega, T_end, dt = 1):
 
     if iteration:
         Y0 = Y_full[-1,:]        
-    else:
-        Y0 *= Omega
+    #else:
+    #    Y0 *= N_cells
 
     #print(Y0)
 

cblb/run_ssa_model_mux.py

@@ -58,10 +58,13 @@ def simulate_stochastic_mux(params, Y0, Omega, T_end, dt = 1):
 
 rho_x = 0 
 params = [delta_L, gamma_L_X, n_y, theta_L_X, eta_x, omega_x, m_x, delta_x, rho_x, gamma_x, theta_x, r_X]
+# reaction space volume for the whole cell population
+# N_cells should be set to 1
 Omega = 10
 
+
 # I0, I1, I2, I3
-I = np.array([0, 1, 0, 0])
+I = np.array([0, 1, 0, 0])*100
 # S0, S1
 S = np.array([1, 0])
 
@@ -73,7 +76,7 @@ def simulate_stochastic_mux(params, Y0, Omega, T_end, dt = 1):
 Y0[4:6] = S
 
 
-T, Y = simulate_stochastic_mux(params, Y0*Omega, Omega, 100)
+T, Y = simulate_stochastic_mux(params, Y0, Omega, 100)
 
 out = Y[:,-1]
 #plt.plot(T,out)

cblb/run_ssa_model_toggle.py

@@ -67,8 +67,11 @@ def simulate_stochastic_toggle(params, Y0, Omega, T_end):
 
 #params = delta_L, gamma_A, gamma_B, n_a, n_b, theta_A, theta_B, eta_a, eta_b, omega_a, omega_b, m_a, m_b, delta_a, delta_b, rho_a, rho_b, r_A, r_B
 params = [delta_L, gamma_B, gamma_B, n_b, n_b, theta_B, theta_B, eta_b, eta_b, omega_b, omega_b, m_b, m_b, delta_b, delta_b, rho_b, rho_b, r_B, r_B]
+# reaction space volume for the whole cell population
+# N_cells should be set to 1
 Omega = 10
 
+
 #Y = L_A, L_B, a, b, N_A, N_B
 Y0 = np.zeros(6)#([0.0]*6)
 # number of cells
@@ -80,7 +83,7 @@ def simulate_stochastic_toggle(params, Y0, Omega, T_end):
 Y0[2] = 1 # a
 #Y0[3] = 0.1 # b
 
-T, Y = simulate_stochastic_toggle(params, Y0*Omega, Omega, 300)
+T, Y = simulate_stochastic_toggle(params, Y0, Omega, 500)
 
 
 L_A = Y[:,0]

robustness/models.py

@@ -22,6 +22,7 @@ def not_cell(state, params):
 
     return dL_X_dt, dx_dt
 
+"""
 def not_cell_stochastic(state, params, Omega):
     L_X, x, y, N_X, N_Y = state
     delta_L, gamma_L_X, n_y, theta_L_X, eta_x, omega_x, m_x, delta_x, rho_x = params
@@ -54,7 +55,7 @@ def not_cell_stochastic(state, params, Omega):
     p[4] = rho_x * x
 
     return p
-
+"""
 
 def yes_cell(state, params):
     x, y, N_X, N_Y = state
@@ -67,7 +68,7 @@ def yes_cell(state, params):
     dx_dt = N_X * gamma_x * (y ** n_y)/(1 + (theta_x*y)**n_y ) - N_Y * (delta_x * x) - rho_x * x
 
     return dx_dt
-
+"""
 def yes_cell_stochastic(state, params, Omega):
     x, y, N_X, N_Y = state
     gamma_x, n_y, theta_x, delta_x, rho_x = params
@@ -92,7 +93,7 @@ def yes_cell_stochastic(state, params, Omega):
     p[2] = rho_x * x
     
     return p
-
+"""
 
 def population(state, params):
     N = state
@@ -101,7 +102,7 @@ def population(state, params):
     dN = r * N * (1 - N)    
 
     return dN
-
+"""
 def population_stochastic(state, params, Omega):
     N = state
     r = params
@@ -112,7 +113,7 @@ def population_stochastic(state, params, Omega):
     p[1] = r * Omega * N**2
 
     return p
-
+"""
 
 def toggle_model(state, T, params):
     L_A, L_B, a, b, N_A, N_B = state
@@ -132,7 +133,7 @@ def toggle_model(state, T, params):
     dN_B_dt = population(N_B, r_B)
 
     return np.array([dL_A_dt, dL_B_dt, da_dt, db_dt, dN_A_dt, dN_B_dt])
-
+"""
 def toggle_model_stochastic(state, params, Omega):
     L_A, L_B, a, b, N_A, N_B = state
 
@@ -213,7 +214,7 @@ def toggle_generate_stoichiometry():
         N[idx_b, r] = -1
 
         return N
-
+"""
 
 # L_A ... intermediate
 # a ... out