Repository Citation This Work Is Licensed under a Creative Commons Attribution 4.0 License. a Deterministic Model Predicts the Properties of Stochastic Calcium Oscillations in Airway Smooth Muscle Cells

A deterministic model predicts the properties of stochastic calcium oscillations in airway smooth muscle cellsA deterministic model predicts the properties of stochastic calcium oscillations in airway smooth muscle cells" (2014). A deterministic model predicts the properties of stochastic calcium oscillations in airway smooth muscle cells Comments This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Abstract The inositol trisphosphate receptor (IP 3 R) is one of the most important cellular components responsible for oscillations in the cytoplasmic calcium concentration. Over the past decade, two major questions about the IP 3 R have arisen. Firstly, how best should the IP 3 R be modeled? In other words, what fundamental properties of the IP 3 R allow it to perform its function, and what are their quantitative properties? Secondly, although calcium oscillations are caused by the stochastic opening and closing of small numbers of IP 3 R, is it possible for a deterministic model to be a reliable predictor of calcium behavior? Here, we answer these two questions, using airway smooth muscle cells (ASMC) as a specific example. Firstly, we show that periodic calcium waves in ASMC, as well as the statistics of calcium puffs in other cell types, can be quantitatively reproduced by a two-state model of the IP 3 R, and thus the behavior of the IP 3 R is essentially determined by its modal structure. The structure within each mode is irrelevant for function. Secondly, we show that, although calcium waves in ASMC are generated by a stochastic mechanism, IP 3 R stochasticity is not essential for a qualitative prediction of how oscillation frequency depends on model parameters, and thus deterministic IP 3 R models demonstrate the same level of predictive capability as do stochastic models. We conclude that, firstly, calcium dynamics can be accurately modeled using simplified IP 3 R models, and, secondly, to obtain qualitative predictions of how oscillation frequency depends on parameters it is sufficient to use a deterministic model. Copyright: ß 2014 Cao et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability: The authors confirm that all data underlying the findings are fully …


Introduction
Oscillations in cytoplasmic calcium concentration (½Ca 2z i ), mediated by inositol trisphosphate receptors (IP 3 R; a calcium channel that releases calcium ions (Ca 2z ) from the endoplasmic or sarcoplasmic reticulum (ER or SR) in the presence of inositol trisphosphate (IP 3 )) play an important role in cellular function in many cell types. Hence, a thorough knowledge of the behavior of the IP 3 R is a necessary prerequisite for an understanding of intracellular Ca 2z oscillations and waves. Mathematical and computational models of the IP 3 R play a vital role in studies of Ca 2z dynamics. However, over the past decade, two major questions about IP 3 R models have arisen.
Firstly, how best should the IP 3 R be modeled? Models of the IP 3 R have a long history, beginning with the heuristic models of [1][2][3]. With the recent appearance of single-channel data from IP 3 R in vivo [4,5], a new generation of Markov IP 3 R models has recently appeared [6,7]. These models show that IP 3 R exist in different modes with different open probabilities. Within each mode there are multiple states, some open, some closed. Importantly, it was found [8] that time-dependent transitions between different modes are crucial for reproducing Ca 2z puff data from [9]. However, it is not yet clear whether transitions between states within each mode are important, or whether all the important behaviors are captured simply by inter-mode transitions.
Secondly, why do deterministic models of the IP 3 R perform so well as predictive models? Deterministic models of the IP 3 R have proven to be useful predictive models in a range of cell types. For example, IP 3 R-based models have been developed to study Ca 2z oscillations in airway smooth muscle cells (ASMC) [10][11][12][13], and these models have made predictions which have been confirmed experimentally. This shows the usefulness of such models in advancing our understanding of how intracellular Ca 2z oscillations and waves are initiated and controlled in ASMC. However, these models are deterministic models which assume infinitely many IP 3 R per unit cell volume, an assumption that contradicts experimental findings in many cell types showing that Ca 2z puffs and spikes occur stochastically, and that intracellular Ca 2z waves and oscillations arise as an emergent property of fundamental stochastic events [9,14,15].
Here, we answer these two fundamental modeling questions using data and models from ASMC. Firstly, we show that a simple model of the IP 3 R, involving only two states with time-dependent transitions, suffices to generate correct dynamics of Ca 2z puffs and oscillations. Secondly, we show that, although Ca 2z oscillations in ASMC are generated by a stochastic mechanism, a deterministic model can make the same qualitative predictions as the analogous stochastic model, indicating that deterministic models, that require much less computational time and complexity, can be used to make reliable predictions. Although we work in the specific context of ASMC, our results are applicable to other cell types that exhibit similar Ca 2z oscillations and waves.

Results
A two-state model of the IP 3 R is sufficient to reproduce function We have previously shown [8] that the statistics of Ca 2z puffs in SH-SY5Y cells can be reproduced by a Markov model of the IP 3 R based on the steady-state data of [5] and the time-dependent data of [4]. In this model the IP 3 R can exist in 6 different states, grouped into two modes, which we call Drive and Park (see In our previous study on calcium puffs [8], we showed that, to reproduce the experimentally observed non-exponential interspike interval (ISI) distribution and coefficient of variation (CV) of ISI smaller than 1, the time-dependent intermodal transitions are crucial. Lack of time dependencies in the Siekmann model leads to exponential ISI distributions and CV = 1, which is not the case for calcium spikes in ASMC. Fig. 2A shows an example of Ca 2z oscillations generated by 50 nM methacholine (MCh, an agonist that can induce the production of IP 3 by binding to a G proteincoupled receptor in the cell membrane) in ASMC. By gathering data from 14 cells in 5 mouse lung slices, we found that the standard deviation of the interspike interval (ISI) is approximately a linear function of the ISI mean, with a slope clearly between 0 and 1 (i.e. CVv1), indicating that the spikes are generated by an inhomogeneous Poisson process (a slope of 1 would denote a pure Poisson process) (see Fig. 2B). This shows the necessity of inclusion of time-dependent transitions for mode-switching.
Using a quasi-steady-state approximation, and ignoring states with very low dwell times, it is possible to construct a simplified two-state version of the full six-state model (see Materials and Methods). In the simplified model the intramodal structure is ignored, and only the intermodal transitions have an effect on IP 3 R behavior. In Fig. 3 we compared the simplified IP 3 R model to the full six-state model. Both models have the same distribution of interspike interval, spike amplitude and spike duration. Moreover, by looking at a more detailed comparison between the two model results (Figs. 4A, C and E) and experimental data (Figs. 4B, D and F), we found the 2-state model not only can reproduce the behaviour of the 6-state model, but can also qualitatively reproduce experimental data. The average experimental ISI shows a clear decreasing trend as MCh concentration increases (although a saturation occurs in the data for high MCh), a trend that is mirrored by the model results as the IP 3 concentration increases. Unfortunately, since the exact relationship between MCh concentration and IP 3 concentration is uncertain, a quantitative comparison is not possible. In both model and experimental results, the average peak and duration of the oscillations are nearly independent of agonist concentration. The quantitative difference in spike duration between the model results and the data in Figs. 4E and F are most likely due to choice of calcium buffering parameters. For example, adding 3 or 5 mM fast Ca 2z buffer (see Materials and Methods) increases the average spike duration to 0.54 s or 0.7 s respectively, which are close to the levels shown in the data.
Thus, the intramodal structure of the six-state model is essentially unimportant, as the model behavior (in terms of the qs are rates of state-transitions between two adjacent states and q 42 and q 24 are transitions between the two modes [7]. doi:10.1371/journal.pcbi.1003783.g001

Author Summary
The inositol trisphosphate receptor (IP 3 R) is one of the most important cellular components responsible for calcium oscillations. Over the past decade, two major questions about the IP 3 R have arisen. Firstly, what fundamental properties of the IP 3 R allow it to perform its function? Secondly, although calcium oscillations are caused by the stochastic properties of small numbers of IP 3 R is it possible for a deterministic model to be a reliable predictor of calcium dynamics? Using airway smooth muscle cells as an example, we show that calcium dynamics can be accurately modeled using simplified IP 3 R models, and, secondly, that deterministic models are qualitatively accurate predictors of calcium dynamics. These results are important for the study of calcium dynamics in many cell types.
statistics of puffs and oscillations) is governed almost entirely by the time dependence of the intermode transitions, particularly the time dependence of the rapid inhibition of the IP 3 R by high ½Ca 2z i , and the slow recovery from inhibition by Ca 2z . The multiple states within each mode are necessary to obtain an acceptable quantitative fit to single-channel data, but are nevertheless of limited importance for function. Hence, even when simulating microscopic events such as Ca 2z puffs it is sufficient to use a simpler, faster, two-state model, rather than a more complex six-state model. In the following, we will use the 2-state IP 3 R model to generate all the simulation results.

Prediction of stochastic Ca 2z behavior by a deterministic model
Although the data (Fig. 2) show that Ca 2z oscillations in ASMC are generated by a stochastic process, not a deterministic one, we wish to know to what extent a deterministic model can be used to make qualitative (and experimentally testable) predictions. Our   . More detailed comparisons between the 2-state and the 6-state IP 3 R models, and a comparison to experimental data. As a function of IP 3 concentration (p), the two models give the same ISI (A), peak ½Ca 2z i (C) and spike duration (E). These results agree qualitatively with experimental data, as shown in panels B, D and F respectively. Quantitative comparisons are generally not possible as the relationship between IP 3 concentration and agonist concentration is not known. Error bars represent mean + SEM. Data for each MCh concentration are obtained from at least three different cells from at least two different lung slices. doi:10.1371/journal.pcbi.1003783.g004 simplified 2-state Markov model of the IP 3 R can be converted to a deterministic model (see Materials and Methods). The result is a system of ordinary differential equations (ODEs) with four variables, which takes into account the increased ½Ca 2z i at an open IP 3 R pore, as well as the increased ½Ca 2z i within a cluster of IP 3 R; the four variables are the ½Ca 2z i outside the IP 3 R cluster (c), the ½Ca 2z i within the IP 3 R cluster (c b ), the total intracellular Ca 2z concentration (c t ) and an IP 3 R gating variable (h 42 ). We refer to the reduced 4D model as the deterministic model for all the results and analyses.
Note that there is no physical or geometric constraint enforcing a high local ½Ca 2z i ; in this case the spatial heterogeneity arises solely from the low diffusion coefficient of Ca 2z . Our use of c b is merely a highly simplified way of introducing spatial heterogeneity of the Ca 2z concentration. Since the IP 3 R can only ''see'' c b (as well as the Ca 2z concentration right at the mouth of an open channel, which we denote by c p ), but cannot be influenced directly by c (the experimentally observed Ca 2z signal), our approach allows for the functional differentiation of the rapid local oscillatory Ca 2z in the cluster, from the slower Ca 2z signal in the cytoplasm, without the need for computationally intensive simulations of a partial differential equation model. Quantitative accuracy is thus sacrificed for computational convenience.
Calcium oscillations in the stochastic and deterministic models are shown in Fig. 5A. According to our previous results [8], the average value of h 42 over the cluster of IP 3 R primarily regulates the termination and regeneration of individual spikes. This can be seen in the stochastic model by projecting the solution on the c b , h 42 phase plane (Fig. 5B). Upon an initial Ca 2z release from one or more IP 3 R, a large spike is generated by Ca 2+ -induced Ca 2z release (via the IP 3 R) during which time a decreasing h 42 gradually decreases the average open probability of the clustered IP 3 R. The spike is terminated when h 42 is too small to allow further Ca 2z release. This phenomenon is qualitatively reproduced by the deterministic model (Fig. 5D). In both the stochastic and deterministic models the decrease in average IP 3 R open probability of a cluster of IP 3 R caused by Ca 2z inhibition is the main reason for the termination of each spike.
According to Figs. 5B and D, regeneration of each spike requires a return of h 42 back to a relatively high value (i.e., recovery of the IP 3 R from inhibition by Ca 2z ). The deterministic model sets a clear threshold for the regeneration, as can be seen in Fig. 5C, where an upstroke in c b occurs when the trajectory creeps beyond the sharp ''knee'' of the white curve. When the trajectory reaches the knees of the white curve it is forced to jump across to the other stable branch of the critical manifold, resulting in a fast increase in c b followed by a relatively fast increase in c (seen by combining Figs. 5C and D).
In contrast, the stochastic model enlarges the contributions of individual IP 3 R so that the generation of each spike is also effectively driven by random Ca 2z release through the IP 3 R, which can be seen in the inset of Fig. 5B where the site of spike initiation (blue bar) exhibits significantly greater variation than that of spike termination (green bar). In spite of this, the essential similarities in phase plane behavior result in both deterministic and stochastic models making the same qualitative predictions in response to perturbations, such as changes in IP 3 concentration (½IP 3 ), Ca 2z influx or efflux. In the following, we illustrate this by investigating a number of experimentally testable predictions. Due to the extensive importance of frequency encoding in many Ca 2zdependent processes, we focus particularly on the change of oscillation frequency in response to parameter perturbations. As a side issue we also investigate how the oscillation baseline depends on physiologically important parameters.

Dependence of oscillation frequency on IP 3 concentration
In many cell types a moderate increase in ½IP 3 increases the Ca 2z oscillation frequency (see Fig. 2A in [11], Fig. 4E in [16] and Fig. 6B in [17]), a result that is reproduced by both model types (Fig. 6A). As ½IP 3 increases, the stochastic model increases the probability of the initial Ca 2z release through the first open IP 3 R and of the following Ca 2z release, thus shortening the average ISI. Although the oscillatory region of the deterministic model is strictly confined by bifurcations which do not apply to the stochastic model, the deterministic model can successfully replicate an increasing frequency by lowering the ''knee'' of the red curve in Fig. 5D and shortening the time spent from the termination point c to the initiation point a (thus shortening the ISI). Hence, although the deterministic model cannot be used to predict the exact values of ½IP 3 at which the oscillations begin and end, as stochastic effects predominate in these regions, it can be used to predict the correct qualitative trend in oscillation frequency.

Dependence of oscillation frequency on Ca 2z influx and efflux
In many cell types, including ASMC, transmembrane fluxes modulate the total intracellular Ca 2z load (c t ) on a slow time scale [16,18], and thereby modulate the oscillation frequency [19]. Experimental data can be seen in Fig. 8 in [16] and Fig. 2 in [18]. Figs. 6B and C show that both stochastic and deterministic models predict the same qualitative changes in oscillation frequency in response to changes in membrane fluxes (through membrane ATPase pumps and/or Ca 2z influx channels such as receptoroperated channels or store-operated channels).

Dependence of oscillation frequency on SERCA expression
The level of sarco/endoplasmic reticulum calcium ATPase (SERCA) expression (or capacity) is important for airway remodeling in asthma [20] and ASMC Ca 2z oscillations [21]. We thus investigated the predictions of the two models in response to changes in SERCA expression (V s ). As V s decreases, the deterministic model exhibits a decreasing frequency, in agreement with experimental data (see Figs. 3 and 4 in [21]). The same trend is seen in the stochastic model with only 20 IP 3 R (see Fig. 6D).

Dependence of oscillation frequency on Ca 2z buffer concentration
Calcium buffers have been shown to be able to change the ISI and spike duration, which in turn change the oscillation frequency [15,22]. We compared the effects on the two models of varying total buffer concentration (B t ) by adding one buffer with relatively fast kinetics to the models (see Materials and Methods for details). In both models the frequency decreases as B t increases (see Fig. 6E), which is consistent with experimental data (Fig. 2B in [18]). This is not surprising, because increasing B t can decrease the effective rates of SR Ca 2z release and reuptake.
Dependence of oscillation baseline on Ca 2z influx and

SERCA expression
Sustained elevations of baseline during agonist-induced Ca 2z oscillations or transients have been observed experimentally, and are believed to be a result of an increase in Ca 2z influx caused by opening of membrane Ca 2z channels [13,16]. Furthermore, there is evidence showing that decreased SERCA expression could also increase the baseline (Fig. 4 in [21]). Those phenomena are successfully reproduced by both models (see Fig. 7).

Discussion
In this paper we address two current major questions in the field of Ca 2z modeling. Firstly, we show that Ca 2z puffs and stochastic oscillations can be reproduced quantitatively by an extremely simple model, consisting only of two states (one open, one closed), with time-dependent transitions between them. This model is obtained by removing the intramodal structure of a more complex model that was determined by fitting a Markov model to single-channel data [7]. We thus show that the internal structure of each mode is irrelevant for function and mode switching is the key mechanism for the control of calcium release. The necessity for time-dependent mode switching is shown not only by the dynamic single-channel data of [4]), but also by the puff data of [9] and our ASMC data.
Secondly, we investigate the role of stochasticity of IP 3 R in modeling Ca 2z oscillations in ASMC by comparing a stochastic IP 3 R-based Ca 2z model and its associated deterministic version, for parameters such that both of the models exhibit Ca 2z spikes but the stochastic model cannot necessarily be replaced by a mean-field model. We find that a four-variable deterministic model has the same predictive power as the stochastic model, in  that it correctly reproduces the process of spike termination and predicts the same qualitative changes in oscillation frequency and baseline in response to a variety of perturbations that are commonly used experimentally. The mechanism for termination of individual spikes is fundamentally a deterministic process controlled by a rapid inhibition induced by the high local ½Ca 2z i in the IP 3 R cluster, whereas spike initiation is significantly affected by stochastic opening of IP 3 R. Hence, repetitive Ca 2z cycling is primarily induced by the time-dependent gating variables governing transitions of the IP 3 R from one mode to another.
Our simplified two-state model of the IP 3 R is identical in structure (although not in parameter values) to the well-known model of [23]. It is somewhat ironic that after 20 years of detailed studies of the IP 3 R and the construction of a plethora of models of varying complexity, the single-channel data have led us around full circle, back to these original formulations. Excitability is arising via a fast activation followed by a slower inactivation, a combination often seen in physiological processes [24]. Encoding of this fundamental combination results directly from the two-mode structure of the IP 3 R. Although similar single-channel data have been used to construct three-mode models [6,25], neither of these models has yet been used in detailed studies of Ca 2z puffs and waves, and it remains unclear whether or not they have a similar underlying structure.
In contrast to previous deterministic ODE models, our fourvariable Ca 2z model includes a more accurate IP 3 R model, as well as local control of clustered IP 3 R by two distinct Ca 2z microdomains; one at the mouth of an open IP 3 R, the other inside a cluster of IP 3 R. Neglect of either of these microdomains leads to models that either exhibit unphysiological cytoplasmic Ca 2z concentrations or fail to reproduce reasonable oscillations. This underlines the importance of taking Ca 2z microdomains into consideration when constructing any model. Our microdomain model is highly simplified, with the microdomain being treated simply as a well-mixed compartment. More detailed modeling of spatially-dependent microdomains is possible, and not difficult in principle, but requires far greater computational resources. It is undeniable that a more detailed model, incorporating the full spatial complexity -and possibly stochastic aspects as well -would make, overall, a better predictive tool. However, our goal is to find the simplest models that can be used as predictive tools.
An important similar study is that of Shuai and Jung [26]. They compared the use of Markov and Langevin approaches to the computation of puff amplitude distributions, compared their results with the deterministic limit, and showed that IP 3 R stochasticity does not qualitatively change the type of puff amplitude distribution except for when there are fewer than 10 IP 3 R. Here, we significantly extend the scope of their study by exploring the effects of IP 3 R stochasticity on the dynamics of Ca 2z spikes, and we do this in the context of an IP 3 R model that has been fitted to single-channel data. Although this is true in a general sense for the Li-Rinzel model, which is based on the DeYoung-Keizer model, which did take into account the opening time distributions of IP 3 R in lipid bilayers, neither model can reproduce the more recent data obtained from on-nuclei patch clamping. When these recent data are taken into account one obtains a model with the same structure, but quite different parameters and behavior.  We find that, in spite of a relatively large variation in spike amplitude which is partially caused by a large variation in ISI (Fig. 5B), the mechanism governing individual spike terminations is the same for both a few or infinitely many IP 3 R, which explains why the one-peak type of amplitude distribution is independent of the choice of IP 3 R number (see Fig. 6A in [26]).
Another important relevant study was done by Dupont et al. [27], who compared the regularity of stochastic oscillations in hepatocytes for different numbers of IP 3 R clusters. They found that the impact of IP 3 R stochasticity on global Ca 2z oscillations (in terms of CV) increases as the total cluster number decreases. Our study here extends these results, and demonstrates how well stochastic oscillations can be qualitatively described by a deterministic system, even when there is only a small number of IP 3 R (which appears to be the case for ASMC, in which the wave initiation site is only 2*4 mm in diameter). Indeed, as we have shown, for the purposes of predictive modeling a simple deterministic model does as well as more complex stochastic simulations.
Ryanodine receptors (RyR) are another important component modulating ASMC Ca 2z oscillations [16,28,29] but are not included in our model. This is because the role of RyR is not fully understood and may be species-dependent; for example, in mouse or human ASMC, RyR play very little role in IP 3 -induced continuing Ca 2z oscillations [17,30], but this appears not to be true for pigs [28]. Our study focuses on the calcium oscillations in mouse and human (as we did in our experiments) where inclusion of a deterministic model of RyR should have little effect. An understanding of the role of RyR stochasticity and how the IP 3 R and the RyR interact needs a reliable RyR Markov model, exclusive to ASMC, which is not currently available. Multiple Markov models of the RyR have been developed for use in cardiac cells [31], but these are based on single-channel data from lipid bilayers, and are adapted for the specific context of cardiac cells. Their applicability to ASMC remains unclear.
Although we have not shown that the deterministic model for ASMC has the same predictive power as the stochastic model in all possible cases (which would hardly be possible in the absence of an analytical proof) the underlying similarity in phase plane structure indicates that such similarity is plausible at least. Certainly, we have not found any counterexample to this claim. However, whether or not this claim is true for all cell types is unclear. Some cell types exhibit both local Ca 2z puffs and global Ca 2z spikes (usually propagating throughout the cells in the form of traveling waves), showing that initiation of such Ca 2z spikes requires a synchronization of Ca 2z release from more than one cluster of IP 3 R [14]. This type of spiking relies on the hierarchical organization of Ca 2z signal pathways, in particular the stochastic recruitment of both individual IP 3 R and puffs at different levels [32], and therefore cannot be simply reproduced by deterministic models containing only a few ODEs. However, Ca 2z oscillations in ASMC, as observed in lung slices, may not be of this type, as IP 3 R-dependent puffs have not been seen in these ASMC. It thus appears that, in ASMC in lung slices, every Ca 2z ''puff'' initiates a wave, resulting in periodic waves with ISI that are governed by the dynamics of individual puffs.

Ethics Statement
Animal experimentations carried out were approved by the Animal Care and Use Committee of the University of Massachusetts Medical School under approval number A-836-12.

Lung slice preparation
BALB/c mice (7-10 weeks old, Charles River Breeding Labs, Needham, MA) were euthanized via intraperitoneal injection of 0.3 ml sodium pentabarbitone (Oak Pharmaceuticals, Lake Forest, IL). After removal of the chest wall, lungs were inflated with *1:1 ml of 1.8% warm agarose in sHBSS via an intratracheal catheter. Subsequently, air (*0:3 ml) was injected to push the agarose within the airways into the alveoli. The agarose was polymerized by cooling to 4 0 C. A vibratome (VF-300, Precisionary Instruments, San Jose, CA) was used to make 180 mm thick slices which were maintained in Dulbecco's Modified Eagle's Media (DMEM, Invitrogen, Carlsbad, CA) at 37 0 C in 10% CO 2 /air. All experiments were conducted at 37 0 C in a custom-made temperature-controlled Plexiglas chamber as described in [17].

The calcium model
Inhomogeneity of cytoplasmic Ca 2z concentration not only exists around individual channel pores of the IP 3 R, where a nearly instantaneous high Ca 2z concentration at the pore (denoted by c p ) leads to a very sharp concentration profile, but is also seen inside an IP 3 R cluster where the average cluster Ca 2z concentration (c b ) is apparently higher than that of the surrounding cytoplasm (c) [33]. This indicates that during Ca 2z oscillations each IP 3 R is controlled by either the pore Ca 2z concentration (when it is open) or the cluster Ca 2z concentration (when it is closed). Neither of these local concentrations influence cell membrane fluxes or the majority of SERCAs, which we assume to be distributed outside the cluster.
The scale separation between the pore Ca 2z concentration and the cluster Ca 2z concentration allows to treat c p as a parameter, providing a simpler way of modeling local Ca 2z events (like Ca 2z puffs) that has been used in several previous studies [8,34,35]. However, evolution of the cluster concentration and wide-field cytoplasm Ca 2z concentration are not always separable, so an additional differential equation for the cluster Ca 2z is necessary.
A schematic diagram of the model is shown in Fig. 8. The corresponding ODEs are where c t~c zc b =c 1 zc s =c 2 representing total intracellular Ca 2z concentration, and thus SR Ca 2z concentration, c s is given by c s~c2 (c t {c{c b =c 1 ). c 1 and c 2 are the volume ratios given in Table 1. J IPR is the flux through the IP 3 R, J leak is a background Ca 2z leak out of the SR, and J serca is the uptake of Ca 2z into the SR by SERCA pumps. J pm is the flux through plasma pump, and J in represents a sum of main Ca 2z influxes including J rocc (receptor-operated Ca 2z channel), J socc (store-operated Ca 2z channel) and J leakin (Ca 2z leak into the cell). J diff coarsely models the diffusion flux from cluster microdomain to the cytoplasm. Calcium concentration at open channel pore (c p ) does not explicitly appear in the equations but is used in the IP 3 R model introduced later. c p is assumed to be proportional to SR Ca 2z concentration (c s ) and is therefore simply modeled by c p~cp0 (c s =100) where c p0 is the value corresponding to c s~1 00 mM. Alternatively, c p can also be assumed to be a large constant (say greater than 100 mM) without fundamentally altering the model dynamics. The choice of c p0 is not critical as long as it is sufficiently large to play a role in inactivating the open channels. All the parameter values are given in Table 1.

The data-driven IP 3 R model
The IP 3 R model used in our ASMC calcium model is an improved version of the Siekmann IP 3 R model which is a 6-state Markov model derived by fitting to the stationary single channel data using Markov chain Monte Carlo (MCMC) [5,7,8]. Fig. 1 has shown the structure of the IP 3 R model which is comprised of two modes; the drive mode, containing three closed states C 1 , C 2 , C 3 and one open state O 6 , and the park mode, containing one closed state C 4 and one open state O 5 . The transition rates in each mode are constants (shown in Table 2), but q 42 and q 24 which connect the two modes are Ca 2z -/IP 3 -dependent and are formulated as   Figure 9. Stationary data and fits of q 24 and q 42 . Stationary transition rates of q 24 and q 42 , q 24? and q 42? , as functions of Ca 2z concentration were estimated and fitted for two ½IP 3 , 1 mM (A) and 10 mM (B). Circles and squares represent the means of q 24 and q 42 distributions computed by MCMC simulation [7]. Note that MCMC failed to determine the values of q 24 and q 42 atĉ c~1, 10 mM for 10 mM IP 3 , as the IP 3 R was almost in the drive mode for these cases. The corresponding fitting curves (solid for q 42 ; dashed for q 24 ) are produced using Eqs.
The expressions of as, V s, ns and ks are chosen as follows so that Eq. 11 and Eq. 12 capture the correct trends of experimental values of q 24 and q 42 (see Fig. 9) and generate relatively smooth open probability curves (see Fig. 10), V 24~6 2z880=(p 2 z4) a 24~1 z5=(p 2 z0:25) V 42~1 10p 2 =(p 2 z0:01) a 42~1 :8p 2 =(p 2 z0:34) k 24~0 :35 k 42~0 :49z0:543p 3 =(p 3 z64) Note that the above formulas are different from the relatively complicated formulas used in [8]. The rates, l m24 , l h24 and l m42 , are constants estimated by using dynamic single channel data [4] and given in Table 2, whereas l h42 is not clearly revealed by experimental data. However we have shown that it should be relatively large for highĉ c but relatively small for lowĉ c for reproducing experimental puff data [8]. By introducing two Ca 2z concentrations, c b and c p , l h42 and the state of the IP 3 R channel become highly correlated, so that we can assume l h42 is a relatively large value H if the channel is open and is a relatively small value L if the channel is closed. Hence, l h42 is modeled by the logic function The IP 3 R model reduction Here we reduce the 6-state model to a 2-state open/closed model. The reduction takes the following steps: N The sum of the probabilities of C 1 , C 3 and O 5 is less than 0.03 for anyĉ c, so they are either rarely visited by the IP 3 R or have a very short dwell time. This implies they have very little contribution to the Ca 2z dynamics. Therefore, we completely remove the three states from the full model. N Transition rates of q 26 and q 62 are about 2 orders larger than that of q 24 and q 42 , which allows us to omit the fast transitions by taking a quasi-steady state approximation. This change will affect two aspects. First, we have O 6~C2 q 26 =q 62 which allows us to combine C 2 and O 6 to be a new state D, which satisfies D~O 6 (q 62 zq 26 )=q 26    N To describe an average rate that infinitely many receptors are rapidly inhibited by high Ca 2z concentration but slowly restored from Ca 2z -inhibition. l h42 is proposed to be N l h 42~( 1{D)LzDH:

Deterministic formulation of the stochastic model
Based on the above changes, the full deterministic model containing 8 ODEs is presented as follows, where q 24 and q 42 are functions of the gating variables given by Eqs. 4

Reduction of the full deterministic model
The full deterministic model contains 8 variables which make the model difficult to implement and analyze. Thus, we reduce the full model to a minimal model that still captures the crucial features of the full model. First of all, l m42 , l m24 and l h24 are sufficiently large so that we can assume they instantaneously follow their equilibrium functions. Therefore, by taking quasi-steady state approximation to m 24 , h 24 and m 42 , we remove the three timedependent variables from the full model.
By now, the full model has been reduced to a 5D model, Second, the rate of change of D approaching its equilibrium, l D~( q 42 q 62 zq 42 q 26 zq 24 q 62 )=(q 62 zq 26 ) (calculated from Eq. 24), is at least one order larger than those of c, c t and h 42 , indicating that taking the quasi-steady state approximation to Eq. 24 could not significantly affect the evolutions of c, c t and h 42 . That is, D~q 42 (q 62 zq 26 ) q 42 q 62 zq 42 q 26 zq 24 q 62 : We emphasize here that the theory of the quasi-steady state approximation has not yet been well established, particularly about the rigorous conditions under which such a reduction is valid. Thus, our criterion of judging the validity of the reduction is checking whether the solutions of the reduced model are capable of qualitatively reproducing that of its original model. For this model, we find the reduction works. Hence, the full model is eventually reduced to a 4D model summarized as follows, dh 42 dt~l h 42 ( where D is given by Eq. 26.

Inclusion of calcium buffers
To check the effect of calcium buffers on oscillation frequency, we introduce a stationary buffer (no buffer diffusion), as mobile buffers are too complicated to be included in the current deterministic model. Since we have two different cytoplasmic Ca 2z concentrations, c and c b , two pools of buffer with the same kinetics should be considered. Hence, the inclusion of a stationary calcium buffer is modeled by the following system, where b (b 1 and b 2 ) and B t represent the concentrations of Ca 2zbound buffer and total buffer respectively. k z and k { are the rates of Ca 2z -binding and Ca 2z -dissociation, indicating how fast the time scale of the buffer dynamics is. Fast buffer refers to the buffer with relatively large k z . In the simulations, we use a fast buffer with k z~1 00 mM {1 : s {1 and k {~1 00 s {1 and vary B t to test if the stochastic model and the deterministic model have a qualitatively similar B t -dependency. Results are given in Fig. 6E.

Numerical methods and tools for deterministic and stochastic simulations
For the stochastic model, Eqs. 1-3 and ODEs of the four gating variables in the IP 3 R model are solved by the fourth-order Runge-Kutta method (RK4) and the stochastic states of IP 3 R determined by the IP 3 R model are solved by using a hybrid Gillespie method with adaptive timing [37]. The maximum time step size is set to be either 10 {4 s (for the 6-state IP 3 R model) or 10 {3 s (for the reduced 2-state IP 3 R model). All the computations are done with MATLAB (The MathWorks, Natick, MA) and the codes are provided in Supporting information (Text S1-S2). For the deterministic model, we use ode15s, an ODE solver in MATLAB. Accuracy is controlled by setting an absolute tolerance of 10 {8 applied to all the variables.

Statistical analysis
Data analysis is performed on the Ca 2z traces with relatively stable baselines and less noise. A moving average of every 3 data points is used to improve the data by smoothing out short-term fluctuations ( Fig. 2A is an improved result). Due to large variations in baseline, amplitude, and level of noise in data, we used two thresholds to get samples: a low threshold, 20% of the amplitude of the largest spike above the baseline, to initially filter baseline noise out; and a relatively high threshold, 50% of the amplitude of the largest spike above the baseline, to further remove small spikes that cannot initiate waves. For simulated stochastic traces of variable c, we first convert it to fluorescence ratio (F =F 0 ) by using F =F 0~c (c 0 zK d )=(c 0 czc 0 K d ) where the dissociation constant of Oregon Green K d~0 :17 mM and resting ½Ca 2z i c 0~0 :1 mM. We then used the same sampling procedure mentioned above to obtain samples. After samples are chosen, ISIs and spike durations are calculated based on the low threshold. Simulated traces used to calculate average frequency are about 200-400 seconds long. All the samplings and linear least-squares fittings are implemented using MATLAB (see Text S3-S4 for Matlab codes).

Supporting Information
Dataset S1 ASMC calcium fluorescence trace data. Text S1 Matlab code for simulation using 6 state IP 3 R model.

(DOCX)
Text S2 Matlab code for simulation using 2 state IP 3 R model.