A new mechanochemical model for apical constriction: Coupling calcium signalling and viscoelasticity

Embryonic epithelial cells exhibit strong coupling of mechanical responses to chemical signals and most notably to calcium. Recent experiments have shown that the disruption of calcium signals during neurulation strongly correlates with the appearance of neural tube defects. We, thus, develop a multi-dimensional mechanochemical model and use it to reproduce important experimental findings that describe anterior neural plate morphogenetic behaviour during neural tube closure. The governing equations consist of an advection-diffusion-reaction system for calcium concentration which is coupled to a force balance equation for the tissue. The tissue is modelled as a linear viscoelastic material that includes a calcium-dependent contraction stress. We implement a random distribution of calcium sparks that is compatible with experimental findings. A finite element method is employed to generate numerical solutions of the model for an appropriately chosen range of parameter values. We analyse the behaviour of the model as three parameters vary: the level of IP₃ concentration, the strength of the stretch-sensitive activation and the maximum magnitude of the calcium-dependent contraction stress. Importantly, the simulations reproduce important experimental features, such as the spatio-temporal correlation between calcium transients and tissue deformation, the monotonic reduction of the apical surface area and the constant constriction rate, as time progresses. The model could also be employed to gain insights into other biological processes where the coupling of calcium signalling and mechanics is important, such as carcinogenesis and wound healing.

1 Introduction

During the early stages of the development of an embryo’s central nervous system (CNS), neuroepithelial cells undergo a shape change via apical constriction (AC), a morphogenetic process controlled by apical actomyosin contraction that is induced by calcium transients (Christodoulou and Skourides, 2015; Suzuki et al., 2017). AC results in the folding of the neural plate and in the formation of the neural tube. It is not fully understood how AC in the neural plate is controlled and how it contributes to tissue morphogenesis but recent experiments have shown that Ca²⁺ plays a crucial role in regulating AC during neural tube closure (NTC) (Christodoulou and Skourides, 2015; Suzuki et al., 2017). Furthermore, pharmacologically inhibiting Ca²⁺ has been shown to lead to neural tube defects (Smedley and Stanisstreet, 1986; Wallingford et al., 2001; Christodoulou and Skourides, 2015), such as Spina Bifida and anencephaly.

Many experiments have documented that intracellular Ca²⁺ release triggers actomyosin-based contractions, in both embryonic and cultured cells (Wallingford et al., 2001; Herrgen et al., 2014; Hunter et al., 2014; Christodoulou and Skourides, 2015; Suzuki et al., 2017). The ability of cells to sense and respond to forces by elevating their cytosolic Ca²⁺ is also well established; mechanically stimulated Ca²⁺ waves have been observed propagating through ciliated tracheal epithelial cells (Sanderson and Sleigh, 1981; Sanderson et al., 1988; Sanderson et al., 1990), rat brain glial cells (Charles et al., 1991; Charles et al., 1992; Charles et al., 1993), developing epithelial cells in Drosophila wing discs (Narciso et al., 2017) and many other cell types (Young et al., 1999; Bereiter-Hahn, 2005; Tsutsumi et al., 2009; Yang et al., 2009). Thus, different types of mechanical stimuli, from shear stress to direct mechanical stimulation, can elicit Ca²⁺ elevation (although the sensing mechanism may differ in each case). Moreover, localisation of stresses or strains within the cells can generate alteration in patterns of Ca²⁺ distribution in a tissue by changing cell displacement magnitude, direction, and velocity (Lecuit and Lenne, 2007; Guiu-Souto and Munuzuri, 2015). This is especially noteworthy since distinct Ca²⁺ signalling patterns differentially modulate AC for efficient epithelial folding. The latter mechanism has a broad range of physiological outcomes (Suzuki et al., 2017).

Since mechanical stimulation elicits Ca²⁺ release and Ca²⁺ elicits contractions which are sensed as mechanical stimuli by the cell, a two-way mechanochemical feedback between Ca²⁺ and contractions should be at play. Motivated by the recent experimental observations (Christodoulou and Skourides, 2015; Suzuki et al., 2017) where, during AC, increasing tension in the contracting cells yields Ca²⁺ release which, in turn, elicits contractions in the cells which are sensed as mechanical stimuli by the neighbouring cells, we develop a new mechanochemical model that captures the interplay of Ca²⁺ signalling and mechanical forces during AC.

This paper extends the mechanochemical model in (Kaouri et al., 2019), which describes the coupling of Ca²⁺ signalling with the mechanics of the embryonic epithelial tissue during AC in one spatial dimension; it also extends the multi-dimensional model presented in (Kaouri et al., 2022). In the aforementioned models, following the early mechanochemical models in (Murray and Oster, 1984), where small strains are assumed, the embryonic tissue is assumed to be a linear viscoelastic (Kelvin–Voigt) solid (with one elastic spring and two viscous dashpots), where only after the initial stress has vanished, does the material go back to its original state. Also, In the model proposed here, as in (Murray and Oster, 1984; Banerjee and Marchetti, 2011; Kaouri et al., 2019; Kaouri et al., 2022), we assume that the viscoelastic stress includes an active contraction stress which depends on the cytosolic Ca²⁺ concentration. The models in (Kaouri et al., 2019; Kaouri et al., 2022), as well as the model presented here, employ the well-established Ca²⁺ signalling model from Atri et al. (1993), called the “Atri model” hereafter. The Atri model captures the experimentally verified Ca²⁺-induced Ca²⁺ release (CICR) process. It consists of two differential equations, one PDE for the cytosolic Ca²⁺ concentration and another PDE for the percentage of the non-inactivated IP₃ receptors on the endoplasmic reticulum (ER) which allow release of Ca²⁺ from the ER into the cytosol.

In Figure 1 we show still images from a time lapse recording of the anterior neural plate during the last stage of neural tube closure (stage 16 of Xenopus embryo development). For live imaging, 4-cell stage Xenopus laevis embryos were injected with the mRNA encoding membrane-GFP and the calcium sensor GECO-RED at the two dorsal blastomeres to target the neural tissue. Subsequently the embryos were allowed to develop until stage 14 and imaged during neural tube closure. The time lapse recordings of neural tube closure that we represent in Figure 1 were generated on a ZEISS LSM 710 confocal microscope with a 30 s time interval. At this stage, the ectoderm of the embryo consists of the neuroepithelium, which is surrounded by the surface ectoderm. The last stage of anterior neural tube closure is controlled by AC of neuroepithelial cells (Christodoulou and Skourides, 2015) and lasts about 40 min; this is the stage we model here. During AC cells reduce their apical surface area and change their shape from columnar to a wedge shape (Christodoulou and Skourides, 2015). These cell shape changes subsequently drive the bending of the neuroepithelium and the formation of the neural tube. Note that the frequency of calcium transients has been quantified in (Christodoulou and Skourides, 2015). The observation there is that the frequency increases as neural tube closure progresses. This information ties up with the data presented also in (Kaouri et al., 2019). Thus, experimental evidence shows a clear correlation between the appearance of calcium transients, and the reduction of the apical surface area during neural tube closure. Even though the process of AC is three-dimensional we focus attention to the stage where the apical surface area reduces (in a ratchet-like manner). This is the active driver of the process and it is sufficient to describe it with a two-dimensional model, as we do below. For stage 16 of Xenopus embryo development studied here we can assume small strains; hence the tissue can be modelled as a linear Kelvin–Voigt viscoelastic solid. This linear viscoelastic material is 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 (Nguyen et al., 2020). In the present mechanochemical model, we assume that the viscoelastic stress includes a contraction stress which depends on Ca²⁺ concentration, following the formulation in (Murray and Oster, 1984; Banerjee and Marchetti, 2011; Kaouri et al., 2019; Kaouri et al., 2022).

FIGURE 1

FIGURE 1. Snapshots of the closure of the neural tube over a period of 55 min (from top-left to bottom-right). The green colour indicates the cell membrane and magenta is the calcium sensor.

The new mechanochemical model proposed here, following (Kaouri et al., 2022) is underpinned also by the following fundamental assumptions: a) the equilibrium of the mechanics 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 (Herrmann, 1965)]. The introduction of solid pressure contributes to achieve robustness in the nearly incompressible regime assumed for the tissue. This occurs when the Poisson ratio approaches 0.5, implying that the first Lamé parameter defining the dilation properties of the material is very large. Also, the mechanochemical coupling is modelled directly in the viscoelastic stress through a Hill function that depends on Ca²⁺ and through the modification of the reaction kinetics by volume change. The two-way coupling mechanism we adopt follows the model structure used in (Murray, 2001; Murray, 2003; Neville et al., 2006; Ruiz-Baier et al., 2014; Kaouri et al., 2019; Kaouri et al., 2022).

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, via an implicit, fully coupled finite element method (Ruiz-Baier et al., 2014; Kaouri et al., 2022). Following (Kaouri et al., 2022), we nondimensionalise the model using experimentally verified parameter values from neuroepithelial cells undergoing AC during NTC (Atri et al., 1993; D’Angelo et al., 2019; Benko and Brodland, 2007) to investigate whether our model reproduces important features of NTC observed in the experiments of (Christodoulou and Skourides, 2015).

This paper is organised as follows. In Section 2 we present a new mechanochemical model capturing the coupling of calcium signalling to forces in a deforming embryonic epithelial tissue undergoing AC. We also present the computational implementation of the model, using a Finite Element Method. Next, in Section 3 we present the simulations and discuss how they reproduce important experimental features. Finally, Section 4 includes our conclusions and future research directions.

2 Methods

2.1 A new mechanochemical model for apical constriction, coupling calcium signalling and mechanics

Here, we present a new mechanochemical model coupling calcium and mechanics in AC. We adapt, as previously, the Atri et al. (1993) model to write the governing equations for the cytosolic calcium concentration and the percentage of non-inactivated IPR—for more details on this well-established model see (Atri et al., 1993; Kaouri et al., 2019; Kaouri et al., 2022). We also assume, as previously, that the tissue is represented by a linear viscoelastic, Kelvin–Voigt material (Murray and Oster, 1984; Kaouri et al., 2019; Kaouri et al., 2022). The system is assumed to be in mechanical equilibrium, that is the contraction forces generated by the calcium are in mechanical equilibrium with the viscoelastic forces. The model is as follows:

${\partial }_{t}c+{\partial }_t\mathit{u}\cdot \nabla c-D{\nabla }^{2}c=f\left(c,h,\mathit{u}\right)\text{in\hspace{0.17em}}\mathrm{\Omega }\times \left(0,t_{\text{final}}\right],(1)$

${\partial }_{t}h+{\partial }_t\mathit{u}\cdot \nabla h=g\left(c,h\right)\text{in\hspace{0.17em}}\mathrm{\Omega }\times \left(0,t_{\text{final}}\right],(2)$

$-\mathrm{div}\left(\mathit{ϵ}\left(\mathit{u}\right)-p\mathbf{I}+{\alpha }_1{\partial }_t\mathit{ϵ}\left(\mathit{u}\right)-{\alpha }_2{\partial }_{t}p\mathbf{I}+T\left(c\right)\mathbf{I}\right)=\mathbf{0}\text{in\hspace{0.17em}}\mathrm{\Omega }\times \left(0,t_{\text{final}}\right],(3)$

$p+\frac{\nu }{\left(1-2\nu \right)}\mathrm{div}\mathit{u}=0\text{in\hspace{0.17em}}\mathrm{\Omega }\times \left(0,t_{\text{final}}\right],(4)$

where [Ca²⁺] = c is the cytosolic calcium concentration, h is the percentage of non-inactivated IPR, u is the tissue displacement and p is the Herrmann pressure. ν is the Poisson’s ratio, α₁ and α₂ are the shear and bulk viscosities, respectively, and D is the diffusion coefficient of cytosolic calcium. The Cauchy stress has elastic, viscous, and active calcium-dependent stress components. The active stress and the reaction kinetics are specified as follows:

$T\left(c\right)={\beta }_1\frac{c^n}{{\beta }_2+c^n},f\left(c,h,\mathit{u}\right)=I_{\text{rand}}\left(\mathit{x},t\right)\mu h{K}_1\frac{b+c}{1+c}-\frac{Gc}{K+c}+\lambda \mathrm{div}\mathit{u},g\left(c,h\right)=\frac{1}{1+c^2}-h.(5)$

The function I_rand multiplying the CICR Ca²⁺ flux is a random-in-space distribution of Ca²⁺ sparks which increases in frequency and in amplitude with time, observed in the experiments of (Christodoulou and Skourides, 2015). Since there is a (thin) circumferential layer of epidermal cells surrounding the neuroepithelial cells, the Young’s modulus, E, is discontinuous across the interface of these two regions (Wiebe and Brodland, 2005). Hence, we assume that the Young’s modulus in the domain is given by

$E=E_{\text{in}}{\chi }_{{\mathrm{\Omega }}_{\text{in}}}+E_{\text{out}}{\chi }_{{\mathrm{\Omega }}_{\text{out}}},(6)$

where χ_M denotes the characteristic function on the generic subdomain M, and by E_in and E_out we denote the Young modulus in the inner and outer regions, respectively. The model parameters and their values are discussed in more detail in Section 2.2.

The PDE system (2.1) is complemented with appropriate initial data for c, h, u and p, respectively given by

$c\left(0\right)=c_0,h\left(0\right)=h_0=\frac{1}{1+c_0^2},\mathit{u}\left(0\right)=\mathbf{0},p\left(0\right)=0,\text{in\hspace{0.17em}}\mathrm{\Omega },(7)$

where c₀ is the steady state value of c. We also assume stress-free and zero-flux boundary conditions on the domain boundary, as follows:

$\left(\mathit{ϵ}\left(\mathit{u}\right)-p\mathbf{I}+{\alpha }_1{\partial }_t\mathit{ϵ}\left(\mathit{u}\right)-{\alpha }_2{\partial }_{t}p\mathbf{I}+T\left(c\right)\mathbf{I}\right)\mathit{n}=\mathbf{0}\text{and}D\nabla c\cdot \mathit{n}=0\text{on\hspace{0.17em}}\partial \mathrm{\Omega }\times \left(0,t_{\text{final}}\right].(8)$

These pure-traction boundary conditions necessitate imposing an additional condition to render the system well-defined. We, hence, impose that the displacements are orthogonal with respect to the space of rigid motions, that is

$\mathbb{R}\mathbb{M}\left(\mathrm{\Omega }\right)≔\left\{\mathit{v}\in {\mathbf{H}}^1\left(\mathrm{\Omega }\right):\mathit{ϵ}\left(\mathit{v}\right)=\mathbf{0}\right\}.(9)$

Note that T(c) in Eq. 3 has an opposite sign to that in the models of (Kaouri et al., 2019; Kaouri et al., 2022). There, the opposite sign corresponded to dilation instead of contraction here [see also (Murray and Oster, 1984; Moreo et al., 2010)].

2.2 Model parameter values

The Atri Ca²⁺ signalling model we use (Atri et al., 1993) captures the Ca²⁺ release to the cytosol via the IPR/Ca²⁺ channels, relying on experimental data from the Xenopus laevis oocyte. We nondimensionalised this model in detail in (Kaouri et al., 2019); here we present it in its nondimensional form, choosing the same parameter values.

The values of the mechanical parameters were taken from Xenopus and Drosophila embryos (D’Angelo et al., 2019; Benko and Brodland, 2007; Zhou et al., 2009; Wiebe and Brodland, 2005) and are collected in Table 1. (The parameter values can be taken from two different species since the magnitude of sub-cellular forces are similar across species.) To determine the Young’s modulus and viscosity of neural tissues, in (Zhou et al., 2009) the authors measured the stiffness of dorsal isolate explants of Xenopus laevis embryos over different stages of development, from gastrulation to neurulation. They recorded the values of Young’s modulus and viscosity over five dorsal isolate explants taken from stage 16 embryos. We averaged these five values to obtain E_in. To determine the Young’s modulus for the epidermal cells surrounding the neuroepithelium, E_out, we use the ratio between the stiffness moduli of the neuroepithelium and epidermis, as determined in (Wiebe and Brodland, 2005); we, thus, estimate E_out = 0.55E_in.

TABLE 1

TABLE 1. Parameter values for the mechanochemical model.

We determine the shear and bulk viscosities, ${\tilde{\alpha }}_1$ and ${\tilde{\alpha }}_2$ using data from (Zhou et al., 2009) and (D’Angelo et al., 2019). In (Zhou et al., 2009), measurements were taken from five neurulating embryos and we determined the value of viscosity to be the average of the five values. This viscosity value was then split using the ratio between the shear and bulk viscosity given in (D’Angelo et al., 2019).

For the Poisson’s ratio, we assume, as in (Zhou et al., 2009), that the embryonic tissue is a nearly incompressible material and hence we set ν = 0.4. This value is also consistent with the range of values in (D’Angelo et al., 2019) and in other experimental studies, e.g., (Zhou et al., 2015). The value of the maximum (saturation) traction stress, T₀, is difficult to determine but experiments on zebrafish primordium tissue suggest that the value can range from 50 Pa to 450 Pa (Yamaguchi et al., 2022). This range is supported by (Kraning-Rush et al., 2012) and (Oakes et al., 2014) where traction force microscopy revealed the average traction stresses of cells on 2D substrates to be between 100 Pa and 1,000 Pa. In our model, T₀ was set by tuning its value while keeping all other parameter values constant.

The area of a single neuroepithelial cell at the start of apical constriction is approximately $250\mu {\mathrm{m}}^2$. We have 256 cells (Suzuki et al., 2017) tightly packed in the tissue and, hence, the initial tissue area is approximately 64,000μm². We assume that the tissue is a disc-shaped domain of radius $\approx 143\mu \mathrm{m}$. The spatial variables have been non-dimensionalised using L = 100 μm. Thus, the tissue is represented as a disc of radius R = 1.43, in non-dimensional terms.

We take the non-dimensional parameters from (Kaouri et al., 2019), as follows: D = 0.004, K₁ = 46.28, G = 5.71, and K = 0.14. The three parameters we are going to vary in the simulations are

$\mu ,\lambda \mathrm{a}\mathrm{n}\mathrm{d}{\beta }_1=\frac{T_0\left(1+\nu \right)}{E}.(10)$

As the Young modulus, E, is discontinuous across the interface of the neuroepithelium and the epidermis, so are the parameters α₁, α₂ and β₁. From the nondimensionalisation it arises that ${\alpha }_1=\frac{(1+\nu )}{E{\tau }_j}{\tilde{\alpha }}_1$ and ${\alpha }_2=\frac{(1+\nu )(1-2\nu )}{E\nu {\tau }_j}{\tilde{\alpha }}_2$, where τ_j = 2s and other values as in Table 1. On the neuroepithelium we, hence, have ${\alpha }_1^{\text{in}}=59.94$ and ${\alpha }_2^{\text{in}}=4.35$, whereas on the epidermis we take ${\alpha }_1^{\text{out}}=113.67$ and ${\alpha }_2^{\text{out}}=8.25$.

2.3 Computational implementation of the model using the finite element method

The mechanochemical model Eqs 1–4 has been discretised using the finite element method (FEM). The open-source FEM library FEniCS (Langtangen and Logg, 2017) was used to obtain the numerical approximation of the variational formulation of the governing equations. Due to the nonlinear nature of the model, the Newton–Raphson method was used and at each iteration the linear tangent system was solved with the MUMPS direct solver (Amestoy et al., 2000). Time derivatives were approximated by a fully implicit backward differencing scheme. A mixed finite element formulation based on the MINI element (Cioncolini and Boffi, 2019) was used for the numerical approximation of the displacement and the rescaled Herrmann pressure, and piecewise linear and overall continuous elements for the calcium concentration and IPR. The rigid motions in a finite-dimensional subspace of Eq. 9 were removed from the set of admissible displacements using a Lagrange multiplier approach, as described in (Kaouri et al., 2022). For the Newton–Raphson iterative algorithm, a tolerance of 10^–8 was used. Note that, as long as an appropriate discrete inf-sup condition is satisfied and one can still construct a suitable auxiliary discrete problem to remove rigid motion solutions, high-order elements can also be used, for example generalised Taylor–Hood elements of degree k ≥ 2 for the approximation of displacement and rescaled Herrmann pressure, and piecewise polynomials of order k for the approximation of calcium and IPR concentrations. For the lowest-order Taylor–Hood elements one has an additional order of convergence with respect to the MINI-element. That is, we expect an improvement in model predictions as the mesh is refined, but at the price of solving a larger system at each Newton–Raphson iteration. More details about the code can be found in (Ruiz-Baier et al., 2014; Kaouri et al., 2022).

3 Results and discussion

In this section we present numerical solutions of the model for a range of parameter values and explore the agreement of the model with experimental results. We treat μ, λ, and β₁ as bifurcation parameters and identify the set of values for which the model exhibits agreement with the experimental results in (Christodoulou and Skourides, 2015). The computational domain is a disk of radius R = 1.43, discretised into an unstructured triangular mesh of 34,947 elements. A fixed time-step of Δt = 0.1 is used in all simulations.

In order to generate the random field of calcium sparks, I_rand, we impose a frequency that is linearly increasing from 0.1 to 0.4 and set the amplitude to 1 + a_mpl, with a_mpl increasing from 0.47 to 0.78 quadratically–see Figure 2. A single spark at the domain centre is included in all simulations.

FIGURE 2

FIGURE 2. Frequency and normalised amplitude of calcium sparks versus time used to construct I_rand, fitted to experimental data from (Kaouri et al., 2019).

We proceed to extract transients of the model variables c, h, p, u at 80 points which are located in a square of side 0.25, centred at the origin. We then average these values over space to generate the evolution of the system over time. The simulations are presented in Figures 3, 4. Figure 3 shows the results obtained when ${\beta }_1^{\text{in}}=T_0(1+\nu )/E_{\text{in}}=3.16$, that is when T₀ = 100 Pa. When μ = 0.288 the Atri model (Atri et al., 1993; Kaouri et al., 2019) does not exhibit oscillatory behaviour. The increase of the Herrmann pressure (and of the displacement) is monotonic, which indicates a monotonic contraction and area reduction, as observed in experiments (Christodoulou and Skourides, 2015; Suzuki et al., 2017). Hence, our model reproduces this key experimental feature. This behaviour occurs because the random-in-space calcium sparks (modelled by I_rand) exist elsewhere in the domain and increase in amplitude and frequency as time progresses (see Figure 2).

FIGURE 3

FIGURE 3. Plot of statistics (quartiles, ranges, and average) of the model variables over time, spatially averaged over 80 points near the disk centre (square of side 0.25). The scales for calcium concentration and for the percentage of open IPR are on the right axes, whereas the scales for the displacement magnitude and for the Herrmann pressure are on the left axes. Parameters are ${\beta }_1^{\text{in}}=3.16$, μ = 0.288 and λ ∈ {0.01, 0.1, 0.5, 1.3}, varied from top left to bottom right.

FIGURE 4

FIGURE 4. Plot of statistics (quartiles, ranges, and average) of the model variables over time, spatially averaged over 80 points near the disk centre (square of side 0.25). The scales for calcium concentration and for the percentage of open IPR are on the right axes, whereas the scales for the displacement magnitude and for the Herrmann pressure are on the left axes. Parameters are ${\beta }_1^{\text{in}}=7.91$ (T₀ = 250 Pa), μ = 0.288 and λ ∈ {0.01, 0.1, 0.5, 1.3}, varied from top left to bottom right panels.

Increasing T₀ to 250 Pa gives ${\beta }_1^{\text{in}}=7.91$. In this case we see in Figure 4 that the pressure and displacement approximately double in magnitude compared to those for T₀ = 100 Pa $({\beta }_1^{\text{in}}=3.16)$.

In Figures 3, 4, as time advances the Herrmann pressure increases, the tissue contracts and the area decreases monotonically. However, for any fixed value of μ, the contraction decreases as λ increases. Since λ is a measure of the strength of the coupling between the calcium signalling system and the mechanics of the tissue this result indicates that the stronger the coupling the smaller the contraction.

In Figure 5 we plot the tissue area as time progresses. We compute the area using the expression ∫_Ω det(I + ∇u)dx, and compare it with the experimental data in Figure 2. In the top left plot, for ${\beta }_1^{\text{in}}=7.91$ we plot the area reduction as μ and λ vary. We get a good fit with the area data in Figure 2 for μ = 0.288 and λ = 0.5 (yellow line). For the chosen parameters μ = 0.288, λ = 0.5, in the top right plot, we determine the area for ${\beta }_1^{\text{in}}\in =1.58,3.16,4.74,7.91$ (corresponding to T₀ = 50, 100, 150, 250 Pa and changing accordingly ${\beta }_1^{\text{out}}$); we see that the area reduction is quite sensitive to the choice of the active contractile force parameter, T₀, which confirms the nonlinear nature of the model. In the bottom plot the red dash-dotted curve depicts the constriction rate (rate of area reduction) for the parameters $\mu =0.288,\lambda =0.5,{\beta }_1^{\text{in}}=7.91$. The constriction rate is approximately constant, as identified in the experiments of (Christodoulou and Skourides, 2015).

FIGURE 5

FIGURE 5. Area (in non-dimensional units) with respect to time for different parameter values, plotted against experimental data from (Kaouri et al., 2019). The red, dash-dotted curve in the bottom plot depicts the approximately constant constriction rate, for the parameter set $\mu =0.288,\lambda =0.5,{\beta }_1^{\text{in}}=7.91$.

In Figure 6 we visualise the deformation of the tissue (disc domain) and the associated calcium distribution, at different times. The boundary of the initially non-deformed disc is also shown, for comparison. For λ = 0.01, we observe nucleation of calcium waves—synchronous waves that are sustained for a longer time. In Figure 6 we clearly visualise what has been already noted in Figure 4: that, as time advances the area always decreases monotonically and that the larger λ is the smaller the contraction. In Figure 7, for the same set of parameters, μ = 0.288, ${\beta }_1^{\text{in}}=7.91$, and varying λ, we plot all field variables at time t = 35 (to show a different time snapshot than the ones depicted before). For all cases, a larger displacement is observed near the boundary.

FIGURE 6

FIGURE 6. Snapshots of the contracting domain and the associated calcium distribution, at t = 20, 30, 40. Parameters are μ = 0.288, and λ = {0.01, 0.1, 0.5, 1.3} from the top to the bottom row, respectively. Here we use ${\beta }_1^{\text{in}}=7.91$.

FIGURE 7

FIGURE 7. Snapshots of the contracting domain and the associated calcium distribution (left), Herrmann pressure (centre), and displacement (right), at time t = 35. Parameters are: μ = 0.288, λ = 0.01 (top), λ = 0.1 (second row), λ = 0.5 (third row), λ = 1.3 (bottom). We use ${\beta }_1^{\text{in}}=7.91$.

Testing the code for different viscoelasticity parameters

Although real tissues usually show intermediate degrees of viscoelasticity (von Dassow and Davidson, 2011; Parvini et al., 2022; Brückner and Janshoff, 2015; Lange and Fabry, 2013; Karcher et al., 2003) it is not uncommon that mechanical models of morphogenesis assume one of two extremes in material behaviour: either a purely elastic or a purely viscous fluid. For the specific case of embryonic tissues, in (Davidson, 2011) authors conclude that most embryonic tissues should be considered viscoelastic in order to fully understand the mechanisms behind deformation of a multicellular tissue in response to any force or stress. Along similar lines, in (von Dassow and Davidson, 2011) it is noted that the stability and robustness of specific physical mechanisms of morphogenesis will likely have a strong dependence on the viscoelasticity of the tissue.

Here, for the stage of embryogenesis we consider (stage 16), we can assume small strains; hence we modelled the tissue as a linear Kelvin–Voigt viscoelastic material. However, we emphasize that the formulation we propose along with the finite element method we employ do accommodate for the pure elastic case and also for material constants close to the incompressibility limit. This is tested in the following simple example where we choose the parameters μ = 0.288, λ = 0.5, T₀ = 250, and take a higher Poisson ratio (= 0.4999) with or without shear-bulk viscosities. The results are visualised in Figure 8. They need to be compared with the base-line case shown in Figure 4 (bottom left panel). Both panels use ν = 0.4999 (which corresponds to the slightly higher ${\beta }_1^{\text{in}}=8.47$, since β₁ is proportional to the Poisson ratio). The left panel shows the behaviour when maintaining the base-line viscoelastic parameters. The displacement and pressure exhibit an initial peak and then return to a plateau phase. On the right panel we focus on the case with ${\tilde{\alpha }}_1={\tilde{\alpha }}_2=0$ (on both neuroepithelium and epidermis regions). Both calcium and the percentage of open IPR are quite similar to the base-line case. As in the left panel, both the Hermann pressure and the displacement exhibit an initial peak followed by an undershoot and then reach a plateau phase. However, the displacement magnitude is much lower than in the viscoelastic case (in the base-line case, both mechanical fields were increasing monotonically).We also note that in order to properly capture other stages of NTC, we would require a large-strain viscoelasticity framework in combination with a remodelling approach. This constitutes a direction for future research.

FIGURE 8

FIGURE 8. Testing behaviour in the nearly incompressible purely elastic and viscoelastic cases. Parameters are ν = 0.4999, ${\tilde{\alpha }}_1=3790,{\tilde{\alpha }}_2=550$ (left); and ν = 0.4999, ${\tilde{\alpha }}_1={\tilde{\alpha }}_2=0$ (right).

4 Conclusion

We propose a new mechanochemical model that reproduces important experimental findings on the apical constriction (AC) during the last stage of neural tube closure (NTC). AC is controlled by the complex coupling of calcium signalling to the mechanics of the embryonic epithelial tissue; disruption of calcium signals and consequently of AC leads to significant embryo malformations such as Spina Bifida and anencephaly.

The model builds on other recent mechanochemical models (Kaouri et al., 2019; Kaouri et al., 2022). The calcium-induced calcium release process allowing calcium to get released from the ER into the cytosol has been modulated with a random-in-space distribution of calcium sparks of which the amplitude and frequency increase with time. The distribution has been fitted to agree with experimental data presented in (Christodoulou and Skourides, 2015; Kaouri et al., 2019)–see Figures 1, 2. The embryonic tissue is modelled as a linear viscoelastic material, including three types of stresses: viscous, elastic and an active contraction stress which increases with calcium concentration until it saturates.

We have simulated the model using a finite element method on a disc domain packing 256 epithelial cells—details about the numerical method can be found in (Kaouri et al., 2022). We have studied the behaviour of the model as three parameters vary: μ, the level of IP₃, λ, which measures the strength of the stretch-sensitive activation and ${\beta }_1^{\text{in}}$ which represents the maximum contraction stress.

The model shows that for any value of μ, λ, and ${\beta }_1^{\text{in}}$ the tissue area is decreasing monotonically over time, as observed in experiments (Christodoulou and Skourides, 2015; Suzuki et al., 2017). Furthermore, we have identified that for μ = 0.288, λ = 0.5, and ${\beta }_1^{\text{in}}=7.91$, the monotonic area decrease fits to the experimental curve, generated in (Christodoulou and Skourides, 2015). Also in Figure 5 (bottom plot) we have quantified and plotted the constriction rate (red, dash-dotted curve) which is approximately constant as observed in experiments (Christodoulou and Skourides, 2015).

We also found that as λ increases the contraction effect decreases–see Figures 4, 6. This result could be tested in future experiments.

Data availability statement

The codes used to generate the Figures are provided on GitHub: https://github.com/ruizbaier/viscoelasticCalciumEpithelial

Author contributions

All authors were involved in the conceptualisation and design of the mathematical model, the computational implementation, and in the analysis and interpretation of results. All authors have read, revised and approved the final manuscript.

Funding

This work has been supported by the European Regional Development Fund and the Republic of Cyprus through the Research and Innovation Foundation (POST-DOC/0718/0087); by a Marie Skłodowska-Curie individual fellowship grant (101038073); by the Monash Mathematics Research Fund S05802-3951284; by the Australian Research Council through the Future Fellowship grant FT220100496 and Discovery Project grant DP22010316; 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-2022-304.

Acknowledgments

We thank Prof. Lance Davidson, for fruitful discussions, in particular regarding choosing appropriate parameter values. We also thank Prof. Tim N. Phillips for useful discussions on the modelling of materials.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Amestoy, P. R., Duff, I. S., L’Excellent, J. Y., and Koster, J. (2000). Mumps: A general purpose distributed memory sparse solver. Int. Workshop Appl. Parallel Comput., 121–130. doi:10.1007/3-540-70734-4_16

CrossRef Full Text | Google Scholar

Atri, A., Amundson, J., Clapham, D., and Sneyd, J. (1993). A single-pool model for intracellular calcium oscillations and waves in the xenopus laevis oocyte. Biophys. J. 65, 1727–1739. doi:10.1016/S0006-3495(93)81191-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Banerjee, S., and Marchetti, M. C. (2011). Instabilities and oscillations in isotropic active gels. Soft Matter 7, 463–473. doi:10.1039/c0sm00494d

CrossRef Full Text | Google Scholar

Benko, R., and Brodland, G. W. (2007). Measurement of in vivo stress resultants in neurulation-stage amphibian embryos. Ann. Biomed. Eng. 35, 672–681. doi:10.1007/s10439-006-9250-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Bereiter-Hahn, J. (2005). Mechanics of crawling cells. Med. Eng. Phys. 27, 743–753. doi:10.1016/j.medengphy.2005.04.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Brückner, B. R., and Janshoff, A. (2015). Elastic properties of epithelial cells probed by atomic force microscopy. Biochim. Biophys. Acta 1853, 3075–3082. doi:10.1016/j.bbamcr.2015.07.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Charles, A. C., Dirksen, E. R., Merrill, J. E., and Sanderson, M. J. (1993). Mechanisms of intercellular calcium signaling in glial cells studied with dantrolene and thapsigargin. Glia 7, 134–145. doi:10.1002/glia.440070203

PubMed Abstract | CrossRef Full Text | Google Scholar

Charles, A. C., Merrill, J. E., Dirksen, E. R., and Sandersont, M. J. (1991). Intercellular signaling in glial cells: Calcium waves and oscillations in response to mechanical stimulation and glutamate. Neuron 6, 983–992. doi:10.1016/0896-6273(91)90238-u

PubMed Abstract | CrossRef Full Text | Google Scholar

Charles, A. C., Naus, C., Zhu, D., Kidder, G. M., Dirksen, E. R., and Sanderson, M. J. (1992). Intercellular calcium signaling via gap junctions in glioma cells. J. Cell Biol. 118, 195–201. doi:10.1083/jcb.118.1.195

PubMed Abstract | CrossRef Full Text | Google Scholar

Christodoulou, N., and Skourides, P. A. (2015). 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. doi:10.1016/j.celrep.2015.11.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Cioncolini, A., and Boffi, D. (2019). The MINI mixed finite element for the Stokes problem: An experimental investigation. Comput. Math. Appl. 77, 2432–2446. doi:10.1016/j.camwa.2018.12.028

CrossRef Full Text | Google Scholar

D’Angelo, A., Dierkes, K., Carolis, C., Salbreux, G., and Solon, J. (2019). In vivo force application reveals a fast tissue softening and external friction increase during early embryogenesis. Curr. Biol. 29, 1564–1571. doi:10.1016/j.cub.2019.04.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Davidson, L. A. (2011). Embryo mechanics. Curr. Top. Dev. Biol., 215–241.

PubMed Abstract | CrossRef Full Text | Google Scholar

Guiu-Souto, J., and Munuzuri, A. P. (2015). Influence of oscillatory centrifugal forces on the mechanism of Turing pattern formation. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 91, 012917. doi:10.1103/PhysRevE.91.012917

PubMed Abstract | CrossRef Full Text | Google Scholar

Herrgen, L., Voss, O. P., and Akerman, C. J. (2014). Calcium-dependent neuroepithelial contractions expel damaged cells from the developing brain. Dev. Cell 31, 599–613. doi:10.1016/j.devcel.2014.10.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Herrmann, L. R. (1965). Elasticity equations for incompressible and nearly incompressible materials by a variational theorem. AIAA J. 3, 1896–1900. doi:10.2514/3.3277

CrossRef Full Text | Google Scholar

Hunter, G. L., Crawford, J. M., Genkins, J. Z., and Kiehart, D. P. (2014). Ion channels contribute to the regulation of cell sheet forces during Drosophila dorsal closure. Development 141, 325–334. doi:10.1242/dev.097097

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaouri, K., Maini, P. K., Skourides, P. A., Christodoulou, N., and Chapman, S. J. (2019). A simple mechanochemical model for calcium signalling in embryonic epithelial cells. J. Math. Biol. 78, 2059–2092. doi:10.1007/s00285-019-01333-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaouri, K., Méndez, P. E., and Ruiz-Baier, R. (2022). Mechanochemical models for calcium waves in embryonic epithelia. Vietnam J. Math. 50, 947–975. doi:10.1007/s10013-022-00579-y

CrossRef Full Text | Google Scholar

Karcher, H., Lammerding, J., Huang, H., Lee, R. T., Kamm, R. D., and Kaazempur-Mofrad, M. R. (2003). A three-dimensional viscoelastic model for cell deformation with experimental verification. Biophys. J. 85, 3336–3349. doi:10.1016/S0006-3495(03)74753-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Kraning-Rush, C. M., Califano, J. P., and Reinhart-King, C. A. (2012). Cellular traction stresses increase with increasing metastatic potential. PloS one 7, e32572. doi:10.1371/journal.pone.0032572

PubMed Abstract | CrossRef Full Text | Google Scholar

Lange, J. R., and Fabry, B. (2013). Cell and tissue mechanics in cell migration. Exp. Cell Res. 319, 2418–2423. doi:10.1016/j.yexcr.2013.04.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Langtangen, H. P., and Logg, A. (2017). Solving PDEs in Python: The FEniCS tutorial I. Oslo: Springer.

Google Scholar

Lecuit, T., and Lenne, P. F. (2007). Cell surface mechanics and the control of cell shape, tissue patterns and morphogenesis. Nat. Rev. Mol. Cell Biol. 8, 633–644. doi:10.1038/nrm2222

PubMed Abstract | CrossRef Full Text | Google Scholar

Moreo, P., Gaffney, E. A., García-Aznar, J. M., and Doblaré, M. (2010). On the modelling of biological patterns with mechanochemical models: Insights from analysis and computation. Bull. Math. Biol. 72, 400–431. doi:10.1007/s11538-009-9452-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Murray, J. D. (2001). Mathematical biology II: Spatial models and biomedical applications. New York: Springer.

Google Scholar

Murray, J. D. (2003). On the mechanochemical theory of biological pattern formation with application to vasculogenesis. C. R. Biol. 326, 239–252. doi:10.1016/s1631-0691(03)00065-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Murray, J. D., and Oster, G. F. (1984). Generation of biological pattern and form. IMA J. Math. Appl. Med. Biol. 1, 51–75. doi:10.1093/imammb/1.1.51

PubMed Abstract | CrossRef Full Text | Google Scholar

Narciso, C. E., Contento, N. M., Storey, T. J., Hoelzle, D. J., and Zartman, J. J. (2017). Release of applied mechanical loading stimulates intercellular calcium waves in Drosophila wing discs. Biophys. J. 113, 491–501. doi:10.1016/j.bpj.2017.05.051

PubMed Abstract | CrossRef Full Text | Google Scholar

Neville, A., Matthews, P., and Byrne, H. (2006). Interactions between pattern formation and domain growth. Bull. Math. Biol. 68, 1975–2003. doi:10.1007/s11538-006-9060-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Nguyen, T. L., Polanco, E. R., Patananan, A. N., Zangle, T. A., and Teitell, M. A. (2020). Cell viscoelasticity is linked to fluctuations in cell biomass distributions. Sci. Rep. 10, 7403–7411. doi:10.1038/s41598-020-64259-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Oakes, P. W., Banerjee, S., Marchetti, M. C., and Gardel, M. L. (2014). Geometry regulates traction stresses in adherent cells. Biophys. J. 107, 825–833. doi:10.1016/j.bpj.2014.06.045

PubMed Abstract | CrossRef Full Text | Google Scholar

Parvini, C. H., Cartagena-Rivera, A. X., and Solares, S. D. (2022). Viscoelastic parameterization of human skin cells characterize material behavior at multiple timescales. Commun. Biol. 5, 17. doi:10.1038/s42003-021-02959-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Ruiz-Baier, R., Gizzi, A., Rossi, S., Cherubini, C., Laadhari, A., Filippi, S., et al. (2014). Mathematical modelling of active contraction in isolated cardiomyocytes. Math. Med. Biol. 31, 259–283. doi:10.1093/imammb/dqt009

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanderson, M. J., Charles, A., and Dirksen, E. R. (1990). Mechanical stimulation and intercellular communication increases intracellular Ca2+ in epithelial cells. Cell Regul. 1, 585–596. doi:10.1091/mbc.1.8.585

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanderson, M. J., Chow, I., and Dirksen, E. R. (1988). Intercellular communication between ciliated cells in culture. Am. J. Physiol. 254, C63–C74. doi:10.1152/ajpcell.1988.254.1.C63

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanderson, M. J., and Sleigh, M. A. (1981). Ciliary activity of cultured rabbit tracheal epithelium: Beat pattern and metachrony. J. Cell Sci. 47, 331–347. doi:10.1242/jcs.47.1.331

PubMed Abstract | CrossRef Full Text | Google Scholar

Smedley, M. J., and Stanisstreet, M. (1986). Calcium and neurulation in mammalian embryos: II. Effects of cytoskeletal inhibitors and calcium antagonists on the neural folds of rat embryos. Development 93, 167–178. doi:10.1242/dev.93.1.167

PubMed Abstract | CrossRef Full Text | Google Scholar

Suzuki, M., Sato, M., Koyama, H., Hara, Y., Hayashi, K., Yasue, N., et al. (2017). Distinct intracellular ca2+ dynamics regulate apical constriction and differentially contribute to neural tube closure. Development 144, 1307–1316. doi:10.1242/dev.141952

PubMed Abstract | CrossRef Full Text | Google Scholar

Tsutsumi, M., Inoue, K., Denda, S., Ikeyama, K., Goto, M., and Denda, M. (2009). Mechanical-stimulation-evoked calcium waves in proliferating and differentiated human keratinocytes. Cell Tissue Res. 338, 99–106. doi:10.1007/s00441-009-0848-0

PubMed Abstract | CrossRef Full Text | Google Scholar

von Dassow, M., and Davidson, L. A. (2011). Physics and the canalization of morphogenesis: A grand challenge in organismal biology. Phys. Biol. 8, 045002. doi:10.1088/1478-3975/8/4/045002

PubMed Abstract | CrossRef Full Text | Google Scholar

Wallingford, J. B., Ewald, A. J., Harland, R. M., and Fraser, S. E. (2001). Calcium signaling during convergent extension in Xenopus. Curr. Biol. 11, 652–661. doi:10.1016/s0960-9822(01)00201-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Wiebe, C., and Brodland, G. W. (2005). Tensile properties of embryonic epithelia measured using a novel instrument. J. Biomech. 38, 2087–2094. doi:10.1016/j.jbiomech.2004.09.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamaguchi, N., Zhang, Z., Schneider, T., Wang, B., Panozzo, D., and Knaut, H. (2022). Rear traction forces drive adherent tissue migration in vivo. Nat. Cell Biol. 24, 194–204. doi:10.1038/s41556-022-00844-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, W., Chen, J. Y., and Zhou, L. (2009). Effects of shear stress on intracellular calcium change and histamine release in rat basophilic leukemia (RBL-2H3) cells. J. Environ. Pathol. Toxicol. Oncol. 28, 223–230. doi:10.1615/jenvironpatholtoxicoloncol.v28.i3.30

PubMed Abstract | CrossRef Full Text | Google Scholar

Young, S., Ennes, H., McRoberts, J., Chaban, V., Dea, S., and Mayer, E. (1999). Calcium waves in colonic myocytes produced by mechanical and receptor-mediated stimulation. Am. J. Physiol. 276, G1204–G1212. doi:10.1152/ajpgi.1999.276.5.G1204

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, J., Kim, H. Y., and Davidson, L. A. (2009). Actomyosin stiffens the vertebrate embryo during crucial stages of elongation and neural tube closure. Development 136, 677–688. doi:10.1242/dev.026211

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, J., Pal, S., Maiti, S., and Davidson, L. A. (2015). Force production and mechanical accommodation during convergent extension. Development 142, 692–701. doi:10.1242/dev.116533

PubMed Abstract | CrossRef Full Text | Google Scholar