A Conservative Finite Element ALE Scheme for Mass-Conserving Reaction-Diffusion Equations on Evolving Two-Dimensional Domains

Mass-conservative reaction-diffusion systems have recently been proposed as a general framework to describe intracellular pattern formation. These systems have been used to model the conformational switching of proteins as they cycle from an inactive state in the cell cytoplasm, to an active state at the cell membrane. The active state then acts as input to downstream effectors. The paradigm of activation by recruitment to the membrane underpins a range of biological pathways - including G-protein signalling, growth control through Ras and PI 3-kinase, and cell polarity through Rac and Rho; all activate their targets by recruiting them from the cytoplasm to the membrane. Global mass conservation lies at the heart of these models reflecting the property that the total number of active and inactive forms, and targets, remains constant. Here we present a conservative arbitrary Lagrangian Eulerian (ALE) finite element method for the approximate solution of systems of bulk-surface reaction-diffusion equations on an evolving two-dimensional domain. Fundamental to the success of the method is the robust generation of bulk and surface meshes. For this purpose, we use a moving mesh partial differential equation (MMPDE) approach. Global conservation of the fully discrete finite element solution is established independently of the ALE velocity field and the time step size. The developed method is applied to model problems with known analytical solutions; these experiments indicate that the method is second-order accurate and globally conservative. The method is further applied to a model of a single cell migrating in the presence of an external chemotactic signal.

by guanine nucleotide dissociation inhibitors (GDI), which act as inhibitors of activation. The active forms, which are bound to guanosine triphosphate (GTP), bind strongly to effectors (for example the kinase PAK). As the effectors are recruited to the membrane they too become activated. Together the Rho GTPases and their effectors trigger downstream signalling events, such as actin polymerisation via the Arp2/3 complex, which lead to cell migration. Over the time scales of polarisation and cellular migration, the Rho family of GTPases are conserved quantities in that the total amount of active (membrane bound) and inactive (cytoplasmic and membrane bound) species remains constant. This suggests that the correct mathematical framework for modelling GTPase activity is mass-conservative reaction-diffusion (McRD) systems [25,26].
Early McRD models were based on the assumption that the cell could be represented as a onedimensional slice through the cell from front to back and that the different physical locations of the protein conformations were accounted for by the use of different diffusivities, i.e the membrane bound conformations having a smaller diffusivity compared to the cytoplasmic conformations. Examples of models of this type include Narang et al. [54], Subramanian and Narang [68], Otsuji [57], and Mori et al. [49]. Recent research has tried to more faithfully account for the different spatial locations of the inactive and active GTPases. So-called bulk-surface models consist of a combination of equations posed on the cell membrane (surface region), equations for cytoplasmic diffusion (bulk region), and Robin-type flux boundary conditions coupling the membrane bound and cytoplasmic species. For example, Cusseddu et al. [6] recently proposed a bulk-surface extension of the classical wave-pinning (WP) model [49] to investigate cell polarisation in three dimensions. Giese et al. [21,22] constructed a two-dimensional model to study the effect of cell size and cell shape on the signalling pathway from the cell membrane to the nucleus. Diegmiller et al. [10] used a bulk-surface model to look at pattern formation on the surface of a spherical geometry. The effect of geometry on a pattern formation in a two-dimensional setting was also investigated using a bulk-surface model by Thalmeier et al. [69]. All of the models have become amenable to numerical computation using recent developments in numerical techniques to approximate the solution of PDEs on surfaces [12,13,16,17,15,44,43].
One area where bulk-surface McRD models have not been applied extensively is the study of cell migration. The main reason for this is the additional computational challenges which arise from the time dependence of the bulk and surface domains and the need for the numerical solution to be massconservative. In [56] the authors presented a finite volume method using a fixed background mesh. As their method is based on a finite volume spatial discretisation, by construction it is both locally and globally conservative. However, the spatial rate of convergence was seen to be less than second order on a range of test cases with known analytical solutions. Strychalski et al. [65,66] used a cut-cell approach, again with a fixed background mesh, and applied their method to a conservative reaction-diffusion system on an evolving domain. However, conservation was not maintained at the discrete level although the error in conservation could be reduced by grid refinement. Other popular numerical methods include: immersed boundary methods (see e.g. [70,3,4]); and phase-field methods (see e.g. [36,1,52]).
One approach to ensure global conservation is to use domain fitted meshes which follow the evolving cell membrane and cytoplasm. In [37] we introduced a moving finite element method for the solution of bulk-surface PDEs on evolving two-dimensional bulk domains. An adaptive moving mesh based on the solution of moving mesh PDEs [29,31] was coupled to a novel aproach for moving a mesh on an evolving curve. The moving mesh method was further extended to higher degree of temporal accuracy and applied to forced curve shortening problems in [42]. The main aim of the work presented here is to investigate the mass conservation property of a fully discrete finite element approximation of an ALE reformulation of a McRD bulk-surface system on an evolving domain. For efficiency and accuracy it is also important that the grid generation procedure is robust. We therefore also present a refinement of the finite element technique in [37] to solve the MMPDEs which results in a much more robust procedure.
The layout of the rest of this paper is as follows: In Section 2 we introduce a class of bulk-surface reaction-diffusion systems on evolving domains and in Section 3 we present a weak ALE reformulation and show that the continuous system is globally conservative independent of the ALE velocity. A conservative FEM discretisation is introduced in Section 4 and a second-order two-stage time integration scheme is given in Section 5 for the numerical solution of the coupled ODE system obtained after spatial discretisation. The grid generation procedure is discussed in Section 6. In Section 7 we prove the discrete FEM scheme is globally conservative independently of the ALE velocity and time step size. In Section 8 some numerical experiments verifying the global conservation properties of the FEM scheme, as well as spatial convergence rates, are given. We finally apply the method to a bulk-surface form of the WP model [6] for cell polarisation and chemotaxis on a moving domain. Figure 1: We consider the simulation of a motile cell through a fixed lab frame of reference Λ. The exterior region is denoted by Ξ(t) and the interior region is denoted by Ω(t). The exterior and interior regions are separated by the curve Γ(t). After a time interval of size ∆t, the material point located at X p (t) on Γ(t) evolves to the new location X p (t + ∆t) on Γ(t + ∆t).

Reaction-diffusion systems on evolving domains
The physical layout for the cell migration simulations considered later is shown in Fig. 1. We assume that the cell moves through a fixed frame of reference Λ. Let the time interval be denoted by I = (0, T ] ⊂ R such that the closure is denoted I = [0, T ]. For each t ∈ I, let D(t) ⊂ R 2 be a simply connected domain containing a smooth, closed curve Γ(t) which separates D(t) into an interior region Ω(t) ⊂ R 2 and an exterior region Ξ(t) ⊂ R 2 such that D(t) = Ω(t) ∪ Ξ(t). Note that ∂Ω(t) = Γ(t) and ∂Ξ(t) = Γ(t) ∪ ∂D(t), where ∂D(t) denotes the outer boundary of the exterior domain Ξ(t). To improve computational efficiency, the governing equations for the extracellular region will only be computed over Ξ(t), which is centred on the centroid of Ω(t).
We will assume that a material particle P located at X P (t) ∈ Γ(t) has velocityẊ P (t) ∈ R 2 . Therefore, we assume that there exists a velocity field u Γ : Γ(t) × I → R 2 so that material points on Γ(t) evolve with a velocity fielḋ X P (t) = u Γ (X P (t), t).
Similarly, material points on the outer boundary ∂D(t), in the interior region and exterior region will be assumed to have material velocities u ∂D : ∂D(t) × I → R 2 , u Ω : Ω(t) × I → R 2 and u Ξ : Ξ(t) × I → R 2 , respectively.
Let n = (n 1 , n 2 ) denote the unit outward normal to Γ(t) and let N (t) be any open subset of R 2 containing Γ(t). For any function ζ which is differentiable in N (t), we define the tangential gradient on Γ(t) by ∇ Γ ζ = ∇ζ − (∇ζ · n)n, where · denotes the usual scalar product and ∇ζ denotes the usual gradient on R 2 . For a vector function ζ = (ζ 1 , ζ 2 ) ∈ R 2 , the tangential divergence is defined by The Laplace-Beltrami operator on Γ(t) is defined as the tangential divergence of the tangential gradient ∆ Γ ζ = ∇ Γ · (∇ Γ ζ). Furthermore, let us define The closure of the sets defined in (2.1) will be denoted Q Ξ T , Q Ω T , Q Γ T and Q ∂D T , respectively, where it is understood that the closure occurs in both space and time.
In the exterior domain Ξ(t), we consider N l ∈ N freely diffusing chemical species whose concentrations are given by l (i) : Q Ξ T → R, i = 1, . . . , N l . The vector of concentrations of the N l exterior bulk species is given by Excluding cross-diffusion and coupling the reaction-diffusion systems by the reaction-kinetics only, these concentrations evolve according to the general advection-reaction-diffusion equation l is the diffusion coefficient of the i-th exterior chemical species; the material velocities of Ξ(t), Γ(t) and ∂D(t) are denoted by u Ξ , u Γ and u ∂D , respectively; the function f (i) l : Q Ξ T × R N l → R describes the reaction kinetics between the interacting exterior bulk species l; the function g (i) : Q Γ T × R N l → R describes the boundary conditions for the interacting exterior bulk species; and the functionĝ (i) : Q Γ T × R N l × R N ls → R describes the interaction between the exterior bulk (l) and surface (l s ) species. The normals n D and n Ξ are the outward unit normals to Ξ(t) from ∂D(t) and Γ(t), respectively.
On the surface Γ(t), we consider N ls ∈ N chemical species that are free to diffuse tangentially along the surface Γ(t) and whose concentrations are given by l (j) s : Q Γ T → R, j = 1, . . . , N ls . The vector of concentrations of the surface species is given by These concentrations evolve according to the general surface advection-reaction-diffusion equation ls is the diffusion coefficient of the j-th surface chemical species; the function f (j) ls : Q Γ T × R N ls → R describes the reaction kinetics between the interacting surface species; and the function g (j) s : Q Γ T × R N l × R N ls → R describes the interaction between the exterior bulk (l) and surface (l s ) species.
Similarly, in the interior domain we consider N c ∈ N freely diffusing chemical species whose concentrations are given by c (p) : Q Ω T → R, p = 1, . . . , N c . The vector of concentrations of the N c interior bulk species is given by These concentrations evolve according to the general advection-reaction-diffusion equation c is the diffusion coefficient of the p-th interior chemical species; the material velocity of Ω(t) is denoted by u Ω ; the function f (p) c : Q Ω T × R Nc → R describes the reaction kinetics between the interacting interior bulk species c; the function r (p) : Q Γ T × R Nc → R describes the boundary conditions for the interacting interior bulk species; and the functionr (p) : Q Γ T × R Nc × R Nc s × R N ls → R describes the interaction between the interior bulk (c) and surface (c s ) species under stimulation from the exterior surface species (l s ). The outward unit normal to Ω(t) from Γ(t) is denoted by n Ω .
Finally, on the surface Γ(t) we consider N cs ∈ N chemical species that are free to diffuse tangentially along the surface Γ(t) and whose concentrations are given by c (q) s : Q Γ T → R, q = 1, . . . , N cs . The vector of concentrations of the surface species is given by These concentrations evolve according to the general surface advection-reaction-diffusion equation cs is the diffusion coefficient of the q-th surface species; the function f (q) cs : Q Γ T × R Nc s → R describes the reaction kinetics between the interacting surface species; and the functionr (q) s : Q Γ T ×R Nc × R Nc s × R N ls → R describes the interaction between the interior bulk (c) and surface (c s ) species under stimulation from the exterior surface species (l s ).

Weak ALE formulation
When the domain D(t) is moving it is conventional to adopt a common frame of reference for computational purposes. A popular choice is the Arbitrary Lagrangian Eulerian (ALE) frame. Let A(·, t) : D c → D(t) be a family of bijective mappings, which at each t ∈ I map points in a reference or computational configuration D c ⊂ R 2 with coordinates ξ = (ξ, η) to points in the current physical configuration Here we assume that the computational configuration is the initial configuration D c = D(0).
Let ψ : D(t) × I → R be an arbitrary function defined on the fixed Eulerian frame and φ : D c × I → R be the corresponding function defined on the ALE frame, such that φ(ξ, t) = ψ(A(ξ, t), t). Then the temporal derivative of ψ with respect to the ALE frame is defined as A standard application of the chain rule yields where the time derivative of the ALE mapping defines the ALE velocity w as Note that in general, the ALE velocity will differ from the material velocities. The reformulation of (2.2a) and (2.3) in terms of the ALE reference frame takes the form for each i = 1, . . . , N l , and for each j = 1, . . . , N ls . Similarly, the reformulation of (2.4a) and (2.5) in terms of the ALE reference frame takes the form for each p = 1, . . . , N c , and for each q = 1, . . . , N cs . For full details of the ALE formulation, the reader is referred to MacDonald et al. [37].
To construct a weak formulation of (3.1), we assume Dirichlet boundary conditions on the outer boundary ∂D(t) and consider the space of admissible test functions defined on the reference domain made of functionsv : Hence, the ALE mapping defines a set H(Ξ(t)) of test functions v : Ξ(t) → R as follows: For each i = 1, . . . , N l , multiply (3.1) by a test function v ∈ H(Ξ(t)) and integrate over Ξ(t) to give the weak form: find l (i) ∈ H(Ξ(t)) such that For each p = 1, . . . , N c , multiply (3.3) by a test function v ∈ H(Ω(t)) and integrate over Ω(t) to give the weak form: find c (p) ∈ H(Ω(t)) such that where we have applied the flux boundary conditions (2.4b). Finally, to construct the weak formulation of (3.2) and (3.4), consider the space of admissible test functions defined on the reference curve made of functionsv s : Γ c → R,v s ∈ H 1 (Γ c ). Once again, the ALE mapping defines a set H s (Γ(t)) of test functions v s : Γ(t) → R, as follows: For each j = 1, . . . , N ls and q = 1, . . . , N cs , multiply (3.2) and (3.4) by a test function v s ∈ H s (Γ(t)) and integrate over Γ(t) to give the weak form: find l (j) s ∈ H s (Γ(t)) and c (q) and s (c, c s , l s )v s ds, ∀v s ∈ H s (Γ(t)), (3.8) where we have assumed tangential surface fluxes only.

Global conservation of the continuous formulation on evolving domains
When Dirichlet boundary conditions are used on the exterior outer boundary, global conservation cannot be expected. Global conservation can be obtained if zero flux boundary conditions are employed in (2.2b). However, in order to maintain a homogeneous distribution of exterior chemical species under movement, one would require a non-zero exterior material velocity u Ξ , which may be unphysical depending on the problem under consideration. Finally, one could employ a non-zero flux boundary condition on the exterior outer boundary. However, such a boundary condition will in general be unknown a-priori and conservation would only be achieved if Therefore, in the results presented in §8.2 we will employ Dirichlet boundary conditions on the exterior outer boundary and, in the studies of global conservation throughout this article, focus our attention on the interior bulk-surface environment only; that is, global conservation of c and c s . We require the following assumption: Assumption 3.1. We assume that on the boundary Γ(t), the normal components of the ALE velocity and material velocity are identical; that is, w · n = u Γ · n.
The shape of an evolving curve is determined purely by its normal velocity and therefore, Assumption 3.1 is necessary to ensure that the evolution of the curve with respect to the ALE frame matches the evolution of the curve prescribed by the material velocity u Γ . We are now ready to prove global conservation of the reaction-diffusion system given by (2.4) and (2.5).  are satisfied, then the reaction-diffusion system given by (2.4) and (2.5) is globally conservative independent of the ALE velocity; that is, the following holds true: Proof. In (3.6) we set v = 1, apply the divergence theorem and sum over all p = 1, . . . , N c to give Similarly, in (3.8) we set v s = 1 and sum over all q = 1, . . . , N cs to give s (c, c s , l s ) ds. (3.13) Due to Assumption 3.1, the normal component of the ALE and material velocites are the same. Hence, the first term on the right-hand side of the above is zero. Also, due to the curve Γ(t) being closed there is no boundary of Γ(t) and therefore, there cannot be any tangential flux; in other words, Therefore, (3.13) reduces to The result is then immediate following the application of (3.9) and (3.10).

Finite element spatial discretisation
We will assume that the reference D c and physical D(t) domains are approximated by polygonal domains D c,h and D h (t), ∀t ∈ I, respectively. Let D c,h be covered by a fixed triangulation T D h,c with straight edges, so that The curve Γ(t) separates D(t) into exterior Ξ(t) and interior Ω(t) regions. Therefore, we construct the triangulation T D h,c so that it is fitted to the curve Γ(t). To do this, we let Ξ c,h and Ω c,h be covered by fixed triangulations T Ξ h,c and T Ω h,c , respectively, with straight edges so that and Let E denote the set of all edges of the triangulation T Ω h,c and define Thus, the curve Γ c is approximated using the boundary of Ω c,h so that Γ c,h := ∂Ω c,h . Note that due to T D h,c being fitted to Γ c , we could also define Γ c,h using the inner boundary of Ξ c,h . Similarly, the approximation to the outer boundary ∂D c is obtained using the outer boundary of Ξ c,h . The number of elements of T D h,c , T Ξ h,c and T Ω h,c will be denoted by N D e , N Ξ e and N Ω e , respectively; the total number of vertices of T D h,c , T Ξ h,c and T Ω h,c will be denoted by N D , N Ξ and N Ω , respectively; and the number of vertices on the boundaries will be denoted by N Γ s and N ∂D s . The finite element spaces on T D h,c can then be defined as wherev h : D c,h → R and P 1 (K) is the space of linear polynomials on K. Similarly, the finite element spaces on T Ξ h,c and T Ω h,c can be defined as respectively. Furthermore, we define The ALE mapping is interpolated using piecewise linear finite elements giving rise to a spatially discrete denotes the position of the µ-th node at time t, and {φ µ (ξ)} N D µ=1 are the nodal basis functions in L 1 (D c,h ). The discretised ALE velocity therefore takes the form The physical triangulation is defined as the image of the reference triangulation under the discrete ALE mapping: for example, , t ∈ I. The finite element spatial discretisation of the ALE formulation (3.5) therefore takes the form: find l Similarly, the ALE mapping defines the finite element test space on Ω h (t) as the set The finite element spatial discretisation of the ALE formulation (3.6) therefore takes the form: find c Finally, the ALE mapping defines the finite element test space on Γ h (t) as the set , t ∈ I. The finite element spatial discretisation of the ALE formulations (3.7) and (3.8) take the form: find l (j) s,h ∈ H s,h (Γ h (t)), for each j = 1, . . . , N ls , and c For each i = 1, . . . , N l , j = 1, . . . , N ls , p = 1, . . . , N c and q = 1, . . . , N cs , the finite element approximation of the bulk and surface species can be expressed as such that, for each i = 1, . . . , N l , j = 1, . . . , N ls , p = 1, . . . , N c and q = 1, . . . , N cs , Then we can express (4.1) and (4.2) as systems of ordinary differential equations whose time-dependent matrices are given by The time-dependent matrices for the interior bulk environment Ω h (t) which are given in (4.7) are defined similarly to the above with the appropriate substitutions for the interior region. The time-dependent boundary and load vectors are given by Similarly, we can express (4.3) and (4.4) as systems of ordinary differential equations and whose time-dependent matrices are given by and time-dependent interaction and load vectors are given by

Temporal discretisation
To obtain a temporal discretisation of (4.6), (4.10), (4.7) and (4.11), we subdivide the time domain [0, T ] into N T equal time intervals of size ∆t = T /N T and denote t n = n∆t, n = 0, . . . , N T . Note that the superscript n will be used to denote the time level t n . Following [37], we discretise the ALE mapping using linear interpolation between time levels: where A h (·, t) is the piecewise linear map at time t. The mesh velocity is therefore piecewise constant in time and is given by Following [37], the temporal discretisation of (4.6), (4.10), (4.7) and (4.11) is obtained using a modified Crank-Nicolson semi-implicit approach. First, for each j = 1, . . . , N ls and q = 1, . . . , N cs , the surface approximations (4.10) and (4.11) are predicted using an implicit-explicit Euler method: respectively, where the diffusion and advection terms are treated implicitly whilst the reaction terms are dealt with explicitly. For each i = 1, . . . , N l and p = 1, . . . , N c , the bulk approximations (4.6) and (4.7) are then updated using a Crank-Nicolson step and Note that in general (5.2) and (5.3) result in nonlinear systems. In order to guarantee that the surface approximations are second-order, we perform the correction steps 6 Mesh generation 6

.1 Surface domain
Following [37], we consider boundary movements given by the following evolution law for the normal velocity where α and β are given functions and κ is the curvature. In this article, the curve Γ(t) is represented by the position of a discrete set of nodal points on the curve. These mesh nodes evolve according to the equatioṅ where n and t are the outward unit normal and unit tangent vectors, respectively. The tangential motion B is necessary to maintain mesh quality. Following [42], the normal and tangential velocities can then be expressed in terms of the parameter- where F σ = ∂F/∂σ, |a| denotes the l 2 norm and τ is the mesh relaxation time. In this article, we do not consider adaptivity and therefore, the spatial balancing operator P = 1 and monitor function M = 1 but are included above for completeness. In the interests of brevity, we omit the derivation of (6.3) as they are covered in detail in [42]. To solve (6.2), we follow [37,42] by using second-order central finite differences to approximate the spatial derivatives and a first-order backward Euler scheme to approximate the temporal derivatives. As the mesh is evolved from time level t n to t n+1 , n = 0, . . . , N T , the surface is evolved first and then the nodal positions are used as Dirichlet boundary conditions for the bulk mesh movement which we discuss in the next section.

Bulk domain
In the interests of brevity, we discuss the motion of the interior region Ω h (t) only, as the procedure is identical for the extracellular environment with the minor addition, that after the surface Γ h (t) has evolved, the outer boundary is translated so that the interior region remains in the centre.
The evolution of the bulk mesh is assumed to satisfy a Moving Mesh Partial Differential Equation (MMPDE) [29,30,31,28]. The MMPDE is derived for the inverse ALE mapping A −1 (x, t) = ξ(x, t) as this prevents mesh crossings or foldings [11]. The mapping ξ(x) = (ξ(x), η(x)), corresponding to a fixed t, is chosen so that it minimises the functional [31] where G is a 2 × 2 monitor matrix. Following [37], we choose the monitor matrix to be where M is the monitor function. Rather than minimising the functional directly, the mapping is evolved according to the modifed gradient flow equations [31,37] ∂ξ ∂t where τ denotes the mesh relaxation time and P denotes the spatial balancing operator. To obtain an evolution equation for the physical mesh points {x i (t)} N Ω i=1 , we interchange the dependent and independent variables to give the MMPDE [37] τ p(x, t) ∂x ∂t where p(x, t) = 1/ P (x, t). The coefficients of the MMPDE are defined by [37] a = 1 M where J = x ξ y η − x η y ξ is the Jacobian of the mapping. To construct the initial bulk mesh, the equidistributed surface mesh is used as a set of fixed points for a Fortran implementation of the Matlab toolbox Distmesh [58]. The initial computational mesh is then defined as: Ω c,h = Ω h (0). As mentioned in §6.1, the nodal positions of the surface mesh evolution are used as Dirichlet boundary conditions for the bulk MMPDE (6.4).
For the numerical solution of (6.4), we discretise in space using piecewise linear finite elements and in time using a backward Euler method. To avoid solving nonlinear algebraic systems, the coefficients (6.5), and spatial balancing operator p, are treated explicitly. The discretisation of the MMPDE (6.4) therefore takes the form: find x n+1 for allv h ∈ L 1 0 (Ω c,h ) 2 . The resulting linear system is solved using the UMFPACK library [8] of the SuiteSparse software collection [9,7]. As the positions {x i (t)} N Ω i=1 are approximated using piecewise linear finite elements, the coefficients are in general piecewise constant and hence, not globally continuous. Therefore, within a standard Galerkin formulation, the coefficients require a numerical recovery procedure so that they can be differentiated accurately. Here, we use a spatial averaging technique; for example, let a : Ω c,h → R be a piecewise constant function, then one can define the function a by where a K = a| K is the local restriction of a to element K of the triangulation T Ω h,c . This procedure results in a ∈ C 0 (Ω c,h ) and therefore, we can approximate a using piecewise linear finite elements. The finite element approximation a h can then be differentiated in the usual way. Coincidentally, we also note that the spatial balancing operator may in general depend on these coefficients and is therefore calculated using the recovered coefficients. The increased regularity of the coefficients, due to the numerical recovery technique, produces a more robust moving mesh procedure and therefore, in general, will require remeshing techniques to be used far less frequently. s,h can be determined in terms of the column sums of the time-dependent mass matrices; that is, and where the time dependent mass matrices M Ω (t) and M s (t) are as defined in (4.8a) and (4.12a), respectively.
Proof. First, as a consequence of the finite element approximation, we note that the bulk and surface nodal basis functions form partitions of unity  s,h (t) over Ω h (t) and Γ h (t), respectively, and using the finite element approximation (4.5). Proof. Making use of the fact that both the gradient and tangential gradient operators are linear and that the finite element basis functions form a partition of unity (7.1), the result follows immediately: Proof. Once again, making use of the linearity of the gradient operator as well as the partition of unity of the finite element basis functions (7.1), we have for the bulk advection matrix Using the fact that the curve Γ h (t) is closed and the partition of unity of the surface nodal basis functions, the result for the surface advection matrix B s (t, w h ; u Γ ) follows similarly.
Similar to the continuous case §3.1, in order to demonstrate global conservation we require the following assumption: Assumption 7.1. We assume that on the boundary Γ h (t), the normal components of the ALE velocity and material velocity are identical; that is, w h · n h = u Γ · n h where n h denotes the outward unit normal to Γ h (t).
We are now ready to prove global conservation of the semi-discrete system given by (4.7) and (4.11).  are satisfied, then the semi-discrete system given by (4.7) and (4.11) is globally conservative independent of the ALE velocity; that is, given Lemma 7.1.1, the following holds true: where the time-dependent mass matrices are given by (4.8a) and (4.12a).
Proof. Following from Assumption 7.1, the matrix A Γ (t, w h ; u Γ , n Ω ) in (4.8c) is zero. We take the column sum of (4.7), apply Lemmas 7.1.2 and 7.1.3, and sum over p = 1, . . . , N c to give Similarly, taking the column sum of (4.11), applying Lemmas 7.1.2 and 7.1.3, and summing over q = 1, . . . , N cs yields The second term on the left hand side of (7.4) is given by where we have applied (7.1) and (7.2). Using similar arguments, it can be shown that where, once again, we have applied (7.1) and (7.2). The result then follows immediately from adding (7.4) and (7.5), then applying (7.6), (7.7) and (7.3).

Fully discrete
Following the semi-discrete case §7.1, in order to establish global conservation for the fully discrete system, we require the following assumption: Assumption 7.2. For any time level t n , n = 0, . . . , N T , we assume that on the boundary Γ h (t n ) = Γ n h , the normal components of the ALE velocity and material velocity are identical; that is, w n h · n n h = u n Γ · n n h where n n h denotes the outward unit normal to Γ n h . Observe that Lemma 7.1.2 and Lemma 7.1.3 hold for any fixed time t and therefore, hold at any time level t n , n = 0, . . . , N T . We are now in a position to prove global conservation for the fully discrete system (5.1), (5.3) and (5.5). where (7.8) has been used. Note that as a consequence of (7.9), we have, for any fixed time level t n , n = 0, . . . , N T , that p e T R (p),n (C n , C n s , L n s ) = q e T s R (q),n s (C n , C n s , L n s ). Adding (7.10) and (7.11), then applying (7.13) and (7.9) gives where, from (7.10) and (7.13d), the term in the square brackets in (7.15) is zero. Therefore, we can deduce that the predicted solution C s (t), ALE velocity and time step size ∆t play no role in global conservation. It is clear that (7.15) holds for any n = 0, . . . , N T and thus, the fully discrete system is globally conservative, as required.
Remark 2. As mentioned in Remark 1, the results presented in §8.1 consider interior bulk species only where we employ an implicit-explicit backward Euler scheme for the temporal discretisation. Therefore, Global conservation of (7.16) can then be obtained by following the steps in the proof of Theorem 7.2: apply Assumption 7.2; take column sums of (7.16); apply Lemmas 7.1.2 and 7.1.3; sum over p = 1, . . . , N c ; and, finally, apply (7.8).

Numerical Experiments
In the examples that follow, the conservation error is defined as the absolute difference between the total concentration at time levels t n+1 and t n ; that is, following Theorem 7.2, the conservation error is given by  The L 2 norm is then approximated using Gaussian quadrature , and the L ∞ norm is approximated by

Diffusion in a moving unit circle
In this example, we do not consider an exterior environment. Therefore, N l = 0 and N ls = 0. We consider diffusion of a single interior bulk species (N c = 1) and no surface species (N cs = 0) in a moving circular domain. We also assume that the bulk reaction term in (2.4a) is zero and employ zero-flux boundary conditions in (2.4b); that is, (1) ) = 0, andr (1) (c, c s , l s ) = 0. Following Novak and Slepchenko [56], we assume that the bulk domain Ω(t) is a unit circle moving with a constant velocity u b = (1, 0). We also assume that the boundary Γ(t) evolves with the same constant velocity u Γ = u b . However, the material velocity of the bulk domain u Ω = 0. With the diffusion coefficient D (1) c = 1, this problem has an exact solution given by where c 0 ∈ R. In [56], the constant value c 0 is calculated on the stationary box and is given by For the moving circular domain, a value of c 0 is not given in [56] and therefore, here we assume the same value for c 0 as for the stationary box. In the results which follow, we set the final time T = 0.5, the time step size ∆t = 5 × 10 −4 (so that N T = 1000) and initial condition c (1) (x, y, 0) = c 0 exp (−x). The initial mesh contains N Ω e = 1499 triangles and N Γ s = 87 boundary nodes. The boundary nodes are translated with the velocity u b and these are then employed as fixed Dirichlet boundary conditions for the MMPDE system with monitor matrix G = I, where I is the identity matrix. Figures 2(a) and 2 (b) show the computed solution at T = 0.5 and a comparison between the computed and exact solution along the cut y = 0, respectively. The computed solution is clearly in very good agreement with the exact solution. Figure 2 (c) confirms that the fully discrete numerical solution is globally conservative to machine zero. The rates of spatial convergence in the L 2 and L ∞ norms are shown in Fig. 2(d). Clearly second-order convergence is demonstrated in both norms. The rate of convergence is an improvement over the convergence rates obtained using the method proposed in [56], where the convergence rates were between 1 and 2. In [56], the impaired rates of convergence are attributed to a reinitialisation procedure of the unknowns at the boundary. Due to the fitted nature of the evolving mesh considered in this article, we do not require such a procedure and therefore, we obtain the expected rates of convergence.

Advection-diffusion in a moving unit circle
As an extension to the previous example, we consider the advection-diffusion of a chemical species in a moving circular domain. Once again, this problem has been considered by Novak and Slepchenko [56] and has an exact solution. We consider only a single interior bulk species (N c = 1), no surface species (N cs = 0) and no exterior environment so that N l = 0 and N ls = 0. We also assume that the bulk reaction term in (2.4a) is zero and employ zero-flux boundary conditions in (2.4b), as before.
Once again, we assume that the bulk domain Ω(t) is a unit circle moving with a constant velocity u b = (1, 0) and that the boundary Γ(t) evolves with the same constant velocity u Γ = u b . However, contrary to the previous example, the material velocity of the bulk domain u Ω = u b . The exact solution is then given by [56] c (1) where λ = 1.841183781340659, J 1 (z) is the first-order Bessel function of the first kind and |x| denotes the l 2 norm. In the results which follow, we set the diffusion coefficient D     Figure 3(d) illustrates the spatial convergence in the L 2 and L ∞ norms. Clearly second-order convergence can be seen in both norms. This is an improvement over the convergence rates shown in [56], where less than second-order convergence was seen for the L ∞ norm. Finally, Fig. 3 (c) confirms that the fully discrete numerical solution is globally conservative to machine zero.

Advection-reaction-diffusion in a moving complex geometry
In this section, we consider an example presented by Strychalski et al. [65,66] where the initial domain is given by We also assume that the material velocity of the bulk domain u Ω = u b . Following [66], the bulk reaction    in (2.4a) and flux boundary conditions in (2.4b) are given by and and where k 1 and k 2 are maximum activation and deactivation rates, respectively; and K m1 and K m2 are Michaelis-Menten constants. It is clear that (8.1) and (8.2) satisfy the conditions (3.9) of Theorem 3.1 and therefore, in light of Remarks 1 and 2, we expect this system to be globally conservative. In this example, c (1) and c (2) denote the concentrations of inactive and active protein conformations, respectively. Initially, the inactive and active concentrations are c (1) (x, 0) = 1 and c (2) (x, 0) = 0, respectively. Due to the complex geometry, this example does not possess a known exact solution. Following [66] we set the final time T = 0.4, the time step size ∆t = 1.25 × 10 −3 (so that N T = 320) and choose a uniform mesh such that where h K denotes the mesh width of element K of the triangulation T D h,c . Thus, the initial mesh contains N Ω e = 18659 triangles and N Γ s = 419 boundary nodes. Analagous to §8.1.1 and §8.1.2, the boundary nodes are translated with the velocity u b and then employed as fixed Dirichlet boundary conditions for the MMPDE system with monitor matrix G = I, where I is the identity matrix. For illustrative purposes, a coarser mesh is depicted in Fig. 4(a). We set the parameters D (2) c = 0.1, k 1 = k 2 = 1.0, K m1 = K m2 = 0.2 and S = 1.0, as in [66]. Figures 4(b)- (d) illustrate the concentration of the active species at times t = 0.15, 0.3, 0.4, respectively. Once again, the solid black line indicates the initial configuration. To plotting accuracy our results agree with those presented in [66]. Figure 4 (e) illustrates that the error in global conservation of the fully discrete numerical solution is machine zero. This is an improvement over the results presented in [66], where the conservation error was of the order 10 −5 to 10 −6 .

Application to directed cell migration
As discussed in §1, Rho GTPase activity is of key importance in directed cell migration and massconservative reaction-diffusion (McRD) equations are an important mathematical framework for their study. Therefore, unsurprisingly, the modelling of Rho GTPase activity has been the focus of many studies; see e.g. [63,47,32,57,24,61]. In resting cells the inactive form of the Rho GTPase is predominate. Some of the inactive form will be cytoplasmic, sequestered by GDI, and some will be membrane bound. The activation of the inactive form (in response to an upstream signal, for example) is facilitated by GDP/GTP exchange factors (GEFs), which stimulate the dissociation of GDP from the inactive form. This process cannot take place within the cytoplasm due to the inactive form being sequestered by GDI. Consequently, activation can only take place on the membrane and therefore, in order for sustainable activation to be seen, the dissociation of GDI is essential. Additionally, the exchange factors (GEFs) are cytoplasmic and therefore, must be recruited to the membrane so that activation can take place. Hydrolysis of GTP to GDP is stimulated by GTPase activating proteins (GAPs) and is responsible for the inactivation of the active form.
Marée et al. [47] considered a multiscale model of a motile cell, where the Rho GTPase activity was coupled to the cytoskeletal dynamics and motion of the cell using the Cellular Potts Model [23]. The GTPase activity demonstrated spatial bistability and, crucially, stable polarisation could be seen outside of the bistable region, meaning that a stable resting cell would be possible. The onset of polarisation was then studied in greater detail in [32], where cooperativity between the interacting Rho GTPases and the fast diffusion of the GDI bound inactive form were seen to be essential in developing robust polarity. The models proposed in [47,32] were analytically intractable due to their complexity. Therefore, building upon the previous work of [47,32], Mori et al. [49,50] introduced a minimalistic model for Rho GTPase activity by considering only a single Rho GTPase which cycles between an inactive (cytoplasmic) and active (membrane bound) state. The model, also known as the wave-pinning (WP) model, consists of two time-dependent reaction-diffusion equations in one spatial dimension and is given by where a, b denote the active and inactive states, respectively and L ∈ R. The reaction kinetics are given by where k 0 , γ, K, δ ∈ R are model parameters and f s is an external stimulus. The external stimulus is assumed to increase the rate of activation. Under the influence of an external stimulus, this model sets up a travelling wave in the active state which is later halted (pinned) due to the interaction with the inactive state. Fundamental assumptions at the heart of this model are the properties of mass conservation and also the large differences in the active and inactive diffusivities. Vanderlei et al. [70] later extended this model to two spatial dimensions, incorporated a mechanical model for the membrane and simulated its movement using the immersed boundary method [59]. Note that following the discussion above, a consequence of this model is that the dissociation of GDI and the action of GEFs to produce the active form, are incorporated into a single step. Also note that the GEFs are assumed to be freely available.
In [47,32,49,70], the active and inactive forms are intermixed and therefore, occupy the same physical domain. One could arrive at this scenario by considering a top-down view of a cell, essentially flattening a  three-dimensional cell into two dimensions. However, as discussed above, the active and inactive (without GDI) forms are membrane bound and therefore, occupy a physically distinct region to the inactive GDI bound state. This compartmentalisation is a key biological component of cell polarisation. Bulk-surface implementations of the wave-pinning model can be found in the literature (see e.g. [60,21,10,6]), where the domain is assumed to be stationary. Although the original formulation of the wave-pinning model has been applied on moving domains (see e.g. [47,70,27]), we are not currently aware of any studies involving the bulk-surface wave-pinning model on moving domains. Section 8.2.1 illustrates the mass-conservative properties of the bulk-surface wave-pinning model on a stationary circular domain. Section 8.2.2 then considers the bulk-surface wave-pinning model on a moving domain, where the motion of the cell membrane is driven by the level of activation and the stimulation is dependent on an extracellular ligand field and membrane bound ligand-receptor concentration.

Bulk-surface wave-pinning on a stationary circular domain
Before we consider the bulk-surface wave-pinning model on a moving domain, we illustrate the mass conserving properties of our algorithm on a stationary, circular domain [6]. Similar to the previous examples, we assume that there is no exterior environment so that N l = 0 and N ls = 0. The initial domain Ω(0) is assumed to be a circle of radius R 0 centred at the origin. As the domain is stationary, Ω(t) = Ω(0), ∀t ∈ I, and the material velocities u Ω = u Γ = 0. In this example, there is a single bulk species (N c = 1) corresponding to the inactive state and a single surface species (N cs = 1) corresponding to the active state. Therefore, in (2.4) f  cs (c s ) = 0. The interaction between the bulk and surface species is given by the reaction-kineticŝ where ω ∈ R is a length scaling parameter. It is clear thatr (1) andr (1) s given above satisfy the conditions (3.10) of Theorem 3.1 and therefore, we expect the system to be globally conservative. Note that in this example, we do not consider an external stimulus; that is, f s = 0. Following [6,62] we assume that ω is defined as the ratio between the bulk volume/area and surface area/length; that is, ω = |Ω(0)|/|Γ(0)|. Following [6], the initial conditions are given by where the (dimensional) parameters are detailed in Table 1 and θ denotes the angle measured anticlockwise from the positive x-axis. The parameters are the same as used in [49,6] and the concentrations s l , s h and c (1) 0 are chosen so that wave-pinning occurs [49,6]. We set the final time T = 6000 and time step size ∆t = 0.1 (so that N T = 60000). The initial mesh is uniform and contains N Ω e = 2097 triangles and N Γ s = 103 boundary nodes.
Parameter Value Units Parameter Value Units Parameter Value Units  The concentration of the active surface state as a function of the scaled arc-length around the circular cell and the inactive cytoplasmic state are illustrated at times t = 5, 45, 105, 145 in Fig. 5. Clearly when t = 5, the surface active state displays a sharp Gaussian-like shape centred around θ = 0 (top image of Fig. 5(a)) which is inherited from the initial condition. Correspondingly, the cytoplasmic inactive  state locally depletes around the highest membane-bound active state concentration (bottom image of Fig. 5(a)). The peak in the membrane bound activated state flattens, and therefore broadens, as time progresses (see e.g. top image of Fig. 5(b)) due to the wave of activation travelling around the membrane. Correspondingly, two locally depleted regions just inside the active wave-front can be seen in the cytoplasmic inactive state (see e.g. bottom image of Fig. 5(b)). This process continues to later times (see e.g. Fig. 5(d)) until eventually the wave is pinned, the process terminates and (assuming no external stimulus) a stable pattern is formed. The results depicted in Fig. 5 are in excellent qualitative agreement with those presented in [6]. Note that in Figs. 5(a)-(d) the colourbar axis changes as time progresses. Due to the high diffusivity inside the cell, the inactive cytoplasmic species is approximately constant and therefore, changing the colourbar axis as time progresses emphasizes the local pattern of the cytoplasmic species. Figure 6(a) confirms that the fully discrete numerical solution is globally conservative to machine-zero. The drop off in conservation error as time progresses is most likely caused by the travelling wave of activation halting as it becomes pinned.

Bulk-surface wave-pinning model on a moving domain with extracellular ligand field
In the exterior environment Ξ(t), we assume there exists an external guidance cue in the form of freely diffusing chemoattractant ligand molecules. At the cell membrane Γ(t), these molecules bind reversibly to membrane bound receptors. These transmembrane receptors have binding domains facing both the exterior and interior environments, making them an essential component in the response of a cell to an external guidance cue. Upon binding with an extracellular ligand, the receptor undergoes a conformational change resulting in an active ligand-receptor complex. This active complex then triggers downstream signalling events which lead to cellular migration. Although important for the activation of intracellular signalling pathways, the action of the active ligand-receptor complex on downstream signalling components is that of a catalyst and therefore, is not used up in subsequent reactions. Although the intracellular signalling pathways have received a large amount of attention in the literature, the coupling of these pathways to an external guidance cue, on a moving domain, has received far less attention. Usually, the external guidance cue (or equivalently, ligand-receptor complex concentration) is prescribed in the activation of membrane bound species (see e.g. [48,36,46,16]) or the interaction is studied on a stationary domain (see e.g. [35,68,19]). Recently, Moure and Gomez [53] proposed a phase-field model for the simulation of obstacle-mediated chemotaxis. The cytoplasmic mechanics of the acto-myosin network were considered to be dependent on a membrane bound activator, which in turn depends on the membrane concentration of chemoattractant. This suggests that the activator approximates the entire signalling pathway from extracellular ligand concentration to the intracellular acto-myosin network. Feng  et al. [18] proposed a lattice Boltzmann-particle method for neutrophil migration, where their model comprised four signalling layers: extracellular chemoattractant and membrane bound receptor; membrane bound G protein dynamics and specific effectors; cytoplasmic and membrane bound Rho family of GTPases; and their downstream effectors. However, the membrane is unfitted to the underlying computational lattice and therefore, an additional lattice-particle map was required to redefine the subset of the computational lattice (which defines the membrane) as the cell migrates. Within the context of conservation, this procedure would require care to ensure information is not lost under motion. We consider a single exterior bulk species (N l = 1) corresponding to the extracellular chemoattractant ligand and a single surface species N ls = 1, which corresponds to the active ligand-receptor complex. We assume that the material velocity of the extracellular ligand field is given by u Ξ = 0. Following [37], we employ Dirichlet boundary conditions on the outer boundary ∂D(t). The exterior bulk reaction term, and boundary conditions for the outer boundary, will be defined shortly. On the inner boundary (cell membrane) Γ(t), the bulk boundary conditions (2.2c) and surface reaction kinetics (2.3) is given by ls (l s ) = 0, andĝ (1) (l, l s ) =ĝ (1) s (l, l s ) = k −1 l (1) respectively, where k 1 is the ligand association rate, k −1 is the disassociation rate and R T is the total concentration of receptors. Following the previous section §8.2.1, in the interior environment Ω(t) we assume there is a single bulk species (N c = 1) and a single surface species (N cs = 1). Furthermore, we assume: that in (2.4) f cs (c s ) = 0. However, contrary to the previous section, the interaction between the interior bulk and surface species are given by the reaction-kineticŝ where ω ∈ R is a length scaling parameter defined in §8.2.1, ε is the rate of stimulation andR = l (1) s /R T defines the local fractional receptor occupancy. We assume (initially) that the intracellular material velocity u Ω = 0. However, due to the dynamic motion of the cell membrane considered in this example, u Γ = 0. In fact, in general the material velocity u Γ is not known a-priori. As the shape of the curve is determined solely by the normal motion, we assume u Γ · t = 0.
The (dimensional) parameters for this model are detailed in Table 1 and Table 2. The parameters are the same as used in [49,6,37,41], with the exception of ε and σ l which are guessed. Furthermore, so that the wave-pinning model has enough time to establish a polarised state, the diffusion parameters D (1) ls and D (1) l have been lowered by an order of magnitude (but maintain the order of magnitude of their ratio) compared with those used in [37,41]. Initially, the extracellular environment Ξ(0) is assumed to be an annulus where the inner boundary Γ(0) is a circle of radius R 0 and the outer boundary ∂D(0) is also a circle of radius R outer . Therefore, initially the intracellular environment Ω(0) is a circle of radius R 0 .
The radius of the outer boundary needs to be sufficiently large so that the Dirichlet boundary condition remains valid. Assuming that there are N R /(4πR 2 0 ) receptors per unit surface area of a spherical cell (number concentration), then the molar concentration of the total number of receptors is calculated using the following formula: where N A = 6.022 × 10 23 is Avogadro's constant and the factor 10 24 arises from the conversion to molar concentration with nanomolar units (nM) and to a length scale with micron units (µm). For the extracellular ligand and membrane bound ligand-receptor complex, we initially prescribe a homogeneous fractional receptor occupancy level (R 0 ) which is given in Table 2. The initial condition for the extracellular ligand field and the membrane bound ligand-receptor complex is then calculated from the initial fractional receptor occupancy level and the steady state [55]:  Table 1. The membrane Γ(t) is evolved according to the procedure summarised in §6.1 where the normal velocity (6.1) is given by [55,37] V(x, t) = K prot c (1) s where K prot is the strength of the protrusive force which is proportional to the concentration of the membrane bound active species and λ(t) is a time-dependent cortical tension factor found according to [55,37]: where υ, β ∈ R + are positive parameters given in Table 2, s,0 and A 0 is the initial area contained in the cell. As a full investigation of parameter space is beyond the scope of this article, K prot , υ and β are chosen such that the cellular migration remains stable. We set the final time T = 24000 and time step size ∆t = 0.1 (so that N T = 240000). The interior and exterior meshes contain N Ω e = 2091 and N Ξ e = 5123 triangles, respectively. The boundary of the interior mesh (and inner boundary of the exterior mesh) contains N Γ s = 101 nodes, whilst the outer boundary of the exterior mesh contains N ∂D s = 68 nodes. For computational efficiency, we employ a graded mesh in the exterior environment, where the elements become larger the further they are from the inner boundary. However, in the interior environment, we employ a uniform triangular mesh. Although the membrane shape is dynamically changing, we fix the outer boundary of the extracellular environment to be a circle which is merely translated so that the cell remains in the centre. The extracellular and intracellular meshes are then evolved according to the procedure outlined in §6.2.
The concentration of the interior active surface species as a function of scaled arc-length around the cell and the intracellular inactive species are illustrated at times t = 25, 105, 225, 300 in Fig. 7. In Parameter Value

Units
Parameter Value Units Parameter Value Units   contrast to the stationary example, at t = 25 the active state displays a broad peak on the side of the cell which is facing the point source and a broad, depleted region in the cytoplasm ( Fig. 7(a)). This broad activation becomes sharper at t = 105 due to the expansion of the cell in the direction of the point source ( Fig. 7(b)) which in turn, produces a localised depleted region in the cytoplasm analagous to the stationary example (Fig. 5a). Due to the sharp peak in the active state, the cell protrudes local to this region trapping the depleted cytoplasmic species thereby creating a cap of cytoplasmic inactivity local to the protrusion (Figs. 7(c) and 7 (d)). Analagous to the stationary case, the peak in the active surface state flattens as time progresses due to the wave of activation travelling around the membrane, eventually pinning to produce a stable pattern. Figure 8 illustrates the extracellular ligand concentration and intracellular inactive concentration at later times than those considered in Fig. 7. In Fig. 8(a) it is clear that the cell is migrating towards the point source located at x p = (15, 0). There is a clear cap of inactivity in the intracellular environment in the direction of the point source, due to the cell protruding in that direction as was seen previously. Once the cell gets sufficiently close to the point source, the point source is moved and the active surface state migrates around the membrane so that it faces the point source. This causes the cytoplasmic cap of inactivity to also migrate so that it is facing the point source and this can be seen clearest in Figs. 8 (b) and 8(d). This type of migration is in contrast with the Meinhardt model [48,55,37], where one expects a new peak of activation to form in the direction of the point source. However, in the wave-pinning model, the activation remains but migrates round the membrane to face the point source (as demonstrated previously in the absence of a freely diffusing external guidance cue [47,46]) which is a testament to the stability of the pattern formed by the wave-pinning mechanism. After turning to face the point source, the cell straightens and migrates persistently towards the point source where the point source is again moved when the cell gets sufficiently near (Figs. 8(c)-(f)). Analagous to the previous section §8.2.1, the colourbar axes in Figs. 7 and 8 change as time progresses to emphasize the local cytoplasmic patterning as the cell migrates. Note that due to the robustness of the MMPDE procedure discussed in §6.2, our algorithm is able to perform the full loop without the need to re-mesh. An illustration of the exterior and interior meshes during migration is given in Fig. 9(a). Figure 6 (b) confirms that the fully discrete numerical solution is globally conservative to machine zero as proven in Theorem 7.2.
For the results depicted in Figs. 7 and 8, we assumed that u Γ · t = 0. This assumption is reasonable as the shape of an evolving curve is solely determined by its normal motion. However, in general u Γ is not known a-priori and therefore, u Γ · t = 0 is entirely possible. In the context of cell migration, the material velocity can be approximated using a mathematical model of the intracellular fluid dynamics and the acto-myosin network; see for example, [1,2,51,53]. However, such an implementation would warrant significant changes to the algorithms presented in this article and is therefore, not considered here. Instead, we assume that material velocities for the intracellular and membrane environments are the same as the ALE velocity; that is, u Γ = w and u Ω = w. Hence, we have assumed that the inside of the cell behaves like an elastic solid. To emphasize the importance of the material velocity, we sharpen the Gaussian point source (8.4) by reducing the value of σ l = 0.1. A comparison of the migration of the cell when the material velocities are given by: u Γ · t = 0 and u Ω = 0; and u Γ = w and u Ω = w, is illustrated in Fig. 9 (b). Both cells successfully migrate towards the initial location of the ligand point source x p = (15, 0) (which is depicted in Fig. 9(b) as a yellow circle). Their migration is slower than depicted in Fig. 8 due to the sharper Gaussian point source. Once close enough, the point source location is moved to x p = (15, 15) (which is illustrated in Fig. 9(b) as the extracellular contour plot). The cell with material velocities u Γ · t = 0 and u Ω = 0, cannot turn sharply enough to be able to migrate towards the new point source location (right-most cell in Fig. 9(b)). Therefore, the cell migrates past the location of the point source. The solid black line in Fig. 9(b) depicts the evolution of the cell centroid in time. However, the cell with material velocities u Γ = w and u Ω = w (left-most cell in Fig. 9(b)) can turn sharply and in fact, gets sufficiently close to the point source that it moves to a new location x p = (0, 15) (depicted as a yellow square in Fig. 9(b)). The black dashed line in Fig. 9(b) depicts the evolution of the cell centroid in time. The material velocity u Ω = w ensures that the inactive cytoplasmic species is transported with the cell as it migrates. This produces a more uniform distribution on the inside of the cell (this can be seen in Fig. 9(b)). Consequently, the inactive species is more readily available to be converted to membrane bound active species, which in turn enables the cell to adapt more efficiently. The depleted regions found in the cytoplasm local to a protrusion (see Figs. 7 and 8) can therefore act as a local inhibitor to migration.

Conclusions
Bulk-surface mass-conservative reaction-diffusion systems are a promising framework for modelling GT-Pase cycling leading to cell polarisation and subsequently cell migration and chemotaxis. We have established theoretically that a recently developed finite element ALE scheme is globally conservative at the fully discrete level and this property holds independently of the chosen time step and of the domain movement. To our knowledge this is the first scheme to be shown to have this desirable property. Furthermore, we have shown that the numerical scheme appears to be second-order accurate when applied to test cases with known exact solutions. The method has been applied to a bulk-surface WP model for cell polarisation in combination with a receptor-ligand model to describe a signalling cascade from external stimulus to downstream cell protrusive activity. The modelled cell has been shown to display efficient chemotaxis even when the source of chemoattractant changes location. Interestingly, numerical experiments suggest that the accuracy of a model cell to direct its motion can be affected by the bulk material velocity responsible for the transport of cytoplasmic proteins.
The main theoretical result of this paper concerns the conservation property of the fully discrete ALE finite element scheme. Future work will focus on the numerical stability and convergence of the scheme. Previous studies on the stability of discretisations of reaction-diffusion equations on evolving and growing domains suggest that discretisations of conservative formulations of the governing equations lead to schemes with superior stability characteristics [39,38,40]. While we have used a relatively simple model for ligand-receptor binding to demonstrate the possibilities of the computational framework, the same techniques could be used to investigate the effect of local extracellular ligand depletion by membrane bound enzymes and intracelluar signalling pathways [41]. The method could also be used to consider the effect of cell shape changes on intracellular signalling from the cell cytoplasm to the cell nucleus [22,5,67]. In future work we also plan to couple the GTPase module to a mechanical model for the bulk material velocity which will be driven by the intracellular forces generated by the interaction of actin  and myosin and substrate adhesion. Models of this type range in sophistication from purely viscous flow models [64,20] to models which take into consideration the orientation of actin filaments [33,34,45].