A non-linear analysis of Turing pattern formation.

Reaction-diffusion schemes are widely used to model and interpret phenomena in various fields. In that context, phenomena driven by Turing instabilities are particularly relevant to describe patterning in a number of biological processes. While the conditions that determine the appearance of Turing patterns and their wavelength can be easily obtained by a linear stability analysis, the estimation of pattern amplitudes requires cumbersome calculations due to non-linear terms. Here we introduce an expansion method that makes possible to obtain analytical, approximated, solutions of the pattern amplitudes. We check and illustrate the reliability of this methodology with results obtained from numerical simulations.

Introduction

Spatiotemporal pattern formation is an all-important and ubiquitous self-organization phenomenon in nature [1]. Processes in Physics, Chemistry, Biology, and even Social Sciences, display the formation of organized structures arising from the interaction of individual constituents [2-6]. In this regard, Alan Turing’s pioneering studies on pattern formation stands out among various mechanisms [7]. To begin with, it explains how one of the main transport mechanisms in nature with homogenizing properties, i.e. diffusion, can, counterintuitively, lead to the formation of inhomogeneous, ordered, structures. Notably, during the last decades different studies have revealed the widespread biological applicability of Turing’s ideas to different processes, such as animal coating [4, 8, 9], vertebrate limb development [10-12], tooth primordium patterning [4, 13], the ruggae spacing in the mammalian palate [14], the Min oscillations in bacteria [15, 16], or epidemiology and ecology problems [17-20].

The basic ingredients of a Turing instability are antagonistic local interactions between species with distinct diffusivity properties, e.g. activators vs. inhibitors in a chemical context or infected vs. susceptible individuals in epidemiology models,. Yet, not all combinations of the parameters result in a Turing instability (i.e. patterning). The regions in the parameter space leading to pattern formation, as well as the spatial periodicity, are usually determined by a linear stability analysis [4]. However, the estimation of other patterning properties, such as the amplitude or the pattern arrangement, require cumbersome calculations, e.g. the amplitude equations formalism [21, 22]. We point out that the pattern amplitude determines key features such as the maximum protein expression levels or the degree of infection during an epidemic. Consequently, having simplified means to estimate amplitude values and to understand how those are modulated by non-linearities is advantageous. Here we address this problem and introduce an expansion method that results in formulas to estimate the pattern amplitudes. Our analysis focuses, following Cross and Hohenberg’s pattern categorization [2], on the common I instabilities (periodic patterns in space that are stationary in time) and we illustrate our findings using the case of two interacting species.

A non-linear analysis of the Turing instability

For the sake of simplicity, we restrict our analysis to two coupled reaction-diffusion equations in one spatial dimension. Thus, let us consider the following set of equations describing the spatiotemporal dynamics of two scalar fields u = u (x, t) and v = v (x, t), The functions f and g (reactive terms) are supposed to have an equilibrium point, , that defines a steady, pattern-free (i.e. homogeneous) solution such that f(u0, v0) = g(u0, v0) = 0. By defining and , where z stands for either the field u or v, the equilibrium point is stable if,

If a stationary pattern (type I) eventually develops we expect the field z (x) to have a stationary solution with a functional form, where stands for the amplitude of the pattern of the Fourier mode, q. By substituting this ansatz, Eq (3), into Eq (1) the system reduces to a set of ordinary differential equations,

We notice that the solution , i.e. the pattern-free solution, is always a stationary solution of Eq (4). Thus, the appearance of stationary patterns rely on the existence of solutions such that,

Alternatively, since with δ = (δ, δ), this condition can be written as a function of the pattern amplitudes, with .

Assuming that f and g are polynomial functions (or admit a polynomial expansion around ), then Eq (6) read, where p + r > 1 and s + t > 1. By defining τ = ϵ/ϵ, if we further assume that one of the non-linearities is dominant (i.e., |τ| < 1), we can perform a perturbative expansion of the pattern amplitudes such that at orden τ0 Eq (8) read, Eq (10) can be interpreted in terms of a phase between u and v such that if then u and v patterns develop in a counter-phase manner. Eqs (9) and (10) can be solved and we obtain the following solutions for and , At order τ1 Eq (8) read, that lead to the solution,

Thus, up to order the solution for the pattern amplitudes is given by with and as prescribed by Eqs (11) and (13) respectively.

In case a pattern develops we expect to do it continuously, that is, with an amplitude growing from zero. Thus, the above solution must contain the trivial solution δ = δ = 0 that defines the patterning separatrix in the parameter space. The latter implies that a pattern develops if there exist a value of q for which δ = 0, that in turn requires that and leads to the condition . We stress that regardless of P(q2) appearing in the denominator of the expression for (see Eq (13)), one can easily check that no singularity develops due to the dependence of the numerator on . We also point out that lim P(q2) → +∞ and, by construction, P(0) > 0 (see Eq (2)). Thus, since P(q2) has an absolute minimum at , the condition of existence of a patterning separatrix can be stated as P(q*2) = 0, that is, That is, patterns develop if . Moreover, the most unstable Fourier mode of the developing pattern can be approximated by q* close to the separatrix. We point out that inside the patterning region a full band of Fourier modes becomes unstable and the fastest growing mode deviates from q*. More importantly, a coupling between Fourier components develops and such phenomenon is incompatible with the hypothesized functional form of the field (see Eq (3)). Consequently, we expect our results to be strictly valid in the vicinity of the patterning separatrix, and to obtain disagreements away from it. As for the role played by non-linear terms, Eq (14) is the same condition that is obtained using the conventional Fourier decomposition method over the linearized system (c.f. [4, 23]). Consequently, non-linear components are irrelevant for determining the boundaries of the patterning region as expected.

Example: An activator-inhibitor model system

In order to test our analytical predictions, we consider the following example where f and g, as prescribed by Eq (7), read, with a > 0. This family of non-linear equations can be mapped into an activator-substrate model proposed to describe pigmentation patterns in sea shells [24] and into an activator-inhibitor model that accounts for the regeneration process in hydra [25]. Here we considered their linearized parts and included phenomenological non-linear terms with p = 3, r = 2, and s = t = 1 (see Discussion).

In this case, there is a single real equilibrium point of the reactive terms f(u0, v0) = g(u0, v0) = 0: , that defines the steady, pattern-free solution and, since , , , and , then the stability conditions of , as given by Eq (2), are automatically satisfied. Patterns develop if the following conditions hold (see Eq (14)),

Under those conditions, close the separatrix the periodicity of the pattern is given by the Fourier mode q*2 = (a − 1/D) / 2 and the amplitudes by with and as prescribed by Eqs (11) and (13) respectively. Thus, the leading order of the amplitudes reads,

In order to check our predictions, we ran numerical simulations using the Method of Lines (MOL) with periodic boundary conditions [26] implemented through a Wolfram’s Mathematica code [27] (S1 Simulation Code). For every value of the parameters a and D we changed the system size, L, such that 6 pattern wavelengths fitted into the spatial domain, i.e., . The initial conditions for u and v were periodic functions with a periodicity λ to facilitate reaching a stationary state faster.

Fig 1 shows that the comparison between the approximated theoretical amplitudes, up to order , and those obtained numerically in the parameter space a and D. We point out that in Fig 1A the theoretical surface grows continuously from the patterning separatrix but such detail was skipped in the plot to perceive better the boundary line. We observed a fair quantitative agreement (see quantification in Fig 2) and the theoretical solutions captured both the behavior and the values of the pattern amplitudes as a function of a and D. We noticed that for this particular choice of equations and parameters the theoretical solutions always underestimated the exact (i.e. numerical) solution. Yet, we found that for some cases of the non-linear exponents the theoretical solution overestimated the exact solution (see Fig 2).

Fig 1

A: Amplitudes of the pattern δ (left) and δ (right) in the parameter space a-D obtained by our approach (continuous surfaces) and in numerical simulations (circles) for the activator, u (green in all panels), inhibitor, v, (red in all panels) system described in the main text (Eq 15). The black dashed line indicates the patterning separatrix, Eq (14), and the black circle indicates the parameter values used to obtain the pattern shown in Fig 2A. The central inset shows, graphically, the regulatory interactions between the species u and v as defined by the sign of the linear components, , , , and . B: Trajectories to sample the theoretical solutions. The symbols indicate simulations results and the lines the approximated solutions for fixed values of D = 9, 20 (solid and dashed lines respectively) and increasing a, and for fixed values of a = 0.4, 0.9 (solid and dashed lines respectively) and increasing values of D.

Fig 2

A: Stationary pattern solutions, u (green) and v (red), obtained by our theoretical approach (continuous lines) and in numerical simulations (dotted lines). Parameters: a = 0.9, D = 9 (black circle in Fig 1A). B: Pattern amplitudes as a function of the nonlinear parameter, τ, for different values of p and r (log-log axes): p = 3 and r = 2 (circles), p = 3 and r = 0 (triangles), and p = 5 and r = 4 (squares). The symbols stand for numerical results and the colored lines for the analytical solutions: p = 3 and r = 2 (solid line), p = 3 and r = 0 (dotted line), and p = 5 and r = 4 (dashed line). The black solids lines are a guide to the eye and show, from top to bottom, the functions , , and . The insets show the relative importance of the first order correction versus the leading order of the expansion by plotting as a function of τ (lines codes as in the main figure). C: Quantification of the relative error, Δ, as a function of different combinations of the values of the exponents p, r, s, and t (the rest of parameters as in panel A). The scale bar (right) is the same for all density plots.

The pattern solutions obtained by the theoretical approximation for q = q*, Eq (3), are indeed in agreement with those obtained in numerical simulations, Fig 2A. We also checked the accuracy of the theoretical approach to capture the power-law functional dependence of the pattern amplitude with ϵ: note that the theoretical solution predicts that the dominant term of the pattern amplitude behaves as, . To that end, we ran simulations keeping constant ϵ = 1/10 and varying ϵ using different values of p and r. In this way, by plotting the values of the amplitudes as a function of the expansion parameter, τ = ϵ/ϵ, we validated both the scaling behavior and the goodness of the theoretical prediction as τ increases (Fig 2B). The results revealed an excellent agreement and showed that, in these particular examples, the estimations were accurate even for values of τ close to one since the first order correction is very small with respect the leading order .

Finally, we investigated the accuracy of the theoretical approach as a function of different combinations of the values of the exponents p, r, s, and t (Fig 2C). To that end, we quantified the relative error by means of . We found that the error was always below 15% (absolute value). While for the different studied cases was not possible to found an overall trend, we found monotonic behaviors as a function of p and r and data revealed that increasing r systematically decreased the error. Also, as mentioned above, in most cases the theoretical solutions underestimate the exact solution (red-ish colors). Yet, some combination of exponents led to an overestimation of the numerical amplitude (blue-ish colors). To evaluate the relevance of truncating the expansion series at order , we repeated our calculations for the case ϵ = 0. We point out that in that case the leading order, , are the exact solution of Eq (8). As shown in Fig 2B, we did not find a perfect agreement between the numerical solution and the theoretical estimate. In fact, the relative errors for some values of the non-linear exponents are slightly larger than those found when ϵ ≠ 0 (cf. columns Δ, 1 and 3, and columns Δ, 2 and 4). This phenomenon can be explained as follows: as mentioned above, the hypothesis about the functional form of the patterning solution, Eq (3), breaks down as soon as we move away from the patterning boundary since a band of Fourier modes, and not just one as we assumed, becomes unstable. In any case, as suggested by Fig 2B, additional terms in the τ-expansion are not expected to increase significantly the accuracy of the theoretical estimation when ϵ ≠ 0.

Discussion and conclusions

Here we have introduced a framework to analyze I patterns driven by a Turing instability. This approach maps reaction-diffusion systems into algebraic problems and makes possible to estimate approximated patterning solutions. Importantly, our method simplifies the cumbersome calculations required to obtain information about the amplitude of patterns. As for the limitations of our approach, the method assumes that the amplitude that corresponds to the most unstable Fourier mode represents accurately the real amplitude of the pattern. As shown here, this assumption leads to quantitative disagreements that are expected to be relevant far away from the patterning boundary. However, by comparing the predicted patterning solution obtained by this methodology with numerical solutions, we have shown that the former reproduces qualitatively, and, to a great extent, also quantitatively, the actual patterning solutions and their functional dependence on the model parameters. As a matter of discussion, our approach accounts for a single non-linear term in each reaction such that one of them is dominant. While this is certainly a limitation, a number of systems fall within that category. Some of them are models about gene networks [28] where the experimental evidence suggests the leading, linear, interactions and non-linearities are phenomenologically included to induce the saturation of the growing pattern, e.g. digit patterning during limb development [10]. Other models are those arising from enzyme kinetics studies where Hill-like regulatory functions are approximated by polynomial functions: e.g. the Selkov model of glycolysis [29, 30]. Finally, in some models, while the requirement about a single non-linearity in each reaction is fulfilled, the non-linear terms are equally dominant, i.e. τ = 1, for example the Schnakenberg model [31]. Still, given that our approach could still be applied under those conditions, see Fig 2B, we argue that even in those cases our methodology may hold. Possible extensions of the method, out of the scope of this manuscript, may include generalizations to other patterning situations, such as non-stationary patterns, or the effect of stochastic perturbations. Altogether, we hope that the methodology introduced here can be useful to analyze Turing patterning systems in different fields.

Mathematica code to perform simulations of Turing patterns.

(NB)

Click here for additional data file.