Abstract
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 IP3 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 Ca2+ plays a crucial role in regulating AC during neural tube closure (NTC) (Christodoulou and Skourides, 2015; Suzuki et al., 2017). Furthermore, pharmacologically inhibiting Ca2+ 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 Ca2+ 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 Ca2+ is also well established; mechanically stimulated Ca2+ 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 Ca2+ 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 Ca2+ 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 Ca2+ 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 Ca2+ release and Ca2+ elicits contractions which are sensed as mechanical stimuli by the cell, a two-way mechanochemical feedback between Ca2+ 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 Ca2+ 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 Ca2+ signalling and mechanical forces during AC.
This paper extends the mechanochemical model in (Kaouri et al., 2019), which describes the coupling of Ca2+ 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 Ca2+ concentration. The models in (Kaouri et al., 2019; Kaouri et al., 2022), as well as the model presented here, employ the well-established Ca2+ signalling model from Atri et al. (1993), called the “Atri model” hereafter. The Atri model captures the experimentally verified Ca2+-induced Ca2+ release (CICR) process. It consists of two differential equations, one PDE for the cytosolic Ca2+ concentration and another PDE for the percentage of the non-inactivated IP3 receptors on the endoplasmic reticulum (ER) which allow release of Ca2+ 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 Ca2+ concentration, following the formulation in (Murray and Oster, 1984; Banerjee and Marchetti, 2011; Kaouri et al., 2019; Kaouri et al., 2022).
FIGURE 1
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 Ca2+ 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:
where [Ca2+] = 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, α1 and α2 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:
The function Irand multiplying the CICR Ca2+ flux is a random-in-space distribution of Ca2+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 bywhere χM denotes the characteristic function on the generic subdomain M, and by Ein and Eout 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 bywhere c0 is the steady state value of c. We also assume stress-free and zero-flux boundary conditions on the domain boundary, as follows: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 isNote 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 Ca2+ signalling model we use (Atri et al., 1993) captures the Ca2+ release to the cytosol via the IPR/Ca2+ 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 Ein. To determine the Young’s modulus for the epidermal cells surrounding the neuroepithelium, Eout, we use the ratio between the stiffness moduli of the neuroepithelium and epidermis, as determined in (Wiebe and Brodland, 2005); we, thus, estimate Eout = 0.55Ein.
TABLE 1
| Parameter list | |||
|---|---|---|---|
| Parameter | Definition | Value | Source |
| Ein | Young’s modulus on the neural plate | 44.26 Pa | Zhou et al. (2009) |
| Eout = 0.55Ein | Young’s modulus on the epidermal layer | 24.34 Pa | Wiebe and Brodland, (2005); Zhou et al., (2009) |
| T0 | Traction stress | 50 − 450 Pa | Kraning-Rush et al., (2012); Oakes et al., (2014); Yamaguchi et al., (2022) |
| ν | Poisson’s ratio | 0.4 | D’Angelo et al., (2019); Zhou et al., (2015) |
| Shear viscosity | 3790 Pa s | D’Angelo et al., (2019); Zhou et al., (2009) | |
| Bulk viscosity | 550 Pa s | D’Angelo et al., (2019); Zhou et al., (2009) | |
Parameter values for the mechanochemical model.
We determine the shear and bulk viscosities, and 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, T0, 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, T0 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 . We have 256 cells (Suzuki et al., 2017) tightly packed in the tissue and, hence, the initial tissue area is approximately 64,000μm2. We assume that the tissue is a disc-shaped domain of radius . 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, K1 = 46.28, G = 5.71, and K = 0.14. The three parameters we are going to vary in the simulations areAs the Young modulus, E, is discontinuous across the interface of the neuroepithelium and the epidermis, so are the parameters α1, α2 and β1. From the nondimensionalisation it arises that and , where τj = 2s and other values as in Table 1. On the neuroepithelium we, hence, have and , whereas on the epidermis we take and .
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 β1 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, Irand, we impose a frequency that is linearly increasing from 0.1 to 0.4 and set the amplitude to 1 + ampl, with ampl 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
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 , that is when T0 = 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 Irand) exist elsewhere in the domain and increase in amplitude and frequency as time progresses (see Figure 2).
FIGURE 3
FIGURE 4
Increasing T0 to 250 Pa gives . In this case we see in Figure 4 that the pressure and displacement approximately double in magnitude compared to those for T0 = 100 Pa .
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 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 (corresponding to T0 = 50, 100, 150, 250 Pa and changing accordingly ); we see that the area reduction is quite sensitive to the choice of the active contractile force parameter, T0, 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 . The constriction rate is approximately constant, as identified in the experiments of (Christodoulou and Skourides, 2015).
FIGURE 5
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, , 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 7
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, T0 = 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 , since β1 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 (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
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 IP3, λ, which measures the strength of the stretch-sensitive activation and which represents the maximum contraction stress.
The model shows that for any value of μ, λ, and 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 , 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.
Statements
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
1
AmestoyP. R.DuffI. S.L’ExcellentJ. Y.KosterJ. (2000). Mumps: A general purpose distributed memory sparse solver. Int. Workshop Appl. Parallel Comput., 121–130. 10.1007/3-540-70734-4_16
2
AtriA.AmundsonJ.ClaphamD.SneydJ. (1993). A single-pool model for intracellular calcium oscillations and waves in the xenopus laevis oocyte. Biophys. J.65, 1727–1739. 10.1016/S0006-3495(93)81191-3
3
BanerjeeS.MarchettiM. C. (2011). Instabilities and oscillations in isotropic active gels. Soft Matter7, 463–473. 10.1039/c0sm00494d
4
BenkoR.BrodlandG. W. (2007). Measurement of in vivo stress resultants in neurulation-stage amphibian embryos. Ann. Biomed. Eng.35, 672–681. 10.1007/s10439-006-9250-1
5
Bereiter-HahnJ. (2005). Mechanics of crawling cells. Med. Eng. Phys.27, 743–753. 10.1016/j.medengphy.2005.04.021
6
BrücknerB. R.JanshoffA. (2015). Elastic properties of epithelial cells probed by atomic force microscopy. Biochim. Biophys. Acta1853, 3075–3082. 10.1016/j.bbamcr.2015.07.010
7
CharlesA. C.DirksenE. R.MerrillJ. E.SandersonM. J. (1993). Mechanisms of intercellular calcium signaling in glial cells studied with dantrolene and thapsigargin. Glia7, 134–145. 10.1002/glia.440070203
8
CharlesA. C.MerrillJ. E.DirksenE. R.SandersontM. J. (1991). Intercellular signaling in glial cells: Calcium waves and oscillations in response to mechanical stimulation and glutamate. Neuron6, 983–992. 10.1016/0896-6273(91)90238-u
9
CharlesA. C.NausC.ZhuD.KidderG. M.DirksenE. R.SandersonM. J. (1992). Intercellular calcium signaling via gap junctions in glioma cells. J. Cell Biol.118, 195–201. 10.1083/jcb.118.1.195
10
ChristodoulouN.SkouridesP. 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. 10.1016/j.celrep.2015.11.017
11
CioncoliniA.BoffiD. (2019). The MINI mixed finite element for the Stokes problem: An experimental investigation. Comput. Math. Appl.77, 2432–2446. 10.1016/j.camwa.2018.12.028
12
D’AngeloA.DierkesK.CarolisC.SalbreuxG.SolonJ. (2019). In vivo force application reveals a fast tissue softening and external friction increase during early embryogenesis. Curr. Biol.29, 1564–1571. 10.1016/j.cub.2019.04.010
13
DavidsonL. A. (2011). Embryo mechanics. Curr. Top. Dev. Biol., 215–241.
14
Guiu-SoutoJ.MunuzuriA. P. (2015). Influence of oscillatory centrifugal forces on the mechanism of Turing pattern formation. Phys. Rev. E Stat. Nonlin. Soft Matter Phys.91, 012917. 10.1103/PhysRevE.91.012917
15
HerrgenL.VossO. P.AkermanC. J. (2014). Calcium-dependent neuroepithelial contractions expel damaged cells from the developing brain. Dev. Cell31, 599–613. 10.1016/j.devcel.2014.10.012
16
HerrmannL. R. (1965). Elasticity equations for incompressible and nearly incompressible materials by a variational theorem. AIAA J.3, 1896–1900. 10.2514/3.3277
17
HunterG. L.CrawfordJ. M.GenkinsJ. Z.KiehartD. P. (2014). Ion channels contribute to the regulation of cell sheet forces during Drosophila dorsal closure. Development141, 325–334. 10.1242/dev.097097
18
KaouriK.MainiP. K.SkouridesP. A.ChristodoulouN.ChapmanS. J. (2019). A simple mechanochemical model for calcium signalling in embryonic epithelial cells. J. Math. Biol.78, 2059–2092. 10.1007/s00285-019-01333-8
19
KaouriK.MéndezP. E.Ruiz-BaierR. (2022). Mechanochemical models for calcium waves in embryonic epithelia. Vietnam J. Math.50, 947–975. 10.1007/s10013-022-00579-y
20
KarcherH.LammerdingJ.HuangH.LeeR. T.KammR. D.Kaazempur-MofradM. R. (2003). A three-dimensional viscoelastic model for cell deformation with experimental verification. Biophys. J.85, 3336–3349. 10.1016/S0006-3495(03)74753-5
21
Kraning-RushC. M.CalifanoJ. P.Reinhart-KingC. A. (2012). Cellular traction stresses increase with increasing metastatic potential. PloS one7, e32572. 10.1371/journal.pone.0032572
22
LangeJ. R.FabryB. (2013). Cell and tissue mechanics in cell migration. Exp. Cell Res.319, 2418–2423. 10.1016/j.yexcr.2013.04.023
23
LangtangenH. P.LoggA. (2017). Solving PDEs in Python: The FEniCS tutorial I. Oslo: Springer.
24
LecuitT.LenneP. F. (2007). Cell surface mechanics and the control of cell shape, tissue patterns and morphogenesis. Nat. Rev. Mol. Cell Biol.8, 633–644. 10.1038/nrm2222
25
MoreoP.GaffneyE. A.García-AznarJ. M.DoblaréM. (2010). On the modelling of biological patterns with mechanochemical models: Insights from analysis and computation. Bull. Math. Biol.72, 400–431. 10.1007/s11538-009-9452-4
26
MurrayJ. D. (2001). Mathematical biology II: Spatial models and biomedical applications. New York: Springer.
27
MurrayJ. D. (2003). On the mechanochemical theory of biological pattern formation with application to vasculogenesis. C. R. Biol.326, 239–252. 10.1016/s1631-0691(03)00065-9
28
MurrayJ. D.OsterG. F. (1984). Generation of biological pattern and form. IMA J. Math. Appl. Med. Biol.1, 51–75. 10.1093/imammb/1.1.51
29
NarcisoC. E.ContentoN. M.StoreyT. J.HoelzleD. J.ZartmanJ. J. (2017). Release of applied mechanical loading stimulates intercellular calcium waves in Drosophila wing discs. Biophys. J.113, 491–501. 10.1016/j.bpj.2017.05.051
30
NevilleA.MatthewsP.ByrneH. (2006). Interactions between pattern formation and domain growth. Bull. Math. Biol.68, 1975–2003. 10.1007/s11538-006-9060-5
31
NguyenT. L.PolancoE. R.PatanananA. N.ZangleT. A.TeitellM. A. (2020). Cell viscoelasticity is linked to fluctuations in cell biomass distributions. Sci. Rep.10, 7403–7411. 10.1038/s41598-020-64259-y
32
OakesP. W.BanerjeeS.MarchettiM. C.GardelM. L. (2014). Geometry regulates traction stresses in adherent cells. Biophys. J.107, 825–833. 10.1016/j.bpj.2014.06.045
33
ParviniC. H.Cartagena-RiveraA. X.SolaresS. D. (2022). Viscoelastic parameterization of human skin cells characterize material behavior at multiple timescales. Commun. Biol.5, 17. 10.1038/s42003-021-02959-5
34
Ruiz-BaierR.GizziA.RossiS.CherubiniC.LaadhariA.FilippiS.et al (2014). Mathematical modelling of active contraction in isolated cardiomyocytes. Math. Med. Biol.31, 259–283. 10.1093/imammb/dqt009
35
SandersonM. J.CharlesA.DirksenE. R. (1990). Mechanical stimulation and intercellular communication increases intracellular Ca2+ in epithelial cells. Cell Regul.1, 585–596. 10.1091/mbc.1.8.585
36
SandersonM. J.ChowI.DirksenE. R. (1988). Intercellular communication between ciliated cells in culture. Am. J. Physiol.254, C63–C74. 10.1152/ajpcell.1988.254.1.C63
37
SandersonM. J.SleighM. A. (1981). Ciliary activity of cultured rabbit tracheal epithelium: Beat pattern and metachrony. J. Cell Sci.47, 331–347. 10.1242/jcs.47.1.331
38
SmedleyM. J.StanisstreetM. (1986). Calcium and neurulation in mammalian embryos: II. Effects of cytoskeletal inhibitors and calcium antagonists on the neural folds of rat embryos. Development93, 167–178. 10.1242/dev.93.1.167
39
SuzukiM.SatoM.KoyamaH.HaraY.HayashiK.YasueN.et al (2017). Distinct intracellular ca2+ dynamics regulate apical constriction and differentially contribute to neural tube closure. Development144, 1307–1316. 10.1242/dev.141952
40
TsutsumiM.InoueK.DendaS.IkeyamaK.GotoM.DendaM. (2009). Mechanical-stimulation-evoked calcium waves in proliferating and differentiated human keratinocytes. Cell Tissue Res.338, 99–106. 10.1007/s00441-009-0848-0
41
von DassowM.DavidsonL. A. (2011). Physics and the canalization of morphogenesis: A grand challenge in organismal biology. Phys. Biol.8, 045002. 10.1088/1478-3975/8/4/045002
42
WallingfordJ. B.EwaldA. J.HarlandR. M.FraserS. E. (2001). Calcium signaling during convergent extension in Xenopus. Curr. Biol.11, 652–661. 10.1016/s0960-9822(01)00201-9
43
WiebeC.BrodlandG. W. (2005). Tensile properties of embryonic epithelia measured using a novel instrument. J. Biomech.38, 2087–2094. 10.1016/j.jbiomech.2004.09.005
44
YamaguchiN.ZhangZ.SchneiderT.WangB.PanozzoD.KnautH. (2022). Rear traction forces drive adherent tissue migration in vivo. Nat. Cell Biol.24, 194–204. 10.1038/s41556-022-00844-9
45
YangW.ChenJ. Y.ZhouL. (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. 10.1615/jenvironpatholtoxicoloncol.v28.i3.30
46
YoungS.EnnesH.McRobertsJ.ChabanV.DeaS.MayerE. (1999). Calcium waves in colonic myocytes produced by mechanical and receptor-mediated stimulation. Am. J. Physiol.276, G1204–G1212. 10.1152/ajpgi.1999.276.5.G1204
47
ZhouJ.KimH. Y.DavidsonL. A. (2009). Actomyosin stiffens the vertebrate embryo during crucial stages of elongation and neural tube closure. Development136, 677–688. 10.1242/dev.026211
48
ZhouJ.PalS.MaitiS.DavidsonL. A. (2015). Force production and mechanical accommodation during convergent extension. Development142, 692–701. 10.1242/dev.116533
Summary
Keywords
mechanochemical model, calcium signalling, linear viscoelasticity, embryogenesis, apical constriction, computational modelling
Citation
Kaouri K, Christodoulou N, Chakraborty A, Méndez PE, Skourides P and Ruiz-Baier R (2022) A new mechanochemical model for apical constriction: Coupling calcium signalling and viscoelasticity. Front. Syst. Biol. 2:962790. doi: 10.3389/fsysb.2022.962790
Received
06 June 2022
Accepted
05 October 2022
Published
28 October 2022
Volume
2 - 2022
Edited by
Jennifer Anne Flegg, University of Melbourne, Australia
Reviewed by
Linda Irons, Yale University, United States
Riccardo Sacco, Politecnico di Milano, Italy
Updates
Copyright
© 2022 Kaouri, Christodoulou, Chakraborty, Méndez, Skourides and Ruiz-Baier.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Katerina Kaouri, kaourik@cardiff.ac.uk
This article was submitted to Multiscale Mechanistic Modeling, a section of the journal Frontiers in Systems Biology
Disclaimer
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.