Turing Patterning Using Gene Circuits with Gas-Induced Degradation of Quorum Sensing Molecules

The Turing instability was proposed more than six decades ago as a mechanism leading to spatial patterning, but it has yet to be exploited in a synthetic biology setting. Here we characterize the Turing instability in a specific gene circuit that can be implemented in vitro or in populations of clonal cells producing short-range activator N-Acyl homoserine lactone (AHL) and long-range inhibitor hydrogen peroxide (H2O2) gas. Slowing the production rate of the AHL-degrading enzyme, AiiA, generates stable fixed states, limit cycle oscillations and Turing patterns. Further tuning of signaling parameters determines local robustness and controls the range of unstable wavenumbers in the patterning regime. These findings provide a roadmap for optimizing spatial patterns of gene expression based on familiar quorum and gas sensitive E. coli promoters. The circuit design and predictions may be useful for (re)programming spatial dynamics in synthetic and natural gene expression systems.

One mechanism by which ensembles of cells could self-organize is the Turing instability [2,12,13] that occurs due to interplay of short-range activation and long-range inhibition.This instability then drives the formation of spatially periodic patterns.The Turing instability has been implicated in morphogenetic processes of amoebae [14], plants [15,16], and animals [9,10,[17][18][19][20].The large number of unknown factors often makes it challenging to elucidate the essential determinants of morphogenesis in biological systems, but the instability has also been directly engineered in low-component chemical reactions (malonic acid [21] and platinum surface [22]).
Despite the ubiquity of Turing patterns at the multi-cellular scale, they have yet to be demonstrated in gene expression sytems.This is surprising given the many alternative pattern forming mechanisms found in natural [3,4,6] and engineered [5,6,[23][24][25][26][27] colonies of cells.For example, researchers have created gene circuits to produce stationary ring patterns in growing colonies of bacteria [5,24,28].However, these stationary patterns required colony expansion and very particular initial conditions to form.Hsia et al. [29] have presented a theoretical model of a quenched oscillator circuit with one diffusible element that produces a Turing-like instability.There however, the temporal dynamics never settled to a stable fixed state due to persistent temporal oscillations.Another group [30] has recently developed a model of Turing pattern circuit using single promoter.However, even with common promoters, the activator and inhibitor genes would result in different production rate functions, especially if the compounds are to have largely differing diffusion rates.The lack of convincing examples raises the question of whether Turing patters can really be produced by genetic circuits.
If it is possible, then Turing patterning of gene expression could be quite useful for basic science and biotechnology applications.These include self-organized spatial sequestration of gene expression and downstream metabolites (as is currently done inside cells to improve yield [31]), spatially-patterned delivery of proteins [32], reprogrammed morphogenesis of cell populations [33], filtering of paracrine inputs from adjacent cells in populations [34], and chemical sensing [27].Biosensing applications are particularly promising inasmuch as the Turing instability is inherently capable of amplifying small differences in initial conditions.
Motivated by this, we propose a realistic synthetic gene circuit that implements repressoractivated degradation using activating homoserine lactone and inhibitory peroxide gas, and demonstrate how this system can be switched between stable uniform states, limit cycles, and Turing patterns.We characterize the effects of experimentally relevant parameter changes on the presence and properties of the expected Turing patterns.We conclude by discussing our plans to experimentally implement the circuit and to refine the mathematical model that describes it.

Methods
The fixed points, eigenvalues and unstable wavenumbers of the differential equations were computed using the symbolic and variable-precision solvers in Matlab [35].Bifurcation analysis was performed using the Matcont continuation package [36], and analysis of the Turing instability was performed as stated in the text, with a custom Matlab script.Spatio-temporal simulations were performed using custom Fortran code [37] using forward Euler with five point discretized Laplacian.Convergence of solutions was validated by simulating successively smaller spatial and temporal step sizes.For most cases we also tested both periodic and no-flux boundary conditions and used dx = 11.5 μm and dt = 10 −5 min.We rescaled space by assuming D L = 5.6 × 10 −3 μm 2 /min.S1

Turing pattern gene circuit with homoserine lactone and hydrogen peroxide gas
Based on our experience with using quorum-sensing signaling in synthetic circuit design [26,27], we decided to use acyl homoserine lactone (AHL) as the short-range activator signal and hydrogen peroxide gas (H 2 O 2 ) as the long-range inhibitor signal (Fig 1).The choice to focus on a pure activator-inhibitor system [38] (where the activator self-activates and inhibitor selfinhibits) was based on the fact that the AHL quorum sensing system has a native positive feedback, and H 2 O 2 does not.Several of the circuit components presented in Fig 1 have already been used to construct gene circuits that produced sustained temporal oscillations of fluorescent proteins [27].The synthase LuxI produces AHL, and the Ndh protein generates H 2 O 2 .AHL binds with LuxR and activates the p lux -like promoters (p 1 and p 2 in Fig 1).The luxR promoter (p 3 in Fig 1) is constitutive (unregulated).Controlled degradation of AHL is mediated by the AHL-lactonase, AiiA [39].The aiiA gene can be activated by H 2 O 2 by putting it under the control of a promoter (p 4 in Fig 1) such as: p topA , which responds to H 2 O 2 via the Fis pathway [40], or a synthetic promoter containing the ArcAB binding site of p lux , which relieves inhibition by transcription factors ArcAB in the presence of the peroxide [41]).The compound chemical reactions (transcription, translation, protein binding processing) underlying the proposed gene circuit are: where L, H, A, L R , P represent AHL, H 2 O 2 , Aiia, LuxR and the AHL-LuxR complex, respectively.Production is controlled by the thresholds for production of each species.Degradation is first order for H 2 O 2 and Aiia, but enzymatic for AHL, with g(L) = A/(k LD + L).Total LuxR is assumed to be constant, and the AHL-LuxR complex (with relatively fast association and dissociation rates k 1 , k 2 ) is assumed to be relatively stable to degradation.Including mass-action dynamics of these reactions along with diffusion, and rescaling x by where P m is the total amount of LuxR.The parameters of Eqs ( 1)-( 4), and their nominal values used are listed in S1 Table.

Reducing Degradation Feedback Generates Limit Cycles and Turing Patterns
Restricting AiiA dynamics to reasonable relative production and degradation rates [26,27], we find parameter regimes of stable fixed points, limit cycles, and Turing patterning.The maximal Aiia production rate, α 3 , can be adjusted with plasmid copy number, ribosome binding site strength, or by using a hybrid promoter tuned with an inducer.Numerical analysis of the eigenvalues of the stable fixed point in Eqs (1)-( 4) shows that slowing this maximal production rate leads to a stable limit cycle oscillation through a supercritical Hopf bifurcation, followed by a Turing instability, and a saddle node bifurcation (Fig 2a).The Turing instability bifurcation point is found by examining the characteristic equation, for eigenvalues λ, where J is the Jacobian, I is the identity matrix, D is the matrix of diffusion coefficients, and k is the wavenumber.This equation is solved for different α  [42][43][44].Lowering the AiiA degradation rate (or the slowing of AiiA dynamics in general) leads to a similar sequence of transitions from a homogeneous steady state to temporal oscillations followed by Turing patterns, but this was not investigated in detail, as such parameters would be more difficult to control experimentally.
Overall, the results indicate that slowing AiiA feedback on AHL degradation can lead to limit cycle oscillations before the onset of Turing patterns.This is consistent with previous studies showing that an increased time delay in a negative feedback signal produces limit cycle oscillations [45][46][47].It is interesting to note that for the cases shown in Fig 2, the fixed point is already unstable at k = 0 (a limit cycle) when undergoing the transition to the Turing instability.This means that the classic notion of Turing patterns arising from a stable fixed point can be extended to include some cases with an unstable point inside a limit cycle.

Controlling Range of Unstable Wavenumbers and Patterning Robustness
We next investigate how additional can be used to experimentally tune spatial patterning in the circuit.In contrast to the single promoter circuit model that was proposed [30], increasing the cooperativity of transcriptional activation in our model does not appear to increase the parameter space leading to Turing patterns, since an increased rate of activation is not automatically balanced by an increased rate of repression.This however is not a major concern since increasing cooperativity would be challenging to implement in experiments.
A more experimentally tuneable parameter is α 1 , affecting the maximal rate of LuxI expression, and experimentally controlled by a regulatable hybrid promoter, plasmid copy number, or the strength of the promoter or ribosome binding site.Covarying this parameter with the AiiA production rate creates basins for stable, Hopf unstable and Turing unstable solutions (Fig 4a).For α 3 <0.16min.−1 there is no α 1 that gives rise to Turing patterning.For the parameter space investigated (Fig 4a) as α 1 is increased, the range of α 3 leading to Turing patterning (α 3 -robustness) also increases.Raising α 3 also tends to increase α 1 -robustness (Fig 4a and 4b).Fig 4b demonstrates that raising the AiiA production rate expands the range of wavenumbers corresponding to the Turing instability, Re(λ(k)) > 0 towards higher wavenumbers (shorter wavelengths).Thus, tuning the maximal AiiA production rate to achieve more α 1 -robust Turing patterning also results in shorter spatial periods in the resultant pattern.
Another parameter to control experimentally is P m , representing the amount of LuxR.At the nominal parameter set (S1 Table ), the Turing instability occurs for P m > 0.18 (Fig 4c).Overall, the P m -robustness of Turing patterning is not very sensitive to α 1 , except near the borders of the Turing space (Fig 4c).Furthermore, the α 1 -robustness of patterning is not much affected by increasing P m (Fig 4c and 4d).It is clear however that a higher amount of LuxR leads to much shorter spatial periods in the Turing pattern (Fig 4d).Both codimension two bifurcation plots (Fig 4a and 4c) show that the Turing regime can be adjacent in parameter space to a regime of limit cycle dynamics (as in Fig 2), but this is not always the case (for example, decreasing α 1 when 0.5 < α 3 < 3.2, in Fig 4a).
Together these results show that Turing patterning with this circuit requires an AHL production rate high enough to prevent a stable focus or limit cycle, but not so high as to lose the fixed point through a saddle node bifurcation.Furthermore, it appears that the AiiA production rate is the highest-priority parameter to adjust in order to achieve Turing patterning that is robust to variations in the other experimental parameters.Once inside the Turing unstable regime, both LuxI production rate and the amount of LuxR can be used to tune the range of spatial periods in the solution.

Discussion
In this paper we have demonstrated the possibility of generating Turing patterns using a longrange gas signal that induces degradation of a short-range quorum sensing molecules.The next important step is to experimentally implement this proposed circuit in and test the predictions made by the model.The gene circuit has features in common with previous designs [26,27], but there are some important differences.The Aiia promoter would likely need to be redesigned to be activated by exclusively by H 2 O 2 instead of LuxR-AHL.Besides the use of ndh, H 2 O 2 generation could also be increased by addition of a lux-induced sodA gene [27].For the lux-like promoters (p 1 and p 3 in Fig 1 ), the model shows that they should be distinct from each other, with leak and maximal production from p 1 being higher than those from p 3 (compare δ 1,2 's and α 1,2 's in S1 Table ).Alternatively, this could be achieved by putting ribosome binding sites on luxI that are stronger than those on ndh.It may also be necessary to remove the ArcAB sites on these promoters to eliminate crosstalk with the H 2 O 2 part of the circuit, but preliminary analysis shows that crosstalk high as 20% can still lead to patterning (see S1 Equations) and S3 Fig.
Pattern formation speed is relevant to producing spatial structures in time scales allowed by experiments.This speed is reflected by the magnitude of Re(λ) > 0, that determines how fast unstable solutions diverge from the fixed point.Although the high dimensionality of our model makes it challenging to derive analytical expressions for the pattern formation speed (as was done in [30,53] for two-dimensional models), computing dispersion relations for different parameters demonstrates that pattern formation speed around the nominal parameter set is inversely related to α Bacteria engineered with this gene circuit could be tested on agar plates or inside microfluidic chambers, both of which are amenable to spatiotemporal data collection with an epifluorescence microscope.In one spatial dimension, one could use microfluidic devices like the long trap [26], or an annular trap [54] to study periodic boundary conditions.Agar could be patterned into one-dimensional spatial structures, but these are likely to exhaust nutrients or lose moisture too quickly.Spatial patterns in two dimensions are easier to generate on agar plates (as in [5,23,24]), because use of microfluidics would likely be complicated by diffusional anisotropies resulting from constant perfusion of liquid media through the channels.If this were not a challenge then one could exploit the gas-permeability of polydimethylsiloxane (commonly used to fabricate microfluidic devices) to allow for selective diffusion of H 2 O 2 (as in [27]), or sequester the gas using vacuum channel in the microfluidic device (another use for the vacuum channel that was originally built for shear-free loading of cells [55]).
Another potential challenge to implementing these results experimentally arises from the assumptions made about diffusion.First, transport through and between the cells is only approximated by the diffusion equation.Limited membrane permeability and cell crowding can slow transport of signaling molecules, resulting in a smaller effective diffusion coefficient [56].The effective diffusion coefficients of AHL and H 2 O 2 have not been measured directly for bacterial monolayers in microfluidic devices.However, Danino and colleagues [26] measured conduction velocity of AHL-mediated activation to be V = 8.5-35 mm min , which combined with the model prediction of V = 0.17 D 1=2 L , results in D L = 2.5-42 × 10 3 mm 2 min .This estimate is comparable to D L = 5.9-29 × 10 3 mm 2 min (lower limit for biofilm; upper limit for water at 25°C) estimated by Stewart [56].For H 2 O 2 , D H = 72-120 × 10 3 mm 2 min [56], implying that the diffusion ratio, D = 4-12.However, Prindle et al. measured H 2 O 2 -mediated fluorescence conduction velocity V % 700 mm min across the biopixels microfluidic device [27], so using the 0.17 prefactor that Danino [26] found demonstrates that the ratio, D = 100, used in our model is reasonable.In addition, decreasing D as much as fourfold from D = 100 can still allow for patterning to occur (S1 Fig) .If required, further disparity between gas and lactone diffusion could be attained by using crowding agents like polyethylene glycol, or engineering circuits with slower diffusing quorum sensing molecules [48].
The mathematical model presented could be further investigated by incorporating several additional layers of detail.This includes a more biophysical description of AHL production (LuxI catalyzes the conversion of SAM to AHL), as well as H 2 O 2 production (NDH-2 catalyzes formation of H 2 O 2 ), and regulation of the aiiA promoter (for example, H 2 O 2 oxidizes ArcA to relieve repression of lux promoter [27]).All of these processes have been lumped together (represented by dashed arrows in Fig 1) in the current model, but may turn out to have a non-negligible effect on dynamics.Our model also lumps the processes transcription, translation and protein maturation.Thus it would also be interesting to include these processes and also compare this system to one reduced to a spatially-extended delay differential equation.Previous work [57,58] suggested that delays can predispose that system to oscillatory patterning, which is consistent with what we found upon the inclusion of slower AiiA dynamics (Figs 2a and 4a).Another feature that the current model ignores is the effect of protein dilution due to cell growth.While this may hinder implementation in exponentially growing cells, there is also work suggesting that an expanding domain can help make the Turing patterning regime more robust [57].Furthermore, the dilution problem can be fundamentally resolved by engineering the chemical reactions in an in vitro transcription-translation system [59].
Considering that protein numbers are low at the steady states analyzed (the nanomolar range corresponds to tens of molecules per E. coli), it would be prudent to further examine this system using stochastic simulations.Hsia and colleagues [29] found that such stochastic effects changed the patterning parameter parameter space, but did not destroy the ability to produce Turing patterning.Although external fluctuations tend to destabilize the Turing parameter regime in a continuum setting [57], it has also been shown that when space is discrete, either intrinsic or extrinsic noise can cause patterning over larger parameter ranges [60,61].Testing these additional considerations is a good direction for future modeling studies.
Despite these theoretical and experimental challenges in further studying the proposed circuit we believe it is well worth the effort, as these types of gas and quorum sensing circuits could be quite useful for studies of bacterial morphogenesis and several synthetic biology applications.The present work provides a conceptual foundation and guidance for the development of gas mediated spatial patterning in multi-cellular and cell-free gene expression systems.

Fig 1 .
Fig 1.The proposed gene circuit for generation of Turing patterns.AHL activates its own production, but is degraded by Aiia.H 2 O 2 is produced via the ndh gene, and activates the transcription of the aiiA gene.Intercellular transport and diffusion of AHL and H 2 O 2 are represented by the thick arrows.The circuit is modeled by Eqs (1)-(4).doi:10.1371/journal.pone.0153679.g001 3 , with the largest positive real parts of the eigenvalues plotted in Fig 2b.Having only two diffusable substances leads to a quadratic equation in k 2 , A (λ) (k 2 ) 2 + B (λ) (k 2 ) + C (λ) = 0. Finding the root of the discriminant of the quadratic, B 2 ðl¼0Þ À 4A ðl¼0Þ C ðl¼0Þ , gives the condition for the Turing instability boundary (α 3 = 3.1 min.−1 in Fig 2a)).Numerical simulations of the system in one spatial dimension with periodic boundary conditions further confirm a stable focus at α 3 = 3.5 min.−1 (Fig 3a), a stable limit cycle at α 3 = 3.2 min.−1 (Fig 3b), and Turing instability at α 3 = 2.5 min.−1 (Fig 3c).In two spatial dimensions the instability leads to formation of spot patterns (for snapshots of simulations with periodic boundary conditions, see Fig 3d, and for movies with no-flux boundary conditions, see S1 Movie).Our observation of Turing instability taking over from Hopf instability has previously seen in other systems

Fig 2 .Fig 3 .
Fig 2. Slowing AiiA production in Eqs (1)-(4) leads to oscillations and Turing patterns.A codimension one bifurcation of the AHL fixed point, L * , losing stability through a Hopf bifurcation and into a Turing instability as Aiia production rate is decreased.(b) Eigenvalue-wavenumber curves at various AiiA maximal production rates, corroborating the bifurcation analysis results above.At the cusp of each curve the eigenvalues become a complex conjugate pair, with each low eigenvalue left off the plot for clarity.doi:10.1371/journal.pone.0153679.g002

Fig 4 .
Fig 4. Dependence of Turing parameter space and wavenumbers on quorum sensing and AiiA parameters.(a,c) Two-parameter sections of Turing space, with curves of limit points (green), Hopf points (red), and Turing points (blue).(b,d) Ranges of Turing unstable wavenumbers for different sections along the maximal LuxI production rate within the Turing parameter space.Parameter values are varied the nominal set (S1 Table).doi:10.1371/journal.pone.0153679.g004

3 (
Fig 2b), but proportional to D (S1 Fig), α 1 and P m (S2 Fig).Although pattern formation speed will typically be very slow near the transition to Turing instability, our simulations within the Turing parameter space demonstrate formation of patterns within two hours (Fig 3c and 3d), making the system quite accessible to experimental observation.
Table lists nominal parameter values used in the model equations.