Stability and Robustness of Unbalanced Genetic Toggle Switches in the Presence of Scarce Resources

Abstract

While the vision of synthetic biology is to create complex genetic systems in a rational fashion, system-level behaviors are often perplexing due to the context-dependent dynamics of modules. One major source of context-dependence emerges due to the limited availability of shared resources, coupling the behavior of disconnected components. Motivated by the ubiquitous role of toggle switches in genetic circuits ranging from controlling cell fate differentiation to optimizing cellular performance, here we reveal how their fundamental dynamic properties are affected by competition for scarce resources. Combining a mechanistic model with nullcline-based stability analysis and potential landscape-based robustness analysis, we uncover not only the detrimental impacts of resource competition, but also how the unbalancedness of the switch further exacerbates them. While in general both of these factors undermine the performance of the switch (by pushing the dynamics toward monostability and increased sensitivity to noise), we also demonstrate that some of the unwanted effects can be alleviated by strategically optimized resource competition. Our results provide explicit guidelines for the context-aware rational design of toggle switches to mitigate our reliance on lengthy and expensive trial-and-error processes, and can be seamlessly integrated into the computer-aided synthesis of complex genetic systems.

1. Introduction

Living organisms synthesize a large collection of complex products relying on a vast array of intertwined processes [1]. Synthetic biology seeks to take advantage of these processes by modifying and rewiring existing connections, and via introducing novel components. This interdisciplinary field thus holds the promise of controlling cellular behavior by combining expertise from a diverse set of domains, including experimental techniques from the life sciences and quantitative tools from engineering disciplines [2,3,4]. As a result of this integrative approach, multiple industries are expected to be transformed and revolutionized, including regenerative medicine, biosensing and bioremediation, as well as sustainable manufacturing and energy production [5,6,7,8].

While synthetic biology bears many similarities to traditional engineering disciplines, designing synthetic gene circuits is often time consuming due to their context-dependent behavior [9,10,11,12,13,14], frequently leading to unexpected and perplexing phenomena [15,16,17]. Thus, the construction of even simple systems typically relies on massive DNA libraries that needs to be iteratively refined, involving high-throughput screening and testing in a lengthy and expensive process [18,19,20]. Although this library-based screening approach can prove successful for modules of modest complexity, the method rapidly becomes infeasible with increasing circuit size. Context-dependence thus poses a critical limitation in synthetic biology by undermining the modular and rational design of large-scale systems.

Unlike the above approach based on combinatorial DNA libraries, quantitative tools and computational techniques facilitate the rational forward-engineering of complex genetic circuits. This is achieved by leveraging mathematical models explicitly accounting for various sources of context-dependence, such as metabolic burden due to the limited availability of shared resources, both internal and external to organisms. As a result, synthesizing heterologous proteins can lead to growth rate reduction, and the expression of two unrelated proteins may become coupled [10,11,13,21,22,23,24,25,26,27,28,29,30]. Tackling this issue requires system-level approaches [31,32,33,34] combining a diverse set of quantitative tools [35,36,37,38,39,40,41,42,43,44,45,46,47]. This data-driven quantitative approach thus holds the promise of transforming the way complex biocircuits are designed by offering a more scalable approach [48].

Considering the central role that multistable switches play in both synthetic and natural systems with application examples ranging from regulating cell fate to autonomous control of maximizing cellular productivity [49,50,51,52,53,54,55,56,57], it is especially troubling that their behavior displays particularly strong dependence on their context [58,59]. Thus, our objective here is to reveal how tunable biophysical parameters of the toggle switch shape its fundamental properties. In particular, we seek to characterize how the interplay of competition for shared cellular resources, parameter asymmetries, and burden from the genetic context affects the stability and robustness of toggle switches [60].

Recent results illustrate that competition for shared cellular resources has profound implications regarding fundamental properties of symmetric genetic switches [61,62,63]. Here, we significantly extend these findings to the case of asymmetric toggle switches [64], thus amplifying the impact of our analysis by considering a much wider and experimentally more realistic set of circuits. Combining nullcline analysis and potential landscape-based robustness analysis, our results reveal that resource competition acts against bistability and the unbalancedness of the toggle switch further exacerbates this effect. Thus, bistability requires reduced parameter asymmetry due to resource competition, that is, the toggle switch needs to be better balanced in the presence of scarce resources [64]. Additionally, we demonstrate that both parameter asymmetries and resource competition reduce the overall robustness of metastable fixed points to noise by decreasing the potential barriers separating these equilibria. As a result, the frequency of random switching between these states increases, thus reducing the long-term reliability of the toggle switch as a memory unit. Illustrating the practical implications of our findings, we further reveal that unbalancedness significantly reduces the critical number of bistable toggle switches that can be simultaneously deployed without their collective behavior becoming monostable due to the additional resource competition they face from each other [64]. All of the above findings highlight that resource competition generally has negative impacts on the stability and robustness properties of the toggle switch. However, our analysis also reveals that by carefully adjusting resource sequestration (e.g., via the introduction of decoy sites [65]), it can also increase the balancedness of the toggle switch, an idea similar to how substrate sequestration can be leveraged for biosensor optimization [66]. As we illustrate, this is crucial for leveraging the toggle switch as a “digital comparator” in genetic optimizer modules [67].

Importantly, the results presented in this paper are underpinned by a mechanistic model capturing the scarcity of transcriptional and translational resources, leading to accurate in vivo and in vitro predictions [26,28]. By collapsing tunable microscopic model parameters (e.g., promoter strength, DNA copy number, and ribosome binding site strength) into lumped constants, we dramatically decrease complexity while preserving interpretability. Our findings thus not only provide explicit guidelines promoting modularity and increased robustness to noise, but they can also be mapped directly to experimental considerations and concrete design choices. Consequently, the stability and robustness of genetic toggle switches can be rationally adjusted by combining a wide variety of experimental tools, such as ribosome binding site and promoter engineering [20,68], introduction of decoy sites [65], and expression of heterologous proteins [25,26].

This paper is organized as follows. After briefly introducing the mathematical model of the toggle switch explicitly accounting for the limited availability of shared cellular resources, we present the quantitative framework underpinning the computational results of the paper (see Supplementary Materials for source code). Following this, we leverage these tools to reveal the role that each parameter plays in determining stability and robustness properties of the toggle switch, and how these results can be translated to explicit design guidelines for the context-aware synthesis of genetic switches.

2. Materials and Methods

Here, we first detail the mathematical model of the toggle switch in the presence of scarce transcriptional/translational resources. Following this, we introduce the main tools and techniques we leverage to uncover how competition for shared resources and parameter asymmetry shape the stability and robustness properties of genetic toggle switches.

2.1. Mathematical Model and Parameters

Comprising the repressor proteins Y and Z, the behavior of the toggle switch [60] evolves according to

$${\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}={\displaystyle \frac{{\alpha }_Y}{1+{\left(Z/K_Z\right)}^{{\theta }_Z}}}-{\gamma }_{Y}Y,{\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}={\displaystyle \frac{{\alpha }_Z}{1+{\left(Y/K_Y\right)}^{{\theta }_Y}}}-{\gamma }_{Z}Z,$$

(1)

where $K_Y$ and ${\gamma }_Y$ are the dissociation and decay rate constants of Y, respectively, with Hill coefficient ${\Theta }_Y$, and $K_Z$, ${\gamma }_Z$, and ${\Theta }_Z$ are defined similarly. In this paper, we consider ${\Theta }_Y={\Theta }_Z=2$ corresponding to the most commonly considered case of repressors binding as dimers [50,52,69,70,71,72,73], but our analysis can be easily extended to other cases as well [20,49,56,59,74,75]. Assuming that protein decay is primarily determined by cell growth, in this manuscript we consider ${\gamma }_Y,{\gamma }_Z\approx \gamma$ where $\gamma$ is the cell growth rate. Note that gratuitous protein expression can negatively affect host growth [11,21,25,76], and the extent of this effect largely depends on experimental conditions [77,78,79]. This in turn can lead to the reallocation of cellular resources [22,23,24], eventually yielding a bidirectional coupling between genetic circuits and their host, further complicating the rational analysis and design of large-scale systems. Motivated by the evidence suggesting that such effects may only be transient in the exponential phase and that they disappear after several generations of exponential growth [26,80], here, we assume that cellular growth rate is constant. In case this assumption does not hold, integrative circuit–host models [10,32,33,34,40] offer a promising avenue of inquiry to accurately predict how cell proliferation and gene expression affect one another in the above bidirectional coupling. Finally, the production rate constants ${\alpha }_Y$ and ${\alpha }_Z$ encompass all transcriptional and translational processes. For instance, considering the mechanistic model detailed in [26,28], we have that

$${\alpha }_Y={\displaystyle \frac{{\lambda }_Y^{\mathrm{TX}}{\lambda }_Y^{\mathrm{TL}}D}{{\delta }_Y{\kappa }_Y{k}_Y}},{\alpha }_Z={\displaystyle \frac{{\lambda }_Z^{\mathrm{TX}}{\lambda }_Z^{\mathrm{TL}}D}{{\delta }_Z{\kappa }_Z{k}_Z}},$$

(2)

where ${\lambda }^{\mathrm{TX}}$ and ${\lambda }^{\mathrm{TL}}$ are transcriptional and translational rate constants, respectively; $\kappa$ and k denote the dissociation constants of RNA polymerase (RNAP) and ribosomes to their targets, respectively; and D and $\delta$ stand for DNA concentration and mrNA decay rate, respectively. Typical values of these parameters are provided in

Table 1. Typical range of parameter values.

Symbol	Meaning	Typical Value	Unit	Reference
$\kappa$	RNAP dissociation constant	1	$\mathsf{\mu }$M	[22,26]
k	ribosome dissociation constant	10	$\mathsf{\mu }$M	[26,81]
K	repressor dissociation constant	0.1	nM	[82]
D	DNA concentration	100–1000	nM	[26,83]
${\lambda }^{\mathrm{TX}}$	transcriptional rate constant	100	1/h	[26,81]
${\lambda }^{\mathrm{TL}}$	translational rate constant	1000	1/h	[26,81]
$\delta$	mRNA decay rate constant	10	1/h	[84]

.

While the above model captures the dynamics of the toggle switch when transcriptional/translational resources are abundant, it fails to account for the competition phenomenon that arises when these resources are scarce [25,26,27,28]. As both repressors, as well as the genetic context of the toggle switch, rely on the same pool of resources (building blocks, energy, RNAP, ribosomes, etc.), the above coupling effects need to be modeled explicitly for predictable system-level behavior.

Accounting for the limited availability of these scarce resources, the dynamics of the toggle switch becomes

$${\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}={\displaystyle \frac{\frac{{\alpha }_Y}{1+{\left(Z/K_Z\right)}^2}}{1+\frac{{\beta }_Y}{1+{\left(Z/K_Z\right)}^2}+\frac{{\beta }_Z}{1+{\left(Y/K_Y\right)}^2}}}-\gamma Y,{\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}={\displaystyle \frac{\frac{{\alpha }_Z}{1+{\left(Y/K_Y\right)}^2}}{1+\frac{{\beta }_Y}{1+{\left(Z/K_Z\right)}^2}+\frac{{\beta }_Z}{1+{\left(Y/K_Y\right)}^2}}}-\gamma Z$$

(3)

according to the works in [26,28], with the lumped constants

$${\beta }_Y={\displaystyle \frac{D_Y}{{\kappa }_Y}}\left(1+{\displaystyle \frac{{\lambda }_Y^{\mathrm{TX}}}{{\delta }_Y{k}_Y}}\right),{\beta }_Z={\displaystyle \frac{D_Z}{{\kappa }_Z}}\left(1+{\displaystyle \frac{{\lambda }_Z^{\mathrm{TX}}}{{\delta }_Z{k}_Z}}\right).$$

(4)

Note that ${\beta }_Y$ and ${\beta }_Z$ decrease the effective production rate constants of Y and Z, respectively, and this effect increases with protein production levels. Therefore, these lumped constants measure resource sequestration associated with the production of Y and Z, respectively, due to the limited availability of shared resources. For instance, ${\beta }_Y={\beta }_Z=0$ in case of abundant resources (i.e., in the absence of competition and unwanted coupling between processes responsible for the expression of Y and Z).

Finally, to simplify further analysis, introduce the dimensionless quantities

$$y={\displaystyle \frac{Y}{K_Y}},z={\displaystyle \frac{Z}{K_Z}},{\alpha }_y={\displaystyle \frac{{\alpha }_Y}{\gamma K_Y}},{\alpha }_z={\displaystyle \frac{{\alpha }_Z}{\gamma K_Z}},t^{\prime }=\gamma t,$$

together with ${\beta }_y={\beta }_Y$ and ${\beta }_z={\beta }_Z$, so that (3) becomes

$${\displaystyle \frac{\mathrm{d}y}{d{t}^{\prime }}}={\displaystyle \frac{\frac{{\alpha }_y}{1+z^2}}{1+\frac{{\beta }_y}{1+z^2}+\frac{{\beta }_z}{1+y^2}}}-y,{\displaystyle \frac{\mathrm{d}z}{d{t}^{\prime }}}={\displaystyle \frac{\frac{{\alpha }_z}{1+y^2}}{1+\frac{{\beta }_y}{1+z^2}+\frac{{\beta }_z}{1+y^2}}}-z.$$

(5)

Based on the typical parameter ranges presented in Table 1, we estimate ${\alpha }_y,{\alpha }_z\approx 1\cdots 100$ and ${\beta }_y,{\beta }_z\approx 0.1\cdots 10$ to be typical from (2) and (4). Naturally, this range can be extended by varying promoter regions, ribosome binding sites, degradation tags, etc.

2.2. Stability Analysis

Considering the model in (1), neglecting the effects of resource competition, it was shown numerically that balancedness of the toggle switch (i.e., ${\alpha }_y\approx {\alpha }_z$) is essential for bistability [60]. To characterize this crucial feature, we introduce $a={\alpha }_y/{\alpha }_z$ together with ${\alpha }_0=\sqrt{{\alpha }_y{\alpha }_z}$ measuring the mean expression strength of the repressors. With this, we write the dynamics (5) as

$$\dot{y}=f_y(y,z),\dot{z}=f_z(y,z)$$

(6)

such that ${\alpha }_y={\alpha }_0\sqrt{a}$ and ${\alpha }_z={\alpha }_0/\sqrt{a}$. In what follows, we assume that $a\ge 1$ without loss of generality (if $a\le 1$ then swapping y and z would result in $a\ge 1$).

The stability profile of (6) can be established using nullcline analysis to reveal the effects of resource competition and parameter asymmetry. Focusing first on balanced toggle switches (i.e., $a=1$), with ${\beta }_0=\sqrt{(1+{\beta }_y)(1+{\beta }_z)}$ the nullclines $0=f_y(y,z)$ and $0=f_z(y,z)$ intersect three times if ${\alpha }_0>2{\beta }_0$ and at a single point otherwise [61] (Figure 1a). Considering the Jacobian of (6) at these intersections, it was also shown in [61] that two of them correspond to stable fixed points, whereas the third one gives rise to an unstable equilibrium. As increased resource sequestration yields greater values of ${\beta }_0$, it pushes the nullclines lower in Figure 1a, eventually causing the transition from bistability (middle panel in Figure 1a) to monostability (right panel in Figure 1a).

While in the case of balanced toggle switches an equilibrium lies either on the $y=z$ or on the $yz=1$ manifolds [61] (middle panel in Figure 1b), in case of unbalanced switches this is no longer true. In particular, fixed points of (6) given by the intersections of the nullclines $0=f_y(y,z)$ and $0=f_z(y,z)$ must also satisfy the constraint

$$0=a{\displaystyle \frac{z}{1+z^2}}-{\displaystyle \frac{y}{1+y^2}},$$

(7)

as illustrated in the right panel in Figure 1b (see Appendix A for more details). As the branches of this constraint move away from each other as a increases, we expect unbalancedness to act against bistability. In particular, as we reveal in Section 3.1, while minor differences between ${\alpha }_y$ and ${\alpha }_z$ (i.e., $a\approx 1$) still yield bistable dynamics with stable fixed points $x_y$ and $x_z$ (y-dominated and z-dominated), exceeding a critical threshold eventually leads to a single stable fixed point, yielding monostability.

2.3. Robustness Analysis

To uncover how the the interplay between unbalancedness and resource sequestration affects the robustness of (6) to noise, here we study the average time trajectories spend near the metastable fixed points $x_y$ and $x_z$ before transitioning towards the other. To this end, we first extend (6) in the form of the overdamped Langevin dynamics as

$$\dot{y}=f_y(y,z)+\sqrt{2ϵ}{\xi }_y,\dot{z}=f_z(y,z)+\sqrt{2ϵ}{\xi }_z,$$

(8)

where $ϵ$ regulates the intensity of the zero-mean $\delta$-correlated Gaussian white noise $({\xi }_y,$${\xi }_z)$ [63]. It is this noise that leads to trajectories (infrequently) leaving the metastable fixed points [85,86,87,88], thus causing unwanted $x_y\leftrightarrow x_z$ transitions. The frequency of these events can be characterized by approximating the underlying dynamics with a Markov jump process [89], where transition rates are parametrized using the Eyring–Kramers formula [85,86,87,88]. Within this framework, transition rates depend on the potential barriers separating metastable fixed points. Therefore, we next define a suitable (quasi) potential.

As $\partial f_y(y,z)/\partial z\ne \partial f_z(y,z)/\partial y$, the dynamics in (6) do not correspond to a gradient system, hence a quasi-potential must be defined [90,91]. Following one of the most common approaches, the quasi-potential $V(y,z)$ changes along trajectories according to

$$\Delta V(y,z)=-\left[f_y^2(y,z)+f_z^2(y,z)\right]\Delta t$$

(9)

for a sufficiently small time step $\Delta t$ [90,92]. As $\Delta V(y,z)\le 0$ and $\Delta V(y,z)=0$ only when $f_y(y,z)=f_z(y,z)=0$, the potential surface behaves like a Lyapunov function [90,92]: system trajectories “flow downhill” towards the metastable fixed points $x_y$ and $x_z$ (for more details on the computation of $V(y,z)$, see Appendix B).

According to the Eyring–Kramers formula [87,88], the average time trajectories spend near a metastable fixed point (mean transition time) is exponentially proportional to the potential barrier required for leaving its neighborhood. To characterize this time and potential barrier, let ${\Omega }_y,{\Omega }_z\in {\mathbb{R}}^2$ denote the regions of convergence of the metastable fixed points $x_y$ and $x_z$, respectively. With this, ${\tau }_y=\mathbb{E}[inf\{t>0:(y,z)\in {\Omega }_y,(y\left(0\right),z\left(0\right))\in {\Omega }_z\}]$ and ${\tau }_z=\mathbb{E}[inf\{t>0:(y,z)\in {\Omega }_z,(y\left(0\right),z\left(0\right))\in {\Omega }_y\}]$ are the average duration trajectories spend near the metastable fixed points $x_y$ and $x_z$, respectively. From the Eyring–Kramers formula, as $ϵ\to 0$ the time ${\tau }_y$ is exponentially proportional to the potential barrier $h_y$ required for leaving ${\Omega }_y$ (Figure 2a), defined as

$${\displaystyle h_y=\underset{\gamma }{inf}\underset{x^{*}\in \chi }{sup}V\left(x^{*}\right)-V\left(x_y\right),}$$

(10)

where $\chi$ denotes continuous paths leading from $x_y$ to ${\Omega }_z$ [85,87,88]. The potential barrier $h_z$ is defined similarly. Therefore, to reveal how resource sequestration and unbalancedness of the toggle switch impact its robustness to noise, in Section 3.2 we characterize how the potential barriers $h_y$ and $h_z$ separating the metastable fixed points depend on these factors.

To further understand the long-term implications of reduced robustness to noise, we consider a two-step process to model both the random transitions between the metastable states $x_y$ and $x_z$ and the doubling of cells (Figure 2b). This way it becomes possible to reveal how the population-level distribution of colonies evolve over time, and how this process is shaped by resource sequestration and unbalancedness of the toggle switch. To this end, let $p_y$ and $p_z$ denote the probability of leaving $x_y$ and $x_z$ during one cell cycle, respectively (STEP 1 in Figure 2b). Furthermore, we assume that switching happens during growth and concludes before doubling takes place, thus cells preserve their states during the deterministic doubling (STEP 2 in Figure 2b). Starting from a single cell in generation 0, the population size after i doublings is $N_i=2^i$.

With this, we are interested in how the population-level distribution of cell fates evolve over time, depending on the state of the seed cell in generation 0. To this end, introduce the random variables $Y_i$ and $Z_i$ to denote the number of cells in the y-dominated and z-dominated states in generation i, respectively. Therefore, $q_i\left(y_i\right)=Pr\left(Y_i=y_i|Y_0=0\right)$ is the probability of observing $y_i$ cells in the former state after i doublings provided that the initial seed cell was in latter state, and similarly, let $r_i\left(y_i\right)=Pr\left(Y_i=y_i|Y_0=1\right)$. As with this we have

$$\begin{array}{cc}\hfill Pr\left(Z_i=z_i|Y_0=0\right)=& Pr\left(Y_i=2^i-z_i|Y_0=0\right)=q_i(2^i-z_i),\hfill \\ \hfill Pr\left(Z_i=z_i|Y_0=1\right)=& Pr\left(Y_i=2^i-z_i|Y_0=1\right)=r_i(2^i-z_i),\hfill \end{array}$$

(11)

the population-level composition of cell fates can be calculated by computing $q_i(·)$ and $r_i(·)$ for $i=0,1,2,\cdots$, as we detail in Section 3.2.

3. Results

Leveraging the tools and results from the previous section, here we reveal how resource competition and parameter asymmetry affect fundamental stability and robustness properties of the toggle switch. In addition to illustrating their unwanted consequences, we also uncover how these tunable biophysical properties can be utilized, for instance, by exploiting resource sequestration to balance asymmetric toggle switches so that they can be employed in genetic optimizer modules [67]. As our approach relies on a mechanistic model of the system dynamics underpinned by biophysical parameters with clear interpretations, the findings presented here are directly translatable to experimental considerations.

3.1. Stability Properties

While the stability analysis of (6) is significantly more complex in the unbalanced case, it is possible to derive sufficient conditions ensuring monostability and bistability (see Appendix A for details). In particular, (6) becomes monostable if

$$1+{\displaystyle \frac{{\beta }_0}{\sqrt{b}}}\ge {\displaystyle \frac{{\alpha }_0}{\sqrt{a}+\sqrt{a-\frac{1}{a}}}}.$$

(12)

In addition to providing explicit design guidelines, the above formula (together with its counterpart in Theorem A1 in Appendix A) also illuminates the role that the parameters ${\alpha }_0=\sqrt{{\alpha }_y{\alpha }_z}$, $a={\alpha }_y/{\alpha }_z$, ${\beta }_0=\sqrt{(1+{\beta }_y)(1+{\beta }_z)}$, and $b=(1+{\beta }_y)/(1+{\beta }_z)$ play in shaping the stability profile of unbalanced toggle switches.

For instance, we have already seen that in case of balanced toggle switches (i.e., $a=1$), the dynamics are bistable if $\alpha >2{\beta }_0$ [61], independent of the value of b (Figure 3a). In case of unbalanced switches, from (12) we expect that greater values of a push the dynamics towards monostability by decreasing the right-hand side of (12), whereas increasing ${\alpha }_0$ would have the opposite effect, confirmed in Figure 3b. Our results also reveal that the repressor with stronger expression can tolerate higher resource sequestration without losing bistability. For instance, assume that $a>1$, thus ${\alpha }_y>{\alpha }_z$ (black circles in Figure 3c). As greater values of ${\beta }_0$ increase the left-hand side of (12), we expect this change to push the dynamics towards monostability. Furthermore, keeping ${\beta }_0$ constant (representing the same overall amount of resource sequestration), while increasing b decreases the left-hand side, thus pushes the dynamics towards bistability, decreasing b has the opposite impact. These effects are confirmed and illustrated in Figure 3c: while low ${\beta }_0$ yields bistable dynamics despite the unbalancedness (first panel), increasing it evenly for both sides of the toggle switch causes a shift to monostability (second panel). Conversely, shifting the same amount of overall resource sequestration exclusively towards the side with higher production rate constant $\alpha$ preserves bistability (third panel), whereas allocating it to the other side has the opposite effect (fourth panel).

3.2. Robustness Properties

Having revealed how resource sequestration and parameter asymmetry affect the stability profile of (6), we next focus on the robustness properties of the metastable fixed points considering (8). As detailed in the previous section, according to the Eyring–Kramers formula, the mean transition time between these points is proportional to the height of the potential barriers separating them. Therefore, here we first focus on how these barriers are shaped by increased competition for shared resources and unbalancedness. To this end, consider a variety of toggle switches with different pairs of $({\alpha }_y,{\alpha }_z)$ and progressively increasing ${\alpha }_0=\sqrt{{\alpha }_y{\alpha }_z}$ (Figure 4, first panel).

To reveal the role of balancedness, assume first that ${\beta }_y={\beta }_z=0$. Note that from a stability perspective, increasing ${\alpha }_0$ pushes the dynamics towards bistability (Figure 3); thus, it is reasonable to expect that robustness to noise also increases as the dynamics lie farther from the monostable/bistable border. While this is certainly the case for balanced switches (toggle variants #1 and #19) as the potential barriers $h_y$ and $h_z$ increase with ${\alpha }_0$ (red and green dotted lines in Figure 4), the relationship in case of unbalanced toggle switches is more nuanced (toggle variants #2–#18). In particular, as ${\alpha }_0$ increases by first increasing ${\alpha }_z$ while keeping ${\alpha }_y$ constant (toggle variants #2–#10), $h_z$ is indeed increasing rapidly but at the expense of $h_y$ decreasing (red and green dotted lines in Figure 4). Therefore, while increasing ${\alpha }_0$ in this case pushes the dynamics farther away from monostability (Figure 3), only the $x_z$ metastable state becomes more robust to noise, the other’s sensitivity to noise instead increases, and a similar trend can be observed when the roles are reversed (toggle variants #10–#18).

The role of resource sequestration can be analyzed similarly. As increasing ${\beta }_0$ pushes the dynamics towards monostability (Figure 3), we expect it to have a negative effect on the robustness of metastable states to noise. This is indeed the case, but parameter asymmetries also play a key role. In particular, the results in Figure 4 illustrate that while increasing ${\beta }_y$ and ${\beta }_z$ fundamentally affect $h_y$ and $h_z$, thus the robustness of $x_y$ and $x_z$ to noise, respectively, cross effects are negligible (Figure 4). That is, increasing loading on one side renders the same side more sensitive to random switchings, but leaves the other side unaffected.

Next, we focus on how the robustness of the metastable states shape the population-level composition of colonies. To this end, note that from the Eyring–Kramers formula the mean transition time ${\tau }_y$ is proportional to $exp(h_y/ϵ)$ where $ϵ$ regulates noise intensity [85,86,87,88]. Furthermore, assuming that random $x_y\to x_z$ transitions are distributed exponentially [89] with parameter $1/{\tau }_y$ (so that the mean wait time is precisely ${\tau }_y$), and considering the doubling time $t_d=ln\left(2\right)/\gamma$ where $\gamma$ is the growth rate, we obtain that the probability of a random $x_y\to x_z$ switching between consecutive doublings is given by $p_y=1-e^{-t_d/{\tau }_y}$. This highlights that $p_y$ increases with the doubling time $t_d$ and decreases with ${\tau }_y$, thus increases with noise intensity and decreases with the potential barrier $h_y$. Similarly, we obtain that $p_z=1-e^{-t_d/{\tau }_z}$ for $x_z\to x_y$ transitions. Having uncovered how the potential barriers are shaped by the interplay between competition for shared resources and balancedness (Figure 4), we next focus on how $p_y$ and $p_z$ affect the evolution of the population-level composition of colonies by computing $q_i(·)$ and $r_i(·)$ from (11). To this end, from Figure 3 with $y_i=2y_i^{\prime }$ and $z_i=2z_i^{\prime }$ it follows that

$$\begin{array}{cc}\hfill q_i\left(y_i\right)=& {\displaystyle Pr\left(Y_i^{\prime }=y_i^{\prime }|Y_0=0\right)=\sum _{y_{i-1}=0}^{2^{i-1}}Pr\left(Y_i^{\prime }=y_i^{\prime },Y_{i-1}=y_{i-1}|Y_0=0\right)}\hfill \\ \hfill =& {\displaystyle \sum _{y_{i-1}=0}^{2^{i-1}}Pr\left(Y_i^{\prime }=y_i^{\prime }|Y_{i-1}=y_{i-1},Y_0=0\right)Pr\left(Y_{i-1}=y_{i-1}|Y_0=0\right)}\hfill \\ \hfill =& {\displaystyle \sum _{y_{i-1}=0}^{2^{i-1}}Pr\left(Y_i^{\prime }=y_i^{\prime }|Y_{i-1}=y_{i-1}\right)q_{i-1}\left(y_{i-1}\right),}\hfill \\ \hfill r_i\left(y_i\right)=& {\displaystyle Pr\left(Y_i^{\prime }=y_i^{\prime }|Y_0=1\right)=\sum _{y_{i-1}=0}^{2^{i-1}}Pr\left(Y_i^{\prime }=y_i^{\prime },Y_{i-1}=y_{i-1}|Y_0=1\right)}\hfill \\ \hfill =& {\displaystyle \sum _{y_{i-1}=0}^{2^{i-1}}Pr\left(Y_i^{\prime }=y_i^{\prime }|Y_{i-1}=y_{i-1},Y_0=1\right)Pr\left(Y_{i-1}=y_{i-1}|Y_0=1\right)}\hfill \\ \hfill =& {\displaystyle \sum _{y_{i-1}=0}^{2^{i-1}}Pr\left(Y_i^{\prime }=y_i^{\prime }|Y_{i-1}=y_{i-1}\right)r_{i-1}\left(y_{i-1}\right),}\hfill \end{array}$$

together with $Pr(Y_i^{\prime }=y_i^{\prime }|Y_{i-1}=y_{i-1})={\sum }_{n_y=0}^{y_{i-1}}\left(\genfrac{}{}{0pt}{}{y_{i-1}}{n_y}\right)p_y^{m_y}{(1-p_y)}^{n_y}\left(\genfrac{}{}{0pt}{}{z_{i-1}}{n_z}\right)p_z^{m_z}{(1-p_z)}^{n_z}$ where $z_{i-1}=2^{i-1}-y_{i-1}$, $m_y=y_{i-1}-n_y$, $m_z=y_i^{\prime }-n_y$, and $n_z=2^{i-1}+n_y-y_i^{\prime }$ from Figure 3. Note that here we used the generalized definition of the binomial coefficients such that $\left(\genfrac{}{}{0pt}{}{n}{k}\right)=0$ if $k<0$ or if $k>n$. Therefore, with the initial conditions $q_0\left(0\right)=r_0\left(1\right)=1$ and $q_0\left(1\right)=r_0\left(0\right)=0$ we can recursively compute $q_i\left(y_i\right)$ and $r_i\left(y_i\right)$.

This result reveals that depending on the transition probabilities $p_y$ and $p_z$, after only a few generations the population-level profile of steady-state distribution can fundamentally differ from the state of the initial seed cell (e.g., red in Figure 5a starting from $Y_0=1$), especially if the distribution mean at the steady state is substantially different from the initial state (Figure 5b). Importantly, the expected composition of the population at steady state does not depend on the state of the initial seed cell (Figure 5c), but the speed at which this state is reached does (Figure 5d).

3.3. Balancing via Optimized Competition

The balancedness of the toggle switch not only fundamentally affects its stability and robustness properties, as we have uncovered in this paper (Figure 3, Figure 4 and Figure 5), it is also crucial when the toggle switch is utilized as a “digital comparator” to optimize cellular performance [67]. We have already highlighted that carefully chosen resource competition can restore bistability (Figure 4c), thus here we explore how it can be leveraged to increase balancedness.

To illustrate this, first consider the case when resource competition is neglected (i.e., ${\beta }_y={\beta }_z=0$) and the toggle switch is balanced (i.e., ${\alpha }_y={\alpha }_z$). In this case, trajectories converge to $x_y$ and $x_z$ if $y\left(0\right)>z\left(0\right)$ and if $y\left(0\right)<z\left(0\right)$, respectively, that is, to the fixed point that corresponds to the dominant initial coordinate. In case of an unbalanced toggle switch (i.e., ${\alpha }_y\ne {\alpha }_z$), this is not true anymore: for instance, if ${\alpha }_y>{\alpha }_z$ then some initial conditions where $y\left(0\right)<z\left(0\right)$ will yield trajectories that converge erroneously to $x_y$ (gray region in Figure 6a). To measure this effect, for the initial condition $(y_0,z_0)$ let $(y,z)\to (y_{\infty },z_{\infty })$ as $t\to \infty$ and define $\Psi =\left\{(y_0,z_0)|(y_0-z_0)(y_{\infty }-z_{\infty })<0\right\}$, that is, the set of initial conditions where the initial and final dominant coordinates are different (Figure 6a). To measure the size of this region, define $e_{\Psi }=\frac{1}{{\alpha }_y{\alpha }_z}{\int \int }_{\Psi }\mathrm{d}A$, so that $e_{\Psi }$ characterizes the fraction of the rectangle $[0,{\alpha }_y]\times [0,{\alpha }_z]$ with incorrect initial/final state pairings (Figure 6a). In particular, $e_{\Psi }=0$ in case of balanced switches and $e_{\Psi }$ increases as the toggle gets increasingly more unbalanced. The data in Figure 6a confirm that while $e_{\Psi }$ increases with unbalancedness, it stays fairly constant for a given level of unbalancedness above a certain threshold value of ${\alpha }_0=\sqrt{{\alpha }_y{\alpha }_z}$.

Importantly, unbalancedness due to ${\alpha }_y\ne {\alpha }_z$ can be mitigated by carefully selecting ${\beta }_y$ and ${\beta }_z$, e.g., via the introduction of decoy sites [65]. This is illustrated in Figure 6b where the optimal choice of $({\beta }_y,{\beta }_z)$ significantly reduces the error $e_{\Psi }$: if $1\le a={\alpha }_y/{\alpha }_z\le 2$ the error decreases from 20% (Figure 6a) to less than 1% (Figure 6b), rendering the toggle switch almost perfectly balanced. Additionally, this optimal loading of the toggle switch also expands the range of $({\alpha }_0,a)$ pairs that yield bistable dynamics, rendering previously monostable dynamics bistable (e.g., orange circle in Figure 6).

3.4. Context Effects

Here, we reveal how competition for shared resources originating in the genetic context of the toggle switch shapes its stability and robustness properties. To this end, let ${\beta }_c$ capture the resource sequestration of modules other than the toggle switch (its context), yielding the dynamics

$$\dot{y}={\displaystyle \frac{\frac{{\alpha }_y}{1+z^2}}{1+\frac{{\beta }_y}{1+z^2}+\frac{{\beta }_z}{1+y^2}+{\beta }_c}}-y,\dot{z}={\displaystyle \frac{\frac{{\alpha }_z}{1+y^2}}{1+\frac{{\beta }_y}{1+z^2}+\frac{{\beta }_z}{1+y^2}+{\beta }_c}}-z.$$

Importantly, with ${\alpha }_w^{\prime }={\alpha }_w/(1+{\beta }_c)$ and ${\beta }_w^{\prime }={\beta }_w/(1+{\beta }_c)$ for $w\in \{y,z\}$ we obtain the dynamics in (6) with these rescaled parameters instead of the original ${\alpha }_w$ and ${\beta }_w$. Thus, the results presented in this paper can be applied in a straightforward manner even in the presence of loading from the context of the toggle switch.

To illustrate the detrimental effects of such competition, consider the collective behavior of N toggle switches, given by the dynamics

$${\dot{y}}_i={\displaystyle \frac{\frac{{\alpha }_{y,i}}{1+z_i^2}}{1+\frac{{\beta }_{y,i}}{1+z_i^2}+\frac{{\beta }_{z,i}}{1+y_i^2}+{\beta }_{c,i}}}-y_i,{\dot{z}}_i={\displaystyle \frac{\frac{{\alpha }_{z,i}}{1+y_i^2}}{1+\frac{{\beta }_{y,i}}{1+z_i^2}+\frac{{\beta }_{z,i}}{1+y_i^2}+{\beta }_{c,i}}}-z_i,$$

where ${\beta }_{c,i}={\sum }_{j\ne i}(\frac{{\beta }_{y,i}}{1+z_i^2}+\frac{{\beta }_{z,i}}{1+y_i^2})$ for $i=1,2,\cdots ,N$. While a single toggle switch alone may display a bistable stability profile, the addition of further switches decreases the separation of the two stable fixed points (Figure 1), eventually leading to the collectively monostable dynamics of individually bistable toggle switches (Figure 7a) as a result of additional loading from each other [61].

To reveal the role that parameter asymmetry and resource competition play in the above phenomenon, consider the simulation data in Figure 7b,c illustrating the critical number $N_{crit}$ of identical toggle switches (i.e., ${\alpha }_w={\alpha }_{w,i}$ and ${\beta }_w={\beta }_{w,i}$ for $w\in \{y,z\}$ and $i=1,2,\cdots ,N_{crit}$) such that one more unit would render the collective behavior monostable. These results reveal two key findings, which follow from the data in Figure 4a,b. First, $N_{crit}$ increases with ${\alpha }_0=\sqrt{{\alpha }_y{\alpha }_z}$ and decreases with ${\beta }_0=\sqrt{(1+{\beta }_y)(1+{\beta }_z)}$. Second, while asymmetry in the parameters measuring resource loading via ${\beta }_y$ and ${\beta }_z$ captured by $b=(1+{\beta }_y)/(1+{\beta }_z)$ does not have any appreciable effect, the balancedness of the toggle switch (i.e., $a={\alpha }_y/{\alpha }_z\approx 1$) is crucial, as increasing the difference between ${\alpha }_y$ and ${\alpha }_z$ significantly decreases $N_{crit}$. These results once again underscore that unbalancedness exacerbates the detrimental effects of resource competition. Naturally, as the number of toggle switches increases and approaches $N_{crit}$, additional resource competition also yields reduced robustness to noise, following directly from Figure 4.

4. Discussion

Given that genetic modules display context-dependent behavior [9,25,28], predictive and quantitative models play a fundamental role in the development of complex genetic circuits [36,37]. One major source coupling the behavior of seemingly unconnected components emerges due to the limited availability of shared cellular resources, thus introducing the “bioenergetic cost” of genes due to their existence and expression [93]. As recent experimental developments enable the precise characterization and separation of this cost into expenses that cells incur at various levels [94,95], these high-throughput technologies illuminate the inner workings of cells at the part-level. As a result, leveraging quantitative modeling and formal analytic tools offer promising avenues for aiding the rational design of synthetic gene circuits.

Therefore, in this paper we focused on revealing how competition for shared cellular resources affects the stability and robustness properties of one of the most widely used genetic modules: the toggle switch [60]. This core building block plays a central role in a vast array of both natural and synthetic systems. One prominent example is checkpoint control enabling the division of complex tasks into independent sub-tasks [96], allowing cells to respond to a wide variety of input signals [97] influencing a diverse set of cellular processes [98,99]. Given their central role, understanding how fundamental dynamic properties of the toggle switch depend on tunable biophysical parameters in real-world applications is thus of considerable interest to promote modularity and predictable system-level performance.

To this end, our combined analytical/numerical approach uncovered explicit guidelines aiding the design and tuning of genetic switches in a rational fashion, even in the presence of competition for shared transcriptional and translational resources. For instance, we revealed that while greater protein expression rates push the dynamics towards bistability and yield increased robustness to noise, resource competition has the opposite effect. Furthermore, our findings highlight that parameter asymmetries play a crucial role in establishing stability and robustness properties: not only can they exacerbate detrimental effects of resource competition, but also restore bistability and balancedness when carefully optimized.

To obtain the results presented in this paper, we explicitly modeled the scarcity of shared resources and the resulting coupling phenomena, both between genes of the toggle switch and arising due to its genetic context. The reduced order model underpinning our findings offers a realistic approximation of the dynamics of the switch as it considers parameter asymmetries and the model assumptions lead to accurate experimental predictions in vitro and in vivo [26,28]. As the parameters all possess clear physical interpretations and correspond to easily tunable properties of standard genetic parts, the design guidelines we uncovered are directly translatable to experimental considerations. For instance, $\alpha$ can be tuned via ribosome binding site engineering [68], $\beta$ via the introduction of decoy sites [65], and ${\beta }_c$ via the expression of heterologous proteins [25,26].

The results presented here are complemented by recent efforts to mitigate the adverse effects of competition for shared cellular resources, for instance, by decoupling resource-coupled gene expression [100], upregulating ribosome production reacting to increased metabolic burden [101], and splitting up multicomponent genetic systems into smaller subcomponents distributed among multiple collaborative cell strains [102]. Additionally, as integrative models [10,32,33,34,40] can accurately capture the bidirectional coupling between genetic circuits and the host harboring them, they offer unique and invaluable insights for the synthesis of large-scale biocircuits, for instance, by revealing how cell proliferation and gene expression affect one another. Our findings together with these tools thus offer promising opportunities for the rational and context-aware design of genetic switches relying on carefully characterized parts [103]. Therefore, we expect the key findings presented in this paper to be incorporated into the computer-aided fabrication of large-scale synthetic circuits [96,104,105], among other effects of context-dependence [9,15,22,23,24].