Cooperativity to increase Turing pattern space for synthetic biology.

It is hard to bridge the gap between mathematical formulations and biological implementations of Turing patterns, yet this is necessary for both understanding and engineering these networks with synthetic biology approaches. Here, we model a reaction-diffusion system with two morphogens in a monostable regime, inspired by components that we recently described in a synthetic biology study in mammalian cells.1 The model employs a single promoter to express both the activator and inhibitor genes and produces Turing patterns over large regions of parameter space, using biologically interpretable Hill function reactions. We applied a stability analysis and identified rules for choosing biologically tunable parameter relationships to increase the likelihood of successful patterning. We show how to control Turing pattern sizes and time evolution by manipulating the values for production and degradation relationships. More importantly, our analysis predicts that steep dose-response functions arising from cooperativity are mandatory for Turing patterns. Greater steepness increases parameter space and even reduces the requirement for differential diffusion between activator and inhibitor. These results demonstrate some of the limitations of linear scenarios for reaction-diffusion systems and will help to guide projects to engineer synthetic Turing patterns.

The self-organization of spatial patterns and structures is a fundamental problem in many fields, especially in developmental biology. There are various forms of patterning required for the development of an organism from a zygote, but a common type of patterning requires regular periodicity, such as in making multiple stripes or spots. The most influential model to explain periodic patterns is the reaction–diffusion model proposed by Alan Turing in the 1950s. Although first suggested by Rashevsky,[2] it was Turing who actually demonstrated the periodic or complex reaction–diffusion patterns that came to be known as Turing patterns.[3] Turing showed that a coupling between diffusion and reaction can lead to instabilities from initially homogeneous paired reagent systems, resulting in the spontaneous emergence of regular patterns. Including seminal contributions from Gierer and Meinhardt in the 1970s,[4] a large body of literature subsequently developed the modeling of self-organizing biological patterns, based on reaction–diffusion theory (see refs (5−8)).

Although other mechanisms of pattern formation exist,[7,9] in this article we will focus on the Turing mechanism, which is unique in terms of self-correction and economy of design, making it a tantalizing engineering target in the emerging field of synthetic biology.[10] The Turing mechanism can be thought of as a competition between self-activation by a slow-diffusing chemical (or biological) activator and inhibition by a faster-diffusing factor.[11] This has been largely applied to explain self-regulated repetitive pattern formation in developing animal embryos.[12−16]

Over the last few decades, research in nonequilibrium chemical systems has led to experimental confirmations of the spatial patterns formed by Turing instabilities.[17−21] However, the experimental verification of Turing-driven developmental phenomena has been elusive in biological systems,[12,22] while other patterning strategies have been more frequently found in development.[23] For example, the formation of the characteristic periodic striped pattern of pair-rule gene expression in Drosophila(7,24) follows a strategy based on a hierarchical regulatory gene network, embedded in a morphogen gradient, rather than a Turing pattern.

Even though many biological patterns resemble Turing patterns in appearance, this alone does not constitute a compelling proof of such mechanisms actually operating in living systems.[12] Alternatively, due to recent advances in biotechnology, the synthetic generation of Turing patterns in cell culture should be possible in the near future.[25,26] This will increase our fundamental understanding of these systems as well as provide tools for biotechnology applications. To implement such engineered biological systems, we need to know a priori the requisite properties of the gene and protein building blocks. These properties mediate the processes that support Turing patterning formation, such as production (through specific regulatory functions), diffusion, and degradation. We have been developing a scaffold for biological network engineering with diffusing activators and inhibitors[1] and decided to implement a computational model to help guide the design of our artificial networks. The key aim, therefore, was to develop a biologically interpretable model that would show us the most flexible parameter relationships for making patterns.

Although Turing instabilities are almost ubiquitous in studies modeling highly repetitive patterns in developmental biology, the regulatory functions used in the modeling are often selected by criteria that simplify the equation analysis rather than being chosen on the basis of their correlation with the actual biological response.[13,16,27−29] These simplifications give us phenomenological descriptions and can lead to model-induced constraints over mathematical parameters, with these constraints being mandatory for patterning to occur. These affected parameters encompass several biophysical processes such as diffusion, regulation, and degradation, all of which underly Turing pattern formation. Unfortunately, such mathematical parameters are not independent of each other and, subsequently, the imposed constraints conform poorly to the actual properties of the available biological building blocks. Despite the extensive literature on Turing patterns in biology, very few studies have considered sigmoidal regulatory functions. In a recent study on hair follicle development, Sick et al. showed that a sigmoidal function for noncompetitive inhibition is also able to display Turing patterns.[30] Despite being phenomenological, sigmoidal functions are much closer to the real behaviors of gene expression systems and are thus more relevant.

In the present article, we study the parameter space where Turing instabilities can occur, in reaction–diffusion systems whose reaction terms involve regulatory functions with greater biological interpretability. To achieve this, we have performed a linear stability analysis to find the constraints on parameters that allow pattern formation. We believe that these findings can be used to guide the engineering of biological systems able to form Turing patterns using synthetic biological scaffolds.[1] Importantly, we find that the cooperativity associated with the regulatory function is a key factor. Together with the differential diffusion of activator/inhibitor, cooperativity determines the size of the parameter region associated with successful patterning. This, in turn, can be be used to predict which properties yield robustness and engineering flexibility.

Results and Discussion

Two-Morphogens Turing Model

In general, a biological system able to present Turing instabilities can be modeled by coupled reaction–diffusion equations of the formwhere a and h denote the spatial concentration of the activator and the inhibitor morphogens, respectively. The functions f and g correspond to the regulatory functions of the genes that encode the morphogens a and h. Such functions describe how the expression (or production) rate depends on the concentration of activator a and inhibitor h (see the sketch of the model in Figure 1A). The last term on the right-hand side of each equation describes the degradation process that is assumed to be linear. In our biological implementation, the activator corresponds to hepatocyte growth factor (HGF), and the inhibitor is a truncated variant of HGF, named NK4.[1] Both the activator and inhibitor are expressed and secreted into the extracellular medium from Madin–Darby canine kidney (MDCK) cells, grown as cysts in 3D collagen cell culture.

Figure 1

Modeling overview. (A) Schematic diagram of the biological system at the cellular and molecular levels: activators (green particles; HGF) and inhibitors (red particles; NK4) are products of genes controlled by the same promoter (MMP-1). These molecules are released into the extracellular medium. Activator and inhibitor molecules compete for binding to membrane receptors (yellow; cMet). Activators bound to membrane receptors regulate the production of morphogens through the regulatory function f. (B) The regulatory function used in this article. The plots show f vs activator concentration a for three different concentrations of inhibitor h: 0.01 (blue), 0.25 (red), and 1.0 (brown). nH = 2, k = k = 0.1.

We will consider here the case that both morphogens are under the control of the same promoter, which corresponds to the matrix metalloproteinase promoter-1 (MMP-1) in our synthetic system.[1] Mathematically, this means that both regulatory functions f and g will be the same, so, hereafter, the regulatory function will be denoted by f. Additionally, we consider that the regulatory functions have sigmoidal forms. In particular, we consider that the regulatory function of the activator obeys a Hill function in the absence of inhibitor. For the inhibitor, there are two alternatives: First, the inhibitor competes with the activator for the same regulatory site, decreasing the effective affinity (see Figure 1B). Thus, the regulatory function is well-described bywhere nH is the Hill exponent that describes the steepness of the response to the activator, while mH is related to the inhibitor. As we are regarding the case that both morphogens are under the control of the same promoter (to simplify downstream biological engineering), we consequently consider that the steepnesses of the responses to the activator and to the inhibitor are the same, i.e., nH = mH. Biologically, this corresponds to HGF or NK4 binding competitively to the extracellular cMet receptor and activating or repressing an intracellular signaling cascade that controls expression of the MMP-1 promoter. The second alternative is that the inhibitor blocks the activation pathway, decreasing the effective production rates (see Figure S1). This type of noncompetitive inhibitor case was used recently for modeling hair follicle development, and Turing patterns were reported for nH = 2, D = 0.005, and D = 0.2.[30]

The Dimensionless Problem

Before determining under which conditions the homogeneous steady state becomes unstable, leading to inhomogeneous Turing patterns, it is convenient to define the following new variables. Thus, the new variables of time t̂, position x̂, and the concentrations â and ĥ can be written in terms of the variables and parameters contained in eqs 1–3 aswhere k and k are related to the effective dissociation constant for the activator and inhibitor, respectively. Rewriting eqs 1 and 2 with these new variables, we getwhere we have introduced the abbreviations D = D/D, μ = μ/μ, r = ρk/μ, and r = ρk/(μk). f̂ denotes the dimensionless form of the regulatory function eq 3 and can be written as . Hereafter, we will write new variables â, ĥ, and the function f̂ without the hat in order to simplify the notation.

These reaction–diffusion equations, with any of the specified regulatory functions, have more parameters and more complex dependences on morphogens than previous ones.[16,27−29] Because it is very difficult to understand the dynamics of such systems directly, a possible approach is to use a linear stability analysis. For example, it is possible to derive a closed expression for the wavenumber dependency with diffusion constants and the derivative of the linear reaction terms.[27] Of course, we can approximate any real reaction terms to a linear one, but such a compromise leaves a gap to bridge, which is essential for engineering a real biological device. It is necessary to find the link between the derivative and the actual biological parameters that fully characterize the reaction term. These include morphogen degradation rates, production rates, cooperativity, and dissociation constants; reconciling them is not straightforward. Despite its complexity, our system is amenable to the typical linear stability analysis, in an analytical fashion, for some particular values of Hill coefficient nH.

The Dispersion Relation

In order to determine the conditions for the development of Turing patterns, we will consider the dimensionless version of the two-morphogen system that evolves in a one-dimensional domain of length L. Given the reduced version of the reaction–diffusion equations, the stationary homogeneous solutions, denoted by (a,h), correspond to the solution of the system in the absence of diffusionThe (a, h) can achieve Turing instabilities depending on the parameter values. Our main goal is to obtain the region in the parameter space where the instability occurs for the regulatory functions f (see eq 3). As is usual, we perturb the solution slightly according to a = a + δa and h = h + δh, where δa = δ e(ω and, similarly, δh = δ e(ω. The wavenumber k is the number of complete waves per 2π units of distance, while the real part of ω can be interpreted as the growth rate. Because morphogens cannot diffuse beyond the boundaries x = 0 and x = L, we need to impose a zero-flux boundary condition, which implies that k takes discrete values k = nπ/L, with n = 0, 1, 2, .... Next, we introduce a + δa and h + δh into eqs 4 and 5 and linearize the system around the point (a,h), which giveswhere f and f are derivatives of the regulatory function f with respect to a and h, respectively, evaluated in the point (a, h).

The perturbation amplitudes δ and δ can be different from zero if and only if the discriminant of is zero. Then, the condition for a nontrivial solution is given by ω2 + αω + β = 0. Thus, we obtain a function, known as dispersion relation, which associates a growth rate ω to each k given bywhere α and β are given by

Because ω–(k) is negative in the region of interest of the parameter space, we will analyze the dependency of the largest root ω+(k) with the wavenumber k. The stationary homogeneous state becomes unstable when the real part of ω+ is positive for some wavenumber k. In the linear approximation, the value of the wavenumber associated with the largest positive amplification rate dominates the evolution of the instabilities; we will denote this wavenumber as kmax. kmax is an important feature of the resulting patterns because of its relation to their size. The typical size of an emergent pattern is given by π/kmax such that higher kmax implies low wave sizes and thus smaller patterns.

If the size of domain is much larger than the typical pattern size (i.e., L ≫ π/kmax), then we can go one step further by considering the discrete variable k as a continuous variable k. In this way, we can obtain an expression for kmax by maximizing ω+Thus, in the framework of the linear stability analysis, it is possible to determine a closed expression for kmax as a function of the parameters setwhere f and f can be written in a closed expression in terms of the parameters μ, r, and r for some values of nH. However, eq 11 alone does not allow us to determine if the Turing pattern associated with kmax will be stable. As we will see later, transient Turing patterns can also occur for some parameter values. The growth of Turing patterns also depends on the size domain L, which, in any case, must be greater than the typical pattern size (i.e., L > π/kmax). This condition imposes a critical length for the domain size L, which implies a constraint by the organism size for Turing patterns to develop.[6]

In the next section, we will analyze some properties of the Turing instabilities.

Parameter Region for Turing Patterns

For further steps in the stability analysis, we need to specify the regulatory function and the Hill coefficient. The system in the absence of diffusion (eqs 6 and 7) has two fixed points for nH > 1; one of these fixed points corresponds to a stable solution, and the other one, to an unstable solution. For nH ≤ 1, the unstable solution disappears, and one cannot expect the occurrence of Turing instabilities. We are able to find closed expressions for the fixed points for nH = 4/3, 3/2, 2, and 3. These closed expressions are large and complicated algebraic expressions in most of cases, with the exception of nH = 2. Thus, we will illustrate the complete analysis for this case, but the same steps can be done for the above-mentioned Hill coefficient values.

To look for an unstable state, we analyze the stability of all fixed points by examining the eigenvalues of the Jacobian matrix. For the case nH = 2, we find that the solution given byhas one positive and one negative eigenvalue, corresponding to an unstable state known as a saddle point. From expressions 12 and 13, we can determine the derivatives f and f, which in turn feed the expression of kmax (11).

The unstable state determined by eqs 12 and 13 allows Turing instabilities in a region of the parameter space spanned by D, μ, r, and r, which are terms related to combinations of diffusion rates, degradation rates, and binding affinities. Because we are dealing with four-dimensional space, we can illustrate only some projections of this region on 3D space. In particular, we are interested in the projections on the space spanned by parameters μ, r, and r for some values of D. In Figure 2, we depict the projection of the region of the parameter space where kmax > 0, i.e., the region where Turing patterns can be developed, for two different values of the ratio of diffusion of activator/inhibitor, D (0.01 and 0.1). Biologically, this corresponds to situations where the inhibitor diffuses 100- and 10-fold faster than that of the activator, respectively. As expected, the volume of this region decreases with D, meaning that it is easier to find Turing patterns when the inhibitor diffuses much faster than the activator (100-fold faster is better than 10-fold faster; see also Table 1).

Figure 2

Parameter region for Turing patterns. Region of the parameter space where Turing instabilities develop (i.e., where kmax is positive) for nH = 2 and D = 0.01 (A) or D = 0.1 (B). The region of pattern-forming parameter space increases when D decreases (i.e., a 100-fold difference in inhibitor diffusion, relative to slower activator diffusion, is more likely to yield Turing patterns than is a 10-fold difference). Turing space also increases with the steepness of the Hill function nH (see Figures 7 and 8).

Table 1

Summary of Parameter Relationships To Engineer Turing Patterns

parameters for Turing patterns	current system value	range of values simulated	preferred parameter value and range	ideas to engineer-implementation	refs
differential diffusion (D = Da/Dh)	1.0	0.01–0.75	0.5, Da < Dh	HGF fusion with collagen binding domains or multimerization domains of any protein.	(38, and 41−44)
differential degradation (μ = μh/μa)	1.0	0.45–1.5	0.5, μa > μh	TEV protease targeted to HGF	(45)
Hill coeff. (nH)	1.2	1.33–6.0	5, nH > 1.33	Having both HGF and cMET-receptor in positive feedback, HGF sequestration on collagen.	(46 and 47)
dimensionless activator prod. rate (ra = ρakan–1/μa)	unknown	2.9–6.0	6.0, ra > 5.0	Greater production of activator (e.g., alternative promoter) for larger, more rapidly forming patterns	(39)
dimensionless inhibitor prod. rate (rh = ρhkan/μhkh)	unknown	5.0–6.0	5.0, 5 < rh < 6	NK4 variants with lower activity for larger more rapidly forming patterns. Greater production of inhibitor (e.g., alternative promoter) for smaller, slower forming patterns	(37 and 39)

Figure 7

The effect of cooperativity. Region of the parameter space where kmax is positive for D = 0.01: nH = 3.0 (A) and nH = 1.333 (B). This result suggests that the presence of cooperativity in the regulatory function will enhance the possible parameter space where Turing patterns occur.

Figure 8

Patterns with low differential diffusion D by using higher Hill coefficients: simulation results obtained for D = 0.5, nH = 5, μ = 0.7, r = 5.5, and r = 6 (A) and for D = 0.75, nH = 6, μ = 0.8, r = 6, and r = 6 (B). The D term reflects the ratio of the diffusion constants of activator/inhibitor, and patterns can be achieved with relatively small differences so long as the response function is steep, i.e., cooperative. Note that different spatial scales are used.

In Figure 3, we show two cross-sections of the volume displayed in Figure 2 (D = 0.01) with two different planes, where the heatmap is related to kmax scale. In panel A, the intersecting plane is r = 5, and panel B corresponds to r = 5. From this plot, we can conclude that the convex side of the volume is related to higher values of kmax and consequently to patterns with smaller typical length. Because μ is the ratio of the degradation rates of inhibitor relative to activator, increasing inhibitor degradation or decreasing activator degradation will increase the pattern size. Alternatively, increasing the production of activator (or decreasing the production of inhibitor) can also increase the pattern size. Figure 4 depicts the numerical solution of the eqs 4 and 5 for the parameter values indicated by the white dots in Figure 3B, i.e., D = 0.01, r = 5, r = 5, and three different values of μ: 0.45, 1, and 1.5. We can see that the size of the pattern increases with μ, as predicted.

Figure 3

Size of Turing patterns. Intersection of the volume of Figure 2A with plane r = 5.0 (A) and r = 5.0 (B); the color scale depicts kmax values. For a higher wavenumber kmax, one expects smaller patterns, due to the typical length L being given by L = π/kmax. Turing patterns for parameter values corresponding to the white dots in panel B are shown in Figure 4.

Figure 4

Turing patterns and morphogen degradation. Striped patterns obtained by numerical integration of eqs 5 and 6, using the same parameters as Figure 3B (white dots), for three different values of μ: 0.45 (A), 1.0 (B) and 1.5 (C). In all cases, the initial condition is a small Gaussian perturbation of the unstable steady state at x = 25. In agreement with Figure 3B, the size of the pattern in the simulation increases with μ. Because μ is the ratio of the degradation rates of inhibitor relative to activator, increasing inhibitor degradation or decreasing activator degradation will increase the pattern size. For μ = 1.5 (C), the typical Turing pattern decays to the stable solution after its initial formation, at t = 75, as a consequence of a field size effect.

For μ = 1.5, the simulation shows that the system, initially in a slightly perturbed unstable homogeneous solution, reaches a stable homogeneous solution after a transitory Turing-like pattern generation. It seems that for small μ there does exist a stable periodic pattern to which the system is attracted once the Turing bifurcation occurs, whereas for larger μ, such a stable pattern does not exist. Interestingly, transient Turing patterns have been observed in chemical systems where the transiency is due to the depletion of chemical species in a closed system.[31,32]

In Figure 5, we can see the growth rate ω+ as a function of the k for different values of μ. The position of local maximum kmax changes with μ (the ratio of the degradation rates of inhibitor relative to the activator). From Figure 5, the maximum growth rate of patterns, denoted by ωmax = ω+(kmax), seems to increase with μ for this particular set of values. In general, the maximum growth rate of patterns depends on the parameter set D, μ, r, and r. Figure 6A,B depicts how ωmax depends on these parameters for cross-sections of the same volume represented in Figure 3 for D = 0.01. These plots show that smaller patterns grow more slowly. This prediction is in agreement with the simulation shown in Figure 6C, which depicts the temporal evolution of the concentration of a at the position x = 25.

Figure 5

Speed of pattern appearance versus pattern size. Behavior of ω+ as a function of k for n = 2, D = 0.01, r = 2.5, r = 7.0, and three different values of μ: 3.75 (green), 4.5 (red), and 6.5 (blue). ω+ presents only one maximum in kmax, whose particular value depends on the parameter set (D, μ, r, and r). For this particular set of parameter values, we can see that speed of pattern appearance, ωmax, increases with the typical size of the pattern.

Figure 6

Growth rates of patterning. (Top) Similar to that in Figure 3, we see the intersection of the volume of Figure 2A, with plane r = 0.5 (A) and with the plane r = 5 (B), but now the color scale depicts ωmax values instead of kmax. ωmax quantifies the speed of pattern appearance; for higher values, one expects faster-forming patterns. (C) Pattern formation dynamics for the concentration of a at x = 25. Parameter values correspond to the white dots in panel B (D = 0.01, nH = 2, r = 5, r = 5, and with three different values of μ: 0.45 (green line), 1.0 (red line), and 1.5 (blue line)).

Turing Patterns and Cooperativity

To further analyze the effects of the D diffusion term and μ, we must also consider the role of cooperativity in the patterning processes. Biological regulatory functions are usually sigmoidal responses that arise by cooperative binding of transcription factors (TFs)[33] or by molecular titration.[34] The steepness of the sigmoidal stimulus–response curve is modeled and controlled by the Hill coefficient nH. Unlike with the other parameters, it is not possible to obtain closed expressions for the fixed points in terms of the exponent nH. However, we obtained closed expressions for some particular values of nH, which allows us to unveil the effects of this important type of nonlinearity.

In Figure 7, we depict projections of the region where the system develops Turing instabilities for different values of the Hill coefficient nH. Again, we are interested in the projections on the space spanned by parameters μ, r, and r for D = 0.01. We can see that the projected volume increases when the cooperativity increases. In particular, for the nH = 3 case (panel A), lower values of μ ensure a large domain in the (r, r) plane, where Turing instabilities can develop. By contrast, for the smaller nH value tested (nH = 4/3, panel B), we find that this region is very small, almost reduced to μ < 5 and a narrow range of the parameter r.

Looking at Figures 2 and 7, we can see that the region of the parameter space where Turing instabilities are expected decreases or increases in size as a function of D and nH. However, there is not much change in the overall position of the region itself. This feature allows us to explore the behavior of the system for higher nH and similar diffusivity for both morphogens, i.e., where D is closer to 1. Figure 8 depicts Turing patterns obtained for the case where the inhibitor diffuses only slightly faster than activator D = 0.5 and 0.75, using higher Hill coefficients. Figure 8A uses D = 0.5 and nH = 5, whereas in Figure 8B, we used D = 0.75 and nH = 6 (see Supporting Information for the code for simulations).

Overall, steeper regulatory functions lead to larger parameter sets compatible with Turing patterns, which should facilitate design decisions for synthetic gene network engineering. Importantly, even situations with small differential diffusion between activator/inhibitor can support Turing patterns in combination with high Hill coefficients (Table 1).

Until now, all of our simulations were done for dimensionless equations (eqs 4 and 5). To translate the results to the spatial and temporal scales of a real biological system, we can transform back to variables with units by multiplying the dimensionless lengths by the factor (D/μ)1/2 and dividing the dimensionless times by μ. Thus, the real typical size of a pattern is given byAssuming that the diffusion constants of inhibitor is similar to our previous experimental estimation for HGF (DHGF = 10–3 mm2/min) and the activator HGF half-life in the extracellular medium is around 21 h,[1] we can predict that the typical repeat size of the pattern is around 4–6 mm. Hence, the time needed to reach fully developed Turing patterns on the domain we have explored (≈50 dimensionless units in Figure 4) is estimated to be around 100 h. These scales are compatible with the system we have described, where MDCK cysts can be grown for up to 2 weeks in glass-bottomed dishes measuring several centimeters in diameter.

A final question was whether altering the relative effective production rates of activator r and inhibitor r would alter pattern formation. In Figure 9, numerical simulations show patterns arising for the effective inhibitor production r = 5.0 and three different values of effective activator production r: 2.9 (A), 3.5 (B), and 4.2 (C). We can see that increasing activator production r increases the pattern size and its speed of formation, in agreement with Figures 3 and 6. This suggests that using different HGF-responsive promoters, with a stronger response for the inhibitor, could tune the system for smaller and slower-evolving patterns.

Figure 9

Turing patterns with parameters to increase potential Turing space in the context of low cooperativity. Simulation results obtained for D = 0.01, nH = 1.333, μ = 0.5, and r = 5 and with three different values of r: 2.9 (A), 3.5 (B), and 4.2 (C).

In summary, we found a series of relationships that qualitatively guide our design decision making for the downstream engineering of Turing pattern systems, using our MDCK system.[1] In the following discussion, we will explore the ways that some of these properties could be engineered.

Discussion

Alan Turing’s hypothesis that an inhibitor and activator can create regular patterns, emerging from uniform but noisy initial conditions, was one of the first attempts at a theoretical explanation of biological pattern formation.[3] Importantly, his system employed faster diffusion and decay of the inhibitor, relative to that of the activator. Today, it is broadly accepted that two requirements, namely, local self-enhancement and long-range inhibition, are essential features of this type of pattern formation.[5] The results presented here suggest one more additional requirement, not so broadly present or acknowledged in previous studies, which is for cooperativity in the regulatory function response.

Several studies have reported Turing-like patterns in developmental biology, but, even now, it is still not absolutely clear whether many real molecular systems employ activators and inhibitors exactly in this way. Excitingly, a recent paper elegantly demonstrated that Nodal/Lefty proteins have differential diffusion properties that underpin a reaction–diffusion patterning system in zebrafish embryogenesis.[35] It is interesting that, in this system, the relative protein half-lives were only a minor contributor to the differences in the range of these factors, with differential diffusivity contributing the most to patterning. For a deeper understanding of these systems that enables bottom-up engineering, it is important to define which properties of the network components are the most critical and which are more forgiving of variation.

In this article, we investigated how the regions of parameter space where Turing instabilities occur are affected by specified parameters that can be engineered biologically. We were motivated to explore the simple diffusing activator and inhibitor system that we recently developed.[1] This system is genetically encoded and is based on MDCK cells secreting an activator, hepatocyte growth factor (HGF), and a truncated variant that acts as an inhibitor (NK4). Both molecules then bind the same membrane receptor (c-Met) extracellularly and thus activate or repress the same matrix metalloproteinase-1 gene promoter (MMP-1), via a signaling cascade, which in turn produces more activator and inhibitor. Although these components have been shown to function in cell culture,[1] they have not yet been connected to form a Turing-like feedback network.[36] To work toward this goal, we were primarily interested in developing a biologically interpretable model to explore whether these components would be at all amenable for supporting Turing pattern formation in a simple, single-promoter implementation. Furthermore, we aimed to identify the parameter properties associated with the largest pattern-forming parameter space, with respect to the other parameters. We thus aimed to define the most robust regions of design space for attempting to engineer artificial Turing patterns in cells.

Our main findings can be summarized as follows (see also Table 1):

(1) We found that putting both morphogens under the control of the same promoter (easier to engineer) still allowed the formation of Turing patterns in the context of our competitive activator–inhibitor system. This could be achieved in principle by expressing both HGF and NK4 with the MMP-1 promoter, resulting in a combined positive–negative feedback system. It should be noted that a previous model[30] also used a single promoter; however, that model was within a framework of noncompetitive inhibition.

(2) As expected, the ratio of diffusion of inhibitor/activator affects the size of the parameter regions associated with Turing patterns. For example, 100-fold faster inhibitor diffusion relative to that of the activator is more likely to achieve patterns than 10-fold faster. We assume that the apparent diffusion constants of HGF and NK4 in collagen are similar. Practically, there are several ways to increase this ratio. First, truncated NK4 variants (NK1 and NK2) may diffuse faster through collagen, although this benefit may be at the cost of reduced cMet repression activity.[37] Second, fusing matrix-binding domains to HGF could retard its diffusion, as was achieved recently with a number of growth factors that were engineered to have superaffinity to extracellular matrix.[38] Fusing multimerization domains could also reduce HGF diffusion, in principle, and increase potential Turing space.

(3) Steeper regulatory function responses, potentially achieved by greater cooperativity, lead to larger potential Turing space. Interestingly, the MMP-1 promoter that we use[1] can be fitted with an apparent Hill coefficient of only nH = 1.2. This is toward the lower end of the range that supports Turing patterns, and the system may benefit from a search for steeper induction responses, with respect to HGF morphogen concentration. Practically, this would mean searching for promoters that respond more steeply to HGF in a dose–response curve and then defining the minimal promoter and enhancer region that captured this function in a synthetic promoter construct. Not only do higher Hill coefficients increase the likelihood of Turing space but also they reduce the requirement for differential diffusion: with nH = 6, we were able to use D = 0.75, which implies that the activator need diffuse only 3/4 as fast as the inhibitor.

(4) Larger patterns (fewer complete waves per spatial interval) develop faster over time, in agreement with ref (27). This can be tuned by varying the relative ratios of inhibitor or activator production and degradation, as follows: (4.1) Increasing the production or activity of activator or decreasing the activity or production of inhibitor can also increase the pattern size and rate of emergence of the patterns. In practice, this might be achieved under the single promoter by adding tandem copies of activator, HGF, linked by internal ribosome entry sites (IRES) or the viral self-cleaving 2A peptide signal or by using NK4 variants with lower activity (e.g., NK1 and NK2). (4.2) Increasing inhibitor degradation or decreasing activator degradation will increase the pattern size and, consequently, the patterns will develop more rapidly. Assuming that NK4 has a similar half-life to that of HGF, an artificially higher degradation function could be achieved by adding site-specific tobacco etch virus (TEV) protease cleavage sites in the linkers between Kringle domains in NK4. TEV protease could then be added to the collagen gel matrix, where the cells are grown, in order to tune this differential degradation function. It should be noted, however, that to facilitate engineering differential degradation is not an absolute requirement.

(5) Using two different HGF-responsive promoters could tune the space and time scale of the system. For example, with a greater production rate of the activator relative to that of the inhibitor, one would achieve larger, more rapidly forming patterns. In practice, the TIMP1 promoter[39] has a reported 200-fold induction with HGF and might be used to drive further HGF production after suitable characterization. By contrast, the MMP-1 promoter fragment we currently use is induced only up to 4-fold with HGF. Assuming similar uninduced production rates, reversing the order of promoter usage would produce more NK4 and lead to smaller, slower-evolving patterns.

Recently, McKane et al. proved that the effect of intrinsic noise translates, in the linear noise approximation scenario at least, into an enlargement of the parameter region yielding Turing mechanisms when compared to a deterministic analysis.[40] This result gives noise a further role in Turing patterning, beyond simply providing the essential initial perturbations that break the homogeneity of the systems. Here, we reported other biological features that could also be essential to support Turing instabilities in a specified biological framework, based on an HGF/NK4 sender–receiver system in Madin–Darby canine kidney (MDCK) cells.[1] We believe that the use of a nonlinear regulatory function that is more biologically interpretable, in combination with stability analyses, offers new insights over the nature of this pattern formation, especially in the context of synthetic biology.

Methods

For the simulation of the dimensionless reaction–diffusion equations (4 and 5), we employed a Mathematica v9.0 script built in-house (see Supporting Information for an example). In all of these cases, the initial condition for the activator consists of the unstable solution plus a small Gaussian perturbation (imperceptible in the plots) in the middle of the spatial domain. The initial condition for the inhibitor was the unperturbed unstable solution h. To explore the effect of the boundary condition, we used both a periodic boundary condition, i.e., a(0) = a(L) and h(0) = h(L), as well as the zero-flux boundary condition, i.e., a(0) = a(L) = 0 and h(0) = h(L). No major differences were observed by using these boundary conditions. For the numerical integration procedure, we use an adaptative temporal step size so that the estimated error in the solution is just within the specified accuracy, 10–6.