This Article Abstract Full Text (PDF) Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Falcke, M. Search for Related Content PubMed PubMed Citation Articles by Falcke, M.

On the Role of Stochastic Channel Behavior in Intracellular Ca²⁺ Dynamics

Hahn Meitner Institute, Glienicker Str. 100, 14109 Berlin, Germany and Max Planck Institute for the Physics of Complex Systems, Nöthnitzer Str. 38, 01187 Dresden, Germany

Correspondence: Address reprint requests to M. Falcke, Hahn Meitner Institute, Glienicker Str. 100, 14109 Berlin, Germany. Tel.: 49-30-80622627; Fax: 49-30-80622098; E-mail: falcke{at}hmi.de.

	ABSTRACT

TOP ABSTRACT INTRODUCTION Materials and Methods RESULTS DISCUSSION REFERENCES

 
I present a stochastic model for intracellular Ca²⁺ oscillations. The model starts from stochastic binding and dissociation of Ca²⁺ to binding sites on a single subunit of the IP₃-receptor channel but is capable of simulating large numbers of clusters for many oscillation periods too. I find oscillations with variable periods ranging from 17 s to 120 s and a standard deviation well in the experimentally observed range. Long period oscillations can be perceived as nucleation phenomenon and can be observed for a large variety of single channel dynamics. Their period depends on the geometric characteristics of the cluster array. Short periods are in the range of the time scale of channel dynamics. Both long and short period oscillations occur for parameters with a nonoscillatory deterministic regime.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION Materials and Methods RESULTS DISCUSSION REFERENCES

 
The beauty of intracellular Ca²⁺ dynamics is that it allows for observation of the build up of periodic global events from local stochastic events. A theoretical analysis on how the stochasticity of the elemental events shows up in the global events is the subject of this article.

Changes in the cytosolic free Ca²⁺ concentration are used by many cells for signaling (Tsien and Tsien, 1990; Berridge et al., 1998). The rise of that concentration is accomplished beside influx through the plasma membrane by release of Ca²⁺ from intracellular stores like the endoplasmic reticulum (ER). Opening and closing of Ca²⁺ channels on the ER membrane controls the release. Channels are closely packed into clusters (Sun et al., 1998; Thomas et al., 1998; Mak et al., 2000; Mak and Foskett, 1998). The clusters are randomly distributed on the ER membrane. Areas with high cluster density are called focal sites (Lechleiter et al., 1991; Callamaras and Parker, 2000; Marchant and Parker, 2001). Channels open and close stochastically. Stochastic behavior manifests itself as spontaneous release events through single channels or several channels in a cluster (Sun et al., 1998; Callamaras and Parker, 2000; Marchant and Parker, 2001; Thomas et al., 1999; Bootman et al., 1997a).

A channel type present in the ER membrane of many cells is the inositol 1,4,5-trisphosphate (IP₃) receptor channel (IP₃R). The open probability of the IP₃R depends on the calcium concentration on the cytosolic side of the channel and the IP₃ concentration (Taylor, 1998; Patel et al., 1999). It increases nonlinearly with the IP₃ concentration and the Ca²⁺ concentration. Hence, Ca²⁺ released by one channel increases the open probability for neighboring channels. That provides a self-amplifying release mechanism. Very high Ca²⁺ concentration inhibits the channel and terminates release.

Another element of intracellular Ca²⁺ handling are buffers. Buffers are proteins binding most of the Ca²⁺ in a cell (up to 99%). They are present in the cytosol as well as the ER. Buffers are considered as mobile or immobile depending on their diffusion characteristics. The rate constants of Ca²⁺ binding and dissociation allow for a distinction between slow and fast buffers.

The dependence of the open probability of the release channels on cytosolic Ca²⁺ creates communication between channels and allows for the formation of spatiotemporal patterns of intracellular Ca²⁺ release. These patterns show a hierarchy of phenomena. The smallest event is the opening of a single channel, called a blip. The next larger event is called puff and is the opening of several closely packed channels. Puffs can cooperate to set off a wave traveling through the cell. Waves would appear as an elevation of the Ca²⁺ concentration engulfing the whole volume in cells smaller than the wave length. If waves occur periodically, they appear as global oscillations. The type of the dominant pattern depends on the IP₃ concentration with puffs at low and waves and oscillations at high values.

Parker et al. investigated the hierarchy of spatiotemporal events in detail in Xenopus oocytes and Bootman et al. in HeLa cells (Marchant et al., 1999; Callamaras et al., 1998; Sun et al., 1998; Marchant and Parker, 2001; Thomas et al., 1999; Bootman et al., 1997a; Bootman et al., 1997b). Typically, a single puff is not enough to initiate a wave (Marchant et al., 1999). Rather, the cooperative action of several puff sites is required. The mechanism by which waves are initiated after a step increase of IP₃ is not an increase in puff amplitude. The amplitude of puffs preceding wave initiation immediately was constant. The way the critical amount of Ca²⁺ is raised is an increasing puff frequency by up to a factor of 10 (Marchant et al., 1999). That causes an elevated average Ca²⁺ concentration at the site of wave initiation before the start of the wave (Marchant et al., 1999).

The necessity for cooperativity of several puffs to initiate a wave allows for the occurrence of single puffs in between consecutive waves. Marchant et al. measured the characteristics of these puffs dependent on wave phase and period (Marchant and Parker, 2001). The measurements demonstrate typical differences between repetitive wave initiation with short and long periods. The peak of the previous wave was chosen as the reference time in the presentation of the results of the experiments (Marchant and Parker, 2001). Puffs did not occur in the first 7 s after a wave passed. Then, puffs occurred with increasing frequency. That increase continued until initiation of the next wave for short period waves (<15 s). The puff amplitude quickly reached a certain level which then stayed constant for the last 30–40% of the wave period. Waves with intermediate periods (15–50 s) exhibit an increase in puff frequency from 7 s after the previous wave until about three quarters through the cycle. The amplitude again soon reaches its steady level which it then keeps for most (60%) of the cycle. Finally, waves with periods longer than 50 s did not show any essential variation in puff amplitude or puff frequency during the 60% of the wave cycle preceding the next wave (Marchant and Parker, 2001).

A noticeable difference between short and long period waves is revealed by looking at base level Ca²⁺ concentration at the site of wave initiation (Marchant and Parker, 2001). It declines until emergence of the next wave between short period waves. Base level Ca²⁺ reaches a minimum at about half the period and then slowly rises until initiation of the next wave between long period waves. That increase in base level Ca²⁺ preceding waves was called pacemaker Ca²⁺.

A wide range of periods was found in Xenopus oocytes with the shortest ones of 10 s and the longest ones up to 2 min (Marchant and Parker, 2001). The time necessary for inhibition of the channels and recovery from it is 10 s in Xenopus oocyte (Ilyin and Parker, 1994) and is in the range of the short periods. Hence, the dynamics of the channel alone cannot account for the long period oscillations of 60 s and more.

Especially for long period waves, the nucleation of a wave by the cooperative action of a few puffs could be demonstrated (Marchant et al., 1999). Activity of a single puff site is typically not sufficient to set off a wave but it needs a supercritical nucleus formed by a few puff sites for the wave to be initiated. Because puffs are stochastic events, the formation of such a supercritical nucleus occurs with a certain probability only. Marchant and Parker (2001) state that the time elapsing between two consecutive waves is determined by two processes: The recovery from inhibition caused by the first wave and the creation of a supercritical nucleus for the second wave (Marchant and Parker, 2001). The probabilistic character of nucleation introduces variability into the wave period. Marchant et al. report a standard deviation up to 40% for long period waves (Marchant and Parker, 2001). That supports the interpretation of long period repetitive waves as nucleation phenomena.

I address the stochastic modeling of nucleation and the complete range of oscillation periods. The purpose of that effort is to check Parker's hypothesis on repetitive wave nucleation during long period oscillations from a theoretical point of view. There are interesting questions linked to that. What determines the nucleation probability? Given the fact that the periods are much longer than any process of the channel dynamics (Sneyd et al., 1995), what is the role of single channel dynamics then? The small number of channels involved into formation of a single puff makes it obvious that puffs should be random events. This is less obvious for nucleation processes involving several clusters—i.e., a few times more channels—because fluctuations decrease like with n being the number of elements involved. Will wave nucleation show stochastic behavior or resemble more a deterministic process? If it is stochastic with a nucleation probability small enough to account for long periods, would the standard deviation of the period be in the right range? Last but not least, the intracellular Ca²⁺ dynamics offers a unique opportunity to investigate the build up of global events from stochastic elemental events.

Starting point of a model for our purposes has to be of course the behavior of a single channel and its subunits as the major source of stochasticity. At the same time the model needs to be extended to length scales of several cluster spacings and simulations lasting several long periods have to be carried out. In technical terms, I intend to model a wide range of time scales (from activation of the subunit by Ca²⁺ in the milliseconds range to several periods i.e., 1000 s) and length scales (single cluster 100 nm to several cluster spacings 100 µm). Recently, we introduced a simulation concept for intracellular Ca²⁺ handling especially suited for that task. It takes into account the spatially discrete arrangement of channel clusters as well as the stochastic behavior of single channel subunits (Falcke et al., 2000b). Single subunit dynamics was represented as a stochastic realization of the DeYoung-Keizer model. I use that simulation concept as a starting point for the work presented here. The previous simulations included buffering as a constant ratio of free to buffered Ca²⁺. Here, I extend the simulations by dropping that approximation and consider state dependent buffer ratios.

I will introduce the modeling concept in the next section. The simulations show oscillation-like behavior with an average period and its standard deviation covering the complete range observed in experiments. Long period oscillations are wave nucleation processes depending most strongly on spatial characteristics of the cluster array and the number of channels per cluster. Short periods arise because fluctuations prevent the dynamics from settling in a stable stationary focus. Finally, I will show for four different cases with low and high [IP₃] that the deterministic limit of the stochastic simulation is nonoscillatory.

	Materials and Methods

TOP ABSTRACT INTRODUCTION Materials and Methods RESULTS DISCUSSION REFERENCES

 
The reaction diffusion equations
I describe the interior of the cell as a spatially two-dimensional area with a regular cluster grid with spacing d as the basic arrangement. Focal sites are mimicked by further randomly distributed clusters in certain regions of the area. The number of channels in an individual cluster was chosen with a uniform distribution between and with being the maximum number of channels per cluster.

The dynamics of the concentrations is determined by release from the endoplasmic reticulum, diffusion, buffering, and uptake of Ca²⁺ by the ER. That can be described by partial differential equations (PDEs):

(1)

The channel clusters enter the PDE for the cytosolic Ca²⁺ concentration c through the factor of the Ca²⁺ source term P_c. It depends on the spatial coordinate r: . P_c(r) vanishes outside a cluster and assumes the spatially constant positive value P inside. The number of open channels inside a cluster is represented by the size of the area where . The radius R_c of that area is determined as , where R_s is the single channel radius and N_O the number of open channels in the cluster. The single channel radius is 35 nm in most of the simulations. That is not meant to be the radius of the pore of an ion channel. Rather, it is supposed to be the size of the area where the Ca²⁺ concentration is essentially spatially homogeneous. Swillens et al. (1998) showed that it is not necessary to resolve smaller length scales for the concentration gradients occurring around a channel. Hence, the exact geometry of the Ca²⁺ source within that volume can be neglected. Furthermore, that choice of R_s describes a channel by at least three points of the numerical grid in radial direction, improving numerical accuracy. The variable E is the free Ca²⁺ concentration in the ER.

Buffers with Ca²⁺ bound in the cytosol are denoted b_i, those in the ER b_E,j. I do not consider any slow buffers for the time being but will deal with them in a later paper (Falcke, 2003). Immobile buffers are modeled by setting their diffusion coefficient equal to 0. Cytosolic buffers comprise endogenous stationary buffer b_s, endogenous mobile buffer b_m and exogenous mobile buffer b_ex. The buffers in the ER are an endogenous immobile buffer b_Es and an endogenous mobile buffer b_Em.

The components of the vector N_O(t) are the numbers of open channels of all clusters. That vector is obtained from the stochastic simulation for the channel dynamics coupled to the solution of Eq. 1 (see below). I assume that the dynamics of the different concentrations are sufficiently fast so that the concentration fields reach the stationary solution before N_O(t) changes. That allows for taking the stationary solution of Eq. 1 belonging to the current N_O(t) as the profile of the concentration field until N_O(t) changes. That approximation based on differential time scales is called adiabatic approximation. The details of how this is done are explained in the last section of this article.

The adiabatic approximation neglects the time it takes to build up the concentration profiles upon opening of a channel. That is well justified for the concentration at the location of the opening channel itself. The approximation neglects the 100–200 ms it takes for Ca²⁺ to diffuse to the neighboring cluster. One consequence is that wave speeds would be overestimated with that approach. Therefore, I do not consider them here. I will suggest a way of mapping the results of our simulations onto results of simulations with fully time dependent diffusion in the discussion section. Similarly, the slow decrease of pacemaker Ca²⁺ after closing of a channel is neglected. I will discuss that as well below.

The stochastic channel dynamics
I adopt the DeYoung-Keizer model for the subunit dynamics (DeYoung and Keizer, 1992; Keizer and DeYoung, 1994). An IP₃R consists of four identical subunits in the framework of that model. There are three binding sites on each subunit: An activating site for Ca²⁺, an inhibiting Ca²⁺ site and an IP₃ binding site. The three binding sites allow for eight different states X_ijk of each subunit. The index i stands for the IP₃ site, j for the activating Ca²⁺ site, and k for the inhibiting Ca²⁺ site. An index is 1 if an ion is bound and 0 if not. Transition probabilities per unit time of transitions involving binding of a molecule are proportional to the concentration of that molecule. Dissociation transition probabilities per unit time are constant. I assume that the channel is open if at least three of the subunits are in X₁₁₀, i.e., they have bound Ca²⁺ at the activating site and IP₃.

The experimental analog to the state "open" is actually a rapid switching between an open and a closed state of the channel with much higher open probability than other states (see e.g., (Moraru et al., 1999; Mak et al., 2001)). However, the consequence for the slow dynamics is the same as that of a current averaged in time over the rapid opening and closing because the transition probabilities depend on the time integral of the Ca²⁺ concentration (dtc(t)) only. On that basis, I use a constant current through the channel when the channel is in the open state.

According to (DeYoung and Keizer, 1992), the transition rates between the states X_0jk and X_1jk (IP₃ binding and dissociation) are two orders of magnitude faster than the other transition rates (see Table 1). Therefore, these pairs of states will reach stationary probabilities before the states of the Ca²⁺ binding sites change. Hence, we can lump the pairs X_0jk and X_1jk into single states .

View larger version (15K): [in this window] [in a new window]   FIGURE 1  Lumped states Xjk of a subunit of the IP3 receptor channel. An index is 1 if an ion is bound and 0 if not. The index j stands for the activating Ca2+ site and k for the inhibiting Ca2+ site. The transition rates are given at the edges of the rectangle.

The kinetic scheme with the lumped states is shown in Fig. 1. The state with no Ca²⁺ ion bound is X₀₀, the activated state is X₁₀ and the inhibited states are X₁₁ and X₀₁. The binding and dissociation of Ca²⁺ at the activating and inhibiting sites are stochastic events rendering the opening and closing of the channel a stochastic process. That stochastic process is coupled to the concentration of cytosolic Ca²⁺ because the binding probabilities per unit time depend on it and vice versa the number of open channels determines the concentration field.

Single cluster profiles
We have a complete description of intracellular Ca²⁺ release with the Eq. 1 and the stochastic scheme in Fig. 1 by now and can start to simulate. I show a snapshot taken from a simulation in Fig. 2. It is the concentration of fast high affinity buffer (b_ex). The total number of open channels for two simulations is shown in Fig. 3. The simulation was carried out with high fast buffer concentration. Single channel openings do not lead to higher activity for a long time but then the system switches to a state with a maintained activity of 10 open channels. The simulations shown in Figs. 2 and 3 were done by solving Eqs. 8 and 9 for each N_O(t). That is computationally expensive. A single run of 100 s real time takes several hundred hours of cpu time. That is one of the motivations to approximate the complete solution of concentration fields by the superposition of single cluster profiles. Besides that, the single cluster profiles are a convenient picture to understand what happens during the simulations.

View larger version (149K): [in this window] [in a new window]   FIGURE 2  Concentration of fast buffer with high affinity as solution of the fully coupled Eqs. 8 and 9. , , , , , , , , , , , , . The integration area is and contains 8 x 8 clusters and a focal site with on a area.

View larger version (20K): [in this window] [in a new window]   FIGURE 3  Global number of open channels No in a simulation solving the fully coupled Eq. 8 and 9 for the concentrations. [IP3] was kept at 0.06 µM until and then set to 0.24 µM. , , , , , , , , , , , , . The integration area is and contains 8 x 8 clusters.

View larger version (31K): [in this window] [in a new window]   FIGURE 4  Characteristics of single cluster profiles of free Ca2+. (A) Free Ca2+ dependent on the distance from the cluster center for a cluster with 10 open channels. (B) Concentration increase of free Ca2+ relative to the base level cb at a distance of 4.32 µm, used in some simulations as cluster spacing d, dependent on the number of open channels in a cluster. (C) Concentration increase of free Ca2+ relative to the base level cb at the center of the cluster. (D) Ratio of the contribution of the open channels (cb base level) at to the contribution at . Circles refer to calculations with and crosses to .

(2)

A single cluster profile of free cytosolic Ca²⁺ is plotted in Fig. 4 A. The concentration of free Ca²⁺ enters the channel dynamics as a factor of the transition probabilities for transitions representing Ca²⁺-binding events. Hence, rather than the absolute values of free Ca²⁺, the change in these transition probabilities per time unit relative to those at the base level concentration c_b is relevant. They are shown in Fig. 4, B–D dependent on the number of open channels in a cluster.

Puffs do not show any essential variation in amplitude when IP₃ is increased (Marchant and Parker, 2001; Sun et al., 1998). That might be due to saturation of the free Ca²⁺ peak value with increasing numbers of open channels. Such a saturation can be caused by slow diffusion in the ER. It was shown that diffusion in the ER hindered by a tubular shape can be taken into account by reducing the diffusion coefficient by about a factor of 1/2 compared to unobstructed space (ölveczky and Verkman, 1998). A further reduction of intraluminal diffusion may be due to a higher viscosity compared to the cytosol. However, high viscosity in the ER is an assumption only because I do not know any experimental indication for that. To obtain saturation of peak Ca²⁺ with an increasing number of open channels I reduced diffusion by a factor of 1/5.5 altogether (see Fig. 4, B and C). The constant puff amplitude observed in experiments may have other causes (e.g., buffer saturation). However, I have chosen to guarantee it by saturating free Ca²⁺. I show in Falcke (2003) that important results do not depend on that saturation.

Stochastic simulations comprising up to 712 clusters were performed on a cluster array like shown in Fig. 5. Clusters were arranged on a hexagonal grid with a few additional clusters scattered in between to mimic focal sites (see Fig. 5 for details).

View larger version (24K): [in this window] [in a new window]   FIGURE 5  Array of clusters used in simulations. The area is quadratic with edge length . Clusters form a hexagonal grid with distance d. NF additional clusters are scattered randomly around the points (L/4, L/4), (L/4, 3L/4), (3L/4, L/4), (3L/4, 3L/4) within a square with edge length 2.42 d mimicking a focal site; here .

	RESULTS

TOP ABSTRACT INTRODUCTION Materials and Methods RESULTS DISCUSSION REFERENCES

 
A set of simulations for different IP₃ concentrations is shown in Fig. 6. Only puffs are found at low concentrations, i.e., release events are localized and not coordinated on a length scale of several cluster spacings. That changes with the onset of global events at a little bit higher [IP₃] (Fig. 6 A). These global events are waves emerging from a nucleation area. They are very rare for low IP₃ concentrations and may travel across the whole system (Fig. 6 A, first and second peak) or fade away before they reach the system boundary (third peak). That parameter regime of abortive wave propagation is characterized by a distribution of the probability for a wave to travel a certain distance before being destroyed by fluctuations rather than steadily propagating waves (Falcke et al., 2000b).

View larger version (40K): [in this window] [in a new window]   FIGURE 6  Simulations and oscillation characteristics for different IP3 concentrations: (A) (B) (C) (D) Tav and its standard deviation dependent on [IP3].

That scenario found while going from low to high [IP₃] agrees with experimental observations (Marchant and Parker, 2001). Especially, the large range of periods of eight times the shortest one reported from the experiments is captured. Short periods could be explained by the longest time scale of the channel dynamics T_d being the transition to the inhibited state and recovery from it. However, long periods last 3–8 times longer and cannot be explained by the channel dynamics alone. That is illustrated in Fig. 7. It shows the number of open channels N_O in panel A and the number of excitable channels N_E in panel B. An excitable channel has at least three subunits in X₀₀. Waves traveling through the whole system appear as spikes in N_O or downward spikes in N_E. N_E reaches again a high level after a wave has passed within . That is far before nucleation of the next wave for at least two of the global events. I have sped up recovery from inhibition (X₀₁ X₀₀) by a factor of 2.5 in the example shown in Fig. 7 to make this difference between T_d and T_av very obvious. That difference demonstrates that the channel dynamics alone cannot account for the long average periods, because the channels are in a state ready to support another wave after time T_d.

View larger version (24K): [in this window] [in a new window]   FIGURE 7  The number of open channels (A) and the number of excitable channels (at least three subunits in X00) (B) during repetitive nucleation of waves. (A) is a blow up of panel B in Fig. 10. The transition from was sped up by a factor of 2.5 compared to the standard parameters of Table 1. See Table 1 for parameters except and .

View larger version (29K): [in this window] [in a new window]   FIGURE 10  Repetitive wave nucleation for three concentrations of exogenous buffer: (A) , (B) , (C) . The number of open channels is shown. The average period is in (A) , (B) and in (C) . The transition from was sped up by a factor of 2.5 compared to the standard parameters of Table 1. See Table 1 for further parameters except .

The nucleation probability is very small just after a wave has traveled across the system because of inhibition of most of the channels by the high Ca²⁺ concentration during the wave. That causes the deterministic part T_d of the time elapsing between two consecutive waves. T_d is determined essentially by the transition rates from the activated state (X₁₀) to the inhibited state (X₁₁) and recovery from inhibition . However, p_n is still small after recovery from inhibition. Hence, the next wave does not emerge immediately but it takes some time before another global event can be set off. I would like to illustrate that with the example of Fig. 10 B. The average period T_av of the complete run over 1000 s is 65 s. The nucleation probability is the inverse of the stochastic part T_st of T_av, i.e., . Hence, the small value of p_n provides for the larger part of T_av.

A nucleation is shown in Fig. 8. The puff activity leading to the supercritical event starts in the area marked in Fig. 8 in the frame . That area is outside but close to a focal site. The area comprises a few clusters and the interaction of these clusters allows for stochastic release to go on for a few seconds. The focal site is activated by the ongoing puff activity in its vicinity (Fig. 8, frame ). It strongly amplifies puff activity and sets off a wave. That wave propagates across the whole simulation area (Fig. 8, frames , , ). The maximum number of open channels in panel A occurs when the wave travels through the area marked in frame . Ca²⁺ release terminates first in the area where the wave started (frame ) due to inhibition. Puff activity is reduced to almost none after the wave ceased (Fig. 8, frame ).

View larger version (70K): [in this window] [in a new window]   FIGURE 8  Nucleation of a wave in a regime with long wave periods. (A) shows the number of open channels in the area where the wave nucleates. The nucleation area is bound by a black rectangle in panel (C) frame and comprises several clusters. (B) shows the global number of open channels. (C) shows six snapshots of the simulation area taken at the times indicated. Shown is Ca2+ bound to exogenous buffer. The color scaling is a minimum–maximum scaling with blue indicating minimal values and red maximal ones.

View larger version (27K): [in this window] [in a new window]   FIGURE 9  The number of excitable channels NE (at least three subunits in X00) (A). The number of open channels in an area of the same size like that marked in panel C of Fig. 8 (B). The number of open channels in the same area as in B plus the eight neighboring areas of same size (C). See Table 1 for parameters except and , .

If the nucleation of a wave is the cooperative action of a few clusters, then the nucleation probability p_n should depend very sensitively on parameters determining spatial coupling. Three of these parameters are the buffer concentrations, the cluster spacing and the number of additional clusters in focal sites. Fig. 10 compares simulations with three different concentrations of exogenous mobile buffer B_ex. An increase of B_ex by 25% extends the average period from s to . Varying the cluster spacing leads to a range of periods again from 17 s up to s (Fig. 11). Abortive waves occur for large cluster spacing. Waves disappear completely, if the spacing is increased even more. These data show that the average period in the nucleation regime does strongly depend on the strength of spatial coupling.

View larger version (48K): [in this window] [in a new window]   FIGURE 11  Repetitive wave nucleation for different cluster spacing (A) , , (B) , , (C) , . The number of open channels NO is shown. The small peaks in panel (C) are abortive waves. I included an abortive wave into the statistics as a global event, if it traveled across an area larger than 65 µm in diameter. That is the size of the area of observation used by Marchant and Parker (2001). (D) shows the period Tav and its standard deviation dependent on the cluster spacing. See Table 1 for further parameters except , (circles), , (diamonds), and , (x).

View larger version (12K): [in this window] [in a new window]   FIGURE 12  Average periods Tav for different numbers of additional clusters in focal sites.

View larger version (27K): [in this window] [in a new window]   FIGURE 13  Simulations using different random realizations of cluster grids with identical average spatial characteristics. (A) , (B) , (C) . See Table 1 for further parameters.

View larger version (42K): [in this window] [in a new window]   FIGURE 14  Simulations using a variety of parameter sets for the local dynamics. Parameters different from those in Table 1 are (top left) , , (increased deactivation), (top right) Tav for different inhibition rate constants a2 (note ) , , (bottom left) , , , , (bottom right) , , .

View larger version (18K): [in this window] [in a new window]   FIGURE 15  Simulations with different . The single channel radius was chosen proportional to . (Left) Simulation with ; (right) Tav dependent on , . See Table 1 for further parameters.

View larger version (24K): [in this window] [in a new window]   FIGURE 16  Simulations with , , and . The single channel radius was chosen proportional to as explained in the text. (A) shows the case corresponding to Fig. 6 A (thick line) and with recovery from inhibition sped up by a factor of 2.5 (thin line). (B) shows NO (thick line) and NE (thin line). The simulation is that of panel A (thick line) but perturbed at . Despite the fact that I increased the number of open channels artificially with the perturbation the consequence is a decrease of NO to almost zero followed by a spike and again stationary behavior.

View larger version (16K): [in this window] [in a new window]   FIGURE 17  Simulations with , in panel A and , in panel B, both . The single channel radius was chosen proportional to as explained in the text. (A) shows the case corresponding to Fig. 6 C. The system reaches a high activity stationary state for large times.

View larger version (15K): [in this window] [in a new window]   FIGURE 18  Comparison of the experimental data of Fig. 2 in Marchant and Parker (2001) with simulations. Shown is the relative standard deviation of Tav/Tav versus Tav.

	DISCUSSION

TOP ABSTRACT INTRODUCTION Materials and Methods RESULTS DISCUSSION REFERENCES

The adiabatic approximation implies that the concentration profile adapts immediately to the number of open channels. That means too that free Ca²⁺ goes to the base level concentration immediately upon closing of the last open channel within a vicinity of two to three cluster spacings. In other words, there is not any pacemaker Ca²⁺ if there are not any open channels around. However, the build up of pacemaker Ca²⁺ was found in experiments. To evaluate the importance of pacemaker Ca²⁺ in experiments when all channels in a wave nucleation area are closed I estimate the probability that this situation is given at an arbitrary moment assuming an average puff duration of 300 ms. That probability is 7 10^-5 for waves with periods between 15 and 50 s and 0.2 for 60-s waves. The large difference between these two values arises from the large difference of puffs per time bin preceding a wave immediately given in (Marchant and Parker, 2001), Fig. 3 C and E. Hence, the typical situation is that there are clusters with open channels in the vicinity of closed clusters. In that situation, the Ca²⁺ coming from open channels is the larger contribution to the concentration value.

The adiabatic approximation may partially compensate for the lack of pacemaker Ca²⁺ with the rapid build up of the concentration profiles. Because of the immediate appearance of an increased concentration at clusters neighbored to one which just opened, the time these channels experience an increased concentration starts earlier compared to simulations with time dependent diffusion but will end earlier too.

Pacemaker Ca²⁺ will most likely increase the opening probability when there are not any open channels in the vicinity but the last open one has closed recently. One effect of a time-dependent simulation of the diffusion and hence a more realistic development of pacemaker Ca²⁺ could be a decrease of the variability of long periods because pacemaker Ca²⁺ would support the emergence of supercritical nuclei. However, if the nucleation probability is constant during the stochastic part of the period, it fixes the standard deviation as well. Therefore, the relative standard deviation could not depend on pacemaker Ca²⁺. Pacemaker Ca²⁺ does not abolish variability of the period as the experimental data suggest. Moreover, pacemaker Ca²⁺ does not render the formation of supercritical nuclei a deterministic process as can be seen in Fig. 3 in (Bootman et al., 1997a). There, a global event does not nucleate even after pacemaker Ca²⁺ has been accumulated to a degree comparable to that one of supercritical nuclei. The role of pacemaker Ca²⁺ is the subject of ongoing simulations.

Ca²⁺ increases the opening probability per unit time p_o. For propagation, the probability that a cluster with open channels opens one channel of its neighbor clusters is relevant. As soon as one channel has opened, it dominates the Ca²⁺ concentration in its cluster. The increase of the opening probability caused by open channels at the neighboring cluster must be considerably larger than the probability of spontaneous opening to observe propagation. Local spontaneous events would dominate otherwise. If c(r,t) denotes the concentration profile caused by a cluster with open channels at and c_b the Ca²⁺ base level, is required for propagation, with l being the average life time of puffs. Hence, if we observe propagation of release, the concentration profile of free Ca²⁺ caused by puffs reaches the neighboring clusters, i.e., it is meaningful to assume single cluster profiles. Moreover, propagation is the first global phenomenon observed while going from low to high IP₃. Therefore, single cluster profiles causing an essential increase of open probabilities at neighboring clusters exist for all regimes we have modeled.

The concentration of free Ca²⁺ at a distance d from a cluster with open channels c(d,t) is of course affected by the adiabatic approximation. It monotonously increases up to a certain maximum upon opening of a cluster and then decreases again monotonously upon closing of the last open channel, if modeled with time dependent diffusion. I replace that smooth function by a stepwise function by the adiabatic approximation. The way to compare results of simulations using the adiabatic approximation with those using time dependent diffusion would be to compare the probability increase caused by a puff at its neighbor's site, i.e., the integrals only. That is of course a weaker requirement than exact agreement of the time courses. The stationary single cluster profiles can be considered typical for puffs in the sense that they generate an increase in open probability comparable to that of puffs.

I could model a range of oscillation periods from 17 s to 2 min. The variability of the period increases with its average value but stays in the experimentally observed range. These findings are in good agreement with experimental results. Long periods are determined by the nucleation probability for a wave and not by a time scale of the single channel dynamics. Hence, processes prolonging inhibition or recovery from it are not necessary to achieve long periods in the case of nucleation oscillations. The parameters the nucleation process is most sensitive to are those influencing the life time of puffs, the number of channels per cluster which can be activated and spatial coupling, i.e., the average period depends on buffer concentration too. Lukyanenko et al. found in rat ventricular cardiac myocytes that wave frequency decreases with increasing EGTA concentration (Lukyanenko and Györke 1999). That supports my results for the wave nucleation regime.

The stochastic part of the period depends on the cell size as long as the number of focal sites increases with the cell size and therefore the probability that one focal site sets off a wave. The standard deviation of the average period decreases with increasing system size in the same way. The dependence on system size does not mean that nucleation oscillations occur in large cells only. It merely means that they occur at lower IP₃ in larger cells or in cells with more focal sites. However, the standard deviation at the same T_av is smaller in large systems with a few focal sites than in small systems with only one focal site.

Dupont et al. suggested another mechanism to prolong oscillation periods (Dupont and Swillens, 1995). That mechanism relies on fast removal of free Ca²⁺ from the channel mouth so that released Ca²⁺ cannot activate other channels of the same cluster. The nucleation mechanism I propose has less stringent requirements because it requires low IP₃ concentration and weak spatial coupling only. Another possibility to prolong Ca²⁺-oscillation periods is a coupling to a phosphorylation–dephosphorylation cycle like shown in (Gall et al., 2000) or receptor phosphorylation (LeBeau et al., 1999). However, most likely both is not the case in Xenopus oocytes and cannot account for the long period oscillations there.

The onset of oscillations by wave nucleation produces the right behavior of the amplitude with increasing IP₃. Oscillations start with large amplitudes and long periods and the amplitude and period decrease with increasing IP₃ (see Marchant and Parker, 2001). Short period oscillations occur at higher IP₃ concentration values. Note that the amplitude of puffs increases only slightly with the number of activatable channels. Hence, the decrease in period is mostly due to the increase of channel number and less to higher Ca²⁺ concentration. Short periods arise under conditions when many groups of clusters are able to set off a wave, not only focal sites. Hence, they depend less on coupling parameters and the system size. The finding that oscillations depend on the number of channels per cluster for high IP₃ rather than the number of focal sites implies that the results in that [IP₃] range are valid for small cells too. Furthermore, the simulations indicate that the decrease of period with increasing IP₃ may be essentially due to increasing channel number rather than dependence of recovery rate on IP₃.

I found long period oscillations for single cluster profiles causing a 100- to 130-fold increase of binding probability for Ca²⁺ at the location of the cluster (r = 0) and approximately twofold at the neighboring cluster. That causes an increase in opening probability of the order 10⁶ and 8 respectively, assuming that all subunits are in X₀₀ and three of them need to bind Ca²⁺. The value at in terms of Ca²⁺ concentration seems rather low compared to simulations of release through a single channel with realistic current values. The computational effort necessary for the (long) simulations presented here does not allow for a comprehensive scan of the parameter space and hence I do not know, if these small values are a necessary condition for long periods to occur. In many simulations with higher concentrations long period oscillations were not observed. Short period oscillations or a high activity state could be found only. Our ongoing research addresses that question. If peak concentrations at the location of clusters as we used should turn out to be a necessary condition, it would mean that the binding probability needs to saturate at very high Ca²⁺ at values around 100 times base level or lower for nucleation oscillations to occur. Such a saturation was assumed in the model by Sneyd and Dufour (2002).

I found that periods depend on the very realization of a cluster field and not on its statistical characteristics only. That may be a simple explanation for different behavior of cells with otherwise similar characteristics. It suggests too that there are optimal configurations for a certain period T_av and T_av/T_av analogously to the finding by Shuai et al. of optimal channel numbers per cluster for quasi periodic behavior of a single cluster (Shuai and Jung, 2002). I have not performed a special search for such an optimal here because I used variable numbers of channels per cluster within a single cluster field.

The finding that the deterministic limit of systems with low and high concentration of IP₃ is a high activity stationary state is the most important result with respect to modeling concepts. That means that for some cells deterministic differential equations do not capture the right mechanism of oscillations. The oscillations in systems with small numbers of channels per cluster which can be activated by the binding of Ca²⁺ can be due to fluctuations. They may be caused by wave nucleation, fluctuations across a threshold or escape from a stationary state of a bistable system. I showed the first mechanism to be responsible for the long period oscillations of the current system.

I do not know yet the way the channel population uses to avoid the high activity stationary state at small channel numbers per cluster. Obviously, that state has a small basin of attraction. If one assumes that the dynamics can be reduced to a plane in phase space, as many deterministic models do (Tang and Othmer, 1996), a phase space structure like shown in Fig. 19 seems plausible. It shows a bistable system with the ingoing separatices of the saddle passing close by the stationary states. The separatrix close to the lower stationary state provides the excitability of the low activity stationary state. The upper stationary state is a focus. The trajectory of the system during a wave in the long period regime is shown in Fig. 19 A. Nucleation corresponds to crossing the separatrix close to the lower stationary state because of fluctuations. That is followed by motion roughly along a deterministic trajectory leading into the vicinity of the upper stationary state. The second ingoing separatrix is very close to that state at low IP₃. Hence, fluctuations push the system across it whenever it approaches the high activity stationary state. Finally, it relaxes to the lower stationary state again in a deterministic-like way.

View larger version (14K): [in this window] [in a new window]   FIGURE 19  Phase space structure which could be imagined as underlying the observed behavior of the stochastic system. Stable stationary point (filled dots); saddle point (shaded dots); outgoing separatrices (thin solid lines) ingoing separatrices (dashed lines); trajectories (thick solid lines) (A) low [IP3], (B) high [IP3].

The stationary high activity state reflects the behavior of the DeYoung-Keizer model. Such a state can be reached in the deterministic model by increasing the Ca²⁺ flux through channels which causes an increase of the cytosolic Ca²⁺ concentration. That situation occurs at the location of open channels in a discrete model. The high local Ca²⁺ concentration suppresses oscillatory behavior which needs intermediate Ca²⁺ concentrations. My findings in simulations approaching the deterministic limit suggest that the oscillatory regime in models is lost during the transition from spatially continuous channel density to spatial discreteness of the channel clusters while keeping the average flux density constant. Fluctuations compensate for that by restoring the ability for oscillatory behavior.

A high activity stationary state as one stable state of a bistable dynamics was found in a model for Xenopus oocytes with energized mitochondria too (Falcke et al., 1999; Falcke et al., 2000a). That state was stabilized by the interaction of the Ca²⁺ handling of mitochondria and the ER. The occurrence of such a state in different models of intracellular Ca²⁺ dynamics supports the idea that the cytosol might be a bistable system for high IP₃. That stationary state of the local dynamics led to the disappearance of spiral waves without abolishing less curved waves in the model with energized mitochondria because of a gap in the dispersion relation (Falcke et al., 2000a). Experimentally as well as in the model, it was found that the cytosol might stay for a rather long time in a high activity stationary state and then leave it by a front of low Ca²⁺ release or as we now know because of fluctuations (Figs. 16 B and 11 A) and Fig. 3, E and F in Falcke et al. (1999). Bistability was assumed as well in a model for the fertilization wave in Xenopus oocyte (Wagner et al., 1998).

The interplay between puffs and waves described in the introduction for experiments and as a result of the simulations is not unique to Xenopus oocytes. Puffs are observed e.g., in HeLa cells as well (Bootman et al., 1997b