4pt 4pt

11institutetext:  A. L. Krause 22institutetext: Department of Mathematical
Sciences, Durham University, Upper Mountjoy Campus, Stockton Road, Durham DH1
3LE, United Kingdom
E. A. Gaffney, T. J. Jewell 33institutetext: Wolfson Centre for Mathematical
Biology, Mathematical Institute, University of Oxford, Oxford, OX2 6GG, United
Kingdom
V. Klika 44institutetext: Department of Mathematics, FNSPE, Czech Technical
University in Prague, Trojanova 13, 120 00 Praha, Czech Republic
B. J. Walker 55institutetext: Department of Mathematical Sciences, University
of Bath, Bath BA2 7AY, United Kingdom; Department of Mathematics, University
College London, London WC1E 6BT, United Kingdom

# Turing instabilities are not enough to ensure pattern formation

Andrew L. Krause111Corresponding author<EMAIL_ADDRESS>Eamonn A.
Gaffney  Thomas Jun Jewell  Václav Klika  Benjamin J. Walker

(Received: date / Accepted: date)

###### Abstract

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.

###### Keywords:

Turing instabilities; pattern formation

##  1 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, 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 pattern-forming 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.

##  2 Models & dispersion relations

We consider four models on the spatial domains $\Omega=[0,L]$ or
$\Omega=[0,L]\times[0,L]$, 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 $\rho_{k}$ of the negative
Laplacian given by

$\displaystyle\nabla^{2}w_{k}(\bm{x})=-\rho_{k}w_{k}(\bm{x})$ (1)

with periodic boundary conditions, ordered via
$0=\rho_{0}<\rho_{1}\leq\rho_{2}\leq\dots$, 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 $\lambda_{k}$, so that
$\Re(\lambda_{0})<0$ and $\Re(\lambda_{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.

Model | $L$ | $D$ | $a$ | $b$ | $c$ | $d$ | $e$
---|---|---|---|---|---|---|---
Reaction–diffusion | 100 | $30$ | $1.75$ | $18$ | $2$ | $5$ | $0.02$
Keller–Segel | $80$ | $1$ | $1$ | $1$ | $3$ | $0.8$ | -
Biharmonic | $100$ | $1.45$ | $5$ | $0.9$ | $1$ | - | -
Nonlocal Advection | $30$ | $1$ | $1$ | $0.45$ | $0.5$ | $20$ | -
Table 1: Base model parameters for the four different models.

###  2.1 Reaction–diffusion system

We first consider a two-component reaction–diffusion system of the form

$\displaystyle\frac{\partial u}{\partial t}=\nabla^{2}u+u-v-eu^{3},$ (2)
$\displaystyle\frac{\partial v}{\partial t}=D\nabla^{2}v+av(v+c)(v-d)+bu-
ev^{3},$

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 LABEL:RD_model around $(u_{0},v_{0})=(0,0)$
gives perturbation growth rates

$\displaystyle\lambda_{k}=\frac{1-acd-\rho_{k}(1+D)+\sqrt{(1-acd-\rho_{k}(1+D)))^{2}-4(\rho_{k}^{2}D-\rho_{k}(D-acd)-acd+b))}}{2}.$
(3)

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)).

###  2.2 Keller–Segel with Allee demographics

We next consider a Keller–Segel (Keller and Segel, 1970; Horstmann, 2003)
model of chemotaxis:

$\displaystyle\frac{\partial u}{\partial t}$
$\displaystyle=\nabla^{2}u-c\bm{\nabla}\cdot(u\bm{\nabla}v)+u(b-u)(u-d),$ (4)
$\displaystyle\frac{\partial v}{\partial t}$
$\displaystyle=D\nabla^{2}v+u-av.$

The system admits three spatially homogeneous equilibria, $v_{0}=u_{0}/a$ with
$u_{0}=0,d,b$. This is bistable in the absence of transport if $b>d>0$, with
stable equilibria $u_{0}=0$ and $u_{0}=b$. Linearising Eq. 4 around
$(u_{0},v_{0})=(b,b/a)$ gives perturbation growth rates

$\displaystyle\lambda_{k}=\frac{T_{\textrm{KS}}+\sqrt{T_{\textrm{KS}}^{2}-4Q_{KS}}}{2},$
(5)

where $T_{\textrm{KS}}=-b(b-d)-a-\rho_{k}(1+D)$, and
$Q_{\textrm{KS}}=\rho_{k}^{2}D-\rho_{k}(cb-a-Db(b-d))+b(b-d)a$. 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)).

###  2.3 Biharmonic instability

Next, we consider a fourth-order model of self-organisation:

$\displaystyle\frac{\partial u}{\partial
t}=-D\nabla^{2}u-\nabla^{4}u+au(c-u)(u-b).$ (6)

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

$\displaystyle\lambda_{k}=D\rho_{k}-\rho_{k}^{2}+ac(b-c).$ (7)

For the base parameters, the equilibrium $u_{0}=c$ is Turing unstable (see the
plot of the dispersion relation in Fig. 1(c)).

\begin{overpic}[permil,width=433.62pt]{figs/1D.png} \put(0.0,320.0){(a)}
\put(255.0,320.0){(b)} \put(505.0,320.0){(c)} \put(760.0,320.0){(d)}
\put(0.0,160.0){(e)} \put(255.0,160.0){(f)} \put(505.0,160.0){(g)}
\put(760.0,160.0){(h)} \put(58.0,330.0){{\footnotesize Reaction--diffusion}}
\put(330.0,332.0){{\footnotesize Keller--Segel}}
\put(590.0,330.0){{\footnotesize Biharmonic}} \put(823.0,330.0){{\footnotesize
Nonlocal advection}} \end{overpic} Figure 1: (a)-(d) Plots of $\lambda_{k}$
against the continuous spatial eigenvalue $\rho_{k}$ in 1D, with orange dots
corresponding to discrete values $\rho_{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.

###  2.4 Nonlocal advection

Finally, we consider an integro-differential model of cell aggregation
(Painter et al., 2015; Jewell et al., 2023):

$\displaystyle\frac{\partial u}{\partial t}$
$\displaystyle=D\nabla^{2}u+au(c-u)(u-b)$ (8)
$\displaystyle-d\bm{\nabla}\cdot\left(u(1-u)\int_{\Omega}\frac{\bm{s}}{\left\lVert\bm{s}\right\rVert}\frac{\mathrm{e}^{-\left\lVert\bm{s}\right\rVert}}{2\pi}u(\bm{x}+\bm{s})d\bm{s}^{N}\right),$

where $d\bm{s}^{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

$\displaystyle\lambda_{k}=-ac(c-b)-D\rho_{k}+\frac{c(1-c)d\rho_{k}}{\pi^{2-N}(1+\rho_{k})^{\frac{N+1}{2}}}.$
(9)

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)).

##  3 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 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_{\bm{x}}(u(10^{4},\bm{x}))-\min_{\bm{x}}(u(10^{4},\bm{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.

\begin{overpic}[permil,width=303.53267pt]{figs/snapshots.png}
\put(-30.0,865.0){(a)} \put(-30.0,645.0){(b)} \put(-30.0,425.0){(c)}
\put(-30.0,205.0){(d)} \put(400.0,-12.0){Time $\longrightarrow$}
\put(-212.0,769.0){{\footnotesize Reaction--diffusion}}
\put(-140.0,549.0){{\footnotesize Keller--Segel}}
\put(-130.0,329.0){{\footnotesize Biharmonic}}
\put(-217.0,109.0){{\footnotesize Nonlocal advection}} \end{overpic} Figure 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.

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.

Model | Final state unpatterned | Turing unstable | Turing conditions insufficient
---|---|---|---
Reaction–diffusion | 1D |  $93.01$% |  $49.76$% |  $86.17$%
2D |  $99.0$% |  $46.5$% |  $98.75$%
Keller–Segel | 1D |  $100$% |  $100$% |  $100$%
2D |  $100$% |  $100$% |  $100$%
Biharmonic | 1D |  $99.43$% |  $100$% |  $99.43$%
2D |  $100$% |  $100$% |  $100$%
Nonlocal advection | 1D |  $100$% |  $100$% |  $100$%
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.

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)222https://visualpde.com/mathematical-biology/Turing-conditions-are-not-
enough. 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.

##  4 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 LABEL:RD_model 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 limitations 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.

###### Acknowledgements.

BJW is supported by the Royal Commission for the Exhibition of 1851.
There is no data presented in this paper. All code associated with the project
can be found on GitHub (Krause et al., 2023). Interactive versions of the
local models can be found at the website https://visualpde.com/mathematical-
biology/Turing-conditions-are-not-enough.

## References

  * Al-Karkhi et al. [2020] T. A. Al-Karkhi, R. Kusdiantara, H. Susanto, and E. A. Codling.  Bloom formation and turing patterns in an infochemical mediated multitrophic plankton model.  _International Journal of Bifurcation and Chaos_ , 30(10):2030028, 2020.
  * Arnoldt et al. [2012] H. Arnoldt, M. Timme, and S. Grosskinsky.  Frequency-dependent fitness induces multistability in coevolutionary dynamics.  _Journal of The Royal Society Interface_ , 9(77):3387–3396, 2012.
  * Bastiaansen et al. [2020] R. Bastiaansen, A. Doelman, M. B. Eppinga, and M. Rietkerk.  The effect of climate change on the resilience of ecosystems with adaptive spatial pattern formation.  _Ecology letters_ , 23(3):414–429, 2020.
  * Bocci et al. [2023] F. Bocci, D. Jia, Q. Nie, M. K. Jolly, and J. Onuchic.  Theoretical and computational tools to model multistable gene regulatory networks.  _Reports on Progress in Physics_ , 2023.
  * Breña-Medina and Champneys [2014] V. Breña-Medina and A. Champneys.  Subcritical turing bifurcation and the morphogenesis of localized patterns.  _Physical Review E_ , 90(3):032923, 2014.
  * Champneys et al. [2021] A. R. Champneys, F. Al Saadi, V. F. Breña-Medina, V. A. Grieneisen, A. F. Marée, N. Verschueren, and B. Wuyts.  Bistability, wave pinning and localisation in natural reaction–diffusion systems.  _Physica D: Nonlinear Phenomena_ , 416:132735, 2021.
  * Doelman and Veerman [2015] A. Doelman and F. Veerman.  An explicit theory for pulses in two component, singularly perturbed, reaction–diffusion equations.  _Journal of Dynamics and Differential Equations_ , 27:555–595, 2015.
  * Feng et al. [2016] S. Feng, M. Sáez, C. Wiuf, E. Feliu, and O. S. Soyer.  Core signalling motif displaying multistability through multi-state enzymes.  _Journal of the Royal Society Interface_ , 13(123):20160524, 2016.
  * Gaffney et al. [2023] E. A. Gaffney, A. L. Krause, P. K. Maini, and C. Wang.  Spatial heterogeneity localizes Turing patterns in reaction-cross-diffusion systems.  _Discrete and Continuous Dynamical Systems-B_ , pages 0–0, 2023.
  * Haas and Goldstein [2021] P. A. Haas and R. E. Goldstein.  Turing’s diffusive threshold in random reaction-diffusion systems.  _Physical Review Letters_ , 126(23):238101, 2021.
  * Höfer et al. [1995] T. Höfer, J. A. Sherratt, and P. K. Maini.  Cellular pattern formation during dictyostelium aggregation.  _Physica D: Nonlinear Phenomena_ , 85(3):425–444, 1995.
  * Holmes [2014] W. R. Holmes.  An efficient, nonlinear stability analysis for detecting pattern formation in reaction diffusion systems.  _Bulletin of mathematical biology_ , 76:157–183, 2014.
  * Holmes et al. [2015] W. R. Holmes, M. A. Mata, and L. Edelstein-Keshet.  Local perturbation analysis: a computational tool for biophysical reaction-diffusion models.  _Biophysical journal_ , 108(2):230–236, 2015.
  * Horstmann [2003] D. Horstmann.  From 1970 until present: the Keller-Segel model in chemotaxis and its consequences i.  _Jahresbericht der Deutschen Mathematiker-Vereinigung_ , 105:103–165, 2003.
  * Jewell et al. [2023] T. J. Jewell, A. L. Krause, P. K. Maini, and E. A. Gaffney.  Patterning of nonlocal transport models in biology: the impact of spatial dimension.  _arXiv:2307.03117_ , 2023.
  * Keller and Segel [1970] E. F. Keller and L. A. Segel.  Initiation of slime mold aggregation viewed as an instability.  _Journal of Theoretical Biology_ , 26(3):399–415, 1970.
  * Klika [2017] V. Klika.  Significance of non-normality-induced patterns: Transient growth versus asymptotic stability.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 27(7), 2017.
  * Klika and Gaffney [2017] V. Klika and E. A. Gaffney.  History dependence and the continuum approximation breakdown: the impact of domain growth on Turing’s instability.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 473(2199):20160744, 2017.
  * Koch and Meinhardt [1994] A.-J. Koch and H. Meinhardt.  Biological pattern formation: from basic mechanisms to complex structures.  _Reviews of Modern Physics_ , 66(4):1481, 1994.
  * Krause et al. [2020] A. L. Krause, V. Klika, T. E. Woolley, and E. A. Gaffney.  From one pattern into another: analysis of turing patterns in heterogeneous domains via WKBJ.  _Journal of the Royal Society Interface_ , 17(162):20190621, 2020.
  * Krause et al. [2021] A. L. Krause, E. A. Gaffney, P. K. Maini, and V. Klika.  Modern perspectives on near-equilibrium analysis of Turing systems.  _Philosophical Transactions of the Royal Society A_ , 379(2213):20200268, 2021.
  * Krause et al. [2023] A. L. Krause, E. A. Gaffney, T. J. Jewell, V. Klika, and B. J. Walker.  Code relating to ‘Turing instabilities are not enough’, 2023.  https://github.com/AndrewLKrause/Turing-instabilities-are-not-enough-code.
  * Landge et al. [2020] A. N. Landge, B. M. Jordan, X. Diego, and P. Müller.  Pattern formation mechanisms of self-organizing reaction-diffusion systems.  _Developmental biology_ , 460(1):2–11, 2020.
  * Laurent and Kellershohn [1999] M. Laurent and N. Kellershohn.  Multistability: a major means of differentiation and evolution in biological systems.  _Trends in biochemical sciences_ , 24(11):418–422, 1999.
  * Marcon et al. [2016] L. Marcon, X. Diego, J. Sharpe, and P. Müller.  High-throughput mathematical analysis identifies Turing networks for patterning with equally diffusing signals.  _eLife_ , 5:e14022, 2016.
  * Meinhardt and Gierer [2000] H. Meinhardt and A. Gierer.  Pattern formation by local self-activation and lateral inhibition.  _Bioessays_ , 22(8):753–760, 2000.
  * Murray [1982] J. Murray.  Parameter space for Turing instability in reaction diffusion mechanisms: a comparison of models.  _Journal of Theoretical Biology_ , 98(1):143–163, 1982.
  * Murray [2003] J. D. Murray.  _Mathematical Biology: II: Spatial Models and Biomedical Applications_ , volume 3.  Springer, 2003.
  * Painter and Hillen [2011] K. J. Painter and T. Hillen.  Spatio-temporal chaos in a chemotaxis model.  _Physica D: Nonlinear Phenomena_ , 240(4-5):363–375, 2011.
  * Painter et al. [2015] K. J. Painter, J. Bloomfield, J. Sherratt, and A. Gerisch.  A nonlocal model for contact attraction and repulsion in heterogeneous cell populations.  _Bulletin ofMmathematical Biology_ , 77:1132–1165, 2015.
  * Potts and Lewis [2016] J. R. Potts and M. A. Lewis.  How memory of direct animal interactions can lead to territorial pattern formation.  _Journal of the Royal Society Interface_ , 13(118):20160059, 2016.
  * Satnoianu et al. [2000] R. A. Satnoianu, M. Menzinger, and P. K. Maini.  Turing instabilities in general systems.  _Journal of mathematical biology_ , 41:493–512, 2000.
  * Scholes et al. [2019] N. S. Scholes, D. Schnoerr, M. Isalan, and M. P. Stumpf.  A comprehensive network atlas reveals that Turing patterns are common but not robust.  _Cell systems_ , 9(3):243–257, 2019.
  * Siegal-Gaskins et al. [2009] D. Siegal-Gaskins, E. Grotewold, and G. D. Smith.  The capacity for multistability in small gene regulatory networks.  _BMC systems biology_ , 3(1):1–14, 2009.
  * Suzuki et al. [2021] K. Suzuki, S. Nakaoka, S. Fukuda, and H. Masuya.  Energy landscape analysis elucidates the multistability of ecological communities across environmental gradients.  _Ecological Monographs_ , 91(3):e01469, 2021.
  * Turing [1952] A. M. Turing.  The chemical basis of morphogenesis.  _Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences_ , 237(641):37–72, 1952.
  * Uecker et al. [2014] H. Uecker, D. Wetzel, and J. D. Rademacher.  pde2path-a Matlab package for continuation and bifurcation in 2D elliptic systems.  _Numerical Mathematics: Theory, Methods and Applications_ , 7(1):58–106, 2014.
  * van de Koppel and Rietkerk [2004] J. van de Koppel and M. Rietkerk.  Spatial interactions and resilience in arid ecosystems.  _The American Naturalist_ , 163(1):113–121, 2004.
  * Vastano et al. [1988] J. A. Vastano, J. E. Pearson, W. Horsthemke, and H. L. Swinney.  Turing patterns in an open reactor.  _The Journal of chemical physics_ , 88(10):6175–6181, 1988.
  * Veerman and Doelman [2013] F. Veerman and A. Doelman.  Pulses in a gierer–meinhardt equation with a slow nonlinearity.  _SIAM Journal on Applied Dynamical Systems_ , 12(1):28–60, 2013.
  * Villar-Sepúlveda and Champneys [2023] E. Villar-Sepúlveda and A. Champneys.  Degenerate Turing bifurcation and the birth of localized patterns in activator-inhibitor systems.  _SIAM Journal on Applied Dynamical Systems_ , 22(3):1673–1709, 2023.
  * Vittadello et al. [2021] S. T. Vittadello, T. Leyshon, D. Schnoerr, and M. P. Stumpf.  Turing pattern design principles and their robustness.  _Philosophical Transactions of the Royal Society A_ , 379(2213):20200272, 2021.
  * Walker et al. [2023] B. J. Walker, A. K. Townsend, A. K. Chudasama, and A. L. Krause.  VisualPDE: rapid interactive simulations of partial differential equations.  _arXiv:2308.01245_ , 2023.
  * Wei and Winter [2013] J. Wei and M. Winter.  _Mathematical aspects of pattern formation in biological systems_ , volume 189.  Springer Science & Business Media, 2013.
  * Woolley et al. [2021] T. E. Woolley, A. L. Krause, and E. A. Gaffney.  Bespoke Turing systems.  _Bulletin of Mathematical Biology_ , 83:1–32, 2021.
  * Wyss and Jorgensen [1998] G. D. Wyss and K. H. Jorgensen.  A users guide to LHS: Sandias Latin hypercube sampling software.  Technical report, Sandia National Lab.(SNL-NM), Albuquerque, NM (United States), 1998.