Turing Instabilities are Not Enough to Ensure Pattern Formation

Symmetry-breaking instabilities play an important role in understanding the mechanisms underlying the diversity of patterns observed in nature, such as in Turing’s reaction–diffusion theory, which connects cellular signalling and transport with the development of growth and form. Extensive literature focuses on the linear stability analysis of homogeneous equilibria in these systems, culminating in a set of conditions for transport-driven instabilities that are commonly presumed to initiate self-organisation. We demonstrate that a selection of simple, canonical transport models with only mild multistable non-linearities can satisfy the Turing instability conditions while also robustly exhibiting only transient patterns. Hence, a Turing-like instability is insufficient for the existence of a patterned state. While it is known that linear theory can fail to predict the formation of patterns, we demonstrate that such failures can appear robustly in systems with multiple stable homogeneous equilibria. Given that biological systems such as gene regulatory networks and spatially distributed ecosystems often exhibit a high degree of multistability and nonlinearity, this raises important questions of how to analyse prospective mechanisms for self-organisation.


Introduction
Nature exhibits diverse structures in the organisation of life across spatial and temporal scales.Elaborate animal coat patterns (Koch and Meinhardt, 1994), emergent territory boundaries between predators (Potts and Lewis, 2016), and complex spatiotemporal arrangements of slime moulds (H öfer et al., 1995) are a few of the patterns researchers have sought to understand.A key mechanism underlying such patterns are symmetry-breaking (Turing) instabilities of spatially uniform equilibria, as explored in Turing's influential Chemical basis of morphogenesis (Turing, 1952).
Typical analysis of these phenomena is often based on linear stability theory, which attempts to ascertain the growth or decay of perturbations to homogeneous equilibria.Due to the nature of the resulting linear equations, such analysis can often be carried out very easily.In addition to its simplicity, a chief advantage to this approach is its generality, as it makes minimal assumptions about the precise form of the underlying system.In turn, this provides reasonably broad statements about the kinds of systems that can exhibit such instabilities, as illustrated by the fact that Turing self-organisation in a two-species reaction diffusion system requires a short-range (self)-activator, and a long-range (self)inhibitor (Meinhardt and Gierer, 2000).The simplicity of linear stability analysis means that, even for many-species systems (Marcon et al., 2016), one can typically classify parameters of the linearised system into those that exhibit pattern-forming instabilities (the so-called 'Turing space'), and those which cannot (Murray, 1982(Murray, , 2003)).There is a large body of work aimed at understanding features of these Turing spaces in various contexts (Klika and Gaffney, 2017;Marcon et al., 2016;Gaffney et al., 2023), but always using some form of linear stability theory, which is a dominant feature of the pattern formation literature.Hence a large number of studies have focused on linear systems exclusively to make general claims about proposed mechanisms (Satnoianu et al., 2000;Krause et al., 2020;Haas and Goldstein, 2021) or to design pattern-forming systems with certain properties (Vittadello et al., 2021;Woolley et al., 2021).
However, when linear analysis identifies a pattern-forming instability, the output is always a local result, in that transient symmetry-breaking patterns are expected to form from perturbations of the homogeneous steady state.Beyond the formation of an initial pattern, linear stability provides no guarantee of a long-time (i.e.stable) patterned state.Notably, the existence of stable patterns can be guaranteed in the case of supercritical Turing bifurcations, but only near the boundary of the Turing space (Vastano et al., 1988).However, the emergence of large-scale, persistent self-organisation is invariably presumed from the linear analysis (including by the authors), often based on intuition and experience with simple examples of minimal complexity (Murray, 2003;Krause et al., 2021).
While this intuition has been seen to be correct for many textbook systems, extensive recent examples highlight that linear stability theory cannot always capture the fundamental dynamics of patternforming systems, such as instabilities due to subcritical bifurcations (Champneys et al., 2021;Villar-Sep úlveda and Champneys, 2023).Unlike in the supercritical case, subcritical bifurcations do not typically admit small-amplitude stable patterned states, even in the weakly nonlinear regime except very near to the codimension-2 point where the criticality of the bifurcation changes (Bre ña-Medina and Champneys, 2014).Such subcritical bifurcations can lead to pattern formation outside of Turing space, as implicated in ecological work on resilience due to patterning (van de Koppel and Rietkerk, 2004;Bastiaansen et al., 2020), among other areas.Subcritical bifurcations can also lead to spatiotemporal oscillations and chaos (Painter and Hillen, 2011).Other secondary bifurcations can eliminate any stable patterned branches, so that systems with multiple spatial homogeneous equilibria may form only transient patterned states; see Figures 8 and 11 in Al-Karkhi et al. (2020) for an example.Non-normality (in the sense of normal matrices/operators) can also lead to different predictions from linear theory, as described by Klika (2017).Thus, classical linear stability conditions are neither necessary nor sufficient for self-organisation.
Here, we demonstrate that this insufficiency of the classical Turing conditions can occur generically in a range of systems.In particular, we exemplify that the presence of multistability can robustly spoil typical predictions of patterning by driving a system to a stable homogeneous equilibrium after the emergence of transient patterns via a Turing instability.Multistability has become an increasingly prominent topic in gene regulatory networks (Laurent and Kellershohn, 1999;Siegal-Gaskins et al., 2009;Feng et al., 2016;Bocci et al., 2023), ecology (Suzuki et al., 2021), and evolutionary biology (Arnoldt et al., 2012), with growing evidence that multistable dynamics are ubiquitous in biological systems.Here, we show that even bistability of reaction kinetics can alter the prospect for pattern formation in a robust way, suggesting a need for better tools to analyze more realistic models of pattern formation in biological systems.
The rest of the paper is organized as follows.In Section 2, we present and perform a linear stability analysis of specific models from four distinct classes of pattern-forming systems.In Section 3, we perform thousands of numerical experiments with random parameters and demonstrate that these models void our long-established intuition for pattern formation relying on linear stability theory, raising important issues regarding the connection between textbook analyses and realistic biological systems, which we discuss in Section 4.

Models & dispersion relations
We consider four models on the spatial domains with periodic boundary conditions.Parameters are assumed to take positive nonzero values, with exemplars given in Table 1 for each model, which we will refer to as the base parameters.
For each model, we perform a linear stability analysis around one of the spatially homogeneous equilibria and record the growth rate of spatial perturbations corresponding to the eigenvalues ρ k of the negative Laplacian given by with periodic boundary conditions, ordered via 0 = ρ 0 < ρ 1 ≤ ρ 2 ≤ . . ., with w k the corresponding eigenfunctions (these are just sinusoidal functions for these domains and boundary conditions).We then write the maximal growth rate of linear perturbations corresponding to eigenfunction w k as λ k , so that ℜ(λ 0 ) < 0 and ℜ(λ k ) > 0 for some k > 0 is our criterion for a Turing instability.Analysing the first three models is standard (Murray, 2003;Krause et al., 2021), whereas the linear stability theory for the nonlocal advection model is given by Jewell et al. (2023).In each case, we will focus on the linear stability of one equilibrium, but each model will also admit one other stable equilibrium.

Reaction-diffusion system
We first consider a two-component reaction-diffusion system of the form which has a homogeneous equilibrium at (u 0 , v 0 ) = (0, 0).For the parameters in Table 1, this equilibrium is stable in the absence of diffusion.There are four further real equilibria, only one of which is stable in the absence of diffusion.Linearising Eq. ( 2) around (u 0 , v 0 ) = (0, 0) gives perturbation growth rates For the base parameters, the equilibrium (u 0 , v 0 ) = (0, 0) is Turing unstable (see the plot of the dispersion relation in Fig. 1(a)).

Biharmonic instability
Next, we consider a fourth-order model of self-organisation: The spatially homogeneous equilibria are u 0 = 0, b, c, which exhibit bistability of u 0 = 0 and u 0 = c in the absence of transport for c > b > 0. Linearising Eq. ( 6) around u 0 = c gives perturbation growth rates For the base parameters, the equilibrium u 0 = c is Turing unstable (see the plot of the dispersion relation in Fig. 1(c)).

Nonlocal advection
Finally, we consider an integro-differential model of cell aggregation (Painter et al., 2015;Jewell et al., 2023): where ds N is the volume element for N = 1 or N = 2 spatial dimensions.We require c > b > 0 for stability of the spatially homogeneous equilibria u 0 = 0, c in the absence of transport, while u 0 = b is unstable.Linearising Eq. ( 8) around the u 0 = c equilibrium gives perturbation growth rates In 1D and 2D with the base parameters, the equilibrium u 0 = c is Turing unstable (see the plot of the 1D dispersion relation in Fig. 1(d)).

Results
Each of these systems admits a Turing instability for the parameters given by Table 1 for one of their equilibria, illustrated in the dispersion plots of Fig. 1(a)-(d).Hence, following commonplace reasoning, one might presume that a pattern (a stationary or spatiotemporal solution bounded away from homogeneous solutions) will form from perturbations of these equilibria.However, numerical simulations of  2023)).The dynamics of (a)-(c) can be explored interactively using VisualPDE (Walker et al., 2023) at https://visualpde.com/ mathematical-biology/Turing-conditions-are-not-enough.
these models in 1D in Fig. 1(e)-(h) and in 2D in Fig. 2 show transient pattern formation that then decays to a different homogeneous equilibrium, all of which are linearly stable.Briefly, the three local models are solved using finite differences, and the nonlocal model using a pseudospectral method combined with finite differences.In all cases, implicit timestepping algorithms are used, with initial data given by normally distributed perturbations of the Turing unstable equilibrium (of standard deviation 10 −2 ), as detailed in the repository Krause et al. (2023).
Importantly, this decay to homogeneity occurs robustly across variation in all parameters.To demonstrate this, we vary parameters and initial conditions as follows.For each model, we multiply every parameter given in Table 1, including the domain length L, by a uniformly random number from the interval [0.95, 1.05] using Latin Hypercube Sampling (Wyss and Jorgensen, 1998).We then simulate the system for t = 10 4 time units from a different random initial perturbation, recording both if there is a Turing instability (by analysing the dispersion relation) and if the system is approaching a homogeneous state (assessed by checking if max x (u(10 4 , x)) − min x (u(10 4 , x)) > 10 −5 ).From this, we can determine the proportion of simulations that exhibited a Turing instability but only patterned transiently from a small random perturbation of the homogeneous equilibrium.We perform 10 4 simulations for all 1D models, and 10 3 simulations for the 2D models.We omit the 2D nonlocal advection system from this analysis due to its numerical complexity.
We give results of this analysis in Table 2.For all but the reaction-diffusion model, an overwhelming majority of cases converged to homogeneous equilibria after transiently patterning via a Turing instability.All of these systems remained within the Turing space of their corresponding equilibria.The reaction-diffusion model also exhibited robust convergence to a homogeneous equilibrium, though as the base parameters of the reaction-diffusion model lie near the boundary of the Turing space, not all simulations were Turing-unstable.A small proportion of the reaction-diffusion and Biharmonic systems were attracted to patterned equilibria, some of which were domain-filling while others appeared spatially localised (Champneys et al., 2021).Both the Keller-Segel and 1D nonlocal advection models only exhibited convergence to a homogeneous equilibrium after transient patterning.Sufficiently changing the parameters of the models can give rise to other behaviours.Rather than detail these observations, we encourage the reader to interactively explore the three local models with VisualPDE (Walker et al., 2023) 1 .For instance, by changing properties of the bistability, one can observe transitions between systems that form no patterns, favour localised solutions, or admit domain-filling patterns.Indeed, exploring the Keller-Segel equation interactively via VisualPDE is how we first observed this behaviour, with the other three models designed to mimic the basic ingredients of bistability and subcriticality.

Discussion
Across a range of models, parameter sets, and different initial conditions, we have robustly observed that possessing a Turing instability is not sufficient for systems to form spatial patterns that persist beyond transient timescales (the timescales observed in the examples in Fig. 1 and Fig. 2 are plausibly too short to be compatible with many examples of biological patterning, though this would depend on the details of nondimensionalisation).Conversely, wave-pinning and other mechanisms can give rise to spatially structured stable states without a Turing-like bifurcation (Champneys et al., 2021).Therefore, while linear theory can have value in detecting self-organisation, it is perhaps not as generally valid as most of the textbook examples (e.g.every reaction-diffusion system in the book (Murray, 2003)) might indicate.
We suspect that the almost ubiquitous association between Turing instabilities and pattern formation is largely because most research on patterning in reaction-transport systems, including our own (Krause et al., 2021), focuses on small systems of at most three or four components with relatively mild nonlinearities.Systems such as Eq. ( 2) are still in this class of relatively simple systems, but the presence of bistability might be more indicative of large and complex reaction networks, which likely exhibit a high degree of multistability.Systematic analyses of such systems are relatively unexplored, and the results we have shown underscore the importance of studying them.Additionally, emphasis on supercritical bifurcations with stable small-amplitude patterns near the bifurcation point can fail to capture both systems exhibiting subcritical bifurcations as well as systems far away from the original Turing bifurcation point.
Among modern tools for far-from-equilibrium analyses, we note that there exist several approaches for studying spike or pulse dynamics (Wei and Winter, 2013;Doelman and Veerman, 2015).Such approaches have shown the importance of even small changes to the nonlinearity on the existence and stability of patterned states (Veerman and Doelman, 2013).Contemporary numerical continuation techniques, such as in the PDE2PATH software (Uecker et al., 2014) can be used to describe the loss of patterned states we explored here, as shown in Figure 11 of Al-Karkhi et al. (2020).These approaches, however, typically focus on studying specific models and parameter sets, and do not lend themselves as easily to studying generic systems, especially those with more than two components.In contrast, linear stability theory has been employed to classify larger reaction-diffusion systems (Marcon et al., 2016;Scholes et al., 2019;Landge et al., 2020).Recent approaches such as Local Perturbation Analysis (Holmes, 2014;Holmes et al., 2015) overcome some of these limitations, at the cost of only strictly applying in particular asymptotic regimes.An important avenue would be the development of more powerful tools to understand complex and nonlinear systems in the context of pattern formation without relying on the 1 https://visualpde.com/mathematical-biology/Turing-conditions-are-not-enoughlimitations of looking only locally in the phase space or in the space of parameters/models.We view this as an exciting and important frontier for future theoretical work.
For the base parameters, the equilibrium (u 0 , v 0 ) = (b, b/a) is Turing unstable (see the plot of the dispersion relation in Fig. 1(b)).
Fig. 1: (a)-(d) Plots of λ k against the continuous spatial eigenvalue ρ k in 1D, with orange dots corresponding to discrete values ρ k from the finite domains of size L given in Table 1, with the equilibrium being perturbed given in the text.(e)-(h) Kymographs of u over time in each model following perturbations from their Turing-unstable equilibria.Columns correspond to the models of Section 2.

Fig. 2 :
Fig. 2: Snapshots of transient 2D dynamics with Turing instabilities.(a)-(d) The evolution of u for the models in Section 2 from initial perturbations.Simulations continue to decay towards homogeneous equilibria for times beyond those shown (see the Videos folder at Krause et al. (2023)).The dynamics of (a)-(c) can be explored interactively using VisualPDE(Walker et al., 2023) at https://visualpde.com/ mathematical-biology/Turing-conditions-are-not-enough.

Table 1 :
Base model parameters for the four different models.

Table 2 :
Column 2: Proportion of simulations unpatterned at the final time t = 10 4 .Column 3: Proportion of simulations that were Turing unstable at the initial homogeneous equilibrium.Column 4: Proportion of the Turing-unstable simulations from Column 3 that decayed to a different homogeneous equilibrium.