added toggle, solver changed

latch/models.py

@@ -0,0 +1,63 @@
+# Urrios 2016: multicellular memory
+
+import numpy as np
+
+def not_cell(state, params):
+    L_X, x, y, N_X, N_Y = state
+    delta_L, gamma_X, n_y, theta_X, eta_x, omega_x, m_x, delta_x, rho_x = params
+
+    f = gamma_X * (y ** n_y)/(1 + (theta_X*y)**n_y )
+    dL_X_dt = f - delta_L * L_X
+
+    dx_dt = N_X * (eta_x * (1/(1+ (omega_x*L_X)**m_x))) - N_Y * (delta_x * x) - rho_x * x
+
+    return dL_X_dt, dx_dt
+
+
+def not_cell_a(state, params):
+    return not_cell(state, params)
+
+
+"""
+    L_A, a, b, N_A, N_B = state
+    delta_L, gamma_A, n_b, theta_A, eta_a, omega_a, m_a, delta_a, rho_a= params
+
+    f_b = gamma_A * (b ** n_b)/(1 + (theta_A*b)**n_b )
+    dL_A_dt = f_b - delta_L * L_A
+
+    da_dt = N_A * (eta_a * (1/(1+ (omega_a*L_A)**m_b))) - N_B * (delta_a * a) - rho_a * a
+"""
+
+def not_cell_b(state, params):
+    return not_cell(state, params)
+
+def population(state, params):
+    N = state
+    r = params
+
+    dN = r * N * (1 - N)
+
+    return dN
+
+
+def toggle_model(state, T, params):
+    L_A, L_B, a, b, N_A, N_B = state
+
+    state_A = L_A, a, b, N_A, N_B
+    state_B = L_B, b, a, N_B, N_A
+
+    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
+
+    params_A = delta_L, gamma_A, n_b, theta_A, eta_a, omega_a, m_a, delta_a, rho_a
+    params_B = delta_L, gamma_B, n_a, theta_B, eta_b, omega_b, m_b, delta_b, rho_b
+
+    dL_A_dt, da_dt = not_cell(state_A, params_A)
+    dL_B_dt, db_dt = not_cell(state_B, params_B)
+
+    dN_A_dt = population(N_A, r_A)
+    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_ODE(T, state, params):
+    return toggle_model(state, T, params)

latch/parameters.py

@@ -0,0 +1,24 @@
+# see Supplementary Table 1 in Urrios 2016
+
+gamma_A = 0.615 #nM/min
+gamma_B = 0.495 #nM/min 
+mu_A = 2 #
+mu_B = 2 # 
+omega_a = 1550 #nM-1
+omega_b = 1550 #nM-1
+eta_a = 0.369#0.0369 #nM/min 
+eta_b = 0.162#0.162 #nM/min 
+r_A = 0.07 
+r_B = 0.07 
+m_a = 2 
+m_b = 2
+delta_L = 0.15 #min-1
+delta_a = 0.05 #min-1
+delta_b = 0.023 #min-1
+theta_A = 0.26 #nM-1
+theta_B = 0.167 #nM-1
+n_a = 0.9
+n_b = 1.2
+rho_a = 5#5 #min-1
+rho_b = 5#5 #min-1
+

latch/test_model.py

@@ -0,0 +1,77 @@
+from scipy.integrate import odeint
+import matplotlib.pyplot as plt
+
+from models import *
+from parameters import *
+
+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_A, gamma_A, n_a, n_a, theta_A, theta_A, eta_a, eta_a, omega_a, omega_a, m_a, m_a, delta_a, delta_a, rho_a, rho_a, r_A, r_A]
+#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]
+
+# simulation parameters
+t_end = 500
+N = 1000
+
+# Y = L_A, L_B, a, b, N_A, N_B
+Y0 = np.array([0]*6)
+Y0[-2] = 1
+Y0[-1] = 1
+Y0[2] = 0
+T1 = np.linspace(0, t_end, N)
+T2 = np.linspace(0, t_end, N)
+T3 = np.linspace(0, t_end, N)
+T4 = np.linspace(0, t_end, N)
+T5 = np.linspace(0, t_end, N)
+
+params[-4] = rho_a
+params[-3] = 0
+
+Y1 = odeint(toggle_model, Y0, T1, args=((tuple(params),)))
+
+params[-4] = 0
+Y2 = odeint(toggle_model, Y1[-1], T2, args=((tuple(params),)))
+
+params[-3] = rho_b
+Y3 = odeint(toggle_model, Y2[-2], T3, args=((tuple(params),)))
+
+params[-3] = 0
+Y4 = odeint(toggle_model, Y3[-2], T4, args=((tuple(params),)))
+
+params[-3] = 0
+Y5 = odeint(toggle_model, Y4[-2], T5, args=((tuple(params),)))
+
+T2 += T1[-1]
+T3 += T2[-1]
+T4 += T3[-1]
+T5 += T4[-1]
+
+Y = np.append(np.append(np.append(np.append(Y1, Y2, axis=0), Y3, axis=0), Y4, axis=0), Y5, axis=0)
+T = np.append(np.append(np.append(np.append(T1, T2, axis=0), T3, axis=0), T4, axis=0), T5, axis=0)
+
+L_A = Y[:,0]
+L_B = Y[:,1]
+a = Y[:,2]
+b = Y[:,3]
+N_A = Y[:,4]
+N_B = Y[:,5]
+
+ax1 = plt.subplot(311)
+ax1.plot(T,L_A)
+ax1.plot(T,L_B)
+ax1.legend(["L_A", "L_B"])
+
+ax2 = plt.subplot(312)
+ax2.plot(T,a)
+ax2.plot(T,b)
+ax2.legend(["a", "b"])
+
+ax3 = plt.subplot(313)
+ax3.plot(T,N_A)
+ax3.plot(T,N_B)
+ax3.legend(["N_A", "N_B"])
+
+#plt.plot(Y)
+#plt.legend(["L_A", "L_B", "a", "b", "N_A", "N_B"])
+#plt.xticks([0, len(T)-1], [0, T[-1]])
+
+plt.show()

latch/test_model_ODE.py

@@ -0,0 +1,183 @@
+from scipy.integrate import ode
+import matplotlib.pyplot as plt
+
+from models import *
+from parameters import *
+
+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_A, gamma_A, n_a, n_a, theta_A, theta_A, eta_a, eta_a, omega_a, omega_a, m_a, m_a, delta_a, delta_a, rho_a, rho_a, r_A, r_A]
+#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]
+
+# simulation parameters
+t_end = 200
+N = t_end*2
+
+# Y = L_A, L_B, a, b, N_A, N_B
+Y0 = np.array([0]*6)
+Y0[-2] = 1
+Y0[-1] = 1
+Y0[2] = 0
+T1 = np.linspace(0, t_end, N)
+T2 = np.linspace(0, t_end, N)
+T3 = np.linspace(0, t_end, N)
+T4 = np.linspace(0, t_end, N)
+T5 = np.linspace(0, t_end, N)
+
+# 1
+
+params[-4] = rho_a
+params[-3] = 0
+
+t1 = t_end
+dt = t_end/N
+T1 = np.arange(0,t1+dt,dt)
+Y1 = np.zeros([1+N,6])
+Y1[0,:] = Y0
+
+r = ode(toggle_model_ODE).set_integrator('zvode', method='bdf')
+r.set_initial_value(Y0, T1[0]).set_f_params(params)
+
+i = 1
+while r.successful() and r.t < t1:
+    Y1[i,:] = r.integrate(r.t+dt)
+    i += 1
+
+# 2 
+
+params[-4] = 0
+params[-3] = 0
+
+t1 = t_end
+dt = t_end/N
+T2 = np.arange(0,t1+dt,dt)
+Y0 = Y1[-1,:]
+Y2 = np.zeros([1+N,6])
+Y2[0,:] = Y0
+
+r = ode(toggle_model_ODE).set_integrator('zvode', method='bdf')
+r.set_initial_value(Y0, T2[0]).set_f_params(params)
+i = 1
+while r.successful() and r.t < t1:
+    Y2[i,:] = r.integrate(r.t+dt)
+    i += 1
+
+T2 += t_end
+
+# 3
+
+params[-4] = 0
+params[-3] = rho_b
+
+t1 = t_end
+dt = t_end/N
+T3 = np.arange(0,t1+dt,dt)
+Y0 = Y2[-1,:]
+Y3 = np.zeros([1+N,6])
+Y3[0,:] = Y0
+
+r = ode(toggle_model_ODE).set_integrator('zvode', method='bdf')
+r.set_initial_value(Y0, T3[0]).set_f_params(params)
+i = 1
+while r.successful() and r.t < t1:
+    Y3[i,:] = r.integrate(r.t+dt)
+    i += 1
+
+T3 += 2*t_end
+
+# 4
+
+params[-4] = 0
+params[-3] = 0
+
+t1 = t_end
+dt = t_end/N
+T4 = np.arange(0,t1+dt,dt)
+Y0 = Y3[-1,:]
+Y4 = np.zeros([1+N,6])
+Y4[0,:] = Y0
+
+r = ode(toggle_model_ODE).set_integrator('zvode', method='bdf')
+r.set_initial_value(Y0, T4[0]).set_f_params(params)
+i = 1
+while r.successful() and r.t < t1:
+    Y4[i,:] = r.integrate(r.t+dt)
+    i += 1
+
+T4 += 3*t_end
+
+# 5
+
+params[-4] = rho_a
+params[-3] = 0
+
+t1 = t_end
+dt = t_end/N
+T5 = np.arange(0,t1+dt,dt)
+Y0 = Y4[-1,:]
+Y5 = np.zeros([1+N,6])
+Y5[0,:] = Y0
+
+r = ode(toggle_model_ODE).set_integrator('zvode', method='bdf')
+r.set_initial_value(Y0, T5[0]).set_f_params(params)
+i = 1
+while r.successful() and r.t < t1:
+    Y5[i,:] = r.integrate(r.t+dt)
+    i += 1
+
+T5 += 4*t_end
+
+
+"""
+Y1 = odeint(toggle_model, Y0, T1, args=((tuple(params),)))
+
+params[-4] = 0
+Y2 = odeint(toggle_model, Y1[-1], T2, args=((tuple(params),)))
+
+params[-3] = rho_b
+Y3 = odeint(toggle_model, Y2[-2], T3, args=((tuple(params),)))
+
+params[-3] = 0
+Y4 = odeint(toggle_model, Y3[-2], T4, args=((tuple(params),)))
+
+params[-3] = 0
+Y5 = odeint(toggle_model, Y4[-2], T5, args=((tuple(params),)))
+
+T2 += T1[-1]
+T3 += T2[-1]
+T4 += T3[-1]
+T5 += T4[-1]
+
+Y = np.append(np.append(np.append(np.append(Y1, Y2, axis=0), Y3, axis=0), Y4, axis=0), Y5, axis=0)
+T = np.append(np.append(np.append(np.append(T1, T2, axis=0), T3, axis=0), T4, axis=0), T5, axis=0)
+"""
+
+Y = np.append(np.append(np.append(np.append(Y1, Y2, axis=0), Y3, axis=0), Y4, axis=0), Y5, axis=0)
+T = np.append(np.append(np.append(np.append(T1, T2, axis=0), T3, axis=0), T4, axis=0), T5, axis=0)
+
+L_A = Y[:,0]
+L_B = Y[:,1]
+a = Y[:,2]
+b = Y[:,3]
+N_A = Y[:,4]
+N_B = Y[:,5]
+
+ax1 = plt.subplot(311)
+ax1.plot(T,L_A)
+ax1.plot(T,L_B)
+ax1.legend(["L_A", "L_B"])
+
+ax2 = plt.subplot(312)
+ax2.plot(T,a)
+ax2.plot(T,b)
+ax2.legend(["a", "b"])
+
+ax3 = plt.subplot(313)
+ax3.plot(T,N_A)
+ax3.plot(T,N_B)
+ax3.legend(["N_A", "N_B"])
+
+#plt.plot(Y)
+#plt.legend(["L_A", "L_B", "a", "b", "N_A", "N_B"])
+#plt.xticks([0, len(T)-1], [0, T[-1]])
+
+plt.show()

mux/check_T_f.py

@@ -3,7 +3,7 @@
 import numpy as np
 import matplotlib.pyplot as plt
 
-from data import *
+from parameters import *
 
 
 for name, max_val in zip(names[:], max_vals[:]):

mux/mux.py

@@ -2,7 +2,7 @@
 import pickle
 
 from transfer import *
-from data import *
+from parameters import *
 
 
 with open('direct_params.pickle', 'rb') as handle: