Home | About | Journals | Submit | Contact Us | Français |

Formats

Article sections

Authors

Related links

eLife. 2017; 6: e23702.

doi: 10.7554/eLife.23702

PMCID: PMC5388541

Fuqing Wu,^{1}^{}^{} Ri-Qi Su,^{1,}^{2}^{}^{} Ying-Cheng Lai,^{2,}^{3,}^{4}^{}^{} and Xiao Wang^{1,}^{*}^{}^{}

Wenying Shou, Reviewing editor

Wenying Shou, Fred Hutchinson Cancer Research Center, United States;

Email: ude.usa@gnawoaix

Received 2016 November 27; Accepted 2017 March 10.

Copyright © 2017, Wu et al

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

The process of cell fate determination has been depicted intuitively as cells travelling and resting on a rugged landscape, which has been probed by various theoretical studies. However, few studies have experimentally demonstrated how underlying gene regulatory networks shape the landscape and hence orchestrate cellular decision-making in the presence of both signal and noise. Here we tested different topologies and verified a synthetic gene circuit with mutual inhibition and auto-activations to be quadrastable, which enables direct study of quadruple cell fate determination on an engineered landscape. We show that cells indeed gravitate towards local minima and signal inductions dictate cell fates through modulating the shape of the multistable landscape. Experiments, guided by model predictions, reveal that sequential inductions generate distinct cell fates by changing landscape in sequence and hence navigating cells to different final states. This work provides a synthetic biology framework to approach cell fate determination and suggests a landscape-based explanation of fixed induction sequences for targeted differentiation.

Cells in animals use a process called differentiation to specialize into specific cell types such as skin cells and liver cells. Proteins called transcription factors drive particular steps in differentiation by controlling the activity of specific genes. Many transcription factors interact with each other to form complex networks that regulate gene activity to determine the fate of a cell and control the whole differentiation process. Some individual gene networks can program cells to become any one of several different cell fates, a feature known as multistability.

In the 1950s, a scientist called Conrad Waddington proposed the concept of an “epigenetic landscape” to describe how the fate of a cell is decided as an animal develops. The cell, depicted as a ball, rolls down a rugged landscape and has the option of taking several different routes. Each route will eventually lead to a distinct cell fate. As the ball moves down the hill, the choice of routes and final destinations becomes more limited. Theoretical approaches have been used to understand how gene regulatory networks shape the epigenetic landscape of an animal. However, few studies have experimentally tested the findings of the theoretical approaches and it is not clear how environmental inputs help to determine which path a cell will take.

Although bacteria cells do not generally specialize into particular cell types, bacteria cells can use multistability in transcription factor networks to switch between different behaviors or “states” in response to cues from the environment. Wu et al. used a bacterium called *E. coli* as a model to investigate how a gene network called MINPA from mammals, which is involved in differentiation and is believed to show multistability, can guide cells to adopt different states. The work combined experimental and mathematical approaches to design, construct and test an artificial version of the MINPA gene network in *E. coli*.

The experiments showed that MINPA could direct the cells to adopt four different stable states in which the cells produced fluorescent proteins of different colors. With the help of mathematical modeling, Wu et al. charted how the landscape of cell states changed when external chemical cues were applied. Exposing the cells to several cues in particular orders guided the cells to different final states.

The findings of Wu et al. shed new light on how the fate of a cell is determined and provide a theoretical framework for understanding the complex networks that control cell differentiation. This could help develop new ways of directing cell fate that could ultimately be used to generate cells to treat human diseases.

Multistability is a mechanism that cells use to achieve a discrete number of mutually exclusive states in response to environmental inputs, such as the lysis/lysogeny switch of phage lambda (Arkin et al, 1998; Oppenheim et al., 2005) and sporulation/competence in *Bacillus subtilis* (Süel et al., 2006; Schultz et al., 2009). In multicellular organisms, multistable switches are also common in the cellular decision-making including the regulation of cell-cycle oscillator during cell mitosis (Pomerening et al., 2003), Epithelial-to-Mesenchymal transition and cancer metastasis (Jolly et al., 2016; Lee et al., 2014a), and the well-known cell differentiation process, which is a manifestation of cellular state determination in a multistable system (Laurent and Kellershohn, 1999; Guantes and Poyatos, 2008). However, loss of multistability can drive cells to acquire metastatic characteristics and stabilize highly proliferative, pathogenic cellular states in cancer (Lee et al., 2014b).

C.H. Waddington hypothesized the ‘epigenetic landscape’ to explain canalization and fate determination mechanism during cell differentiation (Waddington, 1957). In this hypothesis, differentiation is depicted as a marble rolling down a landscape with multiple bifurcating valleys and eventually settles at one of the local minima, corresponding to terminally differentiated cells. More recent theoretical studies further proposed the local minima to be modeled as steady states or attractors of dynamical systems, which can be mathematically described using differential equations (Zhang and Wolynes, 2014; Li and Wang, 2013a). As such, cell differentiation can be interpreted as a state transition process on a multistable dynamic system. A myriad of theoretical analysis have investigated the functioning of such systems and quantified the Waddington landscape and developmental paths through computation of the probability landscape for the underlying gene regulatory networks (Li and Wang, 2013a; Wang et al., 2011; Li and Wang, 2013b; Ferrell, 2012; Bhattacharya et al., 2011; Macarthur et al., 2009; Huang et al., 2007). Recent studies also revealed that the potential landscape and the corresponding curl flux are crucial for determining the robustness and global dynamics of non-equilibrium biological networks (Wang, 2015; Xu et al., 2014; Wang et al., 2008). Furthermore, the multiple stable steady states have been predicted beyond the bistable switches with or without epigenetic effects, which is reflected in slow timescales (Wang, 2015; Xu et al., 2014; Li and Wang, 2013b; Feng and Wang, 2012; Wang et al., 2011; Feng et al., 2011). Experimental researches, however, mostly focus on bistable switches, involving transitions between only two states. And demonstrations, from a combination of experiments and computational modeling, for the existence and operation of such a landscape in a higher dimensional multistable system are still lacking. Moreover, it remains unknown how gene regulatory networks (GRNs), gene expression noise, and signal induction together shape the attractor landscape and determine a cell’s developmental trajectory to its final fates (Schmiedel et al., 2015; Tanouchi et al., 2015; Prindle et al., 2014; Chalancon et al., 2012; Murphy et al., 2010; Balázsi et al., 2011; Kramer and Fussenegger, 2005; Bennett et al., 2008; Maamar et al., 2007).

Complex contextual connections of GRNs have impeded experimentally establishing the shape and function of the cell fate landscape. Rationally designed and tunable synthetic multistable gene networks in *E. coli*, however, could form well-characterized attractor landscapes to enable close experimental investigations of general principles of GRN regulated cellular state transitions. Since the functioning of these principles only requires the most fundamental aspects of gene expression regulation, they would also be applicable for cell differentiation regulations in mammalian cells. Here, we combine mathematical theory, numerical simulations, and synthetic biology to probe all possible sub-networks of mutually inhibitory network with positive autoregulations (MINPA, Figure 1A), which has been hypothesized to have multistability potentials (Guantes and Poyatos, 2008; Huang et al., 2007). Moreover, MINPA and its sub-networks are recurring motifs enriched in GRNs regulating hematopoietic development (Gata1-Pu.1, [Graf and Enver, 2009]), trophectoderm differentiation (Oct3/4-Cdx2, [Niwa et al., 2005]), endoderm formation (Gata6-Nanog, [Bessonnard et al., 2014; Li and Wang, 2013a]), and bone, cartilage, and fat differentiation (RUNX2-SOX9-PPAR-γ, [MacArthur et al., 2008; Rabajante and Babierra, 2015]).

Engineered circuits of MINPA (Figure 1B) and its sub-networks (Figure 1—figure supplement 1A) are designed to use two hybrid promoters, *Para/lac* and *Plux/tet,* which are characterized experimentally to show small leakage and high nonlinearity (Figure 1D–E and Figure 1—figure supplement 1B–D). For MINPA topology, hybrid promoter *Para/lac* drives *AraC* and *TetR* expression, representing the node X in Figure 1A, whereas *Plux/tet* controls *LuxR* and *LacI* transcription, representing the node Y. AraC and LuxR activate *Para/lac* and *Plux/tet* in the presence of Arabinose and AHL (3oxo-C6-HSL) respectively, forming positive autoregulations. IPTG inhibits the repressive effect of LacI on TetR expression, while aTc counteracts TetR repression on LacI. Hence, the two nodes form the topology presented in the conceptual design shown in Figure 1A. Green fluorescent protein (GFP) and mCherry serve as the corresponding readouts of *Plux/tet* and *Para/lac* activities in living cells (Figure 1B).

Topologies of MINPA and all its subnetworks can be divided into four layers, from one- to four-dimensional networks based on the number of regulatory edges (Figure 1C and Figure 1—figure supplement 1E–F) and further categorized into nine groups based on the configurations of activation and inhibition. By computationally searching a large parameter range for each of the nontrivial networks (Faucon et al., 2014), we found that networks with two auto-activations, including A^{2}, RA^{2}, R^{2}A^{2}, have high probability of tristability or quadrastability (Figure 1—figure supplement 1G), defined as having three or four stable steady sates (SSS) under a common induction condition. However, MINPA has broader parameter distributions than the other two (Figure 1—figure supplement 1H–I), which suggests it is more resistant to parameter change and thus likely to achieve multistability in experimental settings.

In order to experimentally evaluate dynamic properties of these networks, we constructed nine circuits including tunable positive feedbacks (T6 and T9), mutual inhibition (T5), dual-positive feedbacks (T10), and their combinations (T7, T11, T13, T14 and T15). One-dimensional networks (T1, T4, T2 and T8) and trivial two-dimensional networks (T3 and T12) are excluded for their low multistability probability. All motifs were constructed using the same set of components (Figure 1).

Probing a circuit’s multistability typically requires thorough hysteresis experiments covering wide ranges of doses for all inducers (Acar et al., 2005; Angeli et al., 2004; Gardner et al., 2000), which becomes infeasible for nine complex networks with four inducers. To improve the efficiency of probing multistability and tunability, we designed a ‘sequential induction’ method to accelerate exploration of unknown high dimensional bifurcation spaces (see Appendix text for details), instead of conventional ‘back and forth’ hysteresis on one parameter dimension. The main concept relies on the fact that multistable gene networks could exhibit discontinuous jump from one state to another in response to changing parameter (inducer) combinations. Taking the classic ‘toggle switch’ as an example (Gardner et al., 2000), the circuit can be tuned by two external inducers and its two-parameter bifurcation diagram has a stretched *S* shape (Figure 2A). Initialized at an arbitrary state A, the cells could reach State C in the bistable region directly when induced with both inducers simultaneously. If the cells are first induced by Inducer I to go to state B, they will also reach State C after Inducer II is added. However, if the same dose of Inducer II is applied first, cells will cross the bifurcation plane to state D on the low-Response surface and then reach state E with addition of Inducer I (Figure 2A). State C and E are two different steady states with the same induction dosages, illustrating hysteresis and verifying multistability.

To test our theoretical analysis, a synthetic toggle switch circuit was constructed (Figure 2—figure supplement 1A). Following experimental design principles (see Appendix text for details), we designed a protocol to show the sequential induction effects. We first employed IPTG to induce the circuit for 5 hr, and then aTc was added. Time course results showed that cells stayed at low-GFP state till 24 hr (Figure 2—figure supplement 1B). However, cells induced with aTc first, and then IPTG mainly stayed at high-GFP state, another stable steady state under this condition. Simultaneous aTc and IPTG induction produced similar cell distributions. These results show that sequential induction can be used as a strategy to quickly explore a multistable potential landscape for complex non-equilibrium systems.

Without knowing the exact bifurcation range beforehand, such ordered sequential inductions could help quickly explore the irregular bifurcation space to reveal multistability for systems with complicated bifurcations, which is typically caused by interfering parameters. Similar sequential induction techniques have been shown to enable access of otherwise hard-to-reach cell death states in breast cancer cells (Lee et al., 2012). This strategy has also been widely employed in directed differentiation of stem cells to specific lineages (Paşca et al., 2015; Pagliuca et al., 2014; Kroon et al., 2008) and reprogramming somatic cells to induced pluripotent stem cells (Liu et al., 2013). Although specific inducer concentrations are required to observe the effects of this strategy in synthetic circuits, sequential induction with pre-selected inducer combinations can help perform a coarse-grained exploration from different directions in the parameter space. Furthermore, stochastic gene expression of the circuits also contributes to cellular population distribution thus leads to pronounced sequential induction effects, given experimentally feasible amount of time, when the system is entering its multistable region from different directions. Therefore, distinct final states, or even different population distributions, under sequential induction strongly suggests the existence of nonlinear dynamics, including multistability (see Appendix text for details).

Using the sequential induction approach, we tested the nine circuits using flow cytometry. Cells were first induced by inducer I, inducer II was then added into the media for another 24 hr. Depending on the network configuration, four different dual-inducer combinations were used. For example, Arabinose and IPTG were applied sequentially and simultaneously to T9, T13, T11 and T15, respectively (Figure 2B). It can be seen only T15 exhibits significant expression difference between three induction patterns, while the others show little change (Figure 2B and Figure 2—figure supplement 2A). It should be noted that T15 also exhibits tri-modality of fluorescence expression, suggesting multistability given the presence of gene expression noise, which is partially consistent with our computational predictions. Similarly, AHL and aTc were applied to T6, T7, T14, and T15, respectively (Figure 2C and Figure 2—figure supplement 2B). Results show that only T15 exhibits significant fluorescence pattern change with different inductions, whereas T6 and T7 exhibit minor uniform shifts of expression. T14, although exhibiting bimodality, only shows a ratio change of two populations between three inductions and no sign of bifurcation. Sequential induction by Arabinose and AHL combinations has little effect on T10, T14 and T11, but T15 displays three notable populations for AHL-then-Arabinose induction (Figure 2D and Figure 2—figure supplement 2C). IPTG and aTc were also tested on T5, T7, T13 and T15, but no notable dynamics were observed (Figure 2—figure supplement 2D and Figure 2—figure supplement 3). Taken together, T15, the full MINPA topology, shows the most variety and complexity in population heterogeneity under sequential inductions, suggesting this circuit has the highest potential to generate complex multistability within our induction range and hence enable us to approach the Waddington landscape.

Next, operating principles and full tunability of T15 (MINPA) were further examined by using four inducers (Arabinose, AHL, aTc, and IPTG) to fine tune the strength of regulations and perturb the system (Figure 3A). Uninduced cells showed low GFP and low mCherry expression (low-low state, LL). In the presence of AHL and aTc, high GFP and low mCherry (GFP state) is observed; low GFP and high mCherry (mCherry state) emerged with induction of Arabinose; and high GFP and high mCherry (high-high state, HH) was achieved when induced with Arabinose and AHL. These results verify that our engineered MINPA circuit is functioning as designed and fully controllable with four distinct states reachable through appropriate inductions, respectively.

To help design experiments to further investigate the circuit’s quadrastability, a detailed mathematical model was developed to describe the system (see Appendix for details). Using parameters derived from hybrid promoter testing experiments, bifurcation analysis was carried out to systematically quantify MINPA’s dynamic behavior (Figure 3B, Figure 3—figure supplement 1 and Figure 3—figure supplement 2A–H). Figure 3B is the three-dimensional bifurcation diagram, where levels of GFP and mCherry represent the states of node X and Y, and ‘*AR/AL*’ is a lumped parameter composed of a fixed ratio of the concentrations of Arabinose and AHL. Overall, it can be seen that the system, initialized without induction, is predicted to be quadrastable (shown as four colored spheres, representing LL (grey), GFP (green), mCherry (rose), and HH (golden) state, respectively) but with the low-low state to have dominant attractiveness (shown as the big gray sphere) when *AR/AL* is low (C1). However, when *AR/AL* level is within an intermediate range, relative stabilities between different states become comparable. When *AR/AL* level increased from C1 to C2, the circuit’s quadrastability becomes well pronounced, illustrated as four similar-sized colored spheres on the same gray plane, which represents the low-low, GFP, mCherry, and high-high state, respectively (Figure 3—figure supplement 1). As *AR/AL* continues to increase from C2 to C3, while the other three SSS remain stable, the stability of the GFP branch disappears. Further increase of *AR/AL* results in only one stable state-the high-high state, shown as the orange sphere with biggest size.

To establish MINPA’s quadrastability and tristability as predicted, hysteresis, a hallmark of multistability (Acar et al., 2005; Wu et al., 2014, 2013), of the network was tested. Initialized at the low-low state, cells were induced by increasing doses of *AR/AL* corresponding to C1 to C4 and measured by flow cytometry (Figure 3C and Figure 3—figure supplement 2I). As predicted, C1_{LL} (cells with initial Low-Low state grown at C1 condition) experiment demonstrates uniform low-low fluorescence profile, due to the low-low state’s dominant attractiveness, and C4_{LL} shows a uniform high-high profile. Interestingly, C3_{LL} indeed illustrates tri-modality, which is the result of predicted tristability. C2_{LL} experiment, on the other hand, exhibits enough heterogeneity to signal high-high, low-low, and mCherry state, but does not illustrate significant trace of GFP state. Given that GFP state is achieved through combinational induction of AHL and aTc (Figure 3A), we hypothesize that the GFP state here is not easily accessible with AHL induction only. Next, cells initialized at high-high states were collected and diluted into fresh media with the same concentrations of *AR/AL* (Figure 3D and Figure 3—figure supplement 2J). As predicted, these cells keep high-high expression profile even with inductions as low as C1, another demonstration that the system is already multistable at C1. Taken together, the two sets of experiments demonstrated clear hysteresis and verified the existence of three of the four predicted SSS.

To further investigate what determines the accessibility of certain SSS in this quadrastable system and how cells navigate this attractor landscape, we take into account gene expression stochasticity (Wu et al., 2013) to sketch out MINPA’s quasi-potential attractor landscape (Figure 4A and Appendix), which is calculated as the negative logarithmic function of stationary distribution density in the phase space of GFP and mCherry. Using the weighted ensemble random walk algorithm (Appendix), the stationary density distribution can be efficiently calculated from the initial uniform distribution. It can be seen that when there is no inducer, MINPA is already quadrastable with four local minima, which is consistent with bifurcation analysis for C1 condition. Furthermore, the much stronger stability of the low-low state (deepest well, Top landscape) and high state-transition barrier explain homogeneous low-low population (C1 experiment in Figure 3C) when cells were initialized with no inductions.

Since Arabinose and AHL combination is not sufficient to enable the cells to reach all four SSS, we chose to add aTc to the mix to further facilitate cell transitions among these four SSS. Using our expanded model, we simulated simultaneous and sequential inductions and computed corresponding quasi-potential landscape (Figure 4A), showing cells harboring the same MINPA network exhibiting distinct landscapes under different inductions. AHL and aTc promote a more stable GFP state (Left center), while Arabinose induction modulates the landscape to be biased toward mCherry state (Right center). When the three inducers were applied simultaneously, the landscape changes and the four states show comparable stabilities (Bottom), suggesting a higher possibility of quadramodal cell population experimentally. Experimental validation is shown as flow cytometry measurements of cells treated with Arabinose, AHL, and aTc simultaneously for 24 hr (Figure 4B, and Figure 4—figure supplement 1A). Such a hybrid induction greatly facilitates the cells’ transition from low-low state to the other three states so that a quadramodal distribution emerges. Single-cell time lapse microscopy results also showed that the initial low-low state cells could differentiate into GFP, mCherry and high-high state cells (Figure 4—figure supplement 1B–D and Appendix 1—Video 1). This also finally verifies predicted quadrastability of MINPA.

There are two other strategies to reach this condition: sequential inductions with AHL-and-aTc and then Arabinose (Figure 4A, Left route) or Arabinose and then AHL-and-aTc (Right route). Even though the initial and final landscapes are the same, the dynamics for each route are quite different, which could lead to distinct outcomes. By comparing state barrier heights (Figure 4A), we hypothesize that cells walking through the left route would start transitioning from low-low state to GFP state upon induction of AHL and aTc. Following Arabinose induction would then make the mCherry state accessible. So some cells with GFP state would transition to high-high state while some low-low state cells transition to mCherry state, resulting in cells in all four states. Experimental testing indeed shows four stable populations (Figure 4C). At 6.5 hrs of AHL and aTc induction, about 12% cells were moving to GFP state while the rest of them still stay ‘undecided’ at low-low state (Figure 4—figure supplement 1A and E). This is consistent with the simulated landscape as these two states are more stable and accessible to each other (Figure 4A, Left). Arabinose induction promoted some cells to transition into mCherry state while some cells continued moving into GFP state, of which some further transitioned to high-high state.

Interestingly, the right route is predicted to generate different results. When first induced with Arabinose, the mCherry valley is so deep that it would be difficult for cells to jump out to high-high state, and low-low state cells are also hardly transit to GFP state due to its low attractiveness, and thus most cells would stay at mCherry and low-low state even with AHL and aTc inductions (Figure 4A, Right). Experimental testing of the right route indeed only produces two populations with low-low and mCherry state (Figure 4D). With 5 hrs of Arabinose induction, most cells still stay at low-low state because of slow transition to the mCherry state (Figure 4—figure supplement 1A), but 84.6% cells transitioned to mCherry state with 15.3% cells at low-low state at 9.5 hr (Figure 4—figure supplement 1A). This is consistent with our model predictions. The high barrier between the mCherry state and high-high state blocks the transition from mCherry state to high-high state, while the low attractiveness and relatively high barrier of the GFP state also decreases the probability of cells transitioning from low-low to GFP state. Hence, when AHL and aTc are applied, cells are predominantly in the mCherry state with a small portion in low-low state with low probability of transitioning out, resulting in a bimodal distribution.

Multistability and the resulting landscape has long been proposed as an underlying mechanism that cells use to maintain pluripotency and guide differentiation (Kauffman, 1993; Laurent and Kellershohn, 1999; Huang et al., 2007; Guantes and Poyatos, 2008; Palani and Sarkar, 2009; Narula et al., 2010; Faucon et al., 2014). Theoretical frameworks have also been established to quantify the Waddington landscape and biological paths for cell development (Li and Wang, 2013a, 2013b; Wang et al., 2011). Experimental validation of this hypothesis and a full understanding of this mechanism will help reveal differentiation dynamics and routes for all cell types, which remains an outstanding problem in biology.

In this study, we engineered the quadrastable MINPA circuit and show that it can guide cell fate choices, represented by fluorescence expression, through shaping the potential landscape. MINPA represents one of the most complicated two-node network topologies and includes four genes to implement a web of regulations. Biological complexity correlates with the number of regulatory connections (Szathmáry et al., 2001), not the number of genes. Hence, dense connectivity and complex dynamics of MINPA may provide a framework to understand similarly densely connected gene regulatory networks.

Combining mathematical modeling and experimental investigation, this study serves as a proof-of-principle demonstration of the Waddington landscape. Furthermore, we used this circuit to demonstrate how different sequential inductions can change the landscape in a specific order and navigate cells to different final states. Such illustrations suggest mechanistic explanations of the need for fixed induction sequences for targeted differentiation to desired cell lineage. Overall, this study helps reveal fundamental mechanisms of cell-fate determination and provide a theoretical foundation for systematic understanding of the cell differentiation process, which will lead to development of new strategies to program cell fate.

All the molecular cloning experiments were performed in *E.coli* DH10B (Invitrogen, USA), and measurements of MINPA and sub-networks were conducted in *E.coli* K-12 MG1655Δ*lacI*Δ*araCBAD* strain as previously described (from Dr. Collins Lab [Litcofsky et al., 2012]). The sequential induction for the toggle circuit was conducted in *E.coli* MG1655Δ*lacI* strain as previously described (Litcofsky et al., 2012). Cells were grown at 37°C in liquid and/or solid Luria-Bertani broth medium with 100 µg/mL ampicillin or kanamycin. Chemicals AHL (3oxo-C6-HSL, Sigma-Aldrich), Arabinose (Sigma-Aldrich, USA), isopropyl β-D-1-thiogalactopyranoside (IPTG, Sigma-Aldrich), and anhydrotetracycline (aTc, Sigma-Aldrich) were dissolved in ddH2O and diluted into indicated working concentrations. Chemical aTc solution was stocked in brown vials, and experiments involving aTc were performed in cabinet without light, and cell cultures were grown in darken incubator at 37°C. Cultures were shaken in 5 mL and/or 15 mL tubes at 220 rotations per minute (r.p.m).

All the plasmids (MINPA and its nine sub-networks) in this study were constructed using standard molecular cloning protocols and assembled by standardized BioBricks methods based on primary modules (Table 1) from the iGEM Registry (www.parts.igem.org). Hybrid promoter *Para/lac* was from Dr. Collins lab and amplified using forward primer: *CGGAATTCGCTTCTAGAGAATTGTGAGCGGATAAC*; and reverse primer: *CGCTGCAGGCACTAGTTTGTGTGAAATTGTTATCCG*. PCR product was purified using GenElute PCR Clean-Up Kit (Sigma-Aldrich), and then cut by restriction enzymes *EcoRI* and *PstI*. The purified product was inserted into pSB1K3 backbone, and finally verified by DNA sequencing. The MINPA circuit was constructed from promoter *Para/lac* and nine other Biobrick standard biological parts: BBa_B0034 (ribosome binding site, RBS), BBa_C0080 (*AraC* gene), BBa_C0040 (*tetR* gene), BBa_K176000 (*Plux/tet* hybrid promoter), BBa_C0062 (*luxR* gene), BBa_C0012 (*lacI* gene), BBa_B0015 (transcriptional terminator), BBa_E0240 (GFP generator), and BBa_J06702 (mCherry generator). The fragment and vector were separated by gel electrophoresis (1% TAE agarose) and purified using GenElute Gel Extraction Kit (Sigma-Aldrich). Then, fragment and vector were ligated together using T4 DNA ligase, and the ligation products were transformed into *E. coli* DH10B and clones were screened by plating on 100 μg/ml ampicillin LB agar plates. Finally, their plasmids were extracted and verified by double digestion (*EcoRI* and *PstI*). The detailed procedures of assembling DNA constructs were described in our previous study (Wu et al., 2014). Restriction enzymes (*EcoRI*, *XbaI*, *SpeI*, and *PstI*) and T4 DNA ligase were purchased from New England Biolabs. All the constructs were inserted into high copy number plasmid pSB1A3 and pSB1K3. All the constructs were verified by DNA sequencing (Biodesign sequencing lab in ASU) step by step.

All the samples were analyzed at the indicated time points on an Accuri C6 flow cytometer (Becton Dickinson, USA) with excitation/emission filters (488/530 nm for GFP, and 610 LP for mCherry). The data were collected in a linear scale and non-cellular low-scatter noise was removed by thresholding. All measurements of gene expression were obtained from at least three independent experiments. For each culture, 100,000 events were collected at a slow flow rate. Data files were analyzed using MATLAB (MathWorks).

For sequential induction, initially uninduced overnight cell culture was diluted into fresh media without or with inducer I, grown at 37°C and 220 r.p.m till OD_{600} is 0.15 ~ 0.25 (the time usually takes 5 ~ 6.5 hr, depends on the inducers and concentrations). For samples induced individually by Ara, or AHL, or IPTG, it is ~5 hr; for samples induced with aTc, it takes ~6.5 hr. According to our experience, gene (GFP) is starting to be partially expressed while steady states are not yet stable. Then inducer II was added into the culture, and grown for another 24 hr. Flow cytometry was performed at 0 hr, 12 hr, and 24 hr after the second inducer was added into the culture. For each set of sequential induction, the first scenario: add inducer I first, then add inducer II; the second scenario: add inducer II first, then add inducer I; the third scenario: add inducers I and II at the same time. As a control, cells without any inducer were also prepared and measured. Inducer I and II were the two of four commercial chemicals: AHL, Arabinose, IPTG, and aTc. All the experiments were repeated for at least three times and only representative results were showed.

For hysteresis experiments, initially uninduced cells were diluted into fresh media and distributed into new 5 ml tubes. Various amounts of Arabinose and AHL (3oxo-C6-HSL) were added into the media, and cells were then grown at 37°C shaker. The initially high-high state cells induced with 2.5 *10^{−3} m/v Arabinose and 1*10^{−4} M AHL were collected with low-speed centrifugation, washed twice, resuspended with fresh medium, and at last inoculated into fresh medium at a 1:100 ratio with the same series of inducer (Arabinose and AHL) concentrations. C1, C2, C3, and C4 (Figure 3B–D) are four increasing concentrations of Arabinose and AHL used for experimental probing, but the ratio of Arabinose and AHL is fixed. Specifically, cells were induced with the Arabinose and AHL at the same time (the third scenario), at concentrations from C1 to C4. C1: no inducers; C2: 2.5*10^{−6}m/v Arabinose and 1*10^{−7} M AHL; C3: 2.5*10^{−5}m/v Arabinose and 1*10^{−6} M AHL; C4: 2.5*10^{−3}m/v Arabinose and 1*10^{−4} M AHL. Flow cytometry analyses were performed at 12 hr and 24 hr to monitor the fluorescence levels. Experiments were repeated two times with three replicates.

Cells with MINPA circuit were grown overnight, which was then re-diluted into 5 mL fresh LB medium with Kanamycin the next day. When OD_{600} of the cells reached about 0.2, cells were spun down with low speed and resuspended in 5 ml of fresh medium and loaded into the device. Detailed description of chip design and device setup could be found from Hasty Lab (Ferry et al., 2011). Two media were prepared: one with inducers and the other without. Cells in the trap were first supplied by the medium without inducer for 6 hr, and then switched to medium with inducers for anther 18 hr, which was controlled by adjusting the heights of the medium syringes relative to one another. Images were taken by using Nikon Eclipse Ti inverted microscope (Nikon, Japan) equipped with an LED-based Lumencor SOLA SE. Phase and fluorescence images were taken every 5 min for 24 hr in total under the magnification 40x. Perfect focus was maintained automatically using Nikon Elements software. Experimental detail can also be found in Appendix.

Ordinary differential equation models were developed to describe and analyze the MINPA system. Details are provided in the Appendix.

We would like to thank Dr. James J Collins for the plasmid Para/lac and the *E.coli* K-12 MG1655Δ*lacI* and MG1655Δ*lacI*Δ*araCBAD* strain. We also thank Dr. Jeff Hasty for the microfluidic chip and setup. We thank Dr. Alexander Green for critical reading of the manuscript and great comments. FW was supported by American Heart Association Predoctoral Fellowship 15PRE25710303. YCL was supported by ARO under Grant No.W911NF-14-1-0504. This study was financially supported by National Science Foundation Grant DMS-1100309, American Heart Association grant 11BGIA7440101, and National Institutes of Health Grant GM106081 (to XW).

Much of the gene regulation in our circuit involves multimerization of protein. These reactions can be described as:

$$\begin{array}{rl}2\left[luxR\right]\phantom{\rule{thinmathspace}{0ex}}& \underset{{k}_{-}^{u}}{\overset{\phantom{\rule{thinmathspace}{0ex}}{k}_{+}^{u}\phantom{\rule{thinmathspace}{0ex}}}{\rightleftharpoons}}\phantom{\rule{thinmathspace}{0ex}}\left[lux{R}_{2}\right]\\ 2\left[tetR\right]\phantom{\rule{thinmathspace}{0ex}}& \underset{{k}_{-}^{t}}{\overset{\phantom{\rule{thinmathspace}{0ex}}{k}_{+}^{t}\phantom{\rule{thinmathspace}{0ex}}}{\rightleftharpoons}}\phantom{\rule{thinmathspace}{0ex}}\left[tet{R}_{2}\right]\\ 2\left[LacI\right]\phantom{\rule{thinmathspace}{0ex}}& \underset{{k}_{-}^{l}}{\overset{\phantom{\rule{thinmathspace}{0ex}}{k}_{+}^{l}\phantom{\rule{thinmathspace}{0ex}}}{\rightleftharpoons}}\phantom{\rule{thinmathspace}{0ex}}\left[Lac{I}_{2}\right]\\ 2\left[Lac{I}_{2}\right]\phantom{\rule{thinmathspace}{0ex}}& \underset{{k}_{-}^{l,d}}{\overset{\phantom{\rule{thinmathspace}{0ex}}{k}_{+}^{l,d}\phantom{\rule{thinmathspace}{0ex}}}{\rightleftharpoons}}\phantom{\rule{thinmathspace}{0ex}}\left[Lac{I}_{4}\right]\\ 2\left[AraC\right]\phantom{\rule{thinmathspace}{0ex}}& \underset{{k}_{-}^{a}}{\overset{\phantom{\rule{thinmathspace}{0ex}}{k}_{+}^{a}\phantom{\rule{thinmathspace}{0ex}}}{\rightleftharpoons}}\phantom{\rule{thinmathspace}{0ex}}\left[Ara{C}_{2}\right].\end{array}$$

From these basic processes, we can calculate the concentrations for all the dimers and *LacI* tetrameter (denoted using subscripts):

$$\begin{array}{rlr}\left[lux{R}_{2}\right]\phantom{\rule{thinmathspace}{0ex}}& =\phantom{\rule{thinmathspace}{0ex}}\frac{{k}_{+}^{u}}{{k}_{-}^{u}}[luxR{]}^{2}& \\ \left[Ara{C}_{2}\right]\phantom{\rule{thinmathspace}{0ex}}& =\phantom{\rule{thinmathspace}{0ex}}\frac{{k}_{+}^{a}}{{k}_{-}^{a}}[AraC{]}^{2}& \\ \left[tet{R}_{2}\right]\phantom{\rule{thinmathspace}{0ex}}& =\phantom{\rule{thinmathspace}{0ex}}\frac{{k}_{+}^{t}}{{k}_{-}^{t}}[tetR{]}^{2}& \\ \left[Lac{I}_{4}\right]\phantom{\rule{thinmathspace}{0ex}}& =\phantom{\rule{thinmathspace}{0ex}}\frac{{k}_{2+}^{I}}{{k}_{2-}^{I}}(\frac{{k}_{+}^{I}}{{k}_{-}^{I}}{)}^{2}[LacI{]}^{4}.& \end{array}$$

We use ${P}_{i,j}^{A}$ to denote promoter *Para/lac* with $i[Ara{C}_{2}]$ and $j[Lac{I}_{4}]$ binding, where $i=0,1$ and $j=0,1,2$ since *Para/lac* has two binding sites for $[Lac{I}_{4}]$. Similarly, we use ${P}_{i,j}^{B}$ to denote promoter *Plux/tet* with $i[lux{R}_{2}]$ and $j[tet{R}_{2}]$ binding, where both $i$ and $j$ can be 0 or 1. The association dynamics for these two hybrid promoters can be described as:

$$\begin{array}{rl}\left[Ara{C}_{2}\right]+{P}_{0,j}^{A}\phantom{\rule{thinmathspace}{0ex}}& \underset{{S}_{-}^{a}}{\overset{\phantom{\rule{thinmathspace}{0ex}}{S}_{+}^{a}\phantom{\rule{thinmathspace}{0ex}}}{\rightleftharpoons}}\phantom{\rule{thinmathspace}{0ex}}{P}_{1,j}^{A},\phantom{\rule{2em}{0ex}}j\in \{0,1,2\}\phantom{\rule{thinmathspace}{0ex}}\\ \left[Lac{I}_{4}\right]+{P}_{i,0}^{A}\phantom{\rule{thinmathspace}{0ex}}& \underset{{S}_{-}^{l}}{\overset{\phantom{\rule{thinmathspace}{0ex}}{S}_{+}^{l}\phantom{\rule{thinmathspace}{0ex}}}{\rightleftharpoons}}\phantom{\rule{thinmathspace}{0ex}}{P}_{i,1}^{A},\phantom{\rule{2em}{0ex}}i\in \{0,1\}\phantom{\rule{thinmathspace}{0ex}}\\ \left[Lac{I}_{4}\right]+{P}_{i,1}^{A}\phantom{\rule{thinmathspace}{0ex}}& \underset{{S}_{-}^{l}}{\overset{\phantom{\rule{thinmathspace}{0ex}}{S}_{+}^{l}\phantom{\rule{thinmathspace}{0ex}}}{\rightleftharpoons}}\phantom{\rule{thinmathspace}{0ex}}{P}_{i,2}^{A},\phantom{\rule{2em}{0ex}}i\in \{0,1\}\phantom{\rule{thinmathspace}{0ex}}\\ \left[lux{R}_{2}\right]+{P}_{0,j}^{B}\phantom{\rule{thinmathspace}{0ex}}& \underset{{S}_{-}^{u}}{\overset{\phantom{\rule{thinmathspace}{0ex}}{S}_{+}^{u}\phantom{\rule{thinmathspace}{0ex}}}{\rightleftharpoons}}\phantom{\rule{thinmathspace}{0ex}}{P}_{1,j}^{B},\phantom{\rule{2em}{0ex}}j\in \{0,1\}\phantom{\rule{thinmathspace}{0ex}}\\ \left[tet{R}_{2}\right]+{P}_{i,0}^{B}\phantom{\rule{thinmathspace}{0ex}}& \underset{{S}_{-}^{t}}{\overset{\phantom{\rule{thinmathspace}{0ex}}{S}_{+}^{t}\phantom{\rule{thinmathspace}{0ex}}}{\rightleftharpoons}}\phantom{\rule{thinmathspace}{0ex}}{P}_{i,1}^{B},\phantom{\rule{2em}{0ex}}i\in \{0,1\},\phantom{\rule{thinmathspace}{0ex}}\end{array}$$

where ${S}_{+}$ and ${S}_{-}$ are protein-DNA association and disassociation rates, respectively, and the superscript represents corresponding repressor or activator, where $a$ represent rates for $\left[Ara{C}_{2}\right]$, $l$ for $\left[Lac{I}_{4}\right]$, $t$ for $\left[tet{R}_{2}\right]$ and $u$ for $\left[lux{R}_{2}\right]$. For simplicity, we have omitted the looping process of plasmid.

Using these binding/unbinding relationships, we can calculate the portion of ${P}_{1,0}^{A}$ and ${P}_{1,0}^{B}$ promoters. Both of them have maximal production rates compared to other promoter binding states, thus they are the dominant states among all possible binding states in determining circuit activities. From the first and second equations we can write down the ratio ${R}_{1,j}^{A}$ for promoter binded by $[Ara{C}_{2}]$ (regardless of the binding of protein $[Lac{I}_{4}]$) is:

$$\begin{array}{c}{R}_{1,j}^{A}=\frac{{P}_{1,j}^{A}}{{P}_{1,j}^{A}+{P}_{0,j}^{A}}=\frac{\frac{{S}_{+}^{a}}{{S}_{-}^{a}}[Ara{C}_{2}]}{1+\frac{{S}_{+}^{a}}{{S}_{-}^{a}}[Ara{C}_{2}]}.\end{array}$$

(1)

We further analyze the ratios of promoters binded by $[Lac{I}_{4}]$. Since there are two binding sites for $[Lac{I}_{4}]$, their dynamics can be expressed as:

$${P}_{i,0}^{A}\underset{{S}_{-}^{l}}{\overset{2{S}_{+}^{l}\cdot [Lac{I}_{4}]}{\rightleftharpoons}}{P}_{i,1}^{A}\underset{2\cdot {S}_{-}^{l}}{\overset{\phantom{\rule{thinmathspace}{0ex}}{S}_{+}^{l}\cdot [Lac{I}_{4}]\phantom{\rule{thinmathspace}{0ex}}}{\rightleftharpoons}}{P}_{i,2}^{A}.$$

(2)

We can write down the dynamical functions for promoter ${P}_{i,0}^{A}$ and ${P}_{i,2}^{A}$ as:

$$\frac{d{P}_{i,0}^{A}}{dt}\phantom{\rule{thinmathspace}{0ex}}=\phantom{\rule{thinmathspace}{0ex}}{S}_{-}^{l}\cdot {P}_{i,1}^{A}-2\cdot {S}_{+}^{l}[Lac{I}_{4}]\cdot {P}_{i,0}^{A},$$

(3)

and

$$\frac{d{P}_{i,2}^{A}}{dt}\phantom{\rule{thinmathspace}{0ex}}=\phantom{\rule{thinmathspace}{0ex}}{S}_{+}^{l}[Lac{I}_{4}]\cdot {P}_{i,1}^{A}-2\cdot {S}_{-}^{l}\cdot {P}_{i,2}^{A}.$$

(4)

So the stable ratio for ${P}_{i,0}^{A}$, or promoter without $[Lac{I}_{4}]$ binding, is:

$${R}_{i,0}^{A}=\frac{{P}_{i,0}^{A}}{{P}_{i,0}^{A}+{P}_{i,1}^{A}+{P}_{i,2}^{A}}=\frac{1}{{\left(1+\frac{{S}_{+}^{l}}{{S}_{-}^{l}}[Lac{I}_{4}]\right)}^{2}}.$$

(5)

Similarly, we can write down equations to describe $[lux{R}_{2}]$ and $[tet{R}_{2}]$ binding activities. Their dynamics can be expressed as:

$$\begin{array}{rl}\frac{d{P}_{1,j}^{B}}{dt}\phantom{\rule{thinmathspace}{0ex}}& =\phantom{\rule{thinmathspace}{0ex}}{S}_{+}^{u}[lux{R}_{2}]\cdot {P}_{1,j}^{B}-{S}_{-}^{u}\cdot {P}_{0,j}^{B},\phantom{\rule{thinmathspace}{0ex}}\\ \frac{d{P}_{i,0}^{B}}{dt}\phantom{\rule{thinmathspace}{0ex}}& =\phantom{\rule{thinmathspace}{0ex}}{S}_{-}^{t}\cdot {P}_{i,1}^{B}-{S}_{+}^{t}[tet{R}_{2}]\cdot {P}_{i,0}^{B},\phantom{\rule{thinmathspace}{0ex}}\end{array}$$

and then solve the stable ratios ${R}_{1,j}^{B}$ and ${R}_{i,0}^{B}$ as:

$$\begin{array}{r}\text{(6)}& {R}_{1,j}^{B}\phantom{\rule{thinmathspace}{0ex}}& =\phantom{\rule{thinmathspace}{0ex}}\frac{\frac{{S}_{+}^{u}}{{S}_{-}^{u}}[lux{R}_{2}]}{1+\frac{{S}_{+}^{u}}{{S}_{-}^{u}}[lux{R}_{2}]},& \text{(7)}& {R}_{i,0}^{B}\phantom{\rule{thinmathspace}{0ex}}& =\phantom{\rule{thinmathspace}{0ex}}\frac{1}{1+\frac{{S}_{+}^{t}}{{S}_{-}^{t}}[tet{R}_{2}]}.& \end{array}$$

So at steady state, the ratios are controlled by association rates and the concentration of activator and inhibitors, as followed:

$$\begin{array}{r}\text{(8)}& {R}_{1,0}^{A}\phantom{\rule{thinmathspace}{0ex}}& =\phantom{\rule{thinmathspace}{0ex}}{R}_{1,j}^{A}\times {R}_{i,0}^{A}=\frac{\frac{{S}_{+}^{a}}{{S}_{-}^{a}}[Ara{C}_{2}]}{1+\frac{{S}_{+}^{a}}{{S}_{-}^{a}}[Ara{C}_{2}]}\times \frac{1}{{\left(1+\frac{{S}_{+}^{l}}{{S}_{-}^{l}}[Lac{I}_{4}]\right)}^{2}},& \text{(9)}& {R}_{1,0}^{B}\phantom{\rule{thinmathspace}{0ex}}& =\phantom{\rule{thinmathspace}{0ex}}{R}_{1,j}^{B}\times {R}_{i,0}^{B}=\frac{\frac{{S}_{+}^{u}}{{S}_{-}^{u}}[lux{R}_{2}]}{1+\frac{{S}_{+}^{u}}{{S}_{-}^{u}}[lux{R}_{2}]}\times \frac{1}{1+\frac{{S}_{+}^{t}}{{S}_{-}^{t}}[tet{R}_{2}]}.& \end{array}$$

Here, we model how the association rates and disassociation rates are affected by two inducers [Arabinose] and [AHL], as followed (Stricker et al., 2008):

$$\begin{array}{r}\text{(10)}& {\displaystyle {S}_{a}\phantom{\rule{thinmathspace}{0ex}}}& =\phantom{\rule{thinmathspace}{0ex}}\frac{{S}_{+}^{a}}{{S}_{-}^{a}}={C}_{min}^{a}+[{C}_{max}^{a}-{C}_{min}^{a}]\cdot \frac{[Arabinose{]}^{{n}_{a}}}{[Arabinose{]}^{{n}_{a}}+{K}_{a}^{{n}_{a}}},& \text{(11)}& {\displaystyle {S}_{u}\phantom{\rule{thinmathspace}{0ex}}}& =\phantom{\rule{thinmathspace}{0ex}}\frac{{S}_{+}^{u}}{{S}_{-}^{u}}={C}_{min}^{u}+[{C}_{max}^{u}-{C}_{min}^{u}]\cdot \frac{[AHL{]}^{{n}_{u}}}{[AHL{]}^{{n}_{u}}+{K}_{u}^{{n}_{u}}},& \end{array}$$

and for the inhibition by aTc and IPTG as:

$$\begin{array}{r}\text{(12)}& {S}_{t}\phantom{\rule{thinmathspace}{0ex}}& =\phantom{\rule{thinmathspace}{0ex}}\frac{{S}_{+}^{t}}{{S}_{-}^{t}}={C}_{min}^{t}+[{C}_{max}^{t}-{C}_{min}^{t}]\cdot \frac{{K}_{t}^{{n}_{t}}}{{K}_{t}^{{n}_{t}}+[aTc{]}^{{n}_{t}}},& \text{(13)}& {S}_{l}\phantom{\rule{thinmathspace}{0ex}}& =\phantom{\rule{thinmathspace}{0ex}}\frac{{S}_{+}^{l}}{{S}_{-}^{l}}={C}_{min}^{l}+[{C}_{max}^{l}-{C}_{min}^{l}]\cdot \frac{{K}_{l}^{{n}_{l}}}{{K}_{l}^{{n}_{l}}+[IPTG{]}^{{n}_{l}}}.& \end{array}$$

There are always leakage for promoters, so we further construct a leakage model to approximate these effects. $p$ is used to denote the basal production rate. Also, we use $a$ to denote the ratio of production rates between non-activator-binding promoters and activator-binding promoters, and $r$ to represent the ratio of production rates between inhibitor-binding promoters and non-inhibitor-binding promoters. Since the fluorescence strength can be directly measured and are proportional to gene expression levels, we use the concentration of mCherry protein, $[mC]$, to represent the activity of promoter *Para/lac* and the concentrations of *LacI* and *AraC*. Similarly, we use $[GFP]$ to represent the concentration of GFP protein and further the concentrations of *tetR* and *luxR*.

When leakage is introduced, we can write down the dynamical functions for MINPA as followed:

$$\begin{array}{r}\text{(14)}& \frac{d[mC]}{dt}\phantom{\rule{thinmathspace}{0ex}}& =\phantom{\rule{thinmathspace}{0ex}}{p}_{m}+{A}_{mC}\left[{a}_{m}(1-{R}_{1,j}^{A})+{R}_{1,j}^{A}\right]\times \left[{r}_{m}(1-{R}_{j,0}^{A})+{R}_{j,0}^{A}\right]-d\cdot [mC]& \end{array}$$

$$\begin{array}{c}\text{(15)}& {\displaystyle \begin{array}{rcc}d[GFP]dt\phantom{\rule{thinmathspace}{0ex}}& =\phantom{\rule{thinmathspace}{0ex}}{p}_{g}+{A}_{GFP}\left[{a}_{g}(1-{R}_{1,j}^{B})+{R}_{1,j}^{B}\right]\times \left[{r}_{g}(1-{R}_{j,0}^{B})+{R}_{j,0}^{B}\right]-d\cdot [GFP].\phantom{\rule{thinmathspace}{0ex}}& \end{array}}\end{array}$$

Here, ${p}_{m}$ and ${p}_{g}$ are basal leakage levels. ${a}_{m}$ and ${r}_{m}$ are reduced ratios of production rate for promoter *Para/lac*, while ${a}_{g}$ and ${r}_{g}$ are for promoter *Plux/tet*. ${A}_{mC}$ is the maximal production rates for *Para/lac* and ${A}_{GFP}$ is the maximal production rate for *Plux/tet*. We can further simplify this model to be:

$$\begin{array}{rlr}\text{(16)}& \frac{d[mC]}{dt}\phantom{\rule{thinmathspace}{0ex}}& =\phantom{\rule{thinmathspace}{0ex}}{p}_{m}+{A}_{mC}^{\mathrm{\prime}}\left[{a}_{m}^{\mathrm{\prime}}+{R}_{1,j}^{A}\right]\times \left[{r}_{m}^{\mathrm{\prime}}+{R}_{j,0}^{A}\right]-d\cdot [mC]& \frac{d[GFP]}{dt}\phantom{\rule{thinmathspace}{0ex}}& =\phantom{\rule{thinmathspace}{0ex}}{p}_{g}+{A}_{GFP}^{\mathrm{\prime}}\left[{a}_{g}^{\mathrm{\prime}}+{R}_{1,j}^{B}\right]\times \left[{r}_{g}^{\mathrm{\prime}}+{R}_{j,0}^{B}\right]-d\cdot [GFP],& \end{array}$$

where ${A}_{mC}^{\prime}={A}_{mC}\cdot (1-{a}_{m})\cdot (1-{r}_{m})$, ${A}_{GFP}^{\prime}={A}_{GFP}\cdot (1-{a}_{g})\cdot (1-{r}_{g})$, ${a}_{m}^{\prime}=\frac{{a}_{m}}{1-{a}_{m}}$, ${a}_{g}^{\prime}=\frac{{a}_{g}}{1-{a}_{g}}$, ${r}_{m}^{\prime}=\frac{{r}_{m}}{1-{r}_{m}}$ and ${r}_{g}^{\prime}=\frac{{r}_{g}}{1-{r}_{g}}$.

In order to fit the parameters for both promoter association and the dosage inductions, we construct and test the promoter expression rates in new constructed gene motifs as shown in Figure 1—figure supplement 1B–C.

As shown in Figure 1D–E, we perform dosage response experiments for *Para/lac* and *Plux/tet*. In each experiment, we grow the cells containing these gene circuits with different dosage combinations of inducer concentrations ([Arabinose] and [IPTG] for *Para/lac*, or [AHL] and [aTc] for *Plux/tet*) and then measure fluorescence strengths.

In both promoter setups, the concentrations of activators and inhibitors are controlled by constitutive promoters. According to the Equation (14) and Equation (15), the stable fluorescence strengths mCherry and GFP are determined by inducer concentrations. The stable concentrations for [mC] and [GFP] are:

$$\begin{array}{rlr}[mC{]}^{\ast}([Arabinose],[IPTG])\phantom{\rule{thinmathspace}{0ex}}& =\phantom{\rule{thinmathspace}{0ex}}\frac{1}{d}\left[{p}_{m}+{A}_{mC}^{\mathrm{\prime}}{H}_{a}([Arabinose])\times {H}_{l}([IPTG])\right]& \\ \text{(17)}& [GFP{]}^{\ast}([AHL],[aTc])\phantom{\rule{thinmathspace}{0ex}}& =\phantom{\rule{thinmathspace}{0ex}}\frac{1}{d}\left[{p}_{g}+{A}_{GFP}^{\mathrm{\prime}}{H}_{u}([AHL])\times {H}_{t}([aTc])\right],& \end{array}$$

where ${H}_{a}([Arabinose])={a}_{m}^{\prime}+{R}_{1,j}^{A}$, ${H}_{l}([IPTG])={r}_{m}^{\prime}+{R}_{i,0}^{A}$, ${H}_{u}([AHL])={a}_{g}^{\prime}+{R}_{1,j}^{B}$ and ${H}_{t}([aTc])={r}_{g}^{\prime}+{R}_{i,0}^{B}$ are dosage response functions that need to be fitted. To isolate response functions for different dosages, we make the assumption that inducers function independently, thus we can isolate and then fit the Hill functions from the dosage response experiment.

Here we use principle component analysis (PCA) to analyze the dosage response data. PCA is a matrix decomposition method which project a matrix $M$ into a set of values of linearly uncorrelated variables called principal components, as $M=U\mathrm{\Sigma}{V}^{*}$. Here $\mathrm{\Sigma}$ is a diagonal matrix with non-negative real number on the diagonal, which represent the variability of the corresponding vectors in $U$ and $V$. Since we assume that the dosage response data are product of two single variable function, we can use the vectors corresponding to the largest variability in $U$ and $V$ as the dosage response functions.

We take the response functions for promoter *Para/lac* as an example to explain our fitting procedure. First we subtract different leakage levels of ${p}_{m}$ from the original data, and perform PCA to obtain response functions for activator and repressor, respectively. The decomposition vectors corresponding to the largest factors are considered as the response functions for ${H}_{a}([Arabinose])$ and ${H}_{l}([IPTG])$, which are plotted as circles in Figure 1—figure supplement 1 D1 and D2 . We calculate the error between the experimental data and the model prediction from response functions, and then decide the optimal value for ${p}_{m}$. We can further determine the value for ${p}_{g}$ and the normalized ${H}_{u}([AHL])$ and ${H}_{t}([aTc])$ function in a similar approach (see circles in Figure 1—figure supplement 1 D3 and D4).

When the response functions are obtained, we can use the theoretical models for dosage response to fit the optimal leakage terms and the parameters in Hill functions (red solid lines in Figure 1—figure supplement 1D).

Using the models and these parameters, we can then determine the quantitative models for MINPA and the other circuits in Figure 1—figure supplement 1A.

We use the exhaustive parameter searching (Faucon et al., 2014; Ma et al., 2009) and simplified ODE models to quantitatively study the multistability of different circuits listed in Figure 1—figure supplement 1A. If a given circuit has larger possibility to be tuned by combination of inducers into multiple states, we consider such circuit has larger multistability likelihood.

First, we simplify the ODE model for circuit T15, the MINPA, and then try to derive simplified models for all the other circuits. The model for circuit T15 are:

$$\begin{array}{rlr}\frac{d[x]}{dt}\phantom{\rule{thinmathspace}{0ex}}& =\phantom{\rule{thinmathspace}{0ex}}{p}_{0}+{A}_{mC}^{\mathrm{\prime}}\left({a}_{0}^{\mathrm{\prime}}+\frac{{S}_{a}[x{]}^{4}}{{S}_{a}[x{]}^{4}+1}\right)\times \left({r}_{0}^{\mathrm{\prime}}+\frac{1}{({S}_{l}[y{]}^{4}+1{)}^{2}}\right)-d\cdot [x]& \\ \frac{d[y]}{dt}\phantom{\rule{thinmathspace}{0ex}}& =\phantom{\rule{thinmathspace}{0ex}}{p}_{0}+{A}_{GFP}^{\mathrm{\prime}}\left({a}_{0}^{\mathrm{\prime}}+\frac{{S}_{u}[y{]}^{4}}{{S}_{u}[y{]}^{4}+1}\right)\times \left({r}_{0}^{\mathrm{\prime}}+\frac{1}{{S}_{t}[x{]}^{4}+1}\right)-d\cdot [y],& \end{array}$$

where we use $[x]$ to represent the activity of promoter *Para/lac*, and $[y]$ for promoter *Plux/tet*, respectively. The leakage levels, which are ${p}_{0}$, ${a}_{0}$ and ${r}_{0}$, are used to denote the leakage in promoter, activation and inhibition of the promoters. We assume that the leakage in promoter are small and symmetrical for both hybrid promoters in all circuits. For simplicity, we use the inducer affinity terms, e.g., ${S}_{a}$, ${S}_{u}$, ${S}_{t}$ and ${S}_{l}$, to represent the actual dosage induction effects. Generally, a larger affinity term for ${S}_{a}$ and ${S}_{u}$ represent higher dosage levels and also higher affinity rates. The meaning of ${S}_{t}$ and ${S}_{l}$ are on the contrary, when larger values represent lower dosage concentration.

When the corresponding genes are removed from circuit T15, we change the value of affinity to address the topology changes. If the corresponding self-activation is removed, e.g., in circuit T5, T7 and T13, there is no activator protein synthesized thus the affinity rates are ignorable. We then set ${S}_{a}$ and/or ${S}_{u}$ value as very small in these circuits. Similarly, if the repressor genes are removed, the affinity rate for ${S}_{t}$ and ${S}_{l}$ will become very small. We also argue that the leakage levels are unchanged since we assume the leakage terms are determined by the hybrid promoter part only and all circuits in Figure 1—figure supplement 1A share the same hybrid promoters. So we can use the same model to represent all synthetic gene circuits in Figure 1—figure supplement 1A, by setting the value of $S$ to a very small number to represent the elimination of certain link.

We use the same method (Faucon et al., 2014; Ma et al., 2009) to quantitatively evaluate the ability to achieve multistability for different circuits. Stronger ability to generate multistability is defined as easier to achieve multistability in arbitrary dosage combinations. So we set a reasonable ranges for the activation/repression strengthes $S$ of the existing links in the circuit to represent the variance of inducer concentrations, and randomly pick up a parameter combination from the given parameter space and calculate whether the parameter combination can generate more than two stable steady states (SSS). We choose ${p}_{0}=0.1$, ${a}_{0}=0.1$, ${r}_{0}=0.1$ and $d=0.05$ in the numerical simulations. The ranges for activation/repression strength are $[0.3,0.8]$. If the link is missing in certain circuit, we set the value for the corresponding activation/repression strengths as very small number, e.g., $0.01$. We repeat this procedure for 2000 times for all sub-networks in Figure 1—figure supplement 1A and then summarize the probability to find a suitable parameter combination and also the distribution for these parameters. As they are shown in Figure 1—figure supplement 1H, only circuits with two auto-activation links, which are ${R}^{2}{A}^{2}$(T15), $R{A}^{2}$(T11 and T14) and ${A}^{2}$(T10), can exhibit tristability (black solid lines) and quadrastability(red solid lines). Only parameter combinations for bistability can be found for circuit with one autoactivation link, i.e., ${R}^{2}A$(T7 and T13). The probabilities to find parameter combination for tristability and quadrastability are given in Figure 1—figure supplement 1G. We also show the parameter distributions for these four circuits in Figure 1—figure supplement 1H, while the distribution for fixed parameters in these circuits are not shown. We can also find that ${R}^{2}{A}^{2}$ has wider ranges of parameters for tristability and quadrastability than the other two circuits.

The scatter plots for all the four parameters in ${R}^{2}{A}^{2}$ that have tristability or quadrastability are shown in Figure 1—figure supplement 1I. The scatter plot between ${S}_{a}$ and ${S}_{u}$ are symmetrical, which suggests that the activation strengths have to be balanced to generate multistability. It is the same for the inhibition strength of ${S}_{t}$ and ${S}_{l}$. We can also find strong correlation between the activation strength and inhibition strength, especially when activation is strong. This is also the reason that the range for activation strength is enlarged. This scatter plots suggest that specific inhibition strength is required to generate multistability when the circuit has strong expression, however, when the activation is weak the multistability is not sensitive to inhibition strength.

We perform a series of hysteresis experiments and bifurcation analysis, including single inducer and dual inducer, to fully explore the multiple stability behaviors of MINPA.

Before performing bifurcation analysis, we verify and adjust the quantitative model Equation (16), which is built on multimerization and dosage response curves using the transition dosage in hysteresis experiments for Arabinose and AHL in Figure 3—figure supplement 2I and Figure 3—figure supplement 2J. We also find that, the system can maintain both LL and HH states as they are shown in panel (C1${}_{LL}$) in Figure 3C and panel (C1${}_{HH}$) in Figure 3D. From these experimental results and modeling analysis, we can hypothesize that the MINPA system is quadrastable when no inducer is added.

The bifurcation analysis are performed using Matlab and bifurcation analysis package matcont. Firstly, we exhaustively search all steady states in the initial conditions. We divide the phase space into $10\times 10$ grids. Within each grid, we solve Equation (16) for the steady states with a random of initial solution located within this grid. Secondly, we initialize from each steady state and using parameter continuation package provided by matcont to go through this branch when the parameter increases. Similarly, we explore all branches initialized from all steady states. Finally, we use Jacobian method to study the stability of all the branch points and then mark the SSS to blue and USS to red.

We use the logarithm of inducer concentrations to perform bifurcation analysis. In Figure 3—figure supplement 2A, initially there are 4 SSSs when all the concentrations of inducers are zero (marked as blue solid lines). When Arabinose increases, the LL state vanishes and become mCherry state, and the GFP state will also become HH state. The HH and mCherry states remain stable even at very high dosage of Arabinose. The bifurcation suggests that we can drive the system from LL state to mCherry state using Arabinose. Also the increase of Arabinose can help the system overcome the barrier between GFP and HH state.

On the contrary to the bifurcation of Arabinose, the bifurcation of AHL suggests that the LL state will not vanish when the concentration of AHL increases, as it is shown in Figure 3—figure supplement 2B. Only the mCherry state will vanish so AHL can be used to drive the mCherry state to HH state. Meanwhile, the GFP state and HH state remain stable. The branches of USS will change as AHL increases, and it suggests that AHL can be used to modify the stability of associated SSSs.

In Figure 3—figure supplement 2C–D we can see that, all branches of SSSs will not disappear when the concentration of inducers of aTc and IPTG increase. So changing the concentration of aTc or IPTG solely cannot drive the system from one stable state to the other one. Similarly to AHL, the USS branches will change their position thus the stability of SSSs will change accordingly.

From these single induction experiment, we can see there are many limitations in driving the initial LL state to all the other three states. The Arabinose can drive the system out of low-mc states. However, AHL, aTc and IPTG can not drive the system out of LL state solely by themselves. They can only change the stability of LL state. For these reasons, we need to combine different dosage to efficiently change the landscape of SSSs and further control the differentiation.

We design a dual induction approach to study the bifurcation behaviors in two dimensional parameter space. Given two changing parameters $[{D}_{1}]$ and $[{D}_{2}]$ and GFP and mCherry, it will be challenging to visualize and interpret results in four dimensions. Here we use a hybrid bifurcation technique, in which the concentration ratio $\beta =[{D}_{1}]/[{D}_{2}]$ between two changing parameters is fixed, to explore the parameter space along different directions. The hybrid bifurcation returns to the conventional one-dimensional bifurcation when $\beta =0$ or $\beta =\mathrm{\infty}$.

In Figure 3C, we can see that, by choosing the dual inducers as Arabinose and AHL and $\beta =10$, the MINPA system initialized from LL state will transit to mCherry and then HH state. In low mix dosage concentration, the system will stay at LL state and start to transit to either mCherry or GFP state when the mix concentration increases. Finally, all of the systems will rest in HH state. Different value of $\beta $ can have slightly different transition behaviors, as they are shown in Figure 3—figure supplement 2E–H. When the value of $\beta $ increases from $0.1$ to $100$, the system would prefer to transit to the GFP state more than the mCherry state. We decide to choose $\beta =10$ because the MINPA will have wider quadrastable region than the others.

We study the effect of sequential inductions using bifurcation analysis and stochastic simulations. Compared to the hybrid bifurcation which apply the inducers at the same time, the sequential induction applies one inducer first and then the other inducer after a duration of $T$. The essence of sequential induction is that the induction effect depends on the SSS which the system stays in. For example, if the system is already at high-mC states, the inhibition from mCherry to GFP would create extra obstacle for the system to transit from GFP low to GFP high state. For this reason, it will require higher concentration of inducers to drive the system from mCherry state to HH state than the process from LL state to GFP state.

In deterministic systems, as it is shown in Figure 2A, different induction sequences may drive the gene network to different terminal states. We perform detailed bifurcation analysis and find that the interaction between inducer one and inducer two on their bistable regions can lead to different terminal states. The two dimensional phase diagram is show in Figure 2A. When the concentration of inducer one increases, the bistable region for inducer two will become broader. Also, its lower boundary increases, which suggest that the toggle will need more inducer two to maintain its high GFP response state. It is similar for the increase of inducer two. We assume that the system initialize at high GFP state and neither concentration of inducer one nor inducer two can cause state transition to low GFP state. If inducer two is introduced first, the bistable region for inducer one will change and the required concentration of inducer two to induce state transition will become smaller because the lower boundary of its bistable region increase. So a state transition will happen because of the sequence of induction. However, if sequence is permuted or two inducers are applied simultaneously, the system remain in high state. We can also learn from the two dimensional bifurcation diagram that it will have larger chance to discover the state transition boundary, comparing to using the hybrid induction method alone (the diagonal induction curve).

We now show how to carefully perform experiment to observe the sequential induction effects based on numerical calculations. First we constructed *E. Coli* plasmids which contain mutual inhibitory circuit and can be induced by aTc and IPTG, as it is shown in Figure 2—figure supplement 1A. Then we chose two dosages of IPTG, ${D}_{I}L$ and ${D}_{I}H$ so that the system could be bistable for aTc under these two different IPTG concentrations. From previous numerical studies, we found that as long as the right boundary of bistable region, ${B}_{a}H$ under high IPTG was shifted to the right of the one ${B}_{a}L$ under low IPTG, the sequential effect can be observed by setting the high aTc concentration ${D}_{a}H$ between $[{B}_{a}L,{B}_{a}H]$. The low aTc concentration ${D}_{a}L$ shall be less than ${B}_{a}L$. In our experiment shown in Figure 2—figure supplement 1B, we chose ${D}_{I}H$ as $8\times {10}^{-5}M$, ${D}_{I}L$ as 0, ${D}_{a}H$ as $100$ ng/ml, and ${D}_{a}L$ as 0.

Different sequences of induction provide more alternatives to explore the parameter space from different directions and can help to generate different population distributions in stochastic systems with appropriate low and high inducer concentrations, given finite amount of time. We can also use this property to design new cell differentiation protocols using this method. The transition rate between two neighboring SSSs depends on the barrier height between them. For example, in order to induce the cells from LL state to GFP state, we need to jointly add aTc and AHL to lower the barrier height between LL and GFP state. Some barrier can be very high thus it is impossible for the system to overcome the barrier under given inducers, e.g., the one between mCherry state to HH state. We can also use the sequential induction to compare the multistability in varying circuits in Figure 1—figure supplement 1A.

Because of the nonlinearity of gene regulation and inducer response functions, the constructed gene networks is highly dissipative thus a conventional potential can not be defined (Wang et al., 2011). For this reason, we construct a quasi-potential landscape to describe the state transition barriers and stability of all SSSs. First, we made the assumption that all the gene network described by these dynamical functions will approach stationary state, where the densities of final state is a constant. We can define the pseudo potential from the stationary densities $P(x)$, and can be defined as $E(x)=-logP(x)$. In the conventional canonical dynamical system, the descending direction of a landscape is always pointing to the minimum of the landscape, however, it is not same in this dissipative system. Also, because of the high dimensional state space and multiple attractor geometry, the convergence to stationary distribution is always slow under the transitional random walk algorithm. In order to address this obstacle, we applied the so-called weight ensemble random walk method to speed the convergence rates (Kromer et al., 2013).

We describe in detail the procedure of calculating the stable density distribution under the dynamical function Equation 16 as follow. First, we need to determine the time step $\tau $ and the space resolution $\mathrm{\Delta}x$ according to reference (Kromer et al., 2013) to ensure convergence. When the noise strength is $D$, the time step of simulation $\tau $ shall satisfy: $\tau <\frac{4D}{max{f}^{2}(\mathbf{\mathbf{X}})}$, and $f(\mathbf{\mathbf{X}})$ is the velocity field. Also, $\mathrm{\Delta}x$ shall satisfy $\mathrm{\Delta}x\le f(\mathbf{\mathbf{X}})\tau $. Then we can discretize the state space, or possible ranges of protein concentration, into $M\times M$ lattices as $M=1/\mathrm{\Delta}x$. We will use the densities of each lattice ${P}_{m,n}$ to approximate densities $P(x)$ in the state space. The initial probability ${P}_{m,n}(0)$ of all gird points are set to be uniform, as ${P}_{m,n}(0)=\frac{1}{M\times M}$. When the simulation is long enough, the initial distribution would not affact the final stationary distribution.

At each step of duration $\tau $, we randomly place $N$ walkers within each grid $[m,n]$, no matter how small its probability is. Each of these walkers carries equal weight ${q}_{m,n}^{i}(t)={P}_{m,n}(t)/N$ and start evloving from its initial position under the system dynamics and noise, and ‘transport’ to nearby grids. Its evolution can be simulated using any numerical integration of Equation 16. Here we choose the stochastic version of second order Runge-Kutta algorithm, or the Heun algorithm. When all of the walkers’ position have been updated, the new probabilities for each grid is:

$${P}_{m,n}(t+\tau )=\sum _{k\in {G}_{m,n}}{q}_{{m}^{\prime},{n}^{\prime}}^{k}(t),$$

(18)

where ${q}_{{m}^{\prime},{n}^{\prime}}^{k}(t)$ is the probability that carried by a walker initialized at grid $[{m}^{\prime},{n}^{\prime}]$ and fell into grid ${G}_{m,n}$. At next time step, probabilities carried by another $N$ new walkers will be updated according to the new updated probability. This procedure repeats until the probability distribution ${P}_{m,n}(t)$ becomes stationary. The landscape can be calculated from ${E}_{m,n}=-log{P}_{m,n}$.

Using potential landscape we can estimate the relative stability of different SSSs and the hardness, or the barrier heights for the transition between two SSSs.

Microfluidics coupled with time-lapse imaging was employed to visualize the state transitions at the single-cell level. Two media were prepared: one with inducers (Arabinose, AHL and aTc) and the other without. After cells were loaded into the trap (1 to 5 cells for a trap), the device was heated up to 37 degree and cells were supplied with LB media without inducers for 6 hr. Sulforhodamine was added as a dye to monitor nutrient transport. Then, the supplied media was switched to the media added with Arabinose, AHL and aTc for another 18 hr. Media switching was controlled by adjusting the heights of the medium syringes relative to one another.

However, cells treated with the three inducers demonstrate symptoms of significant stress and cell death, presumably due to photo toxicity compounded with flow-induced sheer stress and other mechanical stresses in the microenvironment (Kohles et al., 2009; Shen et al., 2014; Shemesh et al., 2015). Lower concentrations of three inducer combinations were also tested but yield no significant improvement of cells viability.

Since the logic of emergence of quadrastability is enhancing the two positive feedbacks of MINPA through adding inducers Arabinose, AHL and aTc, quadrastability could also be achieved through weakening the mutual inhibition using IPTG and aTc. Depending on the basal expression of two hybrid promoters, IPTG and aTc can promote GFP and mCherry expression to a limited extent, which in turn attenuates fluorescent proteins toxicity. Hence, we tried to use IPTG and aTc to induce quadrastability instead of the three inducers tried in flow cytometry. Experimental result showed that the initial low-low state cells could differentiate into GFP, mCherry and high-high state cells with $2\times {10}^{-4}$ M IPTG and 200 ng/ml aTc induction (Figure 4—figure supplement 1D and Appendix 1—video 1). It is interesting that the trajectory for many cells were from GFP to high-high to mCherry state. Altogether, this result further verified the MINPA has the potential to generate quadrastability in living cells.

Click here to view.^{(1.0M, mp4)}

Cells with the MINPA plasmid were cultured overnight at 37 degree and diluted into fresh media with corresponding inducers at 1:100 ratio (O.D. $\approx 0.066$). The four individual inducers are Ara ($2.5\times {10}^{-5}$m/v), AHL ($1\times {10}^{-5}$ M), aTc ($200$ ng/ml), IPTG ($1\times {10}^{-4}$ M), and inducer combinations: AHL and aTc, AHL, aTc, and Ara. Cellular growth rates were measured by using $200\mu $L cultures in a 96-well plate with absorbance at $600$ nm on a plate reader (BioTek, USA). Three replicates were tested for each condition.

Compared to AHL or Arabinose individually, aTc addition (aTc, AHL + aTc, and Ara + AHL + aTc) influenced cells growth and increased the lag phase for about $2.5$ hours (Appendix 1—figure 1). But at about 13 hr, the growth rates are almost the same. Since the timescale in our experiments is much longer ($\ge 24$ h) than 13 hr, we reason the effect of inducers on the cellular growth rates would not change interpretations of experimental observations.

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

This paper was supported by the following grants:

- http://dx.doi.org/10.13039/100000968American Heart Association 15PRE25710303 to Fuqing Wu.
- http://dx.doi.org/10.13039/100000183Army Research Office W911NF-14-1-0504 to Ying-Cheng Lai.
- http://dx.doi.org/10.13039/100000001National Science Foundation DMS-1100309 to Xiao Wang.
- http://dx.doi.org/10.13039/100000002National Institutes of Health GM106081 to Xiao Wang.
- http://dx.doi.org/10.13039/100000968American Heart Association 11BGIA7440101 to Xiao Wang.

The authors declare that no competing interests exist.

FW, Designed the research, Performed the experiments, Wrote the paper.

R-QS, Designed the research, Developed the mathematical modeling, Wrote the paper.

Y-CL, Wrote the paper.

XW, Designed the research, Wrote the paper.

10.7554/eLife.23702.015#### Source code file 1.

**Flow cytometry data analysis file.**

**DOI:**
http://dx.doi.org/10.7554/eLife.23702.015

Click here to view.^{(13K, zip)}

10.7554/eLife.23702.016#### Source code file 2.

**Bifurcation analysis in Figure 3.**

**DOI:**
http://dx.doi.org/10.7554/eLife.23702.016

Click here to view.^{(13K, zip)}

10.7554/eLife.23702.017#### Source code file 3.

**Potential landscape calculation file.**

**DOI:**
http://dx.doi.org/10.7554/eLife.23702.017

Click here to view.^{(13K, zip)}

- Acar M, Becskei A, van Oudenaarden A. Enhancement of cellular memory by reducing stochastic transitions. Nature. 2005;435:228–232. doi: 10.1038/nature03524. [PubMed] [Cross Ref]
- Angeli D, Ferrell JE, Sontag ED. Detection of Multistability, Bifurcations, and hysteresis in a large class of biological positive-feedback systems. PNAS. 2004;101:1822–1827. doi: 10.1073/pnas.0308265100. [PubMed] [Cross Ref]
- Arkin A, Ross J, McAdams HH. Stochastic kinetic analysis of developmental pathway bifurcation in phage lambda-infected Escherichia coli cells. Genetics. 1998;149:1633–1648. [PubMed]
- Balázsi G, van Oudenaarden A, Collins JJ. Cellular decision making and biological noise: from microbes to mammals. Cell. 2011;144:910–925. doi: 10.1016/j.cell.2011.01.030. [PMC free article] [PubMed] [Cross Ref]
- Bennett MR, Pang WL, Ostroff NA, Baumgartner BL, Nayak S, Tsimring LS, Hasty J. Metabolic gene regulation in a dynamically changing environment. Nature. 2008;454:1119–1122. doi: 10.1038/nature07211. [PMC free article] [PubMed] [Cross Ref]
- Bessonnard S, De Mot L, Gonze D, Barriol M, Dennis C, Goldbeter A, Dupont G, Chazaud C. Gata6, nanog and erk signaling control cell fate in the inner cell mass through a tristable regulatory network. Development. 2014;141:3637–3648. doi: 10.1242/dev.109678. [PubMed] [Cross Ref]
- Bhattacharya S, Zhang Q, Andersen ME. A deterministic map of Waddington's epigenetic landscape for cell fate specification. BMC Systems Biology. 2011;5:85 doi: 10.1186/1752-0509-5-85. [PMC free article] [PubMed] [Cross Ref]
- Chalancon G, Ravarani CN, Balaji S, Martinez-Arias A, Aravind L, Jothi R, Babu MM. Interplay between gene expression noise and regulatory network architecture. Trends in Genetics. 2012;28:221–232. doi: 10.1016/j.tig.2012.01.006. [PMC free article] [PubMed] [Cross Ref]
- Faucon PC, Pardee K, Kumar RM, Li H, Loh YH, Wang X. Gene networks of fully connected triads with complete auto-activation enable multistability and stepwise stochastic transitions. PLoS One. 2014;9:e102873 doi: 10.1371/journal.pone.0102873. [PMC free article] [PubMed] [Cross Ref]
- Feng H, Han B, Wang J. Adiabatic and non-adiabatic non-equilibrium stochastic dynamics of single regulating genes. The Journal of Physical Chemistry B. 2011;115:1254–1261. doi: 10.1021/jp109036y. [PubMed] [Cross Ref]
- Feng H, Wang J. A new mechanism of stem cell differentiation through slow binding/unbinding of regulators to genes. Scientific Reports. 2012;2:550 doi: 10.1038/srep00550. [PMC free article] [PubMed] [Cross Ref]
- Ferrell JE. Bistability, bifurcations, and Waddington's epigenetic landscape. Current Biology. 2012;22:R458–R466. doi: 10.1016/j.cub.2012.03.045. [PMC free article] [PubMed] [Cross Ref]
- Ferry MS, Razinkov IA, Hasty J. Microfluidics for synthetic biology: from design to execution. Methods in Enzymology. 2011;497:295–372. doi: 10.1016/B978-0-12-385075-1.00014-7. [PMC free article] [PubMed] [Cross Ref]
- Gardner TS, Cantor CR, Collins JJ. Construction of a genetic toggle switch in Escherichia coli. Nature. 2000;403:339–342. doi: 10.1038/35002131. [PubMed] [Cross Ref]
- Graf T, Enver T. Forcing cells to change lineages. Nature. 2009;462:587–594. doi: 10.1038/nature08533. [PubMed] [Cross Ref]
- Guantes R, Poyatos JF. Multistable decision switches for flexible control of epigenetic differentiation. PLoS Computational Biology. 2008;4:e1000235 doi: 10.1371/journal.pcbi.1000235. [PMC free article] [PubMed] [Cross Ref]
- Huang S, Guo YP, May G, Enver T. Bifurcation dynamics in lineage-commitment in bipotent progenitor cells. Developmental Biology. 2007;305:695–713. doi: 10.1016/j.ydbio.2007.02.036. [PubMed] [Cross Ref]
- Jolly MK, Tripathi SC, Jia D, Mooney SM, Celiktas M, Hanash SM, Mani SA, Pienta KJ, Ben-Jacob E, Levine H, Jolly MK, Tripathi SC, Jia D, Mooney SM, Celiktas M, Hanash SM, Mani SA, Pienta KJ, Ben-Jacob E, Levine H. Stability of the hybrid epithelial/mesenchymal phenotype. Oncotarget. 2016;7:27067–27084. doi: 10.18632/oncotarget.8166. [PMC free article] [PubMed] [Cross Ref]
- Kauffman SA. The Origins of Order: Self-Organization and Selection in Evolution. New York: Oxford University Press; 1993.
- Kohles SS, Nève N, Zimmerman JD, Tretheway DC. Mechanical stress analysis of microfluidic environments designed for isolated biological cell investigations. Journal of Biomechanical Engineering. 2009;131:121006 doi: 10.1115/1.4000121. [PMC free article] [PubMed] [Cross Ref]
- Kramer BP, Fussenegger M. Hysteresis in a synthetic mammalian gene network. PNAS. 2005;102:9517–9522. doi: 10.1073/pnas.0500345102. [PubMed] [Cross Ref]
- Kromer JA, Schimansky-Geier L, Toral R. Weighted-ensemble brownian dynamics simulation: sampling of rare events in nonequilibrium systems. Physical Review E. 2013;87:063311 doi: 10.1103/PhysRevE.87.063311. [PubMed] [Cross Ref]
- Kroon E, Martinson LA, Kadoya K, Bang AG, Kelly OG, Eliazer S, Young H, Richardson M, Smart NG, Cunningham J, Agulnick AD, D'Amour KA, Carpenter MK, Baetge EE. Pancreatic endoderm derived from human embryonic stem cells generates glucose-responsive insulin-secreting cells in vivo. Nature Biotechnology. 2008;26:443–452. doi: 10.1038/nbt1393. [PubMed] [Cross Ref]
- Laurent M, Kellershohn N. Multistability: a major means of differentiation and evolution in biological systems. Trends in Biochemical Sciences. 1999;24:418–422. doi: 10.1016/S0968-0004(99)01473-5. [PubMed] [Cross Ref]
- Lee MJ, Ye AS, Gardino AK, Heijink AM, Sorger PK, MacBeath G, Yaffe MB. Sequential application of anticancer drugs enhances cell death by rewiring apoptotic signaling networks. Cell. 2012;149:780–794. doi: 10.1016/j.cell.2012.03.031. [PMC free article] [PubMed] [Cross Ref]
- Lee J, Lee J, Farquhar KS, Yun J, Frankenberger CA, Bevilacqua E, Yeung K, Kim EJ, Balázsi G, Rosner MR. Network of mutually repressive metastasis regulators can promote cell heterogeneity and metastatic transitions. PNAS. 2014a;111:E364–E373. doi: 10.1073/pnas.1304840111. [PubMed] [Cross Ref]
- Lee J, Tiwari A, Shum V, Mills GB, Mancini MA, Igoshin OA, Balázsi G. Unraveling the regulatory connections between two controllers of breast Cancer cell fate. Nucleic Acids Research. 2014b;42:6839–6849. doi: 10.1093/nar/gku360. [PMC free article] [PubMed] [Cross Ref]
- Li C, Wang J. Quantifying cell fate decisions for differentiation and reprogramming of a human stem cell network: landscape and biological paths. PLoS Computational Biology. 2013a;9:e1003165 doi: 10.1371/journal.pcbi.1003165. [PMC free article] [PubMed] [Cross Ref]
- Li C, Wang J. Quantifying waddington landscapes and paths of non-adiabatic cell fate decisions for differentiation, reprogramming and transdifferentiation. Journal of the Royal Society Interface. 2013b;10:20130787. doi: 10.1098/rsif.2013.0787. [PMC free article] [PubMed] [Cross Ref]
- Litcofsky KD, Afeyan RB, Krom RJ, Khalil AS, Collins JJ. Iterative plug-and-play methodology for constructing and modifying synthetic gene networks. Nature Methods. 2012;9:1077–1080. doi: 10.1038/nmeth.2205. [PMC free article] [PubMed] [Cross Ref]
- Liu X, Sun H, Qi J, Wang L, He S, Liu J, Feng C, Chen C, Li W, Guo Y, Qin D, Pan G, Chen J, Pei D, Zheng H. Sequential introduction of reprogramming factors reveals a time-sensitive requirement for individual factors and a sequential EMT-MET mechanism for optimal reprogramming. Nature Cell Biology. 2013;15:829–838. doi: 10.1038/ncb2765. [PubMed] [Cross Ref]
- Ma W, Trusina A, El-Samad H, Lim WA, Tang C. Defining network topologies that can achieve biochemical adaptation. Cell. 2009;138:760–773. doi: 10.1016/j.cell.2009.06.013. [PMC free article] [PubMed] [Cross Ref]
- Maamar H, Raj A, Dubnau D. Noise in gene expression determines cell fate in Bacillus subtilis. Science. 2007;317:526–529. doi: 10.1126/science.1140818. [PMC free article] [PubMed] [Cross Ref]
- MacArthur BD, Please CP, Oreffo RO. Stochasticity and the molecular mechanisms of induced pluripotency. PLoS One. 2008;3:e3086 doi: 10.1371/journal.pone.0003086. [PMC free article] [PubMed] [Cross Ref]
- Macarthur BD, Ma'ayan A, Lemischka IR. Systems biology of stem cell fate and cellular reprogramming. Nature Reviews Molecular Cell Biology. 2009;10:672–681. doi: 10.1038/nrm2766. [PMC free article] [PubMed] [Cross Ref]
- Murphy KF, Adams RM, Wang X, Balázsi G, Collins JJ. Tuning and controlling gene expression noise in synthetic gene networks. Nucleic Acids Research. 2010;38:2712–2726. doi: 10.1093/nar/gkq091. [PMC free article] [PubMed] [Cross Ref]
- Narula J, Smith AM, Gottgens B, Igoshin OA. Modeling reveals bistability and low-pass filtering in the network module determining blood stem cell fate. PLoS Computational Biology. 2010;6:e1000771 doi: 10.1371/journal.pcbi.1000771. [PMC free article] [PubMed] [Cross Ref]
- Niwa H, Toyooka Y, Shimosato D, Strumpf D, Takahashi K, Yagi R, Rossant J. Interaction between Oct3/4 and Cdx2 determines trophectoderm differentiation. Cell. 2005;123:917–929. doi: 10.1016/j.cell.2005.08.040. [PubMed] [Cross Ref]
- Oppenheim AB, Kobiler O, Stavans J, Court DL, Adhya S. Switches in bacteriophage lambda development. Annual Review of Genetics. 2005;39:409–429. doi: 10.1146/annurev.genet.39.073003.113656. [PubMed] [Cross Ref]
- Pagliuca FW, Millman JR, Gürtler M, Segel M, Van Dervort A, Ryu JH, Peterson QP, Greiner D, Melton DA. Generation of functional human pancreatic β cells in vitro. Cell. 2014;159:428–439. doi: 10.1016/j.cell.2014.09.040. [PMC free article] [PubMed] [Cross Ref]
- Palani S, Sarkar CA. Integrating extrinsic and intrinsic cues into a minimal model of lineage commitment for hematopoietic progenitors. PLoS Computational Biology. 2009;5:e1000518 doi: 10.1371/journal.pcbi.1000518. [PMC free article] [PubMed] [Cross Ref]
- Paşca AM, Sloan SA, Clarke LE, Tian Y, Makinson CD, Huber N, Kim CH, Park JY, O'Rourke NA, Nguyen KD, Smith SJ, Huguenard JR, Geschwind DH, Barres BA, Paşca SP. Functional cortical neurons and astrocytes from human pluripotent stem cells in 3D culture. Nature Methods. 2015;12:671–678. doi: 10.1038/nmeth.3415. [PMC free article] [PubMed] [Cross Ref]
- Pomerening JR, Sontag ED, Ferrell JE. Building a cell cycle oscillator: hysteresis and bistability in the activation of Cdc2. Nature Cell Biology. 2003;5:346–351. doi: 10.1038/ncb954. [PubMed] [Cross Ref]
- Prindle A, Selimkhanov J, Li H, Razinkov I, Tsimring LS, Hasty J. Rapid and tunable post-translational coupling of genetic circuits. Nature. 2014;508:387–391. doi: 10.1038/nature13238. [PMC free article] [PubMed] [Cross Ref]
- Rabajante JF, Babierra AL. Branching and oscillations in the epigenetic landscape of cell-fate determination. Progress in Biophysics and Molecular Biology. 2015;117:240–249. doi: 10.1016/j.pbiomolbio.2015.01.006. [PubMed] [Cross Ref]
- Schmiedel JM, Klemm SL, Zheng Y, Sahay A, Blüthgen N, Marks DS, van Oudenaarden A. Gene expression. MicroRNA control of protein expression noise. Science. 2015;348:128–132. doi: 10.1126/science.aaa1738. [PubMed] [Cross Ref]
- Schultz D, Wolynes PG, Ben Jacob E, Onuchic JN. Deciding fate in adverse times: sporulation and competence in Bacillus subtilis. PNAS. 2009;106:21027–21034. doi: 10.1073/pnas.0912185106. [PubMed] [Cross Ref]
- Shemesh J, Jalilian I, Shi A, Heng Yeoh G, Knothe Tate ML, Ebrahimi Warkiani M, Yeoh GH, Tate MLK, Warkiani ME. Flow-induced stress on adherent cells in microfluidic devices. Lab Chip. 2015;15:4114–4127. doi: 10.1039/C5LC00633C. [PubMed] [Cross Ref]
- Shen F, Li X, Li PC. Study of flow behaviors on single-cell manipulation and shear stress reduction in microfluidic chips using computational fluid dynamics simulations. Biomicrofluidics. 2014;8:014109 doi: 10.1063/1.4866358. [PubMed] [Cross Ref]
- Stricker J, Cookson S, Bennett MR, Mather WH, Tsimring LS, Hasty J. A fast, robust and tunable synthetic gene oscillator. Nature. 2008;456:516–519. doi: 10.1038/nature07389. [PubMed] [Cross Ref]
- Szathmáry E, Jordán F, Pál C. Molecular biology and evolution. can genes explain biological complexity? Science. 2001;292:1315–1316. doi: 10.1126/science.1060852. [PubMed] [Cross Ref]
- Süel GM, Garcia-Ojalvo J, Liberman LM, Elowitz MB. An excitable gene regulatory circuit induces transient cellular differentiation. Nature. 2006;440:545–550. doi: 10.1038/nature04588. [PubMed] [Cross Ref]
- Tanouchi Y, Pai A, Park H, Huang S, Stamatov R, Buchler NE, You L. A noisy linear map underlies oscillations in cell size and gene expression in Bacteria. Nature. 2015;523:357–360. doi: 10.1038/nature14562. [PMC free article] [PubMed] [Cross Ref]
- Waddington CH. The Strategy of the Genes; a Discussion of Some Aspects of Theoretical Biology. London: Allen & Unwin; 1957.
- Wang J, Xu L, Wang E. Potential landscape and flux framework of nonequilibrium networks: robustness, dissipation, and coherence of biochemical oscillations. PNAS. 2008;105:12271–12276. doi: 10.1073/pnas.0800579105. [PubMed] [Cross Ref]
- Wang J, Zhang K, Xu L, Wang E. Quantifying the Waddington landscape and biological paths for development and differentiation. PNAS. 2011;108:8257–8262. doi: 10.1073/pnas.1017017108. [PubMed] [Cross Ref]
- Wang J. Landscape and flux theory of non-equilibrium dynamical systems with application to biology. Advances in Physics. 2015;64:1–137. doi: 10.1080/00018732.2015.1037068. [Cross Ref]
- Wu M, Su RQ, Li X, Ellis T, Lai YC, Wang X. Engineering of regulated stochastic cell fate determination. PNAS. 2013;110:10610–10615. doi: 10.1073/pnas.1305423110. [PubMed] [Cross Ref]
- Wu F, Menn DJ, Wang X. Quorum-sensing crosstalk-driven synthetic circuits: from unimodality to trimodality. Chemistry & Biology. 2014;21:1629–1638. doi: 10.1016/j.chembiol.2014.10.008. [PMC free article] [PubMed] [Cross Ref]
- Xu L, Zhang K, Wang J. Exploring the mechanisms of differentiation, dedifferentiation, reprogramming and transdifferentiation. PLoS One. 2014;9:e105216 doi: 10.1371/journal.pone.0105216. [PMC free article] [PubMed] [Cross Ref]
- Zhang B, Wolynes PG. Stem cell differentiation as a many-body problem. PNAS. 2014;111:10185–10190. doi: 10.1073/pnas.1408561111. [PubMed] [Cross Ref]

- eLife. 2017; 6: e23702. »
- Decision letter

2017; 6: e23702.

Wenying Shou, Reviewing editor

Wenying Shou, Fred Hutchinson Cancer Research Center, United States;

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "Engineering of a synthetic quadrastable gene network to approach Waddington landscape and cell fate determination" for consideration by *eLife*. Your article has been favorably evaluated by Naama Barkai (Senior Editor) and three reviewers, one of whom, Wenying Shou (Reviewer #1), is a member of our Board of Reviewing Editors. The following individual involved in review of your submission has agreed to reveal their identity: Gabor Balazsi (Reviewer #2).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

To quantitatively and experimentally understand Waddington's landscape analogy commonly used to depict multiple cell fates, authors constructed a genetic network called MINPA (mutual inhibition and positive autoactivation). By analyzing various sub-networks of MINPA and the full MINPA network, authors find that the full MINPA is capable of achieving the most fate states of mCherry/GFP expression patterns. A total of four states can be achieved, with low or high mCherry and GFP. Aided by modeling, the authors charted how the "energy" landscape of cell fates changes in the presence of various chemicals, and was able to "guide" cells to transit from one state to another by applying chemicals at the correct orders.

All three reviewers thought that your work was important and interesting. However, there are issues that will need to be clarified and addressed before acceptance. Details are provided in the reviewers' comments below. Please note that depending on your response, reviewers might conclude that additional major experiments could be necessary. To help us assess the likelihood of such an outcome, it would help if you responded to the critiques with a plan of action and a timetable for completion of the work recommended by the reviewers. The editor and reviewers will assess your response and provide a recommendation to assist you in preparing a revised submission.

Reviewer #1:

I enjoyed reading it as someone with an informal interest in landscape.

Figure 3B: Not clear to me that at C2, the mCherry state (red dot) has low GFP. Need to make the 3-dimensional illustration more clear.

Figure 3D: How long can the HH state be maintained? Does your model also predict the time scale of decay?

Needs to better describe how "potential" in Figure 4 is calculated.

Reviewer #2:

This is an important and interesting manuscript that investigates experimentally the possibility of multistability in a library of two-gene networks with increasingly complex connectivity consisting of mutual inhibition and positive autoregulation. This topology has counterparts in mammalian cell fate-regulatory networks. The library consists of several synthetic networks, with a number of regulatory links that systematically increases from 2 to 4. Analyzing the dynamics of these networks reveals that the one with all 4 links present has the highest chance for multistability, with more than 2 stable steady states for a broad range of parameters. These findings are experimentally validated with clever sequential induction experiments and hysteresis using 4 different inducers.

The manuscript represents an important step forward in synthetic biology and beyond, having relevance to many cell fate determination networks, including in higher organisms. Nonetheless, some details are clarifications are still needed.

Therefore, I recommend publication once the following comments can be addressed:

1) Synthetic biology often measures its progress by the number of genes in synthetic gene circuits. However, this view ignores the fact that biological complexity does not really correlate with gene number. In fact, biological complexity correlates with the number of regulatory connections. Some plants have more genes than us, humans. In fact, a 10-gene network may be much simpler from a dynamical or an information-processing perspective than a two-gene network – provided that the latter has more complex connectivity. The arguments on network complexity resulting from links rather than nodes may be a nice addition to the Discussion. See Science 292(5520):1315-6 (2001).

2) From the Methods it seems like these networks are on high-copy plasmids rather than integrated. What is the effect of plasmid copy number variation on the results? Could the plasmid be lost from some cells, generating the impression of multistability? Flow cytometry at multiple time points will provide at least a partial answer.

3) The sequential induction method is quite interesting and important for revealing multistability. However, its application to the T15 network (or other networks) could be explained better. I guess the sequence in which inducers are added determines how the network's dynamics changes from monostable to multistable – but the details may be important. If bistable, tristable, quadristable regions are small in the 4-dimensional space then the goal is to "hit" these regions as inducers are sequentially added. While this is shown for the toggle switch, it is unclear how was this done for the more complex networks. This should be much more clearly described, especially for T15 and possibly for some simpler networks. Overall, going from the toggle switch to T15 and even to simpler networks is a big step, which requires more explanation and investigation. In addition, to toggle switch would be worth adding to the experimentally measured networks as control, considering that its behavior is computationally predictable. Do the experiments validate the predictions in Figure 2A for the toggle switch?

Reviewer #3:

In this work, the authors quantified the cell fate decision making landscape through engineering the synthetic circuits. They have identified multiple states. They also considered the cell fate decision at different conditions. The results are interesting and I recommend its publication after revisions upon the following comments.

1) Recently, the theoretical advances have been able to identify the landscape and flux as the driving force of the global dynamics of the biological circuits (Proc. Natl. Acad. Sci. USA, 105: 12271-12276. (2008).). The quantification of the Waddington landscape has been achieved and is directly related to the underlying gene regulatory network for cell fate decision of stem cell differentiation and development (Proc. Natl. Acad. Sci. USA. 108(20):8257-8262(2011).). These should be reflected in the revision.

2) Furthermore, the multiple states have been predicted beyond the bistable states/switches without and with the epigenetic effects (Proc. Natl. Acad. Sci. USA. 108(20):8257-8262(2011); J. Phys. Chem. B, 115, 1254 (2011); Sci. Rep., 2,550 (2012); J R Soc Interface 10: 20130787 (2013); PLoS ONE. 9(8): e105216 (2014); Advances in Physics, 64:1, 1-137. (2015).). These should be reflected in the revision.

3) The authors should explain more clearly the physical and biological origins of the quadrastable states. Why the triple stable states are less likely than the quadrastable states? What is the origin of the LL state?

- eLife. 2017; 6: e23702. »
- Author response

2017; 6: e23702.

[…] Reviewer #1:

I enjoyed reading it as someone with an informal interest in landscape.

Figure 3B: Not clear to me that at C2, the mCherry state (red dot) has low GFP. Need to make the 3-dimensional illustration more clear.

We thank the reviewer for this great comment. To make the illustration clearer, especially for GFP at C2, we rotated the diagram and got a better view of the states at C2, which is added as Figure 3—figure supplement 1.

Figure 3D: How long can the HH state be maintained? Does your model also predict the time scale of decay?

In our experiments, we found the HH state is very stable. When we tried to use AHL and Arabinose to induce MINPA, the HH state could be maintained up to 24 hours after induction. In the hysteresis results (Figure 3D and Figure 3—figure supplement 2J), the HH state cells were washed and resuspended in fresh media with Arabinose and AHL, and the state is stable for at least 24 hours.

In our deterministic and stochastic models, protein degradation is the main factor determining time scales of state transitions. Once the system reaches its stable steady states, it will not change states until new induction is applied. Therefore, our model could qualitatively represent speed of state transitions by choosing appropriate protein degradation rate.

Needs to better describe how "potential" in Figure 4 is calculated.

We thank the reviewer for prompting us to improve the manuscript. As suggested, we summarized the method to calculate the “potential” in the revised main text (subsection “Experimental demonstration of model-guided quadrastability of MINPA”, first paragraph), and more details can be found in the Appendix (“Construct Quasi-potential Attractor Landscape”).

Reviewer #2:

[…] 1) Synthetic biology often measures its progress by the number of genes in synthetic gene circuits. However, this view ignores the fact that biological complexity does not really correlate with gene number. In fact, biological complexity correlates with the number of regulatory connections. Some plants have more genes than us, humans. In fact, a 10-gene network may be much simpler from a dynamical or an information-processing perspective than a two-gene network – provided that the latter has more complex connectivity. The arguments on network complexity resulting from links rather than nodes may be a nice addition to the Discussion. See Science 292(5520):1315-6 (2001).

We thank the reviewer for this great suggestion. We have added it into the second paragraph of our Discussion section.

2) From the Methods it seems like these networks are on high-copy plasmids rather than integrated. What is the effect of plasmid copy number variation on the results? Could the plasmid be lost from some cells, generating the impression of multistability? Flow cytometry at multiple time points will provide at least a partial answer.

We thank the reviewer for this question. We constructed these networks on high-copy plasmids. It is true that the plasmid copy number in different cells harboring the same network and backbone may not be the same, resulting in varied gene expression levels within a cell population. Moreover, stochasticity in gene expression arising from fluctuations in transcription and translation would also lead to different levels of gene expressions. We think these factors are two of the most important ones causing noisy gene expression. Experimentally, the consequence of plasmid copy number and gene expression stochasticity is reflected on the size of each population. For example, the Low-Low state cells in Figure 3A is a population with green fluorescence ranging from 3*10^{2} to 4*10^{3} (a.u.), and GFP state cells from 10^{4} to 2*10^{5} (a.u.). If all the cells have the non-fluctuating gene expression, we would expect a much smaller cloud on the 2-D fluorescence plots.

As the reviewer commented, it is possible that the plasmid may be lost from some cells, especially for long-term culture in which the active antibiotic concentration decreases over time. The cells lost their plasmid would have no fluorescent protein expression in a long time course, and might be only mixed/hidden with the low-low state cells. So it would not influence cells in the other three states (GFP, mCherry, and high-high states). A time-course induction with AHL, aTc first and then Arabinose (Left route in Figure 4A) were provided in the revision (Figure 4—figure supplement 1E), from which the low-low state cells has a percentage from 79.1% (0 h after inducer II addition) to 21.2% (12 h), to 24.6% (24 h), to 30.6% (36 h). If the multistability were generated by the plasmid lost, then it would not be stable for the four populations till 36 h. Thus, we think the plasmid lost may increase the portion of low-low state cells for long-time cultures, but would not generate the four populations.

3) The sequential induction method is quite interesting and important for revealing multistability. However, its application to the T15 network (or other networks) could be explained better. I guess the sequence in which inducers are added determines how the network's dynamics changes from monostable to multistable – but the details may be important. If bistable, tristable, quadristable regions are small in the 4-dimensional space then the goal is to "hit" these regions as inducers are sequentially added. While this is shown for the toggle switch, it is unclear how was this done for the more complex networks. This should be much more clearly described, especially for T15 and possibly for some simpler networks. Overall, going from the toggle switch to T15 and even to simpler networks is a big step, which requires more explanation and investigation. In addition, to toggle switch would be worth adding to the experimentally measured networks as control, considering that its behavior is computationally predictable. Do the experiments validate the predictions in Figure 2A for the toggle switch?

We thank the reviewer for this great comment. We think sequential induction is a good strategy to quickly evaluate multistable potentials and screen out systems harboring the most interesting dynamics. We added more detail about the sequential induction strategy and the application to T15 and other networks in the revised manuscript (subsection “Systematical multistability evaluation of MINPA and its sub-networks”, fourth paragraph), and Appendix (“Sequential Induction”).

To validate the theoretical analysis for sequential induction in Figure 2A, we constructed a synthetic toggle switch circuit and tested with sequential inductions. We first employed IPTG to induce the circuit for 5 h, and then aTc was added. Time course results showed that cells stayed at low-GFP state till 24 hours (Figure 2—figure supplement 1B). However, cells induced with aTc first, and then IPTG mainly stayed at high-GFP state, another stable steady state under this condition. Simultaneous aTc and IPTG induction produced similar cell distributions. These results show that sequential induction can be used as a strategy to quickly explore a multistable potential landscape for complex non-equilibrium systems. Details can be found in the third paragraph of the subsection “Systematical multistability evaluation of MINPA and its sub-networks”, Figure 2—figure supplement 1, and subsection “Strains, Media, and Chemicals”

Reviewer #3:

[…] 1) Recently, the theoretical advances have been able to identify the landscape and flux as the driving force of the global dynamics of the biological circuits (Proc. Natl. Acad. Sci. USA, 105: 12271-12276. (2008).). The quantification of the Waddington landscape has been achieved and is directly related to the underlying gene regulatory network for cell fate decision of stem cell differentiation and development (Proc. Natl. Acad. Sci. USA. 108(20):8257-8262(2011).). These should be reflected in the revision.

We thank the reviewer for the great suggestions. We have added the discussion and corresponding references in the revision (Introduction, second paragraph; and Discussion, first paragraph).

2) Furthermore, the multiple states have been predicted beyond the bistable states/switches without and with the epigenetic effects (Proc. Natl. Acad. Sci. USA. 108(20):8257-8262(2011); J. Phys. Chem. B, 115, 1254 (2011); Sci. Rep., 2,550 (2012); J R Soc Interface 10: 20130787 (2013); PLoS ONE. 9(8): e105216 (2014); Advances in Physics, 64:1, 1-137. (2015).). These should be reflected in the revision.

We thank the reviewer for the great suggestions. We have added the discussion and corresponding references in the revision (Introduction, second paragraph).

3) The authors should explain more clearly the physical and biological origins of the quadrastable states. Why the triple stable states are less likely than the quadrastable states? What is the origin of the LL state?

We thank the reviewer for the great question. The quadrastability presented in this work comes from coherent regulations in GFP and mCherry expressions. Experimentally, we used two hybrid promoters *Para/lac* and *Plux/tet* to drive the regulations. Each of the two promoters has two protein binding sites. AraC activates *Para/lac* in the presence of arabinose while LacI can bind and block its transcription. LuxR activates *Plux/tet* with induction of AHL while TetR can repress its transcription. Hence, each of the promoters has two possible states: ON and OFF. The interlocked circuit design of MINPA enables it has four possible expression levels in a cell population: *Para/lac* ON and *Plux/tet* OFF, *Para/lac* OFF and *Plux/tet* ON, *Para/lac* OFF and *Plux/tet* OFF, and *Para/lac* ON and *Plux/tet* ON. Since GFP and mCherry are reporters of the two promoters, so they could represent four stable states: GFP, mCherry, low-GFP and low-mCherry, and high-GFP and high-mCherry states.

Physically, the quadrastable states stem from the MINPA design, which is described as ordinary differentiation equations. Our mathematical analysis revealed that there are four stable steady states in the MINPA system. We didn’t compare the emergence probability between tristability and quadrastability. But experimentally, it is shown that tri-modality is easier to emerge, such as in the case of induction of MINPA with Arabinose and AHL (Figure 3C (C2_{LL}, C3_{LL})). However, quadra-modality can only occur with triple inductions (Arabinose, AHL and aTc, Figure 4B-D, and Figure 4−figure supplement 1E).

We also compared our model with previous studies describing the same topology (Proc. Natl. Acad. Sci. USA. 108(20):8257-8262(2011); PLoS ONE. 9(8): e105216 (2014)), we found the main difference is the description of the relationship between positive autoregulation and negative feedback. Previous studies described it as a “plus” (which means the two events are independent and additive), however, we described it as a “multiply” (which means the two events are dependent and multiplicative). Specifically, *Para/lac* activation depends on the relative concentrations of AraC (from node X) and LacI (from node Y) in the cell. LacI binding to the promoter would block RNA polymerase progression and hence repress transcription. Therefore, the LL state origins from the high effective concentrations of LacI and TetR in the cell. We believe incorporating mutual dependence of two binding sites of the hybrid promoter is more appropriate for the MINPA system described in this work, although the “additive” description of hybrid promoter could be suitable for other systems with different regulatory mechanisms.

Articles from eLife are provided here courtesy of **eLife Sciences Publications, Ltd**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |