Team:HZAU-China/Analysis

Simulation and Sensitivity Analysis

After describing the biological processes, we wanted to characterize our designed processing modules by simulation. Before the simulation, we integrated some information about non-induced promoter strength we got from the wet lab. We found that some promoter properties did influence the cell's state directly. To take the intrinsic noise into consideration, we simulated the stochastic time course trajectories of the state of a chemical reaction network using Gillespie algorithm (Gillespie, 1977). We analyzed our two designs respectively. Before the wiki freeze deadline, we got some qualitative and quantitative results of design 2, which verified the validity of our models

4.1 Parameters

The parameters we used in our model of the two designs were mainly from two resources. Parameters in the processing module of design 1 were from the repressilator model in BioModels Database. And parameters used in design 2 came from parameter page collected by ETH Zurich iGEM 2013. And we adjusted some parameters by more reasonable consideration and some promoter tests. We adjusted the protein half life from 10min to 40min and adjusted the relative strength among the promoters we used according to our experimental results (Fig. 1).

Figure 1. The non-induced promoter strength in E. coli DH-5α.

4.2 Design 1

4.2.1 The effect of promoter strength

For those people who seldom pay attention to model, it is natural to believe that the repressilator will exhibit a stable oscillation. But it is not always May. Most of the models about the repressilator retain some assumptions made by Elowitz and Leibler (Elowitz and Leibler, 2000), including \begin{equation} \begin{aligned} \beta_{1(i)}&=\beta_{1}, K_{tl(i)}=K_{tl},\\ K_{R(i)}&=K_{R}, K_{P(i)}=K_{P}. \end{aligned} \end{equation} However the three genes are not identical. Our characterization of the promoters showed that the promoter strengths for these three genes were different. For this reason, we treated them differently by using $i=1$ to index gene cI, $i=2$ to index gene tetR, $i=3$ to index gene lacI. The promoter strength of placI that drives cI was represented by $\beta_{1(1)}$; the promoter strength of pcI that drives tetR in repressilator and drives lacI in toggle switch was represented by $\beta_{1(2)}$; the promoter strength of ptet was represented by $\beta_{1(3)}$.

Figure 2. Dynamics and phase diagrams of repressilator model in two cases. The common parameters we used were $\beta_0=0.03, K_{tl}=6.93, K_{R}=0.347, K_{P}=0.0173$. We emphasized the differences by simulating for a very long time. (A) Dynamics of repressilator model with $\beta_1=(30,30,30)$; (B) Dynamics of repressilator model with $\beta_1=(30,15,20)$; (C) Phase diagram of repressilator model with $\beta_1=(30,30,30)$. The trajectory forms a limit cycle. (D) Phase diagram of repressilator model with $\beta_1=(30,15,20)$. The trajectory converges to a fixed point. (E) Phase diagram of repressilator model with $\beta_1=(30,30,30)$ and different initial protein numbers. The behaviours are robust. (F) Phase diagram of repressilator model with $\beta_1=(30,15,20)$ and different initial protein numbers. The behaviours are sensitive to initial numbers.

It is obvious that such a difference can influence the attractor of this system when we plotted the time response diagram and phase trajectory diagram for the three proteins ($x_1, x_2, x_3$) (Fig. 2). We presented two potential attractors: a fixed point and a limit cycle. If the attractor is a fixed point, the gene expression is easy to converge to a steady state which means the oscillations die out. If the attractor is a limit cycle, the oscillations are stable. We simulated the oscillation by scanning different initial protein numbers to further illuminate the difference between these two attractors. The oscillation behaviour is sensitive to the initial state when the attractor is a fixed point. So we cannot make sure that we will observe the oscillations in this situation. But the oscillation behaviour is robust when the attractor is a limit cycle.

Another question arose: Whether the attractor results from the relative promoter strength or the absolute promoter strength. To find out the answer, we scanned a range of promoter strength and constructed a index named relative amplitude $A_r$, \begin{equation} A_r=\frac{2\cdot(max(x_1)-min(x_1))}{max(x_1)+min(x_1)}. \end{equation}

This index was calculated when the simulation time approached the infinity ($t=120000\ min$). So if the relative amplitude approached 0, we concluded that the system was attracted to a fixed point.

Figure 3. Results of parameter scanning. The left diagram and the right diagram are the same results from different point of views. The colorbar indicates the relative amplitudes of different simulation results. The blank regions reveal that there is a fixed point attractor in the parameter space. The colored regions imply that there is a limit cycle attractor in the parameter space.

The simulation results reveal that there is a critical surface between these two limit states (Fig. 3). An attracting limit cycle emerges as the absolute promoter strengths increase, which is consistent with the previous result (Strelkowa and Barahona, 2011). They also assumed that the promoter strengths for all genes are identical and constructed a lumped parameter $c$. The parameter $c$ is the strength of the repressive coupling, which is calculated by \begin{equation} c=\frac{\beta_1\cdot K_{tl}}{K_R\cdot K_P}. \end{equation}

They found that there is only a fixed point $p_m$ when $c$ is small enough and the uniform solution $p_m$ becomes unstable via a Hopf bifurcation as $c$ increases, resulting in a stable limit cycle.

Then how the promoter strength of each gene influences the final result? We tried to use average promoter strength to simplify this question. Comparing four different mean values, we found that the harmonic mean of the promoter strengths exhibited good competence in distinguishing these two cases. If the harmonic mean of these three genes' promoter strengths is large enough, a stable oscillation can be observed easily in theory. While any one of the promoter strengths is too weak, the harmonic mean will be too small to lead to a stable oscillation, because harmonic mean is easily affected by the minimal value.

Figure 4. Distributions of average promoter strength for two cases. The red areas represent the average promoter strength distribution for those systems with a fixed point attractor. The yellow areas represent the average promoter strength distribution for those systems with a limit cycle attractor. The orange areas represent the overlapping areas of these two distributions

The model for genetic toggle switch, by contrast, is more simple. It has been acknowledged that the relative promoter strength is more important to the steady states. If the promoter strengths are approximately equal, or, in other words, if they are in the same order, the topological structure of toggle switch creates two stable and one unstable steady states. If one promoter strength is considerably larger than the another one, only one stable steady state will be produced.

4.2.2 Situations after rewiring

Not only parameters but also the initial state decides the final steady state of the toggle switch. We set $\beta_1=(30,15,20)$ and constructed an index $S$ to estimate which steady state the toggle switch reaches. $S$ is calculated by \begin{equation} S=\lim_{t\to\infty}\log_{10}\frac{x_1(t)}{x_2(t)}. \end{equation}

The simulation results indicated S had two values which represent the two stable steady states. In our design, the initial state for toggle switch depends on the phase of oscillation when the DNA sequence is inverted. We found that the oscillation with a fixed point attractor that is far from the separatrix of toggle switch was more likely to result in one steady state when the fluctuation of gene expression becomes weaker and weaker and that the oscillation with a limit cycle attractor that crosses the separatrix was more likely to form bistability in population level (Fig. 5). In deterministic model, the steady state of a single cell in design 1 depends on the time when the site-specific recombination happened (Fig. 6). And the stochastic simulation predicts similar behaviors, however it is not always consistent with the deterministic simulation, especially for those cases that oscillations stop at a phase near the separatrix of toggle switch (Fig. 7).

Figure 5. Coupling in design 1. The red trajectory represents the gene expression dynamic in the phase space (lacI, cI). The dash line is a separatrix. The toggle switch with initial state above the separatrix will settle to one stable steady state, whereas a toggle starting below the separatrix will settle to another stable steady state. The left figure shows the situation with $\beta_1=(30,15,20)$. The right figure reveals the situation with $\beta_1=(60,30,40)$.

Figure 6. Deterministic simulation of design 1. The state versus time is simulated by scanning the moment of rewiring from 100 min to 900 min. The interval of the transient state we set is 20 min

Figure 7. Stochastic simulation of design 1. The parameters we used are the same as the deterministic model. The simulation is achieved by standard stochastic simulation algorithm (Direct method).

4.3 Design 2

In design 2, we also paid attention to some characteristics of parts. The pluxR promoter gives weak basal expression independent to the present of AHL. Such a leakage contributes to the formation of positive feedback. The weaker the basal expression of LuxI is, the longer the delay of positive feedback is. If the basal expressive of LuxI is too weak, AHL cannot accumulate in the cells leading to a very low steady state (Fig. 8). The delay may result in bistability when we observe the cells in population level.

Figure 8. Delay of the positive feedback loop. We simulated the positive feedback situations under the different basal expression condition. The basal expression of LuxI $beta_0$ controls the delay of positive feedback loop.

The DNA invertion can lead to a drastic change in gene expression. And the negative feedback loop structure forces the system to converge to a lower steady state. We used an explicit Tau-leaping method (Cao et al., 2006) to perform the stochastic simulation here instead of the classical stochastic simulation algorithm, since the high reaction rates made the classical algorithm ineffecient. The stochastic simulation result reveals that there is a weak oscillation in the negative feedback loop stage. And we ascribed this phenomenon to the delay of gene expression and noise in the systems. The steady states we predicted were consistent with our qualitative and quantitative results from wet lab (Fig. 9).

Figure 9. Result of design 2. In the big picture above, the dark blue line stands for the stochastic simulation result and the red dash line stands for the deterministic simulation result. The small picture is our experiment result. The two bottom pictures are the qualitative results from the fluorescence microscope.

References

Gillespie, D. T. (1977). Exact stochastic simulation of coupled chemical reactions. The journal of physical chemistry, 81(25), 2340-2361.

Elowitz, M. B., & Leibler, S. (2000). A synthetic oscillatory network of transcriptional regulators. Nature, 403(6767), 335-338.

Strelkowa, N., & Barahona, M. (2011). Transient dynamics around unstable periodic orbits in the generalized repressilator model. Chaos: An Interdisciplinary Journal of Nonlinear Science, 21(2), 023104.

Cao, Y., Gillespie, D. T., & Petzold, L. R. (2006). Efficient step size selection for the tau-leaping simulation method. The Journal of chemical physics, 124(4), 044109.