Deterministic and stochastic models of intracellular Ca2+ waves

Intracellular Ca2+ dynamics allows for the observation of wave formation from elemental release events to the global phenomenon. It shows propagating waves with deterministic features and at the same time probabilistic behaviour like variable periods. Formulation and analysis of deterministic wave models sheds new light on experimental observations and contributes to the theory of nonlinear propagating waves by exhibiting new instabilities. That is demonstrated with the example of Ca2+ waves with energized mitochondria. Stochastic models can describe spatio-temporal structures from localized puffs to propagating waves. They explain the origin of long timescales in Xenopus oocytes and wave generation. The implications of the results of both approaches for the theory of intracellular Ca2+ dynamics is discussed.


96.3
). See animation. cells smaller than the wavelength. If waves occur periodically, they appear as global oscillations. The type of the dominant pattern depends on the IP 3 concentration with puffs at low and waves and oscillations at high values.
Mitochondria are involved in energy-yielding metabolism. They contain many enzymes that together catalyze the oxidation of organic nutrients by molecular oxygen. The energy gain from oxidation is used for ATP synthesis. The transport of different chemical species involved in that process across the mitochondrial membrane and the dynamics of the membrane voltage was described in comprehensive models by Magnus et al and Selivanov et al [59,60,61,84]. In the course of these processes, mitochondria take up and release Ca 2+ and thus modulate intracellular Ca 2+ signals [9,20,22,30,34,40,41,45,48,68,77]. Especially (but not only) Friel pointed out that modulation occurs in two ways: a buffer-like action due to Ca 2+ uptake and a slowing of the decay of cytosolic transients due to mitochondrial release [15,16,34].
Most of the theoretical research has focused on deterministic models of intracellular Ca 2+ waves [3,23]- [25,29,89,90,99,100]. Only recently, the stochastic nature of intracellular Ca 2+ release has been considered [27,28,32,52,85]. This development of research reflects the character of Ca 2+ dynamics with its stochastic and deterministic features. The elemental release events can be observed as puffs occurring spontaneously and with an exponential distribution for the interval between puffs. At slightly different parameters, these puffs build up waves with constant velocity and smooth wavefront and -back (see figure 1). The single puffs can still be observed if the observation zooms in and a wave travels through the area of interest [66]. This possibility to observe the elemental events building global phenomena is one of the most fascinating features of intracellular Ca 2+ dynamics.
We discuss an example from deterministic modelling dealing with waves under conditions of energized mitochondria. That example was chosen so as to illustrate basic features of the modelling and because it explains surprising experimental results while revealing a phenomenon new to the theory of pattern formation as well: a gap in the dispersion relation. Then, we look at stochastic modelling, how long timescales can arise from stochasticity and the impact of channel numbers on wave characteristics. In the concluding section we discuss the relation between the models and current problems. 96.4

Mitochondria shape waves
Jouaville et al [48] carried out experiments investigating the role of mitochondria in pattern formation in Xenopus oocytes. Energization of mitochondria by injection of pyruvate/malate abolished spiral patterns but not waves in general. The waves forming the pattern under these conditions showed longer periods, higher amplitudes and higher velocities. Spiral waves could be restored when pump density was increased by overexpression of SERCAs [30,48].
Energization of mitochondria results in an increase of the mitochondrial effective membrane potential from −92 to about −120 mV. The mitochondrial matrix is negative relative to the surrounding cytosol. That leads to increased Ca 2+ uptake through the uniporter. Intuitively, we would expect a decrease of amplitudes and wave velocities upon energization of mitochondria. Indeed, that is what was observed in rat cortical astrocytes [22]. The wave velocity and amplitude increased upon collapsing the mitochondrial membrane potential and thus abolishing mitochondrial Ca 2+ uptake. This action of mitochondria corresponds to their behaviour analogous to buffers. However, we will see that in the experiments in Xenopus oocytes the other characteristic modulation of cytosolic Ca 2+ transients by mitochondria-the biphasic decay [16]-comes into play.
The model by Falcke et al [30] focuses on Ca 2+ uptake and release of mitochondria only to obtain a simple model explaining the impact of increased mitochondrial activity on Ca 2+ wave patterns in Xenopus oocytes. The model neglects the mitochondrial membrane potential dynamics. That can be justified by the results of Magnus and Keizer [59]- [61] who simulated membrane potential dynamics with a more detailed model and found that it changes by about 2% only under normal conditions. Furthermore, Friel et al [16] were able to fit experimental results of mitochondrial fluxes with excellent agreement to expressions for the currents neglecting membrane potential as well.
We extend the Othmer-Tang model [94] to incorporate the mechanisms of mitochondrial Ca 2+ cycling by adding a third equation governing the uptake and release of mitochondrial Ca 2+ (denoted m) and a corresponding term in the differential equation for cytosolic Ca 2+ (c) [30]. The first term of the m-dynamics is Ca 2+ uptake through the uniporter, the second one release through the Na + /Ca 2+ exchanger [30]. These terms appear again in the cytosolic dynamics. The variable c denotes the cytosolic Ca 2+ concentration. The cytosolic dynamics includes diffusion, a leak flux, the c-dependent channel flux and removal of Ca 2+ from the cytosol by pumps. The first rhs term in equation (1) models Ca 2+ diffusion in the cytosol. The second term describes Ca 2+ release from the ER by leak flux (P leak ) and channel flux (P chan ) controlled by the fraction of inhibited channels (1 − n) and the fraction of channels with activating Ca 2+ and IP 3 bound (c/(c + β 1 (1 + β 0 (I ))). The Othmer-Tang model assumes that the receptor channel has a binding site for IP 3 , an activating binding site for Ca 2+ and an inhibiting binding site for Ca 2+ . The channel opens upon binding of IP 3 and Ca 2+ to the activating site. Inhibition is provided by a second Ca 2+ binding site. If Ca 2+ is bound to the inhibiting site, the channel is closed. Binding to the inhibiting site occurs on a slower timescale and with lower affinity. Hence, an adiabatic elimination of the activating Ca 2+ and IP 3 binding processes can be applied leading to the dependence of the channel flux on Ca 2+ and IP 3 given above. The third term (P max ) models the uptake of Ca 2+ into the ER by ATPases (see [58]). The fourth and fifth terms are the contribution of mitochondria. The dynamics of the fraction of inhibited channels n is a relaxation to its Ca 2+ -dependent asymptotic value: Parameters are explained in table 1. The equation for mitochondrial dynamics essentially follows the data presented in the experimental review by Gunter and Pfeiffer [39]. Measured values of the half-maximum value K d for the uniporter from 1 to 189 µM are reported in this review. Values at the upper end of this range would not allow the uniporter to transport Ca 2+ during physiological Ca 2+ transients. Physiological evidence in many cell types shows that mitochondria respond to physiological Ca 2+ transients [2,4,5,15,17,34,36,38,42,48,75,78,80,82,88,93]. That suggests K d values in the range of concentrations occurring during cytosolic Ca 2+ transients to apply in this situation as explicitly stated in [75]. Alternatively, mitochondria might experience Ca 2+ concentrations much larger than the average bulk concentration. There is morphological evidence [21,36,41,79,83,88] showing mitochondria in close proximity to the ER, where they experience Ca 2+ concentrations considerably higher than those in bulk cytoplasm. Rizzuto et al [78,80] estimated that the uptake of Ca 2+ released from internal stores was an order of magnitude faster than that resulting from the average bulk concentration of Ca 2+ . Such a difference to bulk concentration can be included into a model by assuming the concentration at the mitochondrial uniporter always to be by a fixed factor larger than the bulk concentration. That results in re-scaling of K d by this factor. Either way leads to a small half-maximum value of K d for mitochondrial Ca 2+ uptake in the model. Similar conclusions were drawn by Sneyd et al [91]. Mitochondria in Xenopus oocytes have an average distance from release sites of 2.3 µm, which is not very close in terms of the diffusion length of free Ca 2+ in the cytosol and the argument of close proximity does not apply [67]. On the other hand, puff sites spaced less than 2.2 µm apart showed high correlation in their activity [102]. That means mitochondria are within the length scale of Ca 2+ increases at least sufficient to set off a puff. The other parameter of the uniporter dynamics, V (1) max , is varied to model the increase in membrane potential due to injection of pyruvate/malate.
The buffer-like action of mitochondria can best be observed for solitary waves since their velocity is not influenced by the dispersion relation, unlike waves in a periodic wave train. Increasing mitochondrial uptake slows down propagation of solitary waves and decreases wave amplitude (figure 2). That is analogous to the effect of buffers on wave speed [46,89,99]. Waves in astrocytes showed an increase of wave velocity and amplitude when the mitochondrial potential was collapsed (using antimycin A1 with oligomycin) corresponding to a decrease of V (1) max from some positive value to 0 in our terms [22]. This observation is compatible with the results shown in figure 2. Similarly, the observation of mitochondria limiting spatial spread of release activity found in pancreatic acinar cells would be in agreement with the buffer mode [98]. However, reduction of velocity and amplitude is not the major effect observed in Xenopus oocytes.
In Xenopus oocytes, stable spirals are observed under conditions of low mitochondrial Ca 2+ uptake [13,48,54,55]. Spiral waves are periodic patterns, unlike the solitary waves mentioned above. In simulations of this regime, a decrease in rotational frequency of spiral waves is  accompanied by a small increase in wave velocity and a small decrease in wave amplitude with increasing V (1) max ( figure 3). However, this behaviour only holds for small values of V (1) max below a critical value. When mitochondria are energized in Xenopus oocyte experiments, spiral wave patterns become unstable, disappear, and do not reform. In agreement with experiments, spirals cease to exist at a certain critical value of mitochondrial Ca 2+ uptake V (1) max,cr in simulations too. Above this value, it is found that waves emitted from pacemakers form the pattern. Examples of these waves, as observed in experiments and simulations, are shown in figure 4.
increased mitochondrial activity. The amplitudes were normalized to 1. The direction of motion is from right to left. With increased mitochondrial activity the typical biphasic decrease of Ca 2+ in the back of the pulse can be seen. In the refractory area (around 400 µm), Ca 2+ is higher with increased mitochondrial activity than under normal conditions. That prolongs receptor recovery. (B) Pulse amplitude (full curve) and velocity (dashed curve) for periodic pulse trains in dependence on the frequency of the wave train, V (1) max = 8 µM s −1 . (C) Dependence of spiral wave velocity (circles), frequency (crosses) and amplitude (triangles) on V (1) max for low mitochondrial Ca 2+ uptake (V (1) max < V (1) max,cr ) normalized to the value at V (1) max = 0. For parameters see table 1. (Figure from [30].) A simulation was also carried out corresponding to an experiment in which pyruvate/malate was injected into an oocyte [48]. In the course of this simulation, V (1) max was increased to mimic the injection. The changing Ca 2+ concentration at a location about a wavelength off the spiral core is shown in figure 4(C). The Ca 2+ oscillations associated with spiral waves cease as the system moves to a new steady state upon increasing V (1) max . This is followed by the dominance of waves emitted from a pacemaker. Figure 5 shows the transient state of a spiral wave tip when Ca 2+ uptake exceeds the critical value (V (1) max > V (1) max,cr ). When the tip bends in the early stage of spiral formation, another small-amplitude wave emerges from the back of the wave at the highest curvature (indicated by white arrow in figure 5(A)). Mitochondrial Ca 2+ efflux is responsible for this secondary wave, which in turn is responsible for prolonging the refractory state of the IP 3 R and preventing spiral formation. Although efflux plays a fundamental role in the destabilization of the spiral core, it is not the sole determinant. Here, wavefront curvature also contributes to spiral core instability. Near the spiral tip, where the wavefront curvature is the highest, Ca 2+ diffusion down the gradient of the back of the wave is focused. This focal increase in Ca 2+ further prolongs the refractory period of IP 3 Rs. Thus, both curvature and mitochondrial efflux are responsible for the generation of the secondary wave which forces the tip outward, thereby preventing spiral pattern formation ( figure 5(B)). This phenomenon was experimentally observed in the oocyte after energization as shown in figures 5(C) and (D). The free end of a Ca 2+ wave is forced outward by a secondary Ca 2+ wave and the spiral fails to form.
When periodic wave patterns of different frequencies are present in a medium, they compete for space. As time goes on, the pattern with the highest frequency generally gains spatial control of the field. It was shown experimentally that pacemakers are present in oocytes and that spiral waves dominate pacemakers in Xenopus oocytes with normal mitochondrial respiration [56]. Thus, it is only after the spirals have disappeared above V (1) max,cr -when mitochondria are energized-that the lower frequency pacemakers can govern the pattern formation in the oocyte. Hence, energization results in a wave pattern dominated by slow pacemakers. The smaller The local dynamics of the model yield three stationary states, each with different concentrations of cytosolic Ca 2+ . At low mitochondrial Ca 2+ uptake (small V (1) max ), only the state with the lowest cytosolic Ca 2+ is stable and the cytosol behaves as an excitable system (figure 6). The system becomes bistable at the uptake value max,cr ). At this point, the stationary state with the highest cytosolic Ca 2+ concentration is stabilized by increased mitochondrial Ca 2+ cycling.
When both stable stationary states exist at adjacent locations, the front connecting the states moves so that the volume occupied by one of the states grows at the expense of the other. The state which loses volume is metastable. Whether the system switches by a front from low to high cytosolic Ca 2+ or vice versa depends on the degree of mitochondrial energization. In most of the bistable region that we consider here, the state of high cytosolic Ca 2+ is metastable. Note that in the bistable regime, both pulses and fronts occur and below V (1) max,cr spirals form. Above V (1) max,cr , the region of high Ca 2+ can expand if it is surrounded by a pulse, even though it is the metastable state. Thus, a front of transition from low to high Ca 2+ can occur in this parameter range if the front immediately follows a pulse. This occurs when the unstable spiral core expands (figures 5(A)-(D)). If the pulse leading the front is extinguished by collision with another pulse, the front reverses its direction of motion and the patch of high Ca 2+ shrinks and disappears. Another way that a patch of high Ca 2+ in figure 5(B) may disappear is that a pacemaker inside it starts a front that returns the cytosol to a state of low Ca 2+ . Finally, energization of mitochondria creates a pattern in which not spirals but the waves emitted by pacemakers become the dominant structure of the bistable system.  At very high energization of the mitochondria, fronts from low to high cytosolic Ca 2+ continue to exist outside the spiral core. The region of high cytosolic Ca 2+ emerging from the spiral instability continues to expand even if the leading pulse becomes annihilated. Experimental evidence for such a transition in oocytes is shown in figure 7. Fronts in both directions (and pulses) co-exist before at still even higher V (1) max fronts switching the cytosol from high to low Ca 2+ cease to exist. If a front is initiated in this range of V (1) max , the resting Ca 2+ concentration in the oocyte is predicted to switch to the stationary state with high cytosolic Ca 2+ and wave activity stops. That is evocative of the fertilization wave in oocytes which was modelled as a bistable system too [33,100].
The mechanism of mitochondrial-induced spiral instability described above suggests that spirals could be recovered by increasing cytosolic Ca 2+ removal. Experimental studies in Xenopus oocytes show that overexpression of SERCAs permits spiral wave formation even in the presence of energized mitochondria. Consistent with experimental observations, simulations showed that an increase in SERCA density restores spiral formation at high mitochondrial Ca 2+ uptake [30].
In summary, the model consisting of equations (1) could cast the main experimental findings into a unified picture. It is worth asking whether mechanisms can be identified underlying the wave instabilities observed in simulations [31]. For the presentation of these results we change to dimensionless variables and measure space units in √ D/ and time units in 1/ . The Ca 2+ concentration is given in units of C M .
We investigate the existence and stability of a pulse with periodic boundary conditions for a range of V (1) max in the bistable regime (see figure 8). Waves are found as stationary solutions in the co-moving coordinate system. There are four branches of periodic solutions with stable and unstable parts originating from two saddle-node bifurcations at small wavelengths (see figure 8) [31]. The 'inner' pair of branches-formed by the unstable high-velocity branch and the stable low-velocity branch (see figure 8, top)-does not exist at low V (1) max . The high-velocity branches correspond to a range of high Ca 2+ travelling in a medium in the lower stationary state (see figure 9, right). The wave trains forming the outer part of spirals are on the stable branch. The unstable high-velocity pulses are distinguished from the stable ones by a delayed drop of the Ca 2+ concentration to the level of the refractory area and a correspondingly shorter refractory area. The low-velocity branch is a range of low Ca 2+ travelling in a medium in the upper stationary state (see figure 9, left).
At a certain value of V (1) max = V (1) max,tc = 16.9 µM s −1 , the stable high-velocity branch collides with the unstable one and a gap in velocity opens up (figure 8, bottom). The collision of the branches is a transcritical bifurcation at parameters (V (1) max,tc , λ tc ), where the wavelength λ of the 96.12 wave trains has to be considered as an additional bifurcation parameter. The situations above and below V (1) max,tc shown in figure 8 represent the typical unfoldings of this bifurcation with respect to the branches involved in the collision. The dispersion gap creates a range of forbidden periods. In figure 10 we show a simulation of a wave train where V (1) max is shifted from 16.5 to 17.0 µM s −1 with a period in the opening gap. Shortly after the parameter shift, the wave train disappears by flooding the refractory area with Ca 2+ . This is caused by strong Ca 2+ release from the mitochondria and is in agreement with the spiral instabilities described above.
The mechanism destroying the pulse becomes more plausible by considering the pulse profile and the unstable eigenmode ( figure 11). Mitochondria take up a great amount of Ca 2+ during the first part of the excited phase of the pulse. That mitochondrial Ca 2+ is released when cytosolic Ca 2+ decreases in the back of the pulse as can be seen in figure 11, bottom. That slows down the decrease of cytosolic Ca 2+ in the back. With increasing mitochondrial uptake (i.e. increasing V (1) max ) the amount of Ca 2+ taken up in the first part of the pulse increases and so does the amount released in the back. Finally, at the critical value of V (1) max , the large amount of released mitochondrial Ca 2+ causes a transition to the upper stationary state in the back of the pulse instead of the transition to the refractory state. This mechanism is supported by the shape of the eigenvector of the leading eigenvalue ( figure 11, top). A comparison of that eigenmode with the Goldstone mode (not shown), which would cause a simple shift of the pulse, led to the identification of the mechanism in biological terms as prolongation of recovery from inhibition. The eigenmode is identical to the Goldstone mode in the front range of the pulse. However, the unstable mode raises the level of cytosolic Ca 2+ c where the transition from the excited phase of the pulse to the refractory area occurs and hence increases the inhibitor n in the refractory area of the pulse. Free spirals at low V (1) max are stable. With slowly increasing V (1) max , spirals are destroyed at V (1) max,cr < V (1) max,tc . The occurrence of the instability in the spiral core at lower V (1) max is in agreement with the ideas about the mechanism and due to curvature effects close to the spiral core as explained above. If we stabilize the spiral core artificially, spirals are destroyed when the rim of the dispersion gap reaches their intrinsic period [31].
The dispersion gap can be perceived as a nonlinear frequency filter. Given the fact that waves on the lower branch of the dispersion relation cannot exist in two or more spatial dimensions since they shrink from the free ends [31], it would be a low pass filter. Experiments on oscillations in cells of the salivary gland of blowfly by Zimmermann et al [103] very much support this notion.
In those experiments, mitochondria were also energized by pyruvate/malate. That reduced the frequency of oscillations by a factor of 2-3 and increased their amplitude. Evidence for the reduction of frequency by mitochondria has been found in oocytes as well. Marchant et al [67] measured the puff frequency of release sites in dependence on the number and proximity of surrounding mitochondria. They found that puff sites surrounded by larger mitochondrial mass or situated less than 1.25 µm away from mitochondria exhibit puff frequencies on average 1.7 96.14 times smaller than puff sites further away from mitochondria. In each imaging frame the fastest site far away from mitochondria was on average six times faster than the slowest puff site close to mitochondria. Remarkably, the amplitude of puffs was not influenced by mitochondria [67]. These data demonstrate that mitochondria even reduce frequencies of elemental events and without being energized, if they are close enough to the release site.
Mitochondrial Ca 2+ fluxes modulate cytosolic Ca 2+ transients. One action of mitochondria is a buffer-like reduction of amplitudes of cytosolic Ca 2+ spikes and wave velocities [15,22,68]. Another characteristic effect on cytosolic Ca 2+ signals is the appearance of biphasic decay of cytosolic Ca 2+ spikes [5,16,75,79,88,93]. Mitochondrial efflux slows down the decrease of cytosolic Ca 2+ as soon as mitochondrial release becomes faster than uptake through the uniporter [16]. The effect may even lead to a plateau phase of cytosolic Ca 2+ concentration, if mitochondrial release compensates almost completely for all other fluxes removing Ca 2+ from the cytosol. The dispersion gap is the manifestation of this slowing down of the decay of cytosolic transients by mitochondria in wave characteristics, since it is caused by release of Ca 2+ in the back of a wave. Abortion of wave propagation due to the buffer features of mitochondria occurs at even higher values of mitochondrial uptake.

Stochastic simulations
The existence of spontaneous wave creation by pacemakers with long periods is crucial for the explanation of the pattern formation with energized mitochondria given above. It can be explained by stochastic models which we will consider now. However, this is not the only question motivating research. Rather, in certain cells the question for the origin of long periods in general is of interest. For instance in Xenopus oocytes, wave periods of up to 120 s are observed, but the longest timescale of the channel dynamics is only about 10-15 s. Hence, the long periods cannot be explained by the local dynamics in Xenopus oocytes. One possibility to prolong Ca 2+ oscillation periods is coupling to a phosphorylation-dephosphorylation cycle or receptor phosphorylation [35,53]. However, this is probably not the case in Xenopus oocytes. Influx across the plasma membrane does not seem to set the period either, since dynamics on the timescale of a period is not very sensitive to the Ca 2+ concentration in the extracellular medium in Xenopus oocytes [12,47].
Regardless of these considerations, the best reason for stochastic models is the observation of stochastic phenomena in intracellular Ca 2+ dynamics, the most prominent of which are the small localized release events called puffs. Puffs are the elemental events which are used to build up global events like waves and oscillations. Parker et al investigated this hierarchy of spatio-temporal structures in detail in Xenopus oocytes and Bootman et al in HeLa cells [7,8,10,65,66,92,97]. Typically, a single puff is not enough to initiate a wave [65]. Rather, the cooperative action of several puff sites is required. The mechanism by which waves are initiated following a step increase of IP 3 is not an increase in puff amplitude. The amplitude of puffs immediately preceding wave initiation was constant. The way the critical amount of Ca 2+ is raised is an increasing puff frequency by up to a factor 10 [65]. Wave initiation has been investigated for periodic waves as well [66]. 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 [66]. Puffs did not occur in the first 7 s after a wave passed. Then, puffs occurred with increasing frequency. That increase continued till initiation of the next wave for short-period waves (<15 s). The puff  (15-50 s) exhibit an increase in puff frequency from 7 s after the previous wave until about three-quarters of the way 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 [66]. Especially for long-period waves, the nucleation of a wave by the cooperative action of a few puffs can be demonstrated [65]. Since puffs are stochastic events, the formation of a supercritical nucleus occurs only with a certain probability. Marchant et al [66] 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. The probabilistic character of nucleation introduces variability into the wave period. Marchant et al [66] report a standard deviation of up to 40% for long-period waves. That supports the interpretation of long-period repetitive waves as nucleation phenomena.

Stochastic channel models
The processes causing random behaviour in intracellular Ca 2+ dynamics are the transitions between the different states of the channel subunits and the channel. Channels open and close randomly. The opening and closing probability depends on the state of the channel subunits. The opening probability is the highest, if a minimum number of subunits are activated.
The Othmer-Tang model reproduces the basic findings on IP 3 receptor channel dynamics. However, it does not consider certain details of regulation of the IP 3 R by Ca 2+ and IP 3 , as e.g. the number of Ca 2+ ions which need to bind for channel opening to occur. That might become relevant in stochastic simulations dealing with single binding and dissociation events. Hence, we use a more detailed model from now on. The model by DeYoung and Keizer was set up as a deterministic model by the authors [19] and used later on as a stochastic scheme by Falcke et al [28,32]. The model assumes identical, independent subunits. Each subunit has three binding sites: an activating binding site, an inhibiting binding site and a binding site for IP 3 . Subunits are activated, if IP 3 is bound, Ca 2+ is bound to the activating site and is not bound to the inhibiting site. The channel is activated, if at least three subunits are activated. The stochastic events are binding of Ca 2+ and IP 3 to and dissociation from the channel subunits.
We can simplify the model by taking advantage of the timescale separation between IP 3 binding and dissociation on one side and Ca 2+ binding and dissociation on the other. DeYoung and Keizer assumed the IP 3 processes to be two orders of magnitude faster than the other reaction rates. That implies that the binding state of IP 3 will be in a stationary distribution most of the time and pairs of states with identical Ca 2+ binding configuration can be lumped into one state. That leaves four lumped states of a subunit corresponding to Ca 2+ bound or not bound to the activating and inhibiting binding sites. Given one of these four states, the subunit is in one of the substates of IP 3 binding with the probability given by the stationary distribution. Such a four-state scheme is shown in figure 12. The probability for transitions corresponding to binding of Ca 2+ ions depends on the concentration of free Ca 2+ at the location of the channel.

Modelling concentrations
The Ca 2+ in a cell is transported through channels and by pumps, diffuses in the cytosol and the ER and reacts with buffers. These processes are described by reaction diffusion equations for the 96  Living cells and their subcompartments are of course spatially three-dimensional objects. Nevertheless, the reaction diffusion problem is often reduced to two spatial dimensions. The fluxes between the compartments cytosol and ER are then scaled by the volume ratio γ = V cytosol / V ER . That leads to the partial differential equations The first term in the equation for the cytosolic Ca 2+ concentration c is the diffusion term. The second and third terms model the Ca 2+ flux through the membrane of the ER. P l is the coefficient for a leak flux proportional to the concentration difference E − c. The function P c (r, t) describes the location and opening state of channel clusters. P c (r, t) has a positive value P at the location of a cluster, if channels in the cluster are open. We include the number of open channels by the size of the area where P c (r, t) is larger than zero. That area equals N o R 2 S . R S is a radius reflecting the spacing of channels within a cluster and N o is the number of open channels. The third term of the equation for c (P p ) models the action of the pumps transporting Ca 2+ from the cytosol into the ER. We describe that as a spatially continuous flux density depending on c. The last term of the equation describes the reaction of free Ca 2+ with buffers. The equation for the 96.18 dynamics of the lumenal concentration E includes diffusion, the flux through the membrane and the reaction with lumenal buffers. The flux terms are those appearing as well in the equation for cytosolic Ca 2+ only with opposite sign and scaled by the volume ratio γ . The last two equations describe the dynamics of the buffers with Ca 2+ bound. The first term is the buffer diffusion, and the second term the binding of Ca 2+ to the free buffer represented as the difference between the (constant) total concentration and the concentration of Ca 2+ -bound buffer, e.g. (B i − b i ). The last term of buffer dynamics is the rate of Ca 2+ dissociation from the buffer.
When a channel opens, a concentration profile builds up quickly in the cytosol and inside the ER. A stationary state is reached when free Ca 2+ has spread far enough so that removal of Ca 2+ from the cytosol back into the ER balances the flux through the channel.
In order to reach a simulation procedure able to simulate long timescales and large length scales adiabatic approximations were applied [27,28]. The crucial approximation was to replace the time-dependent formation of the concentration profile around a cluster by the stationary solution belonging to the actual number of open channels. This is well justified for the concentrations at the location of the cluster but meaningful on the length scale of cluster spacing too as explained in [28]. Basically, this is based on the localization of the profiles by buffers which reduce the diffusion length by binding free Ca 2+ .
Buffers localize concentration profiles so strongly that the concentration increase due to open channels is about two orders of magnitude smaller at the next neighbour's location than at the cluster with open channels itself. The localization can be used to further simplify simulations. The increase due to release through open channels is the difference between the base level concentrations (when all channels are closed) and the concentration profile (when channels are open) and we call it single-cluster profile. It approaches 0 for large distances from the cluster. In taking advantage of the localization of single-cluster profiles, the complete concentration field can be represented as sum of the base level and a superposition of all single-cluster profiles.
The single-cluster concentration profiles for all possible numbers of open channels can be calculated in advance to the simulation. During the simulation, we appoint a single-cluster profile to a cluster according to its number of open channels. The superposition of all profiles is the global concentration field [28]. The concentration fields and the cluster dynamics are mutually coupled. On the one hand, the transition rates for binding transitions depend on the concentration of free Ca 2+ as shown in figure 12. On the other hand, the concentration fields depend on the configuration of open channels.

Stochastic wave patterns
Stochastic simulations comprising up to 712 clusters were performed on a cluster array like that shown in figure 13. Clusters were arranged on a hexagonal grid with a few additional clusters scattered in between to mimic focal sites (see figure 13 for details). The number of channels in each cluster was drawn from a uniform distribution between N max K /2 and N max K . A set of simulations for different IP 3 concentrations is shown in figure 14. Only puffs are found at low concentrations, i.e. release events are localized and not coordinated on a length scale of several cluster distances. That changes with the onset of global events at a slightly higher IP 3 ( figure 14(A)). These global events are waves emerging from a nucleation area. They are very rare for low IP 3 concentrations and may travel across the whole system (figure 14(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  [28].) a wave to travel a certain distance before being destroyed by fluctuations rather than steadily propagating waves [32].
Increasing IP 3 leads to more frequent waves ( figure 14(B)) and almost every wave travels across the whole system now. The time interval in between waves is not completely regular like the period of a deterministic oscillation, but an average interwave time interval T av and its 96.20  figure 14(C)). That scenario found while going from low to high IP 3 agrees with experimental observations [66]. In particular, the large range of periods of eight times the shortest one reported from the experiments is captured. Short periods can be set essentially by the longest timescale 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.
Long periods can be explained by wave nucleation. Waves emerge from small areas-a nucleus-and than spread through the whole system. Nucleation of global events is probabilistic, because a single puff activates a neighbouring cluster with a certain probability only-not with certainty. The nucleation probability p n is small compared to the puff probability because a large supercritical nucleus of a few clusters is needed. The larger a nucleus, the smaller is the curvature of its boundary. The smaller the curvature of the boundary, the larger is the number of active neighbours of an inactive cluster just outside the nucleus and hence the probability that that cluster is activated. Hence, the larger the nucleus the larger is the probability that it grows. Deactivation, inhibition and fluctuations hinder the growth of a nucleus. In that way, a critical size of a nucleus arises.
The nucleation probability is very small just after a wave has travelled across the system because of inhibition of most of the channels by the high Ca 2+ 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 10 ) to the inhibited state (X 11 ) and recovery from inhibition (X 11 → X 01 → X 00 ). 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. The small value of p n provides for the larger part of T av at low IP 3 [28].
The period increases with decreasing strength of spatial coupling and decreases with increasing system size in agreement with the idea of nucleation processes setting long periods [27,28] (see figure 16). Average periods created by wave nucleation are in the range observed experimentally. The nontrivial information carried by the data is the good agreement between experimentally measured and simulated ratios of the standard deviation of periods to the average period T av /T av (see [28]). That supports the nucleation hypothesis. Simulations with different random realizations of a cluster array with identical average cluster spacing and average number of channels per cluster showed differential T av and a varying relative standard deviation T av /T av as well. That implies that there are optimal arrays having minimal T av /T av [28]. Oscillation periods were sensitive to changes in the number of channels per cluster and the strength of spatial coupling. In particular, the sensitivity for the average number of excitable channels per cluster raises the question for the deterministic limit. The deterministic limit assumes infinitely many channels per cluster while keeping the flux density and N max K R 2 S constant. Simulations showed that T av changes by a factor of 5 while going from 25 to 328 channels per cluster at low IP 3 . Results with even higher values of N max K suggest that the deterministic limit is not oscillatory, but is a stationary state with high activity compared to the activity between waves during long-period oscillations (figure 17). That means that the oscillations are completely due to fluctuations, i.e. stochastic binding and dissociation of Ca 2+ . Obviously, the fluctuations occurring with realistic channel numbers are large enough to leave the attractive region (stable manifold) of the high-activity stationary state, which is observed for channel numbers approaching the deterministic limit. We find again a high-activity stationary state for larger IP 3 and large N max K (figure 17, insets). That means that short-period oscillations with small N max K , too, are due to fluctuations. The stationary high-activity state reflects the behaviour of the DeYoung-Keizer model. Such a state can be reached in the deterministic model by increasing the Ca 2+ flux through channels, which causes an increase of the cytosolic Ca 2+ concentration. That situation occurs at the location of open channels in a discrete model. The high local Ca 2+ concentration suppresses oscillatory behaviour, which needs intermediate Ca 2+ concentrations. That suggests 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 behaviour.
In the simulations presented up to now, we have found a strong dependence of pattern timescales on the number of channels per cluster. Such dependence exists for other characteristics 96.22 of patterns, too. We have merely to consider larger length scales in order to be able to observe it. Patterns on these large spatial scales are shown in figures 18 and 19. The spiral in the left panel of figure 18 was simulated with a larger number of channels per cluster than the right one, but again both systems have the same mean field equations. The difference in the wavelengths of the spirals is obvious. The spiral on the left-hand side rotates with a period of about 45 s and the one on the right-hand side with approximately 75 s. The transition to the deterministic limit is demonstrated in figure 19 on the basis of spatio-temporal patterns. Pulses forming more or less stable spirals are observed at small numbers of channels per clusters. At large numbers, front waves occur. They are the trademark pattern of bistable local dynamics. The bistability disappears in a saddle node bifurcation at very low IP 3 . A very similar dependence of wavelength and periods on noise strength was found by Jung and Mayer-Kress [50] in a two-dimensional array of pulse coupled threshold elements.

Concluding remarks
Intracellular Ca 2+ dynamics exhibits a duality of stochastic and deterministic features. Puffs are clearly stochastic events. Stochasticity shows up in global events as well in the form of period distributions instead of regular oscillations and as the wave creation mechanism [27,28,32,66]. The smooth wavefronts observed in experiment (figure 1) and in stochastic simulations (figure 19) suggest that a deterministic description should be possible whereas the strong dependence of wave characteristics on channel numbers contradicts that assumption. Langevin equations like those used by Shuai et al [49,85,86] for single-cluster models are suggested by the observations in simulations. However, they fail at the small numbers of channels per cluster which occur in the wave nucleation regime.
The results of deterministic and stochastic modelling complement one another. Deterministic models allow the understanding of wave propagation with the theory of nonlinear partial differential equations and of wave instabilities with bifurcation theory. Stochastic descriptions explain local spontaneous events and the generation of waves. In the course of the research on wave instabilities with energized mitochondria, bistability was suggested by theory and confirmed experimentally by the observation of fronts figure 7. It was assumed in  the model for the fertilization wave in Xenopus oocytes as well [100]. An explanation of long timescales in Xenopus oocytes and slow pacemakers was suggested by stochastic modelling. Simulations with large numbers of channels per cluster suggest bistability as the deterministic limit of the stochastic system too (see figure 17). In both systems, pulses-a wave type more typical for excitable or oscillatory dynamic regimes-were found in a bistable regime. Pulses coexist in the deterministic model with fronts. They are due to fluctuations around the higher stationary state, which push the system out of the basin of attraction, in the stochastic model.
The most obvious drawback of the stochastic approach is that it introduces a new parameter which is the number of channels per cluster. Sun et al [92] obtained an estimate of maximal 8 open channels per cluster during a puff. That is close to the estimate by Mak et al of 5 [62] and leads to 15-30 channels per cluster taking the IP 3 binding probability into account. However, these numbers are based on single-channel current estimates which are not known very precisely.
The stochastic simulations presented here and in [27,28,32] suggest that the concept of excitable, bistable or oscillatory dynamic regimes of ordinary differential equations or partial differential equations might be of limited use applied to intracellular Ca 2+ dynamics. That is supported by the results of Shuai and Jung [49,85,86]. They found that a Hopf bifurcation of the deterministic limit of a single-cluster model does not show up in the stochastic puff behaviour. The concept of thresholds separating perturbations which decay from those which are amplified might turn out to be useful for comprehending nucleation phenomena and the transition from long-to short-period oscillations [28]. However, at the current state of research, neither the existence of oscillations nor the wave type can be reliably predicted from just the type of deterministic regime.
Nucleation could be found for a variety of parameters of the local dynamics [28]. One parameter the global behaviour of intracellular Ca 2+ dynamics is very sensitive to is the number of channels per cluster which can be activated by Ca 2+ . That number is controlled by IP 3 . In summary, the picture arises that nature created a robust stochastic mechanism that does not much care about dynamic regimes but can easily be controlled by the IP 3 concentration.
Parameters like spatial coupling and the number of channels, which can be activated by Ca 2+ , can be tuned experimentally. Hence, the experimental access as a prerequisite for fruitful research is given. Given the value and success of both approaches, the relation between stochastic and deterministic features and models of intracellular Ca 2+ dynamics will be a major issue of this research.