Numerical Preservation of Velocity Induced Invariant Regions for Reaction–Diffusion Systems on Evolving Surfaces

We propose and analyse a finite element method with mass lumping (LESFEM) for the numerical approximation of reaction–diffusion systems (RDSs) on surfaces in R3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbb {R}}^3$$\end{document} that evolve under a given velocity field. A fully-discrete method based on the implicit–explicit (IMEX) Euler time-discretisation is formulated and dilation rates which act as indicators of the surface evolution are introduced. Under the assumption that the mesh preserves the Delaunay regularity under evolution, we prove a sufficient condition, that depends on the dilation rates, for the existence of invariant regions (i) at the spatially discrete level with no restriction on the mesh size and (ii) at the fully-discrete level under a timestep restriction that depends on the kinetics, only. In the specific case of the linear heat equation, we prove a semi- and a fully-discrete maximum principle. For the well-known activator-depleted and Thomas reaction–diffusion models we prove the existence of a family of rectangles in the phase space that are invariant only under specific growth laws. Two numerical examples are provided to computationally demonstrate (i) the discrete maximum principle and optimal convergence for the heat equation on a linearly growing sphere and (ii) the existence of an invariant region for the LESFEM–IMEX Euler discretisation of a RDS on a logistically growing surface.


Introduction
Reaction-diffusion systems (RDSs) are a well-known class of partial differential equations that arise from the mathematical modelling of numerous phenomena taking place on a spacetime domain, see for instance [31]. In many applications the spatial domain is a curved surface rather than a flat domain, which may be time-dependent. Among the applications of RDSs on surfaces we mention brain growth [26], cell migration [3], chemotaxis [12], developmental biology [28], electrodeposition [24] and phase field modeling [42]. The growing interest toward PDEs on evolving surfaces has stimulated the development of several numerical methods for such problems, among which we mention (but not limited to) embedding methods [2], kernel methods [18], implicit boundary integral methods [5,35], surface finite element methods (SFEM) [10] and some of their recent variations and extensions [13,16,17,20,23,40].
An interesting property of some RDSs is the existence of invariant regions. An invariant region is a subset of the phase-space such that, if the initial condition takes values in , then the solution of the RDS takes values in at all times. We recall that, for scalar equations, the well-known notion of maximum principle is equivalent to the invariance of all the regions of the form [0, M], M > 0. For RDSs on a stationary surface, sufficient conditions have been found for a region to be invariant at the continuous level, see [39]. To the best of the authors' knowledge, the extension of these results to RDSs on evolving surfaces has not been considered in the literature. In this paper we will focus on surfaces that evolve according to a prescribed material velocity field. More complicated forms of the evolution law are beyond the scope of this manuscript and are a subject of our current studies.
For a numerical method it is interesting to understand if the invariant regions of the continuous problem are preserved under discretisation. In [15,16] we proved that, for RDSs on stationary surfaces, the lumped surface finite element method (LSFEM), combined with an implicit-explicit (IMEX) Euler method, preserves the invariant regions of the continuous problem. The purpose of the present paper is to extend these results to the case of evolving surfaces, in particular (i) we prove a semi-and a fully-discrete maximum principle for the heat equation with a linear source term and (ii) we provide sufficient conditions under which a region is invariant at the semi-and fully-discrete levels when a lumped evolving surface finite element method (LESFEM) and an IMEX Euler timestepping are considered. In particular we quantify the impact of surface evolution (measured through the dilation rate) on the existence of invariant regions and we find that surface growth or contraction respectively fosters or inhibits the invariance of a given region in the phase-space.
Crucial in our analysis is the assumption that the mesh preserves the Delaunay property under evolution, which is not true for an arbitrary surface evolution law. A class of surface evolution laws for which the Delaunay property is automatically preserved is that of isotropic growth [29], which has biological applications [7,27,33,36]. To the best of the authors' knowledge, an adaptive strategy for the preservation of the Delaunay property under a generic evolution law is still an open problem. An attempt in this direction is the work in [22]. For the special case of isotropic growth, we provide fully practical sufficient conditions for the existence of invariant regions at the semi-and fully-discrete levels. As an application of our general theory, we classify some classes of invariant regions, depending on the growth rate of the evolving surface, for two well-known RD models in the literature: the activator-depleted (or Schnakenberg, also known as the Brusselator model) and the Thomas models. Finally, we provide two numerical examples. In the first example we experimentally show that the LESFEM-IMEX Euler method, applied to the heat equation with a linear source term on a linearly growing sphere, exhibits optimal convergence rates in space and time. In the second example we consider the Thomas RDS on an exponentially growing Dupin ring cyclide thereby showing (i) the existence of invariant regions for the fully discretised model and (ii) the violation of this region in the absence of mass lumping.
The present paper is structured as follows. First, in Sect. 2 we recall (i) the derivation of RDSs on evolving surfaces and (ii) some basic notions about invariant regions and we introduce the notion of dilation rate, which is crucial in our analysis. In Sect. 3 we introduce the LESFEM for the space discretisation of RDSs on evolving surfaces and we carry out a fullydiscrete scheme using the IMEX Euler timestepping. Section 4 deals with the characterisation of the dilation rate in terms of the material velocity. In particular, we compute exactly the dilation rate for the class of isotropic growth laws. In Sect. 5 we prove a semi-and a fullydiscrete maximum principle for the linear heat equation on evolving surfaces with a linear source term. We prove, in Sect. 6, sufficient conditions for the existence of invariant regions for RDSs of arbitrarily many equations on evolving surfaces at the semi-and fully-discrete levels. Section 7 presents some classes of invariant regions for the activator-depleted and the Thomas RD models on evolving surfaces, respectively. Numerical examples are presented in Sect. 8. Finally, in Sect. 9 we conclude and discuss our findings with an eye for future extensions of the present work.

Preliminaries and Basic Results
For t ∈ [0, T ], let Γ (t) be a C 2 orientable surface in R 3 , represented as the zero-level set of a signed distance function d Hence, the outward unit normal vector field on Γ (t) is given by where · denotes the Euclidean norm in R 3 . Following [10, Section 5], we assume that there exists a mapping G : Vice-versa, ifṽ ∈ C 1 ([0, T ], C 2 (R 3 )) is an extension of v, i.e.ṽ(G(x 0 , t), t) = v(G(x 0 , t), t) for x 0 ∈ Γ (0) and t ∈ [0, T ], the mapping G (and thus the time-dependent surface Γ (t)) is recovered by solving, for each x 0 ∈ Γ (0), the Cauchy problem For t ∈ [0, T ] and δ > 0, let U δ (t) be the open neighbourhood of Γ (t) defined by We recall from [10] the following property.
Lemma 1 (Fermi coordinates) For any t ∈ [0, T ] there exists δ > 0 such that for any x ∈ Γ (t) there exists a unique a(x) ∈ Γ (t) that fulfils Hence, for t ∈ [0, T ], every point x ∈ U δ (t) can be described by d(x) and a(x), which are called normal coordinates or Fermi coordinates of x. The function a(·, t) : For any sufficiently smooth g : G T → R and g : G T → R 3 , let ∇ Γ (t) g, Δ Γ (t) g and ∇ Γ (t) · g denote the tangential gradient, the Laplace-Beltrami operator and the tangential divergence of g on Γ (t), respectively (see [10] for the details). In the following, we will write Γ instead of Γ (t) to simplify the notation. Furthermore, let ∂ • g denote the material derivative of g defined by where ∇ is the standard gradient in R 3 andg is any differentiable extension of g defined on a neighborhood of G T . Definition (6) is intrinsic, i.e. it does not depend on the choice of the extensiong (see [11] for further details). Let us recall some basic results from [10].

Derivation of the Reaction-Diffusion Model in Strong Form
Suppose we are given r ∈ N species u k : Γ (t) → R, k = 1, . . . , r , and let q k : Γ (t) → R 3 , k = 1, . . . , r , be their fluxes tangent to Γ (t). We recall from [1] the derivation of a system of r equations for u := (u 1 , . . . , u r ) that accounts for (i) the diffusion on the surface, (ii) the flux across the boundary (if non-empty) and (iii) the production rates f k (u), k = 1, . . . , r , of the given species. To this end, let R(0) be a portion of Γ (0) and let R(t) = G(R(0), t) be the portion of Γ (t) corresponding to the initial portion R(0). Notice that R(t) is itself an evolving surface, hence formulae (7) through (10) still hold if Γ (t) is replaced by R(t). We consider a mass balance on R(t) of the form Since the fluxes q k , k = 1, . . . , r , are tangent to Γ (t), we can apply the integration-by-parts formula (8) to the first term on the right hand side of (11). Then (11) becomes By applying the transport formula (9) to the left hand side of (12), we obtain where v is the material velocity defined in Sect. 2.1. Since R(t) is an arbitrary portion, we conclude that We assume q k corresponds to a diffusive flux according to Fick's law as follows: where d k , k = 1, . . . , r , are positive diffusivity constants. Inserting (15) into (14), we end up with the reaction-diffusion system of the form

Invariant Regions and Maximum Principle
In this section we recall basic notions concerning invariant regions for systems of the form (16) and conjecture a sufficient condition under which system (16) possesses an invariant region. To this end, we give the following definitions.

Definition 1 (Dilation rates)
The minimum and maximum instantaneous dilation rates are defined by respectively. When the minimum and maximum global dilation rates coincide, we call μ * := μ * min = μ * max the global dilation rate.
Definition 2 (Invariant regions) For the system (16), a region in the phase-space R r is said to be an invariant region if, whenever the initial condition u 0 is in , the solution u stays in as long as it exists and is unique.
We focus our attention on regions ⊂ R r of hyper rectangular shape, that is to say of the form where σ k ∈ R ∪ {−∞} and σ k ∈ R ∪ {+∞} for all k = 1, . . . , r . For instance, if σ k = 0 and σ k = +∞ for all k = 1, . . . , r , then is the positive orthant in R r , which means that the solution of the RDS stays positive at all times. Consider the (r − 1)-dimensional hyperfaces For k = 1, . . . , r , we define the constants where μ * min and μ * max are defined in (18). Next, we conjecture a criterion under which a hyper-rectangle is invariant for system (16). This criterion holds true in the stationary cases (when μ * = 0): (i) when Γ is a stationary monodimensional domain in R (see [38]), (ii) when Γ is a stationary k-dimensional domain in R k , k ∈ N (see [6]) and (iii) when Γ is a stationary Riemannian manifold without boundary (see [39]). In the case of isotropically evolving flat domains, the invariance of the positive orthant was studied in [41]. To the best of the author's knowledge, the case of evolving surfaces has not been studied at the continuous level. Hence, we introduce at the continuous level the following conjecture. (19) in the phase space of (16) and let f be Lipschitz on . If

Conjecture 1 Let be a hyper-rectangle as in
then is an invariant region for (16). In particular, when σ k = +∞ and σ k = −∞, then k ∩ R r and k ∩ R r are respectively empty, and so (21) and (22) are automatically fulfilled, respectively.
Notice that on stationary surfaces, since μ * k = μ * k = 0, then conditions (21)- (22) reduce to the inward flux condition considered in [39,Chapter 4]. In Sect. 4 next, we will prove the discrete counterpart of Conjecture 1 obtained by discretizing the RDS in space with the lumped version of the ESFEM method [10] and IMEX Euler in time.
We remark that, on stationary domains, some systems are known to possess an invariant region which do not meet the strict inequalities (21)- (22). For instance, for many mass-action laws, the positive orthant is invariant [4,16] even though the flux of f is tangent to this region, instead of strictly inward. For the scalar case k = 1, the notion of invariant region reduces to the well-known concept of maximum principle. Definition 3 (Maximum principle) For k = 1, the scalar equation (16) fulfils the maximum principle if, for any initial condition, i.e. u(·, 0), the solution u fulfils min 0, min In particular, if the initial condition is nonnegative, u(y, 0) ≥ 0 for y ∈ Γ (0), the solution fulfils Notice that the maximum principle corresponds to the fact that every (monodimensional) region of the form = [σ , σ ], with σ ≤ 0, σ ≥ 0, is invariant.

Derivation of the Variational Formulation
Following [1], we derive the variational formulation of system (16). To this end, for each t ∈ [0, T ] we multiply equations (16) by the respective test functions ) and integrate over Γ (t): for all k = 1, . . . , r . For the rigorous definition of the Sobolev spaces H 1 and H −1 on a manifold see [21], while for the Bochner spaces L 2 ([0, T ]; B), with B being any Banach space, see [34]. By applying the Green formula (10) to the right hand side of (25) we obtain for all k = 1, . . . , r . We assume that either Γ (t) has no boundary, i.e. ∂Γ (t) = ∅, or homogeneous Neumann boundary condition are enforced, i.e. ∇ Γ u k · μ = 0 on ∂Γ (t), so that the last term in (26) vanishes (see Remark 1). Furthermore, by observing that for all k = 1, . . . , r . By applying the transport property to the first term on the left hand side of (27), we have for all k = 1, . . . , r . Therefore, the variational formulation seeks to find u 1 , . . . , ).

Lumped Evolving Surface Finite Element Method
Following the evolving surface finite element method (ESFEM) studied in [1] for the approximation of the variational problem (29), we present its lumped counterpart, the lumped evolving surface finite element method (LESFEM).

Surface Triangulation and Some Definitions
Following We assume that, for each t ∈ [0, T ], Γ h (t) meets the Delaunay condition, defined as follows.
Let e be an edge of Γ h (t) and let Z 1 and Z 2 be the faces of Γ h (t) sharing the edge e. Let α 1 and α 2 be the angles in Z 1 and Z 2 opposite to e, respectively. The triangulation Γ h (t) is said to meet the Delaunay condition if, for all t ∈ [0, T ] and for any edge e of Γ h (t), it holds that see for instance [16]. The space-time triangulated surface G h,T is defined by Let V h be the space of time-dependent, spatially piecewise linear functions defined by Given t ∈ [0, T ] and a function η The discrete material derivative of a sufficiently smooth function U ∈ V h is defined by where v is the material velocity. For our purposes, we define the discrete counterpart of the dilation rates introduced in Definition 1.

Definition 4 (Discrete dilation rates)
The minimum and maximum discrete instantaneous dilation rates are defined by respectively. When the minimum and maximum discrete instantaneous dilation rates coincide, we call H (t) := H min (t) = H max (t) the discrete instantaneous dilation rate. The minimum and maximum discrete global dilation rates are defined by respectively. When the minimum and maximum discrete global dilation rates coincide, we call μ := μ min = μ max the discrete global dilation rate.
For every i = 1, . . . , N , the i-th Lagrange basis function χ i is the unique V h function such that where δ i j is the usual Kronecker symbol. The components U 1 , . . . , U r ∈ V h of the spatially discrete solution may be expressed in the Lagrange basis as

Preliminary Results on Triangulated Surfaces
We recall from [9] the following property of the basis functions.

Lemma 5 (Transport property of the basis functions) The basis functions
Hence, for the functions U k , k = 1, . . . , r defined in (36) it holds that We recall from [10, Lemma 5.6 and Remark 5.7] the following preliminary result.

Lumped Evolving Surface Finite Element Method
The lumped evolving surface finite element method (LESFEM), applied to the variational formulation (29), seeks to find for all k = 1, . . . , r and i = 1, . . . , N . Thanks to the transport property (37) of the basis functions, formulation (40) is equivalent to: find for all k = 1, . . . , r and i = 1, . . . , N . The LESFEM method (41) differs from the evolving surface finite element method (ESFEM) in [1] due to the presence of the interpolant operator on the first and the last terms in (41). By expressing U 1 , . . . , U r according to (36), the matrix form of (41) is where A and M are the (time-dependent) stiffness and lumped mass matrices defined by The above fact, together with the structure of the lumped mass matrix implies that, for any s ≥ 0, See [16] for further details.

Time Discretisation
We are now concerned with the time discretisation of the spatially discrete system (42) arising from the LESFEM. We discretise system (42) by means of the IMEX (IMplicit-EXplicit) Euler method, i.e by treating diffusion implicitly and reaction terms explicitly. To this end, let τ > 0 be a timestep, let t n := nτ for all n = 0, . . . , N T with N T = T τ , let A n and M n be the stiffness and lumped mass matrices at time t n , respectively, let (ξ n 1 , . . . , ξ n r ) be the coefficients of the numerical solution at time t n , and let f n k := f k (ξ n 1 , . . . , ξ n n ) for each k = 1, . . . , r . If (ξ 0 1 , . . . , ξ 0 r ) are the coefficients of the spatially discrete initial datum, the IMEX Euler time discretisation of (42) is or equivalently

Characterisation of Surface Growth
The purpose of this section is to characterise surface growth in terms of the material velocity v, with specific regard to isotropic growth. In fact, the lumped mass M(t) and stiffness A(t) matrices depend on v. In particular: 1. For an arbitrary triangulated surface that evolves with an arbitrary material velocity, we bound the time derivative dM dt of the lumped mass matrix in terms of the constants μ min and μ max defined in (34), i.e. in terms of the divergence ∇ Γ h · I h (v) of the discrete material velocity. We will need this result to prove a sufficient condition for the existence of invariant regions for the semi-and fully-discrete schemes; 2. For an arbitrary smooth or triangulated surface that evolves with an arbitrary material velocity, we characterise the velocity flows ∇ Γ ·v and ∇ Γ h · I h (v) in terms of the mappings G and G h introduced in Sects. 2.1 and 3.1, respectively; 3. For an arbitrary smooth or triangulated surface that evolves isotropically in space, that is where S : [0, T ] → R is an arbitrary smooth function, we compute exactly ∇ Γ · v and This result will yield a fully practical criterion to detect the invariant regions of a given RDS in the case of isotropic surface evolution.

Bounding the Rate of Change of the Mass Matrix in Terms of the Dilation Rates
In this section we bound the time derivative dM dt of M in terms of the discrete dilation rates defined in (34). To this end, we prove the following characterisation of dM dt . Lemma 7 (Transport formula for the lumped mass matrix M) The entries of the lumped mass matrix M fulfil the following property Proof By choosing U = I h (χ i χ j ) for any i, j = 1, . . . , N and V = 1 in the Leibniz formula (39) we have In some proofs we will need the following corollary of Lemma 7.
Corollary 1 (Consequence of the transport formula for the lumped mass matrix M) The diagonal matrix dM dt fulfils the estimates where μ min and μ max are defined in (34).

Characterising Velocity Flows in Terms of the Mappings G and G h
We wish to characterise the continuous and discrete velocity flows Letê i , i = 1, 2, 3, be the standard basis vectors in R 3 . For (θ, t) ∈ B × [0, T ], let J (x, t) ∈ R 2,3 be the (spatial) Jacobian of the function G(X (θ ), t) and let The surface integrals in (52) can be written as integrals on the planar domain B by using the parametrisation G(X (·), t) : where · denotes the Euclidean norm in R 3 , or equivalently, Since (55) holds for any measurable subset B of A, then it holds that d dt det(J (θ, t)) = ∇ Γ · v(G(X (θ ), t)) det (J (θ, t)) .
By applying the chain rule, (56) is equivalent to Given any triangulated surface Γ h (t) (which evolves under the discrete velocity field I h (v)), we notice that every facet Z (t) ∈ Z h (t) is smooth and thus parametrisable. For any facet of the initial triangulated surface, A similar reason as above yields the following discrete counterpart of (57).
Relations (57) and (59) are useful in that they (i) express the velocity flow without tangential derivatives and (ii) can be computed exactly when G and G h are known explicitly, e.g. for the isotropic growth as discussed in the next subsection.

Computing the Dilation Rates for the Isotropic Growth
In this section we compute the dilation rates defined in (18) and (34), respectively, for arbitrary surfaces that evolve with the material velocity (48). The velocity field (48) corresponds to the specific case of isotropic evolution, see for instance [8,29] for the case of evolving planar domains and [1,32] for the general case. In particular, for suitable choices of the function S(t), the growth law (48) admits some specific cases such as uniform, exponential, logistic and periodic growth, see [1,8,29,32]. From (3), it is easy to show that, with the velocity field (48), each x 0 ∈ Γ 0 evolves to the point therefore the evolution induced by an isotropic growth is a time-dependent dilation of the initial surface. The function that appears in (60) is known as the growth function, see for instance [29].
Remark 2 (Properties of isotropic growth) Isotropic growth preserves the angles of triangulated surfaces. This has two consequences: 1. if Γ h (0) meets the Delaunay condition, then Γ h (t) retains the Delaunay condition for all t ∈ [0, T ]; 2. if A(0) and M(0) are the stiffness and the mass matrices at t = 0, then Hence, in implementations, A(t)and M(t) need not be computed at each time step.
In the following result we compute the dilation rates μ min , μ max , μ * min and μ * max on an arbitrary smooth or triangulated surface that evolves with the material velocity (48), in terms of S(t).
Theorem 1 (Velocity flow on isotropically growing smooth or triangulated surfaces) Let Γ (t) be a smooth surface that evolves with the velocity field (48) and let Γ h (t) be the corresponding triangulated surface. Then, the instantaneous dilation rates satisfy Hence, it follows that Proof Let Γ h (t) be an evolving smooth surface and let (A, X ), J (θ, t) andJ (θ, t) as defined in Sect. 4.2. From (60) and (61), J (θ, t) fulfils where J X : A → R 2,3 is the Jacobian of X . It follows that which implies that By differentiating (67), we have By combining (57), (68) and dropping the parameterised coordinates θ , we have which proves the first equality in (63). Analogously, we prove the second equality in (63) by using (59). This completes the proof.
For the uniform, exponential, logistic and periodic growths, the functions φ(t), S(t) and the dilation rates μ min , μ * min , μ max and μ * max are detailed in Table 1, see also [29, Table 1] for the case of evolving planar domains, while the corresponding plots are depicted in Fig. 1.

Linear Heat Equation and Discrete Maximum Principles
We consider, for k = 1, the specific case of linear heat equation on an evolving surface Γ (t): and we prove the semi-and fully-discrete maximum principles for the case when μ min + β ≥ 0. Equation (70) is a special case of the general system (16) that we are interested in. However, we start with this specific case as (i) it provides more insights on the effect of growth on stability, (ii) we are able to prove a better timestep stability condition and (iii) to make the reader familiar with the demonstrative techniques.
The following result, that we proved in [16, Theorem 2.1] for the special case of stationary surfaces, addresses the maximum principle at the semi-discrete level.
The constant r > 0 is the growth rate. For the logistic growth, K > 1 is the carrying capacity, i.e. the square root of the ratio between the asymptotical and the initial area of the surface (see [29])   Table 1 for K = 2 and r = 0.5, 1, 1.5. From left to right: linear, exponential, logistic and periodic growth profiles Theorem 2 (Semi-discrete maximum principle for the linear heat equation (70)) If the velocity field v fulfils with μ min as defined in (34), and the triangulation Γ h (t) meets the Delaunay condition for all t ∈ [0, T ], then the LESFEM solution of (70) fulfils the following discrete maximum principle Proof From (42), the LESFEM spatial discretisation of (70) is By applying the chain rule, (73) becomes By multiplying (74) on the left by M −1 , we havė All we have to prove is that the ODE (75) is dissipative, i.e. −d M −1 A|ξ | − M −1 dM dt |ξ | − β|ξ | ≤ 0. For every t ∈ [0, T ], M is diagonal with positive diagonal entries and A fulfils (43) from the Delaunay condition. Then, it follows that −d M −1 A|ξ | ≤ 0. Hence, it suffices to prove that −(M −1 dM dt + β I )|ξ | ≤ 0, that is true provided By using (51), condition (76) is true if which holds true from assumption (71). This completes the proof.
The next theorem shows the same result for the LESFEM-IMEX Euler full-discretisation of (70), under a timestep restriction. This result holds true for the special case of stationary surfaces (see [16, Theorem 2.2]).

Theorem 3 (Fully-discrete maximum principle for the linear heat equation (70)) If the velocity field v fulfils
with μ min as defined in (34), if the timestep satisfies In particular, there is no timestep restriction if β ≤ 0.

Remark 3 (Interplay between material velocity and source term) Relation (71) implies that
-domain growth (μ min > 0) can enable the discrete maximum principle even for β < 0; -local domain contraction (μ min < 0) can prevent the discrete maximum principle even for β ≥ 0.
This interplay is justified by observing that domain evolution implies a dilution effect, explained as follows. By choosing ϕ = 1 in the variational formulation (29) with k = 1 and f 1 (u) = −βu, we obtain If |Γ (t)| denotes the surface area of Γ (t) and u(t) := 1 |Γ (t)| Γ (t) u denotes the mean value of u, (90) becomes By solving (91) for d dt u(t) , we obtain By choosing g = 1 in the transport formula (9), we have d dt By combining (92) and (93) we have From (94), the dilution effect arising from surface growth can be interpreted as the dampening or uplifting effect of μ * min on u(t) . The estimate (94)

implies that u(t) is non-increasing if
which is the continuous counterpart of (71). We conclude that condition (71) is consistent with the interpretation of surface growth in terms of dilution effect.
Remark 4 (Interplay between timestep restriction and source term) Relation (80) implies that the timestep restriction needed for guaranteeing the discrete maximum principle is independent of the material velocity and it only depends on the stiffness parameter β of the source term. In particular, when the source term is nonnegative (i.e. when β ≤ 0), the LESFEM-IMEX Euler fully-discrete scheme unconditionally fulfils the discrete maximum principle.

Reaction-Diffusion Systems and Invariant Regions
In this section we prove, for the semi-and full-discretisations of RDSs of the form (16), a criterion to test if a hyper-rectangle in the phase-space is invariant. In the case k = 1 of scalar equations, the notion of invariant region collapses to that of minimum-maximum principle, considered in the previous section for the special case of the linear heat equation. We assume that the Delaunay regularity of the mesh is preserved under evolution. For k = 1, . . . , r , we define the constants where μ min and μ max are the dilation rates defined in (34). In the following theorem we prove that, under similar assumptions of Conjecture 1, is an invariant region for the solution obtained from the semi-discrete scheme (42). Hence, the following theorem extends [16,Theorem 3.3] to the case of evolving surfaces. (42)) Let be a hyper-rectangle as in (19) in the phase space of (42), let f be Lipschitz on . If the triangulation Γ h (t) satisfies the Delaunay condition for all t ≥ 0 and

Theorem 4 (Invariant rectangles for
then is an invariant region for (42).
Proof The semi-discrete method (42) can be written, after applying the chain rule to the term d dt (Mξ k ) and multiplying on the left by M −1 aṡ Since M is diagonal with positive diagonal entries and A i j ≤ 0 for i = j from the assumption of Delaunay regularity, proceeding as in the proof of [16,Theorem 3.3], it suffices to verify that, for all (U 1 , . . . , U r ) ∈ , k = 1, . . . , r , and i = 1, . . . , N , where σ k and σ k are as in (19). Using relation (51), conditions (100)-(101) hold if, for all k = 1, . . . , r , with μ k and μ k as in (96), that is true from assumptions (97)-(98). This completes the proof.
The following theorem provides a sufficient condition for regions to be invariant for the LESFEM-IMEX Euler scheme (47) and extends [16,Theorem 3.4]. In contrast to the semidiscrete case, we relax the strict inequalities (21)- (22) with conditions (105)-(106), in which we use the perturbed dilation ratesμ k and μ k given bỹ respectively, and μ min and μ max are defined in (34). Observe thatμ k → μ k and μ k → μ k as τ → 0. (47)) Let be a hyper-rectangle as in (19) in the phase space of (42), let f be Lipschitz on . If the triangulation Γ h (t) meets the Delaunay condition for all t ≥ 0 and

Theorem 5 (Invariant rectangles for
then is an invariant region for (47) if the timestep τ fulfils where, for all k = 1, . . . , r, L k is the Lipschitz constant of f k .
Proof The fully-discrete scheme (47) can be written as n ∈ N ∪ {0}, k = 1, . . . , r . Since the mesh meets the Delaunay assumption at all times, the matrix properties (82)-(83) hold. Then, it suffices to prove that where 1 is the column vector of ones. We will prove the two inequalities in (109) in turn. From (51), the inequality on the right side of (109) holds true if with μ k as defined in (96). Suppose σ k ≥ 0. From assumption (105) we can estimate f n k as follows From (111), condition (110) holds true provided From assumption (107), since ξ n k ≤ σ k 1, then (112) holds true if that is to say which holds true for each τ ∈ R. Suppose, instead, σ k < 0. From assumption (105) we can estimate f n k as follows From (115), condition (110) holds true provided From assumption (107), since ξ n k ≤ σ k 1, then (116) holds true if As (117) always holds with the equality, we conclude that the second inequality in (109) is true under assumptions (105) and (107). Similarly, the inequality on the left side of (109) holds under assumptions (106) and (107). This completes the proof.

Remark 5 (Sharper timestep restriction)
In the specific case of the linear heat equation (70), estimate (80) in Theorem 3 is sharper than estimate (107) in Theorem 5. In fact, since the Lipschitz constant of the source term is L = |β|, the timestep restriction (107) is fulfilled for τ |β| ≤ 1, that is more restrictive than condition (80).

Velocity-Induced Invariant Regions for RD Models
Now, we consider two different RDSs that are well-known in the literature and prove, at the discrete level, the existence of discrete invariant hyper-rectangles for these RDSs, depending on the global discrete dilation rates μ min and μ max defined in (18). The results in this section are confined to the spatially discrete level, but from Conjecture 1, we claim that the same results holds at the continuous level. In the special case of stationary surfaces (i.e. when μ min = μ max = 0), we obtain invariant hyper-rectangles that have been studied in the literature (see [4,6,19]). It is worth remarking that, even though we consider two RD models for illustrative purposes, the following analysis can be easily extended to other types of RDSs. To mention two examples, for the well-known Hodgkin-Huxley model and for the DIB model for electrodeposition considered in [24,25], the invariant region study is carried out in [14].

RDS with Activator-Depleted Kinetics
Let us consider an RDS with the well-known non-dimensional activator-depleted kinetics, also known as Schnakenberg or Brusselator kinetics (see for instance [1,37]), on evolving surfaces where a, b and γ are positive parameters and d is a positive diffusion rate. The model describes a system of two interacting chemicals, in which u 1 ≥ 0 and u 2 ≥ 0 are the respective concentrations. For this reason, we focus our attention on invariant regions contained in the positive ortant. In the following theorem we prove that: (i) the positive orthant is invariant for (118) regardless of μ min and μ max . At the continuous level, the result holds in the specific case of stationary planar domains, see [4]. (ii) when μ min > 0, the model possesses invariant stripes (depending on μ min ) in the positive orthant.
Theorem 6 (Velocity-induced invariant regions for the activator-depleted model (118)) For the LESFEM spatial discretisation of (118), the following statements hold: 1. For any value of the constants μ min , and μ max defined in (34), the positive orthant If μ min > 0 and σ 2 is a constant such that This proves Statement (1). For Statement (2), let μ min > 0 and we assume for the moment that the strict inequality holds in (119). Then the set 1 : This proves Statement (2) when the strict inequality holds in (119). Otherwise, observe that i.e. is the intersection of invariant regions and is thus invariant. This completes the proof of Statement (2).

RDS with Thomas kinetics
Let us consider an RDS with the non-dimensional Thomas kinetics (see for instance [31, p. 78]), on evolving surfaces where α, a, b, γ , K and ρ are positive constants and d is a positive diffusion rate. The model describes a system of two interacting chemicals, in which u 1 ≥ 0 and u 2 ≥ 0 are the respective concentrations. For this reason, we focus our attention on invariant regions contained in the positive orthant.

Theorem 7 (Velocity-induced invariant regions for the Thomas model (121))
For the LESFEM spatial discretisation of (121), the following statements hold: 1. For any value of the constants μ min , and μ max defined in (34), the positive orthant (1, α) and σ 1 and σ 2 are two constants such that Proof To prove Statements (1) and (2), we have to verify conditions (21)- (22). For Statement (1), observe that This proves Statement (1). For Statement (2), let μ min > −γ min(1, α) and we assume for the moment that the strict inequalities hold in (122). Then, observe that • the set 1 := {σ 1 } × [0, σ 2 ] is contained in the region , This proves Statement (2) when the strict inequalities hold in (122). Otherwise, we have that i.e. is the intersection of invariant regions and is thus invariant.

Numerical Examples
The purpose of this section is to provide two numerical examples in which we (i) estimate the experimental order of convergence of the LESFEM and (ii) experimentally show the ability of Theorem 5 to find invariant regions of RDSs on evolving surfaces at the discrete level.

Numerical Example 1: Linear Heat Equation on an Evolving Sphere
In this example, we wish to estimate the experimental order of convergence of the LESFEM. As a test problem, we consider the linear heat equation given by We choose T = 1 to be the final time. The initial domain Γ (0) is the unit sphere S 2 , that evolves under the velocity field v(x, t) : and undergoes linear growth for r = 1, see Table 1. In particular, the domain Γ (t) at time t ∈ [0, T ] is a sphere whose radius is given by the growth function φ(t) = t + 1 and the minimum dilation rate fulfils μ min = 2r r T +1 = 1 (see Table 1). In order to determine the experimental order of convergence, we consider the analytical solution to (124) given by  Fig. 3: the convergence is experimentally optimal in that it is quadratic in the meshsize and linear in the timestep.

Numerical Example 2: RDS with Thomas Kinetics and Invariant Regions
In this example, we show that the LESFEM-IMEX Euler preserves the invariant regions of RDSs on evolving surfaces. Let us consider the following RDS with Thomas kinetics  (124): Error in L ∞ (0, T, L 2 (Γ h (t))) norm (left) and experimental rate of convergence (right). The quadratic convergence in space is optimal considered in Sect. 7.2, with a = 150, b = 100, ρ = 13, K = 0.05 as in [30] and γ = 1 for illustrative purposes. With these parameters, system (127) an orientable surface without boundary that is topologically equivalent to a torus. The surface evolves under the velocity field (48) with S(t) = r K (K −1) K −1+exp(Krt) > 0, t ∈ [0, T ], with r = 0.2, K = 3 and undergoes a logistic growth, see Table 1. As final time we choose T = 100. From (129) It is easy to see that, on a region of the form = [σ 1 , σ 1 ] × [σ 2 , σ 2 ], the Lipschitz constants L 1 and L 2 of the kinetics f 1 and f 2 of (127) fulfil Hence, the timestep restriction (107) becomes τ ≤ 1 max(L 1 ,L 2 ) ≈ 1.676e-3. It follows that is invariant for the LESFEM-IMEX Euler full discretisation under the timestep restriction 0 < τ ≤ min τ , 1 max(L 1 ,L 2 ) ≈ 1.676e-3. We choose τ = 1e-3. The region is smaller than the invariant region provided in Theorem 7, which has a simple analytical expression but is not optimal. The following initial condition u 1,0 (x, y, z) = σ 1 + (σ 1 − σ 1 )ψ(x, y, z); where ψ(x, y, z) := 1 − 25(min(|y|, 1 5 )) 2 fulfils (u 1,0 (x), u 2,0 (x)) ∈ for all x ∈ Γ (0) and is shown in the first snap shot of Fig. 5. We solve the problem with τ = 1e-3 on a sequence of seven meshes Γ h,i (t), i = 1, . . . , 7, with decreasing initial meshsizes h i (0),   Fig. 5. In particular, at the final time T = 100 (see the last snap shot of Fig. 5), the surface is stationary up to machine precision and the numerical solution has reached a stationary pattern. For a given numerical solution (U 1 , U 2 ) on the mesh Γ h,i , i = 1, . . . , 7, consider the following functions These functions are the oriented distances of the numerical solution (U 1 , U 2 ) from the edges of . If the oriented distances η k and η k , k = 1, 2, stay positive at all times, it means that  The minima of η 1 (and thus the minima of U 1 ) coincide up to machine precision

Conclusions
In this paper we have considered a lumped evolving surface finite element method for the spatial discretisation of reaction-diffusion systems on evolving surfaces, by extending substantially the counterpart on stationary surfaces studied in [16]. We have obtained a fullydiscrete scheme by applying the IMEX Euler timestepping to the spatially discrete method. In Sect. 4 we have presented a characterisation of the tangential flow of the continuous and discrete material velocities on evolving smooth and triangulated surfaces, respectively. As a consequence we have obtained, in Theorem 1, an explicit expression for the dilation rates for the special case of isotropic growth. In Theorems 2 and 3 of Sect. 5 we have presented sufficient conditions for the linear heat equation on evolving surfaces to fulfil the maximum principle at the semi-and fully-discrete levels, respectively. In particular, at the fully-discrete level, no timestep restriction is needed if the source parameter β in (70) is nonpositive. In Theorems 4 and 5 of Sect. 6 we provided sufficient conditions under which a hyper-rectangle in the phase space of an RDS on evolving surfaces is invariant at the semi-and fully-discrete levels, respectively. In particular, at the fully-discrete level, a timestep restriction depending on the Lipschitz constants of the kinetics is needed. In Sect. 7 we classified some families of invariant regions for two well-known RD models on evolving surfaces: the activator-depleted and the Thomas models. Two numerical examples are presented in Sect. 8 that experimentally show (i) the optimal convergence (i.e. quadratic in space and linear in time) of the proposed method for the linear heat equation on a linearly evolving sphere and (ii) the existence of a bounded invariant region for an RDS with Thomas kinetics on a logistically growing Dupin ring cyclide at the discrete level. The mathematical and numerical analysis of more complicated settings such as nonisotropic chemically driven is a non-trivial exercise that requires new mathematical tools and techniques. Such a study is beyond the scope of this paper. In order to provide a first step in this direction, we have considered the simplest form of growth, uniform isotropic, which has allowed us to state precisely the conditions necessary for the preservation of invariant regions as well as stating the conditions for positivity of solutions.
In future studies we wish to extend the method by considering an ALE approach as in [13] in order to preserve good mesh properties under evolution. This would involve the study of convection RDSs on evolving surfaces, where the convection term arises from the ALE mapping.
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.