Fire–diffuse–fire model of dynamics of intracellular calcium waves

Traveling waves in living cells can vary greatly in their appearance. For example, the calcium wave in immature Xenopus frog eggs propagates as a sequence of bursts (1–4) (Fig. 1a), whereas the calcium wave that occurs during fertilization in mature Xenopus eggs appears to be continuous (5, 6) (Fig. 1b). It is commonly believed that information is encoded in the time course of the Ca²⁺ signal (7–12). Thus, the distinction between these two modes of propagation likely has physiological significance. Remarkably, the waves in both cell types involve the release of calcium from internal stores via mechanisms that are physically and biologically similar.

We report here a simple physical model that allows us to characterize these behaviors as limiting cases. A prominent feature of this model is the existence of a new time scale in addition to the usual chemical time scales, namely the intersite diffusion time, d²/D. (Here, D is the diffusion coefficient and d is the source separation.) If intersite diffusion is rate-limiting, the wave will exhibit burst-like (or saltatory) behavior (Figs. 1a and 2a). If chemical processes are rate-limiting (i.e., slow compared to intersite diffusion), then the resulting wave will be smooth (Figs. 1b and 2b). Although the intersite diffusion time is likely to be important for other wave phenomena, our analysis of recent experiments suggests that it is an essential feature for intracellular Ca²⁺ waves.

Ca²⁺ is stored intracellularly in the endoplasmic or sarcoplasmic reticulum at 2–3 orders of magnitude greater than its concentration in the cytosol and is released by a process referred to as Ca²⁺-induced Ca²⁺ release (CICR). CICR involves Ca²⁺ release through specialized Ca²⁺ channels that are activated at slightly elevated levels of cytosolic Ca²⁺ and then inactivated as the level of Ca²⁺ rises further. After an open channel closes via inactivation, it cannot reopen for some time, during which it is in a “refractory” state. Calcium release provides a mechanism for wave propagation (10–17), whereas inactivation limits the amount released.

High-resolution imaging of Ca²⁺ in a variety of cell types reveals localized release events called “puffs” (1–4) (Fig. 1a), “sparks” (19), or “quantum emission domains” (8, 9). These events have been correlated with Ca²⁺ release from either individual or small clusters of channels. In muscle cells, sparks are observed at the t-tubule structures in the sarcoplasmic reticulum (20), which are aligned in regular arrays with a characteristic separation of ∼2 μm. In other cell types, inositol trisphosphate stimulates Ca²⁺ release from clusters of inositol triphosphate receptors from the endoplasmic reticulum. Although these inositol triphosphate-dependent sites appear to be randomly distributed, they are characterized by a mean separation (1–4, 21). The release events are the elementary processes underlying Ca²⁺ waves.

Our description of these events is called the “fire–diffuse–fire” model. This model reproduces the full range of wave propagation, from saltatory (Fig. 2a) to continuous (Fig. 2b). In contrast, homogeneous reaction diffusion models predict only continuous propagation with a velocity (v) proportional to \( \begin{equation*}\sqrt{D/{\tau}}\end{equation*}\), with τ the time scale of a chemical process. When applied to Ca²⁺ waves, the fire–diffuse–fire model predicts that the fertilization wave in Xenopus eggs is continuous. Moreover, it shows, in agreement with experiment, that waves in immature Xenopus oocytes (1–4) and in cardiac myocytes are saltatory.

The fire–diffuse–fire model consists of a regular array of point-source release sites with spacing, d, embedded in a continuum in which calcium ions diffuse. A related model, in which the discrete entities are individual cells, has been used for modeling cAMP (22) waves. In cells in which release sites do not occur in regular arrays, d should be thought of as a mean spacing. The pumps that resequester the ions back into the stores operate on a somewhat slower time scale (23) and are neglected here for simplicity. Whenever the cytosolic Ca²⁺ concentration in the vicinity of a release site reaches a threshold value, [Ca²⁺]_T, the site begins releasing Ca²⁺ ions at a rate σ/τ. It remains open for a time τ, and then closes after having released a quantity, σ, of Ca²⁺. Because the binding time for Ca²⁺ with immobile buffers is faster than the other time scales in the problem (15), we treat the action of the Ca²⁺ buffers as instantaneous. To account for the effect of buffering, we define σ as the total amount of Ca²⁺ released divided by a buffering factor of ∼100. We use the data in ref. 24 for the buffered diffusion coefficient in Xenopus oocytes. Note that this means that both σ and the diffusion coefficient, D, inherit the temperature dependence of the buffering factor.

The processes of Ca²⁺ release (or firing) and diffusion can be described by a single evolution equation described in detail in the Appendix and in ref. 25. The dynamics depends on two dimensionless parameters, Γ and β. The first is defined as where [Ca²⁺]_b is the basal Ca²⁺ concentration. The quantity σ/d³ is the release concentration, since it is the amount of Ca²⁺ released per site after buffering (σ) divided by the volume per site (1/d³). So Γ is the ratio of the release concentration to the difference between basal and threshold concentrations. The parameter β is the ratio of the time a site remains open, τ, to the intersite diffusion time, d²/D, i.e. If release is the rate-limiting process, then β ≫ 1, and the waves are continuous. If intersite diffusion is the rate-limiting process, then β ≪ 1, and the waves are saltatory (25).

A traveling saltatory wave solution (β = 0.08 and Γ = 3.05) to the fire–diffuse–fire model is illustrated in Fig. 2a. The vertical axis is the Ca²⁺ concentration, and the horizontal axis is position in site spacing units. The parameters were estimated from the literature on the immature oocyte (1–4, 26). The wave is traveling to the right with a velocity of 20 μm/s. Well behind the front (i.e., for x < ∼2), the concentration is [Ca²⁺](x, t) = σ/d³ + [Ca²⁺]_b = 0.66 μM. The time at which site n begins to fire (i.e., to release Ca²⁺) is denoted by t_n. The first profile in the sequence is the solid line with the short spike above site 5. This profile occurs just after the Ca²⁺ concentration at site 5 crosses threshold at time t = t₅. The second profile in the sequence is the long-dashed line with the tall peak above site 5. This is the profile at time t = t₅ + τ. The peak is wider because of diffusion and higher because the total amount (σ) of Ca²⁺ ions have been released at this time. The third profile is the short-dashed line and corresponds to a time halfway between t₅ and t₆ (after site 5 has completed firing, but before site 6 has reached threshold). The fourth profile is the solid line with the short spike above site 6. This is the profile just after the Ca²⁺ concentration reaches threshold at time t = t₆. We see that saltatory propagation consists of a distinct emission of Ca²⁺ when a particular site fires. Obviously, both Ca²⁺ release and Ca²⁺ diffusion occur simultaneously, but the firing time τ is so short compared to the intersite diffusion time, d²/D, that during the time that a single site is firing, very few ions diffuse to the next site. The shorthand description of these processes, fire–diffuse–fire, is thus an accurate description of the sequence of events shown in Fig. 2a.

Continuous propagation (β = 50 and Γ = 5) is illustrated in Fig. 2b. The parameters in this case are appropriate for the fertilization wave in mature Xenopus eggs (5). For D, we used a higher value (D = 50 μm²/s because the Ca²⁺ concentration is higher, thus implying less buffering in the egg (24). The wave is traveling to the right with a velocity of 5.2 μm/s. The experimental observations provide the concentration in the wave back, from which it is possible only to estimate σ/d³ but not σ and d separately. Here, Dτ = 450 μm². Thus, β = Dτ/d² will be in the range 28–450 for d in the range 1–4 μm. In Fig. 2b, we used d = 3 μm. These values of β are large enough so that the propagation is continuous. In this regime, the velocity does not depend strongly on d provided Γ = 5 is fixed. Decreasing d to 1 μm (while decreasing σ by a factor of 27, so that σ/d³ remains constant) results in a solution that is indistinguishable from the one in Fig. 2b, which shows the wave profile at two different times. It travels nearly without deformation.

For steady propagation, we have derived a relationship (25) of the form where ζ and η are dimensionless velocities and Combining Eqs. 2, 4, and 5, it follows that It can be shown that F ≥ 1 so that there is no solution to Eq. 3 if Γ < 1. This corresponds to the intuition that waves will not propagate if the release concentration is not sufficient to bring the neighboring sites up to threshold. In myocytes, the diffusion constant of Ca²⁺ is of the order of D = 30 μm² s⁻¹, d ≈ 2 μm, τ ≈ 15 ms, and [Ca²⁺]_T ≈ 0.2 μM (24, 28). Under conditions in which isolated sparks are found experimentally, the unitary Ca²⁺ current at a release site is estimated to be 2pA (13, 16), so that including the buffering, we find σ ≈ 1.5 × 10⁻¹⁵ μmol. Using Eqs. 1 and 2, this gives Γ ≈ 1 and β ≈ 0.1, in agreement with the observation that waves do not propagate under these conditions. We have shown for the extreme saltatory limit (τ → 0) that the fire–diffuse–fire model undergoes a series of bifurcations and crises that lead to chaos as propagation failure is approached (i.e., as Γ is decreased) (25, 28). The general feature of complex dynamics should be observed whenever such a system is driven towards propagation failure. On the other hand, waves do propagate in myocytes under conditions of Ca²⁺ overload (20), which increases the unitary current, thus increasing Γ. Simulations and more detailed calculations (27, 28) suggest that Γ ≈ 5. This generates waves with the experimentally observed velocity of about 65 μm/s, and since β ≈ 0.1, these waves are saltatory, in agreement with experiment.

Although the explicit form of F in terms of η and ζ is complicated (25), the resulting function is simple. Fig. 3 is a contour plot of F(η, ζ) upon which are superimposed straight lines ζ = \( \begin{equation*}\sqrt{{\beta}}\end{equation*}\)η for several values of β. Any two of the quantities Γ, β, η, and ζ, when different from zero, determine the other two.

If ζ²/η ≪ 1 (which holds for small enough β) the propagation is saltatory. In this case, F can be approximated by a function of η only so that Eq. 3 reduces to F(η, 0) ≡ g(η) = Γ from which we find that η = g⁻¹(Γ) so that Here, g⁻¹ denotes the inverse function of g. If Γ is large enough, we get an approximate explicit formula for g⁻¹ which gives v ≈ 4D/dlog Γ. If ζ²/η ≫ 1 and β ≫ 1, propagation is continuous. In this case, F can be approximated by a function of ζ only so that Eq. 3 reduces to F(0, ζ) ≡ f(ζ) = Γ, from which we find that ζ = f⁻¹(Γ) so that Here, if Γ is large enough, we get an approximate explicit formula for f⁻¹ that gives v ≈ \( \begin{equation*}\sqrt{{\Gamma}}\end{equation*}\)\( \begin{equation*}\sqrt{D/{\tau}}\end{equation*}\). This is Luther’s equation (15, 16) (v = α\( \begin{equation*}\sqrt{{D/{\tau}}}\end{equation*}\)) with α = \( \begin{equation*}\sqrt{{\Gamma}}\end{equation*}\). For the continuous-fertilization Ca²⁺ wave in Fig. 1b, Γ = 5, and this approximation is accurate to within a few percent. It is, however, completely incorrect for saltatory waves.

Typically, propagation is saltatory if β ≪ 1 and continuous if β ≫ 1. Given a value of Γ, it is clear from Fig. 3 that for small enough β the intersection between the line ζ = \( \begin{equation*}\sqrt{{\beta}}\end{equation*}\)η and the contour will occur at values that satisfy ζ²/η ≪ 1. This is the saltatory case. On the other hand, for large enough β, the intersection will occur at values that satisfy ζ²/η ≫ 1. This is the continuous case. One of the most distinctive physical differences between continuous and saltatory waves is the number of sites that are simultaneously firing. This is given by 1 + vτ/d = 1 + ζ²/η. The curve η = ζ² is shown as the dashed line in Fig. 3. Data points lying below this curve have only a single site firing at a time and those above have two or more. As illustrated by the vertical lines in Fig. 2b, many sites fire simultaneously during continuous propagation. On the other hand, as shown in Fig. 2a, only a single site fires at a time during saltatory propagation. More detailed models will share this and other general features of the fire–diffuse–fire model. Thus, whether a wave is saltatory or continuous will be determined by β. The trend of increasing velocity with increasing release concentration and the approach to propagation failure as the release density is decreased will both hold.

The parameter values used here for oocytes, myocytes, and eggs are within experimentally accepted ranges. However, the specific values were chosen to reproduce the experimental velocities. Other values within the experimentally observed range give velocities that are off by as much as a factor of 5. The match (or mismatch) between experimentally observed and theoretically predicted velocities cannot be used either to validate or to invalidate a particular model. For example, using the Luther equation as in refs. 15 and 16 (v = \( \begin{equation*}\sqrt{{D/{\tau}}}\end{equation*}\)) gives a velocity of 22 μm/s for the waves in the immature oocyte as in Fig. 2a. These values (D = 25 μm²/s, τ = 0.05s, and v = 22μ/s) are all within the experimental uncertainty. However, the Luther equation is manifestly wrong since the wave is saltatory (1,2). Decreasing τ by a factor of 50 in Fig. 2a results in a change in v of only a few percent, whereas the Luther equation yields a velocity of 160 μm/s! The fact that \( \begin{equation*}\sqrt{{D/{\tau}}}\end{equation*}\) is approximately equal to the velocity of the model saltatory wave is purely coincidental.

Whether the mode of propagation is continuous or saltatory is at least as important a fact as the numerical value of the velocity. Our analysis provides a simple, robust prescription for predicting the mode of propagation: estimate β. If estimates of β are all <0.1, it is safe to assume that propagation is saltatory; if they are >10, propagation is continuous. For intermediate values of β (0.1 ≤ β ≤ 10), a more accurate model (namely, ref. 25) would be necessary to determine the mode of propagation. As it is widely held that information is encoded in the Ca²⁺ signal (7–12), the distinction between saltatory and continuous waves likely has considerable physiological significance.

The fire–diffuse–fire model of calcium propagation is described by the single-evolution equation where [Ca²⁺](x, t) is the average concentration of calcium in the directions perpendicular to that of propagation (x), δ(ζ) is the δ function, H(ζ) is the Heaviside step function [H(ζ) = 0 for ζ < 0, H(ζ) = 1 for ζ ≥ 0], D is the diffusion coefficient of calcium, t_i is the first time at which the ith site takes on the threshold value, [Ca²⁺]_T, σ is the total number of ions released by a site, and τ is the “rise time” for the receptor (i.e., the length of time the receptor is open during a release event). The δ function density for the receptors means we are treating the receptor/channels as point sources. As stated in the text, the binding to buffers is assumed instantaneous. Thus, the effect of buffering is reflected in the values used for D and σ. Details can be found at http://www-xdiv.lanl.gov/XCM/pearson/fdf.