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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textit{IP}_3$$\end{document}IP3 concentration, and the “strength” of stretch activation, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document}λ. The calcium system (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda =0$$\end{document}λ=0, no mechanics) exhibits relaxation oscillations for a certain range of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textit{IP}_3$$\end{document}IP3 values. As \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document}λ is increased the range of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textit{IP}_3$$\end{document}IP3 values decreases and oscillations eventually vanish at a sufficiently high value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document}λ. This result agrees with experimental evidence in embryonic cells which also links the loss of calcium oscillations to embryo abnormalities. Furthermore, as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document}λ 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/Ca 2+ 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 Ca 2+ 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 pro-motes 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(Sanderson et al. , 1988Sanderson and Sleigh 1981), rat brain glial cells (Charles et al. 1993(Charles et al. , 1991(Charles et al. , 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 IP 3 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 IP 3 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 IP 3 is a control parameter and oscillations can be sustained at a constant IP 3 concentration. Oscillations exist for a window of IP 3 values; the oscillations are excited at a threshold IP 3 value and they vanish at a suffuciently high IP 3 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 IP 3 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. 1990Sanderson et al. , 1988Sanderson 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 calciumelicited 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.
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.
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 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 t = 0 and the average calcium oscillation frequency was calculated using a 10-minute window (i.e. calcium oscillations for each cell were monitored between t = 0 and t = 10 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.) a b 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 form Here, c is the cytosolic calcium concentration, h is the fraction of IP 3 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), J ER models the flux of calcium from the ER into the cytosol through the IP 3 receptors, μ( p) = p/(k μ + p) is the fraction of IP 3 receptors that have bound IP 3 and is an increasing function of p, the IP 3 concentration. The constant k f denotes the calcium flux when all IP 3 receptors are open and activated, and b represents a basal current through the IP 3 channel. J pump represents the calcium flux pumped out of the cytosol where γ is the maximum rate of pumping of calcium from the cytosol and k γ is the calcium concentration at which the rate of pumping from the cytosol is at half-maximum. J leak models the calcium flux leaking into the cytosol from outside the cell. Note that from now on we will neglect J leak as this is a good approximation for the embryonic epithelial cells we consider.
J SSCC 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 IP 3 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, τ h > 1 in ODE (2). In ODE (3) T D (c) is a contraction stress that expresses how the stress in the cell depends on the cytosolic calcium level. The constants ξ 1 , ξ 2 are, respectively, the shear and bulk viscosities of the cytosol and the constants E = E/(1 + ν) and ν = ν/(1−2ν), 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 J SSCC 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 IP 3 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 , mammalian spermatozoa (Olson et al. 2010), and keratinocytes (Kobayashi et al. 2014(Kobayashi et al. , 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, Sθ , 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 C SSCC the proportion of channels in the closed state, and by O SSCC the proportion of SSCCs in the open state. The calcium flux due to the SSCCs is proportional to the number of open channels so J SSCC = K SSCC O SSCC , where K SSCC is the maximum calcium flux rate going through the SSCCs. As in Vainio et al. (2015), we propose that the evolution of O SSCC is governed by the ODE where k F is the forward rate constant and k B is the backward rate constant. We assume here that O SSCC is quasi-steady, i.e. O SSCC 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 k F k B is of order 1 at most. We obtain Therefore, we have derived, for the first time, an expression for S as a combination of K SSCC , k F and k B , 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 oe is the stress tensor, e = 1 2 (∇u + ∇u T ) is the strain tensor, u the displacement vector, θ = ∇.u is the dilation/compression of the material, and I is the unit tensor. (The subscript t here denotes a partial derivative with respect to time.) E = E 1+v and v = v 1−2v where E and v are the Young's modulus and the Poisson's ratio, respectively. T D (c) is the contraction stress which depends on the cytosolic calcium (Scholey et al. 1980). In one spatial dimension e = e = θ = ∂u ∂ x and therefore (9) becomes, upon integrating with respect to x, The constant of integration A = 0 since when c = 0, T D = 0 , θ = 0 and θ t = 0. Hence, we obtain ODE (3).

Nondimensionalised model
We nondimensionalise the mechanochemical model using c = k 1c and t = τ ht . Dropping bars for notational convenience we obtain In (11) (13), , and in (12) K 2 = k 2 /k 1 . Using the parameter values of Atri et al. (1993) (see Appendix, Table 1), we obtain K 2 = 1, Γ = 40/7, and K = 1/7. Also, taking values of E, ν and of the viscosity from Zhou et al. (2009) (E = 8.5 Pa, ν = 0.4 and ξ 1 + ξ 2 = 100 Pa.s) we find that k θ is 0.4. For simplicity, and since the parameter values for the calcium dynamics are approximate, is nondimensional, and we also fix τ h (ξ 1 +ξ 2 ) T 0D = 1. 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.

The bifurcation diagrams of the Atri model (no mechanics)
The nondimensional Atri system is In 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 0.289 ≤ μ ≤ 0.495, 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 μ = 0.289 and μ = 0.495. 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 μ = 0.511, 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).

Linear stability analysis of the mechanochemical model
The steady states of the system (11)-(12) satisfy  (14)- (15), as μ (IP 3 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 For anyT (c), using (16), we can easily plot the steady states as a function of μ and λ. The Jacobian of (11)-(12) is given by and the characteristic polynomial is conveniently factorised as where ω represents the eigenvalues. As one eigenvalue is always equal to -1, the bifurcations of the system can be studied through the quadratic To identify the μ-λ parameter range sustaining oscillations we seek the Hopf bifurcations, which satisfy Tr(M 2 ) = 0, Det(M 2 ) > 0, Discr(M 2 ) < 0, where Setting and substituting in (16) we obtain Hence, we can easily obtain the Hopf curve, for anyT (c) 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 Det(M 2 ) = 0 Equations (23) and (16) constitute a linear system for μ and λ, so we again easily derive parametric expressions for μ(c) and λ(c). Similarly, to determine parametric expressions for the curve on which Discr(M 2 ) changes sign we set Discr(M 2 ) = 0 which 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 anyT (c).
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.

Hopf curves
We assume that the calcium-induced stressT (c) is the Hill function assuming that when the calcium level increases sufficiently the stress saturates to a constant value, T s = 1. 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,T = 0 when c = 0, i.e. we assume no contraction stress without calcium.T (0) = α is the rate of increase ofT at c = 0 and 1/α is the scale of 'ascent' to the saturation level T s . Therefore, we can call α the 'mechanical responsiveness factor' of the cytoskeleton to calcium. ChoosingT (c) = 10c/(1 + 10c) 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 λ < 4 the equilibrium curve is qualitatively similar to that of the Atri model (see Fig. 3a) but at λ = 4 the curve changes qualitatively and a second non-zero steady state exists for 4 < λ < 40/7, and a part of the curve corresponds to negative values of μ (see Appendix A.3 for details). For λ > 40/7 no steady state exists for positive μ and hence λ ≤ 40/7 is the biologically relevant range of the model, for α = 10 (see Appendix A.3).
In Fig. 5 we plot the Hopf curve and the fold curve. We observe the following: (i) for λ = 0 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, μ min . (iii) The oscillations are suppressed for a critical maximum value of λ, λ max , and the system is in a high calcium state. Overall, we conclude the following from the Hopf curve: -for low IP 3 values the Atri system does not sustain oscillations but there are two possibilities for the mechanochemical model as λ increases: • for μ < μ min no increase in λ will ever elicit oscillations. • for 0.203 = μ min < μ < 0.289 when λ reaches a certain value, λ OSC , oscillations are elicited, and λ OSC decreases as μ approaches 0.289. The oscillations vanish at a larger value of λ. -for IP 3 values for which the Atri system sustains oscillations (0.289 < μ < 0.495) in the mechanochemical model oscillations eventually vanish at a critical λ. This critical λ decreases monotonically as μ increases towards 0.495. -for high IP 3 values (μ ≥ 0.495) no oscillations are sustained in the Atri system and a further increase in λ will never elicit oscillations.
Therefore, for fixed cytoskeletal mechanical responsiveness factor, α = 10, and for fixed parameter values as in Atri et al. (1993) a range of behaviours emerge as μ and λ vary: at low IP 3 levels that do not elicit oscillations in the Atri system mechanical effects can elicit oscillations, for intermediate IP 3 levels that do sustain oscillations in the Atri system increasing mechanical effects always leads to the oscillations vanishing, and for high IP 3 levels that cannot sustain oscillations in the Atri system no amount of stretch activation can ever 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 IP 3 ) directly affect cell contractility (Christodoulou and Skourides 2015). At low levels of IP 3 and hence low levels of calcium, cells are not able to contract and therefore AC does not take place. At a threshold IP 3 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 IP 3 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 IP 3 concentrations. As λ increases this range decreases and eventually vanishes at λ ≈ 0.83, 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 ofT (c) 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 aT (c) that saturates. In other cell types it is possible that the cell can relax back to baseline stress and in such a caseT (c) would not be described by a Hill function, and more experiments should be undertaken to determine the appropriate form ofT (c).
In Fig. 6 we plot the oscillation amplitude as a function of λ, for two selected values of μ, using XPPAUT. For μ = 0.25 (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 μ = 0.3 (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 μ = 0.4 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 λ.
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 λ = λ max = 1.69, which agrees with the Hopf curve in Fig. 5. We also observe that for λ = 0.5 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 λ = 1.5 (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.
In Fig. 8 we plot the frequency of the limit cycles as μ increases, for three values of λ, using XPPAUT. For λ = 0.5 and λ = 1, 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 λ = 1.5 this region vanishes. We see that as λ increases the frequency of oscillations decreases overall.
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 μ = 0.3 and selected values of λ; as expected from the bifurcation diagrams, the oscillations are suppressed when λ is sufficiently increased.

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 α = 1, 2, 10, 100. We observe that the Hopf curve changes qualitatively; for α ≈ 2 it develops a cusp and for smaller values of α there is a "bow-tie". This geometrical change corresponds to yet another Amplitude of calcium oscillations for the system (11)-(12) whenT (c) = 10c 1+10c , as λ is increased, for: a μ = 0.25 b μ = 0.3. 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 bifurcation, with α as a bifurcation parameter 1 . However, as for α = 10, oscillations always vanish for a sufficiently large value of λ, λ max . c λ = 1.5 Fig. 7 Amplitude of calcium oscillations for the system (11)-(12) whenT (c) = 10c 1+10c , 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): a λ = 0.5 b λ = 1 c λ = 1.5. As λ increases, for any fixed μ the amplitude decreases until it becomes zero 1+10c as a function of μ and for λ = 0.5, 1, 1.5 (computed with XPPAUT and exported to Matlab for plotting). For λ = 0.5 and 1 there are stable limit cycles and unstable limit cycles, represented by dots and dash-dotted lines, respectively. For λ = 1.5 there are only stable limit cycles (blue colour online) We also observe that as α increases, λ max , 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 λ max and α as functions of c, and we plot λ max (α) in Fig. 10b. We see that as α increases, λ max decreases monotonically, and hence oscillations are sustained for an increasingly smaller range of λ, which agrees with Fig. 10a. Also, since λ max (α) asymptotes to a positive value as α → ∞ for anyT (c) = αc 1+αc , 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 λ max .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 onT (c), μ min is constant and not zero for any α. Therefore, as we expect, IP 3 is always required in order to obtain oscillations, for any λ and any α but the minimum level of IP 3 does not change with α. Also, for fixed λ, as α, the mechanical responsiveness factor of the cytosol, increases, the IP 3 level required to induce oscillations also decreases. Additionally, for fixed μ, as α increases λ max 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.
Evolution of c(t) with time, solving the system (11)-(12) numerically, whenT (c) = 10c 1+10c , μ = 0.289 a λ = 0 (Atri model): limit cycles b λ = 1: limit cycles with increased frequency and amplitude c λ = 3: decaying solution; limit cycles (oscillations) disappear a b Fig. 10 a Hopf curves for the system (11)-(12) andT (c) = αc 1+αc , α = 1, 2, 10, 100 (see legend) b The maximum value of λ for which oscillations are sustained, λ max , decreases with α. Both plots are drawn using the parametric expressions (21)-(22), in Mathematica. The horizontal line is the asymptote of the λ max curve as a → ∞. It represents the smallest possible λ max in this system and since this is non-zero there are always be calcium oscillations for any value of a

Hopf curves forT(c) a Hill function of order 2
It is instructive to investigate whether a different functional form ofT will change our conclusions. We thus modelT (c) as a Hill function of order 2,T (c) = αc 2 1+αc 2 , which models a cytosol which is less sensitive to calcium for low calcium levels than T (c) = αc 1+αc 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).
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 λ max decreases while μ min 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 IP 3 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 calciumdependent contraction stress in the cytosol as a Hill functionT (c), 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 anyT (c), we have analytically identified the parameter regime in the μ-λ plane corresponding to calcium oscillations and applied this result in two illustrative examples,T (c) = αc/(1+αc) andT (c) = αc 2 /(1+αc 2 ). In both cases, as λ increases, the oscillations are eventually suppressed at a critical λ, λ max -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.
ForT (c) = αc/(1 + αc) we have also shown analytically that as α, the mechanical responsiveness factor of the cytosol, increases, λ max 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 ofT (c) = 10c/(1+10c) 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 caseT (c) cannot be modelled as a Hill function. Experiments could be undertaken also in other calcium-induced mechanical processes to determine the appropriate form ofT (c) 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 IP 3 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. 1994Sneyd et al. , 1998Timofeeva 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. calcium level of cells induced to undergo AC at gastrula stage in order to decouple AC from other morphogenetic movements, like mediolateral junction shrinkage and convergent extension, which take place in later embryogenic stages and which would also influence the cell shape and surface area.
Explanations about Figure 2: (a) Plot of surface area (μm 2 ) and calcium oscillation frequency over time from 10 neural plate cells of stage 16 Xenopus embryo using the mem-GFP +GECO RED sensor molecule. The average surface area of each cell was evaluated for four time intervals; 0-10, 10-20, 20-30 and 30-40 min. For normalization for each cell, the average surface area in each time period was divided by the average surface area of the first period (0-10 min). The calcium oscillation frequency in each cell was calculated by counting the number of calcium oscillations in each time interval. This value was then divided by 10 since there are 10 minutes in each time interval. (b) Plot of average calcium oscillation amplitude of 10 neural plate cells (same cells as in (a)). The signal intensity of the non-ratiometric calcium sensor (GECO-RED) was measured per calcium oscillation in each of the cells over time. For normalization, the values were then divided by the highest intensity value. The average value in each of the four time intervals was plotted.

A.2.1 Linear stability analysis
The steady states (S.S.) of (14)-(15) are the intersections of the nullclines of the system. Setting we obtain which can be cast as a quartic in c. (Note that (16) reduces to (28) for λ = 0, as expected.) In Figure 12 we plot the equilibrium curve (28) in order to visualise the number of steady states and the corresponding value(s) of c, as μ is increased. The qualitative behaviour of the solutions of the system can be determined by plotting the nullclines (26) and (27). When the nullclines cross the system has a steady state, and when they touch the system has a double (degenerate) steady state. Nullcline (26) passes through the origin of the (c,h) plane, has a maximum at h = h M and saturates to the constant value h = Γ μK 1 as c → ∞. h M can be found analytically by differentiating (26):

Fig. 12
The steady states of (14)-(15) as a function of the bifurcation parameter μ, as obtained using the analytical expression (28). As μ increases, from small to large, there is one steady state, a double (degenerate) steady state at μ 1 = 0.28814, then three steady states, then a double (degenerate) steady state at μ 2 = 0.28925, and for values of μ larger than 0.28925 one steady state and setting dh/dc = 0, which leads to a quadratic equation for c. Rearranging, and since c > 0, we discard the negative root, obtaining and, hence, substituting (30) in (26) Therefore, h M scales with Γ /(μK 1 ) and depends on the parameters K and b in a much more complicated manner. For the parameter values from Atri et al. (1993) we have c M = 0.169 and h M = 0.279/μ. Nullcline (27) is a decreasing function of c; it has a maximum at (0,1) and tends to 0 as c → ∞. For μ 1 = 0.28814 and μ 2 = 0.28925 the nullclines touch and there is a double steady state; for values of μ < μ 1 and μ > μ 2 there is one S.S. and there are three S.S. for μ 1 < μ < μ 2 . (μ 1 and μ 2 are obtained by differentiating (28) and finding the roots of dμ dc = 0.) Note that we present the values of μ with five decimal places because the bifurcation analysis depends sensitively on μ, as we will see later.
We then linearise the system near the steady states. We determine the Trace (Tr), Determinant (Det) and Discriminant (Discr) of the Jacobian of the system as follows Taking the parameter values of Atri et al. (1993) we identify the bifurcations of the system as μ increases by determining where the Tr, Det, and Discr change sign. We find a richer bifucation structure as μ increases. This behaviour was not analysed in such detail in Atri et al. (1993) or in later literature. Given the very sensitive dependence on precise values of μ, these details are probably of more mathematical interest than of biological significance. The parameter values are summarised in Table 1.
• 0 < μ < 0.27828: one stable node. • μ = 0.27828: the stable node becomes a stable spiral (bifurcation Discr = 0) • μ = 0.28814: Stable spiral present. Also, a saddle and an unstable node (UN) emerge (bifurcation Det = 0, fold point) • μ = 0.28900: the stable spiral becomes an unstable spiral. The other two S.S. are still a saddle and an unstable node. (Tr=0, Hopf bifurcation) • μ = 0.28924 the unstable spiral becomes an unstable node, and we have two unstable nodes and a saddle (Discr = 0) • μ = 0.28925: one unstable node (Det = 0, fold point) • μ = 0.28950: the unstable node becomes an unstable spiral (Discr = 0) • μ = 0.49500: the unstable spiral becomes a stable spiral. (Tr=0, Hopf bifurcation) From the regimes identified above we are particularly interested in the regime of relaxation oscillations, since their amplitude and/or frequency encodes the information in calcium signals. Relaxation oscillations are sustained for 0.28900 ≤ μ ≤ 0.49500 since at μ = 0.28900 a Hopf bifurcation (HB) arises, the stable spiral becomes unstable, and we expect relaxation oscillations (limit cycles) in the nonlinear system. As μ increases the unstable spiral becomes a stable spiral close to μ = 0.495000 and the limit cycles eventually vanish.

A.3 Linear stability analysis of the mechanochemical model with no IP 3
Here we analyse the case μ = 0. Biologically, this corresponds to treating the cell with thapsigargin so that all calcium is depleted from the ER and there is no flux from the ER. The model (11)-(13) simplifies to Equation (34) is decoupled from Eqs. (32) and (33), and we can thus continue with a two-dimensional analysis for Eqs. (32) Therefore, for anyT (c), Discr(M 1 ) > 0, Tr(M 1 ) < 0 always, and Det(M 1 ) can be negative or positive, the steady states can only be stable nodes or saddles, and oscillations cannot be sustained. For some choices ofT a non-zero steady state may exist and this means biologically that even without calcium flux from the ER a nonzero calcium concentration can be sustained in the cytosol due to the stress-induced calcium release.
We can easily show that the steady state (0,0) loses its stability and becomes a saddle when the second stable S.S. emerges (as a stable node). This means that there is a range of λ values for which the system sustains a non-zero calcium concentration even without a CICR flux. For even larger values of λ there is no steady state and the model ceases to be biologically relevant.