Influence of Curvature, Growth, and Anisotropy on the Evolution of Turing Patterns on Growing Manifolds

We study two-species reaction–diffusion systems on growing manifolds, including situations where the growth is anisotropic yet dilational in nature. In contrast to the literature on linear instabilities in such systems, we study how growth and anisotropy impact the qualitative properties of nonlinear patterned states which have formed before growth is initiated. We produce numerical solutions to numerous reaction–diffusion systems with varying reaction kinetics, manner of growth (both isotropic and anisotropic), and timescales of growth on both planar elliptical and curved ellipsoidal domains. We find that in some parameter regimes, some of these factors have a negligible effect on the long-time patterned state. On the other hand, we find that some of these factors play a role in determining the patterns formed on surfaces and that anisotropic growth can produce qualitatively different patterns to those formed under isotropic growth.


Introduction
Since Turing first proposed reaction-diffusion systems as a model for pattern formation (Turing 1952), much work has been performed to understand the theoretical and biological aspects of this phenomenon (Satnoianu et al. 2000;Green and Sharpe 2015;Marcon et al. 2016;Woolley et al. 2017). Murray (2003) discusses 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. Curved and complex geometries have been of recent interest, as these provide for more realistic regimes B Robert A. Van Gorder within which to test morphogenetic hypotheses (Tse et al. 2010;Trinh and Ward 2016;Núñez-López et al. 2017). Since reaction-diffusion systems are more difficult to analyze on growing domains, pattern formation has been considered on different-sized static domains to simulate very slow growth (Varea et al. 1997). This requires the reaction and diffusion of the chemical species to occur on a much faster timescale than the growth and also be independent of the growth. An alternative simplification would be to assume growth occurs on a much faster timescale than reaction and diffusion and consider the quasi-static approximation of the reaction-diffusion system on the final shape of the domain. Since neither of these approximations are likely to be valid for all systems modeling biological growth, we consider reaction-diffusion systems where growth, reaction, and diffusion occur on comparable timescales. Crampin et al. (1999) explicitly considered uniform and isotropic domain growth in one-dimensional reaction-diffusion systems in the slow and fast growth regimes, demonstrating frequency doubling of the emergent Turing patterns. This approach was discussed in the context of biological patterning problems in Crampin and Maini (2001), Barrass et al. (2006) and extended in Crampin et al. (2002) to consider the influence of non-uniform domain growth on one-dimensional reaction-diffusion systems, including apical or boundary growth. Castillo et al. (2016) considered Turing and Turing-Hopf instabilities in an exponentially and isotropically growing square and suggested that anisotropy and curvature are important considerations for extending their analysis. Following work on characterizing instabilities in nonautonomous reaction-diffusion equations on slowly growing domains in Madzvamuse et al. (2010), Hetzer et al. (2012) derived conditions for diffusion-driven instability on time-dependent manifolds. Plaza et al. (2004) derived a general formulation of reaction-diffusion theory on isotropically evolving one and two-dimensional manifolds, with motivation from biological settings where growth and curvature both play a role in organism development.
Beyond computing linear instability criteria, some authors have analytically explored how patterns change and evolve under growth (Comanici and Golubitsky 2008) by exploiting the framework of amplitude (or Ginzburg-Landau) equations (Cross and Hohenberg 1993;van Saarloos et al. 1994). 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. Additionally, to our knowledge, there does not exist an analytical framework applicable to curved manifolds. In contrast, many more studies have explored such systems numerically, developing a variety of techniques to capture the transient and the final patterned states in reaction-diffusion systems on manifolds (Madzvamuse et al. 2003; Barreira et al. 2011;Macdonald et al. 2013;Madzvamuse and Maini 2007;Tuncer and Madzvamuse 2017). While these numerical methods and their applications will always be isolated to specific instances of reaction-diffusion systems, one can hope to at least explore the kinds of emergent behaviors in candidate systems as a way to understand realistic pattern formation and pattern evolution in more complex settings.
All 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 (Ubeda-Tomás et al. 2008;Corson et al. 2009;Peaucelle et al. 2015). 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 (Menzel 2005;Saez et al. 2007;Bittig et al. 2008;Amar and Jia 2013). Madzvamuse and Barreira (2014) considered anisotropic and concentration-dependent growth in the context of cross-diffusion models, but their emphasis was on cross-diffusion-driven instabilities in complicated settings, rather than on understanding anisotropic growth in simple geometries, which has been suggested as a useful direction (Castillo et al. 2016). Rossi et al. (2016) also studied concentration-dependent growth of a scalar reactiondiffusion equation on a time-dependent manifold. Krause et al. (2018a) investigated the effects of Turing instabilities in two-species reaction-advection-diffusion models on a sphere for several well-known reaction kinetics. They find that advection and the compact geometry of the sphere allow for complex spatiotemporal patterns. Since growth can induce advection-like terms into the morphogen evolution equations (Plaza et al. 2004), it is of interest to consider the effect of anisotropic growth on compact domains, such as the sphere.
Much of the work in reaction-diffusion theory on growing domains has emphasized differences between quasi-static, slow, and fast growth regimes. Growth can lead to a robustness of patterning processes in some scenarios (Crampin et al. 1999), where the period-doubling of pattern splitting from an existing pattern leads to insensitivity to the initial symmetry-breaking perturbation which initiated patterning. Despite such differences, most reaction-diffusion systems give rise to qualitatively similar patterns, with quantities such as the pattern wavelength determined primarily by the final domain size, independent of growth. Determining when different kinds of growth influence the final pattern, and the effects that growth has on patterns, is an important problem to determine the validity of quasi-static approximations. Indeed, other questions arise such as: To what extent does the rate of deformation, type of growth or extremity of growth affect the patterns formed on the domain, if there is any effect at all?
In this paper, we generalize the modeling of reaction-diffusion systems on growing manifolds developed by Plaza et al. (2004) to allow for dilational anisotropic growth. The restriction of dilational growth allows for relatively simple governing equations amenable to simulation, yet still allows us to compare between quasi-static, isotropic, and anisotropic growth regimes. We then systematically explore three kinds of reaction kinetics on planar and curved two-dimensional manifolds. We also vary the kinds of growth, both in terms of average growth rate, functional form of the growth, and exploring growth in different stages. To address the questions of when these factors matter, we restrict attention to the evolution of established patterns far from homogeneous equilibrium states that arise due to Turing instabilities. In other words, all parameters are chosen within the Turing space such that these systems will always develop a heterogeneous solution, and we focus on qualitative features which change in such a solution as the manifold grows. This allows us to determine when each of these different growth scenarios gives rise to qualitatively different patterns. We emphasize qualitative differences and very pronounced quantitative effects, but do not focus on small quantitative differences between simulations, as these can depend sensitively on parameters and the initial data. We find that all of these features (growth rates, anisotropy, curvature) can give rise to qualitatively different patterns in some regimes, but that some of them (e.g., growth rate) are typically more important than others (e.g., the functional form of growth), and that some reaction kinetics, and some parameter regimes, are more or less influenced by growth.
The remainder of this paper is organized as follows. In Sect. 2, we develop our models to study two-species reaction-diffusion systems on a general two-dimensional manifold undergoing isotropic or anisotropic, but dilational growth. In Sect. 3, we discuss a systematic framework by which we study the effects of such growth on a range of reaction kinetics, demonstrating systems where patterning is generally robust (e.g., qualitatively independent of growth). In the following sections we shall employ this systematic approach to study Turing patterns on growing domains in the case where the domain is a flat ellipse (Sect. 4) and in the case where the domain is an ellipsoid surface (Sect. 5), the former demonstrating the role of anisotropy and the latter additionally demonstrating the role played by local curvature of the domain. We then discuss the evolution and numerical stability of stripe and target patterns in the presence of growth in Sect. 6, considering both the ellipse and ellipsoid domains. We finally discuss the implications of our results in Sect. 7.

Mathematical Model for Anisotropic Dilational Growth
We now introduce the model accounting for reaction-diffusion on a growing manifold. We first define the kinds of stationary patterns we are interested in and recall the conditions for diffusion-driven instability. We then introduce the concept of a reactiondiffusion system on a growing manifold and show how to transform this system onto a stationary manifold. This is a generalization of Plaza et al. (2004) where we allow for anisotropic dilational growth. Finally, we apply this general framework to two example manifolds which we will simulate in the following sections: a planar ellipse and the surface of an ellipsoid.

Turing Instability
We consider a general two-species reaction-diffusion system on a stationary manifold Ω given by where u and v are the two interacting chemical species, δ 1 and δ 2 are the diffusion coefficients of u and v, respectively, ∇ 2 is the Laplace-Beltrami operator on Ω, and f and g are nonlinear kinetic functions (reaction terms) to be specified. We will consider manifolds with and without boundaries. When boundaries are present, we augment Eq.
(1) with the conditions where n is the outward unit normal. In this paper, we consider patterns defined as stable, time-independent, spatially heterogeneous solutions to Eq. (1). In general, nonlinear reaction-diffusion systems can give rise to a variety of spatiotemporal behaviors, but we restrict our analysis to stationary patterns, which were first predicted to occur by Turing (1952) in the presence of differential diffusion (e.g., δ 1 = δ 2 ). In the case of stationary manifolds, conditions can be derived for the instability of spatially homogeneous solutions to Eq. (1) to small perturbations, and these will often lead to patterning. For instance, necessary conditions are given in Murray (2003) as where f u , f v , g u , and g v are the partial derivatives of f and g with respect to u and v. Such conditions have been extended to time-dependent domains (Madzvamuse et al. 2010;Klika and Gaffney 2017), as well as to isotropically growing manifolds (Plaza et al. 2004). In these cases, the conditions are more complicated, but the basic idea is the same: differences in the diffusion rates of each species allow for the growth of spatial perturbations to a homogeneous steady state and hence allow for patterning.
Here we are primarily concerned with the interactions between growth and the kinds of patterns formed and so will always study cases where diffusion-driven instability occurs in order to elucidate important differences in the kinds of growth. We note that while growth may influence the stability criteria, we intend to always let patterns emerge on a static manifold before initiating growth. In this way, we will use the above conditions to choose parameters within the Turing space such that patterns form before considering the influence of growth on such patterns.

Reaction-Diffusion Equations on Growing Two-Dimensional Manifolds
We now consider a growing domain in the form of a smooth two-dimensional manifold Ω(t) which is simply connected and compact. The manifold Ω(t) will either have a smooth boundary ∂Ω(t) for all time t ≥ 0 or will be a manifold without boundary for all time t ≥ 0. LetΩ(t) be an area element of the manifold, such thatΩ(t) ⊂ Ω(t). Let (u, v) : Ω(t) → R be a concentration function defined on the manifold Ω(t), where, for example, u, v may describe the concentration of a chemical species, or morphogen, on the manifold Ω(t). We shall assume that u and v are C 1 ([0, ∞)) in time and C 2 (Ω(t)) in spatial coordinates. The conservation of mass equation is where j denotes the flux of concentration u and dΩ is the local area element on the manifold. Using Reynold's transport theorem on the left-hand side of Eq. (4), we have where Q is the local velocity vector field generated by changes in the manifold Ω(t) and dΩ is the local area element on the manifold. An analogous equation holds for v. Note that we denote ∇ Ω(t) · as the divergence operator on Ω(t). Similarly, we define ∇ 2 Ω(t) to be the Laplace-Beltrami operator on Ω(t). By applying Eqs. (4)-(5) and using Fick's law of diffusion (Plaza et al. 2004), we have the reaction-diffusion-advection system for concentrations u and v given by: Noting that the term ∇ Ω(t) · (Qu) can be written as Q · ∇ Ω(t) u + u · ∇ Ω(t) Q, the term Q · ∇ Ω(t) u corresponds to advection due to local growth of the manifold, whereas the u · ∇ Ω(t) Q term corresponds to dilution of the concentration u due to local area changes.

Coordinate Transformation to Stationary Manifolds
In order to solve Eqs. (6a) and (6b), it is often necessary to map the equations to a stationary manifold and solve for u and v numerically in this stationary domain, before transforming the solved problem back to the original time-dependent manifold. Here, we outline a method to map the evolution equations for u and v to the stationary manifold Φ ⊂ R 2 . We follow a similar procedure to Section 2 of Plaza et al. (2004), however, in a more general form, since we do not insist upon the parametrization for Ω(t) to be orthogonal and therefore permit anisotropic growth. First, parameterize the two-dimensional manifold Ω(t) as a surface in R 3 by where α ∈ [α 1 , α 2 ] ⊂ R, β ∈ [β 1 , β 2 ] ⊂ R, and t ≥ 0. As Ω(t) is assumed compact, we have that |α 1 |, |α 2 |, |β 1 |, |β 2 | < ∞; hence, (α, β) is always confined to a rectangle in R 2 . We shall define this set to be the stationary manifold Φ; in other words, Φ = [α 1 , α 2 ] × [β 1 , β 2 ]. One can always normalize variables so that |Φ| = 1, but this is not necessary. We shall assume that except on a set of measure zero, meaning that Ω(t) is a normal surface.
To properly define the Laplace-Beltrami operator, we first define the metric tensor, G, by where The determinant of the metric tensor is while the inverse of the metric tensor, G −1 , is defined by The representation of the Laplace-Beltrami operator ∇ 2 Ω(t) on the stationary twodimensional manifold Φ is then given by In the case of a diagonal metric tensor (as seen in isotropic growth), note that G −1 12 = G −1 21 = 0, and hence, expression Eq. (13) simplifies to a form like that given in Eq. (14) of Plaza et al. (2004).
In order to remove the advection term induced by the growing manifold [the ∇ Ω(t) · (Qu) terms in Eq. (6)], we shall apply a change of variable akin to that used in Crampin et al. (1999) and Plaza et al. (2004). We shall assume that space and time variables in X are separable. Noting that the change of variable will contribute a term −Q · (∇ Ω(t) u), canceling the latter term in Eq. (14). Then, after the coordinate change, we will have a contributing advection term of the form (∇ Ω(t) · Q)u. The form of the divergence of the advection vector generated due to the change in Ω(t) is intrinsic to Ω(t) and may be calculated by the coordinate change X using again the metric tensor. We find In the case where the metric tensor is diagonal (for instance, as occurs for isotropic growth), the formula above reduces to which is equivalent to that given in Eq. (13) of Plaza et al. (2004).

Reaction-Diffusion System on the Stationary Manifold
Using the form of the Laplace-Beltrami operator given by Eq. (13) and the advection term given by Eq. (15), then we can write reaction-diffusion system Eq. (6) on the growing manifold Ω(t) as a spatially and temporally heterogeneous reaction-diffusion system on the stationary manifold Φ, obtaining

Examples of Flat and Curved Domains
We now define a pair of two-dimensional manifolds which we will use as examples for isotropic and anisotropic growth: a planar ellipse and the surface of an ellipsoid. There are many biologically relevant manifolds one could study, but these examples are sufficient to compare the influence of isotropic and anisotropic growth, and compare effects due to (boundary) curvature of the manifold.

Planar Ellipse
We parameterize a growing ellipse in the plane via the semimajor and semiminor axes, given by A(t) and B(t), respectively. These must satisfy We define the surface parametrically via where α ∈ [0, 1] and β ∈ [0, 2π). We calculate The entries of the metric tensor G are given by The determinant of the metric tensor then reads Substituting the above into the formulae for Laplace-Beltrami operator Eq. (13) and for advective growth term Eq. (15), we find after simplifying and where we denote derivatives of We note that in the case where A(t) ≡ B(t) (an isotropically growing circular disk), the system above reduces to This is the same as given for isotropic growth of a disk in Eq. (34) of Plaza et al. (2004).

Growing Tri-Axial Ellipsoid Surface
Another manifold we consider is the surface of a growing tri-axial ellipsoid. Murray (2003) discusses reaction-diffusion models as a mechanism for spiral patterning on amphibian egg shells. Since the shape of these eggs may be approximated by an ellipsoid, studying reaction-diffusion systems on these sorts of domains could help predict the patterns formed due to the calcium waves diffusing over the surface of the shell.
In order to consider reaction-diffusion system Eq. (6) on the surface of the growing tri-axial ellipsoid, we follow the same analysis as outlined in Sect. 2.3. For the growing tri-axial ellipsoid, we define the surface parametrically via where α ∈ [0, π] and β ∈ [0, 2π) are the polar and azimuthal angles, respectively. Then, we calculate Calculating the entries of the metric tensor G, we find these to be Thus, we find the determinant of the metric tensor G to be det G = g 11 g 22 − g 2 12 = A 2 B 2 cos 2 (α) + C 2 A 2 sin 2 (β) + B 2 cos 2 (β) sin 2 (α) sin 2 (α). (35) Substituting these expressions into Laplace-Beltrami operator Eq. (13) we find that where the functions M 1 and M 2 are defined by and Further, substituting the information from the metric tensor G into advective growth term Eq. (15), we have where we have defined for ease of notation. We denote the derivatives of A(t), B(t), and C(t) with respect to time asȦ = d A dt ,Ḃ = dB dt , andĊ = dC dt , respectively. Then, reaction-diffusion system Eq. (6) can be put into the form In the case where where the domain is a growing spherical surface under isotropic growth), the reaction-diffusion system above reduces to In the special case of a spheroid,

Methods for the Study of Growth Rates, Curvature, and Anisotropic Growth
In this section, we will describe the various kinds of reaction kinetics that we have explored, as well as our systematic approach to studying them on growing domains. We are chiefly interested in fully nonlinear and stationary patterns, rather than transient patterning, so we always allow the reaction-diffusion process to settle onto a steady state on a stationary manifold before and after it has undergone growth. In this way, we can compare the kinds of stationary patterns which are arrived at via different growth processes. We will proceed to discuss parameter choices and numerical methods which will be used in our analysis.

Types of Reaction Kinetics
There are a multitude of functions f , g which can represent the interaction between two morphogens. We consider three different reaction schemes: Schnakenberg, Gierer-Meinhardt, and FitzHugh-Nagumo. These are canonical examples of crossand pure activator-inhibitor kinetics, as well as excitable media, respectively.

Schnakenberg Reaction Kinetics
First, we consider Schnakenberg (also known as activator-depleted) kinetics, a wellstudied form of reaction which models activator-inhibitor behavior (Schnakenberg 1979). They were originally developed to explore simple chemical kinetics that gave rise to oscillatory behavior, but they have since been used in a huge amount of literature as a canonical example of activator-inhibitor chemistry (Iron et al. 2004;Liu et al. 2013;Sarfaraz and Madzvamuse 2017).
Taking the non-dimensionalized form of the Schnakenberg kinetic system as in Krause et al. (2018a), the reaction terms f , g in Eq.
(1) are given by where a and b are parameters to be chosen. In order for patterns to form, the parameters a, b, δ 1 , δ 2 must satisfy Turing instability condition Eq.
(3) given in Sect. 2.1. Following the methodology outlined in Sect. 2.1, we first find the homogeneous steady state of the system by setting f = 0 = g. This gives the steady-state u * and the Jacobian of the linearized system J * as Turing instability condition Eqs.
(3) can be written as:

Gierer-Meinhardt Reaction Kinetics
We also consider Gierer-Meinhardt reaction kinetics (Gierer and Meinhardt 1972). These have been used in a variety of theoretical and applied contexts and are particularly used as an example of Turing instabilities which lead to spatially localized patterns, such as spikes in 1-D domains and spots in 2-D domains (Wei and Winter 2013). We note in particular the study by Tse et al. (2010) of these kinetics on the sphere, demonstrating the effect of curvature on the stability and structure of patterns.
In non-dimensional form, these equations are: where the terms au and bv correspond to linear degradation of the morphogens, at rates a and b, respectively, and the terms involving u 2 correspond to the activator's production of both u and v. We note the addition of the constant 1 to the activator kinetics, corresponding to basal production of the activator, which allows for robust spontaneous spot formation (Page et al. 2005;Krause et al. 2018b).
As for the Schnakenberg reaction kinetics, we find the homogeneous steady state of Eq. (47) by setting f = 0 = g. This gives the steady-state u * and the Jacobian of the linearized system J * : The conditions for Turing instability (3) are:

FitzHugh-Nagumo Reaction Kinetics
The final type of reaction scheme we consider are the FitzHugh-Nagumo kinetics, derived independently by FitzHugh (1955FitzHugh ( , 1961 and Nagumo et al. (1962) which is a two-component simplification of the famous Hodgkin-Huxley model (Hodgkin and Huxley 1952), which itself was proposed following experiments with the giant squid axon. The Hodgkin-Huxley equations model the potential across a surface membrane as the movement of potassium and sodium ions (plus a small leakage current of other ions). The FitzHugh-Nagumo reaction kinetics can lead to Turing-Hopf bifurcations (Castillo et al. 2016), as well as a multitude of behavior due to bistability (Keener and Sneyd 1998). These kinetics have also been used to model labyrinthine patterning in neurogenesis (Cartwright 2002). For simplicity, we only consider parameters where there is a unique positive steady state and we use parameters which are not in the Turing-Hopf bifurcation regime, in order to focus on the behavior of stationary patterns under growth. The homogeneous steady state in this case will only be unstable to modes arising from a pure Turing bifurcation. The FitzHugh-Nagumo kinetics can be written as where a, b, c, d > 0 are constant parameters chosen to satisfy the Turing instability conditions given in Sect. 2.1, and in this case we take δ 1 = 1 and d = 0.6. Following the analysis of Castillo et al. (2016), we require to ensure that there is a unique positive steady state for the homogeneous problem. That is, there is only one positive solution to the system f noting that v * = a−u * b and thus that we require 0 < u * < a for v * > 0, to obtain a positive steady state. Thus, the steady-state u * and Jacobian J * are The Turing instability conditions, for the instability of the steady-state u * are,

Numerical Approach
Throughout the rest of the paper, we will demonstrate solutions to equations (25) or (41) for different reaction kinetics f and g, as well as in different growth scenarios. We will use the commercially available finite element solver COMSOL, version 5.3, which will discretize the manifolds using second-order triangular finite elements. We note that these simulations were checked in various static domain cases using the Matlab package Chebfun (Townsend and Trefethen 2013;Townsend et al. 2016), 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 (but let the solver increase the time step freely as the solution evolved). Convergence in time was checked by restricting the maximum time step, and convergence in space was determined via using different numbers of finite elements and comparing the norm of solutions over time and space. For the growing ellipse, we used 25,970 domain elements, and for the ellipsoid we used 19,704 boundary elements, which we now describe in some detail. 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 R n+1 Elliott 2007, 2013;Olshanskii and Xu 2017), so that dilation of a particular coordinate in R n+1 allows a natural construction of the Laplace-Beltrami operator on the surface. We note that the planar ellipse model given by (25) is a flat 2-D manifold, and so the projection in this case is trivial. For Eq. (41), however, we can define both the Laplace-Beltrami operator and the dilution term in the ambient space R 3 as the operators with the γ 's defined by (40), and the variables (x, y, z) from the extrinsic coordinate system restricted to the surface. We can then define (41) in the whole space R 3 using (55) and simply project the equations to the unit sphere in order to obtain the correct dynamics on a growing ellipsoid. We also note that this ambient formulation allows us to see an explicit interaction between curvature and anisotropic growth-namely, the dilution term (55b) contains explicit spatial heterogeneity due to curvature and anisotropy, as these do not appear in the planar ellipse, nor in the isotropically growing sphere described in Plaza et al. (2004).

Robustness of Patterning to Growth
We now outline our systematic study of these three kinds of reaction kinetics on the elliptical and ellipsoidal domains, with dynamics governed by Eqs. (25) and (41), respectively. We vary the functional forms of A(t), B(t), and in the ellipsoidal case, C(t) in order to elucidate effects due to growth rates, anisotropy, and curvature. We consider different forms and rates of growth, in addition to the case of no growth. We note that on many domains, multistability of several non-uniform patterns is generic, especially in manifolds with dimension greater than one (Hunding 1980;Jensen et al. 1993;Borckmans et al. 1995). Because of this, we will emphasize qualitative differences between the kinds of patterns, rather than quantitative differences in the final patterned state. We now describe all of the different growth scenarios, which correspond to different choices in the parameters, manifolds, growth functions, etc. We will always take parameters within the Turing regime, so that the homogeneous steady state is unstable and leads to the formation of stationary patterns. Unless otherwise mentioned, we will always assume initial data of the form u(0) = u * + ξ , where u * is the spatially homogeneous steady state and ξ is a normally distributed spatial random variable with zero mean and standard deviation 10 −3 . We note that the same realization of the spatial noise was used in all simulations. For each scenario, we will consider a fixed static domain Ω * , such as a sphere or an ellipsoid, and compare it to a growing domain Ω(t) such that after a time period T , Ω(T ) = Ω * . We are interested in the equilibrium states on the final manifold, and we are also interested in the evolution of non-uniform states. For these reasons, for each scenario we identify a kinetic timescale T * such that the reaction-diffusion equations have approximately reached their steady-state behavior. We determine T * by simulating the equations on the static manifold Ω * and observe when the patterning process settles into a stationary state, and we check that no further pattern movement or creation occurs in the time period up to 10T * . With such a kinetic timescale in hand, we always let the pattern develop fully before initiating growth, and we always let the pattern equilibrate after growth. Mathematically this means that we take Ω(t) = Ω(0) for t ≤ T * and Ω(t) = Ω(T ) for t ≥ T − T * , so that growth is confined to the interval of time t ∈ (T * , T − T * ). We will always begin in a symmetric manifold, so that Ω(0) is either a circular disk or a sphere, and we will take the initial radius r 0 to be large enough so that patterning occurs.
We use the following functions to characterize the type of growth used and state these in terms of growth of the axis A(t) and assume that this axis has an initial value of A(0) = r 0 and a final value of A(T ) = r T . Note that these functional forms are applied analogously to all growing axes. We also note that these functions are constant at the beginning and end of the simulation time, as described above. We define linear growth, where the domain is growing at a constant rate, by: Linear growth is unlikely to occur in biological circumstances (Plaza et al. 2004); however, it is a simple form of growth frequently used in the literature and so worth investigating. We consider sublinear growth, where the domain is growing at an increasing rate (i.e.,Ä(t) > 0), given by: This growth function is exponential, which is a reasonable model for certain biological growths-e.g., for the initial growth phases of some tissues during development. Superlinear growth, where the domain is growing at a decreasing rate (i.e.,Ä(t) < 0), is defined by: This is a logistic-like growth, which is a realistic approximation of growth for some biological applications, where growth is reliant on limited external resources which can become scarce as they are depleted. We also use a punctuated form of growth defined by alternating piecewise constant and linear regions during the growing phase, but do not explicitly give it for brevity (as we will see, the particular growth function plays almost no role in the qualitative structure of the steady states). To ensure the numerical solver did not encounter any problems due to lack of continuity at t = T f , the growth function is smoothed to be continuous with a transition region of 0.01 at t = T * and t = T − T * . We consider two different growth scenarios: staged growth with a symmetric final manifold, or pure dilational growth with an elliptical or ellipsoidal final manifold. In the former, A(t) and B(t) (and in the ellipsoidal case C(t)) grow in two (three) different stages subsequently and reach the same final value of r T , so that the final domain is a larger disk (sphere). We note in this case that we also considered an analogous isotropic growth scenario for the same manifold (e.g., without stages). In the second scenario, we consider a final manifold which is not symmetric, and so only one of the axes grows leading to an elongated ellipse or ellipsoid at the final time t = T . Additionally, we consider the effect of the timescale over which the domain grows on the patterns formed, denoted as T g = T − 2T * . We consider three different timescales: fast growth given by T g = T * /2, growth comparable to the kinetic timescale so T g = T * , and growth which is much slower given by T g = 10T * . We note that T * is somewhat arbitrary, and in our simulations we take it to be sufficiently large to ensure we are in a patterned state, but we believe these growth timescales at least highlight differences in how the rate of growth influences the evolution of patterns. We note that in staged In all cases, we took δ 1 = 1 and for FitzHugh-Nagumo kinetics took c = 1 and d = 0.6. For Schnakenberg and Gierer-Meinhardt on both domains, and FitzHugh-Nagumo on the ellipsoidal domain, the initial manifold was always taken with r 0 = 10 and the growth was to a maximum of r T = 4r 0 , whereas for FitzHugh-Nagumo on the planar ellipse we used r 0 = 25 and the growth to a maximum of r T = 4r 0 growth, we increase the growth timescale so that growth in each direction occurs over the same period of time T g . We ran simulations of Eqs. (25) and (41) for each of the three reaction kinetics described in Sect. 3.1. In general, each of these models can give rise to a wide variety of behaviors, but for simplicity we consider the 6 parameter sets shown in Table 1. In each case, we explored all combinations of the 4 growth functions, 3 timescales, and both the final manifold scenarios (including an isotropic growth case), as well as 12 simulations on the static domain given by Ω * = Ω(T ), leading to a total of 444 simulations. In each case, we compared the final pattern (in terms of the activator, u, as the inhibitor can be related to it via a phase), and we summarize differences below. We emphasize that many simulations will lead to quantitatively different steady states (e.g., the number of spots may differ slightly), but these quantitative differences will also occur due to different realizations of the initial data, as the precise organization of the final non-uniform state can sensitively depend on the initial data (Maini et al. 2012). Additionally, we ran other simulations with variations in these parameters or growth scenarios to ensure that our results were robust. In particular, we did not systematically vary staged growth at different rates, but did explore some simulations in this case to understand the possible behaviors.
Broadly speaking, we found that for Schnakenberg and Gierer-Meinhardt kinetics, the kinds of patterns observed and their qualitative properties did not change much, except due to the rate of growth, or whether or not the final domain was curved. The functional form of growth is particularly unimportant, only really changing the 'effective' rate of growth in some boundary cases which we will describe below. This leads us to conjecture that for the long-time steady states observed, linear growth with varying rates in growing domains is sufficient to capture all of the qualitative dynamics that can be observed, and that static domain simulations capture many qualitative properties of the final patterned state for these kinds of reaction kinetics. FitzHugh-Nagumo kinetics result in a richer structure, in terms of its interaction with growth, curvature, and anisotropy in the domain or growth. We do note, however, that analogous comments about the functional form of growth can be made. As the functional form of growth only influences the rate of growth, we will only show figures using the Fig. 1 Simulations of Eq. (25) with kinetics given by (44) using parameters S1 for growth only along the x-axis. We plot the final distribution of u on a the static domain and b a domain with linear growth over the timescale T g = 10T * = 50,000 (Color figure online) case of linear growth, but remark that the results with all other growth laws were qualitatively identical throughout. Similarly, we only show cases where noticeable effects occurred, and omit figures with comparable behavior to those shown. Finally, we note that solutions using the Schnakenberg and Gierer-Meinhardt kinetics will always be positive, as these models correspond to morphogen concentrations, whereas the activator u in FitzHugh-Nagumo is the voltage with respect to an arbitrary gauge and hence can locally take positive and negative values.
In the following sections, we shall employ this systematic approach to study Turing patterns on growing domains in the case where the domain is a flat ellipse in Sect. 4 and in the case where the domain is an ellipsoid surface in Sect. 5. We then turn our attention to the study of striped patterns on both surfaces, in Sect. 6.

Reaction-Diffusion Systems on Planar Elliptical Domains
We now describe some of the results obtained in the case of Schnakenberg and Gierer-Meinhardt kinetics in the case of a planar elliptical domain. In Fig. 1, we compare a static ellipse with an ellipse having grown from an initial circle with growth only along the semimajor axis of the final ellipse. We note that growth appears to generate a more symmetric final distribution of spots (and that this result is identical across growth rates and functions in this case). This is consistent with the robustness attributed to growth in 1-D domains (Crampin et al. 1999). We also consider isotropic and staged growth for a circle in Fig. 2, here showing that isotropic growth, staged anisotropic growth, and the static domain all admit a distribution of spots with comparable symmetry. We remark that while some mild robustness under growth is generally independent of growth rates, it is also a small effect that depends on the initial data and the nature of the domain itself, and is in a sense only a weak qualitative effect. Overall, we find that formation of spot patterns in Schnakenberg is qualitatively independent of growth altogether, at least in two spatial dimensions. For all further plots, we omit any simulation which is qualitatively identical to one already presented.
We give analogous plots of the elliptical and circular domains in Figs. 3 and 4, respectively. Similar to the case of spot patterns, the patterns observed in the ellipse in Fig. 3 were essentially insensitive to growth of any rate or functional form, but unlike spots, these patterns were qualitatively different to simulations on a static domain. In particular, we suspect that the thin static ellipse gives a proclivity to vertically Fig. 2 Simulations of Eq. (25) with kinetics given by (44) using parameters S1. We plot the final distribution of u on a the static domain, b a domain with linear isotropic growth over the timescale T g = T * /2 = 2500, c a domain with linear isotropic growth over the timescale T g = 10T * = 50,000, and d staged growth over the timescale T g + T g = 2 × 10T * = 100,000 (Color figure online) symmetric stripes forming, whereas growth at any rate will distort this structure. As before, the functional form of growth did not matter much in any case, but the growth rate and level of isotropy played a much larger role in the circular domain, shown in Fig. 4. Here, isotropic growth which was fast (Fig. 4b) gave rise to symmetrical patterns. Slow isotropic growth (Fig. 4c) gave more circularly symmetric patterning than in the static domain (Fig. 4a), but we note that this is not perfectly circular patterning, so that there is still some sensitivity to the initial data. Finally, anisotropic (staged) growth (Fig. 4d) gave rise to patterns indistinguishable from simulations on a static domain; note in particular that here the pattern is almost a rotation of that shown in Fig. 4a.
Turning now to simulations of the Gierer-Meinhardt kinetics (47), we plot some example simulations in Fig. 5. As shown in Fig. 2, the case of slow isotropic growth gives rise to some symmetry to the final patterned state, but this was not observed for either staged growth or faster growth regimes, which gave the final spot patterns comparably disordered to the static growth case. We note that in the case of the final domain being elliptical, all growth scenarios gave rise to comparably disordered final patterns. We note the inclusion of the constant term (feed rate) in Eq. (47) leads to generally robust dynamics. If this term was not present, then after the initial production of spots, changes to the manifold do not lead to spot splitting or production of new spots. Growth in the case without this feed rate will, in all cases we explored, lead to transport of spot solutions, but no meaningful interactions or new behaviors.
We also consider Gierer-Meinhardt simulations in the case of small isolated spot solutions (parameter set GM2 from Table 1) and plot an example in Fig. 6. We note that the static domain gives rise to many more spots than the case of isotropic growth. In fact, staged growth or isotropic growth at all rates exhibits the same number of spots in analogous configurations as shown in the case of isotropic growth. We suspect this is due to a conservation of spots inherent in the system, as discussed above, despite a positive feed rate. This occurs here as the inhibitor is able to suppress the constant production of activator throughout the domain (due to a much larger diffusion rate of the inhibitor and a smaller degradation rate). New spots are formed as the domain grows, but typically on the boundary and far from existing spots (which generate large regions of inhibition), so that the final number of spots in all growth cases is 12-14, compared with the 31 present in the static domain. This quantitative difference was consistent across 5 different realizations of the initial conditions (with different numbers of spots forming, but always less than 50% in the case of any growth). A similar effect was also observed in the elliptical domain (not shown).
Finally, we consider the FitzHugh-Nagumo kinetics given by (50) on the elliptical domain. For these kinetics, we observe substantial differences due to the growth of the domain, and these depend heavily on the kind of growth and marginally on the rate of growth; as in all other cases, the functional form of the time-dependent growth was equivalent to just using linear growth at a different rate. In Figs. 7 and 8, we show the evolution under isotropic and anisotropic growth of spot solutions to these equations. In both cases, after T * = 2000 units of time, the solutions are identical as shown in Fig. 7a. We note that in both cases, the number of minimal (blue) regions before growth and the number after is the same (7), but that the spots in the static domain are stretched and curved in complex ways, highly dependent on the mechanism of growth, as shown. Compared to the static circular domain in Fig. 7f, this is a strikingly different kind of patterning. We also note that the final patterned state here is stationary and robust against small numerical perturbations. We note that different initial conditions can tend to different long-time distributions, but these have comparable structures (e.g., in the staged growth case we always observe a central curved blue region). Finally, we remark that this conservation of minimal regions is only true for sufficiently slow growth, as the case of isotropic and staged growth shown in Fig. 9 demonstrates. We do note that there is still some qualitative difference between the kinds of growth in this case (namely, the staged growth has patterns more clearly oriented in the y direction), but the effect is weaker due to the introduction of more minimal (blue) regions.
Lastly, we consider the evolution of labyrinthine solutions under different growth regimes for the FitzHugh-Nagumo kinetics. In Figs. 10 and 11, we give examples of this on the circular domain and the elliptical domain, respectively. As in the spot simulations, it is clear that growth can change the qualitative properties of the patterns, and that this can depend on both if the growth is isotropic or anisotropic. In all cases, we observe labyrinthine-like patterns, but see that slow isotropic growth (Fig. 10c) can induce circular striping in these patterns, and that slow staged growth (Fig. 10d) can also introduce alignment along the most recent growth direction (the y axis). All Fig. 7 Simulations of Eq. (25) with kinetics given by (50) using parameters FHN1 for isotropic growth on a circular domain. We plot several distributions of u on a domain with linear isotropic growth over the timescale T g = 10T * = 20,000, at times a T * (fully developed pattern before growth), b 3T * , c 4T * , d 7T * , e T = 12T * (the final distribution), and f the same simulation on a static domain. Also note that the domains are not shown to scale; a has a radius that is 4 times smaller than in e or f (Color figure online) Fig. 8 Simulations of Eq. (25) with kinetics given by (50) using parameters FHN1 for staged growth first in the x-axis only then in the y axis only. Note that at time T * , the solution is identical to Fig. 7a. We plot several distributions of u on a domain with linear growth over the growth timescales of T g = 10T * = 20,000 in each direction. These are shown at times a 11T * (fully grown in the x-direction), b 12.5T * , c 14T * , and d the final patterned state (Color figure online) Fig. 9 Simulations of Eq. (25) with kinetics given by (50) using parameters FHN1 for both staged growth and isotropic growth in a fast growth regime. Note that at time T * , the solution is identical to Fig. 7a. We plot several distributions of u on a domain with linear growth over the growth timescales of T g = T * /2 = 1000 in each direction (a) and isotropically (b) (Color figure online) Fig. 10 Simulations of Eq. (25) with kinetics given by (50) using parameters FHN2 for various kinds of growth. We plot several distributions of u on a domain with linear growth over various growth timescales. These are the a static domain, b isotropic growth over T g = 0.5T * , c isotropic growth over T g = 10T * , and d staged growth over T g = 10T * first in the x direction and then in the y direction (Color figure online) rates of growth changed the final elliptical state in qualitatively similar ways as shown in Fig. 11, in that both maximal and minimal regions align along the semimajor axis of the ellipse, which is also the only axis of growth.

Reaction-Diffusion Systems on Ellipsoidal Domains
We now demonstrate solutions to Eq. (41) for the various reaction kinetics and growth regimes discussed. We note that now with three axes of possible growth, that there are three final possible domains (an oblate and a prolate spheroid or ellipsoid) where two axes are of equal size, depending on if one axis is smaller or larger than the other two. We did not systematically explore all intermediate domains with all variations  (44) using parameters S1 over a growth timescale of T g = 10T * . We plot the distribution of u on a the domain before growth, b the domain after growth along the x axis (T * + T g ), c the domain after growth along the x and y axes (T * + 2T g ), and d the distribution of u on the final manifold (Color figure online) in growth or kinetics, but we did include some in order to determine the influence of anisotropic curvature (as opposed to the isotropic curvature of a perfect sphere) on the final pattern. These effects can be seen during the transient staged growth, despite the shown states not being completely stationary, so we omit these results for brevity. Similarly, in all cases the effects of growth on different kinetics persist in these simulations, so we only show simulations where the influence of curvature (as well as new directions of anisotropic growth) can affect the final pattern. In particular, we do not include simulations of Gierer-Meinhardt, although we did observe the same kinds of patterns as shown in Figs. 5 and 6.
We begin by demonstrating the influence of curvature on Schnakenberg spot solutions in Fig. 12 for slow staged growth, which are analogous to those in Fig. 1. We see that the intermediate domains have spots which are not perfectly circular; as men-tioned before, the patterns do not change much if the growth is stopped and the pattern is allowed to settle into a stationary solution on either of these non-spherical domains. Despite the manifold being locally planar (by definition of a two-dimensional manifold), the spots are large enough to be influenced by local curvature and anisotropy so that they do not resemble their planar analogues. Despite growth, however, the final patterned state on the larger sphere is qualitatively the same on static or grown manifolds, independent of any details of the growth, as in the planar case shown in Fig. 1.
We next consider the evolution of labyrinthine Schnakenberg patterns in Fig. 13. Here we compare staged growth of these patterns in fast and slow regimes. The transient behaviors in the fast and slow growth regimes are quite different, but the final steadystate patterns are comparable with the fast growth case exhibiting some minor defects in the observed final labyrinthine pattern. We also note that the final patterns in the static domain and the slow growth domain, as shown in Fig. 13b, h, respectively, are essentially rotations of one another. Similar behavior is observed in the isotropic growth cases, where the final steady state differs from the static domain in the same way that different random initial data would (not shown).
We now consider the evolution of spot solutions using FitzHugh-Nagumo kinetics, which we show in Fig. 14. We note comparable behavior to the planar case shown in Figs. 7 and 9. Initially, three spots form on the domain of radius r 0 , whereas the static domain with radius 4r 0 has many more spots. Slow growth from the initial sphere, however, preserves these regions topologically (e.g., no new regions are formed), but they are twisted into curved patterns, which depend heavily on the directionality of growth and the domain. As in the planar case, growth which is sufficiently fast can induce the creation of new regions, leading to four minimal regions in staged growth over the growth timescale T g = T * /2 (for each direction) and eight minimal regions in isotropic growth over the same growth period. We recall that the staged growth occurs in each direction for a period of T g time units, so it is reasonable for fast isotropic growth in this case to produce more splitting, as more surface area is created in a shorter period of time. We note, however, that there is some sensitivity to exact geometries and growth rates that are sufficiently fast, so we do not pursue any further quantitative claims here.
Finally, we again consider the evolution of labyrinthine patterns under growth in Fig. 15. We again observe interesting patterns on the intermediate domains, but the effect of anisotropic growth is no longer obvious, as it was in the planar case shown in Figs. 10, 11 and 12. We also note that, independent of the growth rate or the use of static domains or isotropic growth, these figures are indicative of the kinds of patterns found in this case, and these variations do not seem to influence the qualitative behaviors. In all cases, a similar pattern emerges that appears qualitatively as one might expect when studying labyrinthine patterns on curved geometries, and growth only (qualitatively) seems to change one pattern into another without any stark differences due to the growth itself.  (50) using parameters FHN1. We plot the distribution of u on a the domain before growth, b the final static domain, and grow the domain using staged growth over the timescale T g = 10T * in c-e where c is at T * + T g , d is at T * + 2T g , e is the final stationary pattern in the staged growth case, and f is the isotropic growth case over the same growth period as in each staged direction (Color figure online)

Stripe Stability in Isotropic and Anisotropic Growth Regimes
Having systematically explored a wide range of reaction kinetics and growth regimes on planar and curved domains, we now briefly mention some results concerning the numerical stability and evolution of target and stripe patterns under growth. These patterns are special solutions that typically require spatial heterogeneity or specific initial data in order to observe (although we note that growth can do this as well, as shown in Figs. 4b, 13f). There is a large literature on the applications of stripe patterns in real Fig. 15 Simulations of Eq. (41) with kinetics given by (50) using parameters FHN2. We plot the distribution of u and grow the domain using staged growth over the timescale T g = 10T * where a is before growth at time T * , b at time T * + T g , c at time T * + 2T g , and d is the final stationary pattern (Color figure online) organisms and chemical models (Ouyang and Swinney 1991;Nagorcka and Mooney 1992;Kondo and Asai 1995) (and this is something that Turing mentioned in his original 1952 paper), as well as a large literature on their stability and other theoretical properties (Ermentrout 1991;Lyons and Harrison 1991;Moyles et al. 2016).
In the case of the planar ellipse, we will consider initial data of the form where u * is again the spatially homogeneous steady state, α ∈ [0, 1] is the radial coordinate on the unit disk, and k is a positive integer which determines the number of circular stripes in the target pattern. The vector of opposite signs is used as we will only consider the Schnakenberg kinetics with parameters that give rise to stripe solutions (S2), and the activator and inhibitor are always out of phase as is typical for models with cross-kinetics. Depending on the size of the domain, these kinetics and parameter choices lead to a stable pattern very close to the original one, with k − 1 interior stripes, one stripe along the boundary, and a single central spot completing the target pattern. One can readily generate a bifurcation diagram in radius r 0 and k and numerically observe stability of different modes, but for brevity we set k = 2 and r 0 = 20 in the elliptical case, and consider growth up to r T = 4r 0 = 80. We plot the stationary pattern on the initial domain, as well as a simulation using isotropic growth in Fig. 16. We note that simulations using k = 2 on the static final domain become unstable and evolve into the same pattern shown in Fig. 16b, which corresponds to k = 8 without a central spot. These patterns are stable for growth timescales larger than 10 4 , but as shown in Fig. 16c, defects emerge and destroy the perfectly symmetric rings if the growth is too rapid. Likely there is an interplay here between numerical discretization errors, errors in the initial data, and symmetry of growing target patterns, but we leave further exploration of these intertwined effects to future work. We also find that staged or anisotropic growth of any form appears to destroy symmetric stripe patterns, in all simulations we have considered.
Lastly we consider analogous patterns on the ellipsoid. We use a similar form of the initial data as in the case of the ellipse, but project it along one axis, so that the initial pattern develops into targets with surrounding stripes at antipodal points of the y axis. In the extrinsic and stationary Cartesian frame, this is given by where x ∈ [−1, 1]. In Fig. 17, we show the stationary pattern on the initial domain, the final static domain, and the case of slow isotropic growth, as shown in Fig. 16b. We note that the static domain also breaks up into a higher-mode pattern, but growth distorts the pattern by introducing defects throughout the striped patterns. This is likely because isotropic growth of a sphere is not isotropic with respect to the striped pattern, so growth destabilizes these stripes, leading to labyrinthine patterning. This is inherently an effect due to the curvature of the domain not present in planar isotropic growth.
For staged growth, we show some example simulations in Fig. 18 which all correspond to different stages and directions of anisotropy given the same initial state shown in Fig. 17a. Independent of the initial direction of growth, we see staged growth develop into an approximate target pattern for an elongated ellipsoid, but further growth along other axes leads to more complicated labyrinthine patterning, leading to the final patterns comparable to 17c. Compared with Fig. 16, the target patterns on curved domains appear much less stable to growth, likely due to the influence of curvature.

Discussion
Extending the work of Plaza et al. (2004), we have generalized reaction-diffusion systems on growing two-dimensional manifolds to include situations where the growth is anisotropic yet dilational in nature. This allowed us to systematically study the evolution of patterns in a variety of growth scenarios, in order to understand how growth affects the robustness of non-equilibrium patterns. We primarily concentrated on patterns which naturally emerged from Turing instabilities, but also gave some examples of stripe patterns, and their sensitivity to growth and curvature. While the restriction to only dilational anisotropy is somewhat restrictive, it allows for an exploration of the role that anisotropy has on the evolution of patterns under growth in a wide variety of growth scenarios.
Our results show several key findings. Firstly, the growth rate and kind of growth (e.g., isotropic or anisotropic) can significantly change the final pattern, even though we can demonstrate that the final pattern is a steady state. This suggests that the high level of multistability in these systems on fixed domains can be heavily influenced by their history due to growth, as suggested by Klika and Gaffney (2017). However, the functional form of growth was never really important; all simulations using different forms of the growth rates were qualitatively the same with linear growth at a suitable rate, as long as the initial and final manifolds were fixed. Secondly, we have used anisotropy in staged growth as compared to both quasi-static growth and isotropic growth to the same final manifold, as a concrete way of demonstrating the importance of both anisotropy and history in determining patterning. A final key finding is that the qualitative nature of patterns heavily depends on both the nonlinear reaction Fig. 18 Simulations of Eq. (41) with kinetics given by (44) using parameters S2. We plot the distribution of u. In a, c, e, we show staged growth along axes x, y, and then z, whereas in b, d, and f we show staged growth along axes y, x, and finally z. These are at times a, b 6.6 × 10 5 , c, d 1.26 × 10 6 , and e, f 1.92 × 10 6 (Color figure online) kinetics under consideration and the specific parameter regimes studied. While we only explored three candidate choices of reaction kinetics, we believe these are good exemplars of reaction-diffusion systems sufficient to obtain a qualitative picture of the influence of growth on patterned states.
Spot patterns in Gierer-Meinhardt and Schnakenberg were often not influenced by growth (Figs. 1, 2, 5) and manifold curvature only gave rise to slight local distortions in spots (Fig. 12b-c). The conservation of spots in the Gierer-Meinhardt system is something that has not, to our knowledge, been reported in the literature. Specifically, dropping the constant term in the first of Eq. (47) and letting spot patterns stabilize before growth always led to no new spots forming, independent of the manner of growth. With this feed rate term, however, the behavior of the Gierer-Meinhardt system then became heavily dependent on both the parameters used, as well as the growth rate. 'Large' spot solutions behaved as in Schnakenberg, with comparable numbers and distributions of spots in both quasi-static and grown domains. The smaller spot solutions, in Fig. 6, showed a dependence on the growth rate wherein slow enough growth permitted fewer final spots in the final domain. Finally, spot solutions arising from the FitzHugh-Nagumo kinetics exhibited an interesting dependence on both isotropy and curvature (see Figs. 7,8,9,14). Isotropic growth led to labyrinthinelike patterns, whereas staged growth had an obvious effect on the curved regions. For slow growth rates, the total number of such regions was the same between the initial domain and the final one, whereas fast growth could initiate the formation of new regions, akin to spot splitting due to growth. This stark contrast between these reaction kinetics suggests that excitable media might behave quite differently under growth, compared to their dissipative counterparts.
Labyrinthine patterns exhibited a similar richness of behavior. Using Schnakenberg kinetics, stripes could be aligned if the domain was static, but not symmetric, whereas growth would destabilize such alignment (Fig. 3). For some growth rates and kinds of growth, these solutions could spontaneously form target patterns with symmetric stripes (Figs. 4b, 13d). Solutions for FitzHugh-Nagumo kinetics displayed slightly different behaviors, but also showed a strong dependence on anisotropy (Figs. 10,11,15). In this case, we did not see the emergence of target patterns or similarly perfectly symmetric structures, but we suspect that for some parameters this might be possible.
Finally, we explored the evolution of striped patterns in Sect. 6. These findings reinforce the key roles of growth rates, anisotropy, and curvature in determining the properties of the final steady-state patterns. It is interesting that in any growing ellipsoid we do not find any stable target or striped pattern that is not a labyrinthine-like pattern. This may be due to the specific parameter choices we have made, and we leave further investigation of stripes and other specific kinds of patterns to future work.
We have demonstrated the complex influence of growth, curvature, and anisotropy on non-uniform patterns arising from Turing instabilities on planar and curved manifolds. There are many possible extensions to our work, including considering manifolds with holes or cusps, or considering bulk-surface reaction-diffusion systems as in Madzvamuse et al. (2015), under the influence of growth. Additionally, further generalizing the reaction-diffusion framework to account for arbitrary nonuniform growth would allow for more accurate modeling of realistic pattern formation in development and elsewhere. Finally, our work also shows behaviors of non-uniform solutions to nonlinear PDE in complicated domains, and the stability and evolution of such solutions are usually confined to illustrative special cases, as in Trinh and Ward (2016). One could pursue similar analyses in the context of anisotropic growth, in order to deepen our understanding of the effects we have demonstrated here.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.