A simple mechanochemical model for calcium signalling in embryonic epithelial cells

TL;DR

This paper introduces a mechanochemical model linking calcium dynamics and mechanical forces in embryonic cells, supported by experimental data, revealing how stretch-sensitive channels influence calcium oscillations during development.

Contribution

It presents a novel, simplified mechanochemical model incorporating stretch-activated calcium channels and analyzes its bifurcation structure, advancing understanding of calcium-mechanics coupling in embryogenesis.

Findings

Calcium oscillations depend on IP3 levels and stretch activation strength.

Increasing stretch activation reduces oscillation range and amplitude.

High stretch activation suppresses calcium oscillations, correlating with embryo abnormalities.

Abstract

Calcium (Ca2+) signalling is one of the most important mechanisms of information propagation in the body. In embryogenesis the interplay between Ca2+ 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 Ca2+ signals and contractions are coupled via a two-way mechanochemical coupling. We present a new analysis of experimental data that supports the existence of this coupling during Apical Constriction in Neural Tube Closure. We then propose a mechanochemical model, building on early models that couple Ca2+ dynamics to cell mechanics and replace the bistable Ca2+ release with modern, experimentally validated Ca2+ dynamics. We assume that the cell is a linear viscoelastic material and model the Ca2+-induced contraction stress with a Hill function saturating at high Ca2+ levels. We also express, for the first time, the "stretch-activation" Ca2+ flux in the early mechanochemical models as a bottom-up contribution from stretch-sensitive Ca2+ channels on the cell membrane. We reduce the model to three ordinary differential equations and analyse its bifurcation structure semi-analytically as the $IP_3$ concentration, and the "strength" of stretch activation, $\lambda$ vary. The Ca2+ system ($\lambda=0$, no mechanics) exhibits relaxation oscillations for a certain range of $IP_3$ values. As $\lambda$ is increased the range of $IP_3$ values decreases, the oscillation amplitude decreases and the frequency increases. Oscillations vanish for a sufficiently high value of $\lambda$. These results agree with experiments in embryonic cells that also link the loss of Ca2+ oscillations to embryo abnormalities. The work addresses a very important and understudied question on the coupling of chemical and mechanical signalling in embryogenesis.

Tables

Table 1. Table 1 : Parameter values, taken from [ 3 ]

Parameter values
Parameter	Value
b	0.111
$k_1$	$0.7\mu M$
$k_f$	$16.2\mu M/s$
$k_{\mu }$	$0.7\mu M$
$\gamma$	$2\mu M/s$
$k_{\gamma }$	$0.1\mu M$
$\beta$	$0-0.02\mu M/s$
$k_2$	$0.7\mu M$
${\tau }_h$	$2s$

Equations

$$\frac{dc}{dt}$$

$${\tau }_h\frac{dh}{dt}$$

$$\frac{d\theta }{dt}$$

$$\text{where }{J}_{\mathrm{ER}}$$

$$\frac{dc}{dt}$$

$${\tau }_h\frac{dh}{dt}$$

$$\frac{d(O_{\mathrm{SSCC}})}{dt}=k_F\theta -(k_F\theta +k_B)O_{\mathrm{SSCC}},$$

$$O_{\mathrm{SSCC}}\approx \frac{k_F\theta }{k_F\theta +k_B}.$$

$$O_{\mathrm{SSCC}}$$

$$\nabla .\sigma =0\Rightarrow \nabla .(\underset{\text{viscous stress}}{\underbrace{{\xi }_1{\mathbf{e}}_t+{\xi }_2{\theta }_t\mathbf{I}}}+\underset{\text{elastic stress}}{\underbrace{E^{\prime }(\mathbf{e}+{\nu }^{\prime }\theta \mathbf{I})}}-\underset{\text{contraction stress}}{\underbrace{T_D(c)\mathbf{I})}}=0,$$

$$({\xi }_1+{\xi }_2){\theta }_t+E^{\prime }(1+{\nu }^{\prime })\theta -T_D(c))=A.$$

$$\frac{dc}{dt}$$

$$\frac{dh}{dt}$$

$$\frac{d\theta }{dt}$$

$$\frac{dc}{dt}$$

$$\frac{dh}{dt}$$

$$\mu K_1\frac{1}{1+c^2}\frac{b+c}{1+c}-\frac{\Gamma c}{K+c}+\lambda \hat{T}(c)=0.$$

$$M_1=\left[\begin{array}{ccc}{R}_{1c}& \lambda & R_{1h}\\ {\hat{T}}^{\prime }(c)& -1& 0\\ R_{3c}& 0& -1\end{array}\right],$$

$$(1+\omega )(\lambda {\hat{T}}^{\prime }(c)+(R_{1c}-\omega )(1+\omega )+R_{1h}{R}_{3c})=0,$$

$${\omega }^2-\omega (R_{1c}-1)-R_{1c}-R_{1h}{R}_{3c}-\lambda T^{\prime }(c)=0.$$

$$M_2=\left[\begin{array}{cc}{R}_{1c}& \lambda \\ {\hat{T}}^{\prime }(c)& -1\end{array}\right].$$

$$\text{Setting Tr}(M_2)=0⟹\mu (c)$$

$$\lambda (c)=\frac{1}{\hat{T}(c)}\left(\frac{\Gamma c}{K+c}-\frac{(b+c)(1+c)}{1-b}\left(1+\frac{\Gamma K}{(K+c{)}^2}\right)\right).$$

$$⟹\mu \frac{K_1}{(1+c)(1+c^2)}\left(\frac{1-b}{1+c}-2c(b+c)\right)+\lambda T^{\prime }(c)=\frac{\Gamma K}{(K+c{)}^2}.$$

$$⟹(R_{1c}+1{)}^2+4R_{1h}{R}_{3c}+\lambda {\hat{T}}^{\prime }(c)=0,$$

$$\hat{T}(c)=\frac{\alpha c}{1+\alpha c},\alpha >0,$$

$$F=0⟹h$$

$$G=0⟹h$$

$$\mu K_1\frac{1}{1+c^2}\frac{b+c}{1+c}-\frac{\Gamma c}{K+c}=0,$$

$$\frac{dh}{dc}$$

$$c_M(b,K)$$

$$h_M$$

$$\mathrm{Tr}=F_c+G_h=F_c-1=\frac{\Gamma }{(K+c)}\left(-\frac{K}{K+c}+\frac{(1-b)\Gamma c}{(1+c)(b+c)}\right)-1$$

$$\mathrm{Det}=F_c{G}_h-F_h{G}_c=-\frac{\Gamma }{(K+c)}\left(-\frac{K}{K+c}+\frac{(1-b)\Gamma c}{(1+c)(b+c)}\right)+\frac{2\Gamma c^2}{(1+c^2)(K+c)}$$

$$\mathrm{Discr}=(\mathrm{Tr}{)}^2-4Det.$$

$$\frac{dc}{dt}$$

Full text

††journal: Journal of Mathematical Biology

1 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 berridge2000 ; brini2009 ; dupont2016 . 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 lecuit2007cell . The impact of these forces on morphogenesis is directly linked to calcium signalling hunter2014ion . In general, how the mechanics of the cell and tissue are regulated and coupled to the cellular biochemical responseto 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 antunes2013coordinated ; herrgen2014calcium and cancer basson2015 .

Calcium plays a crucial role in every stage of embryonic development starting with fast calcium waves during fertilization deguchi2000spatiotemporal to calcium waves involved in convergent extension movements during gastrulation wallingford2001calcium , to calcium transients regulating neural tube closure christodoulou2015cell , morphological patterning in the brain webb2007ca2+ ; sahu2017calcium and apical-basal cell thinning in the enveloping layer cells zhang2011necessary , either in the form of calcium waves or through Wnt/$\rm Ca^{2+}$ signalling kuhl2000ca2+ ; kuhl2000wnt ; slusarski1997interaction ; slusarski1997modulation ; wallingford2001calcium ; herrgen2014calcium ; hunter2014ion ; christodoulou2015cell ; Narciso2017release ; suzuki2017distinct . Crucially, pharmacologically inhibiting calcium has been shown to lead to embryo defects smedley1986calcium ; wallingford2001calcium ; christodoulou2015cell .

In many experiments actomyosin-based contractions have been documented in response to calcium release in both embryonic and cultured cells wallingford2001calcium ; hunter2014ion ; herrgen2014calcium ; christodoulou2015cell ; suzuki2017distinct and it has become clear that calcium is responsible for contractions in both muscle and non-muscle cells, albeit through different mechanisms cooper2000 . Cell contraction in striated muscle is mediated by the binding of $\rm 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 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 scholey1980regulation . 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 cooper2000 .

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 vijayraghavan2017mechanics . 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 sawyer2010apical ; rohrschneider2009polarity ; christodoulou2015cell . .

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 sanderson1981 ; sanderson1988 ; sanderson1990 , rat brain glial cells charles1991 ; charles1992 ; charles1993 , keratinocytes tsutsumi2009 , developing epithelial cells in Drosophila wing discs Narciso2017release and many other cell types young1999 ; beraheiter2005 ; yang2009 ; tsutsumi2009 . Thus, different types of mechanical stimuli, from shear stress to direct mechanical stimulation, can elicit calcium elevation (although 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 christodoulou2015cell ; we present this analysis in detail in Section 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 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 mechanohemicall feedback loop is in play solon2009pulsed ; saravanan2013local . 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 yao2014force ; petridou2016ligand ; delmas2013mechano .

Given the broad range of critical biological processes involving calcium signaling 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 murray1984generation ; oster1984mechanochemistry ; murray1988 ; murray2001 . 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, and, so, our mechanochemical model includes instead the calcium dynamics of the experimentally verified model in atri1993single , which captures the experimentally observed Calcium-Induced-Calcium Release process and the dynamics of the $IP_{3}$ receptors on the ER. 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 warren2010 , for keratinocytes kobayashi2014 , for white blood cells yao2016 , and for retinal pigment epithelial cells vainio2015computational . 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 Endoplasmic Reticulum (ER), or the Sarcoplasmic Reticulum (SR), according to the well-established nonlinear feedback mechanism of Calcium Induced Calcium Release (CICR)dupont2016 . 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 dupont2016 . 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 in atri1993single is an established Class I model, validated with experimental findings estrada2016cellular . (We will call this model the ‘Atri model’ from now on.) It also has a mathematical structure that allows us to analyse 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 liRinzel1994 ; 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) dupont2016 ; arnadottir2010eukaryotic ; moore2010stretchy ; hamill2006twenty , in this way linking the channel scale with the whole cell scale.

The paper is organised as follows. In Section 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 Section 3 we develop a new mechanochemical model which captures the key ingredients of the two-way coupling. In Section 4 we analyse the model. In Section 4.1 we briefly revisit the analysis of the Atri model atri1993single and show the bifurcation diagrams for the amplitude and frequency of calcium oscillations. In Section 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 Section 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 Section 5.1.2 we study the amplitude and frequency of the calcium oscillations. In Section 5.1.3 we investigate the bifurcation diagrams as the mechanical responsiveness of the cytosol to calcium varies. In Section 5.2 we consider a Hill function of order $2$ and we again identify the parameter range for oscillations. In Section 6 we summarise our conclusions and discuss further research directions.

2 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 sanderson1981 ; sanderson1988 ; sanderson1990 ; charles1991 ; young1999 ; beraheiter2005 ; tsutsumi2009 ; Narciso2017release and that, in turn, contraction of the cytosol is elicited by calcium wallingford2001calcium ; hunter2014ion ; herrgen2014calcium ; christodoulou2015cell ; suzuki2017distinct . Calcium signaling 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 Apical Constriction (AC), where cells contract significantly, such a feedback loop should also be at play martin2014apical ; in this work we present a new analysis of experimental data in christodoulou2015cell which supports this. AC is a calcium-driven morphogenetic movement of epithelial tissues, central in the embryogenesis of both vertebrates and invertebrates vijayraghavan2017mechanics . The apical domain of epithelial cells constricts the apical surface area, contracting, and this leads to changes in the cell geometry that drive tissue bending; in christodoulou2015cell the formation of the neural tube in Xenopus frogs is studied and in solon2009pulsed dorsal closure in Drosophila.

In solon2009pulsed the constriction of mutants that exhibit disrupted myosin activation rescues apical myosin accumulation, suggesting that mechanically stimulating the cell can elicit contractions pouille2009mechanical . 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 saravanan2013local and that this process is driven by calcium pouille2009mechanical ; saravanan2013local ; hunter2014ion . 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 hunter2014ion ; similarly during wound healing antunes2013coordinated .

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 christodoulou2015cell . Here, in order to investigate in detail the relationship between calcium, contraction and mechanical forces we present a new analysis of previously collected data christodoulou2015cell . For a single embryonic epithelial cell (in a tissue), we plot its apical surface area and calcium level over time in Figure 1 and we see that both oscillate, with approximately the same frequency and that the calcium pulse precedes the contraction by 30-40 seconds. (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 Figure 2(a) 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 Figure 2(b), 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 to 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 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.

3 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 cooper2000 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 murray1984generation ; oster1984mechanochemistry ; murray1988 ; murray2001 and here we update those models by replacing the hypothetical bistable calcium dynamics with the experimentally verified calcium dynamics in atri1993single . We also replace the ad hoc stretch activation calcium flux in murray2001 with a “bottom-up” calcium release through the SSCCs, thus linking the channels’ characteristics to the whole cell scale. The model takes the form

[TABLE]

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 $\theta$ is the dilation/compression of the apical surface area of the cell. In ODE (1), $J_{\rm ER}$ models the flux of calcium from the ER into the cytosol through the $IP_{3}$ receptors, $\mu(p)=p/(k_{\mu}+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_{\rm pump}$ represents the calcium flux pumped out of the cytosol where $\gamma$ is the maximum rate of pumping of calcium from the cytosol and $k_{\gamma}$ is the calcium concentration at which the rate of pumping from the cytosol is at half-maximum. $J_{\rm leak}$ models the calcium flux leaking into the cytosol from outside the cell. Note that from now on we will neglect $J_{\rm leak}$ as this is a good approximation for the embryonic epithelial cells we consider.

$J_{\rm 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 dupont2016 ; arnadottir2010eukaryotic ; moore2010stretchy ; hamill2006twenty . The constant $S$ represents the ‘strength’ of stretch activation. In Section 3.1 we will derive a relationship for $S$ as a function of the characteristics of a SSCC.

The inactivation of the $IP_{3}$ receptors by calcium acts on a slower timescale compared to activation dupont2016 and so the time constant for the dynamics of $h$, $\tau_{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 $\xi_{1},\xi_{2}$ are, respectively, the shear and bulk viscosities of the cytosol and the constants $E^{\prime}=E/(1+\nu)$ and $\nu^{\prime}=\nu/(1-2\nu)$, where $E$ and $\nu$ are, respectively, the Young’s modulus and the Poisson ratio.

Our mechanochemical model is also an extension of the Atri model,

[TABLE]

since ODE (1) is ODE (4) but with $J_{\rm SSCC}$ added to the right hand side as an extra source term. The detailed derivation of the Atri model is presented in atri1993single , where it was initially formulated, and the reader is referred there for more details. The parameter values, which we take from atri1993single , 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 keener1998 ; dupont2016 . 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 wilkins1998 , mammalian spermatozoa olson2010 , and keratinocytes kobayashi2014 ; kobayashi2016 . In addition, modified Atri models have been developed in shi2008 ; harvey2011 ; liu2016 . In estrada2016cellular 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 liRinzel1994 . Also, the Atri model has a mathematical structure that allows us to perform a large part of our study analytically. The Atri system 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.

3.1 Stretch-activation calcium flux

In the early mechanochemical models murray2001 the stretch-activation flux $S\theta$ 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 vainio2015computational 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_{\rm SSCC}$ the proportion of channels in the closed state, and by $O_{\rm 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_{\rm SSCC}=K_{\rm SSCC}O_{\rm SSCC}$, where $K_{\rm SSCC}$ is the maximum calcium flux rate going through the SSCCs. As in vainio2015computational , we propose that the evolution of $O_{\rm SSCC}$ is governed by the ODE

[TABLE]

where $k_{F}$ is the forward rate constant and $k_{B}$ is the backward rate constant. We assume here that $O_{\rm SSCC}$ is quasi-steady, i.e. $O_{\rm SSCC}$ remains approximately constant as calcium rapidly evolves. This is a reasonable approximation, as discussed in Section 2.6 of dupont2016 . Therefore,

[TABLE]

We linearise (7) since $\theta$ is small for a linear viscoelastic medium and under the additional assumption that $\frac{k_{F}}{k_{B}}$ is of order 1 at most. We obtain

[TABLE]

Therefore, we have derived, for the first time, an expression for $S$ as a combination of $K_{\rm SSCC},k_{F}$ and $k_{B}$, linking in this way the characteristics of an SSCC to the macroscopic stretch-activation calcium flux.

3.2 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 dassow2010surprisingly show it is reasonable. For a Kelvin-Voigt, linear viscoelastic material sustaining calcium-induced contractions the force balance equation can be written as follows landau1986 ; murray2001 :

[TABLE]

where $\mathbf{\sigma}$ is the stress tensor, $\mathbf{e}=\frac{1}{2}(\nabla\mathbf{u}+\nabla\mathbf{u}^{T})$ is the strain tensor, $\mathbf{u}$ the displacement vector, $\theta=\nabla.\mathbf{u}$ is the dilation/compression of the material, and $\mathbf{I}$ is the unit tensor. $T_{D}(c)$ is the contraction stress which depends on the cytosolic calcium scholey1980regulation . In one spatial dimension $\mathbf{e}=e=\theta=\frac{\partial u}{\partial x}$ and therefore (9) becomes, upon integrating with respect to $x$,

[TABLE]

The constant of integration $A=0$ since when $c=0$, $T_{D}=0$ , $\theta=0$ and $\theta_{t}=0$. Hence, we obtain ODE (3).

3.3 Nondimensionalised model

We nondimensionalise the mechanochemical model using $c=k_{1}\bar{c}$ and $t=\tau_{h}\bar{t}$. Dropping bars for notational convenience we obtain

[TABLE]

In (11) $K_{1}=k_{f}\tau_{h}/k_{1}$, $\Gamma=\gamma\tau_{h}/k_{1}$, $K=k_{\gamma}/k_{1}$, and $\lambda=\tau_{h}S/k_{1}$. In (13), $\displaystyle{k_{\theta}=\frac{\tau_{h}E^{\prime}(1+\nu^{\prime})}{(\xi_{1}+\xi_{2})}}$ and $T(c)=\frac{\tau_{h}}{(\xi_{1}+\xi_{2})}T_{D}(c)$, and in (12) $K_{2}=k_{2}/k_{1}$. Using the parameter values of atri1993single (see Appendix, Table 1), we obtain $K_{2}=1$, $\Gamma=40/7=5.71$, and $K=1/7$. Also, taking values of $E$, $\nu$ and of the viscosity from zhou2009actomyosin ($E=8.5$ Pa, $\nu=0.4$ and $\xi_{1}+\xi_{2}=100$ Pa.s we find that $k_{\theta}$ is between 0.36 and 0.75. For simplicity, and since the parameter values for the calcium dynamics are anyway approximate, we fix $k_{\theta}=1$. Furthermore, $\displaystyle{T(c)=\frac{\tau_{h}}{(\xi_{1}+\xi_{2})}T_{D}(c)=\frac{\tau_{h}}{(\xi_{1}+\xi_{2})}T_{0D}\hat{T}(c)}$ where $\hat{T}(c)$ is nondimensional, and we also fix $\frac{\tau_{h}}{(\xi_{1}+\xi_{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.

4 Analysis of the model

4.1 The bifurcation diagrams of the Atri model (no mechanics)

The nondimensional Atri system is

[TABLE]

In Appendix A.2 we carry out a linear stability analysis of (14)–(15) and a bifurcation analysis with $\mu$ as the bifurcation parameter and we find that the parameter range for relaxation oscillations (limit cycles) is $0.289\leq\mu\leq 0.495$, as in atri1993single . In Appendix A.2 more details on the bifurcation structure of the system are given.

In Figure 3 we plot the bifurcation diagrams for the Atri system. In Figure 3(a) we present the amplitude of oscillations. The left Hopf point (LHP) and the right Hopf point (RHP) are, respectively, at $\mu=0.289$ and $\mu=0.495$. There are stable limit cycles and unstable limit cycles. The amplitude of oscillations increases with $\mu$ except for a small range of $\mu$ values near the RHP. In Figure 3(b) the frequency of the stable and of the unstable limit cycles are shown, respectively. The range of $\mu$ 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 $\mu=0.511$, where the stable and unstable limit cycle branches coalesce, is also visible. The frequency of the stable limit cycles increases slowly as $\mu$ increases and the lower, stable branch approximates the square root of $\mu$, as predicted by bifurcation theory kuznetsov2013 .

4.2 Linear stability analysis of the mechanochemical model

The steady states of the system (11)–(12) satisfy

[TABLE]

For any $\hat{T}(c)$, using (16), we can easily plot the steady states as a function of $\mu$ and $\lambda$. The Jacobian of (11)–(12) is given by

[TABLE]

and the characteristic polynomial is conveniently factorised as

[TABLE]

where $\omega$ represents the eigenvalues. As one eigenvalue is always equal to -1, the bifurcations of the system can be studied through the quadratic

[TABLE]

To identify the $\mu$-$\lambda$ parameter range sustaining oscillations we seek the Hopf bifurcations which satisfy Tr$(M_{2})=0$, Det$(M_{2})>0$, Discr$(M_{2})<0$, where

[TABLE]

[TABLE]

and substituting in (16) we obtain

[TABLE]

Hence, we can easily obtain the Hopf curve, for any $\hat{T}(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 $\mu$-$\lambda$ 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$

[TABLE]

Equations (23) and (16) constitute a linear system for $\mu$ and $\lambda$, so we again easily derive parametric expressions for $\mu(c)$ and $\lambda(c)$.

Similarly, to determine parametric expressions for the curve on which Discr$(M_{2})$ changes sign we set Discr$(M_{2})=0$

[TABLE]

which is quadratic in $\mu$ and linear in $\lambda$. Combining (24) with (16) we can again determine parametric expressions for $\mu$ and $\lambda$. In summary, we have developed a quick method for determining the three key curves characterising the geometry of the bifurcation diagram, for any $\hat{T}(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.

5 Illustrative examples

5.1 Contraction stress is a Hill function $\hat{T}(c)$ of order 1

5.1.1 Hopf curves

We assume that the calcium-induced stress $\hat{T}(c)$ is the Hill function

[TABLE]

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 stefan2008allosteric . Also, $\hat{T}=0$ when $c=0$, i.e. we assume no contraction stress without calcium. $\hat{T}^{\prime}(0)=\alpha$ is the rate of increase of $\hat{T}$ at $c=0$ and $1/\alpha$ is the scale of ‘ascent’ to the saturation level $T_{s}$. Therefore, we can call $\alpha$ the ‘mechanical responsiveness factor’ of the cytoskeleton to calcium.

Choosing $\displaystyle{\hat{T}(c)=10c/(1+10c)}$ as an illustrative example, in Figure 4 we use (16) to plot the steady state as a function of $\mu$, for selected increasing values of $\lambda$ (equilibrium curve). For $\lambda<4$ the equlibrium curve is qualitatively similar to that of the Atri model (see Figure 3(a)) but at $\lambda=4$ the curve changes qualitatively and a second non-zero steady state exists for $4<\lambda<40/7$, and a part of the curve corresponds to negative values of $\mu$ (see Appendix A.3 for details). For $\lambda>40/7$ no steady state exists for positive $\mu$ and hence $\lambda\leq 40/7$ is the biologically relevant range of the model, for $\alpha=10$ (see Appendix A.3).

In Figure 5 we plot the Hopf curve and the fold curve. We observe the following: (i) for $\lambda=0$ we recover the Hopf points and the fold points of the Atri model, as expected. (ii) As $\lambda$ increases the range of $\mu$ that sustains oscillations decreases. There is a global minimum value of $\mu$ that can sustain oscillations, $\mu_{\rm min}$. (iii) The oscillations are suppressed for a critical maximum value of $\lambda$, $\lambda_{\rm\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 $\lambda$ increases:

$\bullet$ for $\mu<\mu_{\rm min}$ no increase in $\lambda$ will ever elicit oscillations.

$\bullet$ for $0.203=\mu_{\rm min}<\mu<0.289$ when $\lambda$ reaches a certain value, $\lambda_{\rm OSC}$, oscillations are elicited, and $\lambda_{\rm OSC}$ decreases as $\mu$ approaches $0.289$. The oscillations vanish at a larger value of $\lambda$.

• •

for $IP_{3}$ values for which the Atri system sustains oscillations ($0.289<\mu<0.495$) in the mechanochemical model oscillations eventually vanish at a critical $\lambda$. This critical $\lambda$ decreases monotonically as $\mu$ increases towards $0.495$.

• •

for high $IP_{3}$ values ($\mu\geq 0.495$) no oscillations are sustained in the Atri system and a further increase in $\lambda$ will never elicit oscillations.

Therefore, for fixed cytoskeletal mechanical responsiveness factor, $\alpha=10$, and for fixed parameter values as in atri1993single a range of behaviours emerge as $\mu$ and $\lambda$ 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 the effect of 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 christodoulou2015cell (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 christodoulou2015cell . At low levels of $IP_{3}$ and hence low levels of calcium, cells are not able to contract and therefore apical constriction 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 $\lambda$ increases this range decreases and eventually vanishes at $\lambda\approx 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 of $\hat{T}(c)$ we may choose, and thus examine quickly the effect of mechanics on calcium oscillations. We note that in the experiments of christodoulou2015cell the calcium-induced stress saturates to a non-zero level as calcium levels increase and hence we chose a $\hat{T}(c)$ that saturates. In other cell types it is possible that the cell can relax back to baseline stress and in such a case $\hat{T}(c)$ would not be described by a Hill function, and more experiments should be undertaken to determine the appropriate form of $\hat{T}(c)$.

5.1.2 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 $\hat{T}(c)=10c/(1+10c)$.

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

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

In Figure 7 we plot the oscillation amplitude as a function of $\mu$, for three selected values of $\lambda$, using XPPAUT. We see that as $\lambda$ increases the amplitude decreases until the oscillations vanish close to $\lambda=\lambda_{\rm max}=1.69$, which agrees with the Hopf curve in Figure 5. We also observe that for $\lambda=0.5$ and $1$, in Figures 7(a) and 7(b) respectively, there are both stable and unstable limit cycles, and the right Hopf point is subcritical. Also, as $\lambda$ increases, the $\mu$-range of unstable limit cycles decreases until it vanishes; for $\lambda=1.5$ (Figure 7(c)) 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 $\mu$ in the mechanochemical system.

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

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

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

5.1.3 Varying the cytosolic mechanical responsiveness factor

We now investigate if the Hopf curve changes qualitatively as the cytosol’s mechanical responsiveness factor, $\alpha$, varies. In Figure 10(a), using the parametric expressions (21)–(22) we plot the Hopf curves for increasing values of $\alpha=1,2,10,100$. We observe that the Hopf curve changes qualitatively; for $\alpha\approx 2$ it develops a cusp and for smaller values of $\alpha$ there is a “bow-tie”. This geometrical change corresponds to yet another bifurcation, with $\alpha$ as a bifurcation parameter111The cusped Hopf curve is the bifurcation curve of a cusp catastrophe surface, according to the catastrophe theory developed by Zeeman zeeman1977 , and subsequently by Stewart and collaborators poston2014 ; stewart2014 .. However, as for $\alpha=10$, oscillations always vanish for a sufficiently large value of $\lambda$, $\lambda_{\rm max}$.

We also observe that as $\alpha$ increases, $\lambda_{\rm max}$, the critical stretch activation value beyond which oscillations vanishes, decreases, i.e. oscillations are sustained for a smaller range of $\lambda$ values. To investigate this more systematically we have determined parametric expressions for $\lambda_{\rm\max}$ and $\alpha$ as functions of $c$, and we plot $\lambda_{\rm\max}(\alpha)$ in Figure 10(b). We see that as $\alpha$ increases, $\lambda_{\rm max}$ decreases monotonically, and hence oscillations are sustained for an increasingly smaller range of $\lambda$, which agrees with Figure 10(a). Also, since $\lambda_{\rm max}(\alpha)$ asymptotes to a positive value as $\alpha\to\infty$, for any $\hat{T}(c)=\frac{\alpha c}{1+\alpha c}$, the system will always sustain some oscillations, irrespectively of the value of $\alpha$. Therefore, we predict that for cytosols that are more responsive to calcium (higher $\alpha$), oscillations vanish at a lower $\lambda_{\rm\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 on $\hat{T}(c)$, $\mu_{\rm min}$ is constant and not zero for any $\alpha$. Therefore, as we expect, $IP_{3}$ is always required in order to obtain oscillations, for any $\lambda$ and any $\alpha$ but the minimum level of $IP_{3}$ does not change with $\alpha$. Also, for fixed $\lambda$, as $\alpha$, the mechanical responsiveness factor of the cytosol, increases, the $IP_{3}$ level required to induce oscillations also decreases. Additionally, for fixed $\mu$, as $\alpha$ increases $\lambda_{\rm 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 $\mu$ and $\lambda$ when the contraction stress is modelled as a Hill function of order 1.

5.2 Hopf curves for $\hat{T}(c)$ a Hill function of order 2

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

Comparing Figure 11 with Figure 10(a) we see that the Hopf curves have the same qualitative behaviour for the two Hill functions. Oscillations can be sustained for any value of $\alpha$ and they always vanish for a sufficiently large value of $\lambda$, Also, as in the Hill function of order 1 as $\alpha$ increases $\lambda_{\rm max}$ decreases while $\mu_{\rm min}$ is constant. Also, a cusp again develops but for the Hill function of order 2 the value of $\alpha$ where 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.

6 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 Apical Constriction (AC), an essential morphogenetic process which, if disrupted, leads to embryo abnormalities christodoulou2015cell .

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 Apical Constriction.

We have then developed a simple mechanochemical ODE model that consists of an ODE for $\theta$, 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 atri1993single . 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 Oster, Murray and collaborators murray1984generation ; oster1984mechanochemistry ; murray1988 ; murray2001 . 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 $\hat{T}(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, $\lambda\theta$, 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 $\lambda$ as a combination of the structural characteristics of a SSC; $\lambda$ 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 $\hat{T}(c)$, we have analytically identified the parameter regime in the $\mu$–$\lambda$ plane corresponding to calcium oscillations and applied this result in two illustrative examples, $\hat{T}(c)=\alpha c/(1+\alpha c)$ and $\hat{T}(c)=\alpha c^{2}/(1+\alpha c^{2})$. In both cases, as $\lambda$ increases, the oscillations are eventually suppressed at a critical $\lambda$, $\lambda_{\rm\max}$—see, respectively, Figures 10(a) and 11. The prediction is in agreement with experiments christodoulou2015cell 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 Apical Constriction 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 $\hat{T}(c)=\alpha c/(1+\alpha c)$ we have also shown analytically that as $\alpha$, the mechanical responsiveness factor of the cytosol increases, $\lambda_{\rm\max}$ decreases but it never becomes zero (see Figure 10(b)). This means that for any $\alpha$, there will always be a $\mu$-$\lambda$ region for which oscillations are sustained. Furthermore, for the illustrative example of $\hat{T}(c)=10c/(1+10c)$ we have determined numerically the oscillation amplitude and frequency as the bifurcation parameters $\mu$ and $\lambda$ vary, using XPPAUT. We found that the behaviour is qualitatively similar to the Atri model (see Figure 3) for lower $\lambda$ values but that it changes for larger $\lambda$ values (see Figure 6). We found that, as $\lambda$ increases the amplitude of oscillations decreases (see Figure 7) but their frequency decreases (see Figure 8). More experiments could be undertaken to test these predictions.

In the experiments of christodoulou2015cell 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 $\hat{T}(c)$ cannot be modelled as a Hill function. Experiments could be undertaken also in other calcium-induced mechanical processes to determine the appropriate form of $\hat{T}(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 brodland2006cell ; luby1999cytoarchitecture ; zhou2009actomyosin . 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 thul2014 ; cao2014 . A multitude of deterministic and stochastic calcium models have been developed atri1993single ; wilkins1998 ; goldberg2010 ; gracheva2001stochastic ; sneyd1994model ; sneyd1998calcium ; timofeeva2003 ; see also the comprehensive reviews sneyd2003 ; thul2014 ; rudiger2014stochastic and the books keener1998 ; dupont2016 , 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 in 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 to models in more realistic geometries. Moreover, we have fixed all parameters here, except $\mu$, $\lambda$ and $\alpha$; 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 them 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.

Acknowledgements.

We thank James Sneyd, Vivien Kirk, Ruediger Thul, Lance Davidson, Abdul Barakat, Thomas Woolley, Iina Korkka and Bard Ermentrout for valuable discussions. Katerina Kaouri also acknowledges support from two STSM grants awarded by the COST Action TD1409 (Mathematics for Industry Network, MI-NET) for research visits to Oxford University.

Appendix A Appendix

A.1 Supplementary information on the presented experimental results

Figure 1: Plot of surface area (measured in $\mu m^{2}$) and calcium oscillation amplitude over time of a single embryonic, epithelial cell undergoing Apical Constriction in Xenopus [13]. Measurements were taken every 10 seconds from a time lapse movie of a stage 9 Xenopus embryo expressing Lulu-GFP + GECO-RED (see [13]). To normalise the surface area, all the area values were divided by the first measured value (which is the largest since the surface area is decreasing with time). For measuring calcium oscillation amplitudes, the signal intensity of the non-ratiometric calcium sensor (GECO-RED) was also measured over time. For normalization, all values were divided by the highest signal intensity value. Note that we tracked the surface area and calcium level of cells induced to undergo AC at gastrula stage in order to decouple AC from other mophogenetic 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.

Figure 2: (a) Plot of surface area ($\mu 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 minutes. 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 Analysis of the Atri model

A.2.1 Linear stability analysis

The steady states (S.S.) of (14)-(15) are the intersections of the nullclines of the system. Setting

[TABLE]

we obtain

[TABLE]

which can be cast as a quartic in $c$. (Note that (16) reduces to (28) for $\lambda=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 $\mu$ 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=\frac{\Gamma}{\mu K_{1}}$ as $c\to\infty$. $h_{M}$ can be found analytically by differentiating (26):

[TABLE]

and setting $dh/dc=0$, which leads to a quadratic equation for $c$. Rearranging, and since $c>0$, we discard the negative root, obtaining

[TABLE]

Therefore, $h_{M}$ scales with $\Gamma/(\mu K_{1})$ and depends on the parameters $K$ and $b$ in a much more complicated manner. For the parameter values from [3] we have $c_{M}=0.169$ and $h_{M}=0.279/\mu$.

Nullcline (27) is a decreasing function of $c$; it has a maximum at (0,1) and tends to [math] as $c\to\infty$. For $\mu_{1}=0.28814$ and $\mu_{2}=0.28925$ the nullclines touch and there is a double steady state; for values of $\mu<\mu_{1}$ and $\mu>\mu_{2}$ there is one S.S. and there are three S.S. for $\mu_{1}<\mu<\mu_{2}$. ($\mu_{1}$ and $\mu_{2}$ are obtained by differentiating (28) and finding the roots of $\frac{d\mu}{dc}=0$.) Note that we present the values of $\mu$ with five decimal places because the bifurcation analysis depends sensitively on $\mu$, 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

[TABLE]

Taking the parameter values of [3] we identify the bifurcations of the system as $\mu$ increases by determining where the Tr, Det, and Discr change sign. We find a richer bifucation structure as $\mu$ increases. This behaviour was not analysed in such detail in [3] or in later literature. Given the very sensitive dependence on precise values of $\mu$, these details are probably of more mathematical interest than of biological significance. The parameter values are summarised in Table 1.

• $\bullet$

$0<\mu<0.27828$: one stable node.

• $\bullet$

$\mu=0.27828$: the stable node becomes a stable spiral (bifurcation Discr=0)

• $\bullet$

$\mu=0.28814$: Stable spiral present. Also, a saddle and an unstable node (UN) emerge (bifurcation Det=0, fold point)

• $\bullet$

$\mu=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)

• $\bullet$

$\mu=0.28924$ the unstable spiral becomes an unstable node, and we have two unstable nodes and a saddle (Discr=0)

• $\bullet$

$\mu=0.28925$: one unstable node (Det=0, fold point)

• $\bullet$

$\mu=0.28950$: the unstable node becomes an unstable spiral (Discr=0)

• $\bullet$

$\mu=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\leq\mu\leq 0.49500$ since at $\mu=0.28900$ a Hopf bifurcation (HB) arises, the stable spiral becomes unstable, and we expect relaxation oscillations (limit cycles) in the nonlinear system. As $\mu$ increases the unstable spiral becomes a stable spiral close to $\mu=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 $\mu=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

[TABLE]

Equation (34) is decoupled from equations (32) and (33), and we can thus continue with a two-dimensional analysis for equations (32)–(33). The steady states satisfy $\displaystyle{\lambda\hat{T}(c)=\frac{\Gamma c}{K+c}}$. The Jacobian of the system (32)–(33) has entries

[TABLE]

Hence

[TABLE]

Therefore, for any $\hat{T}(c)$, ${\rm Discr}(M_{1})>0$, ${\rm Tr}(M_{1})<0$ always, and ${\rm 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 of $\hat{T}$ a non-zero steady state may exist and this means biologically that even without calcium flux from the ER a non-zero calcium concentration can be sustained in the cytosol due to the stress-induced calcium release.

For $\hat{T}(c)$ as given in (25), there is always a steady state ($c^{\star},\theta^{\star}$)=($0,0$). A second steady state

[TABLE]

exists if

[TABLE]

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 $\lambda$ values for which the system sustains a non-zero calcium concentration even without a CICR flux. For even larger values of $\lambda$ there is no steady state and the model ceases to be biologically relevant.