Domain-growth-induced patterning for reaction-diffusion systems with linear cross-diffusion

In this article we present, for the first time, domain-growth induced pat- tern formation for reaction-diffusion systems with linear cross-diffusion on evolving domains and surfaces. Our major contribution is that by selecting parameter values from spaces induced by domain and surface evolution, patterns emerge only when domain growth is present. Such patterns do not exist in the absence of domain and surface evolution. In order to compute these domain-induced parameter spaces, linear stability theory is employed to establish the necessary conditions for domain- growth induced cross-diffusion-driven instability for reaction-diffusion systems with linear cross-diffusion. Model reaction-kinetic parameter values are then identified from parameter spaces induced by domain-growth only; these exist outside the classical standard Turing space on stationary domains and surfaces. To exhibit these patterns we employ the finite element method for solving reaction-diffusion systems with cross-diffusion on continuously evolving domains and surfaces.


1.
Introduction. Understanding of biological processes during growth development is an unresolved issue in developmental biology that is only starting to be addressed recently. Introducing domain growth into the modelling results in nonautonomous systems of partial differential equations whose analytical tractability is not well understood [8,20,27,28,32]. In the area of developmental biology, partial differential equations for pattern formation take the form of reaction-diffusion type [38,47]. On stationary domains, for example, it is well known that one major criticism of reaction-diffusion theory for pattern formation is the tight control of the model reaction kinetic parameter values [2]. Underpinning this theory is the concept of diffusion-driven instability which leads to patterns that are stationary in time and periodic in space. For a two-component reaction-diffusion system, a key requirement for diffusion-driven instability is the concept of long-range inhibition and short-range activation [18]. This implies that one of the species, the inhibitor, must diffuse faster, typically much faster, than the autocatalytic species, the activator, thereby fulfilling one of the necessary conditions for the formation of spatial structure.
Several generalisations of the reaction-diffusion theory for pattern formation have been undertaken in order to relax some of these constraints. One of these generalisations involve the introduction of domain growth [8,27,28,32]. It is well-known that many problems in biology involve growth. In [32] we proved that in the presence of domain growth, it is no longer necessary to restrict reaction kinetics to an activatorinhibitor type; a long-range activation and short-range inhibition and/or activation chemical processes are all capable of giving rise to what we termed domain-growth induced diffusion-driven instability. However, this generalisation still requires that the inhibitor must diffuse much faster than the activator species and therefore equal diffusion coefficients do not give rise to the formation of spatial structure during growth.
To further relax these constraints on the diffusion coefficients, in our most recent studies, we introduced cross-diffusion which further enhanced the pattern formation mechanism as well as proving a larger range of intervals for the diffusion coefficients [33,35]. In many multi-component systems, there are various forms of diffusion depending on the bio-chemical problem at hand [48]. Cross-diffusion is characterised by a gradient in the concentration of one species inducing a flux of another chemical species. In molecular biology, cross-diffusion processes appear in multicomponent systems containing at least two solute components [37,49]. Multicomponent systems containing nanoparticles, surfactants, polymers and other macromolecules in solution play an important role in industrial applications and biological functions [37]. In developmental biology, recent experimental findings demonstrate that crossdiffusion can be quite significant in generating spatial structure [48]. The effects of cross-diffusion on models for pattern formation, i.e. reaction-diffusion type, have been studied in many theoretical papers [5,6,16,17,22,23,41,42,46,50,52,53]. Recently, in Madzvamuse [33] we showed that introducing cross-diffusion to a system of reaction-diffusion equations results in further relaxation of the conditions necessary for the emergency of patterns. In particular, an inhibitor and activator or two activators can diffuse at equal rates, however, the product of the rates of the principal diffusion coefficients must be greater than the product of the crossdiffusion rates. For detailed theoretical analytical and computational results on the effects of domain growth on pattern formation, the interested reader is referred to results published in [8,20,24,25,26,28,32,36]. In this article our focus is to show-case how domain growth induces patterning in the absence and presence of cross-diffusion.
Recently, there has been an increase in the development of numerical methods for approximating solutions of partial differential equations posed on evolving domains and surfaces. Examples include (but are not limited to) the moving grid finite element method [28,29,30], the method of lines [7], finite differences and finite volumes [10], evolving surface finite element methods on triangulated surfaces [3,11,14,34], implicit finite element methods using level set descriptions of the surfaces [12,13,39,45], diffuse interface methods of which phase-fields are an example [9,15], particle methods using level set descriptions of the surface [10,19,21] and closest-point methods [25,26]. In this article we choose to implement the moving and evolving surface finite element method on growing domains and surfaces, respectively. Our approach is based on the finite element method which has the capability of handling complex geometries and shapes [11,12,13,15]. The finite element method is based on the approximation of the domain or surface by a triangulated domain or surface consisting of a union of triangles. For surfaces, the vertices of the triangulation lie on the continuous surface. On evolving domains, given the boundary evolution law, the internal mesh is deformed using a spring analogy [4,28]. On evolving surfaces, we triangulate the surface and approximate the system of partial differential equations using piece-wise linear surface finite element spaces based on the triangulation. In this instance the vertices of the triangulation are moved with a velocity which is either prescribed or is governed by some evolution law. For both numerical methods, we formulate an appropriate conservation law. The finite element methodology exploits the special features of this conservation law when written in an appropriate variational form.
Hence, the focus of this paper is to study the role of domain growth in the generation of pattern formation in the absence and presence of cross-diffusion for systems of reaction-diffusion equations. To the authors' best knowledge, the role of domain growth whereby parameter values are selected from the domain-growth-induced Turing parameter space has not been studied and is a novelty on its own. When coupled with cross-diffusion, patterning enhancement is observed. For analytical purposes, we study uniform isotropic growth, non-uniform growth is beyond the scope of this manuscript, since the analysis for such growth requires more complex mathematical techniques. For illustrative purposes, we exhibit patterns generated by well studied growth functions such as linear, logistic and exponential. These have been shown to be appropriate for describing some growth processes in biolgy [38].
Our paper is therefore structured as follows: in Section 2 we present the model equations posed on evolving domains and surfaces and these consist of a system of reaction-diffusion equations with linear cross-diffusion. Domain and surface evolution terms are modelled through dilution and convective terms. Section 2.3 summarises the conditions for domain-induced cross-diffusion-driven instability on planar domains and surfaces. Time-dependent parameter spaces for the exponential, linear and logistic growth functions are exhibited for the case of evolving domains in Section 3. In Section 4 we present finite element solutions confirming theoretical predictions and these illustrate pattern formation emerging only due to domain and surface evolution. Such patterns do not exist in the absence of domain or surface evolution with or without linear cross-diffusion. Conclusions and discussion on future directions are detailed in Section 5.
2. Models posed on evolving domains and surfaces.
2.1. Notation and preliminaries. In this section, we establish some notation to be used throughout the paper. If Γ ⊂ R 3 is a two dimensional hypersurface with unit normal ν we denote the tangential gradient of a scalar function u by whereū is a smooth extension of u into an open neighbourhood of Γ and ∇ denotes the usual gradient in R 3 . Hence, we define the Laplace-Beltrami operator as the tangential divergence of the tangential gradient ∆ Γ u := ∇ Γ · ∇ Γ u.

ANOTIDA MADZVAMUSE AND RAQUEL BARREIRA
If Γ(t) is an hypersurface in R 3 evolving in time, t, according to a velocity field v, we denote the material derivative of u by 2.2. Reaction-diffusion systems with linear cross-diffusion posed on evolving domains and surfaces.
2.2.1. Reaction-diffusion system with linear cross-diffusion posed on evolving domains. Let Ω(t) ⊂ R m (m = 1, 2, 3) be a simply connected bounded evolving volume for all time t ∈ I = [0, t F ], t F > 0 and ∂Ω(t) be the evolving boundary enclosing Ω(t). Also let u = (u (x, t) , w (x, t)) T be a vector of two chemical concentrations at position x ∈ Ω(t) ⊂ R m . The evolution equations for reactiondiffusion systems with linear cross-diffusion can be obtained from the application of the law of mass conservation in an elemental volume using Reynolds transport theorem. The growth of the volume Ω(t) generates a flow of velocity v to yield the non-dimensional system [33,35] u(x, 0) = u 0 (x), and w(x, 0) = w 0 (x), x on Ω(0), where ∇ 2 is the Laplace operator on domains and volumes, d is the ratio of the diffusion coefficients and d u and d v are the ratios of the cross-diffusion. Here, n is the unit outward normal to Ω(t). Initial conditions are prescribed through nonnegative bounded functions u 0 (x) and w 0 (x). In the above, f (u, w) and g(u, w) represent nonlinear reactions.
2.2.2. Reaction-diffusion system posed on evolving surfaces. Suppose now that the domain is an evolving hypersurface, Γ(t), with velocity v = V ν, where ν is the unit outward pointing normal vector to Γ(t), and V is the normal velocity. In this study, we assume that there are no contributions to the surface evolution law from the tangential component of the flow velocity. Without loss of generality, let u = (u (s, t) , w (s, t)) T be a vector of two chemical concentrations at position s ∈ Γ(t) ⊂ R m+1 (m = 1, 2) . Then adapting the derivation of the reaction-diffusion equations on an evolving surface in [3], taking into account the terms corresponding to linear cross-diffusion, we obtain the following non-dimensionalised system of reaction-diffusion equations [33,34]  u(s, 0) = u 0 (s), and w(s, 0) = w 0 (s), s on Γ(0), (2.2) where we are assuming Γ(t) is a compact smooth connected and oriented evolving hypersurface in R m+1 , with m = 1, 2. As for the boundary conditions, we impose µ · ∇ Γ(t) u = µ · ∇ Γ(t) w = 0, where µ is the co-normal vector to ∂Γ(t), if ∂Γ(t) is nonempty. Otherwise, if the surface has no boundary, then there is no need for boundary conditions. Remark 2.1. Notice that for the planar case the tangential derivative coincides with the classical derivative and in that case the model equations become exactly the same as (2.1). Furthermore, we note that domains are a special case of evolving surfaces with appropriate boundary conditions.

2.2.3.
Activator-depleted reaction kinetics: An illustrative example. For illustrative purposes, we consider a specific type of nonlinear reactions, namely, the activatordepleted model also known as the Brusselator or Schnakenberg model [18,40,44] with a, b ∈ R + .

2.2.4.
Domain and surface evolution law. In many biological problems, volume and surface evolution laws are modelled from experimental observations. In general, the model equation for the flow velocity v is a complex nonlinear function involving the chemical species. In order to carry out some analytical studies, we assume that the volume and surface evolution law is governed by a uniform isotropic growth, independent of the chemical species, and is defined by Under this form of growth, we can then define the flow velocity by As a result, it follows then that we can define Here m denotes the spatial dimension of Ω(t).

Remark 2.2. On surfaces it can be shown that
In the above framework, the growth function ρ(t) is explicitly given and examples include linear, logistic and exponential as well as piece-wise continuous functions. To proceed, we map for all time t, the model equations posed on evolving volumes and surfaces to models posed on stationary reference domains and surfaces as shown below [20,32,35].

2.2.5.
Lagrangian transformation on evolving domains. We summarise how the model system is mapped from a continuously growing domain to a reference or initial stationary domain. A similar derivation suffices for evolving surfaces; we will omit the derivation on evolving surfaces for the sake of saving space. Assuming spatially linear and isotropic evolution of the domain, can can write It can be easily shown that the reaction-diffusion system with linear cross-diffusion (2.1) can be mapped to the following non-autonomous partial differential equations (with time-dependent coefficients) [32,35] where h(t) as defined by (2.6) on planar domains. Similarly for the boundary and initial conditions.
A similar transformation suffices on evolving surfaces, details of which are omitted.

2.3.
Domain-growth induced driven instability for reaction-diffusion systems with linear cross-diffusion. From here onwards, we will drop the hatsf or ease of exposition. Next, we describe the timescales on which the mathematical analysis for the conditions for domain-growth induced parameter spaces is based. These timescales are purely for mathematical analysis and have no bearing on the numerical methods for the generation of these spaces nor on the numerical integrators of the ODE model system. To proceed, it must be noted that in most biological systems, growth is typically driven by cell proliferation, and occurs on the timescale of the cell cycle duration (24 hours or greater). While on the other hand, biochemical reactions generally occur on a much faster time scale of seconds to minutes, and are typically constrained by the diffusion processes. Diffusion and cross-diffusion coefficients are known to be of the orders of 10 −6 cm 2 s −1 to 10 −5 cm 2 s −1 , then the timescale for the diffusive and cross-diffusive processes can be shown to be between 40 seconds to 170 minutes [48]. Given that the length scale associated with a volume of 10 4 to 10 6 cells is about 2 × 10 −4 m to 10 −3 m, then we are in a position to identify a small oder parameter where T g denotes the growth timescale, and T dyn the maximum timescale of the diffusion and cross-diffusion processes and the biochemical reactions. Following [32], it can be shown that for the above range of estimates ∈ [5 × 10 −4 , 0.12].
In the following, we study domain-growth induced driven instability analysis for reaction-diffusion systems in the presence of linear cross-diffusion. We focus on spatially linear, isotropic evolution of the domain. Given that model system (2.9) is a non-autonomous system of partial differential equations with time-dependent coefficients, the model system does not always admit a uniform steady state. To this end, we will consider a time-dependent manifold (i.e. spatially independent) (u S (t), w S (t)) defined as follows.
Definition 2.1. [A spatially independent time-dependent manifold] A spatially independent manifold defined by (u S (t), w S (t)) is a solution of (2.9) if it solves the non-autonomous system of ordinary differential equations (2.11) In the above h(t) is as defined in (2.6) (or (2.7)) on planar domains (or on surfaces). [Domain-growth induced driven instability on evolving domains for reaction-diffusion systems with linear cross-diffusion] Domain-growth induced driven instability on evolving domains occurs when a spatially independent manifold (u S (t), w S (t)), linearly stable in the absence of spatial variations (diffusion and cross-diffusion terms) becomes unstable in the presence of spatial variations.
The analysis of the solutions of the non-autonomous systems of ordinary differential equations is identical to that derived and presented in Madzvamuse et al., [32,35], and for the sake of brevity, we will omit details of the analysis. For detailed analysis of the model system, we refer the interested reader to consult [32,35]. Following Madzvamuse et al., [32], we consider a perturbation at time t 0 , which can be distinct from the initial time, and how it has evolved by time We are now in a position to state the conditions for domain-growth induced driven instability for reaction-diffusion systems with linear cross-diffusion. For this, it will be useful to define t 1 : for use below.
Theorem 2.1. The necessary conditions for a cross-diffusion driven instability (accurate to first order in 1 defined in (2.10)) when considering the system at time t due to a perturbation at time t 0 < t in the presence of spatially linear and isotropic growth as in system (2.9), are given by In the above, the subscripts u, w denote partial differentiation, with the Jacobian components f u , f w , g u and g w evaluated in terms of u S (t 1 ). Strictly, we have h * (t 1 ) ∼ O( ), due to the fact that the growth timescale is much longer than any other timescale so that one may self-consistently neglect O(h 2 * ) corrections in the above.
Proof. See [32,35] Remark 2.5. It must be observed that conditions (2.12)-(2.16) hold true for evolving surfaces with an appropriate h(t). The derivation of these conditions on evolving surfaces is identical to that presented for growing planar domains, the only difference is that on surfaces, we are looking for wavenumbers and eigenfunctions associated with the Laplace-Beltrami operator while on domains, these are associated with the standard Laplace operator.

Numerical generation of domain-growth induced parameter spaces.
To illustrate our theoretical findings, we briefly introduce slow, isotropic growth for exponential, linear and logistic growth profiles of the domain, calculating functions associated with domain evolution such as h(t) required for the subsequent explorations. Our goal is to compute and exhibit parameter spaces induced by domain evolution for a reaction-diffusion system with linear cross-diffusion. We will show that linear cross-diffusion coupled with domain evolution substantially alter classical parameter spaces as well as exhibiting new spaces that emerge only in the presence of domain growth. Next, let us assume that the domain growth is spatially linear and, in higher spatial dimensions, isotropic. Without any loss of generality, let us assume that t 0 = 0. As defined before, let x(t) = ξρ(t) where ρ(t) > 0 is the domain growth function satisfying ρ(0) = 1. As illustrated before, it can be easily shown that where m defines the number of spatial dimensions. In Table 1 we show the corresponding h(t) for linear, exponential and logistic growth functions, with κ = 1. The definition of A in this table ensures that ρ(0) = 1.
Lemma 3.1. Let us assume that the domain growth is spatially linear and, in higher spatial dimensions, isotropic. Then for linear and logistic domain evolution profiles, the function is a monotonically decreasing function whose limit tends to zero, i.e. The exponential growth, h(t) = mr which is constant for all time.
Proof. The proof follows clearly from the definitions of h(t) shown in Table 1.

DOMAIN-GROWTH-INDUCED PATTERNING FOR REACTION-DIFFUSION SYSTEMS 2783
3.1. Exponential, linear and logistic domain evolution. Next we compute parameter spaces generated by exponential, linear and logistic growth functions for two-dimensional evolving domains. For all these growth functions we compute the time-dependent solutions (u S (t), w S (t)) which satisfy the differential equations where h(t) = ∇ · v and given initial values at time t 0 = 0. From the exponential growth, h(t) = mr; for the linear growth function, h(t) = m r rt+1 and similarly for the logistic growth function, h(t) = m r (κ − ϕ(t)) where ϕ(t) = Aκe κrt 1+Ae κrt , with A = 1 κ−1 . Unlike previous studies where parameter spaces where generated for the case of uniform steady states as functions of the growth rate r [32,35], here we want to generate parameter spaces given by the time-dependent manifolds solving the system of non-autonomous differential equations In all our results on parameter space generation, we fix the dimension m = 2 and the scaling parameter γ = 1. For this case, in the absence of domain growth and linear cross-diffusion, the parameter space reduces to the classical Turing diffusiondriven instability parameter space obtained in the case of stationary domains and surfaces [38,47]. However, to generate patterning, there exists a critical scale length γ c associated with a minimal domain length scale favourable for the emergence of pattern formation [32]. In the numerical simulations shown in Section 4, we select γ = 200 > γ c in order for patterning to occur. γ is a scaling parameter reflecting the strengths of the reaction-kinetics and is associated with the minimal domain scale necessary for patterns to emerge. Taking γ large results in an amplification of the continuously evolving parameter spaces with very little changes in topology (results not shown). The interested reader is referred to the work of Murray [38] for further details. In this section, we will compute parameter spaces given by the reaction kinetic parameters (a, b) when varying the diffusion coefficients d, d u , d v and the growth rate r. In the absence of growth, in the case that these spaces exist, they reduce to the well-known standard stationary parameter spaces, otherwise such spaces do not exist.

3.2.
Numerical procedure for the generation of domain-growth induced parameter spaces. To numerically generate the domain-growth induced parameter spaces given by conditions (2.12)-(2.16), for each point (a, b) in the positive twodimensional quadrant, the non-autonomous system of ordinary differential equations (3.1)-(3.2) must be solved numerically in time. For this, we use the MATLAB ode45 solver with some given positive initial condition. In Figure 1, we plot a series of solutions (u S (t), v S (t)) using MAPLE (for illustrative purposes only) to demonstrate the typically observed space-independent solutions in time. The MATLAB algorithm for the numerical generation of the parameter spaces is available on request. In many cases, these are oscillatory solutions giving rise to limit cycles given the activator-depleted model (2.3). The parameter values a and b are selected from the parameter spaces shown in Figures 5 to 12, and their values are as described in the respective captions of these figures. We have chosen these model kinetic parameters only from the spaces generated by domain growth; such parameter values are inadmissible for the case without domain growth (either with or without crossdiffusion). Our primary aim is to demonstrate the role of domain growth to the emergence of pattern formation for reaction-diffusion systems with or without linear cross-diffusion, and that cross-diffusion allows for a wider set of diffusion coefficients for patterning. These u S (t) and v S (t) time-dependent manifolds are stable in the absence of diffusion and cross-diffusion and become unstable when diffusion and/or cross-diffusion processes are introduced, to result in patterns forming as illustrated in Figures 5-14 (for both evolving domains and surfaces). 3.2.1. Parameter spaces generated by the exponential growth of the domain. In this first example, we study the generation of domain-growth induced parameter spaces driven by the exponential growth. Here, we fix the classical diffusion coefficient On the other hand, for each column, we vary the linear cross-diffusion coefficients d u and d v and fix the growth rate r. The colour code for each figure can be decoded as follows: blue represents the classical Turing parameter space in the absence of domain growth (computed at t = 0), green represents a snap-shot of the time-evolution of the exponentially generated parameter space observed at time t = 0.1 and red represents a further snap-shot of the time-evolution of the exponentially generated parameter space observed at time t = 0.3. It must be observed that the blue parameter space is unchanged throughout the individual figures since this represents the space with no growth (exponential or otherwise). We summarise the results shown in Figure 2 as follows: 1. Row 1: We fix the linear cross-diffusion d u = 0 and d v = 0 which entails that the parameter spaces in green and red are generated purely by the exponential growth in the absence of linear cross-diffusion. Here, by varying now the growth rate r = 0.01 (a), r = 0.04 (b) and r = 0.08 (c) we demonstrate the role of domain growth to the generation of parameter spaces during growth development. By varying the growth rate r, the exponentially-generated spaces evolve substantially from one time snap shot to another with an increase in the growth rate (almost vanishing for the growth rate r = 0.08 as indicated in Figure 2  It must be noted that, in general, for slower growth rates there is a substantial overlap of the parameter spaces while for larger values of the growth rate, less overlap of the parameter spaces is observed. The largest parameter spaces are obtained in the absence of growth (at time t = 0) and these are induced by diffusion and cross-diffusion only).
In Figure 3 we exhibit exponentially generated parameter spaces computed by choosing the principal diffusion coefficient to violate the necessary conditions for diffusion-driven instability in the absence or presence of domain growth by picking d = 1. In the absence of cross-diffusion, all the spaces exhibited here are nonexistent, these only emerge due to the presence of cross-diffusion. However, domain growth enhances the formation of these parameters spaces once they are formed as illustrated in Figure 3. 1. Row 1: In Figure 3 (a)-(c) we plot parameter spaces computed by taking d u = 1 and d v = 0.8. As before, the growth rate r is varied accordingly and plots are exhibited at times t = 0 (no growth), t = 0.1 and t = 0.3. Note that the blue parameter space remains unchanged throughout the evolution, this is the space generated by linear cross-diffusion in the absence of domain growth. The green and red spaces are snap-shots of the continuous time-evolution of this blue space once exponential growth is introduced. In all the three simulations in this row, the exponentially generated parameter spaces have very little overlap with the stationary initial cross-diffusion induced parameter space.
2. Row 2: By taking d u = 0.8 and d v = 1, Figure 3 (d)-(f) shows the evolution of the exponentially generated parameter spaces as the growth rate r is varied. The evolving parameter spaces remain almost unchanged as the growth rate is increased, indicating that domain growth plays a relatively less important role in the generation of potential candidate model parameters during growth development of this specific set linear cross-diffusion coefficients. Unlike in Row 1 there is substantial overlap between parameter spaces during growth development. 3. Row 3: In this numerical example, we introduce for the first time in this study, negative linear cross-diffusion in the v-equation of model (2.1) by taking d u = −0.8, and d v = 1, respectively. Figure 3 (g)-(i) shows that in the absence of domain growth (blue spaces), cross-diffusion induces parameter spaces in the form of a thin-stripe and these are much much smaller than those obtained in the presence of domain growth (compare to the green and red spaces). Clearly domain growth coupled with cross-diffusion enhances the mechanism for pattern formation.
Next we illustrate parameter spaces induced by linear cross-diffusion in the presence of linear and logistic growth function (see Figure 4) and compare such parameter spaces to those obtained with the exponential growth as depicted in  [38,47]. The cross-diffusion and domain-growth induced parameter spaces in this section are independent of these random initial conditions.
In order to compute the parameter spaces given by conditions (2.12)-(2.16) for the cases of linear and logistic growth functions, the non-autonomous differential equations (3.1)-(3.2) must be solved for every point (a, b) in the plane at any given time t thereby giving rise to time-dependent parameter spaces. See Section 3.2 for the description of the numerical procedure to generate these domain-growth induced parameter spaces for reaction-diffusion systems with linear cross-diffusion.
It must be noted that parameter spaces generated by the linear growth function are almost identical to those obtained by logistic growth for slower growth rates while these become largely different for faster growth rates. In Figure 4 we superimpose parameter spaces generated in the absence of growth (blue) to those obtained by the exponential (light blue), linear (green) and logistic (red) growth functions with diffusion coefficient d = 1 and cross-diffusion coefficients d u = 1 and d v = 0.8. Snap shots of the evolving parameter spaces, for each set of parameters and growth rates, are captured at three time levels, t = 0, t = 0.01 and t = 0.03, respectively. The results reaffirm our observations that exponential growth gives rise to substantially different parameter spaces to those obtained for the case of logistic and linear growth functions. These spaces become smaller as time passes, in contrast to those obtained for the case of logistic and linear growth profiles which become larger with time.  package deal.ii [1], while the evolving surface finite element solutions on evolving surfaces are computed using the triangulated finite element toolbox ALBERTA [43]. For both the moving and evolving surface finite element methods on evolving planar domains and surfaces, respectively, we employ a modified backward Euler method as described in [31]. The method could also be thought of as a modified implicit-explicit method where we treat the diffusion and all linear terms implicitly and semi-implicitly all the nonlinear terms (a single Picard iterate). For illustrative purposes, we take a unit square and a unit sphere as initial geometries for the computation and mesh generation. In the case of the moving finite element method on planar domains, we take the timestep τ = 10 −3 and approximately 100 2 degrees of freedom. On evolving surfaces, we take the timestep τ = 5 × 10 −4 and approximately 65, 536 degrees of freedom. For growing planar domains and evolving surfaces, the mesh connectivity remains unchanged throughout the numerical solution procedure of growth development.   Figure 5 (b)-(e) where spots transient to stripes and then to spots again. In the absence of domain growth (results not shown), no patterns are exhibited, instead, the system converges to some uniform state ( Figure 5 (a)). Patterns only start to emerge when the domain reaches some critical threshold such that conditions (2.12)-(2.16) are fulfilled. From here onwards, we omit to plot the numerical solution with no patterns (such as that shown in Figure  5  domain. Parameter values have been selected such that without growth, patterning does not occur. It is clear that domain growth enhances the formation of patterns which would otherwise not form in the absence of domain growth. We observe the formation of stripes, spots, and circular patterns during growth development, with transitions between them. In all our simulations, (i) switching off domain growth results in the disappearance of these patterns (results not shown) and (ii) patterning seems to be insensitive to the random initial conditions, patterning is driven and stabilised by growth.
In Figures 13-14, we illustrate pattern formation on evolving surfaces where the emergence of patterns is driven by growth. As before, we select model parameter values from the domain-induced parameter spaces and illustrate that domain growth induces the formation of spatial structure on evolving surfaces. The evolving sphere and saddle show the emergence of spot patterns during growth development. Numerical simulations on evolving surfaces seem to suggest that surface geometry does not play a pivotal role in the generation of cross-diffusion driven patterns. For example, comparing Figures 9 (on growing planar domains) and 14 (on an evolving saddle surface), we observe that by selecting identical parameter values, the same spot patterns are observed, with similar wavelengths and amplitudes. Further studies are necessary to confirm if this is a general phenomenon for the case of linear cross-diffusion.

5.
Conclusion and future directions. In this study, we have carried out a detailed characterisation of domain-growth induced parameter spaces (these are timedependent) for reaction-diffusion systems with linear cross-diffusion in the case of activator-depleted reaction kinetics for the cases of exponential, linear and logistic growth profiles. Domain-growth induces parameter spaces substantially different from those obtained in the absence of domain-growth. By selecting model parameter values a and b only from the domain-growth induced parameter spaces, finite element numerical simulations exhibit the formation of spatial structure during growth development. Such structures do not form when domain growth is absent. These spatial structures are an emergent process due to domain evolution. The emergence of spatial structure during growth development occurs with or without cross-diffusion. Cross-diffusion enhances pattern formation through the emergence of parameter spaces outside of the standard reaction-diffusion without cross-diffusion by allowing for more general choices of the diffusion coefficient d which would not otherwise give rise to patterning in the absence of linear crossdiffusion (e.g. by taking d = 1, say). We have also exhibited patterning on evolving surfaces to reinforce the fact that studies on growing domains generalises to other complex evolving manifolds as well [33] . We have not presented numerical results for linear and logistic growth profiles since these are essentially identical to those already presented in [33,35]. Unlike in previous studies, for exponential growth, we have solved the non-autonomous system of ordinary differential equations (3.3)-(3.4) to obtain numerically the timedependent manifold (u S (t), w S (t)) which is used to compute the domain-induced parameter spaces. In previous studies, a steady state solution was assumed.
Our studies naturally extend to consider nonlinear cross-diffusion where weakly nonlinear analysis can be used to study the dynamical behaviour of the reactiondiffusion systems in the presence of nonlinear cross-diffusion [16,17]. It will be interesting to investigate how nonlinear cross-diffusion coupled with domain growth influences pattern formation during growth development. Mathematical analysis of reaction-diffusions systems in general (and in particular with nonlinear crossdiffusion) involving non-uniform growth is completely open and remains largely unanswered.
Data statement. All data are provided in full in the results section of this paper. The numerical algorithms are available on request.