Team:HZAU-China/Analysis

After describing the biological processes and choosing a set of empirical parameters, we want to simulate our designed processing modules. Before the simulation, we characterized some promoters to estimate the promoter strength, which will sometimes influence the cell's state directly. To take the intrinsic noise into consideration, we simulate the stochastic time course trajectories of the state of a chemical reaction network using Gillespie algorithm (Gillespie, 2001). We will analysis our two designs respectively.

4.1 Design 1

4.1.1 The effect of promoter strength

It is widely believed 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 transcription rates for these three genes are different. For this reason, we treat them differently. We use $i=1$ to indexes gene cI, $i=2$ to indexes gene tetR, $i=3$ to indexes gene lacI. The promoter strength of placI which drives cI is represented by $\beta_{1(1)}$; the promoter strength of pcI which drives tetR in repressilator and drives lacI in toggle switch is represented by $\beta_{1(2)}$; the promoter strength of ptet is represented by $\beta_{1(3)}$.

It is obvious that such a difference can influence the attractor of this system when we draw the time response diagram and phase trajectory diagram for the three proteins ($x_1, x_2, x_3$) (Figure xxx). We present 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 simulate the oscillation by scanning different initial protein numbers to further illuminate the difference between these two attractors. The oscillation behaviours 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 behaviours is robust when the attractor is a limit cycle.

Another question arises: Whether the attractor results from the relative promoter strength or the absolute promoter strength. To find out the answer, we scan a range of promoter strength and construct 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 is calculated when the simulation time approaches the infinity ($t=120000\ min$). So if the relative amplitude approaches 0, we conclude that the system is attracted to a fixed point.

The simulation results reveal that there is a critical surface between these two limiting states (Figure xxx). 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}

Figure 3: Dynamics and phase diagrams of repressilator model in different cases. Parameters we used are $\beta_0=0.03, K_{tl}=6.93, K_{R}=0.347, K_{P}=0.0173$. (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 converge 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.

They found that there is a only 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 every gene influences the final result? We try to use average promoter strength to simplify this question. Comparing four different mean values, we find that the harmonic mean of the promoter strengths exhibits 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.

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.

Figure 4: 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 result. 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.

Figure 5: 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

4.1.2 Situations after rewiring

Not only parameters but also the initial state will decide the final steady state of the toggle switch. We set $\beta_1=(30,15,20)$ and construct 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 indicates S has 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 find that the oscillation with a fixed point attractor that is far from the separatrix of toggle switch is 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 is more likely to form bistability in population level (Figure 10x). In deterministic model, the steady state of a single cell in design 1 depends on the time when the site-specific recombination happened. And the stochastic simulation predict 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.

Figure 6: 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)$.

4.2 Design 2

In design 2, we also pay 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 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 (Figure 13). The delay may result in bistability when we observe the cells in population level.

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 use an explicit Tau-leaping method to perform the stochastic simulation here instead of the classical stochastic simulation algorithm, since the high reaction rates make the classical algorithm ineffecient. The stochastic simulation result reveals that there is a weak oscillation in the negative feedback loop stage. And we ascribe this phenomenon to the delay of gene expression and noise in the systems.

Figure 7: 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 8: Stochastic simulation of design 1. The parameter we used is the same as the deterministic model. The simulation is achieved by standard stochastic simulation algorithm (Direct method).