Quorum-Sensing Synchronization of Synthetic Toggle Switches: A Design Based on Monotone Dynamical Systems Theory

Synthetic constructs in biotechnology, biocomputing, and modern gene therapy interventions are often based on plasmids or transfected circuits which implement some form of “on-off” switch. For example, the expression of a protein used for therapeutic purposes might be triggered by the recognition of a specific combination of inducers (e.g., antigens), and memory of this event should be maintained across a cell population until a specific stimulus commands a coordinated shut-off. The robustness of such a design is hampered by molecular (“intrinsic”) or environmental (“extrinsic”) noise, which may lead to spontaneous changes of state in a subset of the population and is reflected in the bimodality of protein expression, as measured for example using flow cytometry. In this context, a “majority-vote” correction circuit, which brings deviant cells back into the required state, is highly desirable, and quorum-sensing has been suggested as a way for cells to broadcast their states to the population as a whole so as to facilitate consensus. In this paper, we propose what we believe is the first such a design that has mathematically guaranteed properties of stability and auto-correction under certain conditions. Our approach is guided by concepts and theory from the field of “monotone” dynamical systems developed by M. Hirsch, H. Smith, and others. We benchmark our design by comparing it to an existing design which has been the subject of experimental and theoretical studies, illustrating its superiority in stability and self-correction of synchronization errors. Our stability analysis, based on dynamical systems theory, guarantees global convergence to steady states, ruling out unpredictable (“chaotic”) behaviors and even sustained oscillations in the limit of convergence. These results are valid no matter what are the values of parameters, and are based only on the wiring diagram. The theory is complemented by extensive computational bifurcation analysis, performed for a biochemically-detailed and biologically-relevant model that we developed. Another novel feature of our approach is that our theorems on exponential stability of steady states for homogeneous or mixed populations are valid independently of the number N of cells in the population, which is usually very large (N ≫ 1) and unknown. We prove that the exponential stability depends on relative proportions of each type of state only. While monotone systems theory has been used previously for systems biology analysis, the current work illustrates its power for synthetic biology design, and thus has wider significance well beyond the application to the important problem of coordination of toggle switches.


Introduction
In the short period since the pioneering milestones in synthetic biology [1,2], outstanding progress has been made, and considerable practical experience has accumulated, in the construction of genetic circuits that perform various tasks, such as memory storage and logical operations, as well as support biomedical interventions and biotechnological manipulations in response to both exogenous and endogenous stimuli [3][4][5][6][7][8][9]. These circuits often include plasmids or transfected circuits which implement some form of "on-off" binary device, generically referred to as a toggle switch. For example, the expression of a protein used for gene therapy could be triggered by the recognition of some combination of inducers such as antigens, and memory of this event should be maintained across a cell population until a specific stimulus commands a coordinated shut-off [1,3,4]. In this context, as well as in many others, it is desirable for populations of cells to achieve coordinated static and/or dynamic functionalities. However, this coordination is hampered by molecular ("intrinsic") or environmental ("extrinsic") noise, which may lead to spontaneous changes of state in a subset of the population and is reflected in the bimodality of protein expression, as measured for example using flow cytometry.
To achieve robustness across a population, one may implement a "majority-vote" correction circuit that brings deviant cells back into the desired state. Much synthetic biology research focuses on single-cell microorganisms, often bacteria [4,[6][7][8][9]. Bacterial populations are relatively simple, and some aspects of their complex sociality can be rationally understood [10], providing a foundation for building more complex cellular systems. For bacteria, quorumsensing (QS) has been suggested as a way for cells to broadcast their states to the population as a whole so as to facilitate consensus. QS signaling pathways [11] can, for example, regulate bacterial gene expression in response to fluctuations in cell-population density. Bacteria produce and release various signaling molecules called autoinducers (AIs) [11][12][13][14]. The detection of a minimal threshold stimulatory concentration of an AI leads to an alteration in the host's gene expression. Both Gram-positive and Gram-negative bacteria use QS communication to regulate a diverse array of physiological activities. Synthetic biology design has adopted QS communication in its toolbox [15], because natural and artificially engineered QS modules can be used to interface synthetic circuits with exogenous and endogenous cues [4], and a systematic modular approach to standardize engineering toggle genetic circuits that would allow programmed cells to be designed for various specific purposes and to communicate their states to other cells was suggested as a bioengineering "plug-and-play" modular approach [4]. The design of such QS-toggle combinations is the focus of this paper.

A Known Design and Its Drawbacks
We benchmark our design by comparing it to Toggle B2 [4], see Fig SI-1.1 in S1 Text. Despite remarkable properties of design B2, observed experimentally in controllable experimental settings [4], and studied theoretically [16,17], the fact that their functional repertoire includes not only a bistable long-term memory but also the generation of stable oscillations suggests that the environment-toggle system must be tightly controlled in order to avoid spontaneous switching, not merely between different expression states but even between different functions.
To address this challenge, we propose a novel design, which is endowed with mathematically guaranteed properties of stability and auto-correction. Our approach is closely guided by concepts and theory from the powerful framework of monotone dynamical systems pioneered by M. Hirsch and H. Smith [18][19][20][21][22][23].
We employ monotone theory to provide guarantees of global convergence to steady states, thus ruling out unpredictable ("chaotic") behaviors and sustained oscillations. These theorems are valid no matter for all values of parameters and are based only on the network structure. We also provide an extensive computational bifurcation analysis of the corresponding biochemically-detailed and biologically-relevant mathematical models. Our results for homogeneous or mixed populations are valid independently of the number of cells in the population (N ) 1), and depend only on the relative proportions of each type of state.

The Components
As a basic design, we chose a genetic toggle switch consisting of two mutually repressing genes, lacI and tetR [1]. We use two acylated homoserine lactones (Acyl-HSLs), (i) N-butanoyl-lhomoserine lactone (C4-HSL) secreted by Pseudomonas aeruginosa [24], and (ii) N-(3-hydroxy-7-cis-tetradecenoyl)-L-homoserine lactone (3-OH-C14-HSL) produced by Rhizobium leguminosarum [13] as a means of coordinating toggle-host activity. Our design has two QS arms built-in the toggle in such a way that each promoter-repressor pair is controlled by its own QS signaling pathway symmetrically. Because of this "mirror-like" toggle symmetry, we call our design a symmetric toggle or an "S" design.
To benchmark the new S toggle design and the monotone systems approach, we compare the S design to the well-studied asymmetric B2-strain (Fig SI-1.1 in S1 Text) which has one QS arm only [4,16]. In this work, we call the asymmetric B2-strain the "A" design. Our S design cannot be reduced to the A design by removing one QS arm, and, thus, the S design cannot be viewed as a straightforward extension of the A design. From a theoretical standpoint, it is worth remarking that the A design is non-monotone.
Under certain experimentally controllable conditions (e.g, unsaturated levels of AAA+ proteases ClpXP, etc.), the S vs. A toggle comparative results obtained in this work can be summarized as follows: • The S toggle design completely excludes any unpredictable chaotic behaviors, as well as stable oscillations. Typical trajectories converge globally to stable equilibria. This conclusion is valid for all parameter values, and provides a strong theoretical guarantee missing from other synthetic biology designs.
• We refer to mixed states leading to bimodal distributions as spontaneous synchronization errors. We find that the S toggle design is able to self-correct (or, auto-correct) synchronization errors, while the non-monotone A toggle design is not.
• We show how monotone systems theory can predict not only the dynamics of an S toggle population but it also explains certain monotonically increasing or decreasing parametric dependencies of population steady states. Some of these predictions can facilitate self-synchronization and, thereby, reduce the chance for synchronization errors to emerge spontaneously.

Organization of Paper
In Models, the S toggle and A toggle mathematical models are introduced. The basic formalism and fundamental mathematical results of monotone systems theory, including Strong Monotonicity and Hirsch's Theorem [18][19][20][21]23] are as well reviewed. Reference values of dimensionless parameters, scaling, and selection and interpretation of bifurcation parameters are discussed. We also formalize a concept of spontaneous synchronization errors by considering three types of equilibrium populations: One homogeneous population, and two heterogeneous (mixed) populations (bimodal distributions) with both equally (1:1)-mixed and not-equally (N 1 :N 2 )-mixed transcriptional signatures, N 2 ( N 1 , the latter giving rise to spontaneous synchronization errors, where N = N 1 + N 2 , and N is the number of cells in the given population. In Results and Discussion, we proceed to a comparative theoretical and computational analysis of the S and A toggle designs. We begin this section with results on the application of monotone systems theory to the S design, as these results constitute the main conceptual and practical subject motivating this work (Application of Monotone Systems Theory to the S Design). We then explain how monotone systems theory allows one to predict, based on qualitative knowledge only, that generically all solutions converge to equilibria, with no possible oscillations [16] nor chaotic behavior [25], no matter what the kinetic parameters are. This is in contrast to the A design, which may admit oscillations [16]. Next, we analyze single S and A toggles decoupled from the environment (Bistability in Single S and A Toggles), and observe that the S toggle remains bistable even if "redundant" repressor genes are removed from the corresponding plasmids. To show how the S design is more robust than the A design, we carry out a comparative bifurcation analysis of populations consisting of coupled S or A toggles. We select a free (bifurcation) dimensionless parameter which can be interpreted in terms of experimental interventions [6] leading to (a) changes in the membrane permeability, or (b) changes in the half-lives of repressor proteins, or (c) changes in the specific growth rate of the host cell. We finally test the toggle design capabilities to self-correct spontaneous synchronization errors by sampling the basin of attractor of the corresponding equilibrium solutions. We find that the S toggle design is able to self-correct synchronization errors far better than the A toggle design.
The paper also has Supplemental Information (SI) materials S1-S8 which can be found in S1 text (the file S1-text.pdf).
In SI-1 Toggle B2 (see S1 Text), we discuss the relationship between A design and its prototype, the E. coli strain "B2" developed by Kobayashi et al [4] who considered a number of genetic toggle switches, interfaced with a QS signaling pathway.
In SI-2 Model Derivation (see S1 Text), we derive mathematical models and carry out a nondimensionalization procedure, the conclusions of which are used in the main text (Scaling).
In SI-3 Estimation of Parameter Values (see S1 Text), we discuss ranges of biologically meaningful parameter values based on data available in the existing literature. Values of biologically meaningful parameters depend upon experimental conditions and other factors controlled by an experimenter, as reviewed in [6]. Therefore, we provide an example of a concrete estimation of the values of dimensionless parameters, which we interpret in terms of interventions reviewed in [6].
In SI-4 Alternative Definitions of Monotone Systems and Order Preservation (see S1 Text), balanced graphs, relation to graph partitions, and order presentation by flows are explained.
In SI-5 Symmetry (see S1 Text), we formalize symmetry of the S design and discuss interpretation of symmetric results with respect to nonsymmetric perturbations typical for experimental systems.
In SI-6 Exponential Stability of Cellular Populations (see S1 Text), we prove a number of general theorems to analyze exponential stability [26] of both homogeneous and heterogeneous (mixed) population equilibrium states, independently of the number N of cells in the given population, which (i.e., the value of N ! 2) can be a priori unknown.
In SI-7 Additional Figures (see S1 Text), we provide additional bifurcation diagrams.
In SI-8 Modification of the S and A Models to Describe Sequestration of AAA+ protease ClpXP (see S1 Text), additional (modified) mathematical models describing competition of ssrA-tagged protein molecules for AAA+ proteases ClpXP are described.

Models
Although our main objective in this paper is to present a conceptual and general organizing principle for the construction of self-correcting "majority-vote" multistable synthetic systems, we instantiate our ideas through a very concrete set of genes and protein products, all being standard molecular parts in synthetic biology [1,2,4,[7][8][9][27][28][29][30][31]. We do that in order to emphasize the fact that our constructs can be realistically implemented with currently available molecular components. However, replacing these components with others does not change the basic mathematical framework.

Toggle Designs
For the sake of completeness of our description, we begin our discussion of the S toggle and A toggle designs (Fig 1) with two classical orthogonal repressors (Table 1): Next, the communication network among all toggles (Fig 1) is built upon two quorum-sensing (QS) signaling molecules (Table 1): 1. N-(3-hydroxy-7-cis-tetradecenoyl)-L-homoserine lactone (3-OH-C14-HSL); 2. N-butanoyl-l-homoserine lactone (C4-HSL). Monotone-symmetric and nonmonotone-asymmetric toggle designs. S design (top panel): Activation of the expression of gene x (lacI) occurs by binding of autoinducer G (C14-HSL) to promoter P G (P cin ). Inhibition of the expression of both genes x (lacI) and u (cinI) occurs by binding of the gene product Y (TetR) of gene y (tetR) to a single promoter P Y (P tet ). Symmetrically, activation of the expression of gene y (tetR) occurs by binding of autoinducer R (C4-HSL) to promoter P R (P rhl ), while inhibition of the transcription of both genes y (tetR) and w (rhlI) occurs by binding of X (LacI) to a single promoter P X (P lac ). A design (bottom panel): Activation of the expression of gene x (lacI) occurs by binding of autoinducer R (C4-HSL) to promoter P R (P rhl ). Expression of genes y (tetR) and w (rhlI) is driven by a common single promoter P X . Gene products U and W are synthases CinI and Rhil, respectively. Gray horizontal strips correspond to integration plasmids. Genes gfp and rfp correspond to green and red florescent proteins, GFP and RFP, respectively. For the sake of brevity, the QS signaling molecules are called autoinducers G (C14-HSL) and R (C4-HSL). Note that the G-and R-signals (acylated homoserine lactones) are natural biological signals secreted by Gram-negative bacteria, including E. coli, as a means of coordinating cellular activity [4,11].
Finally, to drive the autoinducer concentrations, two synthases are used ( Table 1): 1. CinI, the gene product of cinI, driving the concentration of C14-HSL; 2. RhlI, the gene product of rhlI, driving the concentration of C4-HSL.
Using the above molecular species, we implement and study two different toggle designs called symmetric (S) and asymmetric (A) designs, respectively, (Fig 1): 1. In the S design, each of the two autoinducers activates symmetrically the transcription of the corresponding repressor gene through a single promoter, that is, promoter P cin (P G ) for gene lacI (x) and promoter P rhl (P R ) for gene tetR (y); 2. In the A design, the same repressor genes (as used in the S design) antagonistically repress one another directly, while there is only one autoinducer that asymmetrically facilitates communication between all toggles.
The genetic circuit topology used in the A design (Fig 1) is taken from [16]. In order to keep making a fair comparison with the S design, we have replace the luxR-luxI system considered in [16] by the lacI-tetR system suggested in [1], see SI-1 Toggle B2 in S1 Text. Note that both CinI and RhiI are homologous to LuxI [38].
To host the S and A toggles, we use E. coli, a bacterial cell which has been well-studied in a huge number of relevant experimental and modeling works [32,[39][40][41][42][43][44][45][46][47][48][49][50], and which has been widely used to implement and test various synthetic circuits [1,2,4,[7][8][9]15]. A practical modeling reason for this selection is narrowing-down our search for biologically-meaningful parameters to values known from the E. coli studies. However, our conclusions do not depend in any way on biological properties of the host.
As a readout of the toggle state in individual cells, we further assume that each E. coli cell contains a gene encoding a spectrally distinct fluorescent reporter, GFP for gene lacI, and RFP for gene tetR, driven by promoters that respond to the autoinducers C14-HSL and C4-HSL, respectively. We do not simulate the processes of bio-synthesis and degradation of the fluorescent proteins explicitly, using appropriate cascade models, for two reasons: (i) the "reporter" submodel does not affect the dynamics of the entire model, and (ii) the half-lives of the reporter proteins can be made similar to the half-lives of the repressor proteins [2]. cinI autoinducer gene encodes protein CinI which synthesizes C14-HSL [11,13,14] rhlI autoinducer gene encodes protein RhlI which synthesizes C4-HSL [11,12] LacI x i lactose inhibitor a DNA-binding protein encoded by lacI [1,2,33,37] TetR y i repressor protein a basic element of tetracycline-controlled regulation [1,2,34,36] CinI synthase the gene product of gene cinI [11,13,14] RhlI synthase the gene product of gene rhlI [11,12] C14-HSL g i , g e autoinducer N-(3-hydroxy-7-cis-tetradecenoyl)-L-Homoserine Lactone [13] C4-HSL r i , r e autoinducer N-butyryl-L-Homoserine Lactone [12,24]  Finally, because each toggle can either be in a state where (a) LacI protein is abundant, while TetR protein is scarce, or in a state where (b) TetR protein is abundant, while LacI protein is scarce, we call state (a) a green state or, simply, a G-state and state (b) a red state or, simply, an R-state, respectively.
S model. A mathematical model describing a population of identical S toggles is Here, all state variables and parameters are dimensionless, and are obtained from the corresponding biologically meaningful state variables and parameters describing the lac-tet system ( Table 1) after an appropriate nondimensionalization, see SI-2 Model Derivation in S1 Text.
In the S model Eqs (1)-(6), t is dimensionless time; x i and y i are the dimensionless concentrations (levels) of intracellular repressor proteins LacI and TetR, respectively; g i and r i are the dimensionless concentrations of intracellular autoinducers C14-HSL and C4-HSL, respectively; g e and r e are the dimensionless concentrations of extracellular autoinducers C14-HSL and C4-HSL, respectively.
The dimensionless rate constants a i , i = 1, . . ., 6, depend on the copy numbers of the plasmids that bear the corresponding genes, see relationships Eqs (18) and (19) given in Scaling; n X , n Y , n G , and n R are the corresponding Hill coefficients; d is the dimensionless diffusion coefficient; δ g and δ r are the dimensionless lumped dilution-degradation rates due to the exponential cell growth and degradation of the corresponding species; γ x , γ y , γ g , and γ r are the corresponding leakiness coefficients [16,17].
The degradation rate constants for repressor species x i and y i are scaled out to unity, as it is done in [1,2,16,17], see SI-3 Estimation of Parameter Values in S1 Text; δ e is the dilution rate due to flow in the medium; ρ is a population density; and N is the number of cells in the given population.
A model. A dimensionless mathematical model describing a population of identical A toggles is Here, all state variables and parameters are as defined for the S model Eqs (1)- (6). The A model is mathematically identical to the minimal (simplified) model developed in [16] for the E.coli strain (toggle) B2 [4], shown in Fig SI-1.1 in S1 Text.
Modeling assumptions. Since the intention of our work is to illustrate the application of monotone dynamical systems theory to the S design and, also, to compare the S design with other known designs, we have developed two simplified minimal models, described in S model and A model, respectively. However, because a practical implementation of synthetic toggles is still far from being a routine exercise [6], care should be taken to explain the assumptions used to construct the models.
First of all, we compare the new monotone S design with a well accepted non-monotone toggle design, the E.coli strain B2 [4], which we call the A design, and for which a substantial modeling work has been done [16,17] (SI-1 Toggle B2 in S1 Text). Therefore, for a careful comparison of these two different designs, we have accepted the corresponding modeling simplifications and assumptions [16,17]. We discuss them in SI-1 Toggle B2 and SI-2.1 Mass-Balance Equations in S1 Text.
Models always involve simplifications of reality. The impact of several such simplifications and assumptions, particularly those impacting monotonicity properties and bistabiliy regions, are: (i) a reduced promoter leakiness [2,4,51], (ii) unsaturated levels of AAA+ proteases ClpXP [9], and (iii) small values of Hill coefficients [7][8][9]16]. Because the impact of the variation in the values of the Hill coefficients and leakiness parameters on the bistability region in the A model is carefully analyzed in [16], we validate variations in these parameters for the S model only.
Other assumptions refer to the cells growth conditions, e.g., whether the cells are in solid culture, liquid culture, in a micro-chemostat, on plates, etc. To this end, the A model [16] corresponds to the strain B [4] for which a detailed description of the growth conditions can be found in the corresponding experimental protocol [4]. We can thus assume that both S and A models correspond to the same growth condition, that is, the aerobic growth in LB medium on plates [4].
In both S and A models, there is no time-dependence of dilution or degradation [16]. However, depending on the growth phase of the bacteria, the effects of dilution will not be constant [40,[47][48][49]52]. We use the stationary exponential growth assumption for the sake of simplicity [16].

Model Parameters
Uncertainty about the values of parameters characterizing molecular components of synthetic circuits always presents a significant difficulty in circuit design [2]. Here, we discuss reference values of dimensionless parameters obtained using an appropriate scaling procedure. We also explain how we select and interpret parameters for our bifurcation analysis.
Model parameters. Reference values of all parameters used in our modeling studies are estimated in SI-3 Estimation of Parameter Values in S1 Text, and these correspond to half-lives of all proteins in the range 4-10 min., which are close to a typical mRNA half-life in E. coli [2]. Also, to avoid competition for ribosomes [43], only a few plasmids bearing four promoters P X , P Y , P G , and P R can be used, and we use 1-2 copies per cell, see SI-3 Estimation of Parameter Values in S1 Text. The E. coli replication period is assumed to be around 25 min.
Despite the fact that much is known about E. coli [39-41, 45, 47, 48, 50], it is not possible to model behavior in a quantitatively precise way, since not enough is yet known about molecular interactions between the toggle and the host cell to make such a description realistic [6]. Instead, we hope to identify classes of toggle designs and dynamic behaviors to determine which of the designs could lead to an improved self-synchronization reliability and an improved capability for self-correction of spontaneous synchronization errors, when a small fraction of cells flips to the opposite (undesirable) transcriptional signature state, see Spontaneous Synchronization Errors. We will also make some predictions that might help to facilitate engineering toggles with desired robust traits.
In our computational analysis, the following set of reference parameter values is used: Groups of parameters with identically the same values are used to introduce the toggle mirror (involutive) symmetry into the S model as discussed in SI-5 Symmetry in S1 Text. We find that the working values of parameters estimated in Eqs (11)- (14) are within the range of equivalent parameters (rate constants, Hill coefficients, etc.) used earlier for genetic circuits built from similar (e.g., homologous) molecular entities [1, 2, 4, 7-9, 16, 17, 28-31, 53].
However, there is a variability in the estimation of values for some important parameters in the literature. Specifically, the values of the Hill coefficients (n G and n R ) for the binding of C4-HSL and C14-HSL to the corresponding promoters are equal to 4 (estimated ad hoc) in the model developed in [9]. On the other hand side, the values of the Hill coefficients for C6-HSL and C12-HSL are estimated in the range of values 1-2 in [7,8]. Because it may be unlikely to have high values of the Hill coefficients, and because the analytical and computational analyses of the impact of the Hill coefficients on the bistability regions for the A model have been done in great detail in [16], we have selected a compromising reference values of all Hill coefficients equal to 3 [16].
Another important parameter (γ) is the promoter leakiness [2,4,51]. To this end, very small dimensionless values of the corresponding leakiness parameters are used in the relevant modeling studies [16,17]. Because we compare our S model with the A model [16], and because the zero reference value for the leakiness parameters is used in [16], we also use the zero reference value for the leakiness parameters in both the S and A models. The tightness of the promoter control has been reported in the literature, e.g., see [2,4,12,51], and the reduction in the promoter leakiness is the subject of active ongoing research with the promise to reduce this leakiness even further dramatically [51].
To probe the robustness of the developed theory in the cases when: (i) the Hill coefficients can take on different values, (ii) the promoter leakiness can be allowed to take on nonzero values, and (iii) the S model can lose its "perfect" mirror symmetry property, we additionally analyze the S model with the following modified parameter values while all other parameter values are kept intact as in the reference set Eqs (11)- (14), Specifically, condition n G 6 ¼ n R removes the mirror symmetry from the S model. We will also use the set of modified parameters in the case when the monotonicity property is violated by the sequestration of AAA+ proteases ClpXP (Monotone Parametric Dependencies in the S Design).
Scaling. One of the goals of a model nondimensionalization and scaling is to reduce the number of (correlated) parameters by lumping original parameters into a smaller parameter set. In this case, interpretation of changes in the values of dimensionless parameters should be done carefully, as the set of dimensionless parameters is usually not in one-to-one correspondence with the set of original parameters. For example, mathematical models used for synthetically engineered systems often contain parameters representing multiple biological parts and, so, tuning a dimensionless parameter in the corresponding mathematical model can be implemented experimentally in a number of different ways [6].
The dimensional and dimensionless parameters used in the S and A toggle models are related to one another by the following relationships (see SI-2 Model Derivation in S1 Text): 1. For the dimensionless rate parameters, we obtain: 2. For the dimensionless diffusion and degradation parameters, we obtain: Let us briefly discuss Eqs (18), (19), (20) and (21). Here, the burst parameter b x for the protein X or, equivalently, LacI, depends on the efficiency of translation, controlled by strength of ribosome-binding sites (RBS) [1,6], and the mRNA half-life time [54]; [P X ] is the number of promoters per cell for gene x; k x is an average transcription rate for gene x (lacI); K X is the number of LacI proteins required to half-maximally repress P lac ; k G is the maximum production rate of C14-HSL by CinI, D G is the export rate of C14-HSL; μ is the intracellular specific dilution rate due to the host cell growth, μ = ln2/T, T is the division period. Parameters for other proteins and QS signaling molecules are defined similarly, see SI-2 Model Derivation and SI-3 Estimation of Parameter Values in S1 Text. Based on the fact that N-Acyl Homoserine Lactone Lactonase (AHL-lactonase) hydrolyzes C4-HSL effectively [55], we also assume that specific degradation rate constants for the signaling molecules, C14-HSL and C4-HSL, can be set experimentally [6], corresponding to the parameter values used in our models. We pick these specific promoters and autoinducers in order to be concrete and to justify biologically meaningful values of the model parameters. However, we wish to emphasize that our results are generic for the architectures shown in Fig 1. Bifurcation parameters. In our bifurcation analysis, we use almost all dimensionless parameters given in Eqs (11)-(14) as free parameters allowed to be varied to detect changes in stability of the corresponding solutions. Whenever a new bifurcation point is detected, we provide an appropriate interpretation in terms of interventions reviewed in [6], which can potentially lead to the corresponding effect.
For example, suppose that d is the free parameter used in bifurcation analysis. Due to the relationships Eqs (20)- (21), changes in the values of d may correspond to different and independent experimental interventions [6] leading to: (a) changes in the membrane permeability (i.e., D G and D R ), or (b) changes in the half-lives of repressor proteins, or (c) changes in the specific growth rate of the host cell. As such, both higher values of protein half-lives and diffusion permeability as well as lower values of the specific growth rate (longer replication periods) correspond to higher values of the parameter d. Recall that the value of the parameter d characterizes the strength of the interaction between cells in the given population, which facilitates selfsynchronization [15-17, 25, 56].
More broadly, we can rely upon the fact that all dimensionless parameters are defined via appropriate combinations of the original dimensional parameters Eqs (18), (19), (20) and (21) in our interpretation of results obtained from bifurcation analysis as follows.
The values of dimensionless rate parameters (i.e., a-parameters) can be changed by decreasing or increasing translational efficiency, which depends on the nucleotide sequence of the ribosome binding sites (RBS) located within the upstream noncoding part of the mRNA [1,50]. The values of dimensionless rate parameters can also be changed by decreasing or increasing the lifetime values of appropriate proteins. Indeed, a carboxy-terminal tag, based on ClpX, the ATP-dependent unfoldase/translocase of ClpXP recognizing specific protein substrates bearing ssrA tags [9,57], can be inserted at the 3W end of each repressor gene [2]. Proteases in E. coli recognize this tag and target the attached protein for destruction. Such tags are used to reduce the half-life of the proteins from more than 60 min to around 4 min, which makes it possible and (also convenient) to set the half-life times for all toggle proteins (approximately) equal to one another and close to the half-lives of mRNAs [2,15]. We assume that all ssrA tagged proteins do not compete for AAA+ protease ClpXP [9], in which case the sequestration of AAA+ protease ClpXP is negligible small and is not modeled. To this end, both RBS and carboxy-terminal tags are the principal tools by which the parameters of an engineered gene network can be adjusted experimentally [1,2,6].

Stability and Bifurcation in Cellular Populations
A number of powerful concepts and software tools have been developed to efficiently analyze bifurcations of equilibrium solutions in small-scale and medium-sized dynamical models [58][59][60][61]. To this end, however, the analysis of bifurcations in the A and S models already becomes a formidable task in terms of CPU loads at N = 10. For example, the S model describing 10 coupled toggles includes 42 ODEs. Therefore, special conceptual and computational approaches need to be developed to interpret results of modeling with A and S models for cellular populations consisting of thousands or even millions of cells.
Fortunately, due to the special structure of the Jacobian matrices for the corresponding linearizations of the A and S models, the computation of the characteristic polynomials, which are used to evaluate stability and bifurcation [26,62], can first be conceptually and, then, numerically simplified, by employing Schur's formula [63]. As a result, (i) the stability and bifurcation analyses of homogeneous populations for any N ! 2 can be rigorously reduced to the case of a population consisting of only two toggles, (ii) the analysis of a (1:1)-mixed state for any even N ! 4 can be rigorously reduced to the case of only three toggles, and (iii) the analysis of a (N 1 : N 2 )-mixed state with any N 1 6 ¼ N 2 and N 1 + N 2 = N can be rigorously reduced to the case of only four toggles as described in SI-6 Exponential Stability of Cellular Populations in S1 Text.
Schur's formula [63] also helps to solve another important nontrivial specificity of the Aand S-population models caused by multiplicity of eigenvalues due to the model's symmetry discussed in SI-6 Exponential Stability of Cellular Populations in S1 Text. Computationally, in the case of multiple eigenvalues caused by symmetry, the standard tools [58][59][60][61] cannot be used in a straightforward way, when a special care should be taken. Our theoretic developments can aid in the analysis and interpretation of all such and similar cases arising in modeling of cellular populations, see SI-6 Exponential Stability of Cellular Populations in S1 Text for more rigorous definitions and results.
Indeed, the exact (very large) number of cells, N, in a cell culture is usually unknown, as cells can die or even be washed out. In such cases, the population density parameter ρ is used, and, therefore, stability of and bifurcation in populations with respect to the variability in their densities is done. The corresponding changes in the integer parameter N that reflect changes in ρ assume a formal study of stability with respect to changes in the number of differential equations in the corresponding models. This is an ill-defined perturbation in the number of equations, and we show how it can be avoided by using the stability approach developed in SI-6 Exponential Stability of Cellular Populations in S1 Text.

Spontaneous Synchronization Errors
Capabilities of toggles to fail and recover from spontaneous synchronization errors can be formalized in terms of a multistability concept, that is, as a co-existence of bistable homogeneous populations and various heterogeneous populations (Fig 2), also called mixed states, under the same conditions. Recall that mixed states are known to lead to bistable distributions [4].
Following [4], we call a population heterogeneous or, equivalently, mixed if it comprises toggles with different transcription signatures for the same genes: (i) the repressor gene lacI is active (G-state), while tetR is repressed, and (ii) lacI is repressed, while the repressor gene tetR Quorum-Sensing Synchronization of Toggles is active (R-state), see Toggle Designs. In other words, a homogeneous population is fully characterized by either transcription signature (i) or (ii), while a heterogeneous population is characterized by mixed signatures (i) and (ii) simultaneously present in the population (Fig 2).
Different heterogeneous populations can be characterized by transcription signature "mixtures" with ratio (p:q), p + q = 1, describing the fraction of toggles in the G-state versus the fraction of toggles in the R-state within the same population. For homogeneous populations, we, therefore, have either (1:0) or (0:1) transcriptional signature (Fig 2).
With these concepts, we can formulate more precisely our objective: to find conditions under which heterogeneous (mixed) population equilibrium solutions can loose their stability or can even be eliminated completely.
As a proof of concept, an example of an (9:1)-heterogeneous population (Fig 2) will be used, where the number of toggles in the first, Green-subpopulation (G) (tetR is suppressed) is 9 times bigger that the number of toggles in the second, Red-subpopulation (R) (lacI is suppressed). In this simplest case, the G-subpopulation comprises 9 cells (p = 0.9 or 90%-fraction of all cells), while the R-subpopulation comprises one cell (q = 0.1 or 10%-fraction of all cells).
Note that our analysis of (9:1)-mixed states does not depend on the number of cells N in the entire population, which is usually unknown in experiments. In other words, our results hold for any integers N, N 1 , and N 2 , such that N = N 1 + N 2 , and N 1 : N 2 = 9: 1, where the fractions of cells with different transcription signatures are defined by the numbers p = N 1 /N and q = N 2 /N, respectively, see SI-6 Exponential Stability of Cellular Populations in S1 Text.

Monotone Systems Formalism
The systems considered here are described by the evolution of states, which are time-dependent vectors x(t) = (x 1 (t), . . ., x n (t)). The components x i represent concentrations of chemical species (such as proteins, mRNA, metabolites, and so forth), the dynamics of which are given by a system of ODE's: . . . We also write simply dx/dt = f(x), where f is a differentiable vector function with components f i . The coordinates x i (t) are non-negative numbers. We write φ(t, x 0 ) for the solution of the initial value problem _ xðtÞ ¼ f ðxðtÞÞ with x(0) = x 0 , or just x(t) if x 0 is clear from the context, and assume that this solution x(t) exists and remains bounded for all t ! 0.
Definition of monotone systems. Observe that the definition does not impose any constrains on diagonal entries @f i @x i ðxÞ. These may have arbitrary signs, even depending on x.
Monotone systems [23,64,65] were introduced by Hirsch, and constitute a class of dynamical systems for which a rich theory exists. (To be precise, we have only defined the subclass of systems that are "monotone with respect to some orthant order" but the notion of monotone dynamics can be defined with respect to more general orders.) We assume from now on that our system satisfies the following property: for each pair of distinct nodes i and j, one of these holds: Of course, there are many models for which partial derivatives may change sign depending on the particular point x. With assumptions (1-3), however, the main results that we need from monotone dynamical systems theory will be particularly easy to state.
Monotone systems cannot admit any stable oscillations [19,21,66]. Under a stronger property, described next, only convergence to steady states is generically possible.
Strong monotonicity. The directed species influence graph G associated to a system with n state variables is defined as follows. The graph G has n nodes (or "vertices"), which we denote by v 1 , . . ., v n , one node for each species. If @f i @x j > 0 ðactivationÞ; we introduce an edge labeled "1" from v j into v i . If, instead, @f i @x j < 0 ðinhibitionÞ; we introduce an edge labeled "−1" (or just "−") from v j into v i . Finally, no edge is drawn from node v j into node v i if the partial derivative @f i @x j ðxÞ vanishes identically (no direct effect of the jth species upon the ith species). An alternative is to write a normal arrow "!" or a blunted arrow "a" (or an arrow labeled "−") respectively for the first two cases. The graph G is an example of a signed graph [67], meaning that its edges are labeled by signs.
No self-edges (edges from a node v i to itself) are included in the graph G, whatever the sign of the diagonal entry @f i /@x i of the Jacobian. The sign of this derivative may be positive, negative, or even be state-dependent. Results will not depend on signs of diagonals of the Jacobian of f.
The graph G is said to be strongly connected if, given an arbitrary pair of different indices {i, j}, there is a some, possibly indirect, effect of i on j. Formally, we ask that there is a sequence of indices i = k 0 , k 1 , . . ., k r = j such that @f k sþ1 @x k s 6 ¼ 0 for s ¼ 0; . . . ; r À 1 : A system is said to be strongly monotone if it is monotone and, in addition, its species influence graph G is strongly connected. (As with the definition of monotonicity, one can extend strong monotonicity to far more general classes of systems, but we use a more restrictive notion that makes results less technical to state.) Even when there are multiple steady-states, the Hirsch Generic Convergence Theorem [21,23,64,65] is a fundamental result.
Hirsch's Theorem. Even though they may have arbitrarily large dimensionality, monotone systems behave in many ways like one-dimensional systems: Hirsch's Theorem asserts that generic bounded solutions of strongly monotone differential equation systems must converge to the set of (stable) steady states. "Generic" means here "for every solution except for a measure-zero set of initial conditions." In particular, no nontrivial attractors arise. The genericity qualifier is needed in order to exclude the unstable manifolds of saddles as well as behavior on lower-dimensional sets [18].
The general theory of monotone systems applies to a class of differential equations somewhat larger than the one considered here. What we defined as monotone systems are, to be more precise, "systems monotone with respect to an orthant order''. It is possible to, more generally, define systems that are monotone with respect to orders induced by arbitrary convex proper cones. However, the generality that one obtains in that fashion comes at the cost of conditions which are typically very difficult to verify in examples and, in any event, this generality is not needed for the purpose of analyzing the systems that are the focus of this paper.

Results and Discussion
To carry out computational bifurcation analysis, MatCont [59,68] has been used. A technical description of bifurcation points can be found in [58,59,62,68].

Application of Monotone Systems Theory to the S Design
To apply monotone systems theory to the S toggle model Eqs (1)-(6), we first rewrite the model in the following convenient general form with 4N + 2 variables: dr i dt ¼ h r ðx i ; r i ; r e Þ; dg e dt ¼ H g ðg e ; g 1 ; . . . ; g N Þ; dr e dt ¼ H r ðr e ; r 1 ; . . . ; r N Þ: Here, i = 1, . . ., N, all the functions in the right-hand side are differentiable, and the following signs hold for partial derivatives, everywhere in the state space: @h y @x i < 0; @h y @y i < 0; @h y @r i > 0; @h y @a 2 > 0; @h y @a 6 > 0; ð23Þ @h g @y i < 0; @h g @g i < 0; @h g @g e > 0; @h g @a 5 > 0; @h g @d < 0; ð24Þ @h r @x i < 0; @h r @r i < 0; @h r @r e > 0; @h g @a 6 > 0; @h r @d < 0; ð25Þ @H g @g i > 0; @H g @g e < 0; @H g @d e < 0; ð26Þ @H r @r i > 0; @H r @r e < 0; @H r @d e < 0; i ¼ 1 . . . ; N: ð27Þ Next we observe that the S system is monotone, because we may partition its state variables as follows. One set consists of Moreover, the corresponding graph is strongly connected, as we have the following paths, for each two indices i, j: x j a r j ! r e ! r i ! y i a g i ! g e ! g i ! x i ð30Þ which shows that one can reach any node from any other node by means of a directed path. Thus, the S model Eqs (1)-(6) is strongly monotone. We conclude as follows. Theorem 1 Typical solutions of the S model Eqs (1)-(6) converge to steady states. This fundamental result is robust to parameters as well as to the functional form of the equations. It ensures that our proposed design has theoretically guaranteed global stability properties. No stable oscillations [16] can exist, nor can other (for, example, "chaotic" [25]) solution regimes arise. In addition to these global properties, it is also possible to use the theory of monotone systems in order to make qualitative predictions about bifurcation diagrams as discussed in the next section.
The monotonicity property of the S system has important consequences regarding its transient as well as asymptotic behavior. We discuss in an appendix how Kamke's Theorem characterizes order relations for monotone systems. We explain now what these mean, explicitly, for the S system. Let z i (t) characterize the state of the i-th S toggle at time t ! 0, that is, z i (t) = (x i (t), y i (t), g i (t), r i (t)), i = 1, . . ., N. Let Z(t) characterize the state of the population of cells, Z(t) = (z 1 (t), . . ., z N (t), g e (t), r e (t)). Suppose that we have two initial sets, Z(0) andZð0Þ, of values for the various expression levels of the repressor proteins, LacI and TetR, and we consider the behavior of Z(t) andZðtÞ for t > 0. Now suppose that we wish to understand what is the effect of a perturbation in one of the components of the initial state z i (0) for S toggle i with some fixed i, 1 i N. (A similar argument can be applied to perturbations in other components of the initial state, or even simultaneous perturbations in all the components.) Suppose, for example, that we are interested in understanding the behavior starting from a state in whichx 3 ð0Þ ! x 3 ð0Þ in the 3rd toggle z 3 . This gives rise to a new population-wide solutionZðtÞ, and we use a tilde to denote its coordinates, that is,ZðtÞ ¼ ðz 1 ðtÞ; . . . ;z N ðtÞ;g e ðtÞ;r e ðtÞÞ, wherez i ðtÞ ¼ ðx i ðtÞ; y i ðtÞ; g i ðtÞ; r i ðtÞÞ, i = 1, . . ., N. Then, using the information provided by the partition shown in Eqs (28) and (29), we can predict that, for all t > 0:x i ðtÞ ! x i ðtÞ,ỹ i ðtÞ y i ðtÞ,g i ðtÞ ! g i ðtÞ,r i ðtÞ r i ðtÞ, g e ðtÞ ! g e ðtÞ, andr e ðtÞ r e ðtÞ for all i = 1, . . ., N. As we will see shortly below, a similar conclusion can also be made with respect to perturbations in parameters, not merely initial states.

Monotone Parametric Dependencies in the S Design
As a first step, we can include the eight parameters, a i (i = 1, . . ., 6), δ g , and δ r , as constant state variables by formally adding the corresponding equations da i /dt = 0 (i = 1, . . ., 6), and dδ g /dt = dδ r /dt = 0 to the S-model Eqs (1)- (6). The extended S-model is a monotone system. However, this extended model has no strong monotonicity property, because the nodes corresponding to the parameters cannot be reached from other nodes, as the parametric extension violates the strong connectivity relationships Eq (30). However, this is not of any consequence, as the global stability properties of the S system are determined by constant values of the parameters. We only introduced the extended model in the context of bifurcation analysis. One might add additional constant variables to represent other parameters, such as the d's. These other parameters do not lead to monotonicity, and this lack of monotonicity will have important consequences in bifurcation analysis, as we discuss later.
Dependencies between the S-model state variables and parameters Eqs (22)- (27) are shown in Fig 3 (Top Panel). Here, the set of all molecular entities in the S design is partitioned into two "orthogonal" subsets, S − and S + (Definition of monotone systems). Solid arrows and lines highlighted in light brown color correspond to S − , while solid arrows and lines highlighted in cyan color correspond to S + . Although interactions within each subset contribute to its activate state, the orthogonal subsets repress one another. Here, ClpXP is a pool of AAA+ proteases ClpXP that use the energy of ATP binding and hydrolysis to perform mechanical work during targeted protein degradation within the cell. The corresponding inhibitory (degradation) interactions are shown, using dashed gray lines. If the circuit operates near the saturation condition for the pool of AAA+ proteases ClpXP, the S design may loss its monotone properties.
However, there is a substantial body of literature that gives theorems to the effect that "small" perturbations of monotone systems retain the properties of monotone systems, for example, a smooth regular perturbation of a strongly monotone system also has generic convergence properties [65]. A similar result as well holds for singular perturbations [69].
An example of three identical S toggles interacting via common autoinducers and operating far from the saturation of AAA+ proteases ClpXP is shown in Fig 3 (Bottom Panel).
The monotonicity of the extended model implies that stable loci in bifurcation diagrams depend monotonically on parameter variations. They will increase when the parameter being varied belongs to the component as the variable being analyzed, and will decrease if they are in different components. This property is a consequence of the general order preserving properties of monotone systems, as we explain now.

Suppose that
x 0 is a steady state corresponding to a parameter value p 0 , that is to say, f ð x 0 ; p 0 Þ ¼ 0. Suppose that we now consider p 1 that is very close to p 0 and larger than p 0 , p 1 > p 0 . Suppose in addition that x 1 is a steady state for the parameter value p 1 , f ð x 1 ; p 1 Þ ¼ 0, and that x 1 is stable. Now pick the solution x 1 (t) of _ x ¼ f ðx; p 1 Þ that has initial condition x 1 ð0Þ ¼ x 0 . Suppose that the extended system _ x ¼ f ðx; pÞ and _ p ¼ 0 is monotone. Now, we may consider the following two initial states for the extended system: ð x 0 ; p 0 Þ and ð x 0 ; p 1 Þ. Since the second state is larger (in the sense of Kamke's Theorem as earlier explained) in the monotone order, it follows that the solutions satisfy x 1 ðtÞ !
x 0 for all t > 0, and therefore, taking limits, we conclude that x 1 > x 0 , as desired. Using Fig 3 in conjunction with the dimension analysis in terms of the relationships Eqs (18), (19), (20) and (21), certain qualitative predictions can be made about the parametric  Here, the blue curve connecting the origin (0, 0) and the LP 1 -point corresponds to the stable branch of the (1:1)-mixed state. The green curve connecting the origin (0, 0) with the LP 2 -point corresponds to the stable branch of the (9:1)-mixed state. Because the green curve was plotted after plotting the blue curve, a part of the blue curve is hidden beneath the green curve. Projections of the corresponding plots on 2D-planes often overlap, mixing different colors, which should not lead to any difficulty in recognizing similar monotone ("overlapping") dependencies. Panels (C) and (D). The following color coding schema is used: (i) red plots correspond to stable homogeneous G-states, (ii) violet plots correspond to stable (1:1)-mixed states, and (iii) green plots correspond to stable (9:1)-mixed states. In all cases, blue plots correspond to unstable states. All red filled circles correspond to the LP bifurcations. In the panel (D), the unstable branches for both (1:1) and (9:1)mixed states are not shown because they overlap with the stable ones. d = 0.1 (weak coupling) for all populations and, additionally, we use d = 10 (strong coupling) for the G-homogeneous population only. In the cases shown in Fig 4C and 4D, the mixed populations turn out to be more robust and exist at d = 10.
Using the S-model Eqs (1)-(6) and its sequestration version (SI-8.1) (see SI-8.1 Modification of the S Model in S1 Text) with the values of fixed parameters given in Eqs (11)- (14) and (15)- (17), respectively, we find that Fig 3 predicts monotonically increasing dependencies. The loss of stability and disappearance of the mixed states shown in Fig 4B and 4D as a 5 increases can be interpreted intuitively by the fact that an increase in a 5 leads to an increase in the intracellular levels of the corresponding QS signaling molecules, which, in turn, lead to an increase of extracellular levels of the QS molecules via diffusion, thereby facilitating self-synchronization of the given population of all toggles under conditions corresponding to a stronger interaction among all toggles. In particular, the strong interaction and coupling condition eliminates spontaneous synchronization errors in terms of suppressing the emergence of undesired (9:1)-mixed states.
This result is similar to a well-known fact for oscillators coupled via a common medium that a transition from an unsynchronized to a synchronized regime emerges as the strength of coupling increases [15-17, 25, 56]. Indeed, many microbial species accomplish this via quorum sensing, which entails the secretion and detection of diffusible molecules (autoinducers), whose concentration serves as a proxy for population density [10].
Using the expression for the dimensionless parameter a 5 given in Eqs (18) and (19), see Scaling, we can conclude that the increase in the values of the parameter a 5 leading to the bifurcation point LP 2 (Fig 4) can be achieved by the following experimental interventions: • stabilization of cell division with lower values of the specific growth rate μ (or, equivalently, higher division periods T); • stabilization of relevant proteins, using lower values of r d (or, equivalently, higher half-lives); • an increase in the maximum production rate (k G ) of C14-HSL by enzyme CinI, see SI-3 Estimation of Parameter Values in S1 Text; • an increase in the sensitivity (K G ) of promoter P cin with respect to the number of molecules C14-HSL to half-activate P cin , see Table SI-3.2 given in SI-3 Estimation of Parameter Values in S1 Text.
We have used bifurcation analysis with respect to changes in the values of the parameter a 5 as a way to illustrate predictions from monotone systems theory, and in the process we obtained conclusions regarding improvements of S toggle self-synchronization properties by eliminating the (9:1)-mixed state. To this end, we note that there is no need to further increase values of a 5 to move the system to the bifurcation point LP 1 at which the (1:1)-mixed state loses it stability and disappears, because we do not interpret the (1:1)-mixed state as a spontaneous synchronization error, see Spontaneous Synchronization Errors. Additional parametric dependencies with respect to changes in other parameters are shown in Figs SI-7.1 and SI-7.1 in S1 Text.
We then repeat the analysis of the same parametric dependencies for a (1:1)-mixed state, illustrated in Fig 5 and Fig SI-7.3 in S1 Text. Like in the previous case, we observe that all dependencies are in line with the predictions suggested by Fig 3. To this end, we will not provide here reproduced similar results for the saturated S design (SI-8.1 Modification of the S Model in S1 Text), because in all computationally investigated cases, the parameter monotone dependencies are predicted by the theory and Fig 3. The LP-bifurcation point (Fig 5) can be interpreted as follows. Decreasing values of both parameters δ g and δ r leads to an increase in the intracellular and extracellular levels of the corresponding QS signaling molecules, which, in turn, leads to stronger interactions among all toggles. Indeed, it follows from Eqs (20) and (21) (see Scaling) that the described changes in the values of dimensionless parameters δ g and δ r can be achieved by increasing half-lives of the corresponding QS signaling molecules.
To this end and similarly to the interpretation provided earlier, as the values of the parameters δ g and δ r decrease, the (1:1)-mixed state loses its stability and disappear via an LP-bifurcation (Fig 5), the effect which is similar to the well-known fact that oscillators coupled via common medium synchronize as the strength of coupling increases [15,25,56].
We note that the parametric dependencies for unstable solutions are not described by Fig 3. To explain this observation, we recall that our proof of monotone dependence on parameters applies to stable solutions only, see above.
Finally, the monotone parametric dependencies for (9:1)-mixed states corresponding to spontaneous synchronization errors are illustrated in Figs SI-7.4 and SI-7.5 in S1 Text. In this case, by increasing the strength of interactions between the toggles from the large subpopulation, the spontaneous error can also be eliminated, corresponding to the existence of the LPpoints in panels (A) and (B) of Figs SI-7.4 and SI-7.5 in S1 Text. At the same time, increasing the strength of interactions between the toggles from the small population, the corresponding spontaneous error cannot be eliminated.

Bistability in Single S and A Toggles
Before comparing population properties of our S design to those of the A design, we remark that, even for isolated cells (when the diffusion constant d is zero), there is a larger range of bistability for the S design compared to the A design. Specifically, a bistability region for a single A toggle in the plane (a 1 , a 2 ) at d = 0 is shown in Fig 6(top panel). Similar regions were found in [1,16]. We also observe that the entire quadrant, a 1 ! 0 and a 2 ! 0, spans a bistability region for the S-model at the fixed parameter values given in Eqs (11)- (14). We have computed the bistablility regions for the S design for three different nonzero values of the promoter leakiness parameter γ = 0.01, 0.1, 1.0, respectively, while all other parameter values were kept fixed as in the reference set Eqs (11)- (14), and found that in all the three cases, the entire quadrant, a 1 ! 0 and a 2 ! 0, belongs to the computed regions. In contrast, the bistability region for the A  (14). (Bottom panel). The reduced S R toggle is obtained from the original S toggle (Fig 1) after removal of genes lacI and tetR from the corresponding plasmids bearing promoters P Y and P X , respectively. This reduction procedure corresponds to setting zero values a 1 = a 2 = 0 as discussed in the main text. design depends on the promoter leakiness parameter significantly [16], and we also observed computationally that the bistability region was leaving the domain shown in Fig 6 (top panel) as soon as γ was allowed to take on values larger than 0.5, that is, when γ > 0.5.
Another important observation that follows immediately from Fig 6 (top panel) is that in the case of the S toggle, bistability persists at the origin of the non-negative quadrant in the plane (a 1 , a 2 ), that is, at a 1 = a 2 = 0. The observation remains true even for the nonzero values of the leakiness parameters as discussed earlier. Additionally, the property persists for the saturated S design (SI-8.1 Modification of the S Model in S1 Text) with the updated parameter set Eqs (15)- (17). This simply means that the genes lacI and tetR can be removed from the corresponding plasmids bearing promoters P Y and P X , respectively (Fig 1). In this case (Fig 6) (bottom panel), it is enough to keep the genes on the plasmids bearing the corresponding promoters P G and P R (Fig 1). We view the reduced S toggle as a minimal design that could be implemented experimentally. The fuller construct S is interesting too, in so far as it is based on the well-characterized and studied Cantor-Collins switch, coupled to quorum-sensing components [4]. We find that the full and reduced designs do not differ much in performance, and, so, we do not consider the minimal design in the rest of the paper.

Bistable Homogeneous Populations of S and A Toggles
Bistable homogeneous populations of S toggles persist within large ranges of the model parameters. For example, panels (A) and (B) in Fig SI-7.6 in S1 Text show scaled levels of LacI and C14-HSL, respectively, for a homogeneous population of S toggle in the G-state, depending on the values of the diffusion (membrane permeability) parameter d.
Panels (C)-(F) in Fig SI-7.6 in S1 Text show two stable homogeneous populations of A toggle which coexist while the parameter d is allowed to vary. Because the A toggle design does not have any intrinsic symmetry, the levels of the activated repressor proteins, LacI for the Ghomogeneous population shown in panels (C) and (E), and TetR for the R-homogeneous population shown panels (D) and (F), differ significantly from one another. Recall that the levels of LacI and TetR in the corresponding G-and R-homogeneous populations consisting of S toggles are identically the same due to mirror symmetry.
Our intensive computational studies confirm that the discussed results on the stable homogeneous populations of S and A toggles, as well as their dependencies on the diffusion parameter d, are robust with respect to perturbations in the model parameters, including various combinations in the values of the Hill coefficients, promoter leakiness, and the saturation conditions (SI-8 Modification of the S and A Models to Describe Sequestration of AAA+ protease ClpXP in S1 Text).
The combination of the analyses discussed here can be summarized by saying that under each one of the two designs, S and A, including biological variability in the Hill coefficients, promoter leakiness, and the degradation sequestration conditions, bistable homogeneous stable populations are possible, in either "Red" or "Green" consensus states, and with the same order of magnitude of expression. The difference between these designs, including the sequestration effect for AAA+ proteases ClpXP, become evident, when we study heterogeneous (mixed) populations, as discussed next. We see that as soon as the parameter d takes on larger values, the (1:1)-mixed state loses its stability via a Branch Point (BP) bifurcation [62] (alternatively called "pitchfork" or "symmetry-breaking" bifurcation [70,71]), giving rise to two stable (1:1)-mixed non-symmetric states at d % 1.43. The general symmetry-breaking phenomenon is rigorously studied in [72,73].

Elimination of (1:1)-Mixed Populations of S Toggles
The symmetry-breaking scenario can be described intuitively as follows. Suppose that we start with a mixed population in which 50% of the cells are in "green" state and 50% of the cells are in "red" state, and the nondimensional diffusion coefficient d (which, as we saw, in fact incorporates many of the kinetic parameters in the original system) has a low value. Suppose that we now slowly increase the value of d, and ask what happens to the (1:1)-mixed state. The first event that is observed, at d % 1.43 corresponding to the BP points in all panels of Fig 7, is that this "pure 50-50 mixed state" loses its stability. A new mixed state arises (Fig 8), in which there are two subpopulations, one in which green gene-expression dominates (but with different expression levels of LacI in each of them), and another one which red gene-expression dominates (also with different TetR levels). These two mixed states correspond to the solution branches connecting points marked with labels BP and upper LP, and BP and lower LP, respectively, shown in all panels of Fig 7. Furthermore, as d is increased a bit more (past d % 2.07 corresponding to the two points labeled with LP in all panels of Fig 7, respectively), even these mixed states disappear (Fig 7). Thus, even with moderate diffusion, heterogeneous populations cannot be sustained, emphasizing the consensus-forming character of the S design. This is in marked contrast to the A design, as shown next. The loss of stability by the (1:1)-mixed state increases the robustness of the S toggle design towards its self-synchronization by reducing the number of alternative stable states to which the toggle state can settle.

Robustness of (1:1)-Mixed Populations of A Toggles and Saturated S toggles
In contrast to (1:1)-mixed populations of (unsaturated) S toggles described by the S model Eqs (1)-(6), we observe from Fig SI-7.7 in S1 Text that the original (1:1)-mixed A-population cannot be eliminated (made unstable) by increasing the values of the parameter d within a very large parameter interval. In other words, increasing the strength of interactions between the cells does not help to establish synchronization across the given population of identical A toggles. This is in a total agreement with a similar observation reported in [16], where the A model Quorum-Sensing Synchronization of Toggles is studied in great detail. Specifically, it is found that a strong interaction between A toggles (e.g., high permeability of the membrane to the autoinducer similar to higher values of d) results in the suppression of synchronous oscillations, leading to a transition of the population to a stable heterogeneous state, where individual A toggles are locked in different equilibrium states.
Our computational experiments with (1:1)-mixed populations of (saturated) S m toggles described by the S m model (SI-8.1) (see SI-8.1 Modification of the S Model in S1 Text) led to dependencies qualitatively indistinguishable from those shown in Fig SI-7.7 in S1 Text. Therefore, we can conclude that the degradation saturation (sequestration) effect may prevent the elimination of the undesired mixed states and synchronization.
(9:1)-Mixed Population of S Toggles Next, we consider bistable (9:1)-mixed populations of S Toggles, which as discussed in the introduction, we think of as arising from random synchronization errors. We observe that (9:1)-populations of S toggles become quickly extinct as soon as the values of the nondimensional diffusion parameter d are slightly increased (Fig 9). Robustness of (9:1)-and (1:9)-Mixed Populations of A Toggles and Saturated S Toggles In contrast to the S design, in the A design, the mixed (9:1)-and (1:9)-heterogeneous populations that might arise from random state switching cannot be eliminated by changes in the values of the parameter d (Fig SI-7.8 in S1 Text). This is again in a total agreement with a similar observation reported in [16].
Using simulations carried out with the S m model (SI-8.1 Modification of the S Model in S1 Text), we also observed that the sequestration effect results in stable (9:1)-mixed states for the S design existing for large ranges of the diffusion parameter d.

Probing Capabilities of the S Toggle Design for Self-Correction of Spontaneous Synchronization Errors
To probe and compare capabilities of the S toggle and A toggle designs to correct "spontaneous synchronization errors" caused by a random flip of one toggle (or a small fraction of toggles) from a homogeneous population to the state opposite to the transcription signature adopted by the majority of the cells, we have performed simple random tests. In mathematical and computational terms, these random tests can be interpreted as an elementary numerical procedure to evaluate the size of the basin of attraction for the corresponding equilibrium solutions by sampling the corresponding small neighborhoods of the solutions, using random initial conditions, for each parameter value d 2 {0.01, 10, 100} as follows, (1) find stable G-and R-homogeneous states (for any population size!), (2) flip 10% of population, and (3) explore initial conditions in neighborhood of this state value for the corresponding state variable (for the S design, since symmetric, only the G-homogeneous state needs to be explored).
We can conclude from Fig SI-7.9 in S1 Text that the A toggle does not have any capability for self-correction of spontaneous errors for all tested values of the parameter d (Fig SI-7.9 in S1 Text). The S toggle can self-correct spontaneous synchronization errors for the medium and large values of the parameter d (Fig 10) for all parameters values for which the mixed state becomes unstable, see  first 10 minutes counted from its onset (Fig 10). Unfortunately, theory does not preclude damped oscillations. Thus, all we can do is to make computational estimates in realistic parameter ranges.
To check how the limited availability of AAA+ proteases ClpXP may negatively impact the self-correctness property by the S design, we have developed an additional S m model describing saturation (sequestration) of the AAA+ proteases ClpXP (SI-8.1 Modification of the S Model in S1 Text.) We then conducted additional computations to show that the sequestration, while preserving monotonicity and bistability properties, can lead to the loss of the self-correction of spontaneous errors by the S design. Thus, high levels of these proteases are required to implement successfully the S design.
Finally, we note that the reported results on the (9:1)-mixed states for both S and A designs are independent of the number N of cells in the given population with density ρ and can be applied to any population consisting of thousands or even millions of cells, split into two subpopulations comprising 90% and 10% fractions of all cells with different transcription signatures, respectively (Stability and Bifurcation in Cellular Populations.) Specifically, if the given (9:1)-mixed state is unstable for the S design in the model of 10 (identical) cells, it will be unstable in the model describing a larger population of (identical) cells because its stability is determined form an auxiliary system of four cells only (SI-6 Exponential Stability of Cellular Populations in S1 Text.) The same is true for the stable (9:1)-mixed state for the A design.

Conclusion
In this study, we have shown how synthetic bistable circuits (toggles), and hosting them, programmable cellular populations, can be designed so as to solve a robust molecular task, the maintenance of a coordinated state, and a "majority-vote" auto-correction of deviations, of a binary switch. Our design was guided by concepts from monotone systems theory [18][19][20][21][22][23]. Specifically, we have shown how this concept can be used for the design of a new class of monotone synthetic biological toggles, including predictive capabilities describing both dynamic state variables and monotone parametric tendencies caused by parameter perturbations.
To benchmark the new toggle design, termed the S design, and the monotone systems approach, we have compared the S design with the known (and non-monotone) B2-strain from [4], termed the asymmetric or A design in this work. The B2-strain has been previously studied both experimentally [4] and theoretically [16,17]. Despite a number of remarkable properties of the B2-strain (A design), the A toggle multifunctionality suggests that the design must be tightly controlled to avoid spontaneous switching not only between different expression states, but, as well, between different functions such as a bistable memory and an oscillatory phenotype.
In this respect, modern gene therapy interventions are currently limited to transfected genes to be either in an "on" or "off" state, when the expression of the transfected gene needs to be regulated tightly for the effective treatment of many diseases. To address this challenge, the monotone S toggle design completely excludes any unpredictable chaotic behaviors, as well as undesired stable oscillations. This conclusion is valid (of course, under certain experimentally controllable conditions pointed out in this work) for all parameter values, and provides a strong theoretical guarantee missing from other synthetic biology designs. Some of conditions include: (i) a reduced promoter leakiness [51], and (ii) unsaturated levels of AAA+ proteases ClpXP.
To achieve an in-depth understanding of dynamic properties of the S toggle design, we have developed biochemically-detailed and biologically-relevant mathematical models to test predictions of monotone systems theory by employing computational bifurcation analysis. To have all results biologically grounded, concrete molecular entities have been used, though the results are general and independent of any specific details.
To investigate the effect of a spontaneous toggle switching within cellular populations, leading to bimodal distributions, we have formalized a concept of spontaneous synchronization errors and tested the toggle design capabilities to self-correct spontaneous synchronization errors by sampling the basin of attraction of the corresponding equilibrium solutions. We found that the S toggle design was able to self-correct (or, auto-correct) synchronization errors, while the non-monotone A toggle design was not.
Because the number of cells in populations is a priori unknown, all the above results and conclusions can make sense only if they are made independently of the population size. To justify the above assertion, we have proved a few general theorems on the exponential stability of the equilibrium solutions corresponding to both homogeneous and mixed populations. The simple exponential stability results are independent of the number of cells in the populations and are based on basic first principles of stability analysis resulting from the Schur's formula [63], allowing the characteristic polynomials for the corresponding model linearizations to be computed explicitly.
Using an additional model describing saturation (and sequestration) of AAA+ proteases ClpXP (SI-8 Modification of the S and A Models to Describe Sequestration of AAA+ protease ClpXP in S1 Text), we have observed computationally that even when the above-mentioned conditions (i) and (ii) are violated, the S toggle still demonstrates the monotonicity properties. If proteases are in limited supply, however, the conclusions break down, because of the nonmonotonicity arising from resource competition. Thus, an important consideration when practically implementing our design is to express these proteases at a high enough level.
We remark that our design is based on a bistable design based on deterministic models. This approach is normally used in synthetic biology design of toggle switches, and our goal was to employ ready technology. On the other hand, in gene regulatory networks bistability, or, to be more precise, multi-modality of steady state distributions, may arise in deterministically monostable systems due to low molecule number effects. Intuitively, a slow switch between two promoter states (modeled in simplest terms by a two-state Markov chain) gives rise to a "bimodal" distribution of gene activation (gene is "on" or "off"); this process then may drive a large-molecule number mRNA and protein process, in effect creating a bimodal protein distribution, though this bimodality would be "averaged out" in a deterministic model that considers a large population. In this context, one may mention the work by Thomas et al. [74] which provides a system with two mutually repressing promoters using noncooperative transcriptional regulation but supplemented by a translational control component in which the protein product of one gene binds and degrades the mRNA of the other gene. Because we use cooperative binding (Hill coefficients 2 and larger), our design is specifically geared to bistability even in low noise situations, and the engineered consensus mechanism is designed to correct for noiseinduced transitions. It would be interesting in further work to study consensus designs for toggles based on stochastic bimodality.