STABILITY ANALYSIS OF REACTION-DIFFUSION MODELS ON EVOLVING DOMAINS: THE EFFECTS OF CROSS-DIFFUSION

This article presents stability analytical results of a two component reaction-diffusion system with linear cross-diffusion posed on continuously evolving domains. First the model system is mapped from a continuously evolving domain to a reference stationary frame resulting in a system of partial differential equations with time-dependent coefficients. Second, by employing appropriately asymptotic theory, we derive and prove cross-diffusion-driven instability conditions for the model system for the case of slow, isotropic domain growth. Our analytical results reveal that unlike the restrictive diffusion-driven instability conditions on stationary domains, in the presence of cross-diffusion coupled with domain evolution, it is no longer necessary to enforce cross nor pure kinetic conditions. The restriction to activator-inhibitor kinetics to induce pattern formation on a growing biological system is no longer a requirement. Reaction-cross-diffusion models with equal diffusion coefficients in the principal components as well as those of the short-range inhibition, long-range activation and activator-activator form can generate patterns only in the presence of cross-diffusion coupled with domain evolution. To confirm our theoretical findings, detailed parameter spaces are exhibited for the special cases of isotropic exponential, linear and logistic growth profiles. In support of our theoretical predictions, we present evolving or moving finite element solutions exhibiting patterns generated by a short-range inhibition, long-range activation reactiondiffusion model with linear cross-diffusion in the presence of domain evolution. 2010 Mathematics Subject Classification. Primary: 35K55, 35K57, 35K58; Secondary: 37B55, 37C60, 37C75.


1.
Introduction. Understanding of biological processes during growth development is an unresolved issue in developmental biology that is only starting to be addressed in the last decade. Introducing domain growth into the modelling results in non-autonomous systems of partial differential equations whose analytical tractability is not yet well understood [7,11,18,19,24]. In the area of developmental biology, partial differential equations for pattern formation take the form of reaction-diffusion type [13,30,37]. 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 [3]. 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 reactiondiffusion system, a key requirement for diffusion-driven instability is the concept of long-range inhibition and short-range activation [10]. This implies that one of the species (the inhibitor) must diffuse faster (typically much faster) than the autocatalytic species (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 [7,18,19,24]. It is well-known that many problems in biology involve growth. In [24] 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/activation chemical process 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.
Another generalisation is the introduction of cross-diffusion. In many multicomponent systems, there are various forms of diffusion depending on the biochemical problem at hand [38]. Diffusive processes can be characterised as selfdiffusion, cross-diffusion, mutual-diffusion, tracer-diffusion, intra-diffusion, interdiffusion, uphill-diffusion and negative-or incongruent-diffusion. A detailed review of the different physicochemical interpretations of these forms of diffusion is given by Vanag and Epstein [38]. 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 [28,40]. Multicomponent systems containing nanoparticles, surfactants, polymers and other macromolecules in solution play an important role in industrial applications and biological functions [28]. In developmental biology, recent experimental findings demonstrate that cross-diffusion can be quite significant in generating spatial structure [38]. The effects of cross-diffusion on models for pattern formation (i.e. reaction-diffusion type) have been studied in many theoretical papers [5,6,8,9,12,14,33,34,36,41,42,43]. Recently, in Madzvamuse [26] we showed that introducing cross-diffusion to a system of reactiondiffusion 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 cross-diffusion 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 [7,11,15,16,17,19,24,27]. In this article our focus is to study, both theoretically and computationally, how cross-diffusion induces patterning in the presence of domain evolution.
Hence, our paper is organised as follows: in Section 2 we present the model equations posed on evolving domains and these consists of a system of reaction-diffusion equations with linear cross-diffusion. Domain evolution terms are modelled through dilution and convective terms. For analytical purposes, in this section, we map at all times, the model system from a continuously evolving domain to a reference frame (which can be taken as the initial domain). As a result, a system of non-autonomous partial differential equations is obtained. The bulk of our work is detailed in Section 3 where we derive and prove the conditions for cross-diffusion driven instability on evolving domains. These conditions are a generalisation of the classical Turing and cross-diffusion driven instability conditions on stationary domains. Cross-diffusion induced parameter spaces are computed and exhibited in Section 4 for the special cases of isotropic exponential, linear and logistic growth profiles. It is in this section that we exhibit several parameter spaces induced by cross-diffusion in the presence of domain evolution. By selecting model parameter values from the cross-diffusion induced parameter spaces, in Section 5, we present finite element solutions on twodimensional evolving domains demonstrating the emergence of patterns induced by cross-diffusion and domain growth, for the case of the short-range inhibition, longrange activation. Such patterns can not be formed in the absence of cross-diffusion. In Section 6, we conclude and discuss the implications of our findings to the theory of pattern formation during planar growth development.
2. Reaction-diffusion systems 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) , v (x, t)) T be a vector of two chemical concentrations at position x ∈ Ω(t) ⊂ R m . The evolution equations for reaction-diffusion 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 nondimensional system [11,24,30] where ∇ 2 is the Laplace operator on domains and volumes, d = Dv Du is the ratio of the diffusion coefficients and d u = Duv Du and d v = Dvu Du are the ratios of the crossdiffusion, where D u , D v , D uv and D vu are dimensional diffusion and cross-diffusion coefficients, respectively. Here, n is the unit outward normal to Ω(t). In the above, V (x, t) denotes the flow velocity of the chemical species due to domain evolution. Initial conditions are prescribed through non-negative bounded functions u 0 (x) and v 0 (x). The functions, f (u, v) and g(u, v) represent nonlinear reactions. Note that we have imposed self-organisation on the system through zero-flux (also known as homogeneous Neumann) boundary conditions. Remark 1. Domain growth has the effect of introducing extra terms to the classical model of reaction-diffusion with linear cross-diffusion. For example, for the u equation these are of the form where V · ∇u represents the transport of the chemical species by the flow velocity V and u(∇ · V ) denotes the dilution (or concentration) of the chemical species due to domain growth (or contraction) if ∇ · V > 0 (or ∇ · V < 0), respectively.
For analytical clarity and simplicity, we will write the model system (1) in compact-vector form as where Assumption 2.1 (Flow velocity). We assume that the flow velocity V (x, t) is identical to the domain velocity, i.e., V = dx dt as is standard in the derivation of classical reaction-diffusion models on evolving domains on application of the Reynold's Transport Theorem [1]. Lemma 2.1 (Global existence of solutions of reaction-diffusion systems with cross-diffusion on evolving domains). Under suitable assumptions on the reaction-kinetics f (u, v) and g(u, v), the diffusion matrix D and the domain evolution, problem (1) admits a global classical solution.
Proof. The proof follows directly from Theorem 4.5 in Venkataraman et al., [39] (Theorem 4.5, page 50) under suitable assumptions on the reaction-kinetics f (u, v) and g(u, v), the diffusion matrix D and the domain evolution.

2.2.
Mapping to a stationary domain using a Lagrangian transformation. For analytical convenience we will work with the model system posed on a timeindependent domain. Therefore, we introduce a transformation that maps the model system (2) from a continuously evolving planar domain to a time-independent (stationary) domain (see [2,11,24] for more details of this approach). In order to do this, without too many technical complications, we will restrict our attention to a special class of domain evolution as detailed in the next section. Assumption 2.3 (Isotropic domain evolution). We assume a spatially linear isotropic evolution of the domain Ω(t) of the form where ϕ(t) is a continuously differentiable function with ϕ(0) = 1. In the above, ξ represents the spatial coordinates of the stationary domain. This assumption and Assumption 2.1 imply that whereφ := dϕ dt . Proposition 1. Assumptions 2.1 and 2.3 imply that the divergence of the flow velocity is given by where m defines the number of spatial dimensions.
We are now in a position to state the reaction-diffusion system with linear crossdiffusion posed on a stationary domain. Assuming uniform isotropic evolution of the domain, we can writê As a result, model (2) simplifies to where h(t) := ∇ · V . From here onwards, we will drop the hatsˆfor ease of exposition.
2.3.1. Timescales. For a biological system, growth is typically driven by cell division, and thus occurs on the timescale of the cell cycle duration, which is typically 24 hours or greater. In contrast, the biochemical reaction kinetics are typically considered to occur on a much faster time scale of seconds to minutes, and are typically constrained by the diffusive dynamics. The length scale associated with a volume of 10 4 to 10 6 cells is about 2 × 10 −4 m to 10 −3 m. Given diffusion and cross-diffusion coefficients of between 10 −6 cm 2 s −1 to 10 −5 cm 2 s −1 , the diffusive and cross-diffusive timescales occur between 40 seconds to 170 minutes [38]. We will take advantage of this difference in timescales in our analysis. With T g denoting the growth timescale, and T dyn denoting the maximum of the diffusive timescale and the biochemical kinetics timescale, we will define and utilise the small parameter and therefore we will generally neglect O( 2 ) corrections. Following [24], it can be shown that for the above range of estimates ∈ [5 × 10 −4 , 0.12].
3. Analysis of domain-induced cross-diffusion driven instability.

3.1.
Definitions. In the following, we carry out cross-diffusion driven instability analysis of a spatially linear, isotropic evolution of the domain. Given that model system (5) is a non-autonomous system of partial differential equations with timedependent 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) defined as follows.
Definition 3.1 (A spatially independent time-dependent manifold). A spatially independent manifold u S (t) is a solution of (5) if it solves the non-autonomous system of ordinary differential equations Remark 2. In all our numerical examples presented in Sections 4 and 5, initial conditions u S (0) are taken as small random perturbations around the uniform steady obtained on stationary domains.
Definition 3.2 (Cross-diffusion driven instability on evolving domains). Crossdiffusion driven instability on evolving domains occurs when a spatially independent manifold u S (t), linearly stable in the absence of spatial variations (diffusion and cross-diffusion terms) becomes unstable in the presence of spatial variations.
In the above, J F (t) is the Jacobian matrix corresponding to the linearisation of the non-linear reaction kinetics which is evaluated at u S (t). We define ψ K (ξ) to be the time-independent eigenmode of the transport operator satisfying      A standard analysis shows that the eigenvalue must be real and negative, and thus without loss of generality, we write it as −K 2 in the above. We similarly have that eigenmodes of distinct eigenvalues are orthogonal. Furthermore, we note that the domain Ω(0) is time-independent because ξ is a Lagrangian coordinate system. We proceed by expanding w in terms of the eigenmodes Linearity and the orthogonality of the eigenmodes allow each mode to be considered separately, which is a substantial simplification. In particular, growth of any one of the β K (t) with time is sufficient to drive the full solution away from the time-dependent solution u S (t). Similarly, if all the β K (t) decay with time then a sufficiently small perturbation from the time-dependent solution u S (t) will decay, at least within the resolution of linear theory predictions.
To proceed, we substitute, for each K, a mode w K = β K (t)ψ K (ξ) from expansion (10) into the linearised equation (8) to find Writing where q(t) = exp − t t0 h(τ )dτ and t 0 is the point the perturbation is applied, we have One can immediately deduce that Given initial conditions, b K (t 0 ) at time t 0 we seek a solution at time t > t 0 of the form where λ * K (t) is either the largest time-dependent real eigenvalue of M K (t) or the real part of one of the time-dependent complex conjugate eigenvalues. From equations (13) and (14) we have with c K (t 0 ) = b K (t 0 ).

3.3.
Analysis of solutions of non-autonomous systems of ordinary differential equations. 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., [24]. The only difference now is the inclusion of the cross-diffusion coefficients into the matrix M K (t). We present a summary of the analysis for completeness' sake, the reader is referred to [24] for further details. First, we observe that for a matrix Q, the induced matrix norm is defined as ||Q|| = sup ||x||=1 ||Qx|| and the Lozinskii measure is [29] µ(Q) = lim It can be easily shown that ||Qx|| ≤ ||Q|| ||x||, for any matrix Q and vector x.
Proposition 2. Assume that M K (t) has real and distinct eigenvalues, then with λ * K the largest real eigenvalue of M K .
Proof. First, we prove that c K (t) is bounded from above. The matrix Q K (t) possesses one zero and one negative eigenvalue. Let us denote by e 0 K (t) the unit eigenvector associated with the zero eigenvalue and let e 1 K (t) be the unit eigenvector associated with the negative eigenvalue, such that the angle between e 0 K (t) and e 1 K (t) is acute. These eigenvectors are distinct and are linearly independent and thus form a spanning set. Hence we can write The above satisfies all the axioms required of a norm. It then follows that noting that a matrix with zero and negative eigenvalues has a zero Lozinskii measure. Thus a solution to equation (15) does not grow in time.
Next, we prove that c K (t) does not decay to zero. Taking t fixed, with ∆t = t − t 0 ∼ T dyn , we can find at least one solution of (15) which does not decay on neglecting O( 2 ) corrections. From a Picard iteration we have where The existence of G K (t, t 0 ) is guaranteed from Picard's theorem given the components of Q K (t) are bounded to ensure that Q K (t)c K (t) is Lipschitz. With "dot" denoting a time derivative, we havė since the change in the matrix Q is driven by growth. We set the initial condition c K (t 0 ) = e 0 K (t 0 ) and consider the solution (17) obtained by expanding each Q K (t) in expression (18). Noting Q K (t 0 )e 0 K (t 0 ) is zero and that the n th order time derivatives of Q K (t) scale with n , we have By differentiating Q K (s) e 0 K (s) = 0 with respect to s and then setting s = t 0 we haveQ Let e 0P K (t) be the unit vector which is perpendicular to e 0 K (t) where P denotes perpendicular. Because e 0 K (t) is defined to be a unit vector, we havė where τ (t) has an implicit K dependence. Since the matrix Q K (t) changes on the growth timescale we also have that /τ (t) ∼ O(T −1 growth ) and hence τ (t) ∼ O(T dyn ). We also have the kinematic relation by projecting e 1 K (t) onto the e 0 K (t) and e 0P K (t) directions. Note that ψ(t) has an implicit K dependence and that, without loss of generality, ψ(t) ∈ [0, π/2). Hence . By substituting this into equation (19) we have which is not decaying. From equations (14), (21) and the observation that with λ * K the largest real eigenvalue of M K .
Proposition 3. Assume that M K (t) has complex conjugate eigenvalues, then where are real, distinct, vectors and thus a spanning set of the real plane. Without loss of generality we take f 0 K (t) to be a unit vector and let Ω K def = Im(λ(t)) > 0. We can write any real vector function of time in the form . Following [24] letē K (t) denotes the conjugate of e K (t) then we can writė where A K , B K ,Â K ,B K are defined by the above relations and are time-dependent. The appears as the change in the eigenvector e K (t) is on the timescale of growth. As long as the angle between f 0 K and f 1 K is significantly greater than and the vectors have a ratio of lengths, l, with l −1 , then A K (t) and B K (t) can be treated as O( 0 ) terms for the purposes of asymptotic expansions. This is assumed below.
We also have that and hence, by differentiation, as Ω K changes on the growth timescale T growth = T dyn / . Thus for a general initial condition Calculating the norm reveals K is a phase angle. We can rewrite the above norm in the form where Thus we have an oscillating solution, the amplitude of which is constant at leading order but not at higher orders. One can immediately deduce that, accurate to O( 2 ) 2144 ANOTIDA MADZVAMUSE, HUSSAINI S. NDAKWO AND RAQUEL BARREIRA Lemma 3.3. The important result for investigating stability is that an O( 2 ) accurate estimate of the growth rate between time t 0 and time t is given by with λ * K denoting the real part of the eigenvalues of M K .
Remark 3. The case when the matrix M K (t) has real, repeated eigenvalues requires a mathematically precise parameter fine tuning, and thus is unlikely to occur in most models which is inappropriate for biological models.
3.4. Cross-diffusion driven instability conditions on evolving domains. Following Madzvamuse et al., [24], 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 cross-diffusion driven instability. For this, it will be useful to define t 1 : Theorem 3.4. The necessary conditions for a cross-diffusion driven instability (accurate to first order in 1) 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 (5), are given by In the above, the subscripts u, v denote partial differentiation, with the Jacobian components f u , f v , g u and g v evaluated in terms of u S (t 1 ). Strictly, we have h * ∼ 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. Let λ * K (t 1 ) be the largest real part of any eigenvalue of M K (t 1 ), λ k (t 1 ) which is the root of the quadratic equation Thus, λ K (t 1 ) satisfies the dispersion relation where the scalar variables and kinetic functions. The partial derivatives are evaluated in terms of u s (t 1 ) the solution to the non-autonomous system of ordinary differential equations given by (7) at time t 1 .
Let us define the eigenvalue where b(K 2 * ) and c(K 2 * ) are given in equations (32) and (33), respectively. For convenience's sake, denoting by b h * (K 2 , we obtain the characteristic equation Thus Taking K * = 0 we have the absence spatial variations and thus spatial homogeneity. In the absence of spatial variations, we require u s (t 1 ) to be asymptotically stable to the K * = 0 and this is guaranteed provided This is guaranteed provided b h * (0) > 0 and c h * (0) > 0 if and only if the following conditions hold These first two conditions are domain-induced, and do not reflect the spatial variations as expected. The conditions enforce the requirement that the time-dependent manifold u s (t) is asymptotically stable with respect to spatially homogeneous perturbations. Now, in the presence of spatial variations which includes both diffusion and crossdiffusion, (K 2 since b h * (0) > 0. For (u s (t 1 )) to become unstable, we require that 2146 ANOTIDA MADZVAMUSE, HUSSAINI S. NDAKWO AND RAQUEL BARREIRA thereby requiring that c h * (K 2 * ) < 0 for some K 2 * non-zero. By definition of c h * (K 2 * ) we can further re-arrange to obtain a quadratic polynomial in K 2 * of the form where . Geometrically, c h * (K 2 * ) represents a parabola. For the reaction-diffusion system with cross-diffusion to be well-posed, we require that which implies that the parabola opens upwards. This results in the third condition for cross-diffusion driven instability. It can be easily shown that c h * (K 2 * ) has a minimum value which is negative provided where In all the above, ignoring O(h 2 * ) results in the five conditions for cross-diffusion driven instability conditions (24)- (28).
Remark 4. In addition to inequalities (24)-(28), a cross-diffusion driven instability requires that there exists at least one wavenumber such that K 2 * is contained in the interval where K 2 ± are the roots of c h * (K 2 ± ) = 0 and these are given by 1. In contrast to the autonomous case, if we have an activator and an inhibitor, we do not require short-range activation, long-range inhibition which is the standard mechanism for a diffusively-driven instability. For example, suppose without loss of generality that u is the activator in equation (1) and v is the inhibitor. Thus f u > 0 and g v < 0. From inequalities (24) and (27) we require where γ > 0. Thus, the inhibitor has to diffuse faster than the activator in the absence of growth and cross-diffusion for an instability to be possible. (b) In the absence of growth, i.e. h * = 0, but with cross-diffusion in both components we have Hence, we can take 0 < d ≤ 1 which implies that (i) if f v > 0 and g u > 0, then at least one of the cross-diffusion coefficients must be negative and both must satisfy if f v < 0 and g u < 0, then either at least one of the cross-diffusion coefficients must be negative or both are positive, and for both cases, The presence of growth, i.e. h * = 0 does not alter the fact that 0 < d ≤ 1 for the reaction-diffusion system with linear cross-diffusion in both components, and therefore, the above characterisation of the cross-diffusion and reaction-kinetics still holds. It follows therefore that non-standard reaction-kinetics such as those of the form short-range inhibition, long-range activation, short-range inhibition, longrange inhibition, short-range activation, long-range activation, same-range inhibitor-inhibitor (i.e. d = 1), same-range activator-activator (i.e. d = 1), and same-range activator-inhibitor (i.e. d = 1) can also generate patterns in the presence of cross-diffusion in both components on stationary as well as on evolving domains. 2. In contrast to the autonomous case an activator and inhibitor interaction is not required to generate an instability; for example, two activators can satisfy the instability constraints. Again suppose that u is the activator in equation (1) so that f u > 0. Using inequality (24) we find where γ > 0. Thus, for h * = 0, as with no growth, we have g v < 0, implying that v must be an inhibitor. In contrast for h * > 0 we clearly are not required to have g v < 0 and hence the biochemical corresponding to the concentration v can also be an activator without making a diffusively-driven instability impossible. Therefore an activator-activator (on growing domains) and inhibitor-inhibitor (on contracting domains) mechanisms could give rise to Turing patterns only in the presence of domain growth.

4.
Exhibiting cross-diffusion induced parameter spaces on evolving domains. 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 cross-diffusion. We will show that 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.

Type
Growth 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. Let x(t) = ξϕ(t) where ϕ(t) > 0 is the domain growth function satisfying ϕ(0) = 1. We can compute the domain (mesh) velocity as and hence 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. , m ≥ 1, is a monotonically decreasing function whose limit tends to zero, i.e.
Proof. The proof follows clearly from the definitions of h(t) shown in Table 1.
where a and b are positive parameters to be selected from domain-and crossdiffusion-induced parameter spaces.   For a fixed set of diffusion and cross-diffusion coefficients, varying the exponential growth rate r results in larger and larger parameter spaces emerging during the exponential evolution of the domain. On the other hand, fixing the exponential growth rate r, parameter spaces of different sizes are obtained with the largest parameter space corresponding to the reaction-diffusion system with cross-diffusion in the u-equation only, while the smallest parameter space corresponds to the reactiondiffusion system with cross-diffusion in the v-equation only. This is consistent with theoretical results presented in [26] for the case of stationary domains.
Cross-diffusion driven instability on evolving domains for d = 1 and d u = 1. Here we fix d = 1 and d u = 1 and vary positively only the cross-diffusion coefficient d v such that the inequality d − d u d v > 0 is satisfied. Figure 2 shows the cross-diffusion parameter spaces with cross-diffusion coefficient d v = 0.5 (light blue), d v = 0.6 (red), d v = 0.7 (green) and d v = 0.8 (dark blue), for different exponential growth rates.    As d v −→ 1, we observe larger and larger parameter spaces, while when d v −→ 0, smaller and smaller parameter spaces are exhibited and these vanish as d approaches zero (see Figure 3). In Figure 3 we superimpose parameter spaces when (i) d v is varied positively and (ii) the exponential growth rate r is also varied. We observe that for lower values of d v , distinct parameter spaces are obtained, while for larger values of d v approaching one, the parameter spaces overlap substantially.
Cross-diffusion driven instability on evolving domains for d = 1 and d v = 1. Here we fix d = 1 and d v = 1 and vary positively only the cross-diffusion coefficient d u . Figure 4 shows the cross-diffusion parameter spaces with cross-diffusion coefficient d u = 0.5 (light blue), d u = 0.6 (red), d u = 0.7 (green) and d u = 0.8 (dark blue), for different exponential growth rates.   Similarly, variations in d u results in smaller and smaller parameter spaces as d u −→ 0 and larger and larger parameter spaces as d u −→ 1. Individual parameter spaces are superimposed as illustrated in Figure 5. As the exponential growth rate r increases, for a fixed set of diffusion and cross-diffusion coefficients, substantially different parameter spaces emerge. Exponential evolution of the domain results in substantially topologically different parameter spaces from those obtained in the absence of growth.
Negative cross-diffusion driven instability on evolving domains for d = 1 and d v = 1.
Unlike previous examples, we fix d = 1 and d v = 1 and vary negatively only the cross-diffusion coefficient d u such that the inequality d − d u d v > 0 is satisfied. Figure 6 shows the cross-diffusion parameter spaces with cross-diffusion coefficient d u = −0.8 (light blue), d u = −0.7 (red), d u = −0.6 (green) and d u = −0.5 (dark blue), for different exponential growth rates.  We observe that substantially and topologically different parameter spaces are obtained as the exponential growth rate is varied. However, for a fixed exponential growth rate, variations in the cross-diffusion coefficient do not alter substantially the parameter spaces. This is unlike previous results where variations in crossdiffusion coefficients could alter substantially the parameter spaces (see Figures  1-2

, for example).
Short-range inhibition, long-range activation: Cross-diffusion induced parameter spaces. For the first time, we present a short-range inhibition, long-range activation model that can only give rise to pattern formation in the presence of crossdiffusion. From the cross-diffusion driven instability condition d−d u d v > 0, we take 0 < d = 0.5 < 1, i.e. 0 < D v < D u in dimensional units, and therefore the activator u diffuses much faster than the inhibitor v. To proceed, we further fix d v = 0.8 and vary negatively d u such that the condition d − d u d v > 0 is satisfied. Figure 7 shows the parameter spaces demonstrating the ability of the model system to generate patterns in the presence of cross-diffusion. Without cross-diffusion, such spaces do not exist. We note that d u is constrained by the fact that −1 < d u < 0 in order for the u-equation in system (1) to be well posed. On the other, d v is constrained by the cross-diffusion driven instability condition (26). As the cross-diffusion coefficients d u and d v increase and decrease towards zero values, cross-diffusion driven instability spaces reduce in size and disappear altogether (results not shown). This is consistent with theoretical result that for diffusion-driven instability to occur for

4.2.2.
Linear and logistic domain evolution. Next we compute parameter spaces generated by linear and logistic growth functions for two-dimensional evolving domains. For the linear and logistic growth functions we compute the time-dependent solutions (u S (t), v S (t)) which satisfy the differential equations where h(t) = ∇ · V and given initial values at time t 0 = 0. From 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 . Let us take r = 0.01( 1), m = 2.0, κ = 2 and A = 1. In all our numerical examples, initial conditions (u S (0), v S (0)) are taken as small random perturbations around the uniform steady obtained on stationary domains. 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 (24)-(28) for the cases of linear and logistic growth functions, the non-autonomous differential equations (48)-(49) must be solved for every point (a, b) in the plane at any given time t thereby giving rise to time-dependent parameter spaces. Figure 9-10 exhibit time-dependent parameter spaces generated as a result of linear (first column) and logistic (second column) growth functions of the domain. The first column of Figure 9 corresponds to the linear growth; (c), (e) and (g) are the individual spaces of (a) and the second column corresponds to the logistic growth; (d), (f) and (h) are individual spaces of (b).
Unlike parameter spaces generated by the exponential growth, linear and logistic growth functions generate parameter spaces which overlap substantially. Linear and logistic growth functions generate larger but topologically similar parameter spaces to those obtained in the absence of domain growth. This is in contrast to the exponential growth which generates larger and sometimes topologically different parameter spaces.
In Figure 10 we present parameter spaces corresponding to the linear and logistic growth functions when we fix d = 0.5, d u = −0.5 and vary d v = 0.6, 0.7, 0.8 respectively. As d v −→ 1, larger and larger parameter spaces are observed, while d v −→ 0, smaller and smaller parameter spaces emerge.
Remark 6. It must be noted that all the parameter spaces exhibited in Figures  2-10 are cross-diffusion and domain-growth induced spaces, such parameter spaces do not exist for the classical autonomous system of reaction-diffusion equations.

Domain-and cross-diffusion induced patterns on evolving domains.
In this section we exhibit domain-and cross-diffusion induced patterns on planar evolving square domains. Surface-and cross-diffusion induced patterns on evolving surfaces are presented in Madzvamuse and Barreira [25]. Unlike any previous studies of these types of models, our focus is to present patterns generated from a nonstandard reaction-diffusion system of the form long-range activation, short-range inhibition model that will only give rise to patterning in the presence of crossdiffusion. The effects of domain growth enhances patterning independent of the initial conditions. The patterns exhibited in this section can only be obtained in the presence of cross-diffusion, such a model of the form long-range activation, shortrange inhibition does not give rise to patterning in the absence of cross-diffusion, with or without domain evolution. For illustrative purposes, we take an activatordepleted model known to satisfy the cross/pure kinetics conditions [10,32,35]. Our focus is to study the effects of cross-diffusion on a standard reaction-diffusion model in the presence of domain evolution. We leave studies of non-standard reaction kinetics capable of giving rise to patterning only in the presence of cross-diffusion and domain evolution for future studies. To the best of the authors' knowledge, only models of the form long-range inhibition; short-range activation have the ability to give rise to patterning in the absence of cross-diffusion. We will depart from this framework and select diffusion and cross-diffusion parameter values such that a short-range inhibition; long-range activation can exhibit cross-diffusion induced patterns on evolving domains.l In Madzvamuse and Barreira [25] we presented a wide range of patterns on evolving domains and surfaces. In most of these examples, parameter values were selected outside the standard Turing parameter space (see [25] for details). Comparisons were then made between patterns obtained between the autonomous standard reaction-diffusion model and the non-autonomous nonstandard reaction-diffusion model.
We select two sets of model kinetic parameter values from parameter spaces induced by exponential, linear and logistic growth functions as illustrated in Figures  7 (a) and 10, respectively. These are (i) a = 0.15, b = 0.2 and (ii) a = 0.15 and b = 0.15 respectively. The diffusion and cross-diffusion coefficients are fixed as d v = 0.5, d u = −0.5. and d u = 0.6 (corresponding to case (i)) and d u = 0.8 (corresponding to case (ii)) above respectively. For all the growth functions, we take the growth rate r = 0.01. For illustrative purposes, we take γ = 200. In all the cases, random initial perturbations around the uniform steady state which is calculated on stationary domains are prescribed.
Numerical simulations are carried out using the evolving or moving finite element method; details of the numerical method are given in Madzvamuse et al., [25]. The moving grid finite element method has been applied extensively to study reaction-diffusion systems on stationary and evolving domains with applications to developmental biology. The reader is referred to [4,15,18,19,20,21,25,26] for further references detailing the numerical method and its applications. In Figures 11 -22 we exhibit the results of our simulations, validating theoretical results presented in the previous section. Furthermore, in order to understand the evolution of the bifurcation process during growth development, we plot the evolution of the log of the L 2 -norm of the errors between successive numerical iterate solutions corresponding to patterns shown in Figures 12,14,16,18,20 and 22. Figure 11 shows pattern formation during growth development for the parameter values d = 0.5, a = 0.15, b = 0.2, d u = −0.5, d v = 0.6 and γ = 200. We observe the formation of spots which correlate to the growth of the error norm at the initial stages of growth (see first peak observed in Figure 12); these are rapidly destabilised; with stripe patterns forming (which also now correlate to the second sharp pick observed in Figure 12). The stripe patterns are stable to domain growth for a substantial amount of time (this is manifested by the reduction in the errors which become almost constant as illustrated in Figure 12). Destabilisation of the stripes results in the rapid formation of spot patterns (of different wavelenghts) but these are less stable to domain growth (compare Figure 11 (f) and (g)), thereby giving way to the formation of stripe patterns again. This is reflected in the error graph where there is a rapid growth and decay in errors, but once the spots are formed, there is a rapid, almost period growth and decay in errors. This evolutionary process of spots forming and being destabilised rapidly and then the formation of stripes or spots continues almost periodically. We observe the formation of "peanut" patterns, reminiscent of the activator-depleted model. The evolutionary process becomes substantially different if a different growth function is employed. For example, a linear evolution of the domain results in the formation of spot patterns, similar to the exponential growth at the very early stages of growth development, and these transient into stripes, which are stable for a substantially longer period of time than for the case of the exponential growth (compare Figures 12 and 14). In fact, for the linear growth, the only patterns to form later during growth development are spots which are only stable for a vey short time; these give way to the formation of stripe patterns which continue to be stable (for as long as we evolved the domain).  Figure 11, for the case of exponential growth.
On the other hand, employing a logistic growth function, we observe a remarkably different bifurcation process during domain evolution. The error graph shows only three substantial bifurcation processes (compared to six for the exponential and four or five for the linear growth functions). The logistic growth function is characterised by rapid growth at the early stages (which might imply that lower modes are almost unexcitable during growth development), followed by a linear growth at mid-stages of the domain evolution and then saturation of the growth evolution at later stages. The three phases of logistic growth are therefore characterised by stripes bifurcating into spot-stripe patterns which subsequently transient into circular patterns (assuming one is following the red patterns), and these finally bifurcate into stripes which remain stable as the domain approaches the final size (see Figure 15  To the authors' best knowledge, these patterns are the first to be exhibited and are solely due to cross-diffusion. We also observe, as before, that linear and logistic growth functions seem to stabilise more robustly stripe patterns at later stages of the domain evolution than spot patterns. The exponential evolution favours the continuous formation of stripe, spot, stripe-spot and circular patterns as the domain continues to evolve.  Figure 13, for the case of linear growth. from model system (1) and that without cross-diffusion nor domain growth, the interested reader is referred to Madzvamuse and Barreira [25]. 6. Conclusion and discussion. In this article we have generalised substantially the theory for reaction-diffusion systems with cross-diffusion to consider the effects of domain evolution on cross-diffusion induced patterning. Our results reveal that in the presence of cross-diffusion and domain evolution, the restrictive conditions associated with long range inhibition, short range activation can be relaxed to allow a wide range of model parameter values and reaction-kinetics that would otherwise not give rise to patterning either in the absence of cross-diffusion and/or domain evolution. We have shown that a long range activation, short range inhibition can generate patterns only in the presence of domain evolution and/or cross-diffusion.
These theoretical findings open new research directions, both experimentally and theoretically, where non-standard mechanisms for patterning can be considered and designed.
Theoretically, we have derived and proved the conditions for cross-diffusion driven instability. By using these conditions, we then computed and exhibited a wide variety of parameter spaces for exponential, linear and logistic growth functions. Our results unravel the emergence of parameter spaces with interesting characteristics. For a fixed set of appropriate parameter values, exponential evolution of the domain results in substantially different spaces from those obtained in the absence of domain evolution. Although linear and logistic growth functions yield larger parameter spaces in the presence of domain evolution, these are topologically similar to those obtained in the absence of domain evolution [19,24,25].
Our analytical results reveal that unlike the restrictive diffusion-driven instability conditions on stationary domains, in the presence of cross-diffusion coupled with domain evolution, it is no longer necessary to enforce cross nor pure kinetic conditions. The restriction to activator-inhibitor kinetics to induce pattern formation on a growing biological system is no longer a requirement. Reaction-cross-diffusion models with equal diffusion coefficients in the principal components as well as those of the short-range inhibition, long-range activation and activator-activator form can generate patterns only in the presence of cross-diffusion coupled with domain evolution.  Figure 15, for the case of logistic growth.
In summary, the following characteristics are exhibited by reaction-diffusion systems with linear cross-diffusion in the presence of an exponential domain growth: 1. Large parameter spaces emerge during an exponential evolution of the domain.
In all our simulations, the reaction-diffusion system with cross-diffusion in the u-equation only possesses the largest parameter space in the presence of domain evolution. This is an inherent property of the model and holds true in the absence of domain evolution. 2. On the other hand, the reaction-diffusion system with cross-diffusion in the vequation only possesses the smallest parameter space during domain evolution, in complete agreement with the model system on stationary domains. Again, this is an inherent property of the model system. 3. In all our simulations, the parameter spaces corresponding to the reactiondiffusion system with cross-diffusion in the v-equation only are sub-spaces (for slow domain evolution) of the parameter spaces corresponding to the reactiondiffusion system without cross-diffusion. These, in turn, are sub-spaces of the reaction-diffusion system with cross-diffusion in both the u and v-equations.
Similarly, there parameter spaces are contained fully in the largest parameter spaces corresponding to the reaction-diffusion system with cross-diffusion in the u-equation only. However, for the case of fast domain evolution, distinct and substantially different parameter spaces are obtained. 4. For a given set of diffusion and cross-diffusion values, exponential evolution of the domain results in larger and distinct parameter spaces with or without intersections. The fact that in general there is no overlap between the zero-growth rate parameter spaces and those on evolving domains leads us to conclude that, in general, cross-diffusion driven instability in the absence of growth need not imply diffusively-driven instability in the presence of growth.
To demonstrate the effects of domain evolution on cross-diffusion induced patterning, numerical experiments computed using the evolving or moving finite element method [19,25] are presented for specific choices of parameter values selected from the cross-diffusion induced parameter spaces. An important key observation is that unlike on stationary domains, pattern evolution during growth development is independent of the initial conditions. On stationary domains, patterns are known to depend crucially on the initial conditions which act as a basin of attraction.   Figure 17, for the case of exponential growth.
We have previously reported in the literature that domain growth enhances robustness of certain patterns in a Turing mechanism [27,31]. In this paper, we have shown that growth does not only induce instability in a reaction-diffusion system with cross-diffusion but that cross-diffusion coupled with domain growth expands substantially the range of mechanisms that can give rise to spatial patterns away from the classical short-range activation, long-range inhibition paradigm. Understanding the effects of cross-diffusion to the theory of pattern formation during growth development is crucial in many areas of research such as nanoparticles, surfactants and polymers.
From a theoretical point of view, our studies raise two key research questions to be addressed; namely (i) the derivation of non-standard reaction kinetics that will give rise to cross-diffusion driven-instability only in the presence of cross-diffusion and/or domain evolution and (ii) the derivation of an analytical framework to study bifurcation sequences and processes for non-autonomous reaction-diffusion systems with linear and nonlinear cross-diffusion. It is clear from the L 2 norm graphs of the errors that there exists smooth as well as sharp bifurcation transitions as the transient solutions gain or loose stability during domain evolution. Some of these transitions are reflective of the bifurcation process known as hysteresis [30]; only detailed analytical studies of these models on evolving domains will answer and confirm such observations.