Stochastic Turing patterns in a synthetic bacterial population

Synthetic Biology of a Bacterial Community.

In our synthetic gene network design, which was guided by computational modeling (SI Appendix, section 1), we used two artificial diffusible morphogens: the small molecule N-(3-oxododecanoyl) homoserine lactone, denoted here as A3OC12HSL, and the small molecule N-butanoyl-l-homoserine lactone, denoted here as IC4HSL, from the Pseudomonas aeruginosa las and rhl quorum sensing pathways, respectively, in P. aeruginosa (15). A3OC12HSL serves as an activator of both its own synthesis and that of IC4HSL, while IC4HSL serves as an inhibitor of both signals (Fig. 1 A and B and SI Appendix, section 1). A3OC12HSL activates its own synthesis and synthesis of IC4HSL by binding regulatory protein LasR to form a complex that activates the hybrid promoter PLas−OR1. This promoter regulates expression of LasI, an A3OC12HSL synthase, and RhlI, an IC4HSL synthase. To increase the sensitivity of A3OC12HSL self-activation, LasR is regulated by a second copy of PLas−OR1. IC4HSL inhibits synthesis of A3OC12HSL and itself by forming a complex with the regulatory protein RhlR. This complex activates expression of lambda repressor CI, which in turn, represses transcription of LasI, RhlI, LasR, and RhlR. Pattern formation in our system can be modulated by altering the concentration of isopropyl β-d-1-thiogalactopyranoside (IPTG), a small molecule inducer that binds LacI and alleviates repression of PRhl−lacO. GFP and red fluorescent protein (RFP) are expressed from the rhl and las hybrid promoters, respectively, to aid in experimental observation (SI Appendix, section 2).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Design of a synthetic multicellular system for emergent pattern formation. (A) Abstractly, the system consists of two signaling species A3OC12HSL and IC4HSL. A3OC12HSL is an activator catalyzing synthesis of both species, while IC4HSL is an inhibitor repressing their synthesis, with additional repression by A3OC12HSL via competitive binding. (B) Genetic circuit implementation. Promoter regions are indicated by white boxes, while protein coding sequences are indicated by colored boxes. IPTG is an external inducer modulating system dynamics. (C, Top) Illustration of signaling species concentrations in 1D space. The dashed orange and blue lines correspond to A3OC12HSL and IC4HSL, respectively. (C, Middle) Spatial profiles of reporter proteins. RFP expression (red line) correlates with A3OC12HSL concentrations, while GFP expression (green line) roughly mirrors RFP expression. (C, Bottom) A vertical slice of cell lawn. Cells express fluorescence proteins according to the profiles above and produce a global multicellular pattern.

In our experimental setup, the A3OC12HSL activator diffuses more slowly than the IC4HSL inhibitor (SI Appendix, section 3). The estimated diffusion coefficient for A3OC12HSL is 83 μm2/s, and for IC4HSL, it is 1,810 μm2/s. The experimentally determined ratio of diffusion rates in our system of 21.6 is much higher than the value of 1.5 predicted by Wilke–Chang correlation in water (16), likely due to partitioning of A3OC12HSL in the cell membrane, which slows its diffusion from cell to cell (17). The slower diffusion rate of A3OC12HSL coupled with positive feedback regulating its synthesis allows A3OC12HSL to aggregate in local domains, leading to formation of visible red fluorescent spots (cellular lawn illustration is shown in Fig. 1C). Within these red domains, both A3OC12HSL and IC4HSL are found in high concentrations, but because A3OC12HSL competitively binds RhlR (SI Appendix, section 4 and Fig. S8), GFP is attenuated (18). The faster diffusion rate of IC4HSL allows it to diffuse into regions outside of the red fluorescent domains. Here, IC4HSL is free to bind RhlR, activating GFP expression. Collectively, these processes lead to green regions between red spots.

Experimental Patterns and Controls.

To study pattern-forming behavior, engineered cells are first cultured in liquid media, and then, they are concentrated and plated on a petri dish to form an initially homogeneous “lawn” of cells (Materials and Methods). After plating, the petri dish is incubated for 24 h at 30 °C, and microscope fluorescence images are captured as needed. Before the self-activation of the A3OC12HSL synthase positive feedback loop, the cell lawn exhibits no fluorescence. However, over time, red fluorescent spots emerge with sizes much larger than that of a single cell (10–1,000×). Simultaneously, green fluorescence develops in a pattern with dark voids positioned precisely in the locations of the intense red fluorescence (Fig. 2A). Time series microscopy reveals that patterns begin to emerge after approximately 16 h (SI Appendix, Fig. S12).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Experimental observations of emergent pattern formation measured by green fluorescence intensity (GFI) and red fluorescence intensity (RFI). (A) Representative microscope images (based on six technical replicates) of a typical field of view showing a fluorescent pattern formed by an initially homogeneous isogenic lawn of cells harboring the Turing circuit with no IPTG. Spots and voids appear in the red fluorescence (RF) and green fluorescence (GF) channels, respectively. (Scale bar: 100 μm.) (B) Microscope images of cell lawns with constitutive expression of fluorescent proteins. (Left) Cells expressing RFP, (Center) cells expressing GFP, and (Right) mixed population of red and green cells. (C) Fluorescence density plots computed from the images above (from left to right: red, green, red/green, and Turing). Color intensity is in log scale [arbitrary units (a.u.)].

In control experiments, we show that our patterns are not simply a result of the outward growth of clusters of differentially colored cells (Fig. 2 B and C and SI Appendix, section 5). Also, by performing an experiment with cells that harbor independent bistable green/red toggle switches, we test whether observable patterns would emerge if individual cells autonomously made cell fate decisions at some point after plating (SI Appendix, section 5). The fluorescence fields after 24 h of incubation at 30 °C in both control experiments are uniform, showing no emergence of patterns.

Next, we examine how changes in the strengths of localized interactions lead to different global outcomes in our pattern-forming gene circuit. In our system, IPTG can be used to modulate the inhibitory efficiency of IC4HSL in individual cells by affecting CI expression from PRhl−lacO. Specifically, IPTG relieves LacI repression of CI and GFP reporter. The increased range of CI ultimately increases inhibition of both morphogens, which is expected to decrease activator spot sizes, while causing the field of CI and GFP reporter to be more strongly expressed. Our data show that mean GFP levels increase sigmoidally with inducer concentration, while the overall area of red spots decreases (Fig. 3 A and C). In addition to offering further support that our gene circuit gives rise to emergent patterns, these results show how pattern formation characteristics can potentially be tuned to fit future application needs.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Mathematical modeling and correlation between pattern modulation experiments and simulations. (A) Experimental results for IPTG modulation of pattern formation with microscopy images corresponding to specific IPTG concentrations in B. The same display mappings were used for all images in A. (B) Collectivity, metric parameter Θ is influenced by IPTG modulation. (C) Pattern statistics over IPTG modulation for experimental results. (D) Pattern obtained from simulating a deterministic reaction–diffusion model with Dv/Du=100. (E) Pattern statistics over IPTG modulation for deterministic modeling. (F) Patterns obtained from simulating our deterministic model (Upper) and stochastic spatiotemporal model (Lower) at the measured diffusion ratio of Dv/Du=21.6. (G) Pattern statistics over IPTG modulation for stochastic modeling. a.u., arbitrary unit.

Deterministic Model.

We first developed a detailed deterministic reaction–diffusion model (SI Appendix, section 7). The model explicitly describes chemical reactions for the LasI and RhlI synthases, regulatory protein CI, and synthesis and diffusion of the morphogens A3OC12HSL and IC4HSL. As the overall system involves a large number of reactions with rate constants that span multiple timescales, we made two commonly used simplifying assumptions. First, we assume that operator states of a promoter fluctuate much faster than protein degradation rates. Second, we assume that mRNA half-life is much shorter than protein half-life. These assumptions allow us to eliminate operator fluctuation and mRNA kinetics and model the system at the communication signals and protein levels as follows:∂U∂t=αuIu−γuU+Du∇2U[1]∂V∂t=αvIv−γvV+Dv∇2V[2]∂Iu∂t=αiuF1(X1,C)−γiuIu[3]∂Iv∂t=αivF1(X1,C)−γivIv[4]∂C∂t=αcF2(X2,L)−γcC,[5]where U and V are the concentrations of the two diffusible morphogens A3OC12HSL and IC4HSL, respectively; Iu and Iv are the concentrations of corresponding acyl homoserine lactone (AHL) synthases, respectively; and C refers to CI.

We model the hybrid promoters using the following Hill functions:F1(X1,C)=1+f1X1Kd1θ11+f2−1CKd2θ21+X1Kd1θ11+CKd2θ2[6]F2(X2,L)=1+f3X2Kd3θ31+f4−1LKd4θ41+X2Kd3θ31+LKd4θ4,[7]where F1(X1,C) and F2(X2,L) are the production rates of the promoters PLas−OR1 and PRhl−lacO, respectively; X1 and X2 are the LasR-A3OC12HSL complex and the RhlR-IC4HSL complex, respectively; and L is the concentration of unbound LacI protein. We use the definitionsX1=RuU[8]X2=RvV(1+U/Kc3)[9]L=λl1+f6−1I/Kd6θ61+I/Kd6θ6,[10]where I is the IPTG concentration and Ru and Rv are the regulatory proteins LasR and RhlR, respectively:Ru=λuIu[11]Rv=λv1+f5−1(C/Kd5)θ51+(C/Kd5)θ5.[12]A summary of the variables used in our model is available in SI Appendix, Table S5, and definitions of the rate constants are in SI Appendix, Tables S6 and S7. Hill functions used in this model have a shared form of Y=1+fX/Kθ1+X/Kθ, where X and Y correspond to the input and output of the function, respectively; K is the dissociation constant; θ is the Hill coefficient; and f is the fold change of Y on full induction by X.

We initially ran simulations of this model using a high diffusion rate ratio (DvDu) of 100. These simulations yield patterns of red spots and green voids (Fig. 3D), suggesting that the underlying dynamics of our system are Turing-like, with the potential for Turing instabilities. Deterministic simulations of IPTG modulations also correlate well with the trends of the experimental results (Fig. 3E and SI Appendix, section 6).

While the overall behavior of our system is reminiscent of classical Turing patterns (19), there are key differences. In particular, when we ran simulations at the measured diffusion rate ratio of DvDu≈21.6, patterns did not arise (Fig. 3F). For some two-node implementations of Turing systems, this rate would be sufficient for pattern formation (20). In addition, certain networks with more nodes can allow small or even equal morphogen diffusion rate ratios to generate Turing instabilities (21). However, a practical biological implementation imposes certain dynamics, such as delays associated with protein production, that can strongly impact pattern formation (22, 23). Indeed, our deterministic modeling results suggest that the ratio of diffusion constants for the activator and inhibitor in our system is either barely within the range required for a Turing instability or even outside the range, depending on the precise medium in which signal diffusion is measured (SI Appendix, section 7). In addition, whereas in the deterministic simulation, spots are identical and evenly distributed, those in the experimental systems vary in size, shape, fluorescence intensity, and the intervals between them.

Stochastic Model.

The deterministic modeling results indicate that our system may be beyond the regime where classical Turing patterns are formed but still within the regime where stochastic Turing patterns occur (5⇓⇓–8). Indeed, gene expression in microbes is inherently noisy due to the small volume of cells and the fact that many reactants are present in low numbers, suggesting that stochastic Turing patterns could be present in our system (24).

Noise in stochastic Turing patterns expands the range of parameters in which patterns form, in contrast to the usual expectation that noise serves as a destabilizing agent. The patterns observed in stochastic Turing systems correspond to the slowest decaying mode of the fluctuations. Similar noise stabilization phenomena can be observed in other systems that are out of equilibrium. For example, in predator–prey systems, fluctuations can drive temporal oscillations of populations (25, 26). Noise-driven stabilization has also been recently discovered in the clustering of molecules on biological membranes (27, 28) and in models that exhibit Turing-like pattern formation (7). In particular, whereas spatial symmetry breaking and pattern formation via the original Turing design require two morphogens with diffusion rates that differ by a large factor on the order of 10 or 100 (1), the requirements to form stochastic Turing patterns are less stringent. For example, in a pattern-forming plankton–herbivore ecosystem, the noise associated with discrete random birth and death processes reduces the required ratio of diffusion constants for pattern formation from a threshold of 27.8 for normal Turing patterns to a threshold of 2.48 for stochastic Turing patterns (5⇓⇓–8).

To determine whether noise in the chemical reactions underlying gene expression and morphogen diffusion in our system can cause the emergence of patterns over a wider range of parameters than a deterministic model, we constructed a stochastic spatiotemporal model using the same biochemical reactions, diffusion, and rate constants used in our deterministic model (SI Appendix, section 8). This model captures stochastic effects in the production and degradation of the proteins and morphogens in our system but approximates diffusion as deterministic. Simulations of the stochastic model generically produce patterns with large variability in spot size, shape, intensity, and intervals, which are similar to the patterns observed in our experiments and different from those predicted for the deterministic model (Fig. 3F). We have compared the experimental patterns with stochastic simulations in both real space and in 2D Fourier transform (2DFT) space (SI Appendix, section 8). Neither the experimental 2DFT nor the simulated 2DFT contain pronounced peaks that would be present in a deterministic honeycomb Turing pattern. Moreover, as the IPTG concentration is increased, both experimental and simulated patterns become more regular (Fig. 3G).

To further test the hypothesis that we are observing stochastic Turing patterns, we measured the power spectrum for both of our fluorescent reporters. Theory predicts that the power spectrum will have a power law tail as a function of wavenumber, k, for large wavenumbers, with an exponent characteristic of the noise source (7, 25). The exponent values are −2 and −4 for stochastic Turing patterns and deterministic Turing patterns with additive noise, respectively, and can be interpreted simply as follows. The −2 arises, because at small frequency or wavenumber, the random variable (i.e., concentration) is simply diffusing and therefore, follows the behavior of a random walk, which has a power spectrum that exhibits a −2 power law. The −4 arises, because for a system that is executing deterministic damped periodic motion but driven by additive white noise, the response of the random variable is a Lorentzian, with an asymptotic behavior for the power spectrum that exhibits a −4 exponent.

For the GFP channel, we observe a power law tail with an exponent of −2.3±0.4 (Fig. 4B and SI Appendix, section 9A). For the RFP channel, we also observe a power law tail with an exponent of −3.9±0.4 (SI Appendix, section 9A). To better understand the implications of these tails, we examined our detailed stochastic model of the genetic systems and also developed a reduced stochastic model that explicitly includes only the morphogens (SI Appendix, section 9B). Both models predict that our experimental parameters will produce a stochastic pattern with a power law tail of −2 for both the activator and the inhibitor at asymptotically large wavenumbers (Fig 4C and SI Appendix, section 9C). However, in the range of parameters likely to correspond to the experiments (Fig. 4A and SI Appendix, section 9C), the detailed stochastic model predicts that the exponent of the power law tail for the activator will be −4 over a large range of intermediate wavenumbers before it eventually undergoes a cross-over to a power law with an exponent −2 at high wavenumbers (SI Appendix, Fig. S24). This behavior once again agrees with our experimental data and supports our identification of stochastic Turing patterns. In summary, spectral analysis of the patterns of activator and inhibitor is consistent with a model in which fluctuations in the amount of signaling morphogens drive stochastic Turing patterns.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Spectral analysis and parameter analysis. (A) Pattern-forming regimes in parameter space and estimated parameters for our system. Parameters above the green surface of neutral stochastic stability can form stochastic patterns, and parameters above the blue surface of deterministic neutral stability can form deterministic Turing patterns. The ratio of the diffusion coefficients ν/μ, the ratio of degradation rate to production rate d/p, and the ratio of production rates are estimated for our system by the yellow ellipsoid. The parameters for our system are mostly in the regime where stochastic patterns form and outside the region where deterministic Turing patterns form. Example stochastic simulations are shown for parameters drawn from a deterministic parameter region with Dν/Dμ=100 (Upper Right) and a stochastic region with Dν/Dμ=21.6 (Lower Right). (B) Radial power spectrum of green fluorescence and best fit power law tail with an exponent of −2.3±0.2. (C) Radial power spectrum for eight trials of our stochastic simulation, their mean, and the best fit power law tail.

Our analysis of the stochastic Turing model predicts that stochastic patterns form over a wide range of parameters (SI Appendix, section 8). Indeed, our stochastic model predicts that stochastic Turing patterns are possible at the measured ratio of diffusion rates for A3OC12HSL and IC4HSL (Figs. 3F and 4A). In addition, to determine the sensitivity of the stochastic model to the parameters chosen, we individually varied parameters from 0.5× their nominal value to 1.5× their nominal value while keeping all other parameters fixed at their best estimated value. For each set of parameters, we calculate the analytical power spectrum and the eigenvalues of the Jacobian (linear stability matrix) of the stochastic model evaluated at a fixed point found numerically. Based on this analysis (SI Appendix, section 8), we classify each set of parameters as producing an unstable homogeneous state at wavenumber k=0, a stable homogeneous state, a stochastic Turing pattern, or a deterministic Turing pattern. Specifically, we classify a set of parameters as producing a pattern if they produce a peak in the calculated power spectrum at a nonzero wavenumber. To distinguish between stochastic Turing patterns and deterministic Turing patterns, we examine the eigenvalues of the corresponding Jacobian. If the real part of all of the eigenvalues is negative for all wavenumbers, then the pattern must be due to stochasticity. If there is any range of wavenumbers that have corresponding positive real parts of their eigenvalues, then the pattern is produced by the traditional Turing mechanism. The results of this analysis are shown in SI Appendix, Fig. S22 and illustrate the significant ranges for each parameter that can lead to stochastic Turing patterns. Indeed, the estimated parameter values yield stochastic Turing patterns and variation of Du, Dv, and IPTG, and several other parameters never produce deterministic patterns; therefore, our results are very insensitive to estimation error of these important parameters. Overall, varying the parameters one at a time, 68% of the values yield stochastic Turing patterns.

To quantify the way in which stochasticity enlarges the pattern-forming regime of parameter space, we simultaneously varied all model parameters and performed the classification used above. Specifically, we used Latin hypercube sampling to randomly generate 500 parameter sets, where all of the parameters were allowed to vary between 0.5× and 1.5× their nominal value. For this analysis, we found that 24.8% of parameters produced unstable fixed points, 43.2% produced stable homogeneous states, 13.2% produced stochastic Turing patterns, and 18.8% produced Turing patterns. Thus, over this arbitrarily large range of parameters, pattern formation occurs only 18.8% of the time in the absence of stochasticity but 32% of the time when stochasticity is included. By including stochasticity, the range in which patterns can form has been increased by 70%.