INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORETICAL METHODS RESULTS DISCUSSION AND SUMMARY REFERENCES

Intracellular calcium signals were first observed in medaka eggs and later on in various other cell types (Cornell-Bell et al., 1990; Bezprovanny et al., 1991). Intracellular Ca²⁺ signals are due to release of Ca²⁺ from intracellular stores such as the endoplasmic reticulum (ER) or the sarcoplasmic reticulum (SR) through inositol 1,4,5-trisphosphate receptor channels (IP₃R) or ryanodine receptor channels (RyR). Cytosolic Ca²⁺ signals in intact cells can display spatially and temporally complex patterns (Cornell-Bell et al., 1990; Newman and Zahs, 1997; Harris-White et al., 1998). They can act as second messengers in living cells to regulate multiple cellular functions such as neurotransmitter release, synaptic plasticity, gene expression, and cell death (Golovina and Blaustein, 1997; Koninck and Schulman, 1998; Allen et al., 2000).

Recently, high-resolution recordings enable us to investigate elementary intracellular Ca²⁺ release events. It has been observed that the Ca²⁺ release channels are spatially organized in clusters. The collective opening and closing of several calcium release channels in a cluster causes Ca²⁺ puffs or sparks observed in experiments (Cheng et al., 1993; Callamaras et al., 1998; Melamed-Book et al., 1999; Gonzalez et al., 2000; Mak et al., 2001). Ca²⁺ blips arising from the opening of single release channels have also been observed in experiments (Bootman et al., 1997; Lipp and Niggli, 1998; Sun et al., 1998). Typically, puffs remain spatially restricted at a low concentration of IP₃ stimulus. At high levels of [IP₃], neighboring clusters become functionally coupled by Ca²⁺ diffusion and Ca²⁺-induced Ca²⁺ release (Bezprovanny et al., 1991) to support global Ca²⁺ waves that propagate in a saltatory manner throughout the cell. Therefore, Ca²⁺ puffs serve as elementary building blocks of global Ca²⁺ waves. Moreover, puffs can arise spontaneously before a wave is initiated and can act as the triggers to initiate waves (Bootman et al., 1997). Calcium puffs provide a unique window on the dynamics of local calcium release.

Observations of Ca²⁺ signals of differing magnitudes suggested a hierarchy of calcium signaling events, with the smaller blips representing fundamental events involving opening of single IP₃R and the larger puffs being elementary events resulting from the opening of small groups of IP₃Rs (Lipp and Niggli, 1998; Bootman et al., 1997). Improved spatial and temporal resolutions in recordings reveal that Ca²⁺ release is not functionally quantized into discrete, stereotypical events of clearly separable magnitude. Instead, the amounts of calcium liberated during different events show a continuous distribution over a wide range, even when monitored from a single site (Bootman et al., 1997; Sun et al., 1998; Thomas et al., 1998; Callamaras and Parker, 2000; Marchant and Parker, 2001; Haak et al., 2001). There is not a clear distinction between fundamental and elementary events. Experimental data suggest that the localized calcium release varies in a continuous fashion due to stochastic variation in both numbers of channels recruited and durations of channel openings.

The knowledge about the calcium release mechanism is directly related to the distribution of calcium-puff amplitudes. Generally, the morphology, i.e., spatial extent, duration and amplitude, of puffs or sparks can be used to infer the release flux and the number of release channels involved locally in release (Bootman et al., 1997; Smith et al., 1998; Sun et al., 1998; Thomas et al., 1998; Jiang et al., 1999; Callamaras and Parker, 2000; Marchant and Parker, 2001). The experimental determination of the puff amplitudes is difficult at small amplitudes because the apparatus response function and cutoff can modify the actual amplitudes (Pratusevich and Balke, 1996; Izu et al., 1998; Cheng et al., 1999; Rios et al., 2001). Experimentally, it is not obvious which aspects of Ca²⁺ puffs are originally determined by the dynamics of the Ca²⁺ channels, which properties are determined by the diffusion and Ca²⁺ binding kinetics of both the intrinsic binding sites in the fiber and the Ca²⁺ indicator dye, and which properties are induced from the measurement of confocal line scan image.

It has been shown that the experimental amplitude distributions are far from the true distribution of puff amplitudes due to various factors in the experiment (Pratusevich and Balke, 1996; Smith et al., 1998; Izu et al., 1998; Jiang et al., 1999). Cheng et al. (1999) suggested that the original calcium puffs should have a monotonically decreasing amplitude distribution, regardless of whether the underlying events are stereotyped. In contrast, Rios et al. (2001) reported on either decaying amplitude distributions or distributions with a central peak.

Mathematical and computational models offer another angle to help settle these issues. Such models are directly based on the microscopic kinetics of clustered channels. The small number of calcium release channels in a cluster indicates that deterministic models might be insufficient. Thus, there is an increasing interest for the theoretical discussion on the stochastic dynamics of local Ca²⁺ release (Keizer et al., 1998; Swillens et al., 1999; Dawson et al., 1999; Moraru et al., 1999; Falcke et al., 2000; Bar et al., 2000). However, most of these studies focus on the onset of saltatory propagation of Ca²⁺ waves due to intercluster diffusion of Ca²⁺ (Keizer et al., 1998; Falcke et al., 2000). The stochastic dynamics of clustered IP₃Rs has been studied by Swillens et al. (1999) where the IP₃ receptors are assumed to be spatially distributed and coupled by calcium diffusion. The model requires 17 variables for each IP₃ receptor. They suggested that a typical cluster contained 20-30 channels in close contact to ensure efficient interchannel communication.

In this paper we expand the much simpler two-variable Li-Rinzel model (Li and Rinzel, 1994) to its Markov-stochastic version to simulate stochastic calcium release from small clusters of IP₃Rs. In the model, the channels are assumed to be close enough so that Ca²⁺ concentration can be considered homogeneous throughout the cluster. We neglect Ca²⁺ diffusion between cluster and environment without accounting for spatial aspects of the formation and collapse of localized Ca²⁺ elevations. The amplitude, lifetime, and interpuff interval distribution of calcium puffs are discussed. We show that different numbers of IP₃Rs and different IP₃ stimuli can lead to a variety of different amplitude distributions, including simply decaying distributions and single- and double-peaked distributions. Based on the open channel number distributions we inferconsistent with other independent estimates (Swillens et al., 1999)that the number of channels per cluster is around 20. We also approximate the Markov Li-Rinzel model by a Langevin-type model. It is shown that the Langevin approach is a simple but efficient approximation for the Markov process, even for a cluster with tens of IP₃Rs.

	THEORETICAL METHODS

TOP ABSTRACT INTRODUCTION THEORETICAL METHODS RESULTS DISCUSSION AND SUMMARY REFERENCES

Deterministic Li-Rinzel model

The first theoretical model for agonist-induced [Ca²⁺] oscillations based on microscopic kinetics of IP₃ and [Ca²⁺] gating of the IP₃R was proposed by De Young and Keizer (1992). The model assumes that three equivalent and independent subunits are involved in conduction in an IP₃R. Each subunit has one IP₃ binding site and two Ca²⁺ binding sites, one for activation, the other for inhibition. Thus, each subunit may exist in eight states with transitions governed by second-order and first-order rate constants. Only the state with one IP₃ and one activating Ca²⁺ bound contributes to the subunit's open probability. All three subunits must be in this state for the channel to be open. Although the model is unique in giving detailed gating kinetics, the number of variables is relatively high. It involves eight variables plus the concentration of IP₃ as a control parameter. A simplified version of the model was proposed by Li and Rinzel (1994). It is shown that the full De Young-Keizer model is symmetric in some of the binding processes and that the IP₃ binding is at least 200 times faster than the Ca²⁺ activation binding, while the Ca²⁺ activation binding is at least 10 times faster than the Ca²⁺ inactivation binding and the change rate of [Ca²⁺] during oscillations (Li and Rinzel, 1994). Considering these factors, the De Young-Keizer model can be reduced to the following system of two ordinary differential equations

(1)

(2)

with J_Channel being calcium flux from the ER to the intracellular space through the IP₃R channel, J_Pump being the calcium flux pumped from the intracellular space into the ER, and J_Leak being leakage flux from the ER to the intracellular space. The expressions for the fluxes are given by

(3)

(4)

(5)

with

(6)

The parameters of the model are c₁ = 0.185, v₁ = 6 s¹, v₂ = 0.11 s¹, v₃ = 0.9 µM s¹, k₃ = 0.1 µM, d₁ = 0.13 µM, d₂ = 1.049 µM, d₃ = 0.9434 µM, d₅ = 0.08234 µM, and a₂ = 0.2 µM¹ s¹. Conservation of Ca²⁺ implies the constraint [Ca²⁺]_ER = (c₀  [Ca²⁺])/c₁, with c₀ = 2.0 µM. The concentration [IP₃] is a control parameter.

This simplified model resembles the Hodgkin-Huxley (HH) model for electrically excitable membranes if [Ca²⁺] is replaced by the transmembrane potential. The driving force for Ca²⁺ fluxes is the concentration gradient ([Ca²⁺]  [Ca²⁺]_ER), while the driving force for the ionic currents in the HH equation is the voltage gradient.

Markov-stochastic Li-Rinzel model

Equations 1 and 2 describe the deterministic behavior averaged for a large number of channels. The small number of IP₃Rs in single clusters suggests that a stochastic formulation of these equations is necessary if calcium release from a single cluster is considered. Following the deterministic Li-Rinzel model, we only consider the stochastic opening and closing process for gate h here. Each gate h is an inactivation binding site for Ca²⁺ that is occupied (closing) or non-occupied (open). We describe the binding and unbinding of these three sites by independent two-state Markov processes with opening and closing rates _h and _h, respectively.

Thus, instead of Eq. 2, the stochastic scheme for all three gates is postulated

(7)

Equation 6 shows that the open rate _h is only determined by [IP₃] and the closing rate _h is determined by [Ca²⁺]. A single cluster consists of N IP₃Rs with three stochastic h gates each. They are globally coupled by a common but varying [Ca²⁺].

There are several ways to simulate this stochastic scheme. A widely applied approach is simply to account for the number of channels in each state of the kinetic model (Strassberg and DeFelice, 1993). The IP₃R channel can exist in four different states, and the kinetic scheme describing the behavior of this channel is given by

(8)

where [n_i] is the number of the channels with i open gates, and hence [n₃] labels the open state of the IP₃R channel. The total population of channels in each of their possible states will be tracked with time.

One can also directly simulate the stochastic dynamics of Eq. 7 for each single gate by a two-state Markov process (Jung and Shuai, 2001). This scheme can be expressed directly in terms of a computer algorithm. In detail, the state of the system is updated for every small time step t. Each IP₃R channel has three two-state h gates. If an h gate is closed at time t, then the probability that it remains closed at time t + t is exp(_h · t), and if it is open at time t, then the probability that it remains open at time t + t is exp(_h · t). To determine the state of a gate, random numbers are drawn consistent with these probabilities. Only if all three h gates in an IP₃R channel are open at time t, which means there is no Ca²⁺ bound at each of the three inactivation Ca²⁺ sites, the channel is h disinactivated or h-open. The probability that a h-open channel is conducting Ca²⁺ is given by mn. The expression for the calcium flux through the IP₃R channels replacing Eq. 3 is given by

(9)

where N and N_h-open indicate the total number of IP₃R channels and the number of h-open channels, respectively. N_h-open/N is the h-open fraction, replacing h³ in Eq. 3 of the deterministic model. Note that the h-open fraction in this method becomes a discrete variable with only N + 1 possible values, rather than a continuous variable.

In this paper we discuss the Markov dynamics of calcium release events for different IP₃R numbers N and stimuli [IP₃]. Note that the maximum flux rate c₁v₁ of the IP₃R channels is constant for varying N in our simulation. In experiments, [IP₃] can be stimulated and adjusted by the binding of an extracellular agonist such as a hormone or a neurotransmitter to receptors in the surface membrane (Callamaras et al., 1998; Sun et al., 1998).

In the model, we did not account for spatial aspects of the formation and collapse of localized Ca²⁺ elevations. The role of Ca²⁺ diffusion between the cluster and the environment in the dynamics of Ca²⁺ puff formation is neglected. However, the channels are assumed to be close enough so that the Ca²⁺ concentration within a cluster is homogeneous due to instantaneous Ca²⁺ diffusion (Swillens et al., 1999). In experiment, to study the dynamics of puffs or sparks, the clusters have to be functionally isolated. This is the case if the IP₃ concentration is low (Sun et al., 1998) and the calcium diffusion coefficient is small (Callamaras and Parker, 2000). Experimentally, weak Ca²⁺ diffusion can be achieved by intracellularly loading with the Ca²⁺ buffer EGTA (Mak and Foskett, 1997; Horne and Meyer, 1997; Thomas et al., 1998; Marchant et al., 1999; Cheng et al., 1999; Callamaras and Parker, 2000; Rios et al., 2001). With a large loading of EGTA, the clusters become functionally isolated even at large concentrations of IP₃ (Horne and Meyer, 1997; Thomas et al., 1998; Callamaras and Parker, 2000). Short-range feedback (within 0.2 µm) between individual IP₃Rs in one cluster is still intact. Thus, even with a larger IP₃ concentration, single release sites can be studied. For this experimental design our separated-cluster model with homogeneous calcium is applicable.

Langevin approach

The Markov method is conceptually simple and very accurate, as long as the random number generator is adequate and the time step t is small compared with the speed of fluctuations of the Ca²⁺ signal and channel state. However, this method is inefficient, especially for a large number of channels. It requires a large array to store the state of each gate and the generation of 3N random numbers for each time step. In the following we discuss under what conditions the Markov approach can be approximated by a Fokker-Planck equation, or equivalently by a Langevin equation for the fraction of open inactivation gates.

Because the time scale for [Ca²⁺] in the dynamic equations is the slowest, we consider the gate dynamics with constant [Ca²⁺] during each time step of iteration (0.01 s). For each gate (i = 1, 2, 3) we can write down a master equation for the numbers n_i of IP₃Rs with open gate i (Fox and Lu, 1994)

(10)

For a large number N, this master equation can be approximated by a Fokker-Planck equation, which is a linear partial differential equation for the probability of fraction of h-open gates h_i = n_i/N (Van Kampen, 1976). For every Fokker-Planck equation there is a statistically equivalent set of Langevin equations, i.e., a set of stochastic differential equations (Fox and Lu, 1994). The Langevin equation for the fraction of h-open gate h_i = n_i/N is then expressed as (Fox and Lu, 1994; Fox, 1997)

(11)

where G_hi(t) are zero mean, uncorrelated, Gaussian white-noise terms with

(12)

and i, j = 1, 2, 3. The Langevin approach indicates that the stochastic dynamics of the IP₃R cluster can be treated as a deterministic dynamics disturbed by a Gaussian white noise.

The stochastic equation for [Ca²⁺] flux through the IP₃R is given by

(13)

which replaces the [Ca²⁺] flux through the IP₃R in Eq. 3. To further simplify the problem we replace h₁h₂h₃ by h³, where h obeys Eq. 11. Thus, instead of applying independent fluctuation for each h gate of IP₃R, three identical h gates are assumed. The gain in computational speed is a factor of three. The error calculated from the mean values of [Ca²⁺], [Ca²⁺] through this approximation is <5% for N = 15 and <0.5% for N = 1000 and various [IP₃] (a comparison of [Ca²⁺] between these two approaches is given in Fig. 5 for N = 20, which is a realistic size of a cluster (Swillens et al., 1999)).

In the simulation, the Gaussian noise sources are generated at each integration step by the Box-Muller algorithm. Since h has to be bound between 0 and 1, it is necessary to verify this condition after each iteration step. The approximate nature of Eq. 11 does not automatically maintain h_i in the required interval. We simply disregard an iteration step that leads to a negative value for h. Simulation shows that the results are insensitive to the choice of strategy to keep h in [0,1].

	RESULTS

TOP ABSTRACT INTRODUCTION THEORETICAL METHODS RESULTS DISCUSSION AND SUMMARY REFERENCES

Approaching a deterministic model with large N

The C++ language is used for programming. We simulate the stochastic equations by using an explicit first-order algorithm with a time step of 0.01 s, smaller than the time constant of the h gate, e.g., _h = 1/(_h + _h) > 3 s for [Ca²⁺] <1.0 µM and [IP₃] <1.0 µM. Consistent results are obtained by changing the time step and the total simulation time.

In Fig. 1 we demonstrate that for large numbers of channels, both stochastic methods reproduce the features of the deterministic Li-Rinzel model. As [IP₃] = 0.3 µM (Fig. 1 A), the average Ca²⁺ concentration approaches the deterministic limit (dashed line) at large N. In Fig. 1 B we show that the entire bifurcation diagram is being reproduced by the stochastic model at large N. While the deterministic Li-Rinzel model predicts oscillation for 0.354 µM < [IP₃] < 0.642 µM (minima and maxima plotted), it predicts fixed points anywhere else. For N = 10⁶, this bifurcation diagram is well-reproduced by both stochastic methods (Fig. 1 B). This argument is obvious since Eq. 11 approaches Eq. 2 for N  , and the Langevin equation is the leading order approach to the exact Markov equation (Eq. 10) as the size of the fluctuations is expanded in order of 1/N. The main advantage of the Langevin approach is that the computing times do not depend on the number of channels, as they do with the Markov method.

View larger version (30K): [in this window] [in a new window]   FIGURE 1   (A) Both stochastic models approach the deterministic Li-Rinzel model at large N. In (a) we show the mean value of Ca2+ versus N obtained with the Markov method (solid square) and the Langevin approach (solid star). The deterministic limit of [Ca2+] is also shown by a dotted line. In (B) the bifurcation diagram of the stochastic Li-Rinzel model, obtained with Markov and Langevin methods, is compared with the deterministic result (solid line).

Stochastic oscillation

It is suggested that IP₃R clusters typically contain several or tens of IP₃R channels (Mak and Foskett, 1997; Bootman et al., 1997; Sun et al., 1998; Swillens et al., 1999). In the following we limit our discussion to N  100. In our simulations, different [IP₃] stimuli (i.e., [IP₃] = 0.3, 0.5, and 0.8 µM) are selected to represent the three deterministically distinguished regions (fixed point, oscillation, fixed point).

In Fig. 2 A we show results obtained with the Markov method for N = 20 and three values of [IP₃]. For [IP₃] = 0.3 µM, the deterministic dynamics gives a stable fixed point. The stochastic openings of h gates, however, initiate stochastic Ca²⁺ release, i.e., puffs. The elevated values of [Ca²⁺] in turn lead to large h-gate closing rates and termination of puffs. For [IP₃] = 0.3 µM, we have m = 0.7 and _h = 0.07. If [Ca²⁺] increases from 0.1 to 0.5 µM, m increases by 67% from 0.55 to 0.86, while _h increases by 400% from 0.02 to 0.1. This large increase of Ca²⁺-dependent inhibition (_h) causes the termination of puffs. For [IP₃] = 0.8 µM, the deterministic model gives a spiral fixed point with a pair of complex conjugate eigenvalues. For small IP₃R clusters, fluctuations will initiate large lightly damped stochastic oscillations. As a result, the Ca²⁺ trajectory spends less time at the stable fixed point than that of [IP₃] = 0.3 µM (Fig. 2 A). However, as N gets larger fluctuations become smaller, and so the stochastic Ca²⁺ oscillation will be closer to the fixed point.

View larger version (50K): [in this window] [in a new window]   FIGURE 2   (A) Stochastic oscillation of the calcium release of a cluster of 20 IP3Rs obtained with the Markov method at [IP3] = 0.3, 0.5, and 0.8 µM. The dotted line is for the result obtained with the deterministic model. (B) A calcium trace is compared with the associated trace of the fraction of h-open channels for N = 20 and [IP3] = 0.3 µM.

Most notably, the Ca²⁺ trace in the deterministically oscillatory regime can hardly be distinguished from the other traces with respect to periodicity. They are strongly dominated by the stochastic opening and closing of IP₃Rs. Thus, it is evident that the Ca²⁺ release from clusters of sizes that are believed realistic cannot be described by a deterministic equation. The Ca²⁺ signal trace appears as stochastic in the deterministically oscillating regime as it appears in the deterministically non-oscillating regimes.

The amplitudes and durations of release events (puffs) are determined by the fractions of h-open channels and their opening durations. The importance of the opening duration is substantiated in Fig. 2 B. The calcium amplitude of puff A is smaller than that of puff B, but the corresponding h-open fraction of puff A is larger than that of puff B.

In Fig. 3 A we show the Ca²⁺ signals obtained with the Langevin approach for N = 20. A plot of h-open fraction h(t) is shown in Fig. 3 B for [IP₃] = 0.30 µM and N = 20. Comparing Figs. 2 B and 3 B it can be seen that the Langevin method gives a slightly larger probability for large or small h.

View larger version (56K): [in this window] [in a new window]   FIGURE 3   (A) Stochastic oscillation of the calcium release of a cluster of 20 IP3Rs with the Langevin approach is shown for [IP3] = 0.3, 0.5, and 0.8 µM. (B) The trace of the fraction of h-open channels for N = 20 and [IP3] = 0.3 µM.

Mean value and variance of [Ca²⁺]

In this section we compare statistical properties of the calcium release using the Markov method and the Langevin approach. In Fig. 4 A the mean value [Ca²⁺] is plotted as a function of the cluster size. The Langevin approach yields results that agree qualitatively with Markov simulations even for a few tens of IP₃Rs. Simulations show that for N = 20, which is a realistic size of a cluster (Swillens et al., 1999), the results for [Ca²⁺] agree within 10% error (Fig. 5); for N = 15, the results agree within 12% error with various [IP₃].

View larger version (31K): [in this window] [in a new window]   FIGURE 4   Time average of the calcium signal as a function of the cluster size N at [IP3] = 0.30, 0.50, and 0.8 µM (A). Variance of the calcium signal as a function of the cluster size N at [IP3] = 0.30 µM (B). The results obtained with the Markov method and with the Langevin approach are marked by solid squares and stars, respectively.

View larger version (16K): [in this window] [in a new window]   FIGURE 5   Average of the calcium signal as a function of [IP3] for N = 20. The square indicates results obtained with the Markov method. The star represents results obtained with the single Langevin equation for three identical h gates. The circle represents results obtained with three Langevin equations for three different h gates.

It is interesting to note that [Ca²⁺] decreases with increasing cluster size for [IP₃] = 0.3 or 0.5 µM, while it increases for [IP₃] = 0.8 µM, as shown in Fig. 4 A. At small values of [IP₃], the baseline of [Ca²⁺] is close to the deterministically predicted value, and very small. Thus, fluctuations of the fraction of channel openings mainly yield the increase of [Ca²⁺] (see Fig. 2 A), resulting in the increase of [Ca²⁺] in comparison to the deterministic value. With large [IP₃] (e.g., 0.8 µM), the stochastic channel closings, to the contrary, can lead to a large decrease of [Ca²⁺], compared to [Ca²⁺] in the deterministic case (see Fig. 2 A). Fig. 4 A also shows that the Langevin approach gives a wrong prediction for [IP₃] = 0.8 µM with small N.

In Fig. 4 B the time averaged variance  = ([Ca²⁺]  [Ca²⁺])² is shown as a function of the cluster size N for [IP₃] = 0.3 µM. Similar results are obtained for [IP₃] = 0.5 and 0.8 µM (not shown). The variance decreases with increasing cluster size, which is also predicted by the Langevin approach with Eq. 12. In the oscillatory regime, the variance does not approach zero as N  , the same as the deterministic limit.

Amplitude distribution of puffs

Fig. 2 B shows that the number of IP₃Rs recruited in the Ca²⁺ puffs varies from puff to puff. Thus, the amplitudes and durations of the Ca²⁺ puffs vary. Important, and experimentally recorded, characteristics of these variabilities are amplitude, lifetime, and interpuff-interval distributions (Pratusevich and Balke, 1996; Bootman et al., 1997; Smith et al., 1998; Izu et al., 1998; Sun et al., 1998; Thomas et al., 1998; Cheng et al., 1999; Jiang et al., 1999; Callamaras and Parker, 2000; Marchant and Parker, 2001; Rios et al., 2001). For the analysis of these distributions of puffs, we apply a cutoff filter at a [Ca²⁺] of 0.2 µM to mimic a noise floor.

The shape of the puff amplitude distribution depends on the concentration of [IP₃] and the size of the cluster, characterized by the number of IP₃Rs in the cluster. An approximate phase diagram is shown in Fig. 6 A for N  100 and [IP₃]  1.0 µM obtained with the Markov method. If the [IP₃] stimulus is quite small, the amplitudes of the spontaneous puffs are typically smaller than 0.2 µM and are regarded as noise floor. For clusters with tens of IP₃Rs, monotonically decreasing amplitude distributions are mainly found for small [IP₃] stimulus (region II in Fig. 6 A). Single-peak amplitude distributions are mainly found for large [IP₃] (region IV in Fig. 6 A). In Fig. 6, B and C we show characteristic amplitude distributions in regions II and IV. The transition from regions II to IV is continuous in that a peak that eventually dominates for large enough [IP₃] develops additional to the decay at small amplitudes (region III in Fig. 6 A, or see Fig. 6 D).

View larger version (46K): [in this window] [in a new window]   FIGURE 6   Phase diagram of the puff amplitude distributions is shown in (A) in the [IP3]-N-plane. Samples of puff amplitude distributions obtained with the Markov method are shown in (B)-(E) in various regimes indicated in (A). Panels (F) and (G) show the amplitude distribution obtained with the Langevin approach in comparison to those obtained with the Markov method in (B) and (C).

For clusters with only a few IP₃Rs, single-peak amplitude distributions are observed in region V of Fig. 6 A; two-peak distributions are observed in region VI (see Fig. 6 E). Different from the single-peak puff amplitude distributions in region IV, the single-peak distribution for a few IP₃Rs in region V is strongly asymmetric, with a slow increase and a rapid fall-off.

The amplitude distributions obtained from the Langevin approach are also compared with the Markov method. It is shown that when N is >~15, both distributions exhibit similar shapes. As two examples, one can compare Fig. 6, B and f for N = 50, or Fig. 6, C and G for N = 20. It is also shown that the Langevin approach yields more puffs with larger amplitudes, which leads to (via the normalization) a drop-off of distribution at smaller amplitude.

Fig. 6 A shows that different distributions of puff amplitudes can be observed at different [IP₃] and N. The particular value of N does not influence the different phases very much for 10 < N < 100. The phase boundaries are mostly determined by [IP₃]. For N > 10, the Markov process can be approached by the Langevin approach, suggesting that the stochastic cluster-dynamics can be considered as a deterministic dynamics perturbed by Gaussian noise. In region II (see Fig. 6 A), where the deterministic dynamics approaches a small fixed point, we thus expect a simply decaying distribution of puff amplitude. For increasing [IP₃], the Hopf-bifurcation in the deterministic equation introduces an oscillatory component of the calcium dynamics, giving rise to a characteristic amplitude, manifesting in the peak of the puff amplitude distribution in III and the part of IV with [IP₃] < 0.64 µM. For [IP₃] > 0.64 µM, the deterministic dynamics predicts a large stationary [Ca²⁺]. Channel noise gives rise to a distribution of puff amplitude around this value.

The shapes of the puff amplitude distributions are consistent with observed amplitude distributions in experiments from Xenopus oocytes and HeLa cells (Sun et al., 1998; Thomas et al., 1998; Marchant and Parker, 2001; Haak et al., 2001).

Lifetime distribution of puffs

In addition to a wide distribution of Ca²⁺ puff amplitudes, considerable variations of lifetimes have been observed experimentally (Sun et al., 1998; Thomas et al., 1998; Haak et al., 2001). The lifetimes of puffs are measured as their full width at half-maximal amplitude (FWHM). In other words, FWHM is defined here as the time interval for which the calcium concentration profile is above one-half of the maximum concentration reached during the puff. Because the model presented does not include a spatial aspect, it is impossible to compare to FWHM in the sense of puff width (as done by Smith et al. (1998)).

Compared to various types of amplitude distribution, numerical simulation shows that the shapes of the FWHM distribution turn out to be more uniform. A typical FWHM distribution is shown in Fig. 7 A for [IP₃] = 0.3 µM and N = 20. It exhibits a single peak at ~3 s, i.e., most frequently, the duration of a puff is ~3 s. For N < 10 and IP₃  0.6 µM, a monotonically decreasing distribution is found with the Markov method. The Langevin approach can reproduce the lifetime distribution of puffs satisfactorily (see Fig. 7 B).

View larger version (24K): [in this window] [in a new window]   FIGURE 7   Lifetime distributions of puffs obtained with the Markov method (A) and the Langevin approach (B) are shown at N = 20 and [IP3] = 0.3 µM. The amplitude-lifetime scatter plot (C) demonstrates small correlations between lifetime and amplitudes of puffs.

A broad distribution of lifetimes is also observed in experiments for HeLa or oocyte cells (Sun et al., 1998; Thomas et al., 1998). A difference is that the characteristic time scale of calcium puff for HeLa or oocyte cells is ~100 ms, but the characteristic time scale for the Li-Rinzel model is ~1 s. The time scale for oscillations in the Li-Rinzel model is related to the rate of the IP₃R inactivation process. If a proper inactivation rate is set in the model, shorter lifetimes can then be observed. However, it is also possible that in HeLa and oocyte cells the buffered diffusion of intracellular Ca²⁺ can affect the lifetime of puffs, which is not addressed in the current model.

Correlation between amplitude and lifetime of puffs

Experimental observations indicate small correlations between puff amplitudes and durations (Thomas et al., 1998). Fig. 7 C shows a scatter plot of puff amplitudes versus lifetimes at [IP₃] = 0.3 µM and N = 20. Similar plots are obtained for different values of [IP₃] and N. A large puff amplitude does not correlate with a large lifetime, and vice versa (also see Fig. 2 B). Therefore, one can observe a puff with a large amplitude but short duration, or a puff with a small amplitude but long duration.

To quantitatively discuss the correlation between puff amplitude x₁ and lifetime x₂ of puffs, the correlation  is calculated with the following equation:

(14)

Simulation results show that the correlation values are typically smaller than 0.3 for various N and [IP₃].

Interpuff-interval distribution of puffs

Another important characteristics of calcium puffs is the distribution of times between two consecutive puffs. Recent experimental investigation (Marchant et al., 1999) has revealed interpuff-interval (IPI) distribution that exhibits a single peak mode. The IPIs obtained with the Markov method and Langevin approach for [IP₃] = 0.3, 0.5, and 0.8 µM are shown at N = 20 in Fig. 8, A and B, respectively. The IPIs resulting from the Langevin approach agree with those obtained with the Markov method (Fig. 8, A and B), although the lifetime distributions of puffs are not well described by the Langevin method (Fig. 7, A and B).

View larger version (24K): [in this window] [in a new window]   FIGURE 8   Interpuff-interval distributions of puffs obtained with the Markov method (A) and the Langevin approach (B) are shown at N = 20 and [IP3] = 0.3 µM (solid square), 0.5 µM (solid triangle), and 0.8 µM (open circle).

Size of IP₃R clusters

An important question is the number of IP₃Rs in a cluster. Experiments with Xenopus oocytes suggest that approximately five to eight spontaneously open IP₃Rs can be observed in a cluster (Mak and Foskett, 1997; Sun et al., 1998). Fig. 9 A shows the distributions of h-open IP₃Rs obtained from the Markov method for various cluster sizes N at [IP₃] = 0.3 µM. This value is motivated by experiments that need this [IP₃] to stimulate puffs. For the Langevin approach, the h-open fraction (i.e., h³) is a continuous number between 0 and 1. The corresponding distribution of N_h-open  (integer)(h³N) is shown in Fig. 9 B at [IP₃] = 0.3 µM for various N. It can be seen that both methods yield consistent distributions.

View larger version (33K): [in this window] [in a new window]   FIGURE 9   Distribution of h-open channel numbers for [IP3] = 0.3 µM and various cluster sizes N = 10, 15, 20, 25, 30, 40, and 50 with the Markov method (A) and the Langevin approach (B). For N = 20, e.g., the most likely number of recruited channels of a puff is 7.

Fig. 9 indicates that the larger the cluster, the more unlikely is the opening of a small number of IP₃Rs. Thus, given the typical number of five to eight spontaneously open IP₃Rs, we can conclude from Fig. 9 that the sizes of the clusters are ~1525 IP₃Rs. This result is consistent with the theoretical estimate of the cluster size by Swillens et al. (1999). They predict 20-30 IP₃Rs per cluster based on the requirement of interchannel communication.

	DISCUSSION AND SUMMARY

TOP ABSTRACT INTRODUCTION THEORETICAL METHODS RESULTS DISCUSSION AND SUMMARY REFERENCES

Previous work showed that IP₃-mediated global Ca²⁺ signals can be devolved into localized Ca²⁺ release events due to the clustered distribution of IP₃Rs (Bootman et al., 1997). Furthermore, observations of signals of differing magnitudes suggested a hierarchy of calcium signaling events from small blips to large puffs (Lipp and Niggli, 1998; Bootman et al., 1997). Improved spatial and temporal resolution in recording reveal that there is not a clear distinction between fundamental blips and elementary puffs. A continuum of Ca²⁺ release signals can be achieved due to stochastic variation in both numbers of channels recruited and durations of channel openings within a cluster (Bootman et al., 1997; Sun et al., 1998; Thomas et al., 1998; Callamaras and Parker, 2000; Marchant and Parker, 2001; Haak et al., 2001).

We have put forward a Markov version of the Li-Rinzel (1994) model to study the statistical properties of Ca²⁺ release of clusters of IP₃Rs. In comparison to the other Markov IP₃R channel models (Falcke et al., 2000; Bar et al., 2000; Swillens et al., 1999), our model is relatively simple and is represented by only two variables in the deterministic limit.

The small size of the IP₃Rs introduce stochastic oscillation into the calcium release dynamics. The stochastic oscillations are different from the stochastic excitability discussed by Keizer and Smith (1998). For the stochastic excitability, once [Ca²⁺] randomly becomes larger than a threshold, a fast release (action-potential-like) of [Ca²⁺] followed by a refractory period can be observed. For the stochastic oscillation there is not such a threshold. More broad distributions of puff amplitudes and lifetime are yielded for stochastic oscillation.

Calcium puffs vary in amplitude, lifetime, and interpuff-interval (Figs. 5-7) due to the variability of numbers of recruited channels and their open duration. This result indicates that a fixed puff morphology (i.e., amplitude and lifetime), which is sometimes assumed in literature (Pratusevich and Balke, 1996; Izu et al., 1998; Cheng et al., 1999) is not a good assumption for Ca²⁺ puff analysis.

The shape of puff amplitude distribution was a subject of a recent controversy. Although in experiment single-peak amplitude distributions are typically observed, theoretical study indicated that this feature could be caused by the failure of detecting small-amplitude puffs due to the confocal response function. Cheng et al. (1999) suggested that the original calcium puffs should exhibit an exponentially decaying amplitude distribution, regardless of whether the underlying events are stereotyped. Recently, however, Rios et al. (2001) presented new data with either decaying amplitude distributions or single-peak amplitude distributions. They suggested that if sparks were produced by individual Markovian release channels evolving reversibly, the amplitude distribution should be simply decaying. Channel groups typically give rise to a peak mode distribution in their collective spark.

Our Markov model suggests that both types of puff amplitude distributions are possible, depending on the size of the cluster and the level of [IP₃]. For tens of channels and small [IP₃] the puff amplitude distributions are typically simply decaying due to the small mean value of Ca²⁺ signals, while for a large [IP₃] the stochastic dynamics leads a single peak amplitude distribution around the large mean value.

The amplitude and duration of puffs are related to the fraction and duration of open IP₃Rs in a cluster. Roughly, a large-amplitude puff correlates with a large fraction of h-open IP₃Rs, and a long lifetime of puffs correlates with a long duration of an h-open fraction. However, a large puff may also be caused by a small fraction of IP₃Rs, but with a long duration. To typically observe five to eight spontaneously open channels, we estimate based on the Markov Li-Rinzel model that a single cluster typically includes 15-25 IP₃Rs. This result is consistent with the estimation obtained with independent methods (Swillens et al., 1999).

To shortcut the computationally expensive Markov simulations of the Ca²⁺ release of a cluster of IP₃Rs, we have introduced a Langevin-type description that is analogous to the one put forward by Fox and Lu (1994) for the Hodgkin-Huxley neuron. It is shown that, even for tens of IP₃Rs, the Langevin approximation can be used as a simple but efficient approach for the Markov process.

In this simple stochastic clustered IP₃R model, spatial aspects of the formation and collapse of localized Ca²⁺ elevations are neglected. The Ca²⁺ diffusion between the cluster and the environment is ignored so that an isolated cluster can be discussed. However, the channels are assumed to be close enough and the instantaneous Ca²⁺ diffusion within a cluster is so fast that the calcium concentration within a cluster is always homogeneous.

To observe puffs or sparks in the experiment, calcium diffusion is suppressed by intracellular loading with the Ca²⁺ buffer EGTA (Mak and Foskett, 1997; Horne and Meyer, 1997; Thomas et al., 1998; Marchant et al., 1999; Cheng et al., 1999; Callamaras and Parker, 2000; Rios et al., 2001). With a large loading of EGTA, the clusters become functionally isolated (Callamaras and Parker, 2000). Under these condition, our model is valid.

Experimental and theoretical work (Roberts, 1993; Bertram et al., 1999) suggests that even at steady state the Ca²⁺ diffusion at a Ca²⁺ release site may lead to inhomogeneous profiles, suggesting that the diffusion within a cluster may affect the puff dynamics. However, a more realistic model put forward by Swillens et al. (1999) shows that the simplification applied in our paper affects the main results only insignificantly. Swillens et al. (1999) considered a stochastic clustered IP₃R model within a 3-D Cartesian space. Each channel occupies a certain position inside the cluster on the ER membrane and is in contact with the Ca²⁺ concentration in the adjacent small cubic volume. Their simulation result suggested that a typical cluster could contain 20-30 channels. Their data also showed that the distribution of puff amplitudes is monotonically decreasing for a small [IP₃] and N = 25, and has a single-peak mode for large [IP₃]. Based on the discussion of a simplified model, they indicated that the kinetic behavior of a cluster could be satisfactorily simulated by considering a virtual domain in which all of the channels of a cluster and the calcium concentration were homogeneously distributed.

Three distinct physiological mechanisms have been proposed to underlie Ca²⁺ release and uptake in cells: excitable, oscillatory, and bistable states. All three physiological states can be derived from a Li-Rinzel model with different parameters (Keizer and Smith, 1998). It is interesting to discuss how the stochastic properties of dynamics of Ca²⁺ release depend on the different dynamics states, which is the subject of our current research and will be discussed in a forthcoming paper.

The stochastic oscillation for small [IP₃] facilitates Ca²⁺ signals that may regulate other cell functions. One could hypothesize that the small cluster size serves exactly this purpose. In a recent paper (Shuai and Jung, 2002), this hypothesis is supported by the coherence analysis of Ca²⁺ signals released from a single cluster. The behavior of coherence resonance (Pikovsky and Kurths, 1997) is found for [IP₃] < 0.354 µM, suggesting that the regularity of calcium signaling can be optimized at a certain cluster size.