Abstract
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.
Similar content being viewed by others
Avoid common mistakes on your manuscript.
1 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, leading to complex spatio-temporal patterns. The coupling between rearrangement of tissue, cell migration and active cell contraction to diffusion of morphogens or signalling molecules has been proposed and studied in earlier works [46, 47] but many challenges remain open. 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 [60]. Moreover, localisation of stresses or strains within the cells can generate alteration in patterns of calcium distribution in tissue by changing cell displacement magnitude, direction and velocity [27, 39]. On the other hand, recent experimental evidence [16, 48, 61] shows that tissue mechanics are actively affected by calcium signalling during development.
Hence, more modelling and analysis has recently appeared in order to elucidate the many open questions that surround the healthy evolution of an embryo [14, 33, 48, 61]. A variety of models have described general interactions between chemical species and mechanics; these include descriptions where stress is triggered by chemical signalling [47], or mainly by migration [44]. Here we are interested in the mechanochemical feedback due to the pressure and velocity of the underlying embryonic tissue. Including the coupling with continuum mechanics significantly affects the propagation of the chemical signals. Disrupting calcium signals in apical constriction leads to embryo malformations such as Spina Bifida, which affects millions of babies worldwide [16, 61].
We develop a multi-dimensional extension of the recent mechanochemical model proposed in [33] which describes the interplay of calcium signalling with the mechanics of embryonic epithelial tissue during apical constriction, an active deformation process. We are inspired by the recent, interesting experimental observations in [16, 61] where increasing tension in the contracting cells yields further calcium release and this, in turn, elicits contractions in the cell 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 [20, 21, 51], poroelastic [19, 43, 56, 57], or nonlinear and nonlocal characteristics [42, 65] but more predominantly, viscoelastic effects [2, 8, 10, 24, 34, 53, 66]. 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 [50]. In the present mechanochemical model, we assume that the viscoelastic stress includes an active tension component which depends on the concentration of calcium, following the formulation in [6, 33, 47].
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 or solid pressure (the so-called Herrmann formulation [30] where the introduction of solid pressure contributes to achieve robustness in the nearly incompressible regime. This occurs when the first Lamé parameter defining the elastic properties of the material is very large); b) the spatio-temporal dynamics of calcium concentration and the percentage of IP3 (inositol triphosphate) 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 [33, 45, 46, 49, 59].
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 rescaled hydrostatic 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 [37] (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 [33] 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 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.
2 Model Description and Weak Formulations
2.1 Coupling Calcium Signalling with Solid Mechanics
We assume that the cell/tissue can be macroscopically regarded as a linear, viscoelastic material of Kelvin–Voigt type, occupying the spatial domain \({\Omega }\subset \mathbb {R}^{d}\) with d = 2 or d = 3, and having a smooth boundary ∂Ω with outward-pointing unit normal n. Boldface letters will indicate either vector or tensor-valued fields and differential operators. For instance, ∇ denotes the gradient operator applied to scalar fields and ∇ the gradient operator applied to vectors. Similarly, by div and div we will denote the divergence operator acting on vector and tensor-valued fields, respectively. The symbol I will be used for the second-order identity tensor.
Following, e.g., [33, 44], and assuming that gravitational forces and inertial effects are negligible, one seeks for each time t ∈ (0,tfinal], the displacements of the tissue, \(\boldsymbol {u}(t):{\Omega }\to \mathbb {R}^{d}\), and the dilation \(\theta (t):{\Omega }\to \mathbb {R}^{d}\), such that
where σ is the Cauchy stress tensor, \(\boldsymbol {\varepsilon }(\boldsymbol {u})={\frac {1}{2}(\boldsymbol {\nabla } \boldsymbol {u}+\boldsymbol {\nabla } \boldsymbol {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. c is the concentration of calcium and here, 𝜃 = divu represents the dilation/compression of the apical surface area of the cell. 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 \(\tilde \alpha _{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 an active stress which is dependent on c. As in [33] and [47] we assume a Hill function which saturates to a constant value T0 for high values of c, in line with experimental observations [16]. Therefore,
where n is a positive integer and κ > 0. We, thus, assume that the calcium concentration affects the material’s motion by regulating the active tension and therefore the generation of stress. The effect is an activation and not an inhibition of active tension and so T0 is assumed positive.
As in [44], 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, the other variable of the system is h, which represents the percentage of non-inactivated IPR on the ER. The model is a generalisation of the recent model in [33], which employs the calcium dynamics of the well-known model in [5]. In dimensional form, the model is written as follows
where
Here D is the calcium diffusion coefficient inherent to the tissue, which is positive and assumed constant. The fluxes in (2.3a) are as follows: the term JER models the flux of calcium from the ER into the cytosol through the IPR, μ(p3) = p3/(kμ + p3) is the fraction of IPR that have bound IP3 and is an increasing function of p3, the IP3 concentration. The constant kf denotes the calcium flux when all IP3 receptors are open and activated, and b represents a basal current through the IPR; Jpump 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. We will explore the bifurcations of the system as the IP3 concentration, and correspondingly, μ increases. Note that one could also include calcium fluxes leaking into the cytosol from outside the cell, but we leave those terms out, as they can be assumed small. The term JSSCC encodes the calcium flux due to the activated stretch-sensitive calcium channels (SSCC). These channels have been identified experimentally in recent years [28]—they are on the cell membrane. SSSC are activated when exposed to mechanical stimulation and they close either by relaxation of the mechanical force or by adaptation to the mechanical force [28]. The constant S represents the strength of stretch activation, and it encodes mechanochemical coupling. Since S multiplies 𝜃 it is a measure of the strength of the mechanochemical coupling we are studying here. We will be later treating S as our second bifurcation parameter.Footnote 1
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, is taken as τj > 1 in the ODE (2.3b) governing the dynamics of h. As in [33] and [5] we adopt the value τj = 2s. 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 [41]). In addition, instead of the dilation 𝜃, we use the rescaled Herrmann pressure p. (This does not coincide with the physical pressure, but will be referred to as rescaled hydrostatic pressure from now on.) Doing this leads to the following system, where we have dropped the bars in the dimensionless variables for notational convenience,
where D⋆ = Dτj/L2, K1 = kfτj/k1, G = γτj/k1, K = kγ/k1, and λ = τjS/k1. For the parameter values we have chosen from [1, 5] and also taking D = 20μm2/s and L = 100μ m we obtain the following values for the nondimensional parameters D⋆ = 4 × 10− 3, K1 = 46.28571428, G = 40/7, K = 1/7. Note that the relatively large value of K1 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/Herrmann pressure. 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. [49]). This parameter will be treated as second bifurcation parameter, as in [33]. In [33] it has been shown that as λ increases the system undergoes bifurcations and calcium signals qualitatively change, e.g. going from an action potential to a limit cycle (calcium oscillation). Moreover at a critical value of λ, oscillations were shown to be suppressed in [33]—we expect a similar behaviour in the higher dimensional case here. It is not clear what the maximum λ value should be in higher dimensions. Inspired by the results in [33] we expect the oscillations/waves should still vanish for large enough λ and so in the parametric study for the 2D numerical tests we vary λ from 0 to λ ≈ 2. In Fig. 5 we indeed show that calcium oscillations have been suppressed for λ = 2.
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 [33] 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 (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 [13, Eq. (11.1.7)])
and unique solvability (for a given calcium concentration c) will follow since \(\mathbb {R}\mathbb {M}({\Omega })\) is a null space of the relevant functional space.
2.2 Mixed-primal Weak Formulation
Multiplying the nondimensional governing equations (2.1) 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 u(t) ∈V, p(t) ∈ L2(Ω), c(t),n(t) ∈ H1(Ω), such that
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
These nonlinear terms are the non-dimensional counterparts of the Hill function (2.2), the sum of the reactions JER, Jpump, JSSCC in (2.3a) and the reaction in (2.3b), respectively.
2.3 Fully Discrete Form
Let us consider a shape-regular partition \(\mathcal {T}_{j}\) of \(\bar {\Omega }\) into affine elements (triangles in 2D or tetrahedra in 3D) E of diameter jE, where \(j = \max \limits \{ j_{E}: E\in \mathcal {T}_{j}\}\) denotes the meshsize. Finite-dimensional subspaces for the approximation of displacement, rescaled hydrostatic pressure, and calcium are specified as
where \(\mathbb {P}_{k}(E) \) denotes the space of polynomials of degree less than or equal than k defined locally over the element \(E \in \mathcal {T}_{j}\), and bE := φ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 Vj or \(\widetilde {\mathbf {V}}_{j}\), the pairs Vj × Qj are known as the MINI element and \(\widetilde {\mathbf {V}}_{j}\times Q_{j}\) are known as the Taylor–Hood elements. They are both Stokes inf-sup stable (see, e.g., [11, 55]). On the other hand, an L2-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 \(\mathbb {I}\) encoding rotational kinetic energy [37]. Denoting \(\boldsymbol {x}_{0} = |{\Omega }|^{-1} {\int \limits }_{\Omega } \boldsymbol {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 \(\mathbb {I}\) (see also [36]). As remarked in [37], 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 rescaled hydrostatic pressure the pair \(\mathbb {P}_{2}-\mathbb {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,tfinal] into equi-spaced points tk = kΔt. Applying a backward Euler method for the first-order time derivatives, we can write a semidiscrete form of (2.2) (but now incorporating the orthogonality with respect to rigid motions with a Lagrange multiplier). Now the fully discrete formulation is given by
3 Well-posedness Analysis
We will consider that the initial data are nonnegative and regular enough, hence the numerical scheme starts from discrete initial data \(\boldsymbol {u}_{j}^{0}\), \({p_{j}^{0}}\), \({c_{j}^{0}}\), \({h_{j}^{0}}\), which are the projections of the exact initial conditions onto the finite element spaces.
Let us introduce the trilinear form \(d: \boldsymbol {H}^{1}({\Omega })\times H^{1}({\Omega }) \times H^{1}({\Omega }) \to \mathbb {R}\) and bilinear form \(e: H^{1}({\Omega })\times H^{1}({\Omega }) \to \mathbb {R}\) defined as:
Here we rewrite the nonlinear term u ⋅∇ϕ in its equivalent skew-symmetric form \(\boldsymbol {u}\cdot \nabla \phi +\frac {\phi }{2}\text {div} \boldsymbol {u}\) [25, Chapter IV, Lemma 2.2]. Then, due to the skew-symmetric form of the operator d(⋅;⋅,⋅), we can write
And recalling the Poincaré–Friedrichs inequality [26, Chapter I, Lemma 3.1], we have the coercivity for e(⋅,⋅):
We also assume the nonlinearities in (2.10) satisfy the following conditions:
-
(A1)
J(c) is uniformly bounded, i.e. |J(c)|≤ CJ.
-
(A2)
K(c,h) is uniformly bounded with respect to c and Lipschitz with respect to h, in particular we have: |K(h,c)|≤ Ck|h|.
-
(A3)
β(c) is uniformly bounded, i.e. |β|(c)| < Cβ.
Moreover, we will use the following algebraic relation: Let a and b be two real numbers, then we have
3.1 Continuous Dependence on Data
Theorem 3.1
Let \((\boldsymbol {u}^{k+1}_{j},p^{k+1}_{j},c^{k+1}_{j},h^{k+1}_{j}) \in (\boldsymbol {V}_{j},Q_{j},{\Phi }_{j},{\Phi }_{j})\) be a solution of problem (2.3), with initial data \((\boldsymbol {u}^{0}_{j},{c^{0}_{j}},{h^{0}_{j}})\). Then under the assumption that \(\left (1-\frac {\sqrt {d}}{2}\frac {1-4\nu }{4\nu -3}\right )>0\), the following bounds are satisfied:
where C1, C2 and C3 are positive constants independent of j and Δt.
Proof
On the second equation of problem (2.3) we choose \(q=p^{k+1}_{j}\), to get
which after applying Young’s inequality, implies
hence,
Now, on the first equation in problem (2.3), we choose the test function v = uk+ 1, then multiply by 2Δt, and use the 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 \(\left (1-\frac {\sqrt {d}}{2}\frac {1-4\nu }{4\nu -3}\right )>0\) and (A3), we get
Now we take \(\phi = c^{k+1}_{j}\), use properties (3.1), (3.2) and multiply by 2Δt the third equation in problem (2.3), to deduce:
As before, summing over k from 0 to n ≤ N − 1, and applying Young’s inequality we deduce
By assumption (A2) we have,
Similarly in the fourth equation of (2.3), we take \(\psi = h^{k+1}_{j}\), multiply by 2Δt and use property (3.1) to get
Hence, applying Young’s inequality, assumption (A1) and summing over k, we get
While the first and third results follows directly from (3.6) and (3.8) respectively, we get the second result by substituting the bounds for \(h_{j}^{k+1}\) and \(\boldsymbol {u}_{j}^{k+1}\) into (3.7). □
3.2 Existence Result: Fixed-point Approach
Next, we address the unique solvability of problem (2.3). To that end, we will make use of Brouwer’s fixed-point theorem in the following form (given by [25, Corollary 1.1, Chapter IV]):
Theorem 3.2
(Brouwer’s fixed-point theorem) Let H be a finite-dimensional Hilbert space with scalar product (⋅,⋅)H and corresponding norm ∥⋅∥H. Let Φ: H → H be a continuous mapping for which there exists μ > 0 such that (Φ(u),u)H ≥ 0 for all u ∈ H with ∥u∥H = μ. Then there exists u ∈ H such that Φ(u) = 0 and ∥u∥H ≤ μ.
Theorem 3.3
(Existence of discrete solutions) Under the same assumptions as in Theorem 3.1, problem (2.3) with initial data \((\boldsymbol {u}_{j}^{0}, {p_{j}^{0}}, {c_{j}^{0}},{h_{j}^{0}})\) admits at least one solution.
Proof
To simplify the notation in the proof we introduce the constants:
We proceed by induction on k ≥ 1. We define the mapping
using the relation
Note that this map is well-defined and continuous on Vj × Qj ×Φj ×Φj. On the other hand, if we take
and employ (3.1), (3.2), (3.4), together with assumptions (A1)–(A3), we obtain
Next, using Korn’s inequality (with constant \(\hat {C}_{k}\)) and (3.5), we deduce that
Then, setting
we may apply the inequality \(a+b\leq \sqrt {2}(a^{2}+b^{2})^{1/2}\), valid for all \(a,b\in \mathbb {R}\), to obtain
Hence, the right-hand side is nonnegative on a sphere of radius r := Cr/CR. Consequently, by Theorem 3.2, there exists a solution to the fixed-point problem \({\Psi }(\boldsymbol {u}^{k+1}_{j},p^{k+1}_{j},c^{k+1}_{j},h^{k+1}_{j}) = 0\). □
3.3 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 \(\boldsymbol {u}_{j}^{0}\), \({p_{j}^{0}}\), \({c_{j}^{0}}\), \({h_{j}^{0}}\) (projections of the exact initial conditions onto the finite element spaces), and for ℓ = 1,…, we find \(\boldsymbol {u}_{j}^{\ell +1}\in \mathbf {V}_{j}\), \(p_{j}^{\ell +1}\in Q_{j}\) and \(c_{j}^{\ell +1},h_{j}^{\ell +1}\in {\Phi }_{j}\), 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}:\mathbf {H}^{1}({\Omega })\times \mathbf {H}^{1}({\Omega })\to \mathbb {R}\), \(a^{2}: L^{2}({\Omega })\times L^{2}({\Omega })\to \mathbb {R}\), \(a^{3},a^{4},e^{3},e^{4}: H^{1}({\Omega })\times H^{1}({\Omega })\to \mathbb {R}\), \(b:\mathbf {H}^{1}({\Omega })\times L^{2}({\Omega })\to \mathbb {R}\), \(g: \mathbf {H}^{1}({\Omega })\times \mathbb {R}\mathbb {M} \to \mathbb {R}\), \(d^{1}:H^{1}({\Omega }) \times \mathbf {V}\to \mathbb {R}\), \(d^{3}:\mathbf {H}^{1}({\Omega })\times H^{1}({\Omega }) \to \mathbb {R}\), \(d^{4}:\mathbf {H}^{1}({\Omega })\times H^{1}({\Omega }) \to \mathbb {R}\), where the subscripts explicitly indicate fixed quantities; and linear functionals \(F_{k}:\mathbf {H}^{1}({\Omega })\to \mathbb {R}\), \(G_{k},H_{k},I_{k}: H^{1}({\Omega })\to \mathbb {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.
4 Numerical Results
In this section we solve the system (2.1) and study its behaviour as the two bifurcation parameters, λ and μ vary. All tests here 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.
From now on, and unless otherwise specified, we take values for all model constants as follows: D⋆ = 0.004, ν = 0.4, α1 = 1, α2 = 0.5, β1 = 1.5, β2 = 0.1, K1 = 46.29, K = 0.1429, G = 5.7143, b = 0.111. These parameter values are appropriate for the 2D and 3D models extending the 1D model in [5, 33]. The approximation of displacement – rescaled hydrostatic pressure will be restricted only to the MINI element and to continuous and piecewise linear elements for the calcium concentration and IPR (that is, Vj × Qj ×Φj ×Ψj in (2.11)).
We will first conduct three tests, Test 1, 2, 3. Tests 1A, 1B, 2A, 2B, 3A, 3B, 3C are over a 2D domain and Test 3D is over a cylindrical domain. The boundary conditions for the 2D domains are taken as (2.7), whereas for the 3D domain we will use either (2.6a)–(2.6b) or (2.7). In all cases we use fixed timesteps dictated by the reaction kinetics.
For the 2D cases, we consider as domain a disk centred at (0,0) with radius 2.5. This domain approximates the apical surface area of an epithelial cell, which typically has a diameter of 50μ m. In the initial conditions (2.5), we take h to be at its steady state value and for the calcium concentration we take a thin Gaussian pulse at the domain centre and the steady state value elsewhere, as follows:
Note that cs changes with μ: cs = 0.18572 for μ = 0.288, cs = 0.55633 for μ = 0.3, and cs = 0.9658 for μ = 0.75. The initial disturbance is sufficient to induce a calcium wave, in an appropriate parameter regime. We let the system evolve in time until tfinal = 10.
4.1 Test 1. Calcium Waves on a Fixed Domain
In Test 1, λ is taken to be zero. Then concentration of IP3 in the system (parameter μ) solely dictates the bifurcations in the spatio-temporal patterns of calcium.
The calcium concentration and h are plotted at three different times in Fig. 1, using μ = 0.288,0.3 and 0.75. For μ = 0.288 (panels (a–d)) the Gaussian pulse evolves into a solitary wave (which eventually reaches the boundary of the domain). Panels (e–h) indicates that a wavetrain is emerging. (Compare the effect of a larger domain, in the I, where Fig. 11, shows a fully-fledged periodic wavetrain). For μ = 0.75, panels (i–l) the initial disturbance decays when moving away from the disc centre. While the plots in Fig. 1 highlight the spatial distribution of the calcium concentration, its dynamic behaviour at a fixed point is more clearly exhibited Fig. 2, where we plot the transients of both c and h at the point P : (0,1), for the three values of μ considered above. We see, as expected from [5] that for μ = 0.288 there is an action potential, for μ = 0.3 there are limit cycles and for μ = 075 the transients decay to the steady state in an oscillatory manner.
4.2 Test 2. Effects of the Mechanochemical Coupling
- Test 2A::
-
Increasing λ to 0.35, and leaving all other parameters the same, causes pronounced changes to the calcium pattern, as illustrated in Fig. 3. Going from (a) to (b) we see that the solitary wave becomes a periodic wavetrain. In panel (c), waves also seem to be initiated from the boundary (which is the part of the domain where most pronounced deformation and dilation occurs). These waves are sustained over a longer time, suggesting a long-range stretch response. In panels (d), (e), (f) we plot the rescaled hydrostatic pressure to study the cell’s deformation associated with the observed calcium patterns. We see that there is a pronounced positive dilation of the cell as we go from (d) to (f). Hence, we see that it is important to take the mechanochemical coupling into account. Note that waves from the boundary arise in the application of macroscopic wound healing modelling. Experimental observations and models from, e.g., [31, 38, 63, 64] indicate that active stress in visco-elastic layers onsets more pronounced deformations near the regions adjacent to the wound, that is, near the boundary of the epithelial layer.
- Test 2B::
-
We also select some other parameter regimes that, according to the 1D model [33], we expect that they will produce qualitatively interesting behaviour. In Test 1 (λ = 0), we choose μ = 0.288 and μ = 0.3 which leads to a solitary wave and a periodic wavetrain, respectively. Now, fixing μ = 0.288 we increase λ and we distinguish between the following cases: Case I: λ = 0.1 and cs = 0.4850 leads to a periodic wavetrain. Case II: Setting λ = 0.5 and cs = 0.6824 gives also a periodic wavetrain with a higher frequency. Case III: With λ = 1.3 and cs = 0.9577 we also find a periodic wavetrain that gradually decays, but where also the resting value elsewhere on the domain increases with time. Case IV: With λ = 2 and cs = 1.19817 we see that the periodic wavetrain has a similar period as before but it decays much faster.
For these simulations it suffices to take a timestep of Δt = 0.2. If we set now a larger μ = 0.3, the values λ = 0.1(cs = 0.6034 (Case V) also lead to a periodic wavetrain. Increasing λ further still leads to wavetrains but these decay faster and faster (Cases VI, VII, VIII). For λ = 2,cs = 1.2506 (case VIII), we can also use a larger timestep Δt = 0.2. We select some of the results to show in Fig. 4, where we can see the viscoelastic response of the domain interacting with the calcium dynamics.
Similarly to Fig. 2, in Fig. 5 we also show transients of all field variables, at the fixed point P:(0,1) for some of the parameter values mentioned above. Comparing Fig. 2a with Fig. 5a we see that for μ = 0.288 as λ increases the action potential of calcium is transformed to limit cycles. Then comparing Fig. 5a with Fig. 5b we see that as λ increases further the amplitude and frequency of the calcium limit cycles both increase. When λ = 2 Fig. 5c shows that the limit cycles have been suppressed and oscillations decay. A very similar transition in the qualitative behaviour is observed in Fig. 5b, d, f for μ = 0.3 as λ increases. Of course when μ = 0.3 limit cycles already exist when λ = 0 and their amplitude and frequency increases as we can see by comparing Fig. 2b with Fig. 5b. We see that the dilation undergoes the same transitions as λ increases from λ = 0.1 to λ = 2, for μ = 0.288 and for μ = 0.3. One important feature to note here is that the cytosolic stretches occur a little after the calcium spike - this has been observed experimentally in [16] and [61].
4.3 Test 3. Including Calcium Sparks
Assuming that calcium sparks appear on the cell in a random manner, as experiments show [16, 61], we can study the influence of these sparks 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 If 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.3), 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 deformation 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 Fig. 6 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 [33]), 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 Fig. 7. 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 deformations and sparks could represent single-cell stretch whereas longer-term stretch and propagating fronts would account for synchronous cell movement and rearrangement and collective, larger inter-cellular calcium wavetrains. It is also clear that the deformations are no longer radial, which shows the benefit of incorporating multidimensional models.
Finally, we replicate the behaviour of Test 3B in different geometries. A microscope image of a Xenopus embryonic tissue from [34] is segmented and an unstructured triangular mesh is generated. The sparks are applied now with a lower intensity (75% of cs). In Fig. 8 we show a sample of coarse mesh and a few snapshots of the calcium distribution over the deformed domain. The tissue stretches and contracts back with an initial period of approximately 5 nondimensional units, and then the contraction increases in frequency. We also plot the displacement and rescaled hydrostatic 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 Fig. 9. 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.
4.4 Accuracy Verification
To close this section, 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, rescaled hydrostatic pressure, calcium concentration and non-activated IPR in their natural norms), against the following manufactured exact solutions
with \(p = -\frac {\nu }{1-2\nu }\text {div} \boldsymbol {u}\). For the convergence tests we use the Taylor–Hood element for displacement and rescaled hydrostatic pressure, and continuous and piecewise quadratic elements for h and c. The parameters are ν = 0.49, α1 = α2 = 0.001, β1 = β2 = μ = λ = b = G = K1 = 1, K = 2, D⋆ = 0.1. The exact displacement in the viscoelastic case does not satisfy a homogeneous traction condition, but still a pure-traction boundary condition is imposed by prescribing a synthetic traction computed from the exact stress. In addition, the exact calcium concentration is imposed as a Dirichlet datum everywhere on the boundary. Other combinations of boundary conditions (such as mixed loading) were also tested, and they do not alter the convergence properties. A time step Δt = 0.01 is chosen and we simulate a short time horizon tfinal = 3Δt. Errors es between the approximate and exact solutions are tabulated against the number of degrees of freedom in Table 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 \limits (t)\) and by partitioning the time interval into successively refined uniform steps and computing accumulated errors \(\hat {\texttt {e}}_{s}\). Here
where ∥⋅∥⋆ denotes the appropriate space norm for the generic vector or scalar field s (that is, the L2-norm for rescaled hydrostatic pressure and the H1-norm for the remaining variables). The results are shown in Table 2, confirming the expected first-order convergence. Samples of approximate solutions are depicted in Fig. 10.
5 Concluding Remarks
In this work we have extended to higher spatial dimensions the 1D mechanochemical model in [33] which captures a two-way coupling between calcium signalling and mechanics. We have taken special care of traction boundary conditions for the motion of the cell. The governing equations consist of an advection-reaction-diffusion system for calcium signalling coupled to the mechanical equation for a linear viscoelastic material. We implement the two-way coupling through volume-dependent terms and an active contraction stress which is dependent on calcium dynamics, as in [33, 47]. We have also established the existence of weak solutions, and we have developed and implemented a multidimensional mixed-primal finite element discretisation for coupled chemo-viscoelasticity with pure traction boundary conditions, which is robust with respect to the Lamé parameters, and that captures the salient features of the model solutions in an accurate and efficient manner.
Our model generates detailed insights into the mechanochemical processes, which are very important during embryogenesis. Approximating a single epithelial cell with a circular disc, we have explored the parameter space by varying the level of IP3 in the system (μ parameter), the usual bifurcation parameter in calcium signalling models, and the strength of the mechanochemical coupling (captured here through the parameter λ.) We have shown that the system exhibits a rich and interesting behaviour, which includes solitary waves and periodic wavetrains. We observe that coupling calcium signalling to mechanics generates significant differences in the calcium signals; hence a key takeaway of our work is that including the mechanochemical coupling when studying calcium signals is important. Moreover, a quantitative evaluation of the effects of viscoelasticity and local dilation involved in these mechanochemical processes underlying cytosolic calcium waves has been undertaken.
Additionally, we have explored the presence of randomised calcium sparks in the domain and have shown that these can organise into calcium waves under appropriate conditions. This is a question of interest since different calcium signals have been observed in experiments, ranging from sparks to waves - see for example, [16]. Taking a step closer to experiments we have run simulations over a sample of a deforming Xenopus tissue in 2D and then also considered a cylindrical surface in 3D.
A direct extension of our work here would be to employ the model in order to study calcium signalling and mechanics over a a multicellular domain, comprising hundreds of cells. In this connection, we are in contact with experimentalists working on apical constriction [16]. Elucidating the precise interplay between calcium signalling and mechanics during apical constriction would shed light into embryo malformations, such as Spina Bifida.
As it stands, the present model has a number of limitations. For example, we note that in [34] the authors suggest that cells detect stress rather than stretch or strain because one does not observe transient stretch spreading 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 stress-assisted diffusion models, recently proposed in [15, 40, 54]. This relates also to the observation that in other types of cells that proliferate at a rate dependent on the substrate stiffness [24]. 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 [32, 52, 64], which will strongly depend on the type of phenomenon under consideration. For instance, the adaptation of growth and remodelling theory to the present context could follow the formalism of constrained mixtures developed in, e.g., [12, 17, 18]. The theoretical analysis of stability extended to that case presents a major challenge.
Finally, it would be of interest to study the contraction and elongation of cell-cell 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. Modifications to the theoretical and computational tools presented here could be used to study these and other processes as well.
Evidently, 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, keratinocytes, or white blood cells [23, 35, 67].
Further work is underway to extend the present model and method formulation in order 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 which would capture more effectively the single cell geometries and the boundary contraction and elongation that together with junctional tension comprise the tissue-level deformation [29]. 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 [20] for an application on morphoelasticity.
Change history
07 December 2022
The original online version of this article was revised: Missing Open Access funding information has been added in the Funding Note.
References
Allbritton, N.L., Meyer, T., Stryer, L.: Range of messenger action of calcium ion and inositol 1, 4, 5-trisphosphate. Science 258, 1812–1815 (1992)
Allena, R., Muñoz, J. J., Aubry, D.: Diffusion-reaction model for Drosophila embryo development. Comput. Methods Biomech. Biomed. Eng. 16, 235–248 (2013)
Alnæs, M.S, Blechta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A., Richardson, C., Ring, J., Rognes, M.E., Wells, G.N.: The FEniCS project version 1.5. Arch. Numer. Softw. 3, 9–23 (2015)
Arnold, D.N., Brezzi, F., Fortin, M.: A stable finite element for the Stokes equations. Calcolo 21, 337–344 (1984)
Atri, A., Amundson, J., Clapham, D., Sneyd, J.: A single-pool model for intracellular calcium oscillations and waves in the Xenopus Laevis oocyte. Biophys. J. 65, 1727–1739 (1993)
Banerjee, S., Marchetti, M.C.: Instabilities and oscillations in isotropic active gels. Soft Matter 7, 463–473 (2011)
Barnafi, N., Gatica, G.N., Hurtado, D.E., Miranda, W., Ruiz-Baier, R.: New primal and dual-mixed finite element methods for stable image registration with singular regularization. Math. Models Methods Appl. Sci. 31, 979–1020 (2021)
Bausch, A.R., Ziemann, F., Boulbitch, A.A., Jacobson, K., Sackmann, E.: Local measurements of viscoelastic parameters of adherent cell surfaces by magnetic bead microrheometry. Biophys. J. 75, 2038–2049 (1998)
Berridge, M.J., Lipp, P., Bootman, M.D.: The versatility and universality of calcium signalling. Nat. Rev. Mol. Cell Biol. 1, 11–21 (2000)
Beysens, D.A., Forgacs, G., Glazier, J.A.: Embryonic tissues are viscoelastic materials. Can. J. Phys. 78, 243–251 (2000)
Boffi, D., Brezzi, F., Fortin, M.: Mixed Finite Element Methods and Applications. Springer Series in Computational Mathematics, vol. 44. Springer, Heidelberg (2013)
Braeu, F.A., Seitz, A., Aydin, R.C., Cyron, C.J.: Homogenized constrained mixture models for anisotropic volumetric growth and remodeling. Biomech. Model. Mechanobiol. 16, 889–906 (2017)
Brenner, S.C., Scott, L.R.: The Mathematical Theory of Finite Element Method, 3rd edn. Texts in Applied Mathematics, vol. 15. Springer, New York (2008)
Brinkmann, F., Mercker, M., Richter, T., Marciniak-Czochra, A.: Post-turing tissue pattern formation: Advent of mechanochemistry. PLos Comput. Biol. 14, e1006259 (2018)
Cherubini, C., Filippi, S., Gizzi, A., Ruiz-Baier, R.: A note on stress-driven anisotropic diffusion and its role in active deformable media. J. Theor. Biol. 430, 221–228 (2017)
Christodoulou, N., Skourides, P.A.: Cell-autonomous Ca2+ flashes elicit pulsed contractions of an apical actin network to drive apical constriction during neural tube closure. Cell Rep. 13, 2189–2202 (2015)
Cyron, C.J., Humphrey, J.D.: Growth and remodeling of load-bearing biological soft tissues. Meccanica 52, 645–664 (2017)
Cyron, C.J., Aydin, R.C., Humphrey, J.D.: A homogenized constrained mixture (and mechanical analog) model for growth and remodeling of soft tissue. Biomech. Model. Mechanobiol. 15, 1389–1403 (2016)
De Oliveira Vilaca, L.M., Gómez-Vargas, B., Kumar, S., Ruiz-Baier, R., Verma, N.: Stability analysis for a new model of multi-species convection-diffusion-reaction in poroelastic tissue. Appl. Math. Model. 84, 425–446 (2020)
De Oliveira Vilaca, L.M., Milinkovitch, M.C., Ruiz-Baier, R.: Numerical approximation of a 3D mechanochemical interface model for skin patterning. J. Comput. Phys. 384, 383–404 (2019)
Dillon, R., Othmer, H.G.: A mathematical model for outgrowth and spatial patterning of the vertebrate limb bud. J. Theor. Biol. 197, 295–330 (1999)
Dupont, G., Falcke, M., Kirk, V., Sneyd, J.: Models of Calcium Signalling. Interdisciplinary Applied Mathematics, vol. 43. Springer, Berlin (2016)
Friedl, P., Alexander, S.: Cancer invasion and the microenvironment: plasticity and reciprocity. Cell 147, 992–1009 (2011)
Ghosh, K., Pan, Z., Guan, E., Ge, S., Liu, Y., Nakamura, T., Ren, X.-D., Rafailovich, M., Clark, R.A.F.: Cell adaptation to a physiologically relevant ECM mimic with different viscoelastic properties. Biomaterials 28, 671–679 (2007)
Girault, V., Raviart, P.-A.: Finite Element Methods for Navier-Stokes Equations. Theory and Algorithms. Springer Series in Computational Mathematics, vol. 5. Springer, Berlin (1986)
Girault, V., Rivière, B., Wheeler, M.F.: A discontinuous Galerkin method with nonoverlapping domain decomposition for the Stokes and Navier-Stokes problems. Math. Comput. 74, 53–84 (2005)
Guiu-Souto, J., Muñuzuri, A.P.: Influence of oscillatory centrifugal forces on the mechanism of Turing pattern formation. Phys. Rev. E 91, 012917 (2015)
Hamill, O.P.: Twenty odd years of stretch-sensitive channels. Pflugers Arch – Eur. J. Physiol. 453, 333–351 (2006)
Hara, Y.: Contraction and elongation: Mechanics underlying cell boundary deformations in epithelial tissue. Dev. Growth Differ. 59, 340–350 (2017)
Herrmann, L.R.: Elasticity equations for incompressible and nearly incompressible materials by a variational theorem. AIAA J. 3, 1896–1900 (1965)
Javierre, E., Moreo, P., Doblaré, M., García-Aznar, J.M.: Numerical modeling of a mechano-chemical theory for wound contraction analysis. Int. J. Solids Struct. 46, 3597–3606 (2009)
Jones, G.W., Chapman, S.J.: Modeling growth in biological materials. SIAM Rev. 54, 52–118 (2012)
Kaouri, K., Maini, P.K., Skourides, P.A., Christodoulou, N., Chapman, S.J.: A simple mechanochemical model for calcium signalling in embryonic epithelial cells. J. Math. Biol. 78, 2059–2092 (2019)
Kim, Y., Hazar, M., Vijayraghavan, D.S., Song, J., Jackson, T.R., Joshi, S.D., Messner, W.C., Davidson, L.A., LeDuc, P.R.: Mechanochemical actuators of embryonic epithelial contractility. Proc. Nat. Acad. Sci. USA 111, 14366–14371 (2014)
Kobayashi, Y., Sanno, Y., Sakai, A., Sawabu, Y., Tsutsumi, M., Goto, M., Kitahata, H., Nakata, S., Kumamoto, U., Denda, M., Nagayama, M.: Mathematical modeling of calcium waves induced by mechanical stimulation in keratinocytes. PLoS ONE 9, e92650 (2014)
Kuchta, M., Mardal, K.-A., Mortensen, M.: Characterisation of the space of rigid motions in arbitrary domains. In: Proceedings of 8th National Conference on Computational Mechanics, Barcelona, Spain (2015)
Kuchta, M., Mardal, K. -A., Mortensen, M.: On the singular Neumann problem in linear elasticity. Numer. Linear Algebra Appl. 26, e2212 (2018)
Lång, E., Połeć, A., Lång, A., Valk, M., Blicher, P., Rowe, A.D., Tønseth, K.A., Jackson, C.J., Utheim, T.P., Janssen, L.M.C., Eriksson, J., Bøe, S.O.: Coordinated collective migration and asymmetric cell division in confluent human keratinocytes without wounding. Nat. Commun. 9, e3665 (2018)
Lecuit, T., Lenne, P.-F.: Cell surface mechanics and the control of cell shape, tissue patterns and morphogenesis. Nat. Rev. Mol. Cell. Biol. 8, 633–644 (2007)
Loppini, A., Gizzi, A., Ruiz-Baier, R., Cherubini, C., Fenton, F.H., Filippi, S.: Competing mechanisms of stress-assisted diffusivity and stretch-activated currents in cardiac electromechanics. Front. Physiol. 9, 1714 (2018)
Luu, O., David, R., Ninomiya, H., Winklbauer, R.: Large-scale mechanical properties of Xenopus embryonic epithelium. Proc. Nat. Acad. Sci. USA 108, 4000–4005 (2011)
Mercker, M., Brinkmann, F., Marciniak-Czochra, A., Richter, T.: Beyond turing: mechanochemical pattern formation in biological tissues. Biol. Direct 11, 22 (2016)
Moeendarbary, E., Valon, L., Fritzsche, M., Harris, A.R., Moulding, D.A., Thrasher, A.J., Stride, E., Mahadevan, L., Charras, G.T.: The cytoplasm of living cells behaves as a poroelastic material. Nat. Mater. 12, 253–261 (2013)
Moreo, P., Gaffney, E.A., García-Aznar, J.M., Doblaré, M.: On the modelling of biological patterns with mechanochemical models: Insights from analysis and computation. Bull. Math. Biol. 72, 400–431 (2010)
Murray, J.D.: Mathematical Biology II: Spatial Models and Biomedical Applications. Interdisciplinary Applied Mathematics, vol. 18. Springer, New York (2006)
Murray, J.D.: On the mechanochemical theory of biological pattern formation with application to vasculogenesis. C. R. Biol. 326, 239–252 (2003)
Murray, J.D., Oster, G.F.: Generation of biological pattern and form. IMA J. Math. Appl. Med. Biol. 1, 51–75 (1984)
Narciso, C.E., Contento, N.M., Storey, T.J., Hoelzle, D.J., Zartman, J.J.: Release of applied mechanical loading stimulates intercellular calcium waves in drosophila wing discs. Biophys. J. 113, 491–501 (2017)
Neville, A.A., Matthews, P.C., Byrne, H.M.: Interactions between pattern formation and domain growth. Bull. Math. Biol. 68, 1975–2003 (2006)
Nguyen, T.L., Polanco, E.R., Patananan, A.N., Zangle, T.A., Teitell, M.A.: Cell viscoelasticity is linked to fluctuations in cell biomass distributions. Sci. Rep. 10, 7403 (2020)
Ohayon, J., Tracqui, P.: Computation of adherent cell elasticity for critical cell-bead geometry in magnetic twisting experiments. Ann. Biomed. Eng. 33, 131–141 (2005)
Pandolfi, A., Gizzi, A., Vasta, M.: Visco-electro-elastic models of fiber-distributed active tissues. Meccanica 52, 3399–3415 (2017)
Preziosi, L., Ambrosi, D., Verdier, C.: An elasto-visco-plastic model of cell aggregates. J. Theor. Biol. 262, 35–47 (2010)
Propp, A., Gizzi, A., Levrero-Florencio, F., Ruiz-Baier, R.: An orthotropic electro-viscoelastic model for the heart with stress-assisted diffusion. Biomech. Model. Mechanobiol. 19, 633–659 (2020)
Quarteroni, A., Valli, A.: Numerical Approximation of Partial Differential Equations. Springer Series in Computional Mathematics, vol. 23. Springer, Berlin, Heidelberg (2008)
Radszuweit, M., Engel, H., Bär, M.: An active poroelastic model for mechanochemical patterns in protoplasmic droplets of physarum polycephalum. PLoS ONE 9, e99220 (2014)
Recho, P., Hallou, A., Hannezo, E.: Theory of mechanochemical pattering in biphasic biological tissues. Proc. Nat. Acad. Sci. USA 116, 5344–5349 (2019)
Ruiz-Baier, R.: Primal-mixed formulations for reaction–diffusion systems on deforming domains. J. Comput. Phys. 299, 320–338 (2015)
Ruiz-Baier, R., Gizzi, A., Rossi, S., Cherubini, C., Laadhari, A., Filippi, S., Quarteroni, A.: Mathematical modelling of active contraction in isolated cardiomyocytes. Math. Med. Biol. 31, 259–283 (2014)
Sanderson, M.J., Charles, A.C., Dirksen, E.R.: Mechanical stimulation and intercellular communication increases intracellular Ca2+ in epithelial cells. Cell Regul. 1, 585–596 (1990)
Suzuki, M., Sato, M., Koyama, H., Hara, Y., Hayashi, K., Yasue, N., Imamura, H., Fujimori, T., Nagai, T., Campbell, R.E., Ueno, N.: Distinct intracellular Ca2+ dynamics regulate apical constriction and differentially contribute to neural tube closure. Development 144, 1307–1316 (2017)
Vainio, I., Khamidakh, A.A., Paci, M., Skottman, H., Juuti-Uusitalo, K., Hyttinen, J., Nymark, S.: Computational model of Ca2+ wave propagation in human retinal pigment epithelial ARPE-19 cells. PLoS ONE 10, e0128434 (2015)
Weihs, D., Gefen, A., Vermolen, F.J.: Review on experiment-based two- and three-dimensional models for wound healing. Interface Focus 6, e20160038 (2016)
Wu, M., Ben Amar, M.: Growth and remodelling for profound circular wounds in skin. Biomech. Model. Mechanobiol. 14, 357–370 (2015)
Wyczalkowski, M.A., Chen, Z., Filas, B.A., Varner, V.D., Taber, L.A.: Computational models for mechanics of morphogenesis. Birth Defects Res. 92, 132–152 (2012)
Yamada, S., Wirtz, D., Kuo, S.C.: Mechanics of living cells measured by laser tracking microrheology. Biophys. J 78, 1736–1747 (2000)
Yao, W., Yang, H., Li, Y., Ding, G.: Dynamics of calcium signal and leukotriene c4 release in mast cells network induced by mechanical stimuli and modulated by interstitial fluid flow. Adv. Appl. Math Mech. 8, 67–81 (2016)
Acknowledgements
This work has been partially supported by the Australian Research Council through the Discovery Project Grant DP220103160; by the 730897 HPC-Europa3 Transnational Access Grant HPC175QA9K; by the Monash Mathematics Research Fund S05802-3951284; and by the Ministry of Science and Higher Education of the Russian Federation within the framework of state support for the creation and development of World-Class Research Centers “Digital biodesign and personalised healthcare” No. 075-15-2020-926.
Funding
Open Access funding enabled and organized by CAUL and its Member Institutions.
Author information
Authors and Affiliations
Corresponding author
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Dedicated to Alfio Quarteroni on the occasion of his 70th birthday.
Appendix: Effect of Number of Cells
Appendix: Effect of Number of Cells
Figure 11 illustrates periodic wavetrains appearing on a larger domain (disk of radius 25). This domain approximates the region of 20 neighbouring cells, and representing a cell compound of 1000μ m in diameter.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Kaouri, K., Méndez, P.E. & Ruiz-Baier, R. Mechanochemical Models for Calcium Waves in Embryonic Epithelia. Vietnam J. Math. 50, 947–975 (2022). https://doi.org/10.1007/s10013-022-00579-y
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s10013-022-00579-y
Keywords
- Viscoelasticity
- Advection-reaction-diffusion equations
- Calcium signalling
- Embryogenesis
- Excitability
- Finite element method