Abstract
The inositol trisphosphate receptor (
) 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 
have arisen. Firstly, how best should the 
be modeled? In other words, what fundamental properties of the 
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 
, 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 
, and thus the behavior of the 
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, 
stochasticity is not essential for a qualitative prediction of how oscillation frequency depends on model parameters, and thus deterministic 
models demonstrate the same level of predictive capability as do stochastic models. We conclude that, firstly, calcium dynamics can be accurately modeled using simplified 
models, and, secondly, to obtain qualitative predictions of how oscillation frequency depends on parameters it is sufficient to use a deterministic model.
Citation: Cao P, Tan X, Donovan G, Sanderson MJ, Sneyd J (2014) A Deterministic Model Predicts the Properties of Stochastic Calcium Oscillations in Airway Smooth Muscle Cells. PLoS Comput Biol 10(8):
e1003783.
https://doi.org/10.1371/journal.pcbi.1003783
Editor: Andrew D. McCulloch, University of California, San Diego, United States of America
Received: April 6, 2014; Accepted: June 24, 2014; Published: August 14, 2014
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 available without restriction. All relevant data are within the paper and its Supporting Information files.
Funding: Funding for all authors came from National Heart Lung Blood Institute (USA) RO1 HL103405. http://www.nhlbi.nih.gov/ The funders had no role in study design, data collection and analysis, decisions to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Oscillations in cytoplasmic calcium concentration (
), mediated by inositol trisphosphate receptors (
; a calcium channel that releases calcium ions (
) from the endoplasmic or sarcoplasmic reticulum (ER or SR) in the presence of inositol trisphosphate (
)) play an important role in cellular function in many cell types. Hence, a thorough knowledge of the behavior of the 
is a necessary prerequisite for an understanding of intracellular 
oscillations and waves. Mathematical and computational models of the 
play a vital role in studies of 
dynamics. However, over the past decade, two major questions about 
models have arisen.
Firstly, how best should the 
be modeled? Models of the 
have a long history, beginning with the heuristic models of [1]–[3]. With the recent appearance of single-channel data from 
in vivo [4], [5], a new generation of Markov 
models has recently appeared [6], [7]. These models show that 
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 
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 
perform so well as predictive models? Deterministic models of the 
have proven to be useful predictive models in a range of cell types. For example, 
-based models have been developed to study 
oscillations in airway smooth muscle cells (ASMC) [10]–[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 
oscillations and waves are initiated and controlled in ASMC. However, these models are deterministic models which assume infinitely many 
per unit cell volume, an assumption that contradicts experimental findings in many cell types showing that 
puffs and spikes occur stochastically, and that intracellular 
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 
, involving only two states with time-dependent transitions, suffices to generate correct dynamics of 
puffs and oscillations. Secondly, we show that, although 
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 
oscillations and waves.
Results
A two-state model of the 
is sufficient to reproduce function
We have previously shown [8] that the statistics of 
puffs in SH-SY5Y cells can be reproduced by a Markov model of the 
based on the steady-state data of [5] and the time-dependent data of [4]. In this model the 
can exist in 6 different states, grouped into two modes, which we call Drive and Park (see Fig. 1). The Drive mode (which contains 4 states; 1 open and 3 closed) has an average open probability of around 0.7, while the Park mode (which contains the remaining two states; 1 open and 1 closed) has an open probability close to zero. Transitions between states within each mode are independent of 
and 
; only the transitions between modes are ligand-dependent.
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 
oscillations generated by 50 nM methacholine (MCh, an agonist that can induce the production of 
by binding to a G protein-coupled 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. 
), 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 
behavior. In Fig. 3 we compared the simplified 
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 
concentration increases. Unfortunately, since the exact relationship between MCh concentration and 
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 
fast 
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 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 
by high 
, and the slow recovery from inhibition by 
. 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 
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 
model to generate all the simulation results.
Prediction of stochastic 
behavior by a deterministic model
Although the data (Fig. 2) show that 
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 simplified 2-state Markov model of the 
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 
at an open 
pore, as well as the increased 
within a cluster of 
; the four variables are the 
outside the 
cluster (
), the 
within the 
cluster (
), the total intracellular 
concentration (
) and an 
gating variable (
). 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 
; in this case the spatial heterogeneity arises solely from the low diffusion coefficient of 
. Our use of 
is merely a highly simplified way of introducing spatial heterogeneity of the 
concentration. Since the 
can only “see” 
(as well as the 
concentration right at the mouth of an open channel, which we denote by 
), but cannot be influenced directly by 
(the experimentally observed 
signal), our approach allows for the functional differentiation of the rapid local oscillatory 
in the cluster, from the slower 
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 
over the cluster of 
primarily regulates the termination and regeneration of individual spikes. This can be seen in the stochastic model by projecting the solution on the 
phase plane (Fig. 5B). Upon an initial 
release from one or more 
, a large spike is generated by Ca2+-induced 
release (via the 
) during which time a decreasing 
gradually decreases the average open probability of the clustered 
. The spike is terminated when 
is too small to allow further 
release. This phenomenon is qualitatively reproduced by the deterministic model (Fig. 5D). In both the stochastic and deterministic models the decrease in average 
open probability of a cluster of 
caused by 
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 
back to a relatively high value (i.e., recovery of the 
from inhibition by 
). The deterministic model sets a clear threshold for the regeneration, as can be seen in Fig. 5C, where an upstroke in 
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 
followed by a relatively fast increase in 
(seen by combining Figs. 5C and D).
In contrast, the stochastic model enlarges the contributions of individual 
so that the generation of each spike is also effectively driven by random 
release through the 
, 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 
concentration (
), 
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 
-dependent 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.
Discussion
In this paper we address two current major questions in the field of 
modeling. Firstly, we show that 
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 
in modeling 
oscillations in ASMC by comparing a stochastic IP3R-based 
model and its associated deterministic version, for parameters such that both of the models exhibit 
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 
in the 
cluster, whereas spike initiation is significantly affected by stochastic opening of 
. Hence, repetitive 
cycling is primarily induced by the time-dependent gating variables governing transitions of the 
from one mode to another.
Our simplified two-state model of the 
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 
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 
. 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 
puffs and waves, and it remains unclear whether or not they have a similar underlying structure.
In contrast to previous deterministic ODE models, our four-variable 
model includes a more accurate 
model, as well as local control of clustered 
by two distinct 
microdomains; one at the mouth of an open 
, the other inside a cluster of 
. Neglect of either of these microdomains leads to models that either exhibit unphysiological cytoplasmic 
concentrations or fail to reproduce reasonable oscillations. This underlines the importance of taking 
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 
stochasticity does not qualitatively change the type of puff amplitude distribution except for when there are fewer than 10 
. Here, we significantly extend the scope of their study by exploring the effects of 
stochasticity on the dynamics of 
spikes, and we do this in the context of an 
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 
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 
, which explains why the one-peak type of amplitude distribution is independent of the choice of 
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 
clusters. They found that the impact of 
stochasticity on global 
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 
(which appears to be the case for ASMC, in which the wave initiation site is only 
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 
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 
-induced continuing 
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 
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 
puffs and global 
spikes (usually propagating throughout the cells in the form of traveling waves), showing that initiation of such 
spikes requires a synchronization of 
release from more than one cluster of 
[14]. This type of spiking relies on the hierarchical organization of 
signal pathways, in particular the stochastic recruitment of both individual 
and puffs at different levels [32], and therefore cannot be simply reproduced by deterministic models containing only a few ODEs. However, 
oscillations in ASMC, as observed in lung slices, may not be of this type, as IP3R-dependent puffs have not been seen in these ASMC. It thus appears that, in ASMC in lung slices, every 
“puff” initiates a wave, resulting in periodic waves with ISI that are governed by the dynamics of individual puffs.
Materials and Methods
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 
of 1.8% warm agarose in sHBSS via an intratracheal catheter. Subsequently, air (
) was injected to push the agarose within the airways into the alveoli. The agarose was polymerized by cooling to 
. A vibratome (VF-300, Precisionary Instruments, San Jose, CA) was used to make 
thick slices which were maintained in Dulbecco's Modified Eagle's Media (DMEM, Invitrogen, Carlsbad, CA) at 
in 
/air. All experiments were conducted at 
in a custom-made temperature-controlled Plexiglas chamber as described in [17].
The calcium model
Inhomogeneity of cytoplasmic 
concentration not only exists around individual channel pores of the 
, where a nearly instantaneous high 
concentration at the pore (denoted by 
) leads to a very sharp concentration profile, but is also seen inside an 
cluster where the average cluster 
concentration (
) is apparently higher than that of the surrounding cytoplasm (
) [33]. This indicates that during 
oscillations each 
is controlled by either the pore 
concentration (when it is open) or the cluster 
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 
concentration and the cluster 
concentration allows to treat 
as a parameter, providing a simpler way of modeling local 
events (like 
puffs) that has been used in several previous studies [8], [34], [35]. However, evolution of the cluster concentration and wide-field cytoplasm 
concentration are not always separable, so an additional differential equation for the cluster 
is necessary.
A schematic diagram of the model is shown in Fig. 8. The corresponding ODEs are 
(1)
(2)
(3)where 
representing total intracellular 
concentration, and thus SR 
concentration, 
is given by 
. 
and 
are the volume ratios given in Table 1. 
is the flux through the 
, 
is a background 
leak out of the SR, and 
is the uptake of 
into the SR by SERCA pumps. 
is the flux through plasma pump, and 
represents a sum of main 
influxes including 
(receptor-operated 
channel), 
(store-operated 
channel) and 
(
leak into the cell). 
coarsely models the diffusion flux from cluster microdomain to the cytoplasm. Details of the fluxes are
- • Different formulations of
give different types of models:
- a) For the stochastic model,
where 
is the maximum conductance of a cluster of 
(here 
). 
is the number of open 
determined by the states of 
.
- b) For the deterministic model we set
where 
is the 
open probability, a continuous analogue of 
.
To calculate 
and 
, we use the 
model of [7], [8], with minor modifications described later.
- •
- •
where 
and 
are obtained from [36].
- •
- •
includes a basal leak (
), receptor-operated calcium channel (ROCC, 
), store-operated calcium channel (SOCC, 
). By using the 
concentration (
) as a surrogate indicator of MCh concentration, we assume that 
. SOCC is modeled by 
[13].
- •
Calcium concentration at open channel pore (
) does not explicitly appear in the equations but is used in the 
model introduced later. 
is assumed to be proportional to SR 
concentration (
) and is therefore simply modeled by 
where 
is the value corresponding to 
. Alternatively, 
can also be assumed to be a large constant (say greater than 
) without fundamentally altering the model dynamics. The choice of 
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 
model
The 
model used in our ASMC calcium model is an improved version of the Siekmann 
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 
model which is comprised of two modes; the drive mode, containing three closed states 
, 
, 
and one open state 
, and the park mode, containing one closed state 
and one open state 
. The transition rates in each mode are constants (shown in Table 2), but 
and 
which connect the two modes are 
-/
-dependent and are formulated as 
(4)
(5)where 
, 
, 
and 
are 
-/
-modulated gating variables. 
, 
, 
and 
are either functions of 
or constants and are given later. We assume the gating variables obey the following differential equation, 
(6)where 
is the equilibrium and 
is the rate at which the equilibrium is approached. Those equilibria are functions of 
concentration at the cytoplasmic side of the 
, denoted by 
in the equations, equal to either 
or 
depending on the state of the channel). They are assumed to be 
(7)
(8)
(9)
(10)
Hence, we have stationary expressions of 
and 
, 
(11)
(12)
The expressions of 
s, 
s, 
s and 
s are chosen as follows so that Eq. 11 and Eq. 12 capture the correct trends of experimental values of 
and 
(see Fig. 9) and generate relatively smooth open probability curves (see Fig. 10), 
Note that the above formulas are different from the relatively complicated formulas used in [8]. The rates, 
, 
and 
, are constants estimated by using dynamic single channel data [4] and given in Table 2, whereas 
is not clearly revealed by experimental data. However we have shown that it should be relatively large for high 
but relatively small for low 
for reproducing experimental puff data [8]. By introducing two 
concentrations, 
and 
, 
and the state of the 
channel become highly correlated, so that we can assume 
is a relatively large value 
if the channel is open and is a relatively small value 
if the channel is closed. Hence, 
is modeled by the logic function 
Values of 
and 
are chosen so that simulated 
oscillations in ASMC are comparable to experimental observations.
The 
model reduction
Here we reduce the 6-state model to a 2-state open/closed model. The reduction takes the following steps:
- The sum of the probabilities of
, 
and 
is less than 0.03 for any 
, so they are either rarely visited by the 
or have a very short dwell time. This implies they have very little contribution to the 
dynamics. Therefore, we completely remove the three states from the full model.
- Transition rates of
and 
are about 2 orders larger than that of 
and 
, which allows us to omit the fast transitions by taking a quasi-steady state approximation. This change will affect two aspects. First, we have 
which allows us to combine 
and 
to be a new state 
, which satisfies 
. Although this means 
is a partially open state with an open probability of 
, it can be used as an fully open state in the stochastic simulations by multiplying the maximum 
flux conductance 
by a factor of 
. Secondly, 
needs to be rescaled by 
, i.e., the effective closing rate is 
.
- Due to the combination of
and 
, 
is accordingly modified to
Hence, the reduced two-state model contains one “open” state 
and one closed state 
with the opening transition rate of 
and the closing transition rate of 
.
Deterministic formulation of the stochastic model
Based on the stochastic calcium model and the reduced 2-state 
model, we construct a deterministic model. We need to modify three things that are used in the stochastic model but inapplicable to fast simulations of the deterministic model. The first is the discrete number of open channels; the second is state-dependent use of 
and 
in calculating 
and 
; the last is the logic expression of 
. Details of the modifications are as follows,
- The fraction of open channels (
) is replaced by open probability 
which is 70% of the probability of state 
.
- In the stochastic simulations,
which only controls the 
closing is primarily governed by 
, whereas 
which controls 
opening is mainly governed by 
. Therefore, in the deterministic model, we separate the functions of 
and 
by assuming 
and 
are functions of 
only whereas 
and 
are functions of 
only. That is, 
, 
, 
and 
. Here 
as defined before.
- To describe an average rate that infinitely many receptors are rapidly inhibited by high
concentration but slowly restored from 
-inhibition. 
is proposed to be
Based on the above changes, the full deterministic model containing 8 ODEs is presented as follows,
(13)
(14)
(15)
(16)
(17)
(18)
(19)
(20)where 
and 
are functions of the gating variables given by Eqs. 4 and 5. All the fluxes are the same as those of the stochastic model except 
. All the parameter values of the deterministic model are the same as those of the stochastic model and are therefore given in Tables 1 and 2.
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, 
, 
and 
are sufficiently large so that we can assume they instantaneously follow their equilibrium functions. Therefore, by taking quasi-steady state approximation to 
, 
and 
, we remove the three time-dependent variables from the full model.
By now, the full model has been reduced to a 5D model, 
(21)
(22)
(23)
(24)
(25)
Second, the rate of change of 
approaching its equilibrium, 
(calculated from Eq. 24), is at least one order larger than those of 
, 
and 
, indicating that taking the quasi-steady state approximation to Eq. 24 could not significantly affect the evolutions of 
, 
and 
. That is, 
(26)
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, 
(27)
(28)
(29)
(30)where 
is given by Eq. 26.
Supporting Information
ASMC calcium fluorescence trace data. The data files are in Excel format and compressed in a zip file. Each Excel file has a name showing their information. For example, “S2_SMC6_MCh200nM” means data are from ASMC No. 6 in lung slice No. 2 by using 200 nM MCh. In each file, there are four columns which represent (from left to right) time(s), fluorescence intensity, 
and average 
.
https://doi.org/10.1371/journal.pcbi.1003783.s001
(ZIP)
Acknowledgments
We acknowledge many useful conversations with Martin Falcke.
Author Contributions
Conceived and designed the experiments: JS GD MJS. Performed the experiments: XT. Analyzed the data: PC. Contributed to the writing of the manuscript: PC JS. Conceived and designed the model: PC JS GD. Performed the simulations: PC.
References
- 1.
De Young GW, Keizer J (1992) A single-pool inositol 1, 4, 5-trisphosphate-receptor-based model for agonist-stimulated oscillations in Ca2+ concentration. Proc Natl Acad Sci USA 89: 9895–9899.
- 2.
Dupont G, Goldbeter A (1993) One-pool model for Ca2+ oscillations involving Ca2+ and inositol 1,4,5-trisphosphate as co-agonists for Ca2+ release. Cell Calcium 14: 311–322.
- 3.
Atri A, Amundson J, Clapham D, Sneyd J (1993) A single-pool model for intracellular calcium oscillations and waves in the Xenopus laevis oocyte. Biophys J 65: 1727–1739.
- 4.
Mak DOD, Pearson JE, Loong KPC, Datta S, Fernández-Mongil M, et al. (2007) Rapid ligand-regulated gating kinetics of single IP3R Ca2+ release channels. EMBO Rep 8: 1044–1051.
- 5.
Wagner LE, Yule DI (2012) Differential regulation of the InsP3 receptor type-1 and -2 single channel properties by InsP3, Ca2+ and atp. J Physiol 590: 3245–3259.
- 6.
Ullah G, Mak DOD, Pearson JE (2012) A data-driven model of a modal gated ion channel: the inositol 1,4,5-trisphosphate receptor in insect sf9 cells. J Gen Physiol 140: 159–173.
- 7.
Siekmann I, Wagner LE, Yule DI, Crampin EJ, Sneyd J (2012) A kinetic model for IP3R type i and type ii accounting for mode changes. Biophys J 103: 658–668.
- 8.
Cao P, Donovan G, Falcke M, Sneyd J (2013) A stochastic model of calcium puffs based on single-channel data. Biophys J 105: 1133–1142.
- 9.
Smith IF, Parker I (2009) Imaging the quantal substructure of single IP3R channel activity during Ca2+ puffs in intact mammalian cells. Proc Natl Acad Sci USA 106: 6404–6409.
- 10.
Brumen M, Fajmut A, Dobovišek A, Roux E (2005) Mathematical modelling of Ca2+ oscillations in airway smooth muscle cells. J Biol Phys 31: 515–524.
- 11.
Sneyd J, Tsaneva-Atanasova K, Reznikov V, Bai Y, Sanderson MJ, et al. (2006) A method for determining the dependence of calcium oscillations on inositol trisphosphate oscillations. Proc Natl Acad Sci USA 103: 1675–1680.
- 12.
Wang IY, Bai Y, Sanderson MJ, Sneyd J (2010) A mathematical analysis of agonist- and kcl-induced Ca2+ oscillations in mouse airway smooth muscle cells. Biophys J 98: 1170–1181.
- 13.
Croisier H, Tan X, Perez-Zoghbi JF, Sandrson MJ, Sneyd J, et al. (2013) Activation of store-operated calcium entry in airway smooth muscle cells: insight from a mathematical model. PLoS ONE 8(7): e69598
- 14.
Marchant JS, Parker I (2001) Role of elementary Ca2+ puffs in generating repetitive Ca2+ oscillations. EMBO J 20: 65–76.
- 15.
Skupin A, Kettenmann H, Winkler U, Wartenberg M, Sauer H, et al. (2008) How does intracellular Ca2+ oscillate: by chance or by the clock? Biophys J 94: 2404–2411.
- 16.
Perez JF, Sanderson MJ (2005) The frequency of calcium oscillations induced by 5-ht, ach, and kcl determine the contraction of smooth muscle cells of intrapulmonary bronchioles. J Gen Physiol 125: 535–553.
- 17.
Bai Y, Edelmann M, Sanderson MJ (2009) The contribution of inositol 1,4,5-trisphosphate and ryanodine receptors to agonist-induced Ca2+ signaling of airway smooth muscle cells. Am J Physiol Lung Cell Mol Physiol 297: L347–L361.
- 18.
Bird GSJ, Putney JW (2005) Capacitative calcium entry supports calcium oscillations in human embryonic kidney cells. J Physiol 562(3): 697–706.
- 19.
Sneyd J, Tsaneva-Atanasova K, Yule DI, Thompson JL, Shuttleworth TJ (2004) Control of calcium oscillations by membrane fluxes. Proc Natl Acad Sci USA 101: 1392–1396.
- 20.
Mahn K, Hirst SJ, Ying S, Holt MR, Lavender P, et al. (2009) Diminished sarco/endoplasmic reticulum Ca2+ atpase (serca) expression contributes to airway remodelling in bronchial asthma. Proc Natl Acad Sci USA 106: 10775–10780.
- 21.
Sathish V, Leblebici F, Kip SN, Thompson MA, Pabelick CM, et al. (2008) Regulation of sarcoplasmic reticulum Ca2+ reuptake in porcine airway smooth muscle. Am J Physiol Lung Cell Mol Physiol 294: L787–L796.
- 22.
Zeller S, Rüdiger S, Engel H, Sneyd J, Warnecke G, et al. (2009) Modeling of the modulation by buffers of Ca2+ release through clusters of IP3 receptors. Biophys J 97: 992–1002.
- 23.
Li Y, Rinzel J (1994) Equations for InsP3 receptor-mediated [Ca2+]i oscillations derived from a detailed kinetic model: a hodgkin-huxley like formalism. J Theor Biol 166: 461–473.
- 24.
Keener J, Sneyd J (2009) Mathematical Physiology, Second Edition. Springer, New York.
- 25.
Ionescu L, White C, Cheung KH, Shuai J, Parker I, et al. (2007) Mode switching is the major mechanism of ligand regulation of InsP3 receptor calcium release channels. J Gen Physiol 130: 631–645.
- 26.
Shuai JW, Jung P (2002) Stochastic properties of Ca2+ release of inositol 1,4,5-trisphosphate receptor clusters. Biophys J 83: 87–97.
- 27.
Dupont G, Abou-Lovergne A, Combettes L (2008) Stochastic aspects of oscillatory Ca2+ dynamics in hepatocytes. Biophys J 95: 2193–2202.
- 28.
Kannan MS, Prakash YS, Brenner T, Mickelson JR, Sieck GC (1997) Role of ryanodine receptor channels in Ca2+ oscillations of porcine tracheal smooth muscle. Am J Physiol 272: L659–L664.
- 29.
Tazzeo T, Zhang Y, Keshavjee S, Janssen LJ (2008) Ryanodine receptors decant internal Ca2+ store in human and bovine airway smooth muscle. Eur Respir J 32: 275–284.
- 30.
Ressmeyer AR, Bai Y, Delmotte P, Uy KF, Thistlethwate P, et al. (2010) Human airway contraction and formoterol-induced relaxation is determined by Ca2+ oscillations and Ca2+ sensitivity. Am J Respir Cell Mol Biol 43: 179–191.
- 31.
Soeller C, Cannell MB (2004) Analysing cardiac excitation-contraction coupling with mathematical models of local control. Prog Biophys Mol Biol 85: 141–162.
- 32.
Thurley K, Skupin A, Thul R, Falcke M (2012) Fundamental properties of Ca2+ signals. Biochim Biophys Acta 1820(8): 1185–1194.
- 33.
Dickinson G, Parker I (2013) Factors determining the recruitment of inositol trisphosphate receptor channels during calcium puffs. Biophys J 105: 2474–2484.
- 34.
Rüdiger S, Shuai JW, Sokolov IM (2010) Law of mass action, detailed balance, and the modeling of calcium puffs. Phys Rev Lett 105(4): 048103
- 35.
Rüdiger S, Jung P, Shuai J (2012) Termination of Ca2+ release for clustered IP3R channels. PLoS Comput Biol 8(5): e1002485
- 36.
Chandrasekera PC, Kargacim ME, Deans JP, Lytton J (2009) Determination of apparent calcium affinity for endogenously expressed human sarco(endo)plasmic reticulum calcium-atpase isoform serca3. Am J Physiol Cell Physiol 296: C1105–C1114.
- 37.
Rüdiger S, Shuai JW, Huisinga W, Nagaiah C, Warnecke G, et al. (2007) Hybrid stochastic and deterministic simulations of calcium blips. Biophys J 93: 1847–1857.
concentration
influx and efflux
buffer concentration
influx and SERCA expression
oscillations