Turing conditions for pattern forming systems on evolving manifolds

TL;DR

This paper develops general, practical conditions for predicting pattern formation in reaction-diffusion systems on evolving domains, extending classical Turing instability analysis to dynamic, non-autonomous settings.

Contribution

It introduces versatile differential inequalities for Turing instability onset on time-dependent manifolds, generalizing existing static domain criteria and applicable to various domain evolutions.

Findings

Derived sufficient conditions for diffusion-driven instabilities on evolving domains.

Validated the conditions on different domain growth laws and oscillatory states.

Extended the framework to higher-order spatial systems.

Abstract

The study of pattern-forming instabilities in reaction-diffusion systems on growing or otherwise time-dependent domains arises in a variety of settings, including applications in developmental biology, spatial ecology, and experimental chemistry. Analyzing such instabilities is complicated, as there is a strong dependence of any spatially homogeneous base states on time, and the resulting structure of the linearized perturbations used to determine the onset of instability is inherently non-autonomous. We obtain general conditions for the onset and structure of diffusion driven instabilities in reaction-diffusion systems on domains which evolve in time, in terms of the time-evolution of the Laplace-Beltrami spectrum for the domain and functions which specify the domain evolution. Our results give sufficient conditions for diffusive instabilities phrased in terms of differential inequalities which are both versatile and straightforward to implement, despite the generality of the studied problem. These conditions generalize a large number of results known in the literature, such as the algebraic inequalities commonly used as a sufficient criterion for the Turing instability on static domains, and approximate asymptotic results valid for specific types of growth, or specific domains. We demonstrate our general Turing conditions on a variety of domains with different evolution laws, and in particular show how insight can be gained even when the domain changes rapidly in time, or when the homogeneous state is oscillatory, such as in the case of Turing-Hopf instabilities. Extensions to higher-order spatial systems are also included as a way of demonstrating the generality of the approach.

Equations

$$\frac{d}{dt}{\int }_{\hat{\Omega }(t)}\mathbf{u}d\Omega ={\int }_{\hat{\Omega }(t)}\left(-\nabla \cdot \mathbf{j}+\mathbf{f}(\mathbf{u})\right)d\Omega ,$$

$$\frac{d}{dt}{\int }_{\hat{\Omega }(t)}\mathbf{u}d\Omega ={\int }_{\hat{\Omega }(t)}\left(\frac{\partial \mathbf{u}}{\partial t}+{\nabla }_{\Omega (t)}\cdot (\mathbf{Qu})\right)d\Omega ,$$

$$\frac{\partial \mathbf{u}}{\partial t}+{\nabla }_{\Omega (t)}\cdot (\mathbf{Qu})=D{\nabla }_{\Omega (t)}^2\mathbf{u}+\mathbf{f}(\mathbf{u}).$$

$$\mathbf{X}=\left(r_1(t){\chi }_1(\mathbf{x}),\dots ,r_N(t){\chi }_N(\mathbf{x})\right),$$

$${\nabla }_{\Omega (t)}\cdot \mathbf{Q}=\frac{\partial }{\partial t}\left(\log \left(|\det G{|}^{\frac{1}{2}}\right)\right).$$

$$\frac{\partial \mathbf{u}}{\partial t}=\frac{D}{|\det G{|}^{\frac{1}{2}}}\sum _{i,j=1}^N\frac{\partial }{\partial x_i}\left(|\det G{|}^{\frac{1}{2}}{G}_{ij}^{-1}\frac{\partial \mathbf{u}}{\partial x_j}\right)-\frac{\partial }{\partial t}\left(\log \left(|\det G{|}^{\frac{1}{2}}\right)\right)\mathbf{u}+\mathbf{f}(\mathbf{u}).$$

$${\nabla }_{\Omega (t)}^2\Psi =-\rho \Psi ,$$

$${\psi }_k(\mathbf{x})={\Psi }_k\left(\frac{X_1}{r_1(t)},\dots ,\frac{X_N}{r_N(t)}\right)={\Psi }_k\left({\chi }_1(\mathbf{x}),\dots ,{\chi }_N(\mathbf{x})\right),$$

$$\frac{1}{|\det G{|}^{\frac{1}{2}}}\sum _{i,j=1}^N\frac{\partial }{\partial x_i}\left(|\det G{|}^{\frac{1}{2}}{G}_{ij}^{-1}\frac{\partial {\psi }_k(\mathbf{x})}{\partial x_j}\right)=-{\rho }_k(t){\psi }_k(\mathbf{x}).$$

$$\frac{\dot{\mu }(t)}{\mu (t)}=\frac{\partial }{\partial t}\left(\log \left(|\det G{|}^{1/2}\right)\right)$$

$$\frac{d\mathbf{U}}{dt}=-\frac{\dot{\mu }(t)}{\mu (t)}\mathbf{U}+\mathbf{f}(\mathbf{U}),\mathbf{U}(0)={\mathbf{U}}^{\ast },$$

$$\mathbf{u}(\mathbf{x},t)=\mathbf{U}(t)+ϵ\frac{{\psi }_k(\mathbf{x})}{\mu (t)}{\mathbf{V}}_k(t),$$

$$\frac{d{\mathbf{V}}_k}{dt}=-{\rho }_k(t)D{\mathbf{V}}_k+J(\mathbf{U}){\mathbf{V}}_k,$$

$$\frac{\partial \mathbf{u}}{\partial t}+{\nabla }_{\Omega (t)}\cdot (\mathbf{Qu})=D\mathcal{P}\left({\nabla }_{\Omega (t)}^2\right)\mathbf{u}+\mathbf{f}(\mathbf{u}),$$

$$\frac{\partial \mathbf{u}}{\partial t}=D\mathcal{P}\left(\frac{1}{|\det G{|}^{\frac{1}{2}}}\sum _{i,j=1}^N\frac{\partial }{\partial x_i}\left(|\det G{|}^{\frac{1}{2}}{G}_{ij}^{-1}\frac{\partial \mathbf{u}}{\partial x_j}\right)\right)-\frac{\partial }{\partial t}\left(\log \left(|\det G{|}^{\frac{1}{2}}\right)\right)\mathbf{u}+\mathbf{f}(\mathbf{u}).$$

$$\mathcal{P}\left({\nabla }_{\Omega (t)}^2\right){\psi }_k=\sum _{\ell =1}^p{\left({\nabla }_{\Omega (t)}^2\right)}^{\ell }{\psi }_k=\sum _{\ell =1}^p{\left(-{\rho }_k(t)\right)}^{\ell }{\psi }_k=\mathcal{P}\left(-{\rho }_k(t)\right){\psi }_k$$

$$\frac{d{\mathbf{V}}_k}{dt}=D\mathcal{P}\left(-{\rho }_k(t)\right){\mathbf{V}}_k+J(\mathbf{U}){\mathbf{V}}_k.$$

$$\ddot{Y}+F(t)\dot{Y}+G(t)Y=0.$$

$$G(t)\le -\frac{\ddot{\Phi }}{\Phi }-\frac{\dot{\Phi }}{\Phi }F(t),t\in \mathcal{I}.$$

$$\ddot{Y}+F(t)\dot{Y}-\left(\frac{\ddot{\Phi }}{\Phi }+\frac{\dot{\Phi }}{\Phi }F(t)+H(t)\right)Y=0.$$

$$\dot{Z}=-Z^2-F(t)Z+\frac{\ddot{\Phi }}{\Phi }+\frac{\dot{\Phi }}{\Phi }F(t)+H(t).$$

$$\dot{Z}=-Z^2-F(t)Z+\frac{\ddot{\Phi }}{\Phi }+\frac{\dot{\Phi }}{\Phi }F(t)+H(t)\ge -Z^2-F(t)Z+\frac{\ddot{\Phi }}{\Phi }+\frac{\dot{\Phi }}{\Phi }F(t)={\dot{Z}}_1,$$

$${\dot{Z}}_1=-Z_1^2-F(t)Z_1+\frac{\ddot{\Phi }}{\Phi }+\frac{\dot{\Phi }}{\Phi }F(t),$$

$$Y(t)=Y(t_0)\exp \left({\int }_{t_0}^{t}Z(s)ds\right)\ge Y(t_0)\exp \left({\int }_{t_0}^t{Z}_1(s)ds\right)=\Phi (t),$$

$$G(t)\le -\frac{\ddot{\mu }}{\mu }-2\delta \frac{\dot{\mu }}{\mu }-{\delta }^2-\left(\frac{\dot{\mu }}{\mu }+\delta \right)F(t).$$

$$G(t)<-\frac{\ddot{\mu }}{\mu }-\frac{\dot{\mu }}{\mu }F(t),\text{for all}t\in \mathcal{I}.$$

$$G(t)\le -\frac{\ddot{\mu }}{\mu }-2\delta \frac{\dot{\mu }}{\mu }{t}^{-1}-\delta (\delta -1)t^{-2}-\left(\frac{\dot{\mu }}{\mu }+\delta t^{-1}\right)F(t).$$

$$G(t)<-\frac{\ddot{\mu }}{\mu }-\frac{\dot{\mu }}{\mu }F(t),\text{for all}t\in \mathcal{I}.$$

$$G(t)>-\frac{\ddot{\mu }}{\mu }-\frac{\dot{\mu }}{\mu }F(t),\text{for all}t\in \mathcal{I}.$$

$$G(t)-H(t)=-\frac{\ddot{\mu }}{\mu }-\frac{\dot{\mu }}{\mu }F(t),\text{for all}t\in \mathcal{I},$$

Full text

keywords: pattern formation, Turing instability, evolving spatial domains, reaction-diffusion

AMS Subject Classifications: 35B36, 92C15, 70K50, 58C40, 58J32

1 Introduction

Since Turing first proposed reaction–diffusion systems as a model for pattern formation [112] much work has been done to understand the theoretical, as well as chemical and biological aspects of this mechanism [2, 65]. Many authors have elucidated the importance of domain size and shape on the formation of patterns, and the impact of geometry on the kinds of admissible patterns that can arise due to a Turing instability [75, 97]. While Turing’s original reaction–diffusion theory postulated the existence of a pre-pattern before growth occurs, in the past few decades biological and theoretical evidence has suggested that growth of the patterning field itself influences the pattern forming potential of a system, and modulates the emergent patterns [7, 13, 59, 64, 73, 86]. Since reaction–diffusion systems are more difficult to analyze on growing domains, pattern formation has been typically considered on different–sized static domains to simulate very slow growth [116]. This requires the reaction and diffusion of the chemical species to occur on a much faster timescale than the growth [45, 63], and to be independent of the growth, although this assumption is certainly not valid for all systems in biology and chemistry.

There have been a variety of studies connecting growth and pattern formation in reaction-diffusion systems. Uniform and isotropic domain growth in one–dimensional reaction–diffusion systems in slow and fast growth regimes was considered by [20], and frequency-doubling of the emergent Turing patterns was demonstrated. This frequency-doubling was discussed with the aim to help resolve the (lack of) robustness of pattern formation in reaction-diffusion systems [2, 4]. More recently, it was shown that such frequency-doubling can depend somewhat sensitively on the kind of growth rates involved, even in a 1-D domain growing isotropically [114]. Turing and Turing-Hopf instabilities of the FitzHugh-Nagumo system in an exponentially and isotropically growing square were studied in [11], and this work suggested that anisotropy and curvature are important considerations for extending their analysis. Such instabilities were also studied on isotropically growing spherical and toroidal domains [92]. A general formulation of reaction–diffusion theory on isotropically evolving one and two–dimensional manifolds was given in [83], with motivation from biological settings where growth and curvature both play a role in organism development. Corrections to the classical conditions for Turing instabilities in the case of slow isotropic growth were derived by [63], while [35] considered a type of quasi-asymptotic stability, although large deviations from an approximately homogeneous state at finite time were not considered. In contrast, [45] considered domain growth based on Lyapunov stability, which captures large deviations from the reference state at finite time before growth saturates, thereby capturing some of the history dependence inherent in the growth dynamics. The study of [45] was able to relax a number of restrictive assumptions made in [35] and [63], although the results were obtained for fairly specific special cases. Beyond computing linear instability criteria, [16] analytically explored how patterns change and evolve under growth by exploiting the framework of amplitude equations. While analytical results on mode competition and selection can be valuable, these are often extremely limited as they only apply near the bifurcation boundary in the parameter space, and they become computationally intractable in many cases of interest [54]. While systems with time-dependent diffusion coefficients have been studied via asymptotic and Floquet-theoretic methods [68], such results break down for domain evolution due to the dynamic nature of the base state.

Most of the above models of reaction–diffusion systems on growing domains only analyzed the case of isotropic (or apical) growth, which is unable to recapitulate the full range of complex biological structures found in developing organisms [17, 82, 104, 113]. Investigating arbitrary anisotropic growth in the context of biological patterning is a natural extension to reaction–diffusion theories of pattern formation, and has been considered in biomechanical models of growth across a range of tissues and organisms [1, 8, 69, 91]. Additionally, contraction and other complex tissue movements have been observed in embryogenesis, suggesting the need to further generalize models of growth and domain restructuring over time beyond monotonically increasing domains [1, 77, 111, 120]. The influence of non–uniform domain growth on one–dimensional reaction–diffusion systems, including apical or boundary growth, was considered in [18], while [87] studied concentration-dependent growth of a scalar reaction-diffusion equation on a time–dependent manifold. Anisotropic growth, consisting of independent dilations of an underlying manifold in each orthogonal Euclidean coordinate, was recently studied in [52]. Some experimental models of reaction–diffusion processes on growing domains (using, for instance, photosensitive reactions) have been explored, although these typically involve apical or otherwise spatially-dependent forms of growth [71, 49]. Recent experiments have also explored hydrodynamic instabilities in time-dependent domains, finding important impacts of the growth on the nature of such instabilities and subsequent pattern evolution [31].

Difficulties arise when the local rate of volume expansion or contraction depends on the spatial coordinates (or more generally, the morphogen concentrations themselves), which induces an advection term from mass conservation of the respective chemical species. Such systems are exceptionally difficult to analyze due to spatial heterogeneity in diffusion and advection, in addition to their non-autonomous nature. Still, as we shall later discuss, some forms of anisotropic growth permit volume expansion or contraction which is global, depending only on time and not on the spatial coordinates. It is this class of growing domains for which we provide a general method to compute the instability of a spatially homogeneous equilibrium. This generalizes much of the above literature, and provides a way to compute a time-dependent instability criterion for a large class of growth functions for reaction–diffusion systems posed on smooth, compact, and simply connected, but otherwise arbitrary, manifolds.

The remainder of this paper is organized as follows. In Sec. 2 we outline the derivation of a general reaction-diffusion model on time-evolving domains. We give the precise mathematical formulation of component-wise dilational growth considered throughout the paper, and outline the general spectral problem. We also discuss several difficulties in the analysis of such systems, necessitating the need for new approaches from those often employed in the literature. In Sec. 3, we present the main theoretical results for systems of two reaction-diffusion equations on evolving domains. After first obtaining a general instability result for second-order ODE through a comparison principle, we derive an instability criterion which generalizes the Turing conditions for diffusive instabilities to the non-autonomous case. We also discuss various reductions of these conditions, highlighting that they collapse to the standard Turing conditions on static domains in the appropriate limits. In addition to systems of reaction-diffusion equations, we also discuss the extension of our results to systems with higher-order space derivatives which are also heavily studied in the pattern formation literature, and capture non-local effects in a variety of models. In Sec. 4-7 we provide a number of applications of the theory, and compare our results with direct numerical simulations of various reaction diffusion systems on growing domains in one, two, and three spatial dimensions. We generalize several examples from the literature without restriction to asymptotic regimes, and consider novel classes of domain evolution which have not yet been considered. In particular, applications of our approach for evolving one-dimensional intervals are described in Sec. 4 (including the situation where the evolution is non-monotone), whereas the more complicated applications to evolving manifolds in two or more dimensions are given in Sec. 5. In addition to growing domains, our approach allows us to consider domains which evolve yet preserve area or volume (as discussed in Sec. 6), which has seemingly not been considered previously. We also give some examples related to higher-order systems on evolving manifolds in Sec. 7. We finally discuss and summarize our findings in Sec. 8.

2 General model and diffusive instability framework

2.1 Reaction-diffusion systems on evolving domains

We consider a manifold $\Omega(t)\subseteq\mathbb{R}^{N}$ which grows in a dilational manner along each Euclidean coordinate axis. We also assume that the domain $\Omega(t)$ is compact, simply connected, smooth, and Riemannian, in order to ensure that the spectrum of the Laplace-Beltrami operator is discrete and non-negative for all time. Concentrations on manifolds with boundary will be subject to Neumann data at the boundary. We shall restrict our attention to growth for which the time evolution of $\Omega(t)$ results in spatially homogeneous volume expansion or contraction, and shall make this statement more precise later. The case of locally varying volume expansion or contraction results in strongly non-uniform growth which we do not consider.

Let $\hat{\Omega}(t)$ be a volume element of the manifold, such that $\hat{\Omega}(t)\subset\Omega(t)$. Let $\mathbf{u}(\mathbf{X},t)$, $\mathbf{u}:\Omega(t)\times[0,\infty)\rightarrow\mathbb{R}^{n}$ be a concentration function defined on the manifold $\Omega(t)$. Here $\mathbf{u}$ may describe the concentration of $n\geq 2$ chemical species, or morphogens, on the manifold $\Omega(t)$, though other interpretations such as cells or effective genetic circuits use the same mathematical formulations [48]. We assume that $\mathbf{u}$ is $C^{1}([0,\infty))$ in time and $C^{2}(\Omega(t))$ in the spatial coordinates. We note that formalizing this space of functions is easier to do after mapping back to a static domain, which we will also do for analytic and numerical convenience shortly.

The conservation of mass equation reads

[TABLE]

where $\mathbf{j}$ denotes the fluxes of concentrations $\mathbf{u}$, $\mathbf{f}\in C^{1}(\mathbb{R}^{n})$ is the function denoting the reaction kinetics, and $d\Omega$ is the local volume element on the manifold. Using Reynold’s transport theorem on the left hand side of Eq. (1), we have

[TABLE]

where $\mathbf{Q}$ is the local velocity vector field generated by changes in the manifold $\Omega(t)$. We denote $\nabla_{\Omega(t)}\cdot$ as the divergence operator on $\Omega(t)$ and $\nabla_{\Omega(t)}^{2}$ to be the Laplace–Beltrami operator on $\Omega(t)$.

By applying Eq. (2) to Eq. (1) and using Fick’s law of diffusion, we have the reaction–diffusion–advection system

[TABLE]

Here $D=\text{diag}(d_{1},\dots,d_{n})$ is the matrix of diffusion parameters. The term $\nabla_{\Omega(t)}\cdot(\mathbf{Q}\mathbf{u})$ can be written as $\mathbf{Q}\cdot\nabla_{\Omega(t)}\mathbf{u}+\mathbf{u}\cdot\nabla_{\Omega(t)}\mathbf{Q}$. We note that the term $\mathbf{Q}\cdot\nabla_{\Omega(t)}\mathbf{u}$ corresponds to advection due to local growth of the manifold, whereas the $\mathbf{u}\cdot\nabla_{\Omega(t)}\mathbf{Q}$ term corresponds to dilution of the concentrations $\mathbf{u}$ due to local volume changes. If $\Omega(t)$ is a manifold with boundary $\partial\Omega(t)$, we assume no flux conditions $\frac{\partial\mathbf{u}}{\partial\mathbf{n}}=\mathbf{0}$ for $\mathbf{X}\in\partial\Omega(t)$.

2.2 Domain evolution as dilations in each orthogonal direction

We consider the case where the evolution of the manifold is such that volume expansion or contraction does not vary locally, in other words such that $\nabla_{\Omega(t)}\cdot\mathbf{Q}$ depends strictly on time and never on spatial coordinates. Other kinds of growth, such as apical growth or anisotropic growth of surfaces, may result in spatially dependent volume expansion [52], and while interesting, we do not consider this manner of growth.

Consider moving coordinates $\mathbf{X}$ written as

[TABLE]

for stationary coordinates on the initial manifold $\mathbf{x}=(x_{1},\dots,x_{N})\in\Omega^{*}=\Omega(t=0)$. In addition to covering all cases of general dilational growth (where the dilation may be different along each coordinate), this assumption will ensure that the metric tensor for these coordinates, $G$, will have the property that $\det G$ is multiplicatively separable in time and space. Here each coordinate again has an independent dilation function $r_{j}(t)$, though each $X_{j}$ depends possibly on multiple stationary coordinates. When $r_{j}(t)=r(t)$ for all $j=1,\dots,N$, we have isotropic evolution of the manifold $\Omega(t)$. When at least two $r_{j}(t)$’s differ, then we have anisotropic evolution which is still dilational in the individual orthogonal Cartesian directions.

In order to remove the advection term induced by the growing manifold, $\nabla_{\Omega(t)}\cdot(\mathbf{Q}\mathbf{u})$, we apply a change of variables to a moving coordinate system. As the space and time variables in $\mathbf{X}$ are separable, and noting $\nabla_{\Omega(t)}\cdot(\mathbf{Q}\mathbf{u})=(\nabla_{\Omega(t)}\cdot\mathbf{Q})\mathbf{u}+\mathbf{Q}\cdot(\nabla_{\Omega(t)}\mathbf{u})$, the change of variable will contribute a term $-\mathbf{Q}\cdot(\nabla_{\Omega(t)}\mathbf{u})$, canceling the latter term. After the coordinate change, we will have a contributing advection term of the form $(\nabla_{\Omega(t)}\cdot\mathbf{Q})\mathbf{u}$. We find [52]

[TABLE]

From (5), we see that the manner of growth for which volume expansion or contraction is spatially homogeneous is equivalent to considering a coordinate chart such that the time derivative of $\log(|\det G|)$ is independent of space, i.e. a coordinate chart for which $\det G$ is multiplicatively separable in space and time variables. Considering only such moving coordinates (4), we then have that (3) becomes [52]

[TABLE]

The Laplace-Beltrami operator on the fixed reference manifold $\Omega^{*}$ is time-dependent, as the coordinates (4) depend explicitly on time.

2.3 General linear instability analysis

We now consider diffusion-driven instabilities arising from systems of the form (6). Consider first the eigenvalue problem

[TABLE]

which is held subject to $\frac{\partial\Psi}{\partial\mathbf{n}}=0$ for $\mathbf{X}\in\partial\Omega(t)$ when $\Omega(t)$ has boundary. From the assumptions made earlier on $\Omega(t)$, for any fixed time $t\geq 0$, we have that a non-negative spectrum $\rho_{k}(t)$ exists, where $0=\rho_{0}(t)<\rho_{1}(t)\leq\rho_{2}(t)\leq\cdots\rightarrow\infty$. In general, for manifolds of dimension greater than one, $k$ denotes a multi-index. As the growth functions are assumed smooth, and $\Omega(t)$ is assumed a simply connected Riemannian manifold with smooth boundary for all $t\geq 0$, then we shall assume that $\Omega(t)$ is such that the spectrum can be continued as a smooth function of time, with $\rho_{k}(t)>0$ for all $k\geq 1$. For our purposes, we assume any given $\Omega(t)$ permits such a construction, as we are concerned with dynamics on a prescribed $\Omega(t)$. We denote the corresponding eigenfunctions by $\Psi_{k}(\mathbf{X})$. Constructing such eigenvalues and eigenfunctions can be very difficult, and although our results only require existence rather than explicit construction, we will give examples later for domains where such constructions can be carried out.

If we carry out the change of coordinates (4), and note that the stationary form of each eigenfunction is

[TABLE]

the eigenvalue problem (7) is put into the form

[TABLE]

In the special case where domain evolution is isotropic, that is $r_{1}(t)=\cdots=r_{N}(t)=r(t)$, the Laplace-Beltrami operator simplifies so that we have $\frac{1}{r(t)^{2}}\nabla_{\Omega(0)}^{2}\psi_{k}(\mathbf{x})=-\frac{\rho_{k}(0)}{r(t)^{2}}\psi_{k}(\mathbf{x})$, and hence $\rho_{k}(t)=\frac{\rho_{k}(0)}{r(t)^{2}}$. This is of course not true in general, and for more complicated growth finding $\rho_{k}(t)$ can be involved. However, under growth of the form (4), each eigenfunction is stationary, and each eigenvalue is a function of time as smooth as the dilations $r_{i}(t)$, granting existence of the $\rho_{k}(t)$.

From the form of growth assumed, we have that volume expansion is not dependent on space, so we can write

[TABLE]

for some function $\mu(t)$. As $|\det{G}|=|\mathcal{G}_{1}(t)\mathcal{G}_{2}(\mathbf{x})|$, then $\mu(t)=|\mathcal{G}_{1}(t)|^{1/2}$.

A spatially uniform solution to (6), $\mathbf{u}(\mathbf{x},t)=\mathbf{U}(t)$, is governed by the equation

[TABLE]

where we choose $\mathbf{U}^{*}$ to satisfy $\mathbf{f}(\mathbf{U}^{*})=\mathbf{0}$. We choose this initial condition so that the dynamics will initially agree with those of a time-independent steady state in the absence of growth. In this way, if there is no growth (or, more generally, when there is no net volume expansion or contraction so that $\dot{\mu}\equiv 0$ for all $t\geq 0$), then the exact solution to (11) is $\mathbf{U}(t)\equiv\mathbf{U}^{*}$, which is what is assumed when deriving the standard Turing conditions on a static domain. Therefore, (11) generalizes the static uniform base state to account for dilution due to growth.111We note that for complex reaction-diffusion systems spatio-temporal base states consisting of plane waves are also possible [46]. However, as our concern is with generalising the Turing conditions to account for evolving domains, we only consider spatially uniform base states. Additionally, plane waves do not generalize well to manifolds which are not flat.

We consider a perturbation of this spatially homogeneous solution in order to determine its stability. Although the solution to (11) may tend to a steady state, this is not required, and for many examples will not occur. We choose a general spatial perturbation of the form

[TABLE]

where $\psi_{k}(\mathbf{x})$ is the $k$th scaled eigenfunction in (8) with corresponding Laplace-Beltrami eigenvalue $\rho_{k}(t)$. For each $k=1,2,\dots,$ we then have the following linearized problem for the growth or decay of the perturbation

[TABLE]

which is an ODE for the unknown function $\mathbf{V}_{k}(t)$, the long-time asymptotic behavior of which determines the stability or instability of the perturbation (12). The matrix $J(\mathbf{U})$ denotes the (in general, time-dependent) Jacobian matrix corresponding to the linearization of $\mathbf{f}$ at $\mathbf{u}=\mathbf{U}(t)$.

Equation (13) results in a solution $\mathbf{V}_{k}(t)$, and we say that a perturbation (12) is asymptotically stable for a given $k\in\mathbb{N}$ provided $\mu(t)^{-1}|\mathbf{V}_{k}(t)|\rightarrow 0$ as $t\rightarrow\infty$. We say that the perturbation (12) is asymptotically unstable for a given $k\in\mathbb{N}$ provided $\mu(t)^{-1}|\mathbf{V}_{k}(t)|\rightarrow\infty$ as $t\rightarrow\infty$. If $V_{k}(t)$ satisfies neither of these, then we might say that the perturbation is neutrally stable or unstable, depending upon the context. This is akin to the classical Turing perturbation for which $\mathbf{V}_{k}(t)=\mathbf{C}\exp(\lambda(k)t)$, where $\mathbf{C}\in\mathbb{R}^{N}$ is a constant vector, and $\lambda(k)\in\mathbb{C}$, in which case the perturbation is stable if $\text{Re}(\lambda(k))<0$ and unstable if $\text{Re}(\lambda(k))>0$. We will also discuss conditions for transient stability or instability, wherein a perturbation may decay or grow for some set of time, as this can generate a pattern in the fully nonlinear setting (though such linear analysis cannot guarantee that an instability leads to a patterned state in finite time periods). As we shall see, transient instabilities play a much larger role in evolving domains than any asymptotic stability criterion.

2.4 Systems with higher-order spatial derivatives

There are various applications for which higher-order spatial derivatives are used, with the Swift-Hohenberg equation [105, 21], Cahn-Hilliard equation [9, 43, 78], and Kuramoto-Sivashinsky [55, 102, 27, 39, 88] equations some common examples of pattern forming systems with fourth-order spatial derivatives, with related models of even higher-order arising in applications [50, 81]. Such higher-order derivatives often represent nonlocal interactions, and this has been extensively applied in biological applications such as cellular signalling and ecological interactions; see [14, 79] and Chapter 11 of [74]. We show here how to carry out similar analysis outlined above for such higher-order equations on evolving spatial domains. Similar to what we have done in Section 2.1, we may write such higher-order systems in the form

[TABLE]

where $\mathcal{P}$ is a polynomial of degree $p\geq 1$ satisfying $\mathcal{P}(0)=0$, and the second term involving $\mathbf{Q}$ again arises from a conservation principle on the evolving manifold. On manifolds with boundary, we consider a generalisation of no-flux boundary conditions $\frac{\partial^{2\ell-1}\mathbf{u}}{\partial^{2\ell-1}\mathbf{n}}=\mathbf{0}$ for $\mathbf{X}\in\partial\Omega(t)$, where $\ell=1,2,\dots,p$. All other quantities are as defined earlier, with the only difference between equation (14) and the earlier discussed equation (3) being the more general operator involving the spatial derivatives. Dynamics from the complex Swift-Hohenberg equation, with $u^{3}$ replaced by $|u|^{2}u$, were considered on an evolving domain by [46] and [54], with the domain being a time-dependent interval.

Mapping the problem (14) to the stationary frame (assuming that the manifold $\Omega(t)$ obeys all properties required in Section 2.2), we find that (14) is put into the form

[TABLE]

Note that the equation governing a spatially uniform state is the same as in (11).

Regarding the spectral problem, as the operator $\mathcal{P}$ is a polynomial of the Laplace-Beltrami operator, we have

[TABLE]

for any spatial eigenfunction $\psi_{k}$ satisfying (7) and (8). Therefore, assuming a spatial perturbation $\psi_{k}$ involving the $k$th spatial eigenfunction as in (12), we obtain the problem

[TABLE]

The qualitative analysis for (17) is the same as that discussed for (13).

Of course, depending upon the form of the boundary conditions, there can be other spatial eigenfunctions for higher-order problems; that is to say, the Laplace-Beltrami eigenfunctions $\psi_{k}$ described in (7) and (8) are in general a subset of the eigenfunctions possible for a given higher-order eigenvalue problem such as (16). The study of such cases would involve the classification of eigenfunctions and eigenvalues of the spatial problem $\mathcal{P}\left(\nabla_{\Omega(t)}^{2}\right)\hat{\psi}_{k}=-\zeta_{k}\hat{\psi}_{k}$, where $\hat{\psi}_{k}$ need not be an eigenfunction of the standard Laplace-Beltrami problem $\nabla_{\Omega(t)}^{2}\psi_{k}=\rho_{k}\psi_{k}$. For the purpose of this paper, we shall primarily consider only the relatively simple case of Laplace-Beltrami eigenfunctions, which will be sufficient for studying the instability problem. Of course, given a desired operator, $\mathcal{P}\left(\nabla_{\Omega(t)}^{2}\right)$, one may perform similar calculations and scaling to what we do for the standard Laplace-Beltrami case, in order to find the time-dependent spectrum $\zeta_{k}=\zeta_{k}(t)$. For most higher-order elliptic operators of practical importance, $\zeta_{k}$ will be bounded away from $-\infty$ and non-decreasing, $-\infty<\zeta_{0}\leq\zeta_{1}<\cdots$ with $\zeta_{k}\rightarrow\infty$ as $k\rightarrow\infty$ [15], although for many commonly studied operators on bounded domains the spectrum is non-negative [84, 58, 56]. We avoid a discussion on the classification of such problems, noting that for many higher-order problems the behavior of the spectrum appears to be an open problem. Assuming one can find the spectrum $\zeta_{k}$, we briefly comment on how to make use of this in Section 3.3.

2.5 Difficulties arising in the study of such systems

The study of the asymptotic stability of systems of the form (13) is made difficult for a number of reasons, which we now outline.

The base states governed by (11) depend on the global rate of volume expansion or contraction, and hence are time-varying. Therefore, we expand and linearize the reaction-diffusion system about a spatially uniform yet temporally varying base state, resulting in a non-autonomous Jacobian matrix. This non-autonomy is in addition to the non-autonomy due to the spectrum of the evolving domain, hence the systems for the linearization $\mathbf{V}_{k}(t)$ given in (13) are non-autonomous in all components rather than just in diagonal components. As the Jacobian $J$ depends on the specific form of the nonlinear reaction kinetics $\mathbf{f}$, non-autonomous entries in $J$ may be non-monotone, even for monotone growth functions. Compared to their autonomous counterparts, there is very little general theory for the dynamics of such systems.

We note that many works in this area (see, for instance, [62]) will attempt to overcome this complication by assuming a time-independent steady state solution of the ODE system governing the reaction kinetics in the presence of growth. This would be equivalent to obtaining a fixed point of the right hand side of (11), i.e., to finding an algebraic solution $\mathbf{U}^{*}$ of the algebraic equation $\frac{\dot{\mu}(t)}{\mu(t)}\mathbf{U}=\mathbf{f}(\mathbf{U})$. There are two problems with this, one regarding feasibility and one more philosophical. Regarding feasibility, a time-independent steady state $\mathbf{U}^{*}$ exists only if $\frac{\dot{\mu}(t)}{\mu(t)}$ is identically equal to a constant for all time, which is restrictive of the kinds of growth considered. One exception is to consider a state which is identically zero, provided that the reaction kinetics $\mathbf{f}$ permit this. For a zero state, the volume expansion or contraction will still permit a zero state. Of course, this is then fairly restrictive on the form of the reaction kinetics, particularly in light of the fact that for many physical or biological systems, loss of stability of a positive steady state is useful for applications. Unlike what is done in the static-domain case, shifting an equilibrium to the zero state would influence the dynamics due to the dilution term, as (13) would then no longer be a homogeneous system.

In order to remedy this, one may be tempted to instead consider the limit $t\rightarrow\infty$, for which taking either of these quantities to be constant (at least in the case of growth no more rapid than exponential) is seemingly more sensible. This leads to the second, more philosophical, problem. If one is interested in understanding how both growth and diffusion interact to induce the Turing instability and resulting pattern formation, then as pointed out in [45], history of the domain growth must factor into the Turing conditions in some manner. If one neglects growth in the base state, then one is arriving at the final spatially uniform state after growth has occurred and effectively obtaining Turing conditions for the final configuration of the problem domain. Depending on the properties of the growth function, mass conservation may result in drastic changes in the spatially uniform state over time, and the changes will become more drastic with an increased number of spatial dimensions. As such, we maintain this dependence on growth in the base states despite the added mathematical difficulties and complications. We shall later show that there is indeed one natural scenario for which the base state can be assumed time-independent, corresponding to domains which evolve in such a way that preserves volume and hence mass.

Regarding a second difficulty, we remark that eigenvalues are not the appropriate criterion to employ for determining the long-time dynamics of such non-autonomous systems [41, 63, 70]. Signing the real part of eigenvalues of an appropriate Jacobian matrix is the standard approach for determining the stability of an autonomous ODE system, and is the approach commonly used to deduce conditions for the Turing instability. For non-autonomous systems of the form $\dot{\mathbf{Y}}=A(t)\mathbf{Y}$, this approach is neither informative nor appropriate. For sake of demonstration, [118] provide an example of a time-dependent matrix $A(t)$ with strictly negative eigenvalues admitting a solution which grows without bound as $t\rightarrow\infty$, while [121] give an example of $A(t)$ with one positive eigenvalue that results in bounded solutions. These two counterexamples demonstrate that eigenvalues are predictive of neither stability nor instability of non-autonomous ODE systems. Furthermore, employing time-dependent eigenvalues is perhaps more dubious, and we avoid making use of eigenvalues in this manner.

A final difficulty, which is more prominent in non-autonomous systems, is the transient growth and hence the limitation in the correspondence between the asymptotic stability of the linearized system and the actual long-time evolution. In the simpler autonomous case (a static domain) it can be shown that the significance of this transient effect is limited only to a fine parameter tuning (at the fringe of the classical Turing space) [44]. However, in non-autonomous systems these transient effects can become more frequent, dependent on the wavenumber and initial conditions. In particular, it was shown that for (exponential) growth with a characteristic time-scale comparable to the characteristic time-scale of reaction kinetics, all wavenumbers above certain threshold grow in initial times yielding a breakdown of the continuum description in finite time [45]. Therefore, in order to understand the emergence of patterns when the reaction kinetics are on a compatible timescale with evolution of the domain, it is necessary to consider transient pattern formation, and hence necessary to consider a criterion for the emergence of transient instabilities. This will entail the consideration of spatial modes which may be unstable over some finite interval, before again becoming stable for large time, since by such a large time the pattern has already been selected (the actual patterning process occurs within finite time).

In light of the above, we will develop a criterion for the transient instability of specific spatial perturbations, akin to the Turing criteria for static domains. This will allow us to understand both which spatial modes result in instability and lead to possible patterning, but also the duration over which such instabilities persist. We shall later show that there is good agreement between full numerical simulations of the reaction-diffusion dynamics and the selection of spatial perturbations which are transiently unstable under our criterion.

3 Turing instability criteria on evolving domains

As pointed out in Section 2.5, the perturbation problem is non-autonomous and hence we can no longer rely on eigenvalues. Furthermore, it is transient dynamics, rather than dynamics as $t\rightarrow\infty$, which lead to patterning in reaction-diffusion systems on evolving domains. As such, we consider a comparison principle for the growth of a spatial perturbation, over some interval of time. This will allow us to determine when a perturbation grows with a certain rate and leads to a transient instability. As we shall later show with numerical simulations, these transiently unstable spatial modes do indeed correspond to patterns formed under the full nonlinear dynamics, subject to the restrictions of a linear analysis. Once a nonlinear pattern has developed, our results are no longer formally valid, though they can give some insight into pattern evolution as we will show later.

The long-time behavior of generic non-autonomous systems such as (13) are too complicated to consider in full generality (as even in the autonomous case, one would appeal to the Routh-Hurwitz stability criterion which becomes cumbersome for large systems [94]). In what follows, we restrict our attention to cases commonly arising in the literature. We consider the $n=2$ case for reaction-diffusion systems in detail, as this is the standard case considered in the literature for activator-inhibitor Turing systems. Still, if one is concerned with particular reaction kinetics with $n\geq 3$, then (13) can be solved numerically. We also consider systems which are higher-order in space. Even scalar systems of such a form can permit spatial patterning, and we consider both the scalar case as well as the case of two coupled equations. In all of these results, a dot over a quantity denotes a time derivative, two dots over a quantity denote the second time derivative.

3.1 Comparison principle

Growth rates for a scalar first-order non-autonomous problem are trivial to obtain, and we do not discuss this case. The corresponding second-order problem is not as simple, and we establish some growth bounds on general second-order non-autonomous ODE of the form

[TABLE]

There have been a variety of results for second-order non-autonomous ODE systems [30, 34, 40, 41]. Due to the breakdown of oscillating solutions, determining conditions for stability can be quite involved and can depend strongly on the properties of non-autonomous terms. On the other hand, obtaining conditions which are sufficient for instability can be viewed as somewhat easier. We begin with a result which gives sufficient conditions for a solution $Y(t)$ to (18) to grow on an arbitrary time interval, which we shall denote $\mathcal{I}\subseteq[0,\infty)$. In particular we can choose $\mathcal{I}$ unbounded to satisfy $|Y(t)|\rightarrow\infty$ as $t\rightarrow\infty$ at a prescribed rate of growth, or $\mathcal{I}$ bounded to only denote regions of transient instability. While such results will be sufficient rather than necessary for instability, as we shall later see, these results will provide the most natural generalization of the standard algebraic Turing instability conditions.

Theorem 3.1**.**

Let $\Phi\in C^{2}(\mathbb{R})$ such that $\Phi(t)>0$ for all $t\in\mathcal{I}$. Consider the ODE (18) and suppose that

[TABLE]

Then, (18) has a fundamental solution $Y(t)$ with $|Y(t)|\geq\Phi(t)$ for all $t\in\mathcal{I}$.

Proof.

We begin with the case where equality holds in the bound (19). For this case, one may verify that $\Phi(t)$ is in the fundamental solution set of (13) and hence for general initial data (13) has one solution satisfying $|Y(t)|=\Phi(t)$ for any $t\in\mathcal{I}$, although there may be a second solution which grows faster. Therefore $|Y(t)|\geq\Phi(t)$ for all $t\in\mathcal{I}$, with at least equality holding.

Next, consider $G(t)=-\frac{\ddot{\Phi}}{\Phi}-\frac{\dot{\Phi}}{\Phi}F(t)-H(t)$ for some $H(t)\geq 0$ for all $t\in\mathcal{I}$. Then, (18) takes the form

[TABLE]

There are two fundamental solutions to this equation, and any solution will be a linear combination of these.

We choose initial data $Y(t_{0})=\Phi(t_{0})$ and $\dot{Y}(t_{0})=\dot{\Phi}(t_{0})$, and make the change of variable $Y(t)=Y(t_{0})\exp\left(\int_{t_{0}}^{t}Z(s)ds\right)$, which puts (20) into the form of the Riccati equation

[TABLE]

Note that $Z(t_{0})=\dot{\Phi}(t_{0})/\Phi(t_{0})$, so we have

[TABLE]

where we define $Z_{1}(t)$ as a function satisfying

[TABLE]

with initial data $Z_{1}(t_{0})=\dot{\Phi}(t_{0})/\Phi(t_{0})$. One may verify that the exact solution reads $Z_{1}(t)=\dot{\Phi}(t)/\Phi(t)$. Now, by differential inequality (22) and since $Z(t_{0})=Z_{1}(t_{0})$, we have $Z(t)\geq Z_{1}(t)$ for all $t\in\mathcal{I}$. Integration and exponentiation preserve this ordering, and yield

[TABLE]

since $Y(t_{0})\exp\left(\int_{t_{0}}^{t}Z_{1}(s)ds\right)=Y(t_{0})\exp\left(\int_{t_{0}}^{t}\frac{\dot{\Phi}(s)}{\Phi(s)}ds\right)=\Phi(t)$. Then, for this choice of initial data, $|Y(t)|\geq\Phi(t)$ for all $t\in\mathcal{I}$. This completes the proof. $\blacksquare$ ∎

We note that these inequalities for $G(t)$ are all sufficient conditions for the prescribed time interval including large-time asymptotic behavior if $\mathcal{I}$ is unbounded. There may be specific problems for which these conditions are not necessary. Classifying such dynamics would involve an advanced study of oscillation theory, and we do not address this here, as our goal is to show that these kinds of sufficient conditions are consistent with the standard Turing conditions for static domains. In linear stability theory, one is often interested in the onset of exponential growth of a small perturbation, and for this case we have the following corollary:

Corollary 3.1**.**

Consider $\Phi(t)=\mu(t)\exp(\delta t)$ for some $\delta>0$ (where we scale with $\mu(t)$ since the factor of $\mu(t)^{-1}$ in (12) will moderate any instability). From Theorem 3.1, we have

[TABLE]

Taking $\delta\rightarrow 0^{+}$, and strict inequality, we recover the weakest bound for exponential growth over $t\in\mathcal{I}$,

[TABLE]

While exponential instabilities are the standard for discussing the Turing instability and related patterning (as well as any other stability or instability criteria dependent upon temporal eigenvalues), it is tempting to weaken the strength of the instability, in order to further probe the boundary of stability and instability regions. This difference is only notable for fixed-strength bounds with $\delta>0$ (say, when comparing an instability of rate $\exp(\delta t)$ with an instability of rate $t^{\delta}$), as taking the $\delta\rightarrow 0^{+}$ limit again results in the same bound (26), as seen the the following corollary:

Corollary 3.2**.**

Consider $\Phi(t)=\mu(t)t^{\delta}$ for some $\delta>0$. From Theorem 3.1, we have

[TABLE]

Taking $\delta\rightarrow 0^{+}$ to capture the weakest possible growth rate, and strict inequality, we recover the weakest bound for algebraic growth over $t\in\mathcal{I}$,

[TABLE]

Therefore, we conclude that the bound in equation (26) of Corollary 3.1 is robust in terms of accounting for the possible rates of instability.

In light of Theorem 3.1 along with Corollaries 3.1-3.2 which provide conditions granting a specific growth rate, it is worthwhile to obtain a complementary result on corresponding decay rates.

Theorem 3.2**.**

Consider the ODE (18) and suppose that

[TABLE]

Then, solutions $Y(t)$ to (18) do not grow faster than any exponential or algebraic function of $t\in\mathcal{I}$, namely $|Y(t)|\leq K\mu(t)e^{\delta t}$ for all $\delta>0$ and $t\in\mathcal{I}$ and similarly $|Y(t)|\leq K\mu(t)(1+t)^{\delta}$ where $K=\mathcal{O}(|Y(t_{0})|)$.

Proof.

First note that one can repeat the whole proof of Theorem 3.1 with the opposite inequality yielding that a solution, characterised by the initial condition $Y(t_{0})=\Phi(t_{0}),~{}\dot{Y}(t_{0})=\dot{\Phi}(t_{0})$, satisfies an upper bound $|Y(t)|\leq\Phi(t)$ and then the claim follows for this fundamental solution from the two corollaries 3.1, 3.2.

To finish the proof we need to show a similar estimate for the second fundamental solution of the ODE (18). We consider again a function $H\geq 0$ so that equality in the bound (29) is obtained

[TABLE]

and the same change of variable $Y(t)=Y(t_{0})\exp\left(\int_{t_{0}}^{t}Z(s)ds\right)$, which puts (20) into the form of the Riccati equation

[TABLE]

Note, however, that one cannot capture the second fundamental solution $Y_{2}(t)=\Phi(t)\int_{t_{0}}^{t}\frac{\exp(-\int_{t_{0}}^{\tau}F(s)ds)}{\Phi^{2}(t)}d\tau$ as it corresponds to initial data $Y_{2}(t_{0})=0$ being impossible to be captured by the aforementioned exponential change of variables. Hence we chose instead initial conditions $Y(t_{0})=\Phi(t_{0})$, $\dot{Y}(t_{0})=\frac{1}{\Phi(t_{0})}+\dot{\Phi}(t_{0})$ corresponding to the sum of the two mentioned fundamental solutions.

The initial condition after the change of variables reads $Z(t_{0})=\frac{\dot{Y}(t_{0})}{Y(t_{0})}=\frac{\dot{\Phi}(t_{0})\Phi(t_{0})+1}{\Phi(t_{0})^{2}}$, so we have

[TABLE]

where we define $Z_{2}(t)$ as a function satisfying

[TABLE]

with initial data $Z_{2}(t_{0})=\frac{\dot{Y}(t_{0})}{Y(t_{0})}=\frac{\dot{\Phi}(t_{0})\Phi(t_{0})+1}{\Phi(t_{0})^{2}}$. As we already know a solution to this Riccati equation we can identify a general solution to it

[TABLE]

while at $t_{0}$ its value is

[TABLE]

and hence to satisfy the prescribed initial condition we set $C=1$.

Finally, by differential inequality (31) and since $Z(t_{0})=Z_{2}(t_{0})$, we have $Z(t)\leq Z_{2}(t)$ for all $t\in\mathcal{I}$. Integration and exponentiation preserve this ordering, and yield

[TABLE]

where $\Psi(t)=\frac{1}{\Phi(t)^{2}}\exp\left(-\int_{t_{0}}^{t}F(\tau)d\tau\right)\geq 0$ since $Y(t_{0})=\Phi(t_{0})$, $\exp\left(\int_{t_{0}}^{t}Z_{2}(s)ds\right)=\frac{\Phi(t)}{\Phi(t_{0})}$, and as $\frac{\Psi(t)}{1+\int_{t_{0}}^{t}\Psi(\tau)d\tau}\leq 1$ due to nonnegativity of $\Psi$.

Then, for this choice of initial data, $|Y(t)|\leq e\Phi(t)$ for all $t\in\mathcal{I}$ and thus an arbitrary solution to the ODE (18) has to satisfy $|Y(t)|<K\Phi(t)$ for all $t\in\mathcal{I}$.

To complete the proof it suffices to choose $\Phi(t)=\mu(t)e^{\delta t}$ and $\Phi(t)=\mu(t)(1+t)^{\delta}$ followed by taking the limit $\delta\rightarrow 0^{+}$. $\blacksquare$ ∎

In Section 3.2, we will apply Theorem 3.1 and Corollary 3.1 in order to obtain the natural analogue of the Turing conditions for a system of two reaction-diffusion equations on an evolving domain. In light of the result presented in Corollary 3.2, the bound obtained in these results is sufficiently general to capture transient growth rates leading to instability. Furthermore, in light of Theorem 3.2, we do not expect a weaker bound to be useful, as the reverse strict inequality results in perturbations which are stable. Therefore, Theorem 3.1 and Corollary 3.1 indeed provide the most general bounds one is likely to obtain.

3.2 Instability conditions for systems of two reaction-diffusion equations

Due to the time variability of the base state and the actual growth, the nature of our stability result will be time dependent (rather than for $t\rightarrow\infty$ as is true of the classical Turing conditions), with modes losing and perhaps gaining stability over time. This is exactly along the lines of the history dependence observed in [45]. We shall then phrase the result in terms of a time interval over which the instability is observed. This interval becomes unbounded if the mode remains unstable as $t\rightarrow\infty$. Throughout the time interval on which an instability arises, given by $\mathcal{I}_{k}$ for each $\rho_{k}(t)$, we shall require $J_{12}\neq 0$ for all $t\in\mathcal{I}_{k}$. Otherwise, the equation for the first chemical species would decouple from the second, and either (i) the reaction kinetics would grow without bound for $J_{11}>0$ or (ii) the perturbation (12) can never give instability for $J_{11}<0$ for any arbitrary spatial perturbation, and hence pattern formation would be impossible. Hence, $J_{12}\neq 0$ is a reasonable assumption. Likewise, we shall assume $J_{21}\neq 0$.

We now apply the results of Theorem 3.1 and Corollary 3.1 to obtain conditions on the instability of spatial perturbations of the form (12). As we mentioned above, due to Theorem 3.2, these are the best instability bounds one expects to obtain. In order to invoke Theorem 3.1 and Corollary 3.1, we first convert the non-autonomous first-order system into a scalar non-autonomous second-order scalar ODE. The generalization of the Turing conditions is as follows:

Theorem 3.3**.**

Consider the evolution of a compact, simply connected, smooth Riemannian manifold $\Omega(t)\subset\mathbb{R}^{N}$ as in (4), with Laplace-Beltrami operator spectrum $\rho_{k}(t)\in C^{1}(\mathcal{I}_{k})$ where $\mathcal{I}_{k}\subseteq(0,\infty)$, such that volume expansion or contraction $\mu(t)\in C^{2}(\mathcal{I}_{k})$ given in (10) is independent of space. Assume that $J\in C^{1}(\mathcal{I}_{k})$ is the time-dependent Jacobian matrix of $\mathbf{f}$ evaluated at the spatially homogeneous solution $\mathbf{U}(t)$ to (11), with $J_{12},J_{21}\neq 0$ on $\mathcal{I}_{k}$. For $n=2$ species, the spatially homogeneous state $\mathbf{U}(t)$ for the reaction-diffusion system (3) is linearly unstable under a perturbation of the form (12) corresponding to $\rho_{k}(t)$ for $t\in\mathcal{I}_{k}$, provided that the inequality

[TABLE]

holds for all $t\in\mathcal{I}_{k}$.

Proof.

For $n=2$ species, and for each $k=0,1,2,\dots$, (13) reads

[TABLE]

Recall that $J=J(\mathbf{U})$, where $\mathbf{U}(t)$ is given by (11), hence the components of $J$ are in general time-dependent.

We start with $V_{1}(t)$. Since $J_{12}\neq 0$ for $t\in\mathcal{I}_{k}$, we isolate (33) for $V_{2}(t)$, and use it in (34) to obtain a single second-order ODE for $V_{1}(t)$, finding

[TABLE]

Applying (26) of Corollary 3.1 to (LABEL:v1b), we arrive at the sufficient condition

[TABLE]

which implies exponential growth of the $u_{1}$ component of the perturbation (12). We perform similar calculations using $J_{21}\neq 0$ for $t\in\mathcal{I}_{k}$, in order to obtain a second-order ODE for $V_{2}(t)$. Applying (26) of Corollary 3.1 to this ODE, we arrive at the sufficient condition

[TABLE]

which implies exponential growth of the $u_{2}$ component of the perturbation (12).

As we only require one of (LABEL:T31) or (LABEL:T32) to hold for instability, we take the inequality corresponding to the more extreme inequality in (LABEL:T31) or (LABEL:T32), resulting in the appearance of a max function in (LABEL:T3). This completes the proof. $\blacksquare$ ∎

These are the conditions for the system (6) to exhibit an instability corresponding to the $k$th spatial mode for $t\in\mathcal{I}_{k}$. In practice we shall consider $\mathcal{I}_{k}$ to be the largest interval on which the hypotheses of Theorem 3.3 hold, though for transient or sporadic growth periods there may be distinct intervals. These generalize the standard Turing conditions to corresponding conditions on a smoothly time-evolving manifold, though we remark that we do not yet incorporate a generalization of the standard stability of the homogeneous equilibrium in the absence of diffusion. Akin to what is done for classical Turing conditions, one may choose to group all modes which are unstable at time $t$, and the natural definition for this set will be: $\mathcal{K}_{t}=\{k\in\mathbb{N}|\{t\}\cap\mathcal{I}_{k}\neq\emptyset\}$. Similar generalizations hold when dealing with multi-indices. For higher-dimensional domains with spectra indexed like $\rho_{k_{1},\dots,k_{\ell}}$, we define $\mathcal{I}_{k_{1},\dots,k_{\ell}}$ and $\mathcal{K}_{t}$ accordingly.

Before continuing, we briefly connect our result to the standard Turing condition for instability on a static domain. We remark that in the case where the domain is static in time, the spectrum $\rho_{k}$ is constant, $\dot{\mu}=\ddot{\mu}=0$, and all entries in the matrix $J$ are constant. As such, the condition in Theorem 3.3 reduces to

[TABLE]

which is exactly the standard Turing condition on a static manifold. In the case where the manifold is flat and rectangular, or flat and unbounded, the spectrum $\rho_{k}=|\mathbf{k}|^{2}$ for some wavenumber vector $\mathbf{k}$ with dimension equal to the dimension of the space. For such a case, (38) reduces further to

[TABLE]

and this is the Turing condition most commonly seen in the literature [74].

3.3 Instability conditions for higher-order systems

Returning to the higher-order systems discussed in Section 2.4, we have the following analogue of Theorem 3.3 for coupled pairs ($n=2$) of systems taking the form (14):

Theorem 3.4**.**

Consider the evolution of a compact, simply connected, smooth Riemannian manifold $\Omega(t)\subset\mathbb{R}^{N}$ as in (4), with Laplace-Beltrami operator spectrum $\rho_{k}(t)\in C^{1}(\mathcal{I}_{k})$ where $\mathcal{I}_{k}\subseteq(0,\infty)$, such that volume expansion or contraction $\mu(t)\in C^{2}(\mathcal{I}_{k})$ given in (10) is independent of space. Assume that $J\in C^{1}(\mathcal{I}_{k})$ is the time-dependent Jacobian matrix of $\mathbf{f}$ evaluated at the spatially homogeneous solution $\mathbf{U}(t)$ to (11), with $J_{12},J_{21}\neq 0$ on $\mathcal{I}_{k}$. For $n=2$ species, the spatially homogeneous state $\mathbf{U}(t)$ for the system (14) is linearly unstable under a perturbation of the form (12) corresponding to $\rho_{k}(t)$ for $t\in\mathcal{I}_{k}$, provided that the inequality

[TABLE]

holds for all $t\in\mathcal{I}_{k}$.

The proof of Theorem 3.4 is similar to that of Theorem 3.3, so we omit it.

In the scalar case ($n=1$), it is also possible to have instability, provided that the degree of $\mathcal{P}$ is at least two. In the case where $\mathcal{P}(y)=y$, as for standard reaction-diffusion systems, the scalar reaction-diffusion system admits spatial perturbations (12) which grow or decay like

[TABLE]

Equation (41) is exactly solvable, and we have

[TABLE]

since $d>0$ and $\rho_{k}(t)\geq 0$. Hence, diffusion is always stabilising in the scalar case of one reaction-diffusion equation without any outside forcing. However, in the higher-order case, the structure of the polynomial $\mathcal{P}$ can permit instability in the scalar ($n=1$) case of (14),

[TABLE]

and we give conditions for such an instability in the following Theorem:

Theorem 3.5**.**

Consider the evolution of a compact, simply connected, smooth Riemannian manifold $\Omega(t)\subset\mathbb{R}^{N}$ as in (4), with Laplace-Beltrami operator spectrum $\rho_{k}(t)\in C^{1}(\mathcal{I}_{k})$ where $\mathcal{I}_{k}\subseteq(0,\infty)$, such that volume expansion or contraction $\mu(t)\in C^{2}(\mathcal{I}_{k})$ given in (10) is independent of space. Assume that $f^{\prime}(U)\in C^{1}(\mathcal{I}_{k})$, where $U(t)$ is the spatially uniform state satisfying (11) which in the scalar case reads $\frac{dU}{dt}=-\frac{\dot{\mu}}{\mu}U+f(U)$. For $\mathcal{P}$ of degree at least two, the spatially homogeneous state $U(t)$ for the system (43) is linearly unstable under a perturbation of the form (12) corresponding to $\rho_{k}(t)$ for $t\in\mathcal{I}_{k}$, provided that the inequality

[TABLE]

holds for all $t\in\mathcal{I}_{k}$.

Proof.

For the scalar equation (43), the linear stability of the $k$th spatial perturbation is given by the scalar form of (13),

[TABLE]

and solving (45) we obtain

[TABLE]

There is then growth of the perturbation (12) provided

[TABLE]

for some $\delta>0$, noting that we have assumed an exponential growth rate. (As pointed out in Section 3.1, an exponential growth rate is sufficiently general.) Rearranging, we find

[TABLE]

Both sides of the inequality are equal at $t=t_{0}$, so we may differentiate the inequality since the left hand side must grow faster than the right hand side in order for there to be a transient instability over $t\in\mathcal{I}_{k}$, and we find

[TABLE]

Taking strict inequality in the limit $\delta\rightarrow 0^{+}$, we have the stated inequality (44). This completes the proof. $\blacksquare$ ∎

We will give explicit examples of transient instabilities and pattern formation in fourth-order scalar systems in Section 7.

In Section 2.4, we also commented that the most generqal higher-order spectral problem may not simply involve Laplace-Beltrami eigenfunctions $\psi_{k}$ but also other eigenfunctions $\hat{\psi}_{k}$, depending on the higher-order boundary conditions. In this case, upon considering the stationary coordinates (8), one instead has the spectral problem $\mathcal{P}\left(\nabla_{\Omega(t)}^{2}\right)\hat{\psi}_{k}=-\zeta_{k}(t)\hat{\psi}_{k}$. Through a similar approach to Theorems 3.4-3.5, we find

Theorem 3.6**.**

Consider the evolution of a compact, simply connected, smooth Riemannian manifold $\Omega(t)\subset\mathbb{R}^{N}$ as in (4). Let the differential operator $\mathcal{P}\left(\nabla_{\Omega(t)}^{2}\right)$ have corresponding spectrum $\zeta_{k}(t)\in C^{1}(\mathcal{I}_{k})$ where $\mathcal{I}_{k}\subseteq(0,\infty)$, such that volume expansion or contraction $\mu(t)\in C^{2}(\mathcal{I}_{k})$ given in (10) is independent of space. Assume that $J\in C^{1}(\mathcal{I}_{k})$ is the time-dependent Jacobian matrix of $\mathbf{f}$ evaluated at the spatially homogeneous solution $\mathbf{U}(t)$ to (11), with $J_{12},J_{21}\neq 0$ on $\mathcal{I}_{k}$. Then, the spatially homogeneous state $U(t)$ for the system (43) corresponding to $n=1$ species is linearly unstable under a perturbation of the form (12) corresponding to $\zeta_{k}(t)$ for $t\in\mathcal{I}_{k}$, provided that the inequality

[TABLE]

holds for all $t\in\mathcal{I}_{k}$. Meanwhile, for $n=2$ species, the spatially homogeneous state $\mathbf{U}(t)$ for the system (14) is linearly unstable under a perturbation of the form (12) corresponding to $\rho_{k}(t)$ for $t\in\mathcal{I}_{k}$, provided that the inequality

[TABLE]

holds for all $t\in\mathcal{I}_{k}$.

This is the most general bound on such higher-order problems. However, for the sake of simulations, we restrict our attention to simple higher-order problems with boundary conditions yielding a straightforward collection of eigenfunctions. Still, Theorem 3.6 is more general than Theorems 3.4-3.5, and should be regarded as the primary result for such systems, as for the most general problems it will include spectral contributions that might be ignored by Theorems 3.4-3.5.

3.4 Asymptotic stability results for reaction-diffusion systems

In this section we shall analyse the special case of asymptotic stability, i.e. $\mathcal{I}_{k}=(T,\infty)$.

3.4.1 No pattern in the large-time limit for unbounded growth

Consider first the case of unbounded growth when $\rho_{k}(t)\rightarrow 0$ for all $k$ (which is the case for purely dilatational growth, $\rho_{k}(t)=\rho_{0}(0)/r^{2}(t)$). Then the condition (LABEL:T3) for instability of the $k$-th wavemode reads

[TABLE]

Note the obvious that in this case the condition for instability is independent of the wavemode number. Thus, if this condition is satisfied, then asymptotically all spatial modes lead to instability. Therefore we expect that such a situation would not yield a reasonable pattern (not having arbitrary small lengthscale) and further note that this inequality is also the one satisfied by $\rho_{0}(t)=0$ at finite times, indicating spatially homogeneous instability. Hence, for asymptotically large times, spatial instabilities (of arbitrary mode number) coincide with homogeneous instabilities which are typically precluded.

Therefore, there is no diffusion driven patterning on domains undergoing unbounded growth for asymptotically large time. Any spatial patterning under such growth must therefore occur due to transient dynamics. This is quite distinct from the static, bounded domain case, where diffusive instabilities retain their dominance in the large-time limit.

3.4.2 Saturating growth scenarios

If the growth stopped at a finite time or is saturating at a finite size then the above observations about the asymptotic stability reduce exactly to classical diffusion-driven instability condition. Indeed, in this case we have $\lim_{t\rightarrow\infty}\rho_{k}(t)=\overline{\rho}_{k}>0$, while $\lim_{t\rightarrow\infty}\dot{\mu}=0$ and $\lim_{t\rightarrow\infty}\dot{J}_{ij}=0$, so $\mathbf{U}^{*}$ is the homogeneous steady state solution satisfying $f(\mathbf{U}^{*})=0$ and the sufficient conditions read

[TABLE]

which is exactly the Turing condition which one would derive on the static domain $\overline{\Omega}$ [23] with $\mathcal{I}_{k}$ unbounded. As $d_{1}d_{2}\overline{\rho}_{k}^{2}$ is dominant for large $k$, the spatial modes of high frequency are stable, in line with classical Turing conditions [75].

3.4.3 Asymptotic stability of the base state

It is worth briefly discussing the asymptotic stability of the base state solution of (11), and to do so we consider two cases.

First, suppose that $\lim_{t\rightarrow\infty}\frac{\dot{\mu}}{\mu}=\nu$, a constant. This is true, for example, in the case of exponential growth of a domain. The dynamics of the base state (11) then read

[TABLE]

As this equation is autonomous in the large-time limit, to obtain a steady state we set $-\nu\mathbf{U}+\mathbf{f}(\mathbf{U})=\mathbf{0}$, and denote by $\mathbf{U}^{\dagger}$ a solution of this algebraic equation. Note that $\mathbf{U}^{\dagger}$ is not in general equal to $\mathbf{U}^{*}$, which is the solution of the algebraic equation $\mathbf{f}(\mathbf{U})=\mathbf{0}$ as discussed in Section 2.3. In particular, $\mathbf{U}^{\dagger}=\mathbf{U}^{*}$ when $\nu=0$, i.e., when the change in volume expansion or contraction is zero in the asymptotic limit $t\rightarrow\infty$. The important thing to note here is that in the case where $\lim_{t\rightarrow\infty}\frac{\dot{\mu}}{\mu}$ exists and is equal to some constant $\nu\in\mathbb{R}$, the appropriate spatially uniform base state is $\mathbf{U}^{\dagger}$ satisfying $-\nu\mathbf{U}+\mathbf{f}(\mathbf{U})=\mathbf{0}$ rather than the base state in the absence of domain evolution, $\mathbf{U}^{*}$, which satisfies $\mathbf{f}(\mathbf{U})=\mathbf{0}$.

Assuming a linear perturbation of the form

[TABLE]

equation (54) is linearized as

[TABLE]

where $J=\frac{\partial\mathbf{f}}{\partial\mathbf{U}}$ evaluated at $\mathbf{U}=\mathbf{U}^{\dagger}$. In particular, the Jacobian matrix is now constant as it is evaluated at the steady state $\mathbf{U}^{\dagger}$. In the standard way, for the $n=2$ species case, we have from (56) that the steady state $\mathbf{U}^{\dagger}$ is asymptotically stable provided that the following necessary conditions are satisfied:

[TABLE]

These are natural analogues of the standard necessary conditions for reaction kinetics to be stable on static domains, $\text{tr}(J)<0$ and $\det(J)>0$.

Of course, for more complicated domain evolution, the limit $\lim_{t\rightarrow\infty}\frac{\dot{\mu}}{\mu}$ need not exist, and this leads us to our second case. In this case, we perturb a time-dependent solution of (11) akin to what we did in (55). However, this is equivalent to a perturbation of the form (12) with a spatially homogeneous mode (the $k=0$ mode, which always exists for the Neumann problem on manifolds with boundary, as well as for manifolds without boundary). In light of Theorem 3.2, we anticipate a condition akin to that given in Theorem 3.3 only with a sign reversed. Indeed, carrying out a similar analysis, and invoking Theorem 3.2, we find that a necessary condition for the stability of a spatially uniform yet time varying base state solution to (11) reads:

[TABLE]

This is the complement of the condition given in Theorem 3.3, for the zeroth mode $k=0$ with $\rho_{0}(t)=0$. This condition is also complementary to that given in (52), which makes sense, as (52) was the unrealistic sufficient condition for all modes (even the zeroth mode) to result in an instability.

3.5 Equal diffusion coefficients

To explore whether equal diffusion coefficients permit pattern formation it is advantageous to transform the evolution equations for the perturbation (13) via

[TABLE]

into

[TABLE]

Note that for finite time intervals or when $\rho_{k}(t)\rightarrow 0$ the transformation is such that $\mathbf{A}$ is stable iff $\mathbf{V}$ is stable. Otherwise, i.e. for finite positive limit $\lim_{t\rightarrow\infty}\rho_{k}(t)=\rho_{k}^{*}>0$, stability in the original variable $\mathbf{V}$ is guaranteed only if $\mathbf{A}$ grows at least as fast as the matrix exponential $\exp(\rho_{k}^{*}{D}t)$.

To use Theorem 3.1 and Corollary 3.1, we rewrite the alternative relation for perturbation evolution (60) as a second-order equation

[TABLE]

where the equation for $A_{2}$ is the same with swapped indices $1\leftrightarrow 2$.

In the case where all diffusion coefficients are equal, $D=d_{1}I$, hence $\exp\left(P_{k}(t)D\right)=e^{d_{1}P_{k}(t)}I$, and we have $\dot{\mathbf{A}}=J(\mathbf{U})\mathbf{A}$. Then, $\mathbf{V}_{k}=\exp\left(-P_{k}(t)D\right)\mathbf{A}=e^{-d_{1}P_{k}(t)}\mathbf{A}$, where $\mathbf{A}$ depends on reaction kinetics at the base state through $J(\mathbf{U})$. As $d_{1}>0$ and $P_{k}(t)>0$, the contribution of diffusion is stabilizing, with any instability arising only from a combination of domain evolution (through the $\mu(t)$ term in (13)) and reaction kinetics, precluding spatial patterning due to diffusive instabilities.

In addition, an application of Corollary 3.1 reveals that no instability (pattern) can be expected for equal diffusion coefficients for finite time intervals or for unbounded growth with $\rho_{k}(t)\rightarrow 0$ as $\rho_{k}(t)$ vanishes from both $F(t;k)$ and $G(t;k)$. Finally, in the asymptotic case with $\lim_{t\rightarrow\infty}\rho_{k}(t)=\overline{\rho_{k}}$ the functions $F,~{}G$ are also independent of $k$, however, the requirement of exponential growth at least as fast as the matrix exponential $\exp(\rho_{k}^{*}{D}t)$ results in a dependence of the upper and lower bounds on $k$ only via the exponential bound $\delta=d_{j}\overline{\rho_{k}}$. Because both terms in the bounds depending on wave number $k$ are negative, the boundary between the upper and lower bound (hence the threshold for instability) is more stringent than a sufficient condition for instability without diffusion. Therefore equal diffusion coefficients do not allow the emergence of spatial patterns (not only being diffusion driven) even on evolving domains.

3.6 Transient breakdown of the continuum assumption

As noted above, any understanding of transient dynamics is welcomed, especially on evolving domains. The sufficient condition (the threshold for instability) given in Theorem 3.3 also allows us to study such effects. Indeed, while [45] explored transient growth, and showed that rapid growth can lead to arbitrarily-large wavemode excitation, their results can be found as a special case of the instability results presented here.

We shall use the above transformation (59). From Corollary 3.1, transient exponential growth then happens when

[TABLE]

where $\neg j$ denotes the index not being $j$ (e.g., if $j=1$ then $\neg j=2$). Focusing on large wavenumbers $k$, we find that they become unstable if

[TABLE]

or equivalently, when

[TABLE]

One immediately observes that for sufficiently rapid growth, this inequality is satisfied (for one of the $j$, since we assume $d_{1}\neq d_{2}$), and hence fast growth (measured by $\dot{\mu}/\mu$) always yields transient exponential growth for arbitrarily large wavenumbers, provided that the time-dependent Jacobian entries remain bounded. Such an instability entails a breakdown of the linear analysis as exemplified in [45]. On the other hand, if the Jacobian entries also change rapidly (recall that they depend on the dynamics of the base state governed by (11), and hence on the quantity $\dot{\mu}/\mu$), then this effect will be suppressed. In practice, we do not observe transient exponential growth for arbitrarily large wavenumbers in our simulations, at any time.

4 Applications to reaction-diffusion systems in one space dimension

We illustrate the analytical instability conditions given in Theorem 3.3 by considering various case studies consisting of specific growth functions and domain geometries, some of which extend studies in the literature, and others of which have seemingly never been considered due to their complexity in the face of existing methods. We note that the conditions given in Theorem 3.3 are sufficient for an instability to grow over a specified time interval, but are insufficient to determine if a given time period is sufficient to observe a heterogeneity forming in a simulation of the full system, as this will depend on the specific nonlinearities involved. Nevertheless, we aim to show that the linear stability analysis captures a variety of solution features observed in numerical simulations.

We first consider systems in one spatial dimension in the present section, before moving onto more complicated configurations in the following sections. Before this, we provide a brief discussion of the reaction kinetics and numerical schemes used.

4.1 Reaction kinetics

We consider the Schnakenberg, or activator-depleted, reaction kinetics as a very simple example which is used extensively in the literature [32, 95]. The kinetics and homogeneous equilibria at $t=t_{0}$ read

[TABLE]

where $a,b$ will be taken as non-negative real parameters. We will also consider the FitzHugh-Nagumo kinetics to demonstrate the applicability of our results to an oscillatory base state giving rise to Hopf and Turing-Hopf bifurcations [29, 42, 76]. The kinetics and homogeneous equilibria at $t=t_{0}$ read

[TABLE]

where $a,b,c$ and $i_{0}$ are taken as non-negative constants, and $U_{1}^{*}$ will be the root of $c(U_{1}^{*}-U_{1}^{*3}/3-(a-U_{1}^{*})-i_{0})=0$. For the parameters we will use, this equation will have a unique real root, and so the system will have a unique steady state solution.

4.2 Numerical approach

We simulate (6) with the kinetics (65) using the finite element solver COMSOL, version 5.3, with which we discretize the manifolds using second-order finite elements (which will be triangular and tetrahedral in the higher dimensional examples). We used Matlab to compute the evolution of the homogeneous state, and to generate $\mathcal{I}_{k}$ according to Theorem 3.3. We verified simulations in various static domain cases (1-D intervals, 2-D spheres) using the Matlab package Chebfun [110], in addition to convergence checks in spatial and time discretizations. In all simulations, we used a relative tolerance of $10^{-5}$, and fixed an initial time step of $10^{-6}$. We used COMSOL’s default backward difference formula method which then adaptively updated the time step beyond this initial value. In 1-D we used $10^{4}$ finite elements, and for higher-dimensional simulations used at least $10^{4}$ elements, though this varied for each geometry. Some restrictions were used on the maximum allowable time step to prevent behaviors such as the loss of modes in the initial perturbations. Convergence in time was checked by restricting the maximum time step, and convergence in space was determined via computing solutions across varied numbers of finite elements, and comparing the norm of solutions over time and space.

We emphasize that in the cases of fast or non-monotonic domain growth, extreme care is needed due to the non-autonomous nature of the spatial operator. We note that one advantage of this choice of finite element software, as well as the restriction to dilational growth, is that it allows for simple implementations of growing manifolds where the growth is directed in particular directions in the ambient space. This is because the Laplace-Beltrami operator on a surface of dimension $N$ can be constructed from the Laplace operator in the ambient space $\mathbb{R}^{N+1}$, so that dilation of a particular coordinate in $\mathbb{R}^{N+1}$ allows a natural construction of the time-dependent Laplace-Beltrami operator on the surface. We note that there exist many other choices for numerical methods for such problems [5, 61].

Initial data is taken to be of the form $\mathbf{u}(0,\mathbf{x})=(I+{\zeta}(\mathbf{x}))\mathbf{U^{*}}$, where $I$ is the identity matrix and ${\zeta}=\text{diag}(\zeta_{1},\zeta_{2})$ are normally distributed perturbations which are independent across space and for each morphogen. Specifically, for each $\mathbf{x}\in\Omega^{*}$ and $i=1,2$, we take $\zeta_{i}(\mathbf{x})\sim\mathcal{N}(0,10^{-1})$. We have also considered smaller initial perturbations for each case, and note that whether or not a pattern persists despite transient periods of growth and decay is highly dependent on the size of the perturbation. For this reason, we use this reasonably large perturbation for all simulations, as the finite-time effects we study are intrinsically linked to observing growth of finite perturbations. For each geometry we show simulations using the same realization of the initial data throughout, though for a given size of perturbation (the variance of $\zeta$), we observe qualitatively similar dynamics for different realizations.

We consider two relevant sets to help visualize our instability criterion. We will consider these sets as functions of time. The first is a generalization of a time-dependent generalized Turing space which is the set of all parameters for which Theorem 3.3 predicts an instability for some $k\geq 0$. Here we will consider as an example the non-negative parameters $(a,b)\in\mathbb{R}^{2}_{+}$ for the kinetics given by (65), but of course generalizing these definitions is straightforward. We then define such a space, for a given time $t$ as: $\mathcal{T}_{t}=\left\{(a,b)\in\mathbb{R}^{2}_{+}|\cup_{k\geq 0}(\{t\}\cap\mathcal{I}_{k})\neq\emptyset\right\}$. Of course one can generalize $\mathcal{T}_{t}$ to a set of times, say $\mathcal{S}([t_{1},t_{2}])=\cup_{t_{1}\leq t\leq t_{2}}\mathcal{T}_{t}$, rather than the singleton time, but for our purposes we prefer to think of these as sets parameterized by time. We separately plot the space corresponding to homogeneous instabilities, which are times $t\in\mathcal{I}_{0}$, so that one may consider Turing spaces which exclude these points. Similarly, for fixed parameters, we may be interested in plotting an analogue of the classical dispersion relation which indicates which wavenumbers $k$ are excited as a function of time $t$. We define this dispersion set to be: $\mathcal{K}_{t}=\{k\in\mathbb{N}|\{t\}\cap\mathcal{I}_{k}\neq\emptyset\}$.

We will compare these time-dependent sets to the quasi-static Turing space and dispersion relations. These are given by ignoring the non-autonomous nature of the system, and treating the domain length as a parameter in the classical static Turing conditions. While such quasi-static conditions are not formally valid, we will demonstrate cases where they do seem to capture the qualitative behaviour of the system, in addition to cases where they fail. We will also compare the observed modes from full numerical solutions, deduced via the Fast Fourier transform, with predictions from our linear stability theory in the 1-D setting to provide evidence of the applicability of our analysis. We will only plot the single Fourier mode with the largest absolute power, corresponding to the wavemode with the largest component of an expansion of the full spatial solution. If the variation of the solution across the domain is less than $1\%$ of its mean value, then we set the largest mode to $k=0$, essentially neglecting small variations from the homogeneous base state.

Finally, we again reiterate that the instability criterion given by Theorem 3.3 only tells us if the $k$th mode is growing or not at some time, but not directly the growth rate (though a bound on this can be inferred via Corollary 3.1 which we will use to identify the fastest growing mode at a given time). Additionally, analysis of nonlinearities is necessary to determine conditions for whether or not a pattern fully develops and persists, or undergoes subsequent instabilities, such as peak-splitting. Nevertheless, we have exhaustively explored this condition numerically and confirmed that patterns typically develop if the parameter set is within the Turing space for a sufficiently long time, or equivalently that at least one mode remains unstable for a sufficient period.

4.3 Isotropic evolution of a line segment

The simplest and most commonly studied example in the literature is a uniformly growing line segment. We define $\Omega(t)\subset\mathbb{R}$ by $\Omega(t)=[0,r(t)]$. The moving coordinate is $X=r(t)x$, for $x\in[0,1]$, and we find $\rho_{k}(t)=\frac{\pi^{2}k^{2}}{r(t)^{2}}$ and $\mu(t)=r(t)$. We will use this simple geometric setting to explore various Turing spaces and dispersion relations for a variety of growth functions $r(t)$, to demonstrate how the instability regions change, particularly away from the well-studied case of slow growth. Our main aim is to show that the instability criterion in Theorem 3.3 can effectively capture instabilities in this time-dependent setting, and how it differs radically from either quasi–static approaches [117], or the small corrections due to slow growth previously reported in the literature [45, 63].

Under the criterion given in Theorem 3.3, we find that the $k$th perturbation of the form (12) corresponding to $\rho_{k}(t)=\frac{\pi^{2}k^{2}}{r(t)^{2}}$ is unstable over some interval $t\in\mathcal{I}_{k}$, provided that the inequality

[TABLE]

holds for all $t\in\mathcal{I}_{k}$. In the special case where $r(t)=L$ for constant $L>0$, hence the domain is static, the condition (LABEL:excond1d) reduces to

[TABLE]

which is exactly the classical Turing condition for a static one dimensional domain $[0,L]$.

A specific form of growth which is somewhat popular in the literature is exponential growth, which takes the form $r(t)=r(0)\exp(st)$, $s>0$. In addition to biological plausibility, another reason for the popularity of exponential growth is that it allows for the volume expansion term to take the form $\frac{\dot{\mu}}{\mu}=Ns$ (where $N$ is the dimension of the space domain), a constant, which greatly simplifies the dynamics of the spatially uniform system. Exponential isotropic growth of surfaces in $\mathbb{R}^{3}$ was extensively studied in [109], albeit under the assumption of a time-independent base state for (11).

We choose parameters of the kinetics (65) and an initial domain of size $r(0)=10$ for which the system would be on the boundary of the Turing space for a static domain, only admitting a single unstable wavenumber $k=1$. We then simulate (6) until the domain has grown to $r(t_{f})=30r(0)$. We show our results in Fig. 1. In each row, we plot solutions to the uniform base state $\mathbf{U}$ from Equation (11) in the first column, the PDE solution $u_{1}$ in the second column, and the dispersion set $\mathcal{K}_{t}$ in the third column, with each row demonstrating an increasing growth rate. We observe that the dynamics of the uniform base state plays a substantial role in determining both $\mathcal{K}_{t}$, and consequently the evolution of the pattern.

As exponential growth leads to an autonomous planar system, we observe that the decaying oscillations in Fig. 1(b)i are due to a stable spiral, and that the oscillations in Fig. 1(c)i are due to a Hopf bifurcation which has created a stable limit cycle. These decaying and persistent oscillations have an impact on the timescale over which a pattern can emerge, and we only see the onset of a pattern near the end of the simulation time in Fig. 1(c)ii. The fastest growth rate results in a uniform base state which grows far from the original kinetic equilibrium, and pattern formation is no longer possible. As the set of unstable wavenumbers grows exponentially, there is a hysteresis effect such that if a perturbation has not left the base state sufficiently early on, then a pattern cannot form, whereas a developed pattern persists. The quasi-static Turing space is identical to that shown in Fig. 1(a)iii, and due to the choice of the growth, is independent of the growth rate (up to relabelling time). Hence the qualitative differences in the third column are all manifestations of the non-autonomous nature of the growth.

The largest modal components observed (in blue) roughly follow a peak-splitting mode doubling process, which is most apparent for the slowest growth case in Fig. 1(a)ii, but breaks down for faster growth as in Fig. 1(c)ii, as anticipated by [114]. While the linear analysis does not precisely predict these observed modes, it does give a qualitative insight into the processes leading to these patterned states. Specifically, the numerically observed modes all follow a period of time wherein that specific mode has been unstable and allowed to grow away from the homogeneous base state. Additionally, the lower-frequency solutions seen in the faster growing domains can also be explained as, due to oscillations, the system does not remain in a state admitting a given unstable mode for nearly as long as it does for the slower growth cases.

Next we consider Turing spaces, $\mathcal{T}_{t}$, at different instances in time, in Fig. 2. The first column shows the initial Turing space, which is equivalent to the quasi-static space obtained by just incorporating the growth rate into the kinetics [63], and specifically (i) is equivalent to the static Turing space without growth. As expected, we observe little change from a small kinetic addition at $t=0$, but for larger times we see previously-unstable regions become stable, and regions becoming unstable to homogeneous perturbations, as well as new regions becoming unstable as the Turing space expands around the edges. Such observations are in line with the results of [45], though we remark that these spaces are not equivalent as our approach accounts for discrete wavenumbers, and does not need the assumption of slow growth.

We consider linear growth in Fig. 3, with increasing growth rates in each subsequent row. Other than similar transient effects to before, the final modes observed are similar in each case except as the growth rate surpasses $s=0.16$. Slightly beyond this point, by $s=0.165$, the steady state of the uniform base states is no longer stable, and instead we see in Fig. 3(c)i that $U_{1}$ tends toward [math], and $U_{2}$ diverges to infinity. We remark that this destabilization of the uniform base state’s long-time behavior can be observed in both the dispersion sets and space-time plots. The concentration of $u_{2}$ increases in time uniformly as the domain expands. This phenomenon is inherently non-autonomous, and depends strongly on the initial condition; for other choices of $\boldmath{U}(0)$ we observe different behaviors. In Fig. 3(b)ii,iii, we see sharp oscillations with increasing amplitudes before a pattern is allowed to form, suggesting a kind of excitability inherent in the transient dynamics.

There are a wide variety of more complex kinds of domain evolution one could consider, especially if we allow expansion and contraction rather than monotonic growth. As a simple example of this, we consider periodic growth and contraction given by a sinusoidal function in Fig. 4. Again, the set $\mathcal{K}_{t}$ in the slow case (a)iii is identical to the quasi-static approximation, and the observed modes follow this reasonably well. While the dispersion set only has extremely small changes in (b)iii, we see that the base states in this case slowly oscillate in (b)i, and that the pattern seemingly disappears during the height of contractions in (b)ii, only to reappear later. In the case of more rapid oscillations, the uniform base states oscillate irregularly, and spatial pattern formation is only intermittent (see $t\in[25,40]$) and fails to persist. One key observation that is clear from the plots of $\mathcal{K}_{t}$, is that contraction of the domain is a highly stabilizing effect, as during the contracting period of Figs. 4(b,c)iii, we see substantially fewer unstable modes than during the expanding phase.

While growth has been heavily studied in the reaction-diffusion literature, contraction or other complex forms of domain evolution have not been as thoroughly explored. While these periodically expanding and contracting domains may be somewhat exaggerated from realistic examples, we remark that large contractions have been observed in the blastocyst stage of mice embryos [98]. Such contractions can lead to a decrease in volume of as much as $20\%$, and potentially play important roles in morphogenesis, likely altering local chemical concentrations in addition to mechanical effects. While most biological media are undergoing expansion and growth, we anticipate that a more nuanced and accurate representation of morphogenesis will necessarily involve processes such as contraction. Indeed, if the domain contraction is strong enough, it forces the solution to be spatially heterogeneous (yet still oscillatory in time) and can suppress future patterning. We comment further on this point later.

In many of our plots and simulations, the blue line corresponding to the maximal mode in the FFT lies mostly within the shaded instability region. In other cases, particularly those where the evolution of the instability region is not strictly monotone in time, there is a slight lag in the full nonlinear system in responding to instabilities, which is why the maximal mode will sometimes extend outside of the black region. It is important to note that, at the onset of instability, the maximal mode resulting from the instability lies in the shaded region, since the new pattern is set by the instability. This maximal mode can then persist for a time even as the instability region shifts, due to nonlinear terms stabilizing the fully nonlinear simulations. However, the maximal mode gradually loses stability, and a new dominant unstable mode is selected within the shaded instability region which is valid at that time. This process continues.

For more extreme cases, such as that shown in Fig. 4(c), there is a strong contraction of the domain leading to all modes becoming stable. This stabilizes a uniform solution, and upon the later expansion of the domain, the concentration remains uniform yet oscillates. The reason for this is that upon the second expansion of the domain, the uniform solution is not acted on by a spatial perturbation (as it was going into the first expansion). Without a spatial perturbation, the higher spatial modes are never activated, despite the domain change, and there is hence no Turing instability. This is in particular seen in Fig. 4(c)(ii), where after the first domain contraction the later dynamics are spatially uniform and simply oscillate in time as the density of the chemicals change along with the domain size. Therefore, it appears as those through solutions maintaining some spatial heterogeneity over time are susceptible to later bifurcations leading to new spatial patterns (such as peak splitting leading to a doubling of localized structures at each bifurcation), as the spatial variations provide enough noise to permit successive Turing instabilities as the domain grows. However, those solutions for which domain evolution suppresses spatial heterogeneity resulting in a spatially uniform state, subsequent dynamics associated to domain evolution do not appear sufficient to initiate later Turing instabilities. While our linear instability analysis compares well with full numerical simulations in this regard, a more rigorous nonlinear analysis focused specifically on this behaviour would possibly elucidate this suppression of pattern formation.

4.4 Isotropic evolution of an excitable medium

We now consider the reaction kinetics (66) with parameters corresponding to the Turing (but not Turing–Hopf) space for a static domain (see [92] for bifurcation diagrams). We consider linear and step wise growth functions to demonstrate the impact that an excitable system has on pattern formation. Theorem 3.3 is also useful in determining when spatial modes can destabilize a homogeneous but oscillating base state on a static domain, such as that which occurs generically when the kinetics have undergone a Hopf bifurcation. We remark that linear analysis is insufficient to completely characterize instabilities which involve competition between both unstable Turing and Hopf modes, and generally the behavior can depend on the initial perturbation in addition to the parameters. Nevertheless, we demonstrate here that Theorem 3.3 can give some insight into when these Hopf modes can occur, which is a prerequisite to both purely oscillatory or spatiotemporal dynamics involving the competition of modes from both kinds of instabilities. Additionally we demonstrate how the solution to Equation (11) precisely determines the possibility of oscillatory dynamics.

In Fig. 5 we consider the linear growth case. For very small growth rates we recover the quasi-static dispersion relation (not shown), but as the growth rate is increased we observe transient oscillations as the base state slowly spirals back to its steady state value (Fig. 5(a)). As the growth rate is increased further, the initial disturbance from the kinetic steady state leads to a sustained oscillation (Fig. 5(b)i), which persists even when the growth is no longer substantially influencing the dynamics. The oscillatory base state leads to a dispersion set which is no longer a simply connected set, such that modes oscillate between growing for some time and decaying for others, which prevents the formation of spatial patterns. This occurs because even without growth, the base state dynamics are excitable such that both a stable steady state and a stable limit cycle coexist for these parameters, and growth provides the necessary perturbation to transition between the two attracting sets. We do note that over longer time periods, spatiotemporal patterns appear to form as indicated by the FFT results in Fig. 5(b)iii. There structures do not persist for long, however, with the patterns forming, dissipating, and then forming again, due to the strongly oscillatory state and resulting non-monotone instability region.

Similarly, in Fig. 6 we observe that a short but rapid domain expansion can induce the same type of multistability. If the increase in the size of the domain is sufficiently slow, a connected dispersion set is recovered. In fact, the quasi-static approach would always generate such a continuous set, as it cannot account for the possibility of an oscillatory base state. While stepwise growth is less simple to analyze than that of linear or exponential growth, it has physiological significance in a number of organisms which exhibit pulsatile growth spurts between periods of slow or stagnant growth during development [6, 28] which we model by a rapid smooth expansion.

Fig. 6(b)iii again shows that an oscillatory base state can result in transient patterns rather than a single persistent patter, like what was seen in previous examples. Fig. 6(a)iii, however, shows something new and fairly interesting. The initial pattern (corresponding to a dominant mode of $k\sim 7$) destabilizes, with a new pattern selected during the short but rapid growth near $k\sim 12$. After this, this pattern is locked in, even though the instability region soon after becomes fixed between $14\leq k\leq 20$ once the growth of the domain finishes. Since the pattern was formed during this short growth period, it lies adjacent to but just outside of the instability region for $t\rightarrow\infty$. This highlights an interesting case where the Turing pattern would not have been detected in the asymptotic limit of $t\rightarrow\infty$, even though the dynamics of the perturbation become autonomous in this limit. This again highlights the strong role of hysteresis in forming patterns when domain evolution is involved.

5 Applications to reaction-diffusion systems on isotropically growing manifolds

We next give examples of domains which evolve in more than one spatial dimension.

5.1 Isotropic evolution of a circular disk in $\mathbb{R}^{2}$

Turing conditions for reaction-diffusion systems on a static disk were recently obtained in [93], and taking growth and volume expansion terms to zero we recover their conditions as a special case. Numerical simulations and experimental results for a specific application of Lengyel-Epstein reaction kinetics on a growing disk were given in [85], although external forcing on the reaction-diffusion model was employed. More recent experiments in a radially expanding domain have been performed which show a crucial mode selection phenomenon induced by the speed of the growth [49].

We consider $\mathbf{X}=(r(t)x_{1}\cos(2\pi x_{2}),r(t)x_{1}\sin(2\pi x_{2}))$, $x_{1},x_{2}\in[0,1]$, for the isotropic evolution of a circular disk. Given Neumann data on the circle $|\mathbf{X}|=r(t)$, we have $\rho_{\ell,k}(t)=\frac{\tilde{j}_{\ell,k}^{2}}{r(t)^{2}}$, where $\tilde{j}_{\ell,k}$ denotes the $k$th positive root of the derivative of the Bessel function of the first kind $\mathcal{J}_{\ell}(x)$, $\ell=0,1,2,\dots$, and $\mu(t)=r(t)^{2}$.

Under the criterion given in Theorem 3.3, we find that the $(\ell,k)$th perturbation of the form (12) corresponding to $\rho_{\ell,k}(t)=\frac{\tilde{j}_{\ell,k}^{2}}{r(t)^{2}}$ is unstable over some interval $t\in\mathcal{I}_{\ell,k}$, provided that the inequality

[TABLE]

holds for all $t\in\mathcal{I}_{\ell,k}$. In the special case where $r(t)=R$ for constant $R>0$, the condition (LABEL:exdisk) reduces to

[TABLE]

which is exactly the classical Turing condition for a static circular disk of radius $R$.

We show simulations on a growing disk in Fig. 7. We see that for slower growth rates, the successive instabilities lead to symmetric patterning, though as the growth rate is increased additional unstable modes lead to less-robust patterning, with additional irregularity in the pattern structure. This is analogous to the one-dimensional case where robustness is only attained for certain growth rates [19, 114]. For larger growth rates, there is no Turing pattern and $u_{2}$ grows uniformly while $u_{1}$ decays to zero everywhere. We note that due to the difference in $\mu(t)$ in the setting of a two-dimensional manifold, the specific value at which this instability in the base state occurs is different from the one-dimensional setting shown in Fig. 3. Regarding the dispersion plots, we select a suitably large subset of them and then sort the resulting eigenvalues $\tilde{j}_{\ell,k}$ by magnitude to obtain a one-dimensional dispersion set analogous to those shown in the previous section, though omit any Fourier analysis of the numerical simulations or maximal growth rates. As anticipated, we observe broadly similar curves for different growth rates, though the faster rate admits transient oscillations when the domain is small, as seen earlier for one-dimensional cases.

5.2 Isotropic evolution of an equilateral triangle

We found no studies of Turing patterns on equilaterial triangles (static or growing), yet this is an example for which our results can be easily applied. Consider an equilateral triangle defined by $\Omega(t)=\left\{(r(t)x_{1},r(t)x_{2})~{}|~{}0\leq x_{1}\leq 1,~{}0\leq x_{2}\leq\sqrt{3}x_{1},~{}x_{2}\leq\sqrt{3}(1-x_{1})\right\}$. Due to self-similarity of this domain, growth of a single face is equivalent to growth of all three faces, with the domain remaining an equilateral triangle for all $t\geq 0$. The spectrum for $\Omega(0)$ is given by $\rho_{k_{1},k_{2}}(0)=\frac{16\pi^{2}}{9}\left(k_{1}^{2}+k_{1}k_{2}+k_{2}^{2}\right)$, for $(k_{1},k_{2})\in\mathbb{N}^{2}$, with the spectrum on the evolving domain then given by $\rho_{k_{1},k_{2}}(t)=\frac{16\pi^{2}}{9r(t)^{2}}\left(k_{1}^{2}+k_{1}k_{2}+k_{2}^{2}\right)$. As the domain is flat and planar, we have $\mu(t)=r(t)^{2}$.

Under the criterion given in Theorem 3.3, we find that the $(k_{1},k_{2})$th perturbation of the form (12) corresponding to $\rho_{k_{1},k_{2}}(t)$ is unstable over some interval $t\in\mathcal{I}_{k_{1},k_{2}}$, provided that the inequality

[TABLE]

holds for all $t\in\mathcal{I}_{k}$. In the special case of a static equilateral triangle with constant area $A>0$, we take $r(t)=\sqrt{\frac{4\sqrt{3}}{3}A}$, and the condition (LABEL:extriangle) reduces to

[TABLE]

which is the classical Turing condition for a static equilateral triangle of area $A$.

We consider linear isotropic growth in Fig. 8. We observe a more ordered formation of stripes in the case of slower growth, whereas spots and disordered connections appear for faster growth. We note that the quasi-static dispersion sets are identical to Fig. 8(a)iii, but the inset in Fig. 8(b)iii demonstrates the impact of a faster growth rate.

We also demonstrate the time dependence of sets $\mathcal{K}_{t}$ for more than one index in Fig. 9, where we compare results for the disk and triangle. We observe an increasing band of unstable wavenumbers emanating from the origin. Although dispersion sets $\mathcal{K}_{t}$ agree between both domains in a qualitative sense, we see that the unstable modes in the triangular case are not bounded by lines as in the case of the disk, but instead by circular arcs, due to the form of the Laplace-Beltrami eigenvalues for each respective domain.

5.3 Isotropic evolution of a sphere $S^{N}\subset\mathbb{R}^{N+1}$

There have been various pattern formation studies on static 2-spheres [12, 51, 89, 117]. An exponentially growing 2-sphere was considered by [33], although their analysis was quasi-static, thereby ignoring the role of transients in the dynamics of (11). Similar assumptions were made in [108]. Numerical simulations of pattern formation on growing 2-spheres, as well as anisotropic growth of 2-spheres into ellipsoids, were obtained in [52].

The unit $N$-sphere ${S}^{N}\subset\mathbb{R}^{N+1}$ has spectrum $\rho_{k}(0)=k(k+N-1)$, $k=0,1,2,\dots$, hence $\rho_{k}(t)=\frac{k(k+N-1)}{r(t)^{2}}$ and $\mu(t)=r(t)^{N}$. These eigenvalues will have increasingly large multiplicity corresponding to different eigenfunctions on ${S}^{N}$. That is to say, for a fixed $k$, there can exist multiple distinct $\psi_{k}(\mathbf{x})$ in the general perturbation (12). As such, each of these distinct spatial eigenfunctions can yield patterning when the perturbation corresponding to $\rho_{k}(t)$ is unstable.

Under the criterion given in Theorem 3.3, we find that the $k$th perturbation of the form (12) corresponding to $\rho_{k}(t)=\frac{k(k+N-1)}{r(t)^{2}}$ is unstable over some interval $t\in\mathcal{I}_{k}$, provided that the inequality

[TABLE]

holds for all $t\in\mathcal{I}_{k}$. In the special case where $r(t)=R$ for constant $R>0$, hence the domain is a static $N$-sphere of radius $R$, the condition (LABEL:exsphere) reduces to

[TABLE]

which is exactly the classical Turing condition for a static $N$-sphere.

We first obtain solutions on $S^{2}$ in Fig. 10 for different rates of linear growth. The unstable modes are qualitatively the same as in Fig. 3, up to a rescaling due to domain size. This suggests that the difference in volume expansion in this case does not have a substantial effect on these dispersion sets. In order to better understand the role of volume expansion, in Fig. 11 we compare dispersion sets for spheres of different dimension undergoing linear growth. The spectrum for each is similar, yet due to differences in the volume expansion term, we find that the dispersion sets collapse for large enough $N$, since for such cases volume expansion is far more rapid, resulting in more rapid dilution of the spatially homogeneous state which prevents spatial patterning. Hence, there are indeed differences in patterning due to volume expansion, yet depending on the growth function selected, these differences may manifest only for large $N$.

6 Domain evolution with area or volume conservation

Much literature on evolving domains considers strict growth or expansion of the domain ($\dot{\mu}(t)>0$). However, there are a variety of situations for which area or volume should be preserved while the underlying domain evolves such as in the buckling of intestinal crypts during development [25], and cell shape changes [38, 107] or stationary shapes of vesicles [96]. For such a case, $\dot{\mu}(t)\equiv 0$, and (11) admits a constant exact solution $\mathbf{U}(t)\equiv\mathbf{U}^{*}$, akin to what is considered in the Turing conditions for static domains. For this case, $J$ is a constant matrix, and the condition from Theorem 3.3 reduces to

[TABLE]

This inequality is close to the static domain Turing condition, modified to account for the time dependence of the spectral parameter $\rho_{k}(t)$. We are unaware of any studies on volume-conserving domain evolution presently considered in the literature, and so give two examples.

6.1 Area conserving evolution of a rectangular domain

Regarding asymmetric growth of a rectangular domain, Turing patterning when growth in only one direction with the other direction held fixed was considered in [71]. Consider the evolution of a domain according to the coordinates $\mathbf{X}=(r_{1}(t)x_{1},r_{2}(t)x_{2})$, $x_{1},x_{2}\in[0,1]$, where $r_{1}(t)r_{2}(t)=A$, and consider $r_{1}(0)=A_{1}$ and $r_{2}(0)=A_{2}$ (with $A_{1}A_{2}=A$), along with $(r_{1}(t_{f}),r_{2}(t_{f}))=(A_{1}R^{-1},A_{2}R)$ for some constant $R>0$. The domain is then described by $\Omega(t)=[0,r_{1}(t)]\times[0,r_{2}(t)]$ with $\Omega(0)=[0,A_{1}]\times[0,A_{2}]$, $\Omega(t_{f})=[0,A_{1}R^{-1}]\times[0,A_{2}R]$, and $|\Omega(t)|=A$ for all $t\geq 0$, so this manner of growth does indeed preserve area. The spectrum is $\rho_{k_{1},k_{2}}(t)=\frac{\pi^{2}k_{1}^{2}}{r_{1}(t)^{2}}+\frac{\pi^{2}k_{2}^{2}}{r_{2}(t)^{2}}$, $k_{1},k_{2}\in\mathbb{N}$, which is an example were the spectrum is not necessarily monotone in time. Without loss of generality, we choose $r_{1}(t)=r(t)$ and $r_{2}(t)=\frac{A}{r(t)}$, where $0<r(t)<\infty$, and we then write $\rho_{k_{1},k_{2}}(t)=\frac{\pi^{2}k_{1}^{2}}{r(t)^{2}}+\frac{\pi^{2}k_{2}^{2}}{A^{2}}r(t)^{2}$ for $k_{1},k_{2}\in\mathbb{N}$.

Under the criterion given by (75), we find that the $(k_{1},k_{2})$th perturbation of the form (12) corresponding to $\rho_{k_{1},k_{2}}(t)=\frac{\pi^{2}k_{1}^{2}}{r(t)^{2}}+\frac{\pi^{2}k_{2}^{2}}{A^{2}}r(t)^{2}$ is unstable over some interval $t\in\mathcal{I}_{k_{1},k_{2}}$, provided that the inequality

[TABLE]

holds for all $t\in\mathcal{I}_{k_{1},k_{2}}$. In the special case of a static rectangle $[0,L_{1}]\times[0,L_{2}]$, the condition (LABEL:exarea) reduces to

[TABLE]

which is the classical Turing condition for a static rectangular domain.

We consider the evolution of such a domain in Fig. 12, showing both slow and fast linear evolution of the domain. For both cases, we find that the thinner rectangular domains admit solutions with more spots than when the rectangle passes through a transient square configuration. The final configuration in both cases similarly admits modes only in the $x$ direction, as one might expect from a quasi-static analysis of thin domains. The faster evolution leads to more disordered structures, as it is further from a true quasi-static picture, and the nonlinear reaction kinetics are unable to stabilize ordered spatial patterns in these transient time periods. We also plot dispersion sets $\mathcal{K}_{t}$ which characterize the unstable modes. The collection of unstable modes seen in Fig. 13 corresponds well with the patterns present at each respective time in Fig.12. While the unstable modes are quite similar in each case (with only modest differences present on transient timescales - see $\mathcal{K}_{0.1667t_{f}}$), the combination of growth with nonlinear selection of patterns does play a role in which unstable modes are selected in the final pattern.

6.2 Volume conserving evolution of a solid cylinder

Three-dimensional Turing patterns have been explored in many systems and geometries, though a complete categorization of such structures and criteria for when they emerge does not yet exist as far as we are aware [10, 22, 57, 100]. This is in contrast to the theory in two-dimensions, for which a reasonable classification of patterns exists, at least on rectangular domains [26, 99]. Such emergent structures have even been observed to be suitable for a variety of applications, such as the design of water filters [106]. Pattern formation on growing cylindrical domains is of strong relevance to plant growth [67]. Simulations and experimental observations of various Turing patterns in static cylindrical domains were shown in [3], and it was shown that three-dimensional Turing patterns exhibit an extraordinarily richer set of patterns than in one or two dimensions. Finally we remark that [38] compared linear analysis and simulations on quasi-static cylinders and spheres to argue that the flattening of cells can have an impact on mitosis, and specifically during cytokinesis when cell shape changes regularly occur.

Consider coordinates $\mathbf{X}=(r_{1}(t)x_{1}\cos(2\pi x_{2}),r_{1}(t)x_{1}\sin(2\pi x_{2}),r_{2}(t)x_{3})$, where $x_{1},x_{2},x_{3}\in[0,1]$, which defines a cylindrical domain $\Omega(t)$. We consider $r_{1}(t)^{2}r_{2}(t)=V$ so that the volume of the cylinder is conserved. Choose $r_{1}(0)=V_{1}$, $r_{2}(0)=V_{2}$ such that $V_{1}^{2}V_{2}=V$, and $(r_{1}(t_{f}),r_{2}(t_{f}))=\left(V_{1}R^{-1},V_{2}R^{2}\right)$ for some constant $R>0$. Then $\Omega(t)=D(r_{1}(t))\times[0,r_{2}(t)]$, where $D(r_{1}(t))$ denotes a disk of radius $r_{1}(t)$ centered at the origin, with $\Omega(0)=D(V_{1})\times[0,V_{2}]$ and $\Omega(t_{f})=D(V_{1}R^{-1})\times\left[0,V_{2}R^{2}\right]$. We have that $|\Omega(t)|=\pi r_{1}(t)^{2}r_{2}(t)=\pi V_{1}^{2}V_{2}=\pi V$ for all $t\geq 0$, hence volume is conserved. The spectrum of the Laplace-Beltrami operator over this domain will take the form $\rho_{\ell,k,m}(t)=\frac{\tilde{j}_{\ell,k}^{2}}{r_{1}(t)^{2}}+\frac{\pi^{2}m^{2}}{r_{2}(t)^{2}}$, for $\ell,k,m\in\mathbb{N}$. Recalling $r_{1}(t)^{2}r_{2}(t)=V$, we write $r_{1}(t)=r(t)$ and $r_{2}(t)=\frac{V}{r(t)^{2}}$ for some function $r(t)$ satisfying $0<r(t)<\infty$. The time-dependent spectrum then becomes $\rho_{\ell,k,m}(t)=\frac{\tilde{j}_{\ell,k}^{2}}{r(t)^{2}}+\frac{\pi^{2}m^{2}r(t)^{4}}{V^{2}}$.

Under the criterion given in (75), we find that the $(\ell,k,m)$th perturbation of the form (12) corresponding to $\rho_{\ell,k,m}(t)$ is unstable over some interval $t\in\mathcal{I}_{\ell,k,m}$, provided that the inequality

[TABLE]

holds for all $t\in\mathcal{I}_{\ell,k,m}$. In the special case of a static cylinder of radius $R$ and height $H$, the Turing condition (LABEL:exvolume) reduces to

[TABLE]

We show simulations of such evolving cylinders in Fig. 14, where we threshold the solutions to focus on the regions of high $u_{1}$ concentration. For this choice of parameters, the typical three–dimensional Turing structures are arrangements of spheres of high activator concentration, with partial spheres on the boundary. We see these spheres form quickly, and their number and placement change during slow growth. Eventually these structures flatten into quasi-two-dimensional cylinders near the end of the simulation shown in (a). In contrast, the structures in (b) do not have time to organize into spheres as the domain evolves more rapidly, and so we see multiple regions coalescing and mixing. The flattening and distortion of three-dimensional structures was discussed in [38], but the influence of rapid domain evolution was neglected, and it is clear that it plays a role in the emergent structures for such a domain. Dispersion sets are shown for each case in Fig. 15, and as in Fig. 13, we observe that modes of similar magnitudes are excited throughout the simulation times, but that modes corresponding to the vertical and radial directions are excited at different times. Note that rapid growth can especially influence transient mode selection (consider $\mathcal{K}_{0}$), and hence depending on the nonlinearities involved, the final patterned state.

7 Applications to systems with higher-order spatial derivatives

As a specific example of a high-order spatial system, we focus on the scalar Swift-Hohenberg equation [105]. In addition to being a canonical model for pattern formation on static domains, we remark that in the past few years this equation has been studied in evolving domains [46, 54], in particular on one-dimensional intervals of the form $[0,r(t)]$. We choose to study the following form of the Swift-Hohenberg equation:

[TABLE]

where $d>0$ is the diffusion coefficient, while $q$ and $a$ are positive parameters. We choose this form as it admits a single stable steady state, $u=0$. We again mention that dynamics from the complex Swift-Hohenberg equation, with $u^{3}$ replaced by $|u|^{2}u$ and the sign of $a$ reversed, were considered on an evolving domain by [46] and [54], with the domain being a time-dependent interval. Therefore, the Swift-Hohenberg equation is a natural canonical example of a fourth-order equation which permits pattern formation in the presence of an evolving domain. In order to match (80) with the theory of Sections 2.4 and 3.3, we note that the choice of $\mathcal{P}(y)=-y^{2}+2qy$ corresponds to (80).

As (80) is a good deal different from the reaction-diffusion equations (3) we have studied thus far, we first consider the Turing instability condition for (80) on a static domain. In the static case, $\Omega(t)=\Omega(0)$. The zero steady state is always stable to homogeneous perturbations. A perturbation of the form (12) of this steady state is governed by

[TABLE]

The perturbation corresponding to the $k$th eigenfunction is unstable provided $-d\left(\rho_{k}^{2}-2q\rho_{k}\right)-a>0$, which is in turn possible provided that $\rho_{k}$ satisfies the inequality

[TABLE]

Clearly, we must have that $q>\sqrt{\frac{a}{d}}$, in order for there to be the possibility of a Turing instability. As $\rho_{k}\rightarrow\infty$ with $k\rightarrow\infty$, there are at most finitely many such $\rho_{k}$ satisfying (82), as in the standard results on Turing instabilities in systems of two reaction-diffusion equations with second-order space derivatives.

Returning to the dynamics of (80) on an evolving domain, an application of Theorem 3.5 gives that the condition for transient instability of the $k$th perturbation (12) is

[TABLE]

where the base state $U(t)=0$ no longer depends on time. Therefore, there is a transient instability due to the $k$th spatial mode for $t\in\mathcal{I}_{k}$ provided that $\rho_{k}(t)$ satisfies

[TABLE]

for all $t\in\mathcal{I}_{k}$.

We remark that this analysis is substantially simpler than that for the reaction-diffusion systems in the preceding sections for two reasons. Firstly it is simpler due to the scalar nature of the condition from Theorem 3.5, and secondly because of the simple choice of kinetics which give a trivial base state (and hence no nonautonomous impacts from the base state’s variation). Even in cases with kinetics admitting a nontrivial base state, we remark that such a function would evolve according to a scalar ODE, and hence essentially be slaved to the dynamics of $\frac{\dot{\mu}(t)}{\mu(t)}$.

We further remark that although we have used the spatial perturbation (12) corresponding to a Laplace-Beltrami eigenfunction, and hence invoked Theorem 3.5, we shall later see the linear instability result is in good agreement with the patterns which emerge from numerical simulations. Therefore, the simplified perturbation and resulting instability criteria given in Theorem 3.5 did not discard any useful information about the onset of instability. As such, we did not have to consider the more complicated spectral problem for the fourth-order spatial operator in order to obtain a more general spectrum as used in the instability criteria of Theorem 3.6. That said, for more complicated spatial operators $\mathcal{P}$ or more complicated spatial domains $\Omega(t)$, one should be aware that the instability criteria of Theorem 3.6 may be necessary.

7.1 Swift-Hohenberg on an evolving interval

As in Section 4.3, we first consider the dynamics of the Swift-Hohenberg equation (80) on an evolving line segment, $x\in[0,r(t)]$. As before, we consider the instability of modes with eigenvalues given by $\rho_{k}(t)=\frac{\pi^{2}k^{2}}{r(t)^{2}}$, and volume expansion $\mu(t)=r(t)$. Broadly we find similar behaviours to that found in the reaction-diffusion setting, with the dispersion sets providing an approximation to the observed patterns.

We first consider exponentially growing domains in Fig. 16. As in Fig. 1, we observe mode-doubling behaviour for some growth rates, and the observed modes are broadly characterized by the dispersion sets. However, we note that the dispersion sets and our instability analysis is much less sensitive to the rate of domain evolution than in the coupled reaction-diffusion system setting, which is likely due to the trivial base state. Nevertheless, full numerical simulations depend on the rate of domain evolution due to nonlinear pattern evolution, and hence we cannot capture the variations in unstable modes precisely.

Next we consider a periodically evolving interval in Fig. 17. Qualitatively the dynamics are very similar to the Schnakenberg example on such a domain given in Fig. 4. We do remark that, again, the dispersion sets are more regular than in cases where the base state can exhibit complex dynamics. We also note that the linear analysis does not explain the loss of pattern in Fig. 4b(ii), where seemingly the rapid movement between unstable modes does not leave enough time for a net growth of any particular instability, and hence the transient pattern eventually falls back into the homogeneous equilibrium state.

7.2 Swift-Hohenberg on an evolving sphere

Swift-Hohenberg equations and generalizations thereof have been studied on static spherical surfaces [66, 101], with the model associated with the study of wrinkling and buckling of membranes [103].

In Fig. 18 we give an example of solutions to (80) on the surface of a linearly growing 2-sphere. As anticipated from the one dimensional examples, the dispersion sets for increasing growth rate are very similar, having only small differences in the size of the unstable region at the beginning of the growth interval (this is due to the form of $\frac{\dot{\mu}(t)}{\mu(t)}$ for linear growth, which is much larger for small $t$). The simulations shown indicate that even this small variation can have a marked impact on the timescale at which patterns form; in the case of $s=0.01$, patterns are developed (e.g. $\max(u)\sim O(1)$) as early as $t=7\approx 0.002t_{f}$, whereas for $s=0.15$, we see that the initial $O(10^{-1})$ perturbation has substantially decayed by $t=0.2t_{f}$, though a pattern does eventually form. As in the reaction-diffusion systems on manifolds shown before, these more rapidly grown manifolds lead to less regular patterns, though the wavelength is approximately the same.

8 Discussion

We have reviewed and generalized much of the contemporary work studying diffusion-driven instabilities in reaction-diffusion systems on growing domains. Our extensions include properly accounting for the non-autonomous nature of the base state of the system which is perturbed, allowing for general dilational evolution of the domain (of which a special case is the more commonly studied isotropic growth), and deriving a differential inequality (involving model parameters and the Laplace-Beltrami spectrum) to determine if a specific mode becomes unstable during a given time frame. Theorem 3.3 is a natural generalization of the Turing conditions on a static domain, yet explicitly accounts for the history-dependence due to the non-autonomous nature of the problem [45], and allows for arbitrary growth functions without the need to rely on slow growth or other simplifying assumptions.

To summarize, the instability conditions given in Theorem 3.3 take the form

[TABLE]

Here “reaction kinetic terms” include terms resulting from the dynamics of (11) which have been linearized (including terms involving the Jacobian matrix $J(\mathbf{U}(t))$), while “domain evolution terms” involve terms which are specific to the manner of domain growth, such as terms involving the time derivatives of $\mu(t)$ and $\rho_{k}(t)$. If the instability is diffusion driven, then for the $k=0$ mode $\rho_{k}(t)\equiv 0$, we should have

[TABLE]

If (86) does not hold, then there is some homogeneous instability not due to diffusion. If (86) holds, while (85) holds for a particular index $k=k^{*}>0$ and a particular interval $\mathcal{I}_{k^{*}}$, then the spatial perturbation (12) corresponding to this particular $k^{*}$ inducing an instability for $t\in\mathcal{I}_{k^{*}}$. When the domain is static, the domain evolution terms vanish, and we are left with instability conditions of the form

[TABLE]

while the $k=0$ mode (86) is stable when reaction kinetic terms are positive, which is just the classical Turing condition on a static domain. Therefore, the condition given in Theorem 3.3 is a fairly natural generalization of the classical Turing conditions. Hence, although we have made few assumptions, and have avoided both asymptotic approximations of growth functions or assuming a constant base state, within our general framework we have captured the spirit of the original Turing conditions. In the limit of no growth, our results recover the Turing conditions for a static domain in a completely natural manner, without further effort or appeal to simplifications. Similarly, instability conditions for particular kinds of slow growth known in the literature (such as slow exponential growth), as well as quasi-static approximations for the slow growth regime, fall out of our results in the relevant limits.

Due to the time-dependent nature of growth terms and of the Laplace-Beltrami spectra, we have phrased our results in terms of instabilities present over a given time interval, rather than as $t\rightarrow\infty$ like in the classical Turing instability. This more general approach allows for the understanding of transient instabilities. This is a useful generalization, as in practice Turing patterns are selected in finite time, with patterns then remaining static once formed. As such, there is an interplay between the rate of growth of the domain and the rate at which this mode selection occurs. Hence, pattern formation will rely strongly on when certain modes result in instability on the timescale of the reaction kinetics, rather than simply if such modes ever induce instability at any arbitrary time. In contrast to the standard Turing theory for static domains, we conjecture that it may be the time duration for which a mode remains unstable that matters more than the degree or magnitude of instability at any instantaneous time (which is often considered by comparing the real part of eigenvalues in the static case in the limit $t\rightarrow\infty$). Of course, one would need nonlinear theory to address such issues (which is beyond the scope of the present paper), and even then, resolution is likely only on a case-by-case basis for given reaction kinetics, growth functions, and spatial patterns.

Our numerical examples illustrate both the power of our analytical results (the dispersion sets obtained broadly match simulations, with some exception), and the large variety of phenomena one can expect when studying such problems. We remark that we intentionally only considered a small fraction of the parameter spaces for only two kinds of reaction kinetics, and already have observed behaviours which are qualitatively distinct from what has been commonly observed in related Turing systems in the past. In particular, we note that our main choice of the Schnackenberg kinetics (65) gives a relatively simple Jacobian, and a base state with an uncomplicated evolution equation. More complicated reaction kinetics can likely lead to many new phenomena due to the complexity of non-autonomous phase spaces. Our results will immediately apply in such cases, providing insight into spatial instabilities around even a time-dependent base state, as shown for instance in Sec. 4.4.

Similarly, we have only considered a handful of growth functions and domain geometries, but unlike many results in the literature, the instability condition in Theorem 3.3 applies so long as one can compute derivatives of the growth functions and the time-dependent Laplace-Beltrami spectrum of the domain. We have used the cases of volume-preserving evolution and the Swift-Hohenberg example to show that even when the base state is static, domain evolution can substantially change the unstable modes observed, along with the qualitative behavior of spatial patterns altogether (cf. Sec. 6). Hence, changes in the structure of a domain are sufficient to modify the linear stability properties, as well as the patterns formed, even if the area or volume of the domain is fixed. Indeed, these examples show some of the largest discrepancies between quasi-static or asymptotically large-time approaches common in the literature (which tend to work best when evolution is slow and monotone) and pattern selection which actually occurs at transient timescales, as seen in the strong time dependence of not only the extent but also the shape of the dispersion sets $\mathcal{K}_{t}$ shown in Figs. 12-14.

We have also demonstrated the application of the same basic instability analysis on a canonical higher-order scalar spatial system, the Swift-Hohenberg equation. This demonstrates the generality of our approach, which is to be expected of a linear stability analysis. Additionally, this model is a much simpler test bed for exploring a variety of nonlinear kinetics and complicated evolution scenarios, as one can explore the influence of growth in the absence of a base state which can exhibit complex temporal dynamics. Extensions to planar domains could also be considered as generalizations of planar studies in the literature [60, 36]. These would be especially simple in the case of area-preserving evolution, as the nonautonomous dynamics would then be entirely embedded in the time-dependent spectral parameters.

The framework developed in this paper may also be applied to problems with time-dependent reaction kinetics, on either growing or static domains. In particular, temporal oscillations have been employed in photosensitive reactions to control Turing patterns, and in some cases eliminate them [24, 37, 119]. Spatiotemporal forcing has also been used to mimic domain growth in such systems [49, 71, 72, 90]. As our approach allows for non-autonomous Jacobian matrices, such applications will naturally benefit from our approach. Another natural extension is to consider the problem of domain evolution simultaneous to an imposed flow, to study problems for which the velocity due to dilution is augmented with a velocity field which transports the chemical species or morphogens in some manner. Recently, the flow properties and cross-section geometry for activator-inhibitor systems within a tube were shown to influence emergent Turing patterns [115]. The extension of such results to scenarios where the tube dilates periodically in time would be one such extension incorporating evolving spatial domains which would have biochemical and physiological relevance. Additionally, spatial heterogeneity in reaction-diffusion systems is also of contemporary interest, and incorporating this alongside growth would lead to a much more biologically realistic extension of Turing’s theory [80, 53]. We also mention that many other contemporary problems in a range of fields consist of systems on time-evolving domains [46, 47]. While a general theory of such systems does not exist, understanding the impact of growth for reaction-diffusion systems can provide an important example for these larger classes of systems. We anticipate that the results presented here can be readily generalized to other settings.