Team:SCU China/Model/Repressilator
Genetic regulatory networks are complex dynamic systems, and they can describe the interaction of gene expression of living cells. When designing novel system of wanted function, the intrinsic stabilization and regulation of kinetics activity could determine the outcome of the design, which highlight the crucial role of mathematical and physical modelling. Many works have been done to analyze the behavior of synthetic oscillators. The most well-known synthetic oscillator is the repressilator, which was the first one to be achieved with experiment by Elowitz&Leibler in 2000[7]. Their realized synthetic oscillation in one cell. In 2016, Potvin-Trottier and his coworker successfully measured the repressilator functioning features in a cluster of cells[16]. However, in order to enhance the synchronization in a population of cells, in our project, we coupled quorum sensing system with repressilator to synchronize repressilators oscillation pattern. In addition, to measure the time to reach quasi-steady-state, we included the delay time τ into the equations since time delays are omnipresent in actual complex systems[5][13][15]. The equations that we use to describe this behavior of the system are ordinary differential equations (ODE) and delay differential equations (DDE). At last, the model is shown below.
ExtraFig.1 Scheme of the repressilator network coupling a quorum-sensing system. m1, m2, m3 are the promoters of LacI, TetR and CI. The parameter β and τ correspond to the notation used later. Line with bar represents inhibition. The period of oscillation is synchronized through Ai signal(N-acyl homoserine lactones).
Firstly, let’s take a look at promoter. A promoter is a region of DNA just upstream of the coding region of a gene that acts as a regulator of the gene’s expression. There are generally three types of promoters used in synthetic biology: (1) promoters that are up-regulated by transcriptional activators, as shown in (2) promoters that are down-regulated by transcriptional repressors, as shown in Fig. 1.
Fig1. Taken from [15]. Regulation of promoters through transcriptional activators and repressors. (a) Transcriptional activators initiate transcription of a gene once bound to the promoter. In the absence of a transcriptional activator the gene is silenced, or off. (b) Transcriptional repressors inhibit transcription of a gene by binding to the promoter region. The removal of the transcriptional repressor allows the promoter to become transcriptionally active.
To derive the equation for our repressilator, we need to start with the equation of single protein. If the promoter fires at a rate αu when unbound and αb when bound, then the total rate at which initiation events occur is
On the other hand, the concentration of protein won’t increase without limitation. The most ubiquitous form of protein ‘‘degradation’’ is due to cellular growth and division, i.e. the protein concentration dilutes as the total reaction volume grows despite the real situation may be more complicated. For simplicity, we assume the cells’ volume grow exponentially with a rate μ as is shown in Fig.2, then we can write
When we introduced our project, we attracted people’s attention successfully. Many people showed their interest in our project in the panel session and they also gave us some good suggestions.
But sometimes this approximation of dynamics can be poor, the general way is to include the dynamics of related mRNA. The reason is that it actually requires quite a lot of steps to produce a functional protein. Of course, we can’t write down every step’s equation. It’s hard to do that, and it takes much time for computer to numerically calculate the result. However, generally, ODE including mRNA dynamics is accurate enough for most of the situations since some reactions like oligomerization, TF binding and protein folding occur on a much faster time scale. Hence, we can write a 2-dimensional system as
As for the system of repressilator: LacI, TetR, CI, following the work of Elowitz and Leibler[7], we can naturally expand equations to 3 groups and they are coupling with each other.
Here, m1, m2, m3 are the concentrations of LacI, TetR and CI. Let αi be the maximum production rate. For simplicity, we assume those constants are the same for each protein and mRNA, i.e. α1=α2=α3; n₁=n₂=n₃; K₁=K₂=K₃; λ₁=λ₂=λ₃; αleak1=αleak2=αleak3. By now, we have got the equations to represent the repressilator system. In fact, we can rescale concentrations and time to arrive at dimensionless equations. So that time is rescaled by mRNA lifetime, the protein concentration is written in units of K, the mRNA concentration is rescaled by translation efficiency (number of proteins produced per mRNA).
Here, i represents a gene in the circuit, j represents its predecessor. α and α0 are the dimensionless maximum production rate and leakiness. n here equals to 2. β is the ratio of protein decay rate and mRNA decay rate. Bifurcation analysis is done by [2][3][4][15]. The results showed that a stable oscillation should satisfy the constrain that β>0.07 approximately, α>αbifurcation=r0n+1+r0 and r0=(2/(n-2))1/n. Otherwise, the result will be like Fig.3.
However, the equations above is still too ideal. In real situation, those parameters may be different and even fluctuate through time. If we solve the equations above with no change, the curve will be the same and this system won’t oscillate as is shown in Fig.3. Consider the situation that β is the most marked parameter, which can affect the oscillation frequency markedly. We assume that β satisfies the Gaussian distribution. The standard deviation of β is βΔ. Also, we have to admit that the standard deviation of β is beyond our ability to measure it with experiments, but after all it is still one of simulating the fluctuation. The interesting fact is that the system will always oscillate with even small deviation as is shown in Fig.4.
Fig. 3. (a) is the concentration of protein with time, (b) is the phase diagram of two protein. We can see that the trajectory will go along the line and the turn back to the steady point (red dot). The parameters are chosen as α=216, α₀=α/1000, β=1, n=2 following [7].
(b): The phase diagram of (a).
(c): The concentration of proteins with β=1.01, 0.99, and 1.0.
(d): The phase diagram of (c).
The other parameters are all the same as Fig.2. However, similar oscillation occurs if we add relative error in , α and α₀. Conclusion: the continuous oscillation depends on the relative error of parameters. Further reaserch can focus on bifurcation analysis of β. However, bigger relative error only increases the time to get to the limit cycle and it has little influence on the period of oscillation.
However, the equations above is still too ideal. In real situation, those parameters may be different and even fluctuate through time. If we solve the equations above with no change, the curve will be the same and this system won’t oscillate as is shown in Fig.3. Consider the situation that β is the most marked parameter, which can affect the oscillation frequency markedly. We assume that β satisfies the Gaussian distribution. The standard deviation of β is βΔ. Also, we have to admit that the standard deviation of β is beyond our ability to measure it with experiments, but after all it is still one of simulating the fluctuation. The interesting fact is that the system will always oscillate with even small deviation as is shown in Fig.4.
So far it seems good. The oscillation of the proteins that we want is here. However, this system that we have derived only works in a single cell. We have no reason to believe that clusters of cells’ proteins will have synchronized oscillation as is shown in Fig.5. Here, we introduce another parameter imax represents the cell number.
(b) shows that if we consider the fluctuation of different cells, the total concentration oscillates but the difference between the peak and valley concentration will be much smaller compared to one cell situation. It doesn’t seem good if the concentration of Melatonin won’t change strongly. Of course, we have to admit that the value of βΔ is hard to measure, but it is at least a good fluctuation. A similar case analyzed the influence of βΔ[11], and it showed that larger βΔ comes with bigger difficulty to synchronize. Here, we no longer talk about the influence of βΔ and consider it equals to 0.05 all the time.
Aiming to solve this problem, we introduced another function into the system: Quorum sensing. We can replace the equation in (1.5) with:
The third term is the quorum sensing term, with assumption that it follows a standard Michaelis-Menten kinetics. The concentration of Ai inside each cell is denoted by Sᵢ, which is also scaled by its Michaelis constant. κ is the maximum contribution to LacI transcription in the presence of saturating amounts of AI. On the other hand, the dynamics of intracellular Ai concentration is influenced by cell division’s dilution (degradation term), synthesis and the diffusion with intercellular medium. To model the membrane diffusion between cell and the medium we applied the picture of diffusion through a long, thin channel to membrane transport. Thus we expect that the flux through the membrane will be of the form
Here the permeation constant of solute, Ps, is a number depending on both the membrane and the molecule whose permeation. In simple cases, the value of Ps roughly reflects the width of the pore, the thickness of the membrane (length of the pore), and the diffusion constant for the solute molecules. We can firstly imaging the cell as a perfect ball membrane with surface area A and volume V. The value of A and V is ~π μm2 and ~1 μm3. Data A, V and PS are taken from www.bioNumbers.org. Since
we can get that the concentration of Si’s diffusion obeys the equation (1.12) by solving the equations (1.8-1.11).
After adding the degradation and synthesis term and rewriting those redundant parameters, we get whole equation (1.13):
Here, we already defined a new parameter η called diffusion constant and let η=APs/V. k₀ and k₁ are the rate constant of degradation and synthesis (rescaled dimensionless constant). It contains three terms: degradation, synthesis and diffusion. As for the extracellular concentration of AI, with similar derivation, we can get:
Here ηₑₓₜ represents the diffusion constant if we consider the extracellular medium as the inside. Still for simplicity, we use the well-known quasi-steady-state approximation(QSSA) which is widely considered a good approximation to describe quorum sensing. Besides, we will use QSSA to replace quasi-steady-state approximation in what follows. Equation(1.15) is the main form of QSSA because we assume that dSₑ/dt=0. Equation(1.15) tells us that extracellular AI follows the evaluation of intracellular AI at all time.
Obviously, Q is linearly proportional to the cell density and its value lies between 0 and 1. With equation (1.15), we can reduce one dimension since it is hard to quantitatively describe the concentration of extracellular AI.
This equation system has(7×imax+1) dimensions. Hence, with the increase of imax, equation (1.16) requires more and more time to solve it numerically. Here, we will set imax=imax=100 and 900(later), and this is high enough for us to eliminate the accidental error caused by Gaussian Distribution parameter β as much as possible. Here, the number of equations to be numerically solved is 141. So let’s take a look at what we got solving (1.16) by changing value Q (cell density) as is shown in Fig. 6, 7 and 8. Here, we showed the evaluation of extracellular AI, which is a good way to analysis the stability of quorum sensing. If the extracellular AI oscillates in a stable period, it means that quorum sensing can synchronize the oscillation.
The reason is that we want to determine the time it takes for Se to reach a relative steady sate. By solving equation (1.17) with the rest of the equations in (1.14) and (1.16), we can obtain the evaluation of Se. For simplicity, we assume that there is only one cell in the medium, which means imax=1 in equation (1.14) and (1.16). Hence, our problem becomes much more simple: we only need to calculate the ratio of Se and Si and estimate the time to reach steady value. Stable value means that Q won’t oscillate sharply like the concentration of those proteins. The result of our calculation is shown in Fig.9. Here, we assume the degradation constant of extracellular AI k₂ is equal to the degradation constant of intracellular AI k₁. Besides, Obviously, β here will have no Gaussian Distribution if there is only one cell in the medium.
(a)
(b)
(c)
Fig.6. The parameters are chosen as Q=0.01, κ=20,
k0=1, k1=0.01 and others are the same as Fig.5. (a) is the concentration of
protein CI in different cells.
(b) is the concentration of protein CI
and
the average of every cell’s protein CI.
(c) is the concentration of
extracellular AI. Here, κ, k1, k2 are also hard to measure in our lab, so we
just followed [16].
a
b
c
Fig.7. The parameters are chosen as Q=0.4, and others are the same as Fig.6. (a) is the concentration of protein CI in different cells. (b) is the concentration of protein CI and the average of every cell’s protein CI. (c) is the concentration of extracellular AI.
a
b
c
Fig.8. The parameters are chosen as Q=0.8, and others are the same as Fig.5. (a) is the concentration of protein CI in different cells. (b) is the concentration of protein CI and the average of every cell’s protein CI. (c) is the concentration of extracellular AI. Conclusion: With Q become larger, which also means larger influence of quorum sensing or cell density, the system tends to synchronize. The concentration of Average protein CI will be more like a single cell situation. It is good to see this. We now can expect our Melatonin to express stably and strongly during the time we want.
So far, we have proved that repressilator system plus quorum sensing can synchronize the oscillation of a cluster of cells. However, there are still many aspects that we didn’t consider, since the real situation is much more complex than our simple model. In order to improve the system further, we intend to analyze the influence of QSSA. Before that, we may ask some questions about QSSA. Question 1: how long does the system take to reach quasi-steady-state or how long do intracellular travel through the cell membrane and the extracellular AI becomes synchronized with intracellular AI ? Question 2: If we consider the time that it takes to reach quasi-steady-state, what will happen to the system’s behavior? Through our research, it turned out that QSSA made the system lose some behavior.
To answer the question 1, we can firstly include equation (1.14) into the equation system instead of using QSSA (it ignores the differential equation of Se). Whereas, we need to rewrite equation (1.14) as
With the result of Fig.9 we successfully estimated the time for intracellular AI. Hence, we set the delay time as τ=4.
With the result of the delay time τ, we can introduce delay differential equation to describe the behavior of the system. What is the definition of delay differential equation? In mathematics, delay differential equations (DDEs) are a type of differential equation in which the derivative of the unknown function at a certain time is given in terms of the values of the function at previous times. Therefore, we rewrite the equation of Si (1.13) as:
The term Sₑ₍t-τ₎ means that the concentration of intracellular AI at time t depends on the concentration of extracellular AI at time t-τ. This change of equation is used to describe the time delay for intracellular AI to travel through membrane and reach steady state with extracellular AI. Then we can get the following results with the new version of equations (1.16) by changing the value Q (cell density):
Interesting fact: owing to the delay time to reach QSS, synchronization by quorum sensing becomes less effective. The phase lag of different cells will not decrease like the original equations in chapter 1.4, and fluctuation of concentration is everywhere. The phase lag brought by delay time will hamper the synchronization because the cell only feels the extracellular AI a moment ago. Although the result is less satisfying, but our new equations (1.18) is at least closer to the real situation. On the other hand, if we make the period of oscillation long enough, the fluctuation of Melatonin will be less. According to our result here, our project must make sure the cells that produce Melatonin must gather together to let quorum sensing synchronize their oscillation period.
With the analysis of chapter 1.6, we have already got a basic image of the behavior of the repressilator. Aiming to make our project more reliable, we need to control the period of oscillation. In our project, we intend to make the production of Melatonin follow people’s biological clock. Nevertheless, according to the results in chapter 1.6, the period with β=1.0 is around 12 times mRNA lifetime, which yields around 1 hour if we set mRNA lifetime equals to 5 minutes. Unfortunately, this is still hard to reach our design. Therefore, we are going to analyze the controlling parameters of period. Apparently, β is the marked parameter. The other parameters like α₀, n, κ andη mainly affect the expressing strength (concentration of peak value), while β can change the period of proteins. So, we set β=2.0, 1.0, 0.5, 0.2, 0.1 respectively and the results are shown in Fig.13. Of course, β here means the mean value, and the standard deviation βΔ still equals to 0.05.
From the Fig.13, we can apparently
conclude that the period will become
longer with the decrease of β. However, longer period comes along with
unsteadiness. We can see that if we set β<0.2, the system won’t even
oscillates (some of the cell). We can set βmin=0.2 as the critical value.
Below βmin, where mRNA decay rate is 5 times as much as the protein decay
rate, we’d better avoid that if we want our Melatonin produce steadily.
Approximately, we can make our mRNA decay a little bit faster like 2 to 4
times will be good. The period will increase roughly from one hour to 1.5
hours and 3 hours. Objectively speaking, this is still not up to standard.
But at least we can figure out other ways to increase mRNA lifetime. If we
can lengthen the mRNA lifetime to 4 times as much as the original one, we
may hopefully reach our aim.
What’s more, we can analyze the phase
diagram of protein CI in one cell and
average protein CI in every cell. The results are in Fig.14.
In order to see the oscillation of Melatonin intuitively, our team linked CFP gene to the CI gene. We set imax=900, which suggests that we have a 6301 equations system to be solved. We improved explicit Runge-Kutta method by adding the time delay term into the algorithm so that we can solve a DDE system. Fig.15 is a experimental picture from our team. Those green dots represent E.coli, and they are diluted so that we can see it clearly.
The simulation of repressilator is shown in Fig.16. There are three groups of figures. The left images are arrayplots. Each pixel represents a cell and the intensity of the fluorescence or concentration of Melatonin is represented as green color. The Right images are the average intensity or concentration of Melatonin.
So far, our modeling shows that with high Q value (means high cell density for quorum sensing) the concentration of the target protein or Melatonin will oscillate even there is fluctuation and time delay to the system. To be honest, it’s hard to judge the reliability of the modeling of repressilator, and some of the parameters we used may not be that accurate due to the fluctuations of the cells’ physical property. Biology system is complex, and we still have much room to improve our model like considering the concentration impact of cell division[1][8] or introducing stochastic differential equations[5]. However, we do have the confidence that our modeling can be an effective and reasonable reference to design our project or similar projects.
1. Bierbaum V, Klumpp S. Impact of the cell division cycle on gene circuits.
Phys Biol. 2015;12(6):066003.
2. BUSE O, KUZNETSOV A, ́REZ RAP. Existence of limit cycles in the repressilator
equations. Journal of Bifurcation and chaos. 2009;19.
3. Buse O, Perez R, Kuznetsov A. Dynamical properties of the repressilator
model. Phys Rev E Stat Nonlin Soft Matter Phys. 2010;81(6 Pt 2):066206.
4. BUZZI. CA, LLIBRE. J. HOPF BIFURCATION IN THE FULL REPRESSILATOR EQUATIONS.
2010.
5. Cai X. Exact stochastic simulation of coupled chemical reactions with delays.
J Chem Phys. 2007;126(12):124108.
6. Churilov AN, Medvedev A, Zhusubaliyev ZT. Impulsive Goodwin oscillator with
large delay: Periodic oscillations, bistability, and attractors. Nonlinear
Analysis: Hybrid Systems. 2016;21:171-83.
7. Elowitz MB, Leibler S. A synthetic oscillatory network of transcriptiona.
Nature. 2000;403.
8. Gonze D. Modeling the effect of cell division on genetic oscillators. J Theor
Biol. 2013;325:22-33.
9. Goodwin BC. Oscillatory Behavior in Enzymatic Control Process. Advances in
Enzyme Regulation. 1965.
10. J.D.Dockery, J.P.Keener. A mathematical model for quorum sensing in
pseudomonas aeruginosa. Mathematical Biology. 2001.
11. Jordi Garcia-Ojalvo MBE, and Steven H. Strogatz. Modeling a synthetic
multicellular clock: Repressilators coupled by quorum sensing. PNAS.
2004;101.
12. Kopell DMaN. Synchronizing Genetic Relaxation Oscillators by Intercell
Signaling. PNAS. 2002;99.
13. Ling G, Guan ZH, He DX, Liao RQ, Zhang XH. Stability and bifurcation
analysis of new coupled repressilators in genetic regulatory networks with
delays. Neural Netw. 2014;60:222-31.
14. Muller S, Hofbauer J, Endler L, Flamm C, Widder S, Schuster P. A generalized
model of the repressilator. J Math Biol. 2006;53(6):905-37.
15. O'Brien EL, Itallie EV, Bennett MR. Modeling synthetic gene oscillators.
Math Biosci. 2012;236(1):1-15.
16. Potvin-Trottier. L, Lord1. ND, Vinnicombe. G, Paulsson. J. Synchronous
long-term oscillations in a synthetic gene circuit. Nature. 2016;538.
17. Purcell O, Savery NJ, Grierson CS, di Bernardo M. A comparative analysis of
synthetic genetic oscillators. J R Soc Interface. 2010;7(52):1503-24.
18. Shampine LF, Thompson S. Numerical Solution of Delay Differential Equations.
2009:1-27.
19. Philip Nelson. Biological Physics: Energy, Information, Life.
