Emergence of qualitative states in synthetic circuits driven

General modeling framework for synthetic gene circuit with ultrasensitive growth feedback

Exogenous gene expression causes a metabolic burden to the host cell—leading to a reduced cell growth rate, which in turn affects the gene expression through a dilution effect. These mutual inhibitions create a double-negative feedback loop in the host-circuit system (Fig 1A). To understand the impacts of growth feedback on the synthetic gene circuits, we followed our previous work model construction [21] and used the following general ordinary differential equations (ODEs) model to describe the dynamics of the gene expression in one gene circuit.

Fig 1. Ultrasensitive cell growth feedback could confer bistability in a simple synthetic gene circuit with a constitutive promoter.

(A) Diagram of the mutual inhibition between the host cell growth and synthetic gene expression. (B) Model fitting (solid lines) to experimental data (circle) of the relation between growth rate and exogenous gene expression from different media conditions [14]. (C) Rate balance plot shows the emergence of bistability (top) based on the curves of degradation plus dilution rate versus gene expression (d(x)) with various m values (heatmap) and the curve of production versus gene expression f(x) (horizontal black dashed line). Dash-dotted black line shows the degradation rate (d₀*x) curve, and the dotted black line shows the degradation plus max dilution rate (d₀*x+kg₀*x). Solid and open circles represent stable and unstable steady states, respectively. The black arrows on the x-axis denote the directional field. Vertical dotted lines between the top and bottom panels represent the x values at the stable states. The corresponding growth rate curves versus gene expression are shown at the bottom. (D) Bifurcation diagram of the gene expression x with respect to kg₀. Solid and dashed lines represent stable and unstable steady states with m = 2, respectively. (E) Two-parameter bifurcation diagram with respect to m and kg₀. SN1 (blue line) and SN2 (orange line) are the saddle-node bifurcation points shown in (D). The enclosed green region between the saddle nodes represents the bistability. The black dashed line corresponds to the case shown in (D) with m = 2.

https://doi.org/10.1371/journal.pcbi.1010518.g001

To understand the consequence of this growth feedback on the host-circuit system, we developed a simple mathematical modeling framework,

where the f(x) represents the gene production rate, which could depend on the exogenous gene expression level (x) based on the regulations in the circuits, the d₀ is the degradation rate constant of the circuit product x_i, and the last term GR(x)*x represents the dilution rate mediated by host cell growth. It is noted that the growth rate GR(x) is a function of gene expression level x, given that the expression of the synthetic gene circuit causes a metabolic burden that slows down the host cell growth. Previously, GR(x) was formulated in a Michaelis-Menten format to represent the growth of the cell under a metabolic burden [20,21]. Here, for generality, we use a Hill function:

where kg₀ is the maximal growth rate of the host cell without the exogenous genes, J indicates the expression capacity for the exogenous gene expression that determines the threshold above which significant metabolic burden is induced. For simplicity, it was also assumed that the metabolic burden caused by exogenous gene expression affects only the growth rate. This framework is suitable for analyzing the system’s steady state in the exponential growth phase. The Hill coefficient indicates the sensitivity of metabolic burden. To quantify the sensitivity of metabolic burden, we used the Hill function to fit experimental data obtained from Ref [14]. The data shows the dependence of the growth rate on exogenous gene expression levels for cells in the exponential growth phase under three different growth conditions (Fig 1B). Interestingly, the fitted sensitivity of metabolic burden is greater than one for all the media conditions (Fig 1B), suggesting ultrasensitive growth feedback. However, how the sensitivity of the growth feedback affects synthetic gene circuits’ functions is still less understood.

The emergence of bistability in a simple synthetic gene circuit with a constitutive promoter

To understand the implication of the ultrasensitive growth feedback to the host-circuit system, we first analyzed a simple circuit with a constitutive promoter. For this system, the production rate is expected to be constant and independent of x, that is f(x) = k₀. Applying the general modeling framework gives the following model for this simple system,

First, we performed a rate balance analysis by plotting the production rate (k₀) and dilution plus degradation rate (d(x)) curves together (Fig 1C). The production rate is a constant k₀, indicated as a horizontal curve (Fig 1C, black dashed line). Under the condition without cell growth, d(x) is a linear curve (Fig 1C, black dash-dotted line). One stable steady state with high-level gene expression (Fig 1C, right solid circle) is found at the intersection point with the production curve. When the growth rate is considered but independent of gene expression, the d(x) curve is also linear but with a much greater slope (Fig 1C, black dotted line). One stable steady-state with low-level gene expression (Fig 1C, left solid circle) is found instead. However, when the growth rate depends on the gene expression level, the d(x) curve is nonlinear and constrained within two dotted and dash-dotted lines (Fig 1C). It is worthy to note that the shape of d(x) curve depends on the values of m. Fig 1C shows the d(x) curve (top) and growth rate (bottom) for a range of values of m. For m>1, the d(x) curve can exhibit local maximum and minimum, which allow the system to have three steady states at three intersections of the k₀ and d(x) curves. Two of them (solid circles) are stable, representing the ON state with a high expression but low growth rate and the OFF state with low expression and high growth rate, and the other is unstable (open circle), representing the threshold that separates the two stable states. Intuitively, if circuit gene expression is not above the threshold, the caused burden is not significantly high, thus the cell growth rate is not attenuated much and still leads to a significant dilution that further leads to a decrease of the circuit gene expression and the burden. Under this condition, we can expect that the system stabilizes in one state in which the cell grows fast while circuit gene expression is low. On the other hand, if the burden caused by the circuit is very high, the growth rate is reduced significantly, so the dilution efforts will be significantly reduced, which leads to a further increase in gene expression and metabolic burden. Under this condition, we expect that the system stabilizes in a state in which the cell grows very slow and circuit gene expression is high. This intuitive explanation is indicated by the direction arrows (Fig 1C). Therefore, a simple circuit with a constitutive promoter could be bistable with ultrasensitive growth feedback.

One condition for the emergent bistability is that the d(x) curve is nonmonotonic with both local maximum (d(x_man)) and minimum (d(x_min)). Once the d(x) curve is nonmonotonic, we also need the production rate in the range d(x_min)<k₀<d(x_max) so that the three intersection points could be found (See Note A in S1 Text). Given that , three parameters, kg₀, J, and m can change the shape of d(x) curve. We perturbed these parameters further to understand the underlying mechanism and conditions of emergent bistability. S1A Fig shows the difference between the maximum and minimum of d(x) (Δd = d(x_max)−d(x_min)) in the space of kg₀ and m, where bistability could be found in the upper-right corner. It is noted that increasing kg₀ and m increase Δd, and thus enhance the range of k₀ for the system to be bistable. That is, an increase of kg₀ or m is beneficial for bistability emergence. The bistability could be lost if kg₀ is reduced (S1B Fig, blue to orange curve), and regained with an increase of m (yellow curve), and the bistable range could be increased with an increase of kg₀ (purple curve) or further increased with an increase of both kg₀ and m (green curve). The boundary curve is (See Note A in S1 Text for detail). Interestingly, this condition is independent of parameter J. Hence, our analysis suggests that a system with relatively high values of m and kg₀ sets up the conditions for one steady-state emergence.

To further demonstrate the emergence of bistability in the system due to the growth feedback, we performed a one-parameter bifurcation analysis with respect to kg₀ (Fig 1D). Notice that the system is monostable for small and high values of kg₀ corresponding to high and low gene expression states, respectively. While for moderate values of kg₀, the system becomes bistable. For fixed kg₀ in this range, the system shows bistable in a broad range of k₀, as shown in the one-parameter bifurcation diagram (S1C Fig). Saddle node points SN1, and SN2 (Figs 1D and S1C) set the lower and higher values of kg₀ or k₀ for bistability. We further performed a two-parameter bifurcation analysis with respect to kg₀ and m, which discloses the change in the kg₀ values of the saddle nodes with respect to the m value (Fig 1E). The two saddle-node get closer for decreasing values of m. The merge of the saddle nodes implies the loss of bistability. Notice that the system can display bistability for values on m between 1 and 2. This is also true in the two-parameter bifurcation analysis with respect to k₀ and m (S1D Fig). Remarkably, our analysis suggests that an emergence of a steady-state could occur with m values estimated with experimental data. Therefore, bistability emerges with ultrasensitive growth feedback.

Tristability emergence in the self-activation switch under ultrasensitive growth feedback

We further extend the above analysis to a general self-activation gene circuit (Fig 2A). Applying the general modeling framework to a self-activation circuit gives the following model,

Here, the production rate is modeled with , where k₀ is the basal production rate, k₁ the maximum production rate, K the dissociation constant, and n the hill coefficient. For a cooperative self-activation circuit, the f(x) is ultrasensitive as a function of x, shown as a sigmoid curve in Fig 2B (black curve, n = 2).

Fig 2. Tristability Emergence in the self-activation switch under ultrasensitive growth feedback.

(A) Mutual inhibition interaction between growth rate and self-activation circuit. (B) The rate balance plot shows the emergence of tristability (top) based on the production rate curve (f(x), black line) and degradation plus dilution rate d(x) (red line, top) versus gene expression x. Solid and open circles represent stable and unstable steady states. The black arrows denote the directional field. The corresponding growth rate (black line) and dilution rate (red line) curves are shown in the bottom panel. Vertical dotted lines between the top and bottom panels represent the x values at stable states. (C) Bifurcation diagram of the gene expression x with respect to kg₀. Solid and dashed lines correspond to the stable and unstable steady states. (D) Bifurcation diagram of growth rate with respect to kg₀. The solid circles in C-D mark Saddle-Nodes SN1~4. (E) Two-parameters bifurcation diagram with respect to m and kg₀ shows the dependence of the Saddle-Nodes SN1 (blue line), SN2 (orange line), SN3 (yellow line), and SN4 (purple line) on m. The region encloses between SN1 and SN2 corresponds to the bistability from the self-activation circuit. The region between SN3 and SN4 represents the bistability from ultrasensitive growth feedback. The overlap region of the two bistable regions gives tristability.

https://doi.org/10.1371/journal.pcbi.1010518.g002

With a strong self-activation, the system could be bistable without considering the metabolic burden. Three steady-states can be found at the intersections between the production curve and the degradation curve. However, when the metabolic burden is considered, the degradation d(x) curve bends down and could intersect the function f(x) curve up to five times (Fig 2B, m = 2). Thus, the system can have five steady-states; three are stable, and two are unstable. The three stable steady-states are low gene expression (OFF-state) with a high growth rate, moderate gene expression (ON-state) with a slightly reduced growth rate, and high gene expression (super ON-state) with a low growth rate and dilution rate (Fig 2B). Hence, a bistable system from a cooperative self-activation circuit becomes tristable under high-order growth feedback.

The bifurcation diagram shows the steady-state values of the gene expression (x_ss) and growth rate (GR) as a function of the maximum growth rate kg₀ (Fig 2C–2D). It is noted that four saddle-node points (SN1~4) were found, indicating the emergence of more stable qualitative states in the system. The three stable states are shown as solid branches. The first switch (indicated by SN1 and SN2) shows a similar growth rate (GR) but a considerable difference in gene expression (x) (Fig 2C–2D), which results from the positive feedback loop of the SA gene circuit. The second switch (indicated by SN3 and SN4) displays a significant difference in growth rate and gene expression produced by the growth feedback loop. For kg₀ values lower than the smallest SN point (SN1 solid blue circle) or higher than the largest SN point (SN4, solid purple circle), the system is monostable with a high expressed gene and a low growth rate or a low expression gene expression and a high growth rate, respectively. Between the two left SN points (SN1<kg₀<SN3) or two right SN points (SN2<kg₀<SN4), the system is bistable with two of the three states. Interesting, the coexistence of all the three states (i.e., tristability) is found between two middle SN points (SN3<kg₀<SN4).

We further performed a two-parameter bifurcation analysis with respect to kg₀ and m. Fig 2E illustrates the positions of the saddle-node points with respect to kg₀ by varying the m. The m−kg₀ space is divided into five regions. The white region represents the system in the monostable OFF state for low values of kg₀, or monostable ON state for low values of m and high values of kg₀, respectively. In the red region, the system is bistable, which results from the gene circuit, while in the green region, the system is bistable but from growth feedback. The overlapped region with both red and green is tristable. It is worthy to note that the tristability area decreases as the SN3 (yellow line) and SN4 (purple line) get closer until merging with decreasing values of m, similar to the system with a constitutive promoter (Fig 1E). Therefore, reducing the values of m makes the system lose one switch mediated by growth feedback. Contrastingly, SN1 (blue line) and SN2 (orange line) exist for all values of m between 1 and 2. These results suggest that the bistability enclosed of SN3 and SN4 is determined by the high-order growth feedback (Fig 2E green region), and the bistability marked by SN1 and SN2 is governed by the cooperative self-activation circuit (Fig 2E red region). In summary, the underlying mechanism of tristability lies in the coupled cooperative self-activation feedback in the circuit and the growth feedback between the circuit and the host cell.

Tristability emergence or enhancement of bistability depends on the synergism between the circuit feedback and growth feedback

To further understand the synergy between the circuit feedback and growth feedback on circuit dynamics, we performed a bifurcation analysis with varying parameter values in the self-activation circuits. We choose two key parameters that affect the self-activation circuit’s bistability capacity, the hills coefficient n representing the sensitivity of the circuit feedback, and K representing the threshold of the self-activation.

Fig 3 shows the two-parameters bifurcation with respect to kg₀ and m, for different values of n and K. It is noted that with K fixed at 0.1, decreasing the value of n shifts the SN1 and SN2 curves right until they disappear at n = 1 where the system loses the circuit bistability (first rows in Fig 3 and first/second rows in S2 Fig). In contrast, the SN3 and SN4 curves do not move with changes in n values indicating that the growth bistability is unaffected. On the other hand, with fixed n at 2 and K at 0.1, varying m value does not change SN1 and SN2 curves significantly but could collapse SN3 and SN4, leading to the loss of growth bistability (first two panels in the first row of Fig 3, and first/second columns in S2 Fig). This analysis reveals that the two switches are independently regulated by changing the circuit sensitivity and growth feedback sensitivity.

Fig 3. Tristability emerges from synergism between the self-activation circuit feedback and growth feedback.

Two-parameter bifurcation diagram with respect m and kg₀, for different values of hills coefficient n representing the sensitivity of the circuit feedback and K representing the threshold of the self-activation. Solid blue, orange, yellow, and purple lines represent saddle-node SN1~SN4, respectively.

https://doi.org/10.1371/journal.pcbi.1010518.g003

However, the interdependence between circuit and growth bistability differed with varying K values. SN3 and SN2 were lost with an increase of K values, representing the loss of the middle branch, which belongs to both circuit and growth switches (Fig 3 first/second columns, S3 Fig first/second rows). That is, the system loses tristability, but the two feedback are now synergistic to make the system bistable. Decrease of m shrinks the bistable range of this synergetic switch (second/third rows in Fig 3, second/third columns in S3 Fig). Thus, growth feedback could enhance the robustness of the circuit bistability in the region where tristability was lost. Taken together, the emergence of the tristability or enhancement of bistability depends on the synergism between the circuit feedback and growth feedback.

Tristability appears in the synthetic toggle switch circuit under ultrasensitive growth feedback

Our previous work showed how circuits interact with growth feedback depending on their network topologies [21,22]. Here, we study how the qualitative state of the toggle switch circuit is affected under the ultrasensitive growth feedback. The toggle switch circuit consists of mutual inhibition of two genes, represented by x₁ and x₂ (Fig 4A). Applying the general modeling framework to a toggle switch circuit, we built the following model,



where k_0,i is the basal expression rate for gene-i, k_1,i is the maximum inducible expression rate of gene-i, K_i is the concentration of its repressor producing half inhibition of the inducable expression rate of gene-i, n_i is the repression sensitivity of gene-i.

Fig 4. Tristability appears in the synthetic toggle switch circuit under the ultrasensitive growth feedback.

(A) Diagram of interaction between the toggle switch circuit and growth rate. (B) Nullcline analysis shows the emergence of tristability. Red and blue lines represent the nullcline (dx₁/dt = 0, dx₂/dt = 0) for x₁ and x₂, respectively. Solid and open circles represent stable and unstable steady states. Black arrows represent the directional field of the system. (C) The bifurcation diagram of x₁ with respect to kg₀. Solid and dashed lines correspond to stable and unstable steady state points. Solid circles denote the saddle-nodes (blue and orange for SN1-2, respectively) and pitchfork bifurcation point (yellow for PF). (D) Two-parameter bifurcation diagram with respect to m and kg₀ shows the dependence of the saddle-nodes (SN1-2) and pitchfork bifurcation point (PF) on m. SN1 (blue line) and SN2 (orange line) correspond to the saddle-node bifurcation points. Yellow line represents the PF point. The left red region of the blue (SN1) and orange (SN2) lines represents the bistability of the toggle switch circuit. The green region enclosed between the orange lines corresponds to the system bistability due to the growth feedback. Tristability is found where both bistable regions intersect.

https://doi.org/10.1371/journal.pcbi.1010518.g004

Fig 4B shows the nullclines (blue and red curve) and direction field (black arrows) for the toggle switch under ultrasensitive growth feedback (m = 2). Each intersection between the nullcline represents a steady-state (circles). Notice that there are five steady states, of which three are stable (solid circle), and two are unstable (non-solid circle). The stable steady states are or with a low growth rate and with a high growth rate. Hence, tristability could also appear in a toggle switch circuit under ultrasensitive growth feedback.

Fig 4C shows the bifurcation diagram with respect to the kg₀ (m = 2), where solid and dash lines symbolize a stable and unstable point, respectively. Our analysis reveals the existence of three bifurcations points (solid circle Fig 4C), two saddle-node (SN) points, and one pitchfork bifurcation (PF) point. For low values of kg₀ (<SN1), the system shows the native bistability, which results from the double-negative feedback in the circuit. However, increasing of kg₀ value enhances the dilution of circuit gene expression and weakens the feedback strength of the circuit. The native bistability from the circuit disappears at the PF point, above which the two stable states merge into one state with two genes coexpressed at moderate levels. SN1 and SN2 mark the switch that results from the growth feedback, similar to the constitutively expressing circuit (Fig 1D) and the self-activation circuit (Fig 2C). Notice that the tristability occurs for values of kg₀ between SN1 and PF.

We performed additional bifurcation analysis with various combinations of Hill coefficients in the toggle switch and growth feedback. As shown in S4 Fig, we found that increasing the value of the metabolic sensitivity (m) not only gives a larger bistable range for the toggle switch but also enhances the chances of tristability. With m = 1, we only have a narrow range of kg0 for the system to be bistable and no tristability, but the increase of m to 3 significantly increases the bistable range and shows the emerged tristability. On the other hand, with an increase of the Hill coefficient for toggle switch (n) value, both the bistable and tristability ranges increase. Fig 4D shows the two-parameter bifurcation with respect to kg₀ and m, which reveals that the saddle nodes merge for decreasing value of m, consistent with the constitutively expressing circuit (Fig 1E) and the self-activation circuit (Fig 2E). On the other hand, the kg₀ value corresponding to the PF point increases with m, resulting in an increased bistable region and the emergence of a tristability region in the top-right corner. Thus, our analysis suggests that tristability could result from a synergy between the growth feedback and the toggle switch.

The existence of pitchfork bifurcation requires a perfectly symmetric parameter setting. An asymmetric parameter setting with k_1,1 ≠ k_1,2 indeed makes the pitchfork bifurcation disappears. However, another saddle-node bifurcation (SN3) emerges, and tristability can also be seen for a range of kg₀ (S5A Fig). The bifurcation diagram in terms of the production rate of gene x₁ (k_1,1) (S5B Fig) shows a stepwise ternary switch, which also confirms the existence of tristability with k_1,1 values different from k_1,2. Hence, the emergent tristability does not require a symmetrical parameter setting.

We further tuned the circuit parameters to see whether there is a topology dependence effect under ultrasensitive growth feedback. We found that with lower mutual inhibition thresholds (K₁ and K₂ values), the system could have four stable steady states, which is different from the case of self-activation. Fig 5A shows the nullclines and directed field of this scenario, where seven steady states were found, four of which are stable, and three are unstable. While three of the stable states are the same as in the case of tristability (Fig 4A), one additional state with moderate coexpression of both genes emerges. That is, the toggle switch circuit could show quadstability under ultrasensitive growth feedback. Furthermore, we performed a one-parameter bifurcation with respect to kg₀ to determine the mechanism of quadstability (Fig 5B). Our analysis reveals the existence of two additional saddle nodes (SN3 and SN4), which squeeze the coexpressing state within the original tristable region. Therefore, a toggle switch circuit with a relatively strong mutual inhibition could be quadstable under ultrasensitive growth feedback.

Fig 5. Quadstability emerged in the toggle switch circuit under ultrasensitive growth feedback.

(A) Nullcline analysis shows the emergence of quadstability. Red and blue lines represent the nullcline for x₁ and x₂, respectively. Solid and non-solid circles correspond to stable and unstable steady states. Black arrows represent the directional field. (B) The bifurcation diagram of x₁ with respect kg₀. The solid circles reveal four saddle-nodes (blue, orange, purple, and green for SN1~4) and pitchfork (yellow, PF) bifurcation points.

https://doi.org/10.1371/journal.pcbi.1010518.g005