Turing's diffusive threshold in random reaction-diffusion systems

Turing instabilities of reaction-diffusion systems can only arise if the diffusivities of the chemical species are sufficiently different. This threshold is unphysical in most systems with $N=2$ diffusing species, forcing experimental realizations of the instability to rely on fluctuations or additional nondiffusing species. Here we ask whether this diffusive threshold lowers for $N>2$ to allow"true"Turing instabilities. Inspired by May's analysis of the stability of random ecological communities, we analyze the probability distribution of the diffusive threshold in reaction-diffusion systems defined by random matrices describing linearized dynamics near a homogeneous fixed point. In the numerically tractable cases $N\leqslant 6$, we find that the diffusive threshold becomes more likely to be smaller and physical as $N$ increases and that most of these many-species instabilities cannot be described by reduced models with fewer species.

In 1952, Turing described the pattern-forming instability that now bears his name [1]: diffusion can destabilize a fixed point of a system of reactions that is stable in well-mixed conditions.Nigh on threescore and ten years on, the contribution of Turing's mechanism to chemical and biological morphogenesis remains debated, not least because of the diffusive threshold inherent in the mechanism: chemical species in reaction systems are expected to have roughly equal diffusivities, yet Turing instabilities cannot arise at equal diffusivities [2,3].
It remains an open problem to determine how much of a diffusivity difference is required for generic systems to undergo this instability, yet this diffusive threshold has been recognized at least since reduced models of the Belousov-Zhabotinsky reaction [4,5] only produced Turing patterns at unphysically large diffusivity differences.
For this reason, the first experimental realizations of Turing instabilities [6][7][8] obviated the threshold by using gel reactors that greatly reduced the effective diffusivity of one species [9,10].(Analogously, biological membrane transport dynamics can increase the effective diffusivity difference [11].)Later work showed how binding to an immobile substrate, or more generally, a third, nondiffusing species, can allow Turing instabilities even if the  = 2 diffusing species have equal diffusivities [12][13][14].Such nondiffusing species continue to permeate more recent work on the network topology of Turing systems [15,16].
Moreover, Turing instabilities need not be deterministic: fluctuation-driven instabilities in reaction-diffusion systems have noise-amplifying properties that allow their pattern amplitude to be comparable to that of deterministic Turing patterns [17], with a lower diffusive threshold than the deterministic one [18][19][20][21].A synthetic bacterial population with  = 2 species that exhibits patterns in agreement with such a stochastic Turing instability, but does not satisfy the conditions for a deterministic instability [22], was reported recently.
These experimental instabilities relying on fluctuations or the dynamics of additional nondiffusing species and the nonlinear instabilities arising from finite-amplitude perturbations [2] are not however instabilities in Turing's own image.Can such instabilities be realized instead in systems with  > 2 diffusing species?Equivalently, is the diffusive threshold lower in such systems?These questions have remained unanswered, perhaps because, in contrast to the textbook case  = 2 and the concomitant picture of an "inhibitor" out-diffusing an "activator" [23,24], the complicated instability conditions for  > 2 [25] do not lend themselves to analytical progress.
Here, we analyze the diffusive threshold for Turing instabilities with 2  6 diffusing species.Inspired by May's work on the stability of random ecological communities [26], we analyze random Turing instabilities by sampling random matrices that represent the linearized reaction dynamics of otherwise unspecified reaction-diffusion systems.A semianalytic approach shows that the diffusive threshold is more likely to be smaller and physical for  = 3 compared to  = 2, and that two of the three diffusivities are equal at the transition to instability.We extend these results to the remaining numerically tractable cases of reaction-diffusion systems with 4  6 and two different diffusivities: their Turing instabilities are still more likely to have a smaller and physical diffusive threshold, but most of them cannot be described by reduced models with fewer species.
We begin with the simplest case,  = 2, in which species ,  obey The conditions for Turing instability in this system [24] only depend on the four entries of the Jacobian the partial derivatives of the reaction system at a fixed point ( * ,  * ) of the homogeneous system.This fixed point is stable to homogeneous perturbations iff  ≡ det J > 0 and  1 ≡ tr J < 0. A stable fixed point of this kind is unstable to a Turing instability only if  ≡ −     > 0 [24].Defining the diffusion coefficient ratio This diffusivity difference  * 2 , which is required mathematically for instability, is unphysical [Fig.1(a)] if it exceeds the diffusivity difference D > 1 of the physical system: D ≈ 1 for similarly sized molecules in solution, but, e.g., D ≈ 20 for the stochastic Turing instability observed in Ref. [22].Hereinafter we take D = 5 arbitrarily (but have checked that the value of D does not affect results qualitatively).To quantify  * 2 , we introduce the range  of kinetic parameters, Equivalently,   ,   ,   ,   ∈  ≡ [−, −1] ∪ [1, ] up to scaling, with one parameter equal to ±1 and one equal to ±.
One deduces [27] that as shown in Fig. 1(a).Note that  max 2 → 1 as  → 1; in this limit, there is no diffusive threshold:  ≈ 1 is a particular instance of the converse fine-tuning problem for the reaction kinetics that allows Turing instabilities at nearly equal diffusivities more generally [3].If  1, then  max 2 =   2 .This does not imply the existence of a threshold, for it does not preclude most systems with range  having  * 2  max 2 .The existence of a diffusive threshold therefore relates to the distribution of  * 2 for systems with range .To understand this distribution, we draw inspiration from May's statistical analysis of the stability of ecological communities [26], which studies random Jacobians, corresponding to equilibria of otherwise unspecified population dynamics.By analogy, we study random Turing instabilities, sampling uniformly and independently random Jacobians corresponding to equilibria of otherwise unspecified reaction kinetics, and analyze the criteria for them to be Turing unstable.There is of course no more reason to expect the kinetic parameters to be independent or uniformly distributed than there is reason to expect the linearized population dynamics in May's analysis [26] to be independent or normally distributed.Yet, in the absence of experimental understanding of what the distributions of these parameters should be (in either context), the potential of the random matrix approach to reveal stability principles has been amply demonstrated in population dynamics [34][35][36][37][38][39][40][41][42][43][44].
We sample the kinetic parameters in Eq. ( 2) independently and uniformly from , set one of them equal to ±1 and one equal to ±, and thus estimate the probability distribution ( * 2 ) for fixed  [Fig.1(b)].The threshold is quantified by the probability of a Turing instability being physical, Both from the estimates in Fig. 1(b) and by evaluating the integral in closed form [27], we find that P( * 2 < D) is tiny [Fig.1(c)], except if  is small, which is the fine-tuning problem again.In other words, the required diffusivity difference is very likely to be unphysical.This expresses Turing's diffusive threshold for  = 2.
To investigate how this threshold changes with  we consider next the  = 3 system where we have rescaled space to set   = 1.We introduce the matrix of diffusivities and the reaction Jacobian, in which the entries of J are again the partial derivatives evaluated at a fixed point ( * ,  * ,  * ) of the homogeneous system.This fixed point is unstable to a Turing instability if it is stable but, for some eigenvalue − 2 < 0 of the Laplacian, J  2 = J −  2 D is unstable [3], i.e. has an eigenvalue  such that Re() < 0.More precisely, a Turing instability arises when a real eigenvalue of J  2 crosses zero, i.e. when J  2 ≡ det J  2 = 0, and therefore arises first at a wavenumber  =  * with J  2 * = J / 2  2 * = 0 [3].Hence J , a cubic polynomial in  2 , has a double root at  2 =  2 * > 0, so its discriminant [29] vanishes.This discriminant, Δ(  ,   ), is a polynomial in   ,   .We denote by  (  ,   ) the double root of J corresponding to a point (  ,   ) on the curve Δ(  ,   ) = 0.
Determining the diffusive threshold for Turing instability in Eqs.(7) thus requires solving the problem in which the diffusion coefficient ratio is With the aim in mind of obtaining statistics for the minimal value  * 3 , direct numerical solution of this constrained optimization problem is obviously not a feasible approach.In the Supplemental Material [27], we therefore show how to reduce solving problem (9) to polynomial root finding.This semianalytic approach reveals a particular class of minima, attained at the vertices of the contours of  3 (  ,   ) [Fig.2(a)], i.e. at   = 1,   = 1, or   =   .In these cases, Δ(  ,   ) = 0 is a (sextic) polynomial in the single variable   ,   , or  =   =   , respectively.We call these minima "binary", since the corresponding systems have only two different diffusivities.We implement this approach numerically [27], and sample random systems similarly to the case  = 2, drawing the entries of J in Eq. ( 8) uniformly and independently at fixed range .
Remarkably, all global minima we found numerically were binary [27].This means that the minimizing systems come in two flavors: those with two "fast" diffusers and one "slow" diffuser, and those with one "fast" diffuser and two "slow" diffusers.Systems with a nondiffusing species are a limit of the former; this point will be discussed below.The latter arise in models of scale pattern formation in fish and lizards [45,46], in which short-range pigments respectively activate and inhibit a long-range factor.
The distribution of  * 3 , shown for different values of  in Fig. 2 is more likely to be physical for  = 3 than for  = 2: the diffusive threshold is lowered.
The proportion   of random kinetic Jacobians that have a Turing instability (be it physical or unphysical) is smaller for  = 3 than for  = 2 [Fig.2(d), inset].This is not surprising, because a random Jacobian is less likely to correspond to a stable fixed point (which, we recall, is a necessary condition for Turing instability) for  = 3 than for  = 2, essentially because its entries have to satisfy more conditions for stability if  = 3.It is therefore striking that the threshold is reduced sufficiently for  = 3 compared to  = 2 for the proportion   (D) =   P( *  < D) of random Jacobians that have a physical Turing instability to be larger for  = 3 than for  = 2 [Fig.2(d)], even though a Turing instability of any kind is more likely if  = 2.
To extend these results to  > 3 diffusing species, we consider the (linearized) reaction-diffusion system where J is a random kinetic Jacobian, and D is a diagonal matrix of diffusivities.Even with our semianalytic approach, this cannot be analyzed for general D: not even for  = 4 were we able to obtain closed forms of the required polynomials.To make further progress, we therefore restrict to binary D in which the  diffusivities take two different values only, since we showed above that  * 3 is attained for such binary D. As in the case  = 3, this reduces the discriminant condition Δ(D) = 0 to polynomial equations in one variable that determine the minimum diffusivity difference  *  for Turing instability in these binary systems [27].
Figure 3(a) shows that the diffusive threshold lowers further for 4  6 in these systems.At the same time, the fact that most stable random kinetic Jacobians undergo such a binary Turing instability [Fig.3(b)] suggests that these provide a useful picture of the diffusive threshold.However,   decreases further for  4 [Fig.3(c), inset], and the widening of the bottleneck is not sufficient to prevent   (D) from decreasing for  4. Nonetheless, since both P( *  < D) and the proportion   of stable random Jacobians that are Turing unstable increase [Figs.3(a) and 3(b)], so does the proportion of stable random Jacobians that have a physical Turing instability.
How then to realize "true" Turing instabilities experimentally?Our analysis shows that the diffusive threshold of a Turing instability is more likely to be physical the more species there are, but how to find an experimental Turing instability in the first place?Turing instabilities remain rare in random reaction systems even as the number of species is increased, but the above shows that this rareness mainly results from the rareness of stable equilibria in such systems.The proverbial search for the needle in a haystack can therefore be avoided by exploring biochemical systems that admit a stable equilibrium, and evolving them towards a "true" Turing instability.
This analysis does not however reveal whether these instabilities lead to patterns that are observable at the physical scale of the system.Analysis of the wavenumber at which the linear instability first arises [27] suggests that we can extend our conclusions: Turing instabilities with more species are more likely to have physical diffusivity differences and to be observable.However, our statistical, linearized analysis cannot fully answer this question of observability, because it fundamentally depends on the system through details of the nonlinearities of its reaction kinetics, which set the precise nature and scale of the Turing patterns that develop beyond onset of the instability; this is why we have relegated this discussion to the Supplemental Material [27].
The different species in the systems with 3  6 analyzed above separate into "fast" and "slow" diffusers.The diffusion of these "slow" species is often ignored in the analysis of systems of many chemical reactions [30], such as the full Belousov-Zhabotinsky reaction [47].Corresponding reduced models are obtained by substituting the steady-state kinetics of the "slow" species into the remaining equations, thereby eliminating them from the system [30].The conditions for Turing instability in these reduced models are (almost) equivalent to those for the full model with nondiffusing "slow" species [30].However, the diffusion of the "slow" species cannot in general be ignored: up to reordering species and rescaling space, where  < 1 is the common diffusivity of the slow diffusers.Results of Ref. [30] imply that there is a Turing instability with nondiffusing "slow" species, i.e. with  = 0, only if J 11 − J 12 J −1 22 J 21 has a positive (real) eigenvalue [27].Although the proportion of Turing unstable systems that have  2 fast diffusers (and hence could a priori still undergo a Turing instability with  = 0) is large [Fig.4(a)], the proportion of systems that do undergo such an instability is small, even if we restrict to those systems with physical diffusivity differences [Fig.4(b)].Hence most of these Turing instabilities with  > 2 species require all species to diffuse.These many-species Turing instabilities, although binary, are thus more general than the instabilities of systems with nondiffusing species realized experimentally in gel reactors [6][7][8] and analyzed theoretically in Ref. [30].In particular, this shows that reduced models give but an incomplete picture of Turing instabilities.Together with our main result, that the diffusive threshold lowers as  increases, this implies that the failure of a reduced model to produce a physical Turing instability cannot be taken as an indication that a Turing instability cannot exist in the full system that the reduced model seeks to describe.
In this Letter, we have analyzed random Turing instabilities to show how the diffusive threshold that has hampered experimental efforts to generate "true" Turing instabilities in systems of  = 2 diffusing species lowers for systems with  3, most of whose instabilities cannot be described by reduced models with fewer species.All of this does not, however, explain the existence of a "large" threshold in the first place: even though Turing instabilities at equal diffusivities are impossible [2,3], this does not mean that the threshold needs to be "large".In this context, we prove an asymptotic result in the Supplemental Material [27]: for a Jacobian J to allow a Turing instability at almost equal diffusivities D ≈ I, J must be even closer to a singular matrix J 0 , i.e.J − J 0 D − I.In this sense, the threshold D − I is asymptotically "large".Understanding how a large threshold arises more generally outside this asymptotic regime and lowers as  increases remains an open problem, as do extending the present analysis to include the nonlocal interactions [48,49] that arise for example in vegetation patterns [50] and extending previous work [16,51] on the robustness of Turing patterns to  3. The latter in particular may help to identify those chemical or biological pattern forming systems with  3 in which the "true" Turing instabilities discussed here can be realized experimentally.

C. Calculation of P P P(𝑫
There are 48 ways of assigning values ±1 and ± to two of the entries   ,   ,   ,   of J. Integrating the conditions for Turing instability of the remaining entries in each of these cases using M (Wolfram, Inc.) gives the area of parameter space in which a Turing instability arises, where we use the shorthand dJ = d   d   d  d  .To analyze the condition  * 2 < , we note that the expression for  * 2 in Eq. ( 3) shows that we may swap   ,   and   ,   .Hence the 48 cases reduce to 4 cases (corresponding to the entries ±1 or ± being on the the same or on different diagonals): Moreover, since  > 0, we may take   > 0 and   < 0 without loss of generality.We now discuss these cases separately.
From a more physical point of view, as discussed in our Letter, it is more natural to consider the probability P  * 2 < D , for some constant D > 1.Since "small" values  D require fine-tuning of the reaction kinetics, we restrict to D , so that  * 2 < D is only possible in case (2) above.We consider again the case   > 0,   < 0. Similarly to the derivation of conditions (S9), we find In particular, in which, since   1, We notice that The area of parameter space described by conditions (S12) is therefore for  > D, and, as above, we conclude that, for  > D,

D. Nondimensionalization
We close by remarking on the (absence of) nondimensionalization of the reaction system.Indeed, up to rescaling time, one among   ,   ,   ,   can be set equal to ±1.Moreover, one more parameter can be set equal to ±1 by rescaling ,  differently.However, if we made those choices, we could no longer sample from a fixed interval.

II. SEMIANALYTIC METHOD FOR 𝑵 = 3
A. Derivation of the semianalytic method

Preliminary observations
Before deriving the semianalytic method, we need to make two preliminary observations.First, the necessary and sufficient (Routh-Hurwitz) conditions for the homogeneous system to be stable include  1 ≡ tr J < 0 and  ≡ det J < 0 [S1].By definition, J  2 * has one zero eigenvalue.The other two eigenvalues are either real or two complex conjugates ,  * .In the second case, they are both stable (i.e. have negative real parts) since Hence Eqs.(7) are not unstable to an oscillatory (Turing-Hopf) instability at ( *  ,  *  ), so, by minimality of ( *  ,  *  ), the system destabilizes to a Turing instability there.
Moreover, since J , viewed as a polynomial in  2 * , has leading coefficient −    and constant term J (0) =  < 0, the double root  (  ,   ) varies continuously with   ,   and cannot change sign on a branch of Δ(  ,   ) = 0 in the positive (  ,   ) quadrant.

Reduction of problem (9) to polynomial equations
The discriminant of J , viewed as a polynomial in the two variables   ,   , is where  00 =  10 =  01 =  34 =  43 =  44 = 0 and (complicated) expressions for the 19 non-zero coefficients can be found in terms of the entries of J using M (Wolfram, Inc.).
The contours of  3 (  ,   ) are drawn in Fig. 2(a) of our Letter and show that tangency to a contour in case (i) requires d  = 0 or d  = 0 or d  /d  =   /  .(S18) Since Δ(  ,   ) = 0, the chain rule reads Hence there are two subcases: In subcase (a), Δ viewed as a polynomial in   or   has a double root, and so its discriminant [S3] must vanish.On removing zero roots, this discriminant of a discriminant is found to be a polynomial of degree 20 in   or   , respectively; complicated expressions for its coefficients in terms of the non-zero coefficients   in Eq. (S17) are obtained using M .Similarly, in subcase (b), the resultant [S3] of Δ and   Δ/  +   Δ/  , viewed as polynomials in   or   must vanish.This resultant is another polynomial of degree 20 in   or   .
Thus, we have reduced finding candidates for local minima in (9) to solving polynomial equations: this defines our semianalytic approach.The global minimum is found among those local minima with  (  ,   ) > 0; in case (i), the roots only correspond to local minima if additionally   ,   > 1 or   ,   < 1 in subcase (a) and   < 1 <   or   < 1 <   in subcase (b) [Fig.2(a)].

Extension to binary systems with 𝑁 > 3
For binary systems, the diagonal entries of D take two different values,  1 ,  2 only.Up to rescaling space,  1 = 1 and  2 = , which turns the condition Δ(D) = 0 into 2  −1 − 1 different polynomial equations in the single variable , corresponding to the different combinatorial ways of assigning diffusivities  1 ,  2 to the  species (in such a way that not all species have the same diffusivity).Determining the minimum value  *  of   = max {, 1/} for these binary systems is thus reduced, again, to solving polynomial equations.
The argument we used above to show that coexistence of Turing and Turing-Hopf instabilities is not possible for  = 3 does not, however, carry over to  > 3. Numerically, it turns out, however, that systems in which Turing and Turing-Hopf instabilities coexist are rare.We therefore treat these systems in the same way as we treat systems for which the numerics fail (as discussed below).

B. Numerical implementation
Implementing the semi-analytical approach for  = 3 and its extension to binary systems with 4  6 numerically takes some care as the coefficients of the polynomials that arise can range over many orders of magnitude.Our python3 implementation therefore uses the mpmath library for variable precision arithmetic [S4].
To determine the positive real roots of the polynomials that arise in the semi-analytical approach, we complement the Durand-Kerner complex root finding implemented in the mpmath library [S4] with a test based on Sturm's theorem [S3], to ensure that all positive real roots are found.Those systems in which root finding fails-either because the Durand-Kerner algorithm fails to converge or because it finds an incorrect number of positive real roots-are discarded, but included in error estimates where reported.
For  = 3, we ran both a search for general, non-binary systems and a (larger but numerically less expensive) search for binary systems only.Since the first search only yielded binary global minima (as stated in our Letter), we used the results of the second, larger search for Figs. 3 and 4.

III. WAVENUMBER STATISTICS
In this Section, we discuss the wavenumber  *  at which a Turing instability first arises at   =  *  .In particular, as discussed in our Letter, we must ask whether a Turing instability is "observable at the system size".This observability requires the lengthscale 1/ *  of the linear instability to be (a) smaller than the system  and (b) larger than /ℓ, for some scale difference ℓ > 1.We are thus led to consider the probability P( <  *  < ℓ), where  = 1/.It is instructive to start by considering the case  = 2.For the reaction-diffusion system in Eq. ( 1), a Turing instability arises for  2 =  * 2 at a wavenumber  * 2 = (/    ) 1/4 [S1].We stress that this value depends on   ,   not only through their ratio  =   /  .To absorb the dependence on the dimensional system scale, it is natural to consider as the maximal probability of a Turing instability being observable at some inverse system scale  over a fixed scale difference ℓ.We denote by  2 (ℓ) the corresponding maximizing inverse system size.
For  > 2, we correspondingly ask: what is the probability of a Turing instability being observable at this inverse system size?We therefore define Figure S1 plots   (ℓ) against , for fixed values of  and ℓ, but the qualitative behaviour is independent of  and ℓ.We notice that   (ℓ) increases slightly with .If we restrict the analysis to those Turing unstable systems with  *  D, the probability is reduced somewhat for  > 2 compared to the case  = 2.This merely reflects the "fine-tuning problem": the wavenumber is strongly constrained for those very rare systems that have a "small" diffusive threshold at  = 2.Moreover, about three quarters of the Turing instabilities at  > 2 do arise at physical wavenumbers, so we can extend the observations in Figs.2(d) and 3(c) to note that random kinetic Jacobians are still more likely to be unstable to an observable Turing instability with small diffusive threshold for  > 2 than for  = 2.

IV. DIFFUSION OF "SLOW" SPECIES
In the notation of Eq. ( 12) of our Letter, Ref. [S5] shows that Turing instability at  = 0 requires J 22 to be stable (i.e.all its eigenvalues to have negative real part): if it is not, instabilities arise at arbitrarily small and therefore unphysical lengthscales.In particular, det J 22 ≠ 0, and so, using another result of Ref. [S5], where j = J 11 − J 12 J −1 22 J 21 .Hence a Turing instability occurs at  = 0 only if j has a positive real eigenvalue, as claimed in our Letter.and  2 = (1) since Re() < 0. Hence J and J −  2 I have a zero eigenvalue at leading order.Additionally, the eigenvalue correction from − 2 d =   2 must be   2 at least, which occurs iff the (leading-order) zero eigenspaces of J −  2 I and J are defective [S7]; this final implication is discussed in more detail in Ref. [S8].

V. THE ASYMPTOTIC DIFFUSIVE THRESHOLD
The generic case is therefore J = J 0 +  (), where  1 and J 0 has a defective double zero eigenvalue.

FIG. 3 .
FIG. 3.Results for "binary" systems with 4  6.(a) P( *  < D) against  for 3  6, revealing further lowering of the diffusive threshold compared to the case  = 3.(b) Proportion   of random stable kinetic Jacobian that have a (binary, if  > 3) Turing instability, averaged over , and plotted against .(c) Proportion   (D) of random Jacobians that have a physical Turing instability plotted against , for 3  6. Inset: proportion   of random Jacobians that have a (physical or unphysical) Turing instability, averaged over , for 3  6[33].

FIG. S1 .
FIG. S1.Wavenumber statistics.Probability   (ℓ) of a Turing instability being "observable" at a scale difference ℓ plotted against ; see text for further explanation.Larger markers:   (ℓ) estimated from all Turing unstable systems; smaller markers:   (ℓ) estimated from only those Turing unstable systems with  *  < D. Parameter values:  = 10, ℓ = 10, D = 5.Asymmetric error bars again correspond to 95% confidence intervals larger than the plot markers, corrected for systems for which the numerics failed.

TABLE S1 .
Number of random Turing unstable systems used to estimate distributions, averages, and probabilities for the different values of , and corresponding figures.
a Maximal number of Turing unstable systems.b Figure (if any) in which results are shown.