-
Notifications
You must be signed in to change notification settings - Fork 0
/
Repressilator_model.py
114 lines (79 loc) · 2.43 KB
/
Repressilator_model.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sun Jun 6 22:43:47 2021
@author: savannah
"""
# ## simulate the Repressilator using an ODE model
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import math
# ### Define repressilator simulation
def sdot_repressilator(s,t,params):
km, km0, kdm, kp, kdp, K, n = params
m_tetR, m_lacI, m_cI, p_tetR, p_lacI, p_cI = s
rate_m_tetR_prod = (km*((K**n)/((K**n)+(p_lacI**n))))+km0
rate_m_lacI_prod = (km*((K**n)/((K**n)+(p_cI**n))))+km0
rate_m_cI_prod = (km*((K**n)/((K**n)+(p_tetR**n))))+km0
rate_p_tetR_prod = kp*m_tetR
rate_p_lacI_prod = kp*m_lacI
rate_p_cI_prod = kp*m_cI
rate_m_tetR_loss = kdm*m_tetR
rate_m_lacI_loss = kdm*m_lacI
rate_m_cI_loss = kdm*m_cI
rate_p_tetR_loss = kdp*p_tetR
rate_p_lacI_loss = kdp*p_lacI
rate_p_cI_loss = kdp*p_cI
dm_tetR = rate_m_tetR_prod - rate_m_tetR_loss
dm_lacI = rate_m_lacI_prod - rate_m_lacI_loss
dm_cI = rate_m_cI_prod - rate_m_cI_loss
dp_tetR = rate_p_tetR_prod - rate_p_tetR_loss
dp_lacI = rate_p_lacI_prod - rate_p_lacI_loss
dp_cI = rate_p_cI_prod - rate_p_cI_loss
ds = [dm_tetR, dm_lacI, dm_cI, dp_tetR, dp_lacI, dp_cI]
return ds
# ### Parameter values
km = 30
km0 = 0.03
kdm = 0.3466
kp = 6.931
kdp = 0.06931
K = 40
n = 2
params = [ km, km0, kdm, kp, kdp, K, n ]
# ### Initial conditions
m_tetR0 = 5
m_lacI0 = 0
m_cI0 = 0
p_tetR0 = 0
p_lacI0 = 0
p_cI0 = 0
s0 = [ m_tetR0, m_lacI0, m_cI0, p_tetR0, p_lacI0, p_cI0 ]
# ### Time observations
t_max = 1000
t_obs = np.linspace(0,t_max,t_max+1)
# ### Run simulation
s_obs = odeint(sdot_repressilator,s0,t_obs,args=(params,))
m_tetR_obs = s_obs[:,0]
m_lacI_obs = s_obs[:,1]
m_cI_obs = s_obs[:,2]
p_tetR_obs = s_obs[:,3]
p_lacI_obs = s_obs[:,4]
p_cI_obs = s_obs[:,5]
# ### Figure 3: Repressilator ODE simulation results
# Time vs Proteins per cell
vals =[3,4,5]
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
legends = ['tetR', 'lacI', 'cI']
colours = ['r','b','y']
ax.set_xlabel('Time (min)')
ax.set_ylabel('Proteins per cell')
ax.set_ylim(0,8000)
## CREATE TIMESERIES PLOT FULL BEHAVIOUR
ax.set_title('Protein production vs time')
ax.plot(t_obs, p_tetR_obs, '-',label='tetR',color='r')
ax.plot(t_obs, p_lacI_obs, '-',label='lacI', color='b')
ax.plot(t_obs, p_cI_obs, '-',label='cI', color='y')
ax.legend()