Controlling spatiotemporal pattern formation in a concentration gradient with a synthetic toggle switch

Abstract The formation of spatiotemporal patterns of gene expression is frequently guided by gradients of diffusible signaling molecules. The toggle switch subnetwork, composed of two cross‐repressing transcription factors, is a common component of gene regulatory networks in charge of patterning, converting the continuous information provided by the gradient into discrete abutting stripes of gene expression. We present a synthetic biology framework to understand and characterize the spatiotemporal patterning properties of the toggle switch. To this end, we built a synthetic toggle switch controllable by diffusible molecules in Escherichia coli. We analyzed the patterning capabilities of the circuit by combining quantitative measurements with a mathematical reconstruction of the underlying dynamical system. The toggle switch can produce robust patterns with sharp boundaries, governed by bistability and hysteresis. We further demonstrate how the hysteresis, position, timing, and precision of the boundary can be controlled, highlighting the dynamical flexibility of the circuit.


Introduction
Synthetic biology aims to engineer living organisms with standardized and modular circuits that perform their functions in a programmable and predictable way (Brophy and Voigt, 2014, Cameron et al., 2014, Purcell and Lu, 2014. In addition to the promise of providing new technologies for medical and industrial applications (Gilbert and Ellis, 2018, Kitada et al., 2018, Nielsen and Keasling, 2016, Xie and Fussenegger, 2018, recapitulating biological processes synthetically provides a route to understanding the basic necessary mechanisms underpinning biological functions and dissect their properties and limitations Collins, 2018, Li et al., 2018).
Formation of spatiotemporal patterns of gene expression, a crucial process during the development of multicellular organisms, lends itself to be studied by such a synthetic biology approach. During development, pattern formation is achieved through a set of inter-connected gene regulatory programs encoding different non-linear responses to spatial chemical cues. This multiscale complexity makes the elucidation of the core principles of spatial patterning very challenging in living embryos, calling for alternative approaches capable of interrogating and comparing different pattern formation mechanisms.
The rise of synthetic biology has successfully allowed to build synthetic systems able to explore core patterning principles (reviewed in (Santos-Moreno and Schaerli, 2019b, Luo et al., 2019, Davies, 2017, Ebrahimkhani and Ebisuya, 2019). In addition, synthetic pattern formation is also an attractive technology for the engineering of living materials (Gilbert and Ellis, 2018, Nguyen et al., 2018, Moser et al., 2019, Cao et al., 2017 and tissues (Davies and Cachat, 2016, Healy and Deans, 2019, Webster et al., 2016. One ubiquitous strategy of patterning during embryogenesis is positional information, in which signaling molecules -the morphogens -diffuse and generate concentration gradients to specify positions. Specific gene regulatory programs are able to translate the spatiotemporal information provided by the local concentration of morphogen gradients into robust gene expression patterns (Wolpert, 1996, Green and Sharpe, 2015, Rogers and Schier, 2011. This mechanism has been extensively used in synthetic systems to generate spatial patterns, especially stripe patterns, which were produced through synthetic feed-forward loops (Basu et al., 2005, Schaerli et al., 2014, inducible promoters (Grant et al., 2016) and AND-gates (Boehm et al., 2018).
One of the gene regulatory subnetworks able of interpreting positional information is the bistable genetic switch (Zhang et al., 2012, Srinivasan et al., 2014, Balaskas et al., 2012, Kraut and Levine, 1991, Lopes et al., 2008b, Perez-Carrasco et al., 2016, Sokolowski et al., 2012, Zagorski et al., 2017, Alon, 2007, known as toggle switch (TS) ( Figure 1A). The topology of this circuit consists of two cross-repressing nodes that result in the binary mutually exclusive stable expression of one of the nodes. If the expression is influenced by an external signal, the TS provides a mechanism to convert a concentration gradient of this signal into stripes of gene expression (Perez-Carrasco et al., 2016, Sokolowski et al., 2012. Examples of TS-controlled pattern formation have been identified in the mesoderm formation in Xenopus (Saka and Smith, 2007), Drosophila blastoderm gap gene segmentation (Verd et al., 2019, Clark, 2017, and neural specification in vertebrate neural tube Small, 2015, Perez-Carrasco et al., 2018).
The non-linearity of gene regulatory networks such as the TS impedes to intuitively understand the effect that different kinetic parameters have on the dynamics of gene expression. For this reason, during the last decade, gene regulatory networks have been analyzed using tools from dynamical systems theory, associating the stable steady states of the dynamical system with the attainable gene expression states of a cell. Similarly, the change in the availability of cellular states as a consequence of perturbations of kinetic parameters of the network can be associated with the bifurcations of the dynamical system, providing information of the constraints of the possible cellular states. Indeed, the dynamical system of the TS has been thoroughly analyzed in silico, both in single cells as well as in population-level patterning scenarios (Perez-Carrasco et al., 2016, Ferrell, 2002, Guantes and Poyatos, 2008, showing that two possible stable states can coexist for a region of parameters inside which the expression state of the cell will depend on the initial gene expression -a phenomenon known as bistability. Under the control of an external signal, this bistability leads to hysteresis in which the state of the system is robust to signal changes, thus providing memory to gene expression patterns (Wang et al., 2009). Interestingly, the analysis of the steady states of the dynamical system not only provides static information of the cellular states, but also information on the transient dynamics of gene expression by which the states are attained (Verd et al., 2014). Therefore, a map of the underlying dynamical system is critical to fully understand the dynamics of the TS network.
The first synthetic version of the TS network was built almost 20 years ago and was a milestone of synthetic biology (Gardner et al., 2000). Since then, it has been built multiple times, extensively studied and used for its memory, bistability or hysteresis properties (Padirac et al., 2012, Purcell and Lu, 2014, Chen and Arkin, 2012, Lou et al., 2010, Andrews et al., 2018, Yang et al., 2019, Nikolaev and Sontag, 2016, Zhao et al., 2015, Sokolowski et al., 2012, Lebar et al., 2014, for stochasticity fate choice (Axelrod et al., 2015, Perez-Carrasco et al., 2016, Sekine et al., 2011, Wu et al., 2013, Lugagne et al., 2017 and to tune threshold activation (Gao et al., 2018). Nevertheless, its patterning capabilities controlled with a morphogen-like signal have not been studied in a synthetic system. Here, we constructed a "morphogen"-inducible synthetic TS network and studied its ability to produce spatial patterns -governed by bistability and hysteresis -in an Escherichia coli (E. coli) population. A combination of experiments and mathematical modelling allowed us to characterize the underlying bifurcation diagram, unveiling the possible dynamical regimes of the circuit. This enabled us to demonstrate how the inducible TS allows to control hysteresis, precision, position and timing of the pattern boundary.

Inducible toggle switch topology and its spatial patterning behavior.
The inducible TS network consists of two mutually repressing nodes ( Figure 1A). This mutual inhibition architecture ensures that only one of the nodes can be maintained at high expression. The dichotomous response of the cell depends on the asymmetry of the repression strengths and on the production rates of each node. An array of cells under a concentration gradient in charge of controlling any of these two properties will generate a binary spatial pattern. We built a synthetic inducible TS circuit starting from a characterized TS (Litcofsky et al., 2012). The first node of the network is composed of the TetR repressor and the mCherry reporter and will be referred hereon as the red node. The second node contains the LacI repressor and GFP reporter and will be referred as green node ( Figure 1B). Accordingly, the two expression states that the circuit can maintain will be referred as green and red states. TetR and mCherry are regulated by the hybrid pLuxLac promoter (BBa_I751502), which is activated by the LuxR-AHL (N-(β-Ketocaproyl)-L-homoserine lactone) complex and repressed by LacI, whose repression strength can be regulated by isopropyl β-d-1-thiogalactopyranoside (IPTG). LacI and GFP are controlled by the pLtetO promoter, which is repressed by TetR. LuxR is constitutively expressed from a pLuxL promoter on a second plasmid ( Figure 1B). In the absence of AHL (inducer), TetR and mCherry are not expressed, but LacI and GFP are, resulting in the green state. In the presence of AHL, TetR and mCherry can be expressed and repress LacI and GFP expression. Consequently, the system switches to the red state, provided that the concentration of the regulator (IPTG) is high enough.
We studied the capability of the inducible TS circuit to pattern an E. coli population exposed to different concentrations of AHL and IPTG and combinations thereof ( Figure 1C-D). To this end, we grew the cells on a hydrophobic grid placed on an agar plate to give a defined spatial organization to the cells (Grant et al., 2016), while small molecules can freely diffuse between grid squares. AHL and IPTG were pipetted at the left and at the top edges of the grid, respectively, forming gradients of AHL and IPTG by diffusion. We then measured the bacterial fluorescence after overnight incubation. When the starting cells were in the green state (reached by previous incubation in the absence of AHL and IPTG), the switch to the red state was observed only in presence of both AHL and IPTG, in the top left corner of the grid ( Figure 1C). In contrast, the same experiment performed with cells initially in the red state (reached by previous incubation in the presence of AHL and IPTG), showed that the red state is maintained above a certain concentration of AHL, mostly independent of the concentration of IPTG, leading to a green domain at the right ( Figure 1D). An overlay of the two grid patterns allows to highlight the possible stable states available for different concentrations of AHL and IPTG showing two monostable regions for the red and green states, and a bistable region for low values of AHL and high values of IPTG that grant hysteresis to the system ( Figure 1E). In addition to the tuning of the LacI repression by IPTG as explored here, the TetR repressing strength can be regulated by the addition of

Characterization of the toggle switch as a dynamical system
To quantitatively analyze the network behavior with single cell resolution, we grew the TS bacteria in liquid culture for 10 h with defined concentrations of inducer (AHL) and regulator (IPTG) molecules and used flow cytometry to measure the green and red fluorescence of individual bacteria ( Figure 2A). We gated the bacteria to quantify the percentage of cells in the red or green state for different inducer and regulator concentrations and different initial states. The observed behavior is consistent with the spatial patterning on the grid: Starting with an initial population of bacteria in the green state, the entire population (>90% of events) switched to the red state at high concentrations of both AHL and IPTG, whereas they stayed green otherwise. On the other hand, starting in the red state, the entire population stayed red above AHL concentrations of ~11 nM, but switched to the green state at lower AHL concentrations.
At the boundary between the red and green areas (white squares in Figure 2A), we observed cells in both flow cytometry gates. At high IPTG concentrations (≥0.5mM), a single population was located at intermediate red and green fluorescence values, whereas at IPTG concentrations ≤0.5mM the population split into two subpopulations displaying either the red or the green state. A bimodal distribution at the boundary between the two stable states is a known phenomenon of bistable circuits (Axelrod et al., 2015, Wu et al., 2013, Sekine et al., 2011 and can be caused by cell state switching in response to intrinsic gene expression noise (Perez-Carrasco et al., 2016, Sekine et al., 2011. The transition from red to green is less affected by the concentrations of IPTG than the transition from green to red, because of the asymmetry of the network. In the transition from red to green, the ratelimiting step is the degradation and dilution of TetR in the cells, since the clearance of TetR is required before changes in AHL or IPTG concentrations can induce the switch. Therefore, for the rest of the manuscript, we will focus on initial green state populations, while the data for the initial red state is in the supplementary data ( Supplementary Figures 4 and 5).
In order to fully characterize the dynamical capabilities of the circuit, we described the network with a mathematical model composed of ordinary differential equations (ODE) capturing the concentration of all the chemical species over time (see Methods for details). To unveil the dynamical landscape compatible with the synthetic circuit, we parameterized the model by fitting it to our experimental data.
In order to do so, we made use of multitry Markov chain Monte Carlo (MCMC) inference (Laloy andVrugt, 2012, Shockley et al., 2018), with a likelihood function based on the experimentally measured levels of expression and number of steady states found for the different concentrations of AHL and IPTG tested. We obtained a narrow probability distribution for all 11 parameters (Supplementary Figure 2).
The predicted outputs from the parametrized equation were able to recapitulate the bistability and hysteresis observed in the experimental data ( Figure 2B and Supplementary Figure 3). This allowed us to reconstruct the multidimensional bifurcation diagram of the system ( Figure 2C) which we can use as a map to predict the effect that different dynamical gradients will have on the observed gene expression patterns.
Analysis of the bifurcations of the system shows a scenario compatible with other TS, in which a continuous variation of AHL or IPTG can change the number of stable states via saddle-node bifurcations. Starting from a bistable condition, a saddle-node bifurcation occurs through the collision of a stable and an unstable state of the system, leaving only one possible stable state left ( Figure 2C). For a cell in the state about to disappear, a small perturbation of IPTG or AHL induces a sharp switch to the opposite state, producing a sharp transition between cell states.
There are two different saddle-node lines in the phase plane (AHL, IPTG) of the bifurcation diagram that collide at a certain concentration of the inducers (around 10 -3 µM AHL and 10 mM IPTG). Called a cusp bifurcation, this point separates the regions of inducer in the phase plane in which the switch between states occur through a saddle-node bifurcation (bistability) or through the continuous change of expression of a monostable state. Thus, the availability of a cusp point allows to explore the properties of two different patterning strategies (bistability versus sigmoidal response) for the same circuit topology.

Controlling the hysteresis, position, and sharpness of the boundary
Next, we analyzed how the IPTG concentration influences the hysteresis of the circuit, characterized by the range of inducer concentrations for which the circuit shows bistability ( Figure 2D). Since IPTG is in control of the repression of the green node over the red node, it is a perfect candidate to control the hysteresis of the TS. Based on the bifurcation diagram, we expected that the lower the IPTG concentration (i.e. the stronger the repression on the red node), the bigger the range of bistability.
Testing the bistability of the circuit at different IPTG concentrations confirmed this hypothesis, showing a response without bistability to an AHL gradient at high values of IPTG (10 mM) and an increasing range of bistability as IPTG decreases. The highest amount of hysteresis is observed in the absence of IPTG. Here, cells are not able to switch from green to red, even in presence of high AHL. However, AHL is enough to preserve the red state once reached. We thus observed three different dynamic regimes with distinct ranges of hysteresis: 1) no hysteresis (>1 mM IPTG) in a sigmoidal regime. 2) hysteresis, in a bistable regime with the possibility to switch between both states (around 0.125 mM IPTG) and 3) hysteresis in a bistable regime with irreversibility of the green state (around 0 mM IPTG). Those three regimes were also observed as spatial patterns on the grid, when placed on agar plates with uniform IPTG concentrations ( Figure 2D, right).
In addition to controlling the transition between a sigmoidal and a bistable regime and consequently the hysteresis, the IPTG concentration also affects the position of the boundary between the red and green states. Higher concentrations of IPTG allow the system to switch from green to red at lower AHL concentrations (Figure 2A and D), thus controlling the boundary position for a given AHL gradient.
Finally, IPTG also tunes the transition sharpness between the red and green states ( Figure 2E and Supplementary Figure 4). Starting from the green state, at high concentrations of IPTG beyond the cusp point, we observed a sigmoidal expression response in an AHL concentration gradient, similar to the one expected from a circuit with a single repressor (TetR) active. At lower IPTG concentrations, below the cusp point, the system displayed a sharper transition as expected from the saddle-node bifurcation.
Moreover, these different transitions are consistent with the distinct population behaviors during cell state transitions along the gradient: a smooth transition of a single population for the sigmoidal behavior and a bimodal population transition (Figure 2A) for the bistable switch behavior. In summary, the regulator IPTG allows us to choose between a bistable and a sigmoidal regime (at high IPTG concentrations) and thus to control the hysteresis, the position and the sharpness of the boundary.  Figure 3A). Similarly, for a constant amount of IPTG, the switching time depends on the AHL concentration, switching at times as slow as 6-10 hours for an IPTG concentration corresponding to the bistable region (0.125 mM Figure 3B). On the other hand, inducing a pattern in the sigmoidal regime (10 mM IPTG), we observe a more consistent switching time (4 h) across different AHL concentrations ( Figure 3C).
Interestingly, the behaviors of the switch at the population level over time are identical to the ones observed across different AHL concentrations: at high IPTG (beyond the cusp bifurcation, in the sigmoidal regime) we measured one population moving as a whole, while at lower IPTG concentrations (bistable regime) we observed the cells splitting into two divergent subpopulations ( Figure 3A,B and Supplementary Movie 1), suggesting that the dynamics of the sigmoidal patterning are a result of the population relaxing to the unique possible steady state, whereas in the bistable regime the transition is controlled through the stochastic switching between the two possible states, with a time that is determined by the intrinsic noise and can therefore be slower than the degradation rate of the proteins composing the TS.

Patterning in the sigmoidal regime is faster than in the bistable regime
Combining our model with 2D diffusion allowed us to reproduce the patterns observed in the grid assay shown in Figure 1 (Figure 4A). In addition, integrating the diffusion of the inducer with the dynamical properties of the bifurcations of the system can shed light on the different patterning dynamics observed.
In particular, consistent with our flow cytometry data (Figure 3), the model suggests that the switching slows down close to the saddle node bifurcation, a phenomenon called critical slow down. Thus, for different constant values of IPTG, a gradient of AHL is expected to create a moving boundary at different speeds for cells switching from the green to the red state ( Figure 4B) (Perez-Carrasco et al., 2016). To test this prediction, we measured the position of the boundary over time in the grid assay at two different IPTG concentrations corresponding to two different dynamical regimes (sigmoidal, 1 mM and bistable, 0.125 mM) ( Figure 4B). As expected, we observed that the transition to the production of mCherry starts earlier and advances faster in the sigmoidal regime, equipping the TS with time control of the pattern formation through the cusp bifurcation. This was confirmed through simulation of the diffusion of the inducers on the grid, were the model predicts the same trend than the experimental data, with no more than three squares of difference. On the other hand, we did not observe clear differences in the precision of the boundary in the grid assay, probably due to a lack of precision in the concentration gradients accessible by the grid.

Spatiotemporal control of the pattern by using spatially homogenous signals
In addition to the pattern formation through a bistable or monostable sigmoidal regime, the TS offers the possibility to move between both different regimes during the pattern formation, allowing for different patterning strategies that can exploit the properties of both regimes. In particular, by manipulating the homogeneous levels of IPTG in time for a given gradient of AHL, we can control the pattern formation process. An initial population of cells in the green state in the absence of IPTG (bistable irreversible regime) is unaffected by the gradient of AHL. Adding IPTG homogeneously at the desired time point, brings the cells to the sigmoidal monostable regime able to respond to the gradient of AHL and forming a boundary. Once the boundary is located at the desired position, removing IPTG takes the system back to the bistable zone, thus freezing the boundary and making it robust to changes in AHL ( Figure 5A). Consequently, the system is able to maintain a pattern in the absence of IPTG which removes the requirement of maintaining a precise AHL gradient to keep the pattern boundary. Therefore, a pulse of IPTG can combine advantages of two distinct regimes: of the sigmoidal monostable regime for a fast establishment of a pattern ( Figure 4) and of the irreversible regime to make the pattern robust to changes in the AHL gradient. To demonstrate this effect, we grew cells initially in the green state at different concentrations of AHL and exposed them to pulses (2 h, 3.5 h ,5 h, 6.5 h) of 10 mM IPTG ( Figure 5B).
We then grew them for further 6 h in the absence of IPTG, but keeping the AHL concentrations used during the pulse. Indeed, cells receiving enough AHL (>= 0.01 µM), were able to maintain the red state.
Interestingly, this was the case even for short pulses were not the whole population did have enough We used the grid assay to further demonstrate this memory property of the TS and to show that the pattern in the bistable regime is indeed robust to changes in the AHL gradient. We patterned a grid as in Figure 1 and transferred the cells onto a new agar plate where the positional information provided by AHL and IPTG was removed (0 mM IPTG, homogeneous concentration of 5 µM AHL). As predicted, the pattern was maintained in the absence of any spatial information (Supplementary Figure 6). This result demonstrates that an inducible TS network is capable of interpreting spatiotemporal gradients by making use of the memory properties of the circuit.

Discussion
Systems displaying hysteresis, and in particular the bistable switch, allow for a sharp threshold response that can turn a graded input into a binary output (Perez-Carrasco et al., 2016, Sokolowski et al., 2012. Interestingly, the mutual repressing motif of the toggle switch is widely found in natural pattern-forming systems, for example, in networks responsible for patterning the neural tube (Balaskas et al., 2012), the dorsal telencephalon (Srinivasan et al., 2014), the Drosophila embryo segments (Verd et al., 2019, Clark, 2017 and the Xenopus mesoderm (Saka and Smith, 2007). However, in addition to forming sharp boundaries, the successful formation of a pattern requires the control of the position and timing of the boundary formation. In order to explore dynamic patterning properties of the TS gene regulatory network, we have successfully built and characterized a synthetic inducible toggle switch with tunable repression (Figure 1).
We quantified the gene expression at the single cell level over a wide range of inducer and regulator concentrations and fit the data to parameterize a mathematical model of gene regulation. The resulting bifurcation diagram of the dynamical system provided the mechanistic understanding required to interpret the different spatiotemporal patterns observed and experimental design guidance for pattern formation with the synthetic TS. In particular, this approach allowed us to characterize the mechanism by which the system can transit from a bistable regime to a sigmoidal unimodal response to the inducer, determined by the cusp bifurcation of the dynamical system (Figure 2). From the exploration of the differences in boundary precision of both regimes, we observed that while the bistable regime provides the sharper boundary (as previously reported (Lopes et al., 2008a, Isalan et al., 2005), the sigmoidal regime allows for a faster response (Figures 3 and 4). This reveals a trade-off between the timing and precision of the boundary. These observations are in consonance with dynamical systems predictions, in which the sharp transition of the bistable regime -through a saddle-node bifurcation -comes at the price of the critical slowing down close to the bifurcation (Narula et al., 2013, Perez-Carrasco et al., 2016).
In addition, for different levels of a spatially homogeneous regulator, we have been able to control the range of hysteresis and the position of the pattern boundary for a given inducer gradient changed accordingly ( Figure 2). This mechanism is analogue to the one proposed by Cohen et al. (Cohen et al., 2014) to control the boundary position of patterns regulated by morphogen gradients by changing the affinity of one of the elements to a spatially homogeneous transcription factor. This result highlights the evolutionary potential of the circuit, allowing for kinetic mechanisms of controlling the position of the boundary without compromising other aspects of the boundary or requiring an upstream regulation of the morphogen gradient.  (Balaskas et al., 2012). In this latter scenario, the adapatation to the gradient of sonic hedgehog (through Gli signalling) results in a combination of temporal profiles of activation and repression acting on the patterning circuit (Junker et al., 2014, Cohen et al., 2015. This highlights the importance to understand bistable switches inside a dynamical scenario. While such dynamics is challenging to measure in the developing embryo, our synthetic biology framework allowed for an alternative route towards understanding the dynamics of the TS in pattern formation. Our improved understanding of the TS for pattern formation opens up new avenues of research. In particular, for many bistable parameter conditions, single cell expression showed the coexistence of two subpopulation of cells at each one of the available stable states. This provides evidence of the relevance of intrinsic noise in the establishment of the pattern. While previous in silico research shows that intrinsic noise can determine the position and precision of morphogen driven boundaries (Perez-Carrasco et al., 2016, Weber andBuceta, 2013), the actual role of noise in the dynamics of pattern formation in living systems and the possibility of optimal dynamical strategies based on the stochasticity of gene expression remain still a conundrum. In addition, it poses new challenges to dynamical system inference in which different sources of intrinsic noise must be disentangled from measurement noise in order to obtain an accurate characterization of the circuit (Dharmarajan et al., 2019).
Overall, our results underscore the relevance of studying the dynamical context of a gene regulatory network in order to understand patterning processes, not only of synthetic circuits but also of developmental systems (Ferrell, 2002, Sagner and. Future work will reveal if the TS properties are conserved when incorporated in more complex (synthetic) gene regulatory networks, for example when combined with the repressilator (Elowitz and Leibler, 2000, Potvin-Trottier et al., 2016, Santos-Moreno and Schaerli, 2019a to yield the AC-DC network (Perez-Carrasco et al., 2018, Verd et al., 2019, Balaskas et al., 2012, Panovska-Griffiths et al., 2013. Moreover, the here established engineering guidelines on how to control patterning with a synthetic TS will be valuable for future synthetic pattern formation, for example in the context of engineered living materials based on bacterial biofilms (Gilbert and Ellis, 2018, Nguyen et al., 2018, Moser et al., 2019, Cao et al., 2017.

Cloning of the inducible toggle switch circuit
We cloned the morphogen inducible toggle switch plasmid ("TS_pLuxLac") from an already functional toggle switch plasmid (pKDL071) (Litcofsky et al., 2012), kindly provided by Jeong Wook Lee. The two Ptrc-2 promoters upstream of TetR and mCherry coding sequences were replaced by the hybrid promoter pLuxLac (BBa_I751502) with the plug-and-play method initially used to assemble pKDL071 using the restrictions sites NcoI and SalI upstream of TetR and XmaI and MfeI upstream of mCherry.
For the "pCDF_luxR" plasmid, the pLuxL promoter (BBa_R0063) and the LuxR gene (synthesized by GenScript) were introduced into a pCDF plasmid with a customized multiple cloning site (Schaerli et al., 2014) between the KpnI and BamHI sites. The sequences and plasmids will be available through Addgene.

Strain and growth condition
pCDF_luxR and TS_pLuxLac were transformed into the E. coli strain MK01 (Kogenaru and Tans, 2014).
The absence of the lacI gene in this strain avoids unexpected cross talk between the synthetic circuit and the host.
Bacteria were turned into red state by inoculating single colonies into 4 ml M9 medium in presence of 1 mM IPTG (isopropyl β-D-1-thiogalactopyranoside) and 10 µM AHL (N-(β-Ketocaproyl)-L-homoserine lactone). They were grown overnight at 37° C and 200 RPM shaking. The same procedure was used to turn the bacteria into the green state, with the difference that the medium did not contain IPTG and AHL.
These bacteria were plated out on LB agar plates (supplemented with 10 µM AHL for the red state) and incubated overnight at 37° C. The plates were stored at 4° C and single colony from them were picked for the experiments.

Flow cytometry
Single colonies in the red or green state were cultured in M9 for 4-6 h and put at 4° C before entering in stationary phase (below 0.8 OD). The following day, these cultures were diluted to 0.01 OD (NanoDrop 2000, Thermofisher) and added into the wells of a 96-well plate (CytoOne) to a total volume of 120 µL including indicated concentrations of AHL and IPTG. The plate was incubated in a BioTek Synergy H1 plate reader at 37° C with 548 cpm (2 mm) double orbital shaking speed. Absorbance (600 nm) was monitored every 10 min to check that cells did not enter stationary phase (below 0.3 in the plate reader) as cells in stationary phase can no longer switch between the two states. After 5 h, cells were diluted 100 times before incubating them again under the same conditions for a total of 10 h. For the IPTG time pulse experiment, we started the experiment with a culture of OD 0.1 in a medium containing 10 mM IPTG and the indicated concentration of AHL. The samples were incubated as described above until the first time point (2 h) and an aliquot was stored at 4° C. Then, the cells were diluted 1:5 into fresh medium containing the same IPTG and AHL concentrations and further incubated until the next time point (3.5 h). This procedure was repeated for all time points. An aliquot of each sample was analyzed by flow cytometry once all samples were collected. The next day, the cells where diluted 1'000 times into fresh medium with the indicated concentration of AHL and no IPTG. The cells were grown for 6 h and directly measured by flow cytometry.

Grid assay
Single colonies in the red or green state were cultured in M9 for 4-6 h and put at 4° C before entering in stationary phase (below 0.8 OD). The following day, these cultures were washed from IPTG and AHL by centrifugation for 1 min at 13'000 RPM and resuspended in fresh medium without IPTG or AHL. Then cells were diluted to an OD of 0.1 before being pipetted onto the grid (ISO-GRID from NEOGEN) (Grant et al., 2016) (20 μl of cells were loaded for 10 lines, approximatively 0.05 μl per square) that was placed on a M9 agar plate with 0.2% (w/v) glucose and appropriate antibiotics (25 µg/mL kanamycin and 25 µg/mL spectinomycin). The inducers (100 mM IPTG and 100 μM AHL, 5 μl each) were added as indicated in the figures. The plate was incubated overnight at 37° C before green and red fluorescence were measured with a Fusion FX (VILBER) imaging system. We used 0.2 ms exposure with blue light (480nm) and a F-565 filter for the GFP measurement and 1 min exposure with red light (530nm) and F-740 filter for mCherry. ImageJ software (R Core Team, 2017) was used to analyze the picture.
For the time course experiment, bacterial solution with an OD of 3 were loaded on the grid in order to be able to quantify early timepoints. Fluorescence intensity values for each square of the grids were extracted with a custom Fiji imageJ (Schindelin et al., 2012) macro script and the maximum value for each square was normalized to the highest measured fluorescence in all the conditions and replicates.
The data was plotted with R software (R Core Team, 2017).

Model derivation and parametrization
We model the expression dynamics of the TS by describing the change in time of the concentration of the [LacI] and [TetR] as a balance between their regulated production and degradation. In addition, we made the assumptions that the dynamics of the transcripts and promoter binding/unbinding is faster than the dynamics production and degradation of the repressor proteins, and that the reporters follow the same dynamics as their associated repressors.
The expression of LacI is regulated by TetR, which can inhibit the production of LacI by binding to the TetO promoter. Modelling this interaction as a repressive Hill function we can describe the evolution in time of LacI as, 1 Where is the maximum production of LacI in the absence of the repressor TetR, KTetR the TetR concentration producing half repression, nTetR the Hill coefficient and its degradation rate.
On the other hand, the expression of TetR is regulated at the same time by the repression of free LacI 1 In order to parametrize the mathematical model Eqs. (1-3), we compared the experimentally observed cellular states with the stable steady states of the theoretical dynamical system for different sets of concentrations of AHL and IPTG. The experimental steady states were obtained by using the gated expression in cellular populations (see Flow Cytometry for details). For each gated population, the state was accepted when it contained at least 15% of the cellular population.
The stable states for the mathematical model ( In order to compare the computational steady states with fluorescence measurements we assumed a linear relationship between fluorescence and gene expression, * , * Thus, the parametrization of the problem is reduced to the inference of 11 identifiable parameters , , , , , , where the superindex presents the different theoretically predicted stable states in the case of bistability. Finally, parameter penalizes that for the condition given by state (concentrations of AHL and IPTG) the number of different states detected experimentally is not the same than the number of stable states predicted in the mathematical model. Thus, promotes parameter sets that match monostable and bistable regions between experimental and computational steady states. In particular 0 if the number of stable steady states match, and 0 otherwise. For the inference used in the manuscript, a value of 10 was used. The inferred parameters are given in Table 1.  Figure 3. For each parameter, the median and the 95% credibility interval of each marginal distribution are indicated. The median of each parameter corresponds with the value used in the mathematical model for the rest of simulations of the manuscript.

Morphogen diffusion
The diffusion of the morphogens (AHL, IPTG, and aTc) on the agar plate is assumed to be a 2-D free diffusion determining the concentration of each chemical , , at different positions of the plate. This is modelled through the partial differential equations , . This equation was solved by discretizing the space using the square experimental grid. The boundary conditions imposed by the experiment are given by the constant concentration of morphogen along one side of the square 0, , , and assuming sinks in the rest of the perimeter of the square , , , 0, , , 0 , where L is the length of the side of the square. In order to take into account the dilution effect of the initial concentration pipetted on the agar before stablishing the gradient, we set . Parameters used for the grid assays are 10 and 2.5 10 cm /h (Basu et al., 2005, Miyamoto et al., 2018. Finally, in order to provide a time scale to the dynamics of protein, not available from the MCMC fitting of the steady states of the system, the degradation rate of the LacI and TetR was set to