Mechanochemical Models for Calcium Waves in Embryonic Epithelia

In embryogenesis, epithelial cells acting as individual entities or as coordinated aggregates in a tissue, exhibit strong coupling between mechanical responses to internally or externally applied stresses and chemical signalling. One of the most important chemical signals in this process is calcium. This mechanochemical coupling and intercellular communication drive the coordination of morphogenetic movements which are characterised by drastic changes in the concentration of calcium in the tissue. In this paper we extend the recent mechanochemical model in Kaouri et al. (J. Math. Biol. 78, 2059–2092, 2019), for an epithelial continuum in one dimension, to a more realistic multi-dimensional case. The resulting parametrised governing equations consist of an advection-diffusion-reaction system for calcium signalling coupled with active-stress linear viscoelasticity and equipped with pure Neumann boundary conditions. We implement a finite element method in perturbed saddle-point form for the simulation of this complex multiphysics problem. Special care is taken in the treatment of the stress-free boundary conditions for the viscoelasticity in order to eliminate rigid motions from the space of admissible displacements. The stability and solvability of the continuous weak formulation is shown using fixed-point theory. Guided by the bifurcation analysis of the one-dimensional model, we analyse the behaviour of the system as two bifurcation parameters vary: the level of IP3 concentration and the strength of the mechanochemical coupling. We identify the parameter regions giving rise to solitary waves and periodic wavetrains of calcium. Furthermore, we demonstrate the nucleation of calcium sparks into synchronous calcium waves coupled with deformation. This model can be employed to gain insights into recent experimental observations in the context of embryogenesis, but also in other biological systems such as cancer cells, wound healing, keratinocytes, or white blood cells.

model. We also demonstrate the nucleation of calcium sparks into syn-

Introduction
Embryogenesis is a remarkable example of a complex process where different sub-mechanisms involving mechanical and chemical effects closely interact in a self-organised manner, forming complex spatio-temporal patterns. The coupling between rearrangement of tissue, cell migration, active cell contraction to diffusion of morphogens or signalling molecules has been proposed and studied in earlier works [40,41]. It is well known that a variety of responses in the cell are driven by the transduction of mechanical stimulation into chemical signals such as calcium oscillations and waves [54]. In turn, localisation of stresses or strains within the cells can generate patterns of motion distribution in tissue by changing their displacement magnitude, direction and velocity [23,33]. Recent experimental evidence [17,42,55] shows that tissue mechanics are actively coupled to chemical patterns during development.
Hence, more modelling and analysis has recently appeared in order to elucidate the many open questions that directly impact the healthy evolution of an embryo [13,28,42,55]. A variety of models have described general interactions between chemical species concentration and mechanics; these include descriptions where stress is triggered by chemical signalling [41], or mainly by migration [38], but we are also interested more generally in the mechanochemical feedback due to pressure and velocity of the underlying embryonic tissue. In all cases, including the coupling with continuum mechanics significantly affects the propagation of the chemical signals.
Here, we develop a multi-dimensional extension of the recent mechanochemical model proposed in [28] which describes the interplay of calcium signalling with the mechanics of embryonic epithelial tissue during apical constriction, an active deformation process. Inspired by the recent, interesting experimental observations in [17,55] where increasing tension in the contracting cells yields further calcium release and this, in turn, elicits contractions which are sensed as mechanical stimuli by the neighbouring cells. The model we propose here takes a step further in elucidating this important mechanochemical coupling.
Mechanical properties of different cell types indicate diverse behaviour, including elastic [15,18,45], poroelastic [14,37,50,51], or nonlinear and nonlocal characteristics [36,56] but more predominantly, viscoelastic effects [2,8,10,20,29,47,57]. The specific constitutive rheological model to adopt in a tissue depends on the characteristics of each constituent cell, on the properties inherent to distinct biological states, on the nature and intensity of the stresses and strains that are to be applied, and on the spatio-temporal scales involved. Here, we restrict the description to the regime of small strains and model the cell and tissue as a modification of the simple Kelvin-Voigt viscoelastic solid (with one elastic spring and two viscous dashpots), where only after the initial stress has vanished, the material goes back to its original state. These linear viscoelastic materials are completely defined by the stiffness and viscosity, which can be determined using diverse measuring approaches such as pipette suction, optical laser tweezers, microrheology tools, particle tracking, or even contact-free techniques [44]. In the present mechanochemical model, we assume that the viscoelastic stress has an active tension component which depends on the concentration of calcium, following the formulation in [6,28,41].
We adopt the following fundamental assumptions: a) the equilibrium of forces in the system is established by a quasi-static balance of linear momentum using displacements and hydrostatic pressure (the so-called Herrmann formulation [26] where the introduction of solid pressure ensures that the system is robust with respect to the modulus of dilation of the solid); b) the spatio-temporal dynamics of calcium concentration and the percentage of IP 3 receptors (IPR) on the Endoplasmic Reticulum (ER) that have not been inactivated are governed by an advection-reaction-diffusion system; and c) mechanochemical coupling is modelled directly in the viscoelastic stress through an active stress Hill function that depends on calcium and through the modification of the reaction kinetics by volume change. The two-way coupling mechanism we adopt follows the model structure used in [28,39,40,43,53].
Finding closed-form solutions to this inherently highly nonlinear and multidimensional problem is only possible in very restricted scenarios and simplified settings. We, hence, resort to solving the governing equations numerically. The numerical framework undertaken here uses the method of lines, adopting a backward Euler scheme for the discretisation in time and a primal-mixed formulation consisting of mixed approximations for the viscoelasticity in terms of displacement and pressure using the classical Taylor-Hood and MINI-elements [4], and piecewise quadratic or piecewise linear approximations for the calcium concentration and the fraction of non-inactivated IPR. Methods of this type are known to perform well in a variety of scenarios. As we assume that the cell or tissue is not attached to any supporting structure, we consider pure traction boundary conditions. In this case the viscoelasticity problem is not well-posed unless we incorporate a constraint to eliminate the rigid motions from the set of admissible solutions. To achieve this we employ additional vector Lagrange multipliers for the coupled problem following [32] (see also [7]). The overall problem is treated as a monolithic system, so at each time iteration we solve a set of nonlinear equations with the Newton-Raphson method and the tangent system at each step is inverted with a direct solver.
We have organised the contents of the paper as follows. In Section 2 we lay out the details of the mechanochemical 1D model in [28] and construct its multidimensional extension. We also explain the coupling mechanisms and then perform an appropriate nondimensionalisation. In the same section we also state the weak form of the governing equations and introduce a suitable finite element discretisation. In Section 3 we address the continuous dependence on data of the weak formulation, as well as the existence of weak solutions using Brouwer's fixed-point theory. A number of illustrative numerical computations are then presented in Section 4, where we also perform a simple verification of convergence. We specifically explore different regimes of wave propagation of practical interest, such as solitary waves and periodic wavetrains of calcium and the the nucleation of calcium sparks into calcium waves. We conclude in Section 5 with a summary of our findings and a discussion on the limitations of the model, addressing also possible extensions and future directions.

Model description and weak formulations 2.1 Coupling calcium signalling with mechanics
We assume that the cell/tissue can be macroscopically regarded as a linear, viscoelastic material of Kelvin-Voigt type, occupying the spatial domain Ω ⊂ R d with d = 2 or d = 3.
Following, e.g., [28,38], and assuming that gravitational forces and inertial effects are negligible, one seeks for each time t ∈ (0, t final ], the displacements of the tissue, u(t) : Ω → R d , and the dilation θ(t) : Ω → R d , such that where σ is the Cauchy stress tensor, ε(u) = 1 2 (∇u + ∇u T ) is the tensor of infinitesimal strains, and E, and ν are the Young modulus and Poisson ratio associated with the constitutive law of the material, respectively. Equation (2.1c) represents the force equilibrium for the system in the absence of inertia, whereas both (2.1a) and (2.1b) are constitutive equations describing properties of the viscoelastic material. The parametersα i in the constitutive relation (2.1b) are the shear viscosity and the bulk viscosity related to the total Cauchy stress exerted in the cell/tissue, that characterises the viscoelastic response to deformations. In addition, the last term in (2.1b), T (c) describes the active contraction stress which is dependent on the calcium concentration, c. There are various ways to model this term. As in [28] we assume a Hill function which saturates to a constant value T 0 for high values of c, in line with experimental observations from [17]. Therefore, where n is a positive integer and κ > 0. We, thus, assume that the calcium concentration controls entirely the material's motion by regulating the active tension and therefore the generation of stress. Such control is indeed an activation and not an inhibition of active tension, so T 0 is assumed positive. As in [38], we do not incorporate external body loads or restoring displacement-dependent body forces on the right-hand side of (2.1c), since such terms are only relevant to substrate-on-substrate or in tissue-on-substrate configurations.
Focusing on the spatiotemporal behaviour of calcium, we denote its concentration by c and h represents the percentage of non-inactivated IPR on the ER. The model is a generalisation of the recent model in [28], which employs the calcium dynamics of the well-known model in [5]. In dimensional form of the model is written as follows Here, div u = θ represents the dilation/compression of the apical surface area of the cell and D is the calcium diffusion coefficient inherent to the tissue (or to the cytosol in the single-cell case), which is assumed positive and constant. The fluxes in (2.3a) are as follows: the term J ER models the flux of calcium from the ER into the cytosol through the IPR, µ(p 3 ) = p 3 /(k µ + p 3 ) is the fraction of IPR that have bound IP 3 and is an increasing function of p 3 , the IP 3 concentration. The constant k f denotes the calcium flux when all IP 3 receptors are open and activated, and b represents a basal current through the IPR; J pump represents the calcium flux pumped out of the cytosol where γ is the maximum rate of pumping of calcium from the cytosol and k γ is the calcium concentration at which the rate of pumping from the cytosol is at half-maximum. One could also include calcium fluxes leaking into the cytosol from outside the cell, but we leave those terms out, as they are assumed small. The term J SSCC encodes the calcium flux due to the activated stretch-sensitive calcium channels (SSCC). These channels have been identified experimentally in recent years [24]-they are on the cell membrane. SSC are activated when exposed to mechanical stimulation and they close either by relaxation of the mechanical force or by adaptation to the mechanical force [24]. The constant S represents the strength of stretch activation. (In [28] the authors derive a relationship for S as a function of the characteristics of a SSCC, under suitable assumptions.) The inactivation of the IPR by calcium acts on a slower timescale compared to activation [5] and so the time constant for the dynamics of h, τ j > 1 in the ODE (2.3b) governing the dynamics of h. As in [28] and [5] we take τ j = 2 s. The values of all other calcium signalling parameters are also taken as in [5].
We nondimensionalise the set of governing equations (2.1a)-(2.1c) and (2.3a)-(2.3b) using where L is the length of an embryonic epithelial tissue (or the maximal length of the cell if considering the single-cell case, see [35]). In addition, instead of the dilation θ, we use the rescaled pressure p, which leads to the following system, where we have dropped the bars in the dimensionless variables for notational convenience For the parameter values we have chosen from [1,5] and also taking D = 20µm 2 /s and L = 100 µm we obtain the following values for the nondimensional parameters D = 4 × 10 −3 , K 1 = 46.29, G = 40/7, K = 1/7. In this context, a large value of K 1 captures the fact that calcium is a fast messenger [9]. We have assumed here that mechanics modify the behaviour of calcium only through the advection term and an additional calcium flux which is linearly dependent on the local dilation. The latter flux is modulated by λ > 0, thus representing a source for c if the solid volume increases, and a calcium sink otherwise (see e.g. [43]). This parameter will be treated as a bifurcation constant, as in [28]. In [28] it has been shown that for a critical value of λ, oscillations were suppressed, and we expect a similar behaviour in the higher dimensional case here. We will vary λ from 0 (calcium dynamics not affected by mechanics) to λ = 5. However it is not clear what the maximum λ value should be -inspired by the results in [28] we expect the oscillations/waves to vanish for large enough λ. Note that α 1 , α 2 determine the magnitude of the viscous effects. Also, the smaller β 2 is, the faster T (c) in (2.2) (actually its nondimensional form) will tend to its saturation value β 1 . In [28] the authors have explored this variation of the Hill function T (c); here we take the values α 1 = 1, α 2 = 0.5, β 2 = 0.1.
The system composed by equations (2.4a)-(2.4b) and (2.4c)-(2.4d) is complemented with appropriate initial data at rest Homogeneous boundary conditions can be incorporated for normal displacements, calcium fluxes, and traction, in the following manner where the boundary ∂Ω = Γ ∪ Σ is disjointly split into Γ and Σ where we prescribe slip boundaries and zero traction, respectively. This case assumes that the tissue is allowed to slip along the substrate on Γ, while it is free to deform on Σ. However, in most of our numerical tests we will consider instead the pure traction boundary conditions In this case an additional condition is required to make the system well-defined. For instance, we can impose orthogonality to the space of rigid motions defined as (see [12, Eq. (11.1.7)])

RM(Ω)
and unique solvability (for a given calcium concentration c) will follow since RM(Ω) is a null space of the relevant functional space.

Mixed-primal weak formulation
Multiplying the nondimensional governing equations (2.4) by adequate test functions and integrating by parts (in space) whenever appropriate, we end up with the following variational problem, here restricted to the case of pure traction boundary conditions: For a given t > 0, find (2.9d) Here we have defined the additional space meaning that all rigid motions are discarded as feasible displacement solutions. Alternatively, one can weakly enforce orthogonality to the space of rigid motions by a Lagrange multiplier, which is what we will do at the discrete level. Furthermore, we use an equivalent skew-symmetric form for advection terms and define β(c) = β1c n β2+c n , K(h, c) := µhK 1 b+c 1+c − Gc K+c and J (c) = 1 1+c 2 .

Fully discrete form
Let us consider a shape-regular partition T j ofΩ into affine elements (triangles in 2D or tetrahedra in 3D) E of diameter j E , where j = max{j E : E ∈ T j } denotes the meshsize. Finite-dimensional subspaces for the approximation of displacement, pressure, and calcium are specified as where P k (E) denotes the space of polynomials of degree less than or equal than k defined locally over the element E ∈ T j , and b E := ϕ 1 ϕ 2 ϕ 3 is a polynomial bubble function in E, and ϕ 1 , ϕ 2 , ϕ 3 are the barycentric coordinates of E.
Depending on whether the displacement is approximated with V j or V j , the pairs V j × Q j are known as the MINI element and V j × Q j are known as the Taylor-Hood elements. They are both Stokes inf-sup stable (see, e.g., [11,49]). On the other hand, an L 2 -orthonormal basis for the space of rigid motions is constructed in terms of translations and rotations associated with the principal axes of the inertial tensor I encoding rotational kinetic energy [32]. Denoting x 0 = |Ω| −1 Ω x the centre of mass of the domain, such basis (having dimension 6 in 3D) assumes the form where the (λ i , ν i ) are the eigenpairs of I (see also [31]). As remarked in [32], we emphasise that not all Stokes inf-sup stable pairs will lead to j−robust bounds when approximating the singular Neumann problem of linear elasticity (for instance, using for displacement and pressure the pair P 2 − P 0 does not allow for constructing a well-posed mixed Poisson auxiliary problem needed in establishing orthogonality with respect to the kernel).
Next we discretise the time interval (0, t final ] into equi-spaced points t k = k∆t. Applying a backward Euler method and an implicit centred difference discretisation of the first-and second-order time derivatives, respectively; we can write a semidiscrete form of (2.9) (but now incorporating the orthogonality with respect to rigid motions with a Lagrange multiplier). Now the fully discrete formulation is given by (2.11d)

Well-posedness analysis
We will consider that the initial data are nonnegative and regular enough, hence the numerical scheme starts from discrete initial data u 0 j , p 0 j , c 0 j , h 0 j , which are the projections of the exact initial conditions onto the finite element spaces.
Let us introduce the trilinear form d : Then, due to the skew-symmetric form of the operator d(·; ·, ·), we can write And recalling the Poincaré-Friedrichs inequality [22, Chapter I, Lemma 3.1] , we have the coercivity for e(·, ·): We also assume the nonlinearities satisfy the following conditions: is uniformly bounded with respect to c and Lipschitz with respect to h, in particular we have: Moreover, we will use the following algebraic relation: Let a and b be two real numbers, then we have

Continuous dependence on data
where C 1 , C 2 and C 3 are positive constants independent of j and ∆t.
Proof On the second equation of problem (2.11) we choose q = p k+1 j , to get which after applying Young's inequality, implies Now, on the first equation in problem (2.11), we choose the test function v = u k+1 multiply by 2∆t and use algebraic relation (3.3) in combination with equations (3.4), (3.5). We obtain: By Korn's and Young's inequalities, we get Summing over k from 0 to n ≤ N −1, applying the Poincaré-Friedrichs inequality and assuming 1 −

Existence result: Fixed-point approach
Next, we address the unique solvability of problem (2.11). To that end, we will make use of Brouwer's fixed-point theorem in the following form (given by [   Proof To simplify the proof we introduce the constants: We proceed by induction on k > 2. We define the mapping using the relation Note that this map is well-defined and continuous on  Next, using Korn's Inequality (with constantĈ k ) and (3.5), we deduce that Then, setting we may apply the inequality a + b ≤ √ 2(a 2 + b 2 ) 1/2 , valid for all a, b ∈ R, to obtain

Linearisation
At each time iteration we are left with a nonlinear system, and proceeding with a Newton linearisation we finally obtain the following scheme: Starting from discrete initial data u 0 j , p 0 j , c 0 j , h 0 j (projections of the exact initial conditions onto the finite element spaces), and for = 1, . . ., we find u +1 as the converged solutions of the iteration for k = 0, . . .
where the discrete Newton increments δ(·) j (the iteration superscript is discarded) solve the non-symmetric Jacobian linear problem We have used the bilinear forms a 1 : where the subscripts explicitly indicate fixed quantities; and linear functionals F k : H 1 (Ω) → R, G k , H k , I k : H 1 (Ω) → R where the subscript k denotes that they depend on quantities associated with the state around which one performs linearisation. The forms satisfy the following specifications If one uses (2.6a)-(2.6b) then the third column and the third row (3.9c) in the tangent matrix are not needed.

Preliminaries and convergence verification
All tests in this section have been implemented with the open-source finite element library FEniCS [3]. After matrix assembly, the resulting linear systems at each linearisation step are solved with the direct solver MUMPS. The stopping criterion on the nonlinear iterations of the Newton-Raphson algorithm is based on a weighted residual norm dropping below the fixed tolerance of 1 · 10 −7 . Before computing numerical solutions pertaining to the application of calcium wave trains in 2D and 3D, we briefly address the verification of convergence for the space-time discretisation. For this we simply consider a unit square domain, partitioned uniformly into triangles of successive refinement. On each resolution level we compute approximate solutions and compare (using errors for displacement, pressure, calcium concentration and non-activated IPR in their natural norms), against the following manufactured exact solutions For the convergence tests we use the Taylor-Hood element for displacement and pressure, and continuous and piecewise quadratic elements for h and c. The parameters are ν = 0.49, α 1 = α 2 = 0.001, β 1 = β 2 = µ = λ = b = G = K 1 = 1, K = 2, D = 0.1. The exact displacement in the viscoelastic case does not satisfy a homogeneous traction condition, so a synthetic traction is imposed computed from the exact stress. In addition, the exact calcium concentration is imposed as a Dirichlet datum everywhere on the boundary. A time step ∆t = 0.01 is chosen and we simulate a short time horizon t final = 3∆t. Errors e s between the approximate and exact solutions are tabulated against the number of degrees of freedom in Table 4.1. For the spatial convergence verification we have used f (t) = t. The expected convergence behaviour (quadratic rates for all fields) is observed. The Newton-Raphson algorithm takes, in average, four iterations to reach the prescribed residual tolerance. The convergence in time is verified with f (t) = sin(t) and by partitioning the time interval into successively refined uniform steps and computing accumulated errorsê s . Here   where · denotes the appropriate space norm for the generic vector or scalar field s (that is, the L 2 −norm for pressure and the H 1 −norm for the remaining variables). The results are shown in Table 4.2, confirming the expected firstorder convergence. Samples of approximate solutions are depicted in Figure 4.1.

Test 1. Calcium waves on a fixed domain
The initial conditions are as in (2.5), in particular, for Tests 1-3 the initial condition for h is the steady state and for calcium consists of a Gaussian on the domain centre and the steady state elsewhere, as follows , whereas for the 3D cases we will use either (2.6a)-(2.6b) or (2.7). In all cases we use fixed timesteps dictated by the reaction kinetics. Let us consider as domain the disk centred at (0, 0) with radius 2.5. The system evolves in time until t final = 10. For this case the coupling constant λ is taken simply as zero. We observe that a single calcium pulse propagates from the centre of the domain towards the boundary. The amount of IPR in the system also has an influence on the spatio-temporal patterns of calcium. Examples of concentrations at three different times are depicted in Figure 4.2,

Test 2. Effects of mechanochemical coupling
Setting now a different value of the coupling constant λ = 0.35 results in the behaviour illustrated in Figure 4.3. No other parameters have been changed and we can see significant effects from the mechanochemical coupling. In particular, we observe that additional secondary waves are initiated from the boundary (which is the part of the domain where most pronounced deformation and dilation occurs), and these are sustained over a longer time, suggesting a long-range contraction response. We also select some other parameter regimes that, according to the 1D theory, would produce qualitatively interesting behaviour such as oscillatory calcium transients. In connection with Test 1 above, note that in the uncoupled scenario of λ = 0, the cases µ = 0.288468 and µ = 0.3 generate a solitary wave and a periodic wave train, respectively. On the other hand, if we fix µ = 0.288468 and vary λ we distinguish between the following cases: Case I: λ = 0.1 and an equilibrium calcium state predicted by the 1D theory of c s = 0.4850 produces a periodic wavetrain. Case II: Setting λ = 0.5 and c s = 19817 we see that the periodic wavetrain has a similar period as before but it decays much faster. For these rounds it suffices to take a timestep of ∆t = 0.2. If we set now a larger µ = 0.3 we see that a smaller ∆t = 0.1 is needed to capture the faster domain velocity. The values λ = 0.1, c s = 0.6034 (Case V) also produce a periodic wavetrain, whereas increasing λ further will still show wavetrains but decaying faster and faster (Cases VI, VII, VIII). For λ = 2, c s = 1.2506 (case VIII), we can also use a larger timestep ∆t = 0.2. A few samples of these findings are shown in Figure 4.4.

Test 3. Including calcium sparks
Assuming that calcium sparks appear on the continuum randomly, as experiments show, we can study their influence on the calcium propagation through the modification of (2.4c) as follows where F (·, ·, ·) is a generic reaction term specified by the sought model (e.g. [5]), and I f encodes a collection of flashes of different amplitudes, having a uniform random distribution. From a phenomenological perspective, if the impulses are applied sufficiently close to each other, a simple superposition phenomenon explains why the intensity of individual local pulses contributes to form a synchronous wave propagating front. These impulses are implemented as point sources (Dirac delta functions on a linear variational form) localised at randomly distributed points in Ω, with amplitude equal to the constant calcium equilibrium found for each parameter configuration, and switched on at a given frequency. They are projected onto the appropriate finite element space (the one used for calcium approximation) and applied after assembling of the residual that constitutes the right-hand side of the tangent system (3.9), at each Newton iteration.
In the case of λ = 0, although the advection effect is still turned on, the calcium transients do not exhibit substantial differences with respect to their counterparts on a fixed domain (that is, when also the calcium-dependent contraction is turned off, with β 2 = 0). Also, the structure of the waves (when regarded locally) does not change with respect to the observed behaviour in the case of a single calcium stimulus (meaning that the nature of solitary waves and periodic wavetrains also occurs in the case of spark phenomena). However, the nucleation does plays a role, as anticipated above. We show in Figure 4.5 snapshots of the calcium concentration at four time instants and compare the effect of increasing the frequency of the applied stimuli, and see that isolated point sources initiate a propagating wave but eventually decay when the frequency is low enough (this is the expected trend observed also in the one-dimensional model from [28]), while as soon as the frequency is increased from one each 20 time steps to one each 5 time steps, the solitary waves form synchronous wavefronts that are sustained in time.  Then we simulate for larger values of λ and examine the differences in body deformation and calcium distribution. Samples of these computations are presented in Figure 4.6. We see that even with lower frequencies of calcium sparks (here taking one each 10 time steps) one is able to generate propagating fronts that in turn induce a periodic (but not uniform) dilation of the cell. In addition, if we regard the domain as a multi-cell aggregate instead as of a single cell, the short-time contractions and sparks could represent single-cell contractions whereas longer-term contractions and propagating fronts would account for synchronous cell movement and rearrangement and collective, larger intercellular calcium wavetrains. It is also clear that the deformations are no longer radial, which shows the benefit of incorporating multidimensional models.
We finish by replicating the behaviour of Test 3B in different geometries. A microscope image of a Xenopus embryonic tissue from [29] is segmented and an unstructured triangular mesh is generated. The sparks are applied now with a lower intensity (75% of c s ). In Figure 4.7 we show a sample of coarse mesh and a few snapshots of the calcium distribution over the deformed domain. The tissue expands and contracts with an initial period of approximately 5 nondimensional units, and then the contraction increases in frequency. We also plot the displacement and pressure distribution, as well as the magnitude of the strain tensor, at t = 10. We also consider a thin cylinder of radius 2.5 and height 0.4. The boundary is split into the bottom disk Γ where we impose the slip condition (2.6a) and on the remainder of the boundary we prescribe zero traction (2.6b). Orthogonality with respect to the space (2.8) is therefore not required in this case. These mechanical boundary conditions represent a tissue that slips on a substrate, and in order to prevent it from slipping away we also prescribe zero displacement on the origin. Apart from a larger diffusion D = 0.02 and a higher Poisson ratio ν = 0.45, all other mechanochemical parameters are maintained as before. We compare qualitatively the deformation patterns with respect to imposing pure traction boundary conditions as before and show the outcome in Figure 4.8. The former boundary treatment leads to mechanochemical patterns similar to those encountered in some of the 2D cases reported in Test 3B, whereas for the latter boundary conditions the deformations show a slightly more pronounced displacement in the z-direction.

Concluding remarks
In this work we have extended to higher spatial dimensions the mechanochemical model in [28] representing a two-way coupling between calcium signalling and mechanics in one dimension. We have taken special care of traction boundary conditions for the motion of the cell. The governing equations couple an advection-reaction-diffusion system for calcium signalling with a linear viscoelastic material. We include the two-way coupling mechanism through volume-dependent terms and an active contraction stress which is dependent on calcium dynamics as in [28,41]. We have also established the existence of weak solutions, and we have proposed a finite element discretisation.
Our model enables us to gather detailed insights in the mechanochemical processes during embryogenesis. More precisely, a quantitative evaluation of the effects of viscoelasticity and local dilation as mechanochemical processes underlying cytosolic calcium waves was undertaken. Regarding future directions, we note that in [29] the authors suggest that cells detect stress rather than stretch or strain because one does not observe transient stretch as contractions spread from cells exposed to ATP. This indicates that probably an appropriate formalism for further elucidating calcium signalling coupled to mechanics of cells and tissues is the framework of stressassisted diffusion models, recently proposed in [16,34,48]. This relates also to the observation that in other types of cells that proliferate at a rate dependent on the substrate stiffness [20]. Furthermore, at certain stages of morphogenesis, the assumption of small strains is no longer adequate and one needs to describe the motion through nonlinear elasticity. One could explore growth models using multiplicative decompositions of deformation gradients [27,46], which will strongly depend on the type of phenomenon under consideration. Finally, it would be of interest to study the contraction and elongation of cellcell gap junctions on the apical end of the tissue, since these regions contain a clustering of epithelia, and an important part of the overall morphogenesis occurs therein.
The two-way coupling between mechanical forces and calcium signalling at both individual and collective cell levels is of course not unique to embryogenesis-it is a phenomenon shared by many other biological systems such as cancer cells, wound healing, keratinocytes, or white blood cells [19,30,58]. Modifications to the theoretical and computational tools presented here could be used to study these other processes too.
Work is underway to extend the present model and method formulation to simulate the contact of the cell with a surrounding fluid (using an immersed boundary finite element method with distributed Lagrange multipliers) as well as the cell-to-cell interactions using a virtual element discretisation capturing more effectively the single cell geometries and the boundary contraction and elongation that together with junctional tension comprise the tissue-level deformation [25]. Moreover, since calcium signalling evolves rapidly it would be natural to explore time-adaptive schemes based on local error estimators such as those advanced in [15] for a similar application.