A simple mechanochemical model for calcium signalling in embryonic epithelial cells.

Calcium signalling is one of the most important mechanisms of information propagation in the body. In embryogenesis the interplay between calcium signalling and mechanical forces is critical to the healthy development of an embryo but poorly understood. Several types of embryonic cells exhibit calcium-induced contractions and many experiments indicate that calcium signals and contractions are coupled via a two-way mechanochemical feedback mechanism. We present a new analysis of experimental data that supports the existence of this coupling during apical constriction. We then propose a simple mechanochemical model, building on early models that couple calcium dynamics to the cell mechanics and we replace the hypothetical bistable calcium release with modern, experimentally validated calcium dynamics. We assume that the cell is a linear, viscoelastic material and we model the calcium-induced contraction stress with a Hill function, i.e. saturating at high calcium levels. We also express, for the first time, the "stretch-activation" calcium flux in the early mechanochemical models as a bottom-up contribution from stretch-sensitive calcium channels on the cell membrane. We reduce the model to three ordinary differential equations and analyse its bifurcation structure semi-analytically as two bifurcation parameters vary-the [Formula: see text] concentration, and the "strength" of stretch activation, [Formula: see text]. The calcium system ([Formula: see text], no mechanics) exhibits relaxation oscillations for a certain range of [Formula: see text] values. As [Formula: see text] is increased the range of [Formula: see text] values decreases and oscillations eventually vanish at a sufficiently high value of [Formula: see text]. This result agrees with experimental evidence in embryonic cells which also links the loss of calcium oscillations to embryo abnormalities. Furthermore, as [Formula: see text] is increased the oscillation amplitude decreases but the frequency increases. Finally, we also identify the parameter range for oscillations as the mechanical responsiveness factor of the cytosol increases. This work addresses a very important and not well studied question regarding the coupling between chemical and mechanical signalling in embryogenesis.

Introduction

Calcium signalling is one of the most important mechanisms of information propagation in the body, playing an important role as a second messenger in several processes such as embryogenesis, heart function, blood clotting, muscle contraction and diseases of the muscular and nervous systems (Berridge et al. 2000; Brini and Carafoli 2009; Dupont et al. 2016). Through the sensing mechanisms of cells, external environmental stimuli are transformed into intracellular or intercellular calcium signals that often take the form of oscillations and waves.

In this work we will focus on the interplay of calcium signalling and mechanical forces in embryogenesis. During embryogenesis, cells and tissues generate physical forces, change their shape, move and proliferate (Lecuit and Lenne 2007). The impact of these forces on morphogenesis is directly linked to calcium signalling (Hunter et al. 2014). In general, how the mechanics of the cell and tissue are regulated and coupled to the cellular biochemical response is essential for understanding embryogenesis. Understanding this mechanochemical coupling, in particular when calcium signalling is involved, is also important for elucidating a wide range of other body processes, such as wound healing (Antunes et al. 2013; Herrgen et al. 2014) and cancer (Basson et al. 2015).

Calcium plays a crucial role in every stage of embryonic development starting with fast calcium waves during fertilization (Deguchi et al. 2000) to calcium waves involved in convergent extension movements during gastrulation (Wallingford et al. 2001), to calcium transients regulating neural tube closure (Christodoulou and Skourides 2015), morphological patterning in the brain (Sahu et al. 2017; Webb and Miller 2007) and apical-basal cell thinning in the enveloping layer cells (Zhang et al. 2011), either in the form of calcium waves or through Wnt/ signalling (Christodoulou and Skourides 2015; Herrgen et al. 2014; Hunter et al. 2014; Kühl et al. 2000a, b; Narciso et al. 2017; Slusarski et al. 1997a, b; Suzuki et al. 2017; Wallingford et al. 2001). Crucially, pharmacologically inhibiting calcium has been shown to lead to embryo defects (Christodoulou and Skourides 2015; Smedley and Stanisstreet 1986; Wallingford et al. 2001).

In many experiments actomyosin-based contractions have been documented in response to calcium release in both embryonic and cultured cells (Christodoulou and Skourides 2015; Herrgen et al. 2014; Hunter et al. 2014; Suzuki et al. 2017; Wallingford et al. 2001) and it has become clear that calcium is responsible for contractions in both muscle and non-muscle cells, albeit through different mechanisms (Cooper 2000). Cell contraction in striated muscle is mediated by the binding of to troponin but in non-muscle cells (and in smooth muscle cells) contraction is mediated by phosphorylation of the regulatory light chain of myosin. This phosphorylation promotes the assembly of myosin into filaments, and it increases myosin activity. Myosin light-chain kinase (MLCK), which is responsible for this phosphorylation, is itself regulated by calmodulin, a well-characterized and ubiquitously expressed protein regulated by calcium (Scholey et al. 1980). Elevated cytosolic calcium promotes binding of calmodulin to MLCK, resulting in its activation, subsequent phosphorylation of the myosin regulatory light chain and then contraction. Thus, cytosolic calcium elevation is an ubiquitous signal for cell contraction which manifests in various ways (Cooper 2000).

In some tissues these contractions give rise to well defined changes in cell shape. One such example is apical constriction (AC), an intensively studied morphogenetic process central to embryonic development in both vertebrates and invertebrates (Vijayraghavan and Davidson 2017). In AC the apical surface of an epithelial cell constricts, leading to dramatic changes in cell shape. Such shape changes drive epithelial sheet bending and invagination and are indispensable for tissue and organ morphogenesis including gastrulation in C. elegans and Drosophila and vertebrate neural tube formation (Christodoulou and Skourides 2015; Rohrschneider and Nance 2009; Sawyer et al. 2010).

On the other hand, the ability of cells to sense and respond to forces by elevating their cytosolic calcium is well established. Mechanically stimulated calcium waves have been observed propagating through ciliated tracheal epithelial cells (Sanderson et al. 1990, 1988; Sanderson and Sleigh 1981), rat brain glial cells (Charles et al. 1993, 1991, 1992), keratinocytes (Tsutsumi et al. 2009), developing epithelial cells in Drosophila wing discs (Narciso et al. 2017) and many other cell types (Beraeiter-Hahn 2005; Tsutsumi et al. 2009; Yang et al. 2009; Young et al. 1999). Thus, different types of mechanical stimuli, from shear stress to direct mechanical stimulation, can elicit calcium elevation (the sensing mechanism may differ in each case). So, since mechanical stimulation elicits calcium release and calcium elicits contractions which are sensed as mechanical stimuli by the cell, it is clear that a two-way mechanochemical feedback between calcium and contractions should be at play.

This two-way feedback is supported by our work here with a new analysis of data from earlier experiments conducted by two of the authors (Christodoulou and Skourides 2015); we present this analysis in detail in Sect. 2. The analysis shows that in contracting cells, in the Xenopus neural plate, calcium oscillations become more frequent and also increase in amplitude as the calcium-elicited surface area reduction progresses. This suggests that as the increased tension around the contracting cell is sensed, it leads to more calcium release and in turn to more contractions, and so on. In addition, experiments in Drosophila also support the hypothesis that a mechanochemical feedback loop is in play (Saravanan et al. 2013; Solon et al. 2009). Thus, data from these two model systems clearly show that mechanical forces generated by contraction influence calcium release and the contraction cycle. The mechanosensing takes place via, as yet undefined, mechanosensory molecules which could be mechanogated ion channels, mechanosensitive proteins at adherens junctions like alpha catenin, or even integrins which have recently been shown to become activated by plasma membrane tension in the absence of ligands (Delmas and Coste 2013; Petridou and Skourides 2016; Yao et al. 2014).

Given the broad range of critical biological processes involving calcium signalling and its coupling to mechanical effects, modelling this mechanochemical coupling is of great interest. Therefore, we develop a simple mechanochemical model that captures the essential elements of a two-way coupling between calcium signalling and contractions in embryonic cells. The first mechanochemical models for embryogenesis were developed by Oster, Murray and collaborators in the 80s (Murray 2001; Murray et al. 1988; Murray and Oster 1984; Oster and Odell 1984). Calcium evolution in those early models was modelled with a hypothetical bistable reaction-diffusion process in which the application of stress can switch the calcium state from low to high stable concentration. We now know that the calcium dynamics are more complicated, so our mechanochemical model includes instead the calcium dynamics of the experimentally verified model in Atri et al. (1993), which captures the experimentally observed Calcium-Induced-Calcium Release (CICR) process and the dynamics of the receptors on the Endoplasmic Reticulum (ER). In this way we update the early mechanochemical models for embryonic cells in line with recent advances in calcium signalling. Note that there are many recent models of calcium signalling induced by mechanical stimulation, for example for mammalian airway epithelial cells (Warren et al. 2010), for keratinocytes (Kobayashi et al. 2014), for white blood cells (Yao et al. 2016), and for retinal pigment epithelial cells (Vainio et al. 2015). However, these models do not include a two-way coupling between calcium signalling and mechanics.

Calcium is stored and released from intracellular stores, such as the ER, or the Sarcoplasmic Reticulum (SR), according to the well-established nonlinear feedback mechanism of CICR (Dupont et al. 2016). There are many models for calcium oscillations, all capturing the CICR process. Many of them model the receptors on the ER in some manner, and they can be classified as Class I or Class II models (Dupont et al. 2016). In all Class I models is a control parameter and oscillations can be sustained at a constant concentration. Oscillations exist for a window of values; the oscillations are excited at a threshold value and they vanish at a suffuciently high value. The Atri et al. model (1993) is an established Class I model, validated with experimental findings (Estrada et al. 2016). (We will call this model the ‘Atri model’ from now on.) It also has a mathematical structure that allows us to investigate our mechanochemical model semi-analytically and easily identify the parameter range sustaining calcium oscillations. Such an analysis cannot be done for other qualitatively similar, minimal Class I models as, for example, the more frequently used Li-Rinzel model (Li and Rinzel 1994); this is one of the contributions of this work.

Another contribution of our work is that we interpret the “stretch-activation” calcium flux from the outside medium, introduced in an ad hoc manner in the early mechanochemical models, as a “bottom-up” contribution from recently identified, stretch sensitive (stretch-activated) calcium channels (SSCCs) (Árnadóttir and Chalfie 2010; Dupont et al. 2016; Hamill 2006; Moore et al. 2010), in this way linking the channel scale with the whole cell scale.

The paper is organised as follows. In Sect. 2 we present a new analysis of experimental data which shows that calcium and contractions in embryonic cells must be involved in a two-way mechanochemical feedback mechanism. In Sect. 3 we develop a new mechanochemical model which captures the key ingredients of the two-way coupling. In Sect. 4 we analyse the model. In Sect. 4.1 we briefly revisit the analysis of the Atri model and show the bifurcation diagrams for the amplitude and frequency of calcium oscillations. In Sect. 4.2 we perform the bifurcation analysis of the mechanochemical model, varying the concentration and the strength of stretch activation, and we identify the parameter range sustaining calcium oscillations. In Sect. 5.1 we model the calcium-induced contraction stress with a Hill function of order 1, and we plot the parameter range for which oscillations are sustained. In Sect. 5.1.2 we study the amplitude and frequency of the calcium oscillations. In Sect. 5.1.3 we investigate the bifurcation diagrams as the mechanical responsiveness of the cytosol to calcium varies. In Sect. 5.2 we consider a Hill function of order 2 and we again identify the parameter range for oscillations. In Sect. 6 we summarise our conclusions and discuss further research directions.

Calcium and contractions are involved in a feedback loop in apical constriction: a new analysis of experimental data

There is ample experimental evidence that mechanical stimulation of cells leads to calcium elevation (Beraeiter-Hahn 2005; Charles et al. 1991; Narciso et al. 2017; Sanderson et al. 1990, 1988; Sanderson and Sleigh 1981; Tsutsumi et al. 2009; Young et al. 1999) and that, in turn, contraction of the cytosol is elicited by calcium (Christodoulou and Skourides 2015; Herrgen et al. 2014; Hunter et al. 2014; Suzuki et al. 2017; Wallingford et al. 2001). Calcium signalling would therefore, at least in part, be regulated by a mechanochemical feedback loop whereby calcium-elicited contractions mechanically stimulate the cell, lead to more calcium release, then to more contractions and so on. In embryogenesis, and in particular during AC, where cells contract significantly, such a feedback loop should also be at play (Martin and Goldstein 2014); in this work we present a new analysis of experimental data in Christodoulou and Skourides (2015) which supports this. AC is a calcium-driven morphogenetic movement of epithelial tissues, central in the embryogenesis of both vertebrates and invertebrates (Vijayraghavan and Davidson 2017). The apical domain of epithelial cells constricts the apical surface area, and this leads to changes in the cell geometry that drive tissue bending; in Christodoulou and Skourides (2015) the formation of the neural tube in Xenopus frogs is studied and in Solon et al. (2009) dorsal closure in Drosophila is investigated.

In Solon et al. (2009) the constriction of mutants that exhibit disrupted myosin activation rescues apical myosin accumulation, suggesting that mechanically stimulating the cell can elicit contractions (Pouille et al. 2009). In addition, experiments using laser ablation, and other methodologies that reduce cell contractility, reveal that mechanical feedback non-autonomously regulates the amplitude and spatial propagation of pulsed contraction during AC (Saravanan et al. 2013) and that this process is driven by calcium (Hunter et al. 2014; Pouille et al. 2009; Saravanan et al. 2013). Therefore, reducing contractility reduces local tension and this suppresses contraction in the control cells. This suggests that mechanical feedback is important during AC.

Moreover, experimental evidence suggests that sensing of mechanical stimuli involves mechanogated ion channels; in Drosophila such ion channels are required for embryos to regulate force generation after laser ablation (Hunter et al. 2014); similarly during wound healing (Antunes et al. 2013).

Previously, two of the co-authors have shown that cell-autonomous, asynchronous calcium transients elicit contraction pulses, leading to the pulsed reduction of the apical surface area of individual neural epithelial cells during neural tube closure (NTC) in Xenopus (Christodoulou and Skourides 2015). Here, in order to investigate in detail the relationship between calcium, contraction and mechanical forces we present a new analysis of previously collected data (Christodoulou and Skourides 2015). For a single embryonic epithelial cell (in a tissue), we plot its apical surface area and calcium level over time in Fig. 1 and we see that both oscillate, with approximately the same frequency and that the calcium pulse precedes the contraction by 30–40 s. (Note that calcium oscillations emerge spontaneously without any periodic external stimulation.) More information about how Figure 1 is produced is found in Appendix A.1.

Fig. 1

Normalised apical surface area and amplitude of calcium oscillations in a single cell undergoing apical constriction. We see that calcium elevation always precedes the initiation of a contraction pulse. At calcium begins to rise and at  s the surface area starts decreasing. The surface area reduction is succeeded by relaxation and stabilization of the cell at a smaller surface area. (This happens repeatedly, leading to significant reduction of the surface area over time.) See Appendix A.1 for further details on how the figure is produced

In Fig. 2a we plot the frequency of calcium transients and the apical surface area over time, averaged over 10 cells. The frequency of calcium oscillations is clearly correlated with the reduction in the surface area - cells with a smaller surface area exhibit more frequent calcium oscillations. Also, in Fig. 2b, for the same 10 cells and in the same timeframe, we plot the calcium oscillation amplitude, which increases with time. Therefore, the reduction in the surface area correlates also with an increase in the amplitude of the calcium oscillations. Therefore, increased surface area reduction (i.e. increased tension and hence increased mechanical stimulation) correlates with increased frequency and increased amplitude, i.e. overall increased calcium release.

Fig. 2

a The normalised surface area reduction is correlated with increasing oscillation frequency (10 cells). b The amplitude of oscillations increases with time (10 cells). We used time lapse sequences from which the surface area of each cell was measured at and the average calcium oscillation frequency was calculated using a 10-minute window (i.e. calcium oscillations for each cell were monitored between and minutes). A 10-minute window was selected so that the typical cell undergoes at least one calcium pulse. (See Appendix A.1 for further details.)

Summarising, our analysis shows that calcium oscillations trigger contraction pulses that lead to pulsed reduction in the apical surface area over time. It also shows that the increasing localized tension in a contracting cell correlates with calcium pulses of higher frequency and larger amplitude, confirming the presence of a mechanochemical feedback loop.

A new mechanochemical model for embryonic epithelial cells

We develop a simple mechanochemical model that captures the essential components of a two-way coupling of contractions and calcium signals in embryonic epithelial cells undergoing AC. Since the cell machinery involved in the mechanochemical coupling is similar in most cell types (Cooper 2000) our model, with some modifications, can also be applicable to other cell types. The essential features of our model are a component modelling the cell mechanics and a component modelling calcium dynamics, coupled through a two-way feedback. Such models have been proposed by Oster, Murray and collaborators in the 80s (Murray 2001; Murray and Oster 1984; Murray et al. 1988; Oster and Odell 1984) and here we update those models by replacing the hypothetical bistable calcium dynamics with the experimentally verified calcium dynamics in Atri et al. (1993). We also replace the ad hoc stretch activation calcium flux in Murray (2001) with a “bottom-up” calcium release through the SSCCs, thus linking the channels’ characteristics to the whole cell scale. The model takes the formHere, c is the cytosolic calcium concentration, h is the fraction of receptors on the ER that have not been inactivated by calcium, and is the dilation/compression of the apical surface area of the cell. In ODE (1), models the flux of calcium from the ER into the cytosol through the receptors, is the fraction of receptors that have bound and is an increasing function of p, the concentration. The constant denotes the calcium flux when all receptors are open and activated, and b represents a basal current through the channel. represents the calcium flux pumped out of the cytosol where is the maximum rate of pumping of calcium from the cytosol and is the calcium concentration at which the rate of pumping from the cytosol is at half-maximum. models the calcium flux leaking into the cytosol from outside the cell. Note that from now on we will neglect as this is a good approximation for the embryonic epithelial cells we consider.

is the calcium flux due to the activated SSCCs. SSCCs have been identified experimentally in recent years - they are on the cell membrane and allow calcium to flow into the cytosol from the extracellular space. They are activated when exposed to mechanical stimulation and they close either by relaxation of the mechanical force or by adaptation to the mechanical force (Árnadóttir and Chalfie 2010; Dupont et al. 2016; Hamill 2006; Moore et al. 2010). The constant S represents the ‘strength’ of stretch activation. In Sect. 3.1 we will derive a relationship for S as a function of the characteristics of an SSCC.

The inactivation of the receptors by calcium acts on a slower timescale compared to activation (Dupont et al. 2016) and so the time constant for the dynamics of h, in ODE (2). In ODE (3) is a contraction stress that expresses how the stress in the cell depends on the cytosolic calcium level. The constants are, respectively, the shear and bulk viscosities of the cytosol and the constants and , where E and are, respectively, the Young’s modulus and the Poisson ratio.

Our mechanochemical model is an extension of the Atri model,since ODE (1) is ODE (4) but with added to the right hand side as an extra source term. The detailed derivation of the Atri model is presented in Atri et al. (1993), where it was initially formulated, and the reader is referred there for more details. The parameter values, which we take from Atri et al. (1993), are summarised in the Appendix, Table 1. The Atri model is one of the minimal Class I models, in which relaxation oscillations can be sustained at constant concentration (Dupont et al. 2016; Keener and Sneyd 1998). It was developed as a model for intracellular calcium oscillations in Xenopus oocytes but it has been subsequently used to model calcium dynamics in other cell types including glial cells (Wilkins and Sneyd 1998), mammalian spermatozoa (Olson et al. 2010), and keratinocytes (Kobayashi et al. 2014, 2016). In addition, modified Atri models have been developed in Shi et al. (2008), Harvey et al. (2011) and Liu and Li (2016). In Estrada et al. (2016) the Atri model was compared to seven other calcium dynamics models and it exhibited the best agreement with experiments along with the more frequently used Li-Rinzel model (Li and Rinzel 1994). The Atri model has a mathematical structure that allows us to perform a large part of our study analytically. The Atri model is also mathematically interesting because its relaxation oscillations have a different asymptotic structure to that of the well-known Van der Pol oscillator and similar excitable systems. We will present an asymptotic analysis of the Atri model and of our mechanochemical model in future work.

Now, we describe our modelling assumptions and the remaining components of the model in more detail.

Stretch-activation calcium flux

In the early mechanochemical models (Murray 2001) the stretch-activation flux, , was introduced in a somewhat ad hoc manner. Here, we derive it in a bottom-up manner, from the contribution of the SSCCs to the cytosolic calcium concentration.

A model for the opening and closing of SSCCs was developed in Vainio et al. (2015) for retinal pigment epithelial cells; we adapt it here for embryonic epithelial cells for which no such modelling has been performed. We denote by the proportion of channels in the closed state, and by the proportion of SSCCs in the open state. The calcium flux due to the SSCCs is proportional to the number of open channels so , where is the maximum calcium flux rate going through the SSCCs. As in Vainio et al. (2015), we propose that the evolution of is governed by the ODEwhere is the forward rate constant and is the backward rate constant. We assume here that is quasi-steady, i.e. remains approximately constant as calcium rapidly evolves. This is a reasonable approximation, as discussed in Section 2.6 of Dupont et al. (2016). Therefore,We linearise (7) since is small for a linear viscoelastic medium and under the additional assumption that is of order 1 at most. We obtainTherefore, we have derived, for the first time, an expression for S as a combination of and , linking in this way the characteristics of an SSCC to the macroscopic stretch-activation calcium flux.

Derivation of ODE (3)

We can obtain ODE (3) from the full force balance mechanical equation for a linear viscoelastic material. Linear viscoelasticity, at first glance, might not be an appropriate approximation for embryogenic tissue undergoing drastic changes, but recent experiments (Von Dassow et al. 2010) show it is reasonable. For a Kelvin-Voigt, linear viscoelastic material sustaining calcium-induced contractions the force balance equation can be written as follows (Landau and Lifshitz 1970; Murray 2001):where is the stress tensor, is the strain tensor, the displacement vector, is the dilation/compression of the material, and is the unit tensor. (The subscript t here denotes a partial derivative with respect to time.) and where E and v are the Young’s modulus and the Poisson’s ratio, respectively. is the contraction stress which depends on the cytosolic calcium (Scholey et al. 1980). In one spatial dimension and therefore (9) becomes, upon integrating with respect to x,The constant of integration since when , , and . Hence, we obtain ODE (3).

Nondimensionalised model

We nondimensionalise the mechanochemical model using and . Dropping bars for notational convenience we obtainIn (11) , , , and . In (13), and , and in (12) . Using the parameter values of Atri et al. (1993) (see Appendix, Table 1), we obtain , , and . Also, taking values of E, and of the viscosity from Zhou et al. (2009) ( Pa, and Pa.s) we find that is 0.4. For simplicity, and since the parameter values for the calcium dynamics are approximate, we fix . Furthermore, where is nondimensional, and we also fix . To our knowledge, there are no measured properties for SSCCs and therefore we take the ‘strength’ of stretch activation as a bifurcation parameter, to explore the behaviour of the model for a range of values.

Table 1

Parameter values, taken from Atri et al. (1993)

Parameter values
Parameter	Value
b	0.111
k_1	0.7μM
k_f	16.2μM/s
k_μ	0.7μM
γ	2μM/s
k_γ	0.1μM
β	0–0.02 μM/s
k_2	0.7μM
τ_h	2 s

Analysis of the model

The bifurcation diagrams of the Atri model (no mechanics)

The nondimensional Atri system isIn Appendix A.2 we carry out a linear stability analysis of (14)–(15) and a bifurcation analysis with as the bifurcation parameter and we find that the parameter range for relaxation oscillations (limit cycles) is , as in Atri et al. (1993). In Appendix A.2 more details on the bifurcation structure of the system are given.

In Fig. 3 we plot the bifurcation diagrams for the Atri system. In Fig. 3a we present the amplitude of oscillations. The left Hopf point (LHP) and the right Hopf point (RHP) are, respectively, at and . There are stable limit cycles and unstable limit cycles. The amplitude of oscillations increases with except for a small range of values near the RHP. In Fig. 3b the frequencies of the stable and of the unstable limit cycles are shown, respectively. The range of for which both a stable and an unstable limit cycle are sustained is clearly visible as the double-valued part of the curve. The limit point of oscillations at , where the stable and unstable limit cycle branches coalesce, is also visible. The frequency of the stable limit cycles increases slowly as increases and the lower, stable branch approximates the square root of , as predicted by bifurcation theory (Kuznetsov 2013).

Fig. 3

Bifurcation diagrams for the ODEs (14)–(15), as ( level) increases: a amplitude of calcium oscillations (limit cycles). The dots represent stable limit cycles and the dash-dotted part corresponds to unstable limit cycles (respectively blue and green colour online). The left Hopf point (LHP) and the right Hopf point (RHP) are indicated. b Frequency of the limit cycles

Linear stability analysis of the mechanochemical model

The steady states of the system (11)–(12) satisfyFor any , using (16), we can easily plot the steady states as a function of and . The Jacobian of (11)–(12) is given byand the characteristic polynomial is conveniently factorised aswhere represents the eigenvalues. As one eigenvalue is always equal to -1, the bifurcations of the system can be studied through the quadraticTo identify the - parameter range sustaining oscillations we seek the Hopf bifurcations, which satisfy Tr, Det, Discr, whereand substituting in (16) we obtainHence, we can easily obtain the Hopf curve, for any by parametrically plotting (21) and (22), with c as a parameter. The interior of the Hopf curve corresponds to an unstable spiral and approximates the - parameter space sustaining oscillations (limit cycles) in the full nonlinear system.

We also determine parametric expressions for the fold curve. Inside the fold curve there are three steady states: on the fold curve two of states coalesce, and outside the fold curve there is one steady state. Setting DetEquations (23) and (16) constitute a linear system for and , so we again easily derive parametric expressions for and .

Similarly, to determine parametric expressions for the curve on which Discr changes sign we set Discrwhich is quadratic in and linear in . Combining (24) with (16) we can again determine parametric expressions for and . In summary, we have developed a quick method for determining the three key curves characterising the geometry of the bifurcation diagram, for any.

It is of course, a fortunate accident of construction that we can obtain these analytical expressions for this particular model. Since our model is qualitatively similar to any other mechanochemical model that is based on Class I calcium models, the analytical progress we make here is very useful since the insights gained from it can be applied to other mechanochemical models. A different model would have a more complex set of linear stability equations, that look quite different, but that are fundamentally saying the exact same thing. Crucial to the behaviour is the shape of the manifolds rather than the detail of the algebraic expressions.

Illustrative examples

Contraction stress is a Hill function of order 1

Hopf curves

We assume that the calcium-induced stress is the Hill functionassuming that when the calcium level increases sufficiently the stress saturates to a constant value, . This is a reasonable assumption since the cells reach a point at which they stop responding mechanically to calcium since the molecules involved in contraction, calmodulin and myosin light chain kinases, saturate for sufficiently high calcium levels (Stefan et al. 2008). Also, when , i.e. we assume no contraction stress without calcium. is the rate of increase of at and is the scale of ‘ascent’ to the saturation level . Therefore, we can call the ‘mechanical responsiveness factor’ of the cytoskeleton to calcium.

Choosing as an illustrative example, in Fig. 4 we use (16) to plot the steady state as a function of , for selected increasing values of (equilibrium curve). For the equilibrium curve is qualitatively similar to that of the Atri model (see Fig. 3a) but at the curve changes qualitatively and a second non-zero steady state exists for , and a part of the curve corresponds to negative values of (see Appendix A.3 for details). For no steady state exists for positive and hence is the biologically relevant range of the model, for (see Appendix A.3).

Fig. 4

Steady states of the system (11)–(12) when as is increased, for selected –from right to left, the thick (red) curve is for and the dashed curve is for . (Plot done with Mathematica.)

In Fig. 5 we plot the Hopf curve and the fold curve. We observe the following: (i) for we recover the Hopf points and the fold points of the Atri model, as expected. (ii) As increases the range of that sustains oscillations decreases. There is a global minimum value of that can sustain oscillations, . (iii) The oscillations are suppressed for a critical maximum value of , , and the system is in a high calcium state. Overall, we conclude the following from the Hopf curve:Therefore, for fixed cytoskeletal mechanical responsiveness factor, , and for fixed parameter values as in Atri et al. (1993) a range of behaviours emerge as and vary: at low levels that do not elicit oscillations in the Atri system mechanical effects can elicit oscillations, for intermediate levels that do sustain oscillations in the Atri system increasing mechanical effects always leads to the oscillations vanishing, and for high levels that cannot sustain oscillations in the Atri system no amount of stretch activation can ever elicit oscillations.

Fig. 5

Geometry of the bifurcation diagrams of the system (11)–(12) when , plotted using the analytical, parametric expressions we derived for and . The Hopf curve is dashed (blue colour online) and the fold curve is drawn with a solid line (red colour online). The horizontal and vertical dashed lines correspond, respectively, to the maximum value of , , and to the minimum value of , , for which oscillations can be sustained. (Plot done with Mathematica.)

for low values the Atri system does not sustain oscillations but there are two possibilities for the mechanochemical model as increases:

for no increase in will ever elicit oscillations.

for when reaches a certain value, , oscillations are elicited, and decreases as approaches 0.289. The oscillations vanish at a larger value of .

for values for which the Atri system sustains oscillations () in the mechanochemical model oscillations eventually vanish at a critical . This critical decreases monotonically as increases towards 0.495.

for high values () no oscillations are sustained in the Atri system and a further increase in will never elicit oscillations.

Overall, we conclude that in this case mechanics can significantly affect calcium signalling. A very important prediction of the model is that oscillations vanish for sufficiently large stretch activation. This prediction agrees with the experiments reported in Christodoulou and Skourides (2015) (Figure 5D); when the cells were forced to enter a high, non-oscillatory calcium state they monotonically reduced their apical surface area without oscillations. Interestingly, although the loss of oscillations did not affect the reduction of the apical surface on average, it led to the disruption of the patterning involved in AC and neural tube closure failed, leading to severe embryo abnormality.

In fact, the model also agrees, qualitatively, with other experimental observations. Intracellular calcium levels (which are regulated by ) directly affect cell contractility (Christodoulou and Skourides 2015). At low levels of and hence low levels of calcium, cells are not able to contract and therefore AC does not take place. At a threshold value the system changes behaviour and calcium oscillations/transients appear (mathematically this corresponds to a bifurcation). The calcium oscillations enable the ratchet-like pulsating process of the AC to progress normally. At high levels of the cell has been shown to enter a high-calcium state with no oscillations, as mentioned above. (This corresponds to another bifurcation since the system changes its qualitative behaviour.)

Regarding bistability, note that the fold curve consists of two branches very close to each other since the Atri system is bistable for a very small range of concentrations. As increases this range decreases and eventually vanishes at , where the two fold curve branches merge.

Summarising, the parametric method we have developed allows us to easily plot the Hopf curve, and the two other important curves of the bifurcation diagram, for any functional form of we may choose, and thus examine quickly the effect of mechanics on calcium oscillations. We note that in the experiments of Christodoulou and Skourides (2015) the calcium-induced stress saturates to a non-zero level as calcium levels increase and hence we chose a that saturates. In other cell types it is possible that the cell can relax back to baseline stress and in such a case would not be described by a Hill function, and more experiments should be undertaken to determine the appropriate form of .

Amplitude and frequency of the calcium oscillations

We now determine numerically the amplitude and frequency of oscillations (limit cycles) of the system (11)–(12) when .

In Fig. 6 we plot the oscillation amplitude as a function of , for two selected values of , using XPPAUT. For (Fig. 6a) the Atri system has no oscillations but stable limit cycles arise in the mechanochemical model as is increased, which agrees with the Hopf curve in Fig. 5. For (Fig. 6b) the Atri system has a stable limit cycle and as increases, stable and unstable limit cycles emerge for a finite -interval, and oscillations eventually vanish for sufficiently large . For the Atri system has a stable limit cycle and as increases, stable and unstable limit cycles emerge for a finite range of values, and oscillations eventually vanish for a large enough value of .

Fig. 6

Amplitude of calcium oscillations for the system (11)–(12) when , as is increased, for: ab. The stable limit cycles are represented by dots and the unstable limit cycles by the dash-dotted parts (respectively with blue and green colour online). The plots are computed with XPPAUT and exported to Matlab for plotting

The oscillation amplitude changes slowly with for a fixed , that is the oscillation amplitude is robust to changes in stretch activation.

In Fig. 7 we plot the oscillation amplitude as a function of , for three selected values of , using XPPAUT. We see that as increases the amplitude decreases until the oscillations vanish close to , which agrees with the Hopf curve in Fig. 5. We also observe that for and 1, in Figs. 7a, b respectively, there are both stable and unstable limit cycles, and the right Hopf point is subcritical. Also, as increases, the -range of unstable limit cycles decreases until it vanishes; for (Fig. 7c) there are only stable limit cycles, and the right Hopf point has become supercritical. We see that as in the Atri model, the oscillation amplitude changes quite rapidly with in the mechanochemical system.

Fig. 7

Amplitude of calcium oscillations for the system (11)–(12) when , as is increased, for selected values of (computed with XPPAUT and exported to Matlab for plotting). The LHP and the RHP are indicated. The stable limit cycles are represented by dots and the unstable limit cycles by the dash-dotted parts (respectively with blue and green colour online): abc. As increases, for any fixed the amplitude decreases until it becomes zero

In Fig. 8 we plot the frequency of the limit cycles as increases, for three values of , using XPPAUT. For and , the frequency increases rapidly close to the LHP and the RHP and there is an ‘intermediate’ region where the frequency varies slowly with , as in the Atri system (see Fig. 3). The ‘intermediate’ region becomes smaller as increases, and for this region vanishes. We see that as increases the frequency of oscillations decreases overall.

Fig. 8

Frequency of calcium oscillations for the system (11)–(12) when as a function of and for (computed with XPPAUT and exported to Matlab for plotting). For and 1 there are stable limit cycles and unstable limit cycles, represented by dots and dash-dotted lines, respectively. For there are only stable limit cycles (blue colour online)

Summarising, for any value of and we can determine the range for oscillations using the parametric expressions (21) and (22), and then use XPPAUT (Ermentrout 2002) or other continuation software to obtain their amplitude and frequency.

In Fig. 9 we plot the evolution of c(t), solving (11)–(12) numerically, for and selected values of ; as expected from the bifurcation diagrams, the oscillations are suppressed when is sufficiently increased.

Fig. 9

Evolution of c(t) with time, solving the system (11)–(12) numerically, when , a (Atri model): limit cycles b: limit cycles with increased frequency and amplitude c: decaying solution; limit cycles (oscillations) disappear

Varying the cytosolic mechanical responsiveness factor

We now investigate if the Hopf curve changes qualitatively as the cytosol’s mechanical responsiveness factor, , varies. In Fig. 10a, using the parametric expressions (21)–(22) we plot the Hopf curves for increasing values of . We observe that the Hopf curve changes qualitatively; for it develops a cusp and for smaller values of there is a “bow-tie”. This geometrical change corresponds to yet another bifurcation, with as a bifurcation parameter1. However, as for , oscillations always vanish for a sufficiently large value of , .

Fig. 10

a Hopf curves for the system (11)–(12) and , (see legend) b The maximum value of for which oscillations are sustained, , decreases with . Both plots are drawn using the parametric expressions (21)–(22), in Mathematica. The horizontal line is the asymptote of the curve as . It represents the smallest possible in this system and since this is non-zero there are always be calcium oscillations for any value of a

We also observe that as increases, , the critical stretch activation value beyond which oscillations vanish, decreases, i.e. oscillations are sustained for a smaller range of values. To investigate this more systematically we have determined parametric expressions for and as functions of c, and we plot in Fig. 10b. We see that as increases, decreases monotonically, and hence oscillations are sustained for an increasingly smaller range of , which agrees with Fig. 10a. Also, since asymptotes to a positive value as for any , the system will always sustain some oscillations, irrespective of the value of . Therefore, we predict that for cytosols that are more responsive to calcium (higher ), oscillations vanish at a lower .To test this experimentally the responsiveness of the cytosol to calcium should be manipulated whilst monitoring whether oscillations appear. The contractility of the cytosol could be manipulated by inhibiting Myosin II contractility using the ROCK inhibitor (Y-27632).

However, since (21) does not depend on , is constant and not zero for any . Therefore, as we expect, is always required in order to obtain oscillations, for any and any but the minimum level of does not change with . Also, for fixed , as , the mechanical responsiveness factor of the cytosol, increases, the level required to induce oscillations also decreases. Additionally, for fixed , as increases reduces.

Summarising, we conclude that as the cytosol’s mechanical responsiveness increases a lower level of stretch activation is sufficient to sustain oscillations. Also, there will always be oscillations for some values of and when the contraction stress is modelled as a Hill function of order 1.

Hopf curves for a Hill function of order 2

It is instructive to investigate whether a different functional form of will change our conclusions. We thus model as a Hill function of order 2, , which models a cytosol which is less sensitive to calcium for low calcium levels than but which saturates faster. In Fig. 11 we plot the Hopf curves of the system (11)–(13) for increasing , the cytosolic mechanical responsiveness factor, using again the parametric expressions (21)–(22).

Fig. 11

Hopf curves for the system (11)–(13) when , for drawn using the parametric expressions (21)–(22), in Mathematica

Comparing Fig. 11 with Fig. 10a we see that the Hopf curves have the same qualitative behaviour for the two Hill functions. Oscillations can be sustained for any value of and they always vanish for a sufficiently large value of , Also, as in the Hill function of order 1, as increases decreases while is constant. Also, a cusp again develops but for the Hill function of order 2 the value of at which this occurs increases. We conclude that the conclusions are robust to the change of the Hill function. In future work Hill functions of higher order or other functional forms of T can be investigated.

Summary, conclusions and future research directions

A wealth of experimental evidence has accumulated which shows that many types of cells release calcium in response to mechanical stimuli but also that calcium release causes cells to contract. Therefore, studying this mechanochemical coupling is important for elucidating a wide range of body processes and diseases. In this work we have focused attention on embryogenesis, where the interplay of calcium and mechanics is shown to be essential in AC, an essential morphogenetic process which, if disrupted, leads to embryo abnormalities (Christodoulou and Skourides 2015).

We have presented a new analysis of experimental data that supports the existence of a two-way mechanochemical coupling between calcium signalling and contractions in embryonic epithelial cells involved in AC.

We have then developed a simple mechanochemical ODE model that consists of an ODE for , the cell apical dilation, derived consistently from a full, linear viscoelastic ansatz for a Kelvin-Voigt material, and two ODEs governing, respectively, the evolution of calcium and the proportion of active receptors. The two latter ODEs are based on the well-known, experimentally verified, Atri model for calcium dynamics (Atri et al. 1993). An important feature of our model is the two-way coupling between calcium and mechanics which was proposed for the first time in models by Murray (2001); Murray et al. (1988); Murray and Oster (1984) and Oster and Odell (1984). However, in those models hypothetical bistable calcium dynamics were considered whereas here we have updated those models with recent advances in calcium signalling, as encapsulated by the Atri model. We have also modelled the calcium-dependent contraction stress in the cytosol as a Hill function , since experiments indicate that the mechanical responsiveness of the cytosol to calcium saturates for high calcium levels.

The early mechanochemical models included an ad hoc stretch activation calcium flux, , in the calcium ODE. Here, we have also derived, for the first time, this “stretch-activation” flux as a “bottom-up” contribution from stretch sensitive calcium channels (SSCCs), thus expressing the parameter as a combination of the structural characteristics of an SSCC can also be thought of as a coupling parameter between calcium signalling and mechanics. Despite an extensive literature search we could not find experimental measurements for SSCCs; this could be a direction for future experiments.

For any, we have analytically identified the parameter regime in the – plane corresponding to calcium oscillations and applied this result in two illustrative examples, and . In both cases, as increases, the oscillations are eventually suppressed at a critical , —see, respectively, Figs. 10a and 11. The prediction is in agreement with experiments (Christodoulou and Skourides 2015) where a high, non-oscillatory calcium state is associated with a very high stress in the cytosol and continuous contraction (Figure 5D). This high-calcium, high-stress state is associated with failure of AC and consequently with defective tissue morphogenesis. This makes sense since calcium oscillations are the ‘information carrier’ in cells so we indeed expect that if they vanish the task at hand, in this case AC, will not be performed correctly. In summary, we have shown that there are scenarios where mechanical effects significantly affect calcium signalling and this is a key result of this work.

For we have also shown analytically that as , the mechanical responsiveness factor of the cytosol, increases, decreases but it never becomes zero (see Fig. 10b). This means that for any , there will always be a - region for which oscillations are sustained. Furthermore, for the illustrative example of we have determined numerically the oscillation amplitude and frequency as the bifurcation parameters and vary, using XPPAUT. We found that the behaviour is qualitatively similar to the Atri model (see Fig. 3) for lower values but that it changes for larger values (see Fig. 6). We found that, as increases the amplitude of oscillations decreases (see Fig. 7) but their frequency increases (see Fig. 8). More experiments could be undertaken to test these predictions.

In the experiments of Christodoulou and Skourides (2015) the calcium-induced stress saturates to a non-zero level as calcium levels increase but in other cell types it is possible that the cell can relax back to baseline stress and in such a case cannot be modelled as a Hill function. Experiments could be undertaken also in other calcium-induced mechanical processes to determine the appropriate form of and the model could then be modified appropriately.

Another approximation we have made is that the mechanical properties of the cell (Young’s modulus, Poisson ratio, viscosity) are constant. However, their values can change significantly with space and also with embryo stage (Brodland et al. 2006; Luby-Phelps 1999; Zhou et al. 2009). One of the next steps in the modelling would be to take these variations into account.

Due to the complexity of calcium signalling all models introduce approximations. One important approximation in this work is that we neglect stochastic effects, even though the opening and closing of receptors and of the SSCCs is a stochastic process. However, the deterministic models still have good predictive power, whilst being more amenable to analytical calculations (Cao et al. 2014; Thul 2014). A multitude of deterministic and stochastic calcium models have been developed (Atri et al. 1993; Goldberg et al. 2010; Gracheva et al. 2001; Sneyd et al. 1994, 1998; Timofeeva and Coombes 2003; Wilkins and Sneyd 1998); see also the comprehensive reviews (Rüdiger 2014; Sneyd and Tsaneva-Atanasova 2003; Thul 2014) and the books (Dupont et al. 2016; Keener and Sneyd 1998), among others. Future work could involve developing stochastic mechanochemical models.

The interplay of mechanics and calcium signalling in non-excitable cells is important in processes occurring not only in embryogenesis but also in wound healing and cancer, amongst many others, and more efforts should be devoted to developing appropriate mechanochemical calcium models that would help elucidate the currently many open questions. In this connection, the insights we have obtained from the simple model we have developed here are a first step in this direction. We will aim to extend our models to more realistic geometries. Moreover, we have fixed all parameters here, except , and ; and the variation of other parameter values may lead to other bifurcations and biologically relevant behaviours.

Finally, the newly discovered SSCCs merit much more experimental investigation and modelling; in this work we have adopted a simple model for their behaviour, assuming that they are quasisteady and also made restricting assumptions about their opening and closing rates. In further experimental work, the parameter values associated with SSCCs should be measured and perhaps more sophisticated models for SSCCs should be developed.