Interpretation of morphogen gradients by a synthetic bistable circuit

During development, cells gain positional information through the interpretation of dynamic morphogen gradients. A proposed mechanism for interpreting opposing morphogen gradients is mutual inhibition of downstream transcription factors, but isolating the role of this specific motif within a natural network remains a challenge. Here, we engineer a synthetic morphogen-induced mutual inhibition circuit in E. coli populations and show that mutual inhibition alone is sufficient to produce stable domains of gene expression in response to dynamic morphogen gradients, provided the spatial average of the morphogens falls within the region of bistability at the single cell level. When we add sender devices, the resulting patterning circuit produces theoretically predicted self-organised gene expression domains in response to a single gradient. We develop computational models of our synthetic circuits parameterised to timecourse fluorescence data, providing both a theoretical and experimental framework for engineering morphogen-induced spatial patterning in cell populations.


Supplementary Figures
Supplementary Figure 1: Circuit variants. A circuit diagram of the exclusive receiver circuit with genetic parts labelled in black that were included in the final circuit and in red that were evaluated but discarded. Not all combinations of parts were tested. nM C12 then exposed to the indicated concentrations of C6 and C12. YFP fluorescence is plotted on the X-axis while CFP fluorescence is plotted on the Y-axis. b, Cells conditioned in 500 nM C6 then exposed to the indicated concentrations of C6 and C12. Square indicates the position of untreated cells. c, Gating strategy. Cells also constitutively express mRFP1 via a genomic transgene. Only RFP + cells were used for analysis and electronic noise, cell debris and cell clusters were excluded sequentially. Source data are provided as a Source Data file. Figure 6: Microfluidics measurements of bistability. Single-cell data used to compute ratios in Figure 2c. Cells were grown in microfluidic chips for 3 hours in the presence of either 37 nM C6 (rows 1 and 2) or 100 nM C12 (rows 3 and 4). Then media was changed to 100 nM C12 + 37 nM C6 (rows 1 and 3) or 100 nM C12 (row 2) or 37 nM C6 (row 4) . Cells were imaged with a frame rate of (1 frame/10 minutes). Left panels in each column are kymographs of the CFP (left column) or YFP (right column) expression per-cell, and fraction of cells as a heat map. Histograms represent the populations at 3 hours (red) and 8 hours (blue). Lines and shaded region represent the mean and standard deviation, respectively, over n = 4 biological replicates performed over 4 different days. Figure 7: Bistability and switching of single cells is robust to high C6 signal concentration. Cells were grown in microfluidic chips for 3 hours in the presence of either 1 µM C6 (rows 1 and 2) or 100 nM C12 (rows 3 and 4). Then media was changed to 100 nM C12 + 1 µM C6 (rows 1 and 3) or 100 nM C12 (row 2) or 1 µM C6 (row 4). Cells were imaged with a frame rate of (1 frame/10 minutes). Left panels are kymographs of the log-ratio of CFP expression per-cell to YFP expression per-cell, and fraction of cells as a heat map. Histograms represent the populations at 3 hours (red) and 8 hours (blue). Lines and shaded region represent the mean and standard deviation, respectively, n = 4 biological replicates performed over 4 different days. Right panels are sample montages of cells switching state (rows 3 and 4) or exhibiting bistablity (rows 1 and 2); phase contrast and fluorescence channel ranges chosen for display.

Supplementary Figure 8: Microfluidics measurements of bistability with high C6 signal concentration.
Single-cell data used to compute ratios in 7. Cells were grown in microfluidic chips for 3 hours in the presence of either 1 µM C6 (rows 1 and 2) or 100 nM C12 (rows 3 and 4). Then media was changed to 100 nM C12 + 1 µM C6 (rows 1 and 3) or 100 nM C12 (row 2) or 1 µM C6 (row 4) . Cells were imaged with a frame rate of (1 frame/10 minutes). Left panels in each column are kymographs of the CFP (left column) or YFP (right column) expression per-cell, and fraction of cells as a heat map. Histograms represent the populations at 3 hours (red) and 8 hours (blue). Lines and shaded region represent the mean and standard deviation, respectively, over n = 4 biological replicates performed over 4 different days.
Supplementary Figure 9: Boundaries summarized in Figure 3c Endpoint fluorescence microscopy of Exclusive Receiver cells grown in transient gradients of signals (C12 diffusing from the left, C6 diffusing from the right) at the spatial average concentrations indicated and in the context of 10 µM IPTG throughout. Figure 10: Swapping primary and secondary morphogens also produces patterning. a, A circuit diagram of exclusive reporter cells co-transformed with the P81-LuxI relay device that responds to C12 by producing C6. Cells that experience high levels of C12 (central cells) will express YFP, TetR, and LuxI, causing them to produce C6 but be unable to sense it. Neighbouring cells (outer cells) that do not experience C12 will sense C6 and express CFP and LacI. b, Isogenic cells transformed with the circuit shown in (a) grown for 24 hours in the presence of a gradient of C12 diffusing from the centre express CFP and YFP in mutually exclusive domains of gene expression.c, A simulation in C6-C12 space over time (t1-t2) labelling points in physical space by their CFP and YFP expression (cyan and yellow points), and showing the production of C6 as vectors (red arrows) that move the spatial average (black point) toward increasing C6. The bistable region is outlined in red.

Differential Equation Models and Parameter Inference
In this section, we derive ordinary differential equation (ODE) models for the reaction kinetics underlying the Exclusive Receiver circuit. These derivations broadly follow the derivations of the Receiver circuit in [3] and [4]. Importantly, we introduce differences in that original derivation that lead to changes in the location of bifurcations in (C12,C6) space, when those derivations are extended to incorporate dynamics of the repressor proteins LacI and TetR, and their chemical inhibitors IPTG and ATC.

Dynamic characterization with inference graphs
In order to infer the parameters of the Exclusive Receiver circuit, we adopt the strategy described in [4], evaluating parameters of sub-circuits first, and propagating their inferred values to larger circuits that embed those same parameters. This results in an inference graph, where we infer parameters over a sequence of models and corresponding datasets.  Blue rectangular nodes are inference problems for collections of synthetic gene circuit models, compared with CFP and YFP measurements, while green rectangular nodes are inference problems for the growth models of the circuits in the downstream blue nodes, compared with OD 600 measurements. Internal to the coloured nodes are white elliptical nodes, which correspond to individual synthetic gene circuits and one or more associated dataset(s). White rectangular nodes are the sets of inferred parameters that are propagated between nodes. To simplify the notation, φ has been used to denote the set of growth model parameters, which are culture-specific (local) values for r, K and t lag , and are propagated as maximum likelihood estimates. All other parameters are global (non-culture-specific) values and propagated as marginal posterior estimates.
The simplest possible cell line to characterize is one in which there is no synthetic gene circuit at all. Applying dynamic characterization in this context enables us to quantify autofluorescence, and so we name this circuit auto. Therefore, we measured cells under a range of conditions to explore how gene expression capacity influenced time-series measurements at fluorescence wavelengths corresponding to CFP and YFP. Subsequently, we characterize a circuit (prpr) in which CFP and YFP are driven by constitutive promoters (PR), enabling us to characterize the rates of degradation of the fluorescent proteins. Next, we use four variants of simple HSL receiver (PCat, R100S32, R33S32 and R33S175) to characterize the genetic parts associated with LuxR and LasR receiver proteins and their interactions with HSL molecules 3OC6HSL (C6) and 3OC12HSL (C12). Finally, having obtained parameter estimates for parts associated with CFP, YFP, LuxR and LasR, we characterize the Exclusive Receiver circuit (exrep) itself, establishing quantitative estimates for the parts associated with TetR and LacI repressor proteins.
In the following, we introduce the models for each of the circuits just mentioned, define their parameters and present results of the inference.

Autofluorescence model
The model we used for autofluorescence assumes that the rate of autofluorescence is constant, and that the fluorescent material dilutes with cell growth. As such, the equations for intracellular autofluorescence corresponding to CFP and YFP are is the specific growth rate of the cell culture with density ρ. To compare with experimental measurements, we consider the bulk fluorescence given by where c 480 and c 530 are modelled as in (1). Here, the quantities b 480 and b 530 represent background fluorescence at 480 nm and 530 nm, corresponding to CFP and YFP respectively.

Inference.
We use the data in [4] to infer the parameters of the auto circuit. Cells were treated with EtOH to perturb cell growth, enabling us to determine how autofluorescence changes with different cellular growth rates. The priors used are detailed in the following table. Supplementary The marginal posterior estimates are shown in Figure 15. Simulation of the maximum likelihood estimate is shown in Figure 16.

Constitutive (prpr) model
The prpr circuit described above uses the constitutive PR promoter to drive CFP and YFP expression, in two separate operons. Following the derivation in [4], we arrive at a system of equations that describe the time-evolution of the intracellular concentrations of CFP and YFP as where a CFP and a YFP are aggregated parameters that incorporate the rate of transcription and translation of CFP and YFP respectively. To compare with bulk culture fluorescence data, we use the observer model where the dynamics of c 480 and c 530 are governed by equations (1) above.

Inference.
We use the data in [4] to infer the parameters of the prpr circuit. Cells were treated with chloramphenicol to perturb cell growth, enabling us to establish how constitutively expressed proteins are altered with different cellular growth rates. While the model above does not explicitly describe any explicit functional response to chloramphenicol, our general strategy of allowing the cell growth parameters to vary across different measurements enables the effect of chloramphenicol on growth rates to be implicitly captured. The quantification of autofluorescence was reused from the auto circuit (upstream in the inference graph), but the background fluorescence parameters were re-inferred. The priors used are detailed in Table 3. The marginal posterior estimates are shown in Figure 17. Parameter Description Unit Distribution Scaling

Supplementary
Simulation of the maximum likelihood estimate is shown in Figure 18. The simulated CFP and YFP largely agree with the measured fluorescence at the culture level. As the effect of chloramphenicol is not explicitly modelled here, this comparison indicates that the majority of the effect of chloramphenicol can be described via its effect on cell growth. Any additional direct effect on CFP and YFP expression directly, is relatively minor.

Receiver model
We consider the dynamic characterization of the HSL Receiver circuit introduced in [3] and modelled dynamically in [4,5]. In this circuit, which we refer to in the main text as the Receiver, two variations of the wild-type PLux promoter, PLux76 and PLas81, were engineered to bind preferentially to activated LuxR and LasR complexes respectively. As LuxR favours binding of C6 and LasR favours binding of C12, optimized expression of LuxR and LasR can lead to near-orthogonal intracellular detection of C6 and C12. The Receiver device was originally measured with PLux76 upstream of the coding sequence for CFP, and PLas81 upstream of the coding sequence for YFP.
Version 1 -Uniform degradation. The first version of the model we introduce is derived in [4]. It is based on the assumption that all degradation processes are of the same order, e.g. LuxR monomers are degraded at a similar rate as LuxR Dimers bound to HSLs. Subsequently, we will introduce an alternative derivation where we assume that complexes are protected from degradation, i.e. that degradation mainly occurs on the monomer level. We start by repeating some of the derivation from [4].
We denote by C k the HSL molecule with length k carbon chain, and by G i the PLux76 and PLas81 promoters. Then similar to the derivation in [3], we can specify all of the reactions between the HSLs, LuxR and LasR, and eventual binding of transcriptional regulators to PLux76/PLas81.
Constitutive expression of LuxR and LasR is described by Degradation of LuxR and LasR is described by where d 1 R and d 2 R are distinguishable from d R to describe the effect of HSL molecules protecting receiver proteins from degradation.
Inducible expression of CFP and YFP by P OLux and P OLas respectively is described by To produce a simplified ODE model amenable to parameter inference, we start with the equations describing LuxR and LasR protein, their complexes involving C 6 and C 12 , and the bound-/unbound promoters. Crucially, in this first derivation, we make a rapid equilibrium assumption for the binding reactions (5), and obtain the following relationships . Therefore (also symmetry of LuxR and LasR), where the new K's are defined as above. To reduce the system to fewer variables, we consider the evolution of total LuxR, and seek to co-ordinate this with the rapid equilibrium relationships above. By denoting the total concentration of LuxR as c R , we can write down its time-evolution as If there are differences between d R , d 1 R and d 2 R , then further reduction is complicated. However, by assuming that HSL is not protective of receiver protein (d 1 , we obtain the simplification Now using a conservation relationship for LuxR, we can obtain When C k is low, total LuxR is closely approximated by free LuxR, c R ≈ [LuxR]. But when C k is high, c R should be partitioned between the [LuxR-C k -LuxR-C k ] and [G i -LuxR-C k -LuxR-C k ] species. Therefore, to simplify the analysis, we propose modelling this by using the assumption which still captures the saturation of LuxR by C k , using the approximations where K (i) GR = K GiDk K Dk is assumed to be independent of which signal is bound (k), equivalent to the derivation in [3]. By symmetry, we immediately obtain equivalent expressions for interactions between LasR, HSL and PLux promoters. We denote the total concentration of LasR as c S .
In addition to the saturation of LuxR and LasR, our reduced model also allows for saturation of G i . By taking advantage of the conservation law where the fractions of bound LuxR and LasR are defined by Here, we have introduce alternative exponents n R and n S , analogous to the usage of n in [3]. Accordingly, we obtain the following system of equations

Inference for version 1 (uniform degradation).
To characterize the LuxR and LasR signalling components, we used measurements of the response of four Receiver circuits from [3] to treatment with C6 and C12 over 3-fold dilutions. The maximum LuxR and LasR production rates were normalized to the values corresponding to the Pcat promoters, as done previously [3], thus setting the scale for unobserved concentrations of LuxR and LasR. We used (uninformative) uniform priors on the previously uncharacterized parameters, and (informative) truncated Gaussian priors on f 480 , f 530 , d CFP and d YFP with mean and standard deviation taken from the marginal posteriors of the prpr circuit characterization. We did not propagate the marginal posteriors of a CFP and a YFP as the promoter involved differed between the prpr circuit and Receiver circuits. The priors used are detailed in the following table.
The marginal posterior estimates are shown in Figure 19. Simulation of the maximum likelihood estimate is shown in Figure 20.

Version 2 -Protected degradation.
Here we provide an alternative derivation for the model reduction based on the assumption that HSL complexes and dimers are protected from degradation, i.e. that degradation predominantly occurs on the level of LuxR/LasR monomers. This is supported by in vitro analysis of purified LuxR suggesting that LuxR protein is unstable in the absence of 3OC6-HSL [6]. Furthermore, we assume that all dilution effects occur on a slow time scale as compared to the kinetic. With this derivation, we obtain the same functional structure of the promoter activities of PLux76 and PLas81 from equation (15), but with the definitions of the bound fraction of LuxR and LasR given instead by where e R12 = K R12 K R6 and e S6 = K S6 K S12 result from dividing by K R6 and K S12 . Consequently, the parameters K  (15) are rescaled by K R6 and K S12 respectively.

Inference for version 2 (Protected degradation).
We carried out parameter inference for version 2 of the model directly equivalent to that done for version 1. The complete list of prior distributions for the parameters is as follows.

Exclusive Receiver model
To model the Exclusive Receiver, we consider the inhibition of PTet by TetR and PLac by LacI, and the mechanism of chemical inhibition by IPTG and ATC. For promoter regulation, we use the inhibition Hill function Typically, these functions would include a parameter for the half-saturation concentration, but we omit that here because, without loss of generality, [LacI] and [TetR] can be arbitrarily scaled by those half-saturation concentrations. In such a rescaling, the half-saturation constants get embedded within the maximal production rates, a L and a T . For the chemical inhibitors, we assume a reaction of the form inhibitor + repressor → sequestered product (20) Correspondingly, the action of IPTG and ATC is proportional to the product of its concentration and its target repressor protein concentration.
where P 76 and P 81 are defined in (15).

Inference for version 1 (uniform degradation).
The inference procedure was less robust for the Exclusive Receiver, as compared with upstream circuits in the inference graph. In particular, we found it was not possible to identify a unique mode within the parameter space when all parameters were allowed to be flexible. Specifically, there was strong interdependency within the triplet {a L , d L , i I } and the triplet {a T , d T , i A }. Our interpretation is that we are unable to fully recover the time-scales of variations in c L and c T , as they are likely to not vary much during the experiments we used for characterization. It's likely that c L and c T quickly stabilise to equilibria when the cells are transferred to the media containing the treatments (specific concentrations of C6, C12, IPTG and ATC). In which case, only those equilibrium values will be identifiable, and not the production and degradation rates separately. Therefore, in the final version of the inference results presented here, we have fixed the degradation rates of LacI and TetR to 1 h −1 . Even when running the inference with d L and d T fixed, we found that chain convergence was not perfect, and some chains got stuck in local optima. Therefore, in our marginal posterior estimates, we have only included chains that converged to relatively good likelihood scores Figure 23. The marginals clearly indicate some additional flexibility in the inferred parameter values, possibly resulting from the larger parameter space being navigated, which includes uninformative priors for the parameters listed in Table 6, but also some flexibility in all of the parameters associated with the Receiver module, despite them having a strong prior.

Inference for version 2 (Protected degradation).
The inference procedure was also not completely robust for version 2 of the Exclusive Receiver model. We used the same uninformative priors for the parameters specific to the Exclusive Receiver model as in version 1, including fixing d L and d T . Again, in our marginal posterior estimates, we have only included chains that converged to relatively good likelihood scores Figure 24.

Bistability Analysis
In this section we outline computations used to create Figure 2b in the main text which compares regions of bistability indicated by hysterisis experiments in flow cytometry to that given by the differential equation model for the exclusive receiver circuit.
To characterize the region in the (c 12 , c 6 ) plane where bistability occurs, we used numerical continuation to calculate a co-dimension two limit curve. The code in get bifurcations.jl in our repository takes advantage of the Julia package PseudoArcLengthContinuation.jl [7]. To calculate the bifurcations, we only need consider the steady states of the model. The auto-fluorescence equations are independent of the others and CFP and YFP are simply readouts of c L and c T respectively, leaving only four coupled equations to solve, defined by state vector c = (c R , c S , c L , c T ). For simplicity of analysis, we assume that cell density ρ is constant, and consequently the specific growth rate γ(ρ) =: γ 0 is constant. Later, we check this assumption in Figure 26. With this simplification, the model given by (21) can be represented compactly as where u = (c 6 , c 12 , c I , c A ) are the experimental control parameters, θ is the vector containing the inferred parameters from Section S1 and γ 0 .
To improve the numerical stability of numerical continuation, we transform the model into log 10 coordinates via the element-wise transformations c → 10 c and u → 10 u yielding dc dt = F θ (10 c , 10 u ) 10 −c ln(10) The steady states are defined by zeros of the numerator of the right-hand side. We can immediately see that this transformation induced a zero at c → ∞ which we are not interested in and thus simply seek to solve F θ (10 c , 10 u ) = 0 (24) To further increase numerical stability of finding 24 we explicitly calculate the Jacobian in the log coordinate system. Luckily the Jacobian in the new coordinates can be expressed in terms of the matrix product between the Jacobian in the original coordinates and the Jacobian of the coordinate transformation where Diag[v] is a diagonal matrix with the vector components of v along the diagonal. The Jacobian in the original coordinates is where partials of inhibitory hill functions and promoter activities are and partial bound molecules for uniform and protected degradation models are respectively ∂B X ∂X = 2c X (K X6 c 6 ) n X + (K X12 c 12 ) n X (1 + K X6 c 6 + K X12 c 12 ) n X ∂B X ∂X = 2c X c n X 6 + (E R12 c 12 ) n X X = R c n X 12 + (E S6 c 6 ) n X X = S Given the rate function F θ and its Jacobian ∂F θ ∂c in log coordinates we can perform a co-dimension one parameter continuation for a fixed value of c 12 along the c 6 direction to find a limit point. Then the solution can be continued along a limit curve in the (c 6 , c 12 ) plane along both directions until the limits of the observation region are met. Figure 25 reveals these curves for different models and inferred maximum likelihood parameter sets θ and Figure 26 reveals that the chosen model for the main text is insensitive to changes in growth γ 0 and therefore we can safely assume that the qualitative behaviour of the model will not change if the cell density is constant. For simplicity, this is what is done in the spatial simulations in Sections 2.3 and 2.4.

Boundary Experiments Simulation and local equilibria
For spatial simulations a simple forward-Euler method is implemented in get movie.py. The bacterial colonies for the spatial experiments were placed on top of agar with no signalling molecules inside it. Then additional volumes of agar were attached either side of the width of the experiment with different concentrations of c 6 and c 12 . The signalling molecules then diffuse in the agar and established a cross-gradient felt by the bacterial colonies. These are governed by diffusion where D 6 = 1.8 · 10 −6 m 2 .h −1 and D 12 = 0.9 · 10 −6 m 2 .h −1 . The initial conditions for spatial simulations are zero everywhere except for c 6 (x, t)| t=0 and c 12 (x, t)| t=0 initialised in small regions widths w on opposite sides of the experiment of width W at concentrations such that the homogeneous equilibrium after diffusion would be c 6 (x, t)| t→∞ = C6 and c 12 (x, t)| t→∞ = C12 with zero-flux boundary conditions. Therefore where H(x) is a unit step function.
Supplementary Figure 27: Initial conditions c 6 (x, t)| t=0 and c 12 (x, t)| t=0 with W = 1.6 cm and w = 0.4 cm Each location x experiences concentrations c 6 , c 12 which define a local equilibrium for the remaining state variables. These local equilibria evolve over time and are chased by the actual concentrations of protein in the cells at that location. As shown by Figure 25, at some concentrations c 6 , c 12 there may two stable equilibria rather than one. Therefore it becomes useful to not only display the dynamics in the one dimensional spatial domain but also in the (c 6 , c 12 ) plane. How and when local equilibria bifurcate reveals the eventual fate of the spatial pattern. Figure 28 shows a snapshot of the dynamics in the spatial domain x and the (c 6 , c 12 ) plane. A sharp boundary in space forms when the state density in (c 6 , c 12 ) moves into the bistable region enclosed by the limit point curve. This means that, given the cross-gradient initial conditions 29, if the homogeneous equilibrium C6, C12 lies within the bistable region, a sharp stationary boundary will form. If C6, C12 lies below the cusp of the limit point curve, only soft boundaries will form. In other regions outside the bistable cone the sharp boundary has a finite velocity and will eventually leave the experimental region. See supplementary movies 2-5 for examples of the above. This motivated the experimental exploration of the space of C6, C12 and measurements of boundary velocity, which are described in the following section.

Computation of the boundary velocity
The velocity of the boundary is determined using get movement.py from the TIFF image stack M[t, x, y, s] obtained by the fluorescence microscope from one experiment, set up with a chosen C6, C12 combination. Here t indexes the time point, x and y index the width and height and s indexes the three channels: CFP,YFP and RFP.
First the data are normalised by the RFP channel. This way the location of the boundary can be defined by comparing the pixel values of one channel against the other. The pixels are masked for the colony grid squares (x, y) ∈ Ω which are otherwise surrounded by hydrophobic ink filter paper, on which no colonies grow. The grid squares are detected by thresholding the constituent RFP channel at the end time point.
With the continuous estimate f φ * (x) for each [t, s], the location of the boundary can be obtained even if it was estimated to lie between two colonies as seen in Figure 29. Since this continuous estimate is obtained for each point in time, the position of the boundary can be tracked in a smooth kymograph as shown in Figure 30. We define the boundary location β[t] to be where the estimate of the CFP channel • is equal to that of the YFP channel •, The distance travelled by the boundary ∆β from its formation time t * to the end as a fraction of the size of the experiment W can now be computed. The formation time t * was judged by eye and seems to lie between 3 − 5 h, at which fluorescence values are sufficiently steep to form a sharp boundary. The boundary should have travelled at least 10% along the width -which is the approximate size of one colony grid square -in order to be classified as moving. The distance travelled ∆β can be investigated for different equilibrium concentrations C6 and C12. Figure 31 shows results for a two dimensional dilution between 5 nM and 25000 nM and fixed 10 µM IPTG. The subsequent classification using ∆β is shown as Figure 3c in the main text.

Use of IPTG to influence bifurcation curve
Solid culture experiments on boundary movement were initially performed in the absence of IPTG. We observed stationary boundaries at the concentrations labelled as grey points ( Figure  S24). The shape of the region encompassing the grey points was qualitatively similar to the conical bistability region we observed in previous experiments and in our models, but was quantitatively shifted such that even very low concentrations of C6 enabled bistability. We hypothesized that this was due to minor differences in culture conditions between solid and liquid cultures. We hypothesized that the addition of a low concentration of IPTG would partially derepress LacI, thereby making the YFP-dominant region larger and more like that seen in liquid culture. We used IPTG to shift the bistability region back ( Figure S25) to coincide with the region in liquid culture. This allowed us to perform the solid culture experiments in a regime in which we could observe the transition from stable boundary to moving boundary with YFP dominance, by using higher concentrations of C12 while keeping C6 constant. Experiments with the relay circuit (Fig. 4b) performed as expected without the need for addition of IPTG or ATC. Experiments with the relay circuit (4a) performed as expected without the addition of IPTG or ATC.