Activation of Store-Operated Calcium Entry in Airway Smooth Muscle Cells: Insight from a Mathematical Model

Intracellular dynamics of airway smooth muscle cells (ASMC) mediate ASMC contraction and proliferation, and thus play a key role in airway hyper-responsiveness (AHR) and remodelling in asthma. We evaluate the importance of store-operated entry (SOCE) in these dynamics by constructing a mathematical model of ASMC signaling based on experimental data from lung slices. The model confirms that SOCE is elicited upon sufficient depletion of the sarcoplasmic reticulum (SR), while receptor-operated entry (ROCE) is inhibited in such conditions. It also shows that SOCE can sustain agonist-induced oscillations in the absence of other influx. SOCE up-regulation may thus contribute to AHR by increasing the oscillation frequency that in turn regulates ASMC contraction. The model also provides an explanation for the failure of the SERCA pump blocker CPA to clamp the cytosolic of ASMC in lung slices, by showing that CPA is unable to maintain the SR empty of . This prediction is confirmed by experimental data from mouse lung slices, and strongly suggests that CPA only partially inhibits SERCA in ASMC.


Introduction
Ca 2z is a ubiquitous cellular messenger, controlling a wide range of biological functions. These include ASMC contraction and proliferation, which are associated with airway hyperresponsiveness (enhanced contractility) and airway remodelling (structural changes) in asthma. The main trigger for cytoplasmic Ca 2z (½Ca 2z i ) increase in ASMC is agonist stimulation at the cell membrane (e.g., by histamine released from mast cells or acethylcholine released from nerves). Binding of agonist to Gprotein coupled receptors induces the production of IP 3 , a second messenger which diffuses into the cytosol and binds to IP 3 receptor Ca 2z channels (IPR) on the sarcoplasmic reticulum (SR) membrane (Fig. 1). This causes the IPR to open and release Ca 2z from the SR into the cytosol (the SR being the main Ca 2z store in ASMC). As ½Ca 2z i exerts a positive feedback on IPR, this results in so-called Ca 2z -induced Ca 2z release (CICR). The release is terminated by the inhibition of the IPR at large ½Ca 2z i , and Ca 2z is pumped back into the SR by Ca 2z ATP-ases (SERCA). Hence, for sufficient IP 3 concentration, cycling of Ca 2z through IPR can occur, and give rise to the repetitive propagation of ½Ca 2z i waves through the cytosol. These appear as ½Ca 2z i oscillations at the whole-cell level. Importantly, airway contraction increases with the frequency of these ½Ca 2z i oscillations [1,2]. Ca 2z dynamics are also involved in ASMC proliferation [3][4][5], and in the assembly of myosin thick filament and actin thin filament [6][7][8], which form the contractile machinery of ASMC. In addition, several Ca 2z channels and pumps in ASMC are regulated by inflammatory mediators present in asthma (e.g., [4,[9][10][11][12]). Ca 2z dynamics therefore appear to be involved in multiple interrelated aspects of asthma at the cellular level. In the present work, we use mathematical modelling to investigate the important Ca 2z pathways at play in Ca 2z dynamics of ASMC and thus improve our understanding of airway hyper-responsiveness and remodelling in asthma.
Store-operated Ca 2z entry (SOCE) is one important Ca 2z entry mechanism, in which plasma membrane (PM) Ca 2z channels open in response to Ca 2z store depletion. These are called store-operated Ca 2z channels (SOCC). Although the concept of SOCE was proposed 25 years ago [13], the mechanism of its activation has been identified only recently [14]. The process is mediated by stromal interaction molecules (STIM), proteins embedded in the SR membrane which are sensitive to SR Ca 2z . Upon dissociation of Ca 2z from their SR binding site, they oligomerise and translocate within the SR membrane to the plasma membrane. Here, STIM proteins bind to Orai and/or TRP, the proteins forming the pore of SOCC, and trigger their opening (Fig. 1). Although SOCE has been identified in many cells, it is generally stimulated by artificial emptying of the Ca 2z store, as there is unfortunately no specific pharmacological SOCC blocker. Hence, the importance of store depletion, and therefore of SOCE, during physiological conditions such as Ca 2z oscillations, remains largely unknown. This may explain why SOCE has been included only in a few mathematical models of Ca 2z dynamics [15][16][17][18]. In particular, no prior modelling work on Ca 2z dynamics in ASMC [19][20][21][22][23] has taken SOCE into account, even though there is evidence that SOCE is up-regulated by inflammatory mediators found in asthma (TNF-a and IL-13) [9,11,24], and is associated with ASMC proliferation [3,5].
In this paper, we develop a mathematical model to evaluate the importance of SOCE in Ca 2z dynamics of ASMC. While there is much evidence that SOCE occurs upon SR depletion in cultured ASMC in vitro (e.g., [3,[25][26][27]), these cultured cells often lose their contractile phenotype, and rarely display agonist-induced Ca 2z oscillations. Hence, ASMC in lung slices, which retain most of their physiological and morphological characteristics, are a more reliable preparation to study ASMC ½Ca 2z i dynamics. Moreover, the available data from lung slices reflect ½Ca 2z i dynamics in individual ASMC, while the majority of works with cultured cells provide only global imaging of ½Ca 2z i over wells containing thousands of ASMC. Therefore, we base our model on data from lung slices. SOCE has not been studied directly in lung slices, but a treatment with ryanodine-caffeine (Rya-Caf) has previously been used to clamp the cytosolic Ca 2z of ASMC [2,28,29], which relies on emptying the SR from Ca 2z . The results of these experiments therefore provide invaluable information about SOCE. Because agonist stimulation was systematically performed before Rya-Caf treatment to ensure that the lung slice is viable, i.e., that ASMC exhibit normal ½Ca 2z i oscillations and contraction, we can construct a mathematical model of ½Ca 2z i dynamics informed by these data that accounts for both physiological and non-physiological conditions. The model is then used to i) evaluate the effect of SOCE up-and down-regulation on agonist-induced ½Ca 2z i oscillations, and (ii) explain the inability of the SERCA pump blocker CPA to clamp the ½Ca 2z i , in contrast with Rya-Caf treatment.

Ethics Statement
The experimental study followed the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. The protocol was approved by the Institutional Animal Care Committee of the University of Massachusetts Medical School (Docket Number: A-836-12). Animals were euthanized with sodium pentobarbital before tissue collection.

Experimental data
Data consist of fluorescence recordings of ½Ca 2z i dynamics in ASMC within intact lung slices. All the materials and methods have been previously described (e.g., [2,28]). Essentially, ½Ca 2z i imaging was performed from regions of about 4mm 2 within ASMC (Fig. 2), using two-photon laser scanning microscopy. The fluorescent indicator employed was Oregon Green BAPTA-1-AM, which has a high affinity for Ca 2z (K d^0 :2mM). We use published data [2] to develop the mathematical model, and new experimental results to test the model predictions (see Results). The latter data can be made freely available upon request for academic, non-commercial use.

Mathematical model
Intracellular Ca 2z dynamics are modelled at the whole-cell level, via the following system of ordinary differential equations (e.g., [30]): where c~½Ca 2z i is the free cytosolic Ca 2z concentration, and c s~½ Ca 2z SR is the free SR Ca 2z concentration. The term J in represents the total influx of Ca 2z into the cytosol through PM channels; J PMCA , the Ca 2z efflux through the PM Ca 2z ATP-ase pumps (PMCA); J rel , the Ca 2z flux of Ca 2z from the SR into the cytosol, and J SERCA , the flux of Ca 2z from the cytosol into the SR through the SR/ER Ca 2z ATP-ases  (SERCA). The factor c represents the ratio of cytoplasmic volume to SR volume, and implicitly incorporates the relative effect of fast, linear (e.g., low affinity) Ca 2z buffers in the SR compared to the effect of similar buffers in the cytosol. Indeed, the effect of fast, linear buffers amounts to a global rescaling of the Ca 2z fluxes in the corresponding compartment (e.g., [30]). The other buffers are assumed to have a negligible effect on Ca 2z dynamics at the whole-cell level (see also Discussion).
We assume that where J leakin is a constant Ca 2z leak through unspecified channels, J ROCC is the Ca 2z influx through receptor-operated Ca 2z channels (ROCC) and J SOCC the influx through SOCC. We neglect the Ca 2z influx through voltage-operated Ca 2z channels (VOCC) because membrane depolarisation plays little role during agonist-induced ½Ca 2z i signalling and contraction in ASMC (in contrast to other types of muscle cells, including vascular smooth muscle cells, where action potentials are crucial to contraction) [1,31]. The Ca 2z influxes are modelled by: where a 0 and a 1 are constants, p is the agonist concentration, V s is the maximum SOCC flux, and P so represents the fraction of STIM proteins bound to Orai/TRP proteins, i.e. the fraction of activated SOCC. This fraction adapts slowly to changes in c s , because the diffusion of STIM within the SR membrane is a slow process [32]. We model this phenomenologically by The steady-state function P ? so can be interpreted as the fraction of STIM proteins dissociated from SR Ca 2z (as a consequence of store depletion), and thus able to oligomerise and move toward the PM to bind with Orai and/or TRP (see also Discussion). P ? so is therefore a decreasing function of ½Ca 2z SR , which we model by the reverse Hill function Eq. (0b), assuming affinity K s for ½Ca 2z SR and Hill coefficient n s~4 [33].
The total Ca 2z flux from the SR into the cytosol is given by where J IPR is the Ca 2z flux through IP 3 receptors (IPR), J RyR the Ca 2z flux through ryanodine receptors (RyR), and J leakSR an unspecified Ca 2z leak out of the SR. We use the formulation (e.g., [30]): where k IPR (resp. k RyR ) is the maximum rate of Ca 2z flow through IPR (resp. RyR). Following [23], the IPR opening probability P IPR is modelled using the Li-Rinzel/Tang et al. reduction of the De Young-Keizer (DYK) model [34][35][36]: where p is the IP 3 concentration, and y is the fraction of inhibited IPR. The latter obeys with The parameters K i~k . . ,4) are equilibrium constants for IP 3 and Ca 2z binding/unbinding to the IPR; we use the original values from the DYK model [34]. The value of k { 2 is scaled so that the range of ½Ca 2z i oscillation frequencies matches the experimental range, with the ratio k 42~k { 4 =k { 2 kept constant to the value in ref. [30] (see also Table 1). In the experiments modelled in this work, either RyR play a negligible role, or they are locked open by Rya-Caf treatment (see Results). Hence, we neglect their dynamics and set the fraction of open RyR, P RyR , either to 0 or 1 depending on the experiment considered.
The Ca 2z ATP-ases are modelled using the usual expressions (e.g., [30]): We do not model Ca 2z pumping into mitochondria explicitly, but acknowledge that a portion of the extrusion process attributed to PMCA might actually be performed by mitochondria uniporters, as these might be activated at average ½Ca 2z i as low as 1mM [20].
Gathering all expressions, the model is described by: In addition to Eq. (11), we use the following expressions to account for the time needed by drugs to reach full effect: These equations describe respectively agonist stimulation, Rya-Caf treatment, and SERCA block by CPA (see Results).
Unless otherwise mentioned, parameter values were freely adapted (within physiological ranges when they are known) to account for the experimental results. The values retained are listed in Table 1. The fitting was performed ''by hand'' (i.e., no algorithmic method was used) within the Mathematica ''Manipulate'' environment (a useful framework for fitting an ODE model to several experimental results as it enables visualisation of the effect of a parameter change on several ODE integrations almost instantaneously). The code can be made freely available upon request for academic, non-commercial use.
All simulations were run from the same initial condition as in the experiment, which is usually the physiological equilibrium.
Bifurcation diagrams were computed using the numerical continuation software AUTO [37,38].

Results
Accounting for ½Ca 2z i dynamics of AMSC in lung slices Fig. 3A-C shows representative ½Ca 2z i dynamics of an ASMC in a human lung slice in response to a three-step experimental protocol [2]. This protocol was originally designed to clamp the ½Ca 2z i of ASMC, in order to study independently the effects of agonist and ½Ca 2z i on airway contraction [28]. The slice is first stimulated with agonist (histamine), to verify its viability (Fig. 3A). This induces ½Ca 2z i oscillations. Agonist is then washed from the slice, and a Rya-Caf treatment is applied (Fig. 3B). This creates a permanent Ca 2z leak through RyR, because caffeine opens RyR and ryanodine locks them open irreversibly. If this Ca 2z leak is large enough, it keeps the SR empty and prevents any further change in ½Ca 2z i , unless extracellular Ca 2z is modified. The effectiveness of the treatment is confirmed by the second application of agonist (Fig. 3C): no further ½Ca 2z i increase is triggered, showing that ½Ca 2z i is clamped. It is important to emphasise that these results are not specific to histamine stimulation of human lung slices: similar results have been obtained in mouse and rat lung slices with methylcholine ( Fig. 6 in ref. [29], Figs. 5B and 6C-D in ref. [28]). For p p~0 and k RyR~0 , the equilibrium Ca 2z concentrations are c Ã~6 8:1nM and c Ã s~1 58mM, which are in the physiological ranges [40,41]. doi:10.1371/journal.pone.0069598.t001 The mathematical model enables the deduction of valuable information from the experimental results. First, from Eq. (14), the new, elevated, ½Ca 2z i equilbrium reached after Rya-Caf treatment satisfies: where c Ã and c Ã s are respectively the equilibrium ½Ca 2z i and ½Ca 2z SR . An important consequence of (18) is that, in the absence of SOCE, c Ã depends only on the Ca 2z fluxes through the PM. This may seem surprising, as any increase in Ca 2z flux out of the SR (J rel in Eq. (1)) is expected to increase ½Ca 2z i . However, the equilibrium equation (18) tells us that such an increase would only be transient (because the PMCA pumping rate is an increasing function of ½Ca 2z i ), unless there is a concomitant permanent increase in Ca 2z influx through the PM. Hence, the persistence of an elevated ½Ca 2z i means that a permanent SOCE has been elicited (as SOCE is the only Ca 2z influx capable of increase upon Rya-Caf treatment). Moreover, the model indicates that ROCE is negligible after Rya-Caf treatment. Indeed, if it was not, the addition of agonist would increase c Ã via the increase in J ROCC . Hence, we assume that the ROCE rate a 1 is small (see Table 1 and Discussion).
Results of ''hand-fitting'' the model to the experimental results are shown in Figs. 3D-F and Fig. 4, with the corresponding parameter values listed in Table 1. The model reproduces (i) the agonist-induced Ca 2z oscillations, (ii) the similar magnitudes of the new equilibrium ½Ca 2z i in Fig. 3B and the amplitude of the oscillations in Fig. 3A, and (iii) the negligible effect of agonist stimulation after Rya-Caf treatment. Agonist-induced Ca 2z oscillations were simulated with P RyR~0 because RyR appear to play a negligible role during agonist-induced Ca 2z oscillations [2,39]. On the other hand, the response to Rya-Caf was simulated with P RyR~1 since the treatment locks open the RyR. We did not attempt to reproduce the magnitude of the initial spike response to Rya-Caf treatment relative to that of the subsequent ½Ca 2z i plateau (Fig. 3B) because the fluorescent dye used in the experiments saturates rapidly with ½Ca 2z i . Parameter values were also adjusted to yield physiological Ca 2z equilibrium concentrations (c*0:1mM [40] and c s *500mM [41]), realistic ½Ca 2z i oscillation amplitude ( * ; 1mM), and to reproduce the range of ½Ca 2z i oscillation frequencies observed in human lung slices as a function of agonist (0.5-11/min [2]). More detail on the parameter estimation procedure is given in Supporting Information S1. Fig. 4 shows the bifurcation diagram of the model as a function of agonist concentration. Periodic solutions (i.e., Ca 2z oscillations) arise through a Hopf bifurcation, and disappear through a saddle-node bifurcation of limit cycles. A second Hopf bifurcation is present on the steady-state branch, and is associated with a region of bistability between the steady-state and the periodic solution at the right of the bifurcation diagram. It is not known whether such bistability occurs in reality. It should also be noted that the steady-state ½Ca 2z i increases with agonist concentration, as is expected (e.g., [30]). This increase is provided by SOCE in our model. Indeed, the Ca 2z flux through IPR increases with agonist, so that store depletion increases as well.

Effect of SOCE regulation on agonist-induced Ca 2z
oscillations SOCE is the main Ca 2z influx in the model, as ROCE is negligible (see above) and the Ca 2z leak influx is (by definition) small. Fig. 3D shows that while SOCE is almost zero at physiological equilibrium (initial condition), it substantially increases during agonist-induced Ca 2z oscillations (final condition; see also Fig. 4), due to significant SR Ca 2z depletion. Therefore, changes in SOCE can be expected to have a substantial effect on Ca 2z oscillations. This is quantified in Fig. 5, where the amplitude and frequency of Ca 2z oscillations are plotted as a function of (a) the maximum SOCE rate, V s , and (b) STIM affinity for Ca 2z , K s (the ½Ca 2z SR at which half SOCC are open). It is found that the ½Ca 2z i oscillation frequency varies as much with V s and K s at fixed agonist concentration (Fig. 5) as it varies with agonist concentration at fixed SOCE parameters (Fig. 4). Moreover, a too big departure from the ''normal'' values (dotted lines, Table 1) leads to the extinction of the Ca 2z oscillations (via a Hopf bifurcation to the left, and a saddle-node to the right, of the bifurcation diagrams in Figs. 5A-B). These results are not very surprising to the extent that Ca 2z oscillations are expected to depend crucially on Ca 2z influx (e.g., [42]). However, they suggest that SOCE could play a role in AHR since (i) ½Ca 2z i oscillations mediate ASMC contraction, and (ii) SOCE upregulation (which increases ½Ca 2z i oscillation frequency) can be triggered by inflammatory mediators commonly found in asthma [9,11,24].

Partial inhibition of SERCA by CPA
We now apply the model to experimental data from mouse lung slices showing an attempt to clamp ½Ca 2z i with the SERCA blocker CPA, instead of Rya-Caf treatment (Fig. 6). After inducing Ca 2z oscillations with agonist, CPA is applied in the presence of agonist (for faster emptying of the SR than CPA alone) and causes a gradual damping of the Ca 2z oscillations, together with a rise of the ½Ca 2z i baseline, until the oscillations become undistinguishable from fluctuations around an elevated steady ½Ca 2z i mean. Because CPA is believed to inhibit SERCA, the assumption, at this stage of the experiment, is that the SR is empty and SOCE fully shows the frequency of the Ca 2z oscillations on the main stable segment (from the upper blue dot to the black cross), which fits the experimental range in human [2]. The stable solutions are represented as thick lines and unstable solutions as thin lines. The green diamonds represent Hopf bifurcations, the black cross, a saddle-node bifurcation, and the blue dots, period-doubling points. Period-doubled branches are not shown because they extend only over a tiny range of p values; moreover it is likely that the deterministic description of Ca 2z oscillations fails at these low agonist concentrations (see Discussion). The vertical dotted line indicates the value of p used in Fig. 3 (Table 1). doi:10.1371/journal.pone.0069598.g004 active. However, when agonist is removed (CPA remains), ½Ca 2z i falls. When agonist is reapplied, ½Ca 2z i increases. These ½Ca 2z i responses to agonist addition and removal are not observed when SOCE is evoked by Rya-Caf treatment. According to our model (Eq.(18)), the decrease in ½Ca 2z i upon agonist removal indicates that SOCE does not remain activated, i.e. that the SR refills with Ca 2z . This suggests that the SERCA are not completely blocked by CPA, as illustrated by the simulations in Fig. 6B-D. If CPA was to fully block the SERCA (Fig. 6B), ½Ca 2z i would not decrease upon agonist removal. If ½Ca 2z i falls, it must be because either CPA requires a longer time than that used in the experiment to fully block the SERCA (Fig. 6C), or CPA achieves only partial block of the SERCA (Fig. 6D).
Experiments of longer duration were performed to test the model predictions. Fig. 7A shows that if CPA is applied in the presence of agonist for 5 minutes, followed by CPA only for a further 10 minutes, ½Ca 2z i still returns to the original equilibrium level when agonist is removed, and remains low until agonist is reintroduced. This suggests that the explanation in Fig. 6C can be rejected, otherwise the longer exposure to CPA should yield a result similar to Fig. 6B. The inability of CPA to fully empty the SR of Ca 2z is confirmed by Fig. 7B, where extracellular calcium is removed before agonist is applied a second time, to prevent any potential ROCE. The Ca 2z response induced can thus be unambiguously attributed to Ca 2z release from the SR.
Hence, our combined modelling and experimental study indicates that CPA blocks only partially the SERCA of ASMC in lung slices (scenario simulated in Fig. 6D). This is a potentially important result given the wide use of CPA in cell biology to study SOCE. We note that Figs. 6A and 7 could also be explained by a model assuming that ROCE, instead of SOCE, is the main Ca 2z influx (e.g., [23]). However, such a model would fail to explain the outcome of Rya-Caf treatment in human and mouse lung slices (both the persistent elevated ½Ca 2z i in the absence of agonist, and the absence of effect of agonist on this elevated ½Ca 2z i ). In contrast, our model, constructed to account for both agonistinduced oscillations and Rya-Caf treatment, explains the CPA results without requiring any modification. Its prediction holds provided CPA is not a 100% efficient SERCA blocker, and this hypothesis is supported by the experimental data in Fig. 7.

Modelling SOCE
Our mathematical model accounts for the two main properties of SOCE: 1) SOCE is an increasing function of Ca 2z store depletion, and 2) it activates slowly upon store depletion. While the mechanisms of SOCE activation are rather well understood [14,32], the mechanisms of SOCE termination remain less clear [43,44]. Hence, we do not explicitly distinguish between SOCE activation and inactivation in the model, and use a single parameter K s for STIM affinity for SR Ca 2z and a single time constant t s for the slow adaptation to changes in ½Ca 2z SR . This is also justified by the fact that most experimental data available on SOCE come from a category of SOCC called CRACC (Ca 2zrelease-activated Ca 2z channels), which are highly selective to Ca 2z , while there is evidence that SOCE in ASMC (and in other cells) occurs at least in part through non-selective Ca 2z channels (NSCC). It could be that the latter operate somewhat differently from CRACC in response to store depletion or refilling.
Our description of SOCE slow activation upon store depletion is continuous, which is easy to handle computationally, and compatible with experimental knowledge. Indeed, it is reasonable to assume that a small fraction of STIM proteins reside in close proximity to the PM, and may thus bind Orai quickly upon store depletion. Hence, a weak SOCE is likely to occur almost instantaneously upon store depletion, rendering unnecessary to introduce a finite activation delay in the model via a delaydifferential equation.
We are aware of only few prior works on Ca 2z dynamics that include a mathematical description of SOCE, all of which are ODE models [15][16][17][18]. The first two were published before the molecular basis for SOCE was established. The latter two works include more realistic descriptions of SOCE, but none of them accounts for the slow translocation of oligomerised STIM to the PM, while it is recognised as the rate-limiting event for SOCE activation [32]. Ong et al. however assume a slow diffusion of Ca 2z between internal SR and superficial SR (modelled as distinct compartments exchanging Ca 2z ), with SOCE being triggered by peripheral SR depletion [17]. Liu et al. explicitly model both SR Ca 2z dissociation from STIM and binding of STIM to Orai. Both models are used to study transient ½Ca 2z i responses only; ½Ca 2z i oscillations are not considered. Prior models of Ca 2z dynamics specific to ASMC did not include SOCE, while we have shown that this is necessary to account for several experimental results obtained with lung slices. The work of Haberichter et al. [19] focused on the influence of the different IPR isoforms on Ca 2z signalling in ASMC. Brumen et al. studied the influence of the total Ca 2z content on the nature (damped or sustained) and frequency of agonist-induced Ca 2z oscillations [21]. Roux et al. did not model Ca 2z oscillations, but transient Ca 2z responses to caffeine [20]. Finally, the model by Wang et al. [23] addressed the different contributions of IPR and RyR to agonist-induced and KClinduced Ca 2z oscillations in ASMC.
From the mathematical point of view, the fact that SOCE is an explicit function of store Ca 2z renders the models of Ca 2z dynamics including this influx qualitatively different from those which do not, as SOCE couples the homogenous steady-state ½Ca 2z i to ½Ca 2z SR (Eq. (18)). This property is essential for the predictions of our model (in particular, the persistence of an elevated ½Ca 2z i upon sustained store depletion in the absence of agonist). On the other hand, whether SOCE is an instantaneous or delayed function of ½Ca 2z SR appears to have little effect on our results.

SOCE vs. ROCE
While Fig. 3C (as well as Fig. 6 in ref. [29], Figs. 5B and 6C-D in ref. [28]) shows that no ROCE is elicited by agonist following Rya-Caf treatment, it does not imply that ROCE cannot play a substantial role during other, more physiological, conditions, such as agonist-induced Ca 2z oscillations. It could be that ROCE is inhibited at the large ½Ca 2z i levels induced by SOCE activation following Rya-Caf treatment. Instead of assuming the existence of an inactivation process at large ½Ca 2z i , we assumed, for simplicity, that ROCE is negligible in the model. This approach enabled us to show that Ca 2z influx through SOCC is sufficient to sustain agonist-induced ½Ca 2z i oscillations, and to explain the experimental results obtained with CPA, although the latter could be interpreted as evidence for ROCE at first sight. The fact that there appears to be no selective blocker for SOCE and ROCE makes it difficult to evaluate experimentally the respective contributions of the two Ca 2z influxes during physiological conditions. These magnitudes are probably also cell-type dependent. Such issues explain the persistence of the controversy regarding SOCE and ROCE [45][46][47][48]. An informative experiment would be to stimulate ASMC using flash photolysis of caged IP 3 instead of agonist stimulation. Indeed, as IP 3 does not induce ROCE, SOCE should be the essential Ca 2z influx left. By comparing the responses to IP 3 stimulation in the presence and in the absence of extracellular calcium, one could then deduce the importance of SOCE in physiological conditions.

Efficacy of CPA
CPA is widely used as a SERCA blocker, having the advantage over Thapsigargin (Tg) of being reversible, and probably less toxic.
Both have been used extensively to study SOCE in different cell types (e.g., [25][26][27]43,49]). Although our work indicates that CPA does not fully block the SERCA in intact tissue such as lung slices, it does not imply that CPA should not be used experimentally to induce SOCE. Indeed, CPA might still cause substantial SOCE activation in the presence of agonist. However, our results indicate that CPA is not a good mean to fully empty Ca 2z stores, and care should be taken in interpreting the experimental results of its application. We suggest that a combined Rya-Caf treatment is a more reliable way to induce a permanent large SR depletion (Fig. 3B, C). There is evidence that Tg is an efficient SERCA blocker in cell lines such as Hela cells [43], but we have not addressed the effect of Tg on ASMC in lung slices in this study.

Modelling IPR
In this work, we followed the approach of Wang et al. [23], in that we have used one of the simplest models of IPR Ca 2z release, namely the Li-Rinzel/Tang et al. reduction of the DYK ODE model [34][35][36]. This category of IPR model produces agonistinduced Ca 2z oscillations characterised by significant SR Ca 2z depletion ( Fig. 3D and [23]), hence the possibility of SOCE being activated during such Ca 2z oscillations. This property might be model-dependent, however there is evidence that the SR is actually depleted to some extent during agonist-induced Ca 2z oscillations in ASMC. Indeed, the absence of effect of ryanodine during agonist-induced oscillations can be explained by the average level of ½Ca 2z SR being too low for RyR activation [1,23]. However, the respective ½Ca 2z SR ''thresholds'' for SOCE and RyR activation are experimentally unknown. In this work, the SOCE activation threshold was deduced from fitting the model simultaneously to Fig. 3A and Figs. 3B-C.
Finally, we note that our whole-cell ½Ca 2z i model would likely not benefit from using a recent Markov model of an IPR (e.g., [50][51][52]), because these models are based on steady-state data only (i.e., single-channel opening and closing times in stationary Ca 2z and IP 3 ) and typically miss the long inactivation timescale which was included ''ad hoc'' in the first IPR models to reproduce the observed behavior at the cell level (i.e., ½Ca 2z i oscillations upon agonist stimulation).

Limitations of the whole-cell model
As we are essentially interested in ½Ca 2z i responses of ASMC at the cell level, we have described Ca 2z dynamics via a deterministic ODE model. The scope of this model is, however, somewhat limited for the following reasons.
First, there is evidence that IPR are not homogeneously distributed on the SR membrane of cells, but are found as dense clusters. This channel clustering is especially patent upon stimulation by low agonist concentrations, for which local, stochastic Ca 2z releases may not propagate to neighboring clusters, resulting in spatially isolated, unsynchronised Ca 2z releases, called ''puffs''. At higher agonist concentrations, the frequency of these puffs increases, allowing Ca 2z releases from close sites to accumulate and propagate further away. This triggers, via CICR, the firing of more distant clusters, and results in Ca 2z waves propagating repeatedly throughout the cytosol. These waves usually appear as Ca 2z oscillations at the whole-cell level. While Ca 2z waves are indeed associated with ½Ca 2z i oscillations in ASMC [1], it has, so far, been impossible to detect Ca 2z puffs. This could arise from a less clustered distribution of IPR in ASMC, compared to the larger cells (ooycytes and Hela cells) where puffs have been characterised. On the other hand, Ca 2z ''sparks'', the equivalent of Ca 2z puffs but mediated by RyR, have been detected in ASMC [1], which supports a clustered distribution of RyR. In this study, we did not attempt to consider these spatial/stochastic aspects of the Ca 2z signals. Our model is thus less reliable at low agonist concentrations.
Second, cytoplasmic microdomains often exist between cell organelles (e.g., between peripheral SR and the plasma membrane, between the SR and mitochondria), out of which ½Ca 2z i cannot diffuse easily. These have consequences for SOCE dynamics. Indeed, it has been reported that upon store depletion, SERCA can colocalise with STIM proteins, in proximity to the PM [49,53]. As a consequence, if SOCE is slow enough, the SR can refill with Ca 2z without a concomitant increase in bulk ½Ca 2z i [49]. Upon large SOCE, this is no longer the case; however, mitochondria prevent the local Ca 2z increase to become too large by pumping Ca 2z from the subplasmalemmal space and releasing it deeper in the cytoplasm, where it can be absorbed by other SERCA [49]. These spatial effects cannot be accounted for by our current non-compartmentalised model.
Finally, Ca 2z dynamics are modified by Ca 2z buffers in the cytosol and SR, which bind 99% of the free Ca 2z . While the effect of fast, linear buffers can be taken into account by a global rescaling of Ca 2z fluxes (see Methods), this is not the case for high affinity buffers, in particular fluorescent dye indicators. Including such buffers in an ODE model of ½Ca 2z i dynamics leads to suppression of ½Ca 2z i oscillations, because the buffer affinity is close to the amplitude of whole-cell ½Ca 2z i oscillations. In reality, ½Ca 2z i reaches much higher levels locally upon IPR opening, so that the buffers become saturated and cannot prevent Ca 2z oscillations. Again, this would have to be accounted for by a spatial model of Ca 2z dynamics.

Future work
Although RyR dynamics play a role only during the initial phases of agonist-induced Ca 2z oscillations and Rya-Caf treatment, the interaction between RyR and IPR may become important in other situations, such as drug-induced RyR sensitisation. We plan to extend our model to these dynamics.
Since our work is part of a broader effort to improve the understanding of airway hyper-responsiveness and remodelling via mathematical modelling [54][55][56][57], we also intent to model the interaction of ASMC Ca 2z signalling with other aspects of lung dynamics. Although mathematical models of ASM contraction have previously been developed [54,55,58], modelling of other signalling pathways, such as inflammation and proliferation, is, to our knowledge, still in its infancy.
Additionally, experimental studies of ASMC inflammation and proliferation in conjunction with ½Ca 2z i imaging in lung slices would be desirable. While such studies have been carried out with cultured ASMC [3][4][5][9][10][11][12], they do not provide individual ½Ca 2z i dynamics; moreover, cultured ASMC often exhibit a different phenotype from ASMC in intact tissues.

Conclusions
The inclusion of SOCE in our mathematical model of Ca 2z dynamics in ASMC enables a better understanding of the experimental physiology of lung slices. It shows that the different abilities of CPA and Rya-Caf treatment to clamp the ½Ca 2z i of ASMC can be explained by their different ability to invoke SOCE. The model predicts that CPA, in contrast with Rya-Caf treatment, is unable to empty the SR because of its inefficiency to fully inhibit the SERCA. Furthermore, by accounting for both agonist-induced Ca 2z oscillations and SOCE activation by SR Ca 2z depletion, the model shows that SOCE can be a major determinant of the frequency of agonist-induced Ca 2z oscillations. Because this frequency of the Ca 2z oscillations regulates airway contraction, the model suggests a role for increased SOCE in AHR, a correlation consistent with SOCE up-regulation under inflammatory conditions typical of asthma. These predictions underscore the synergistic role for mathematical modeling in medical research.

Supporting Information
Supporting Information S1 Details of the parameter estimation procedure. (PDF)