INTRODUCTION

Top Abstract Introduction Discussion Appendix References

Fluorescence imaging of Ca²⁺ in living cells has revealed localized events referred to variously as "puffs" (Parker and Yao, 1991), "quantum emission domains" (Llinas et al., 1992), "sparks" (Cheng et al., 1993), and "elementary calcium-release units" (Horne and Meyer, 1997). These events are associated with Ca²⁺ flux into the cytosol through individual or small clusters of Ca²⁺ channels (Berridge, 1997). Ca²⁺ sparks, first characterized in cardiac myocytes (Cheng et al., 1993), also have been seen in skeletal (Schneider and Klein, 1996) and smooth muscle (Nelson et al., 1995). In myocytes sparks are associated with t-tubule structures and ryanodine receptor (RyR) Ca²⁺ channels in the sarcoplasmic reticulum (SR) (Shacklock et al., 1995; Parker et al., 1996). Ca²⁺ sparks are essential unitary events in excitation-contraction coupling (Cannell et al., 1995), and coronary defects in rats have been shown to correlate with a decreased occurrence of sparks (Gomez et al., 1997).

In myocytes sparks originate from submicron-sized sites, have a spatial extent of several microns, and a peak Ca²⁺ concentration and duration of ~0.3 µM and 100 msec, respectively (Cannell et al., 1995). In low external Ca²⁺ sparks are isolated random events, but after external Ca²⁺ is increased, sparks can serve as sites for initiation and propagation of Ca²⁺ waves. The saltatory nature of these waves (Cheng et al., 1996) and their speed (60-80 µm s¹) suggest that their initiation and propagation is different from other cytosolic Ca²⁺ waves, which can be described by continuous reaction-diffusion equations (Murray, 1989; Jaffe, 1993; Atri et al., 1993; Jafri and Keizer, 1995).

Here we use computer simulations to investigate how a regular array of release sites influences the propagation of Ca²⁺ waves in cardiac myocytes. We introduce a kinetic model of a release site that generalizes an earlier model of adaptation of the ryanodine receptor (Keizer and Levine, 1996) and that mimics the behavior of isolated sparks observed in cardiac myocytes. Simulations with equally spaced release sites in one spatial dimension lead to saltatory propagation of Ca²⁺ waves. We find that the saltatory wave speed is proportional to the diffusion constant of calcium, rather than its square root, as would be expected for a continuum wave. By using a simplified caricature of release sites coupled via Ca²⁺ diffusion (the "fire-diffuse-fire" model), we explore the nature of the saltatory wave. Analysis of the fire-diffuse-fire model defines the parameter range for successful wave propagation and gives a simple criterion for distinguishing saltatory and continuous propagation modes. Suggestions are made for how to distinguish saltatory and continuous propagation experimentally.

	SIMULATION OF SPARKS IN MYOCYTES

We have carried out computer simulations of spark-induced waves to explore the influence of the regular array of release sites on their propagation. The simulations, which are described in the Appendix, combine Ca²⁺ diffusion with a simple kinetic model of the release site (J_site) and a Ca²⁺ leak (J_leak) and re-uptake into the SR via SERCA pumps (J_serca). The model release site reproduces the average rise and refractory times of a spark and includes adaptive behavior that mimics measurements on isolated RyRs in bilayers (Györke and Fill, 1993; Keizer and Levine, 1996). The simulations are deterministic, rather than stochastic, as described elsewhere (Keizer and Smith, 1998), since here we focus on wave propagation rather than initiation and termination. Fig. 1 shows a space-time (or "waterfall") plot for a typical wave (speed v = 67 µm s¹) initiated by a localized increase in the myoplasmic Ca²⁺ concentration, [Ca_i²⁺], around x = 0 that simulates the opening of several release sites. Ca²⁺ diffuses in both directions, triggering release of additional Ca²⁺ from neighboring sites (separated by d = 2.0 µm) via Ca²⁺-induced Ca²⁺ release (CICR). The wave is composed of a regular sequence of sparks, evident as regions of elevated Ca²⁺ that last for ~120 ms. The shape and duration of the sparks and the wave speed, v, are comparable to that found in cardiac myocytes (Cheng et al., 1996).

View larger version (115K): [in this window] [in a new window]   FIGURE 1   Simulated line-scan image of myoplasmic calcium ([Cai2+] = c, represented by color) of a cardiac myocyte as a function of space (horizontal axis) and time (vertical axis). The reciprocal of the slope of the wave front gives the wave speed, v = 67 µm s1. The simulation includes an array of 50 spatially discrete Ca2+ release sites, two spatially homogeneous fluxes (Ca2+ leak and SERCA pumps), and Ca2+ diffusion. See Appendix for methods.

If Ca²⁺ release from the isolated sites in Fig. 1 were replaced with a continuous, uniform rate of the same magnitude per unit length, the local medium would be excitable, i.e., increasing [Ca_i²⁺] above a threshold (~0.14 µM) would cause an action-potential-like spike of Ca²⁺. An excitable medium would support a traveling wave pulse with a speed proportional to the square root of the diffusion constant (Murray, 1989; Tyson and Keener, 1988), as predicted for Ca²⁺ waves in immature Xenopus oocytes (Jafri and Keizer, 1995). Thus we carried out additional simulations in which either the diffusion constant of calcium, D, or the distance between release sites, d, was varied. The main graph in Fig. 2 shows that v is approximately proportional to D/d, rather than (Jaffe, 1993). We find that if d is too large or D is too small, waves do not propagate.

View larger version (14K): [in this window] [in a new window]   FIGURE 2   Main graph: Wave speed, calculated as in Fig. 1, is approximately a linear function of D/d. Open squares are simulations with d = 2.0 µm and D varied; waves do not propagate for D < 10 µm2 s1. Filled circles are simulations for D = 30 µm2 s1 and d varied; waves do not propagate for d > 3.0 µm. Slope of full line is 4.5. Inset: Open squares, axes, and line as in main graph; filled circles for D = 30 µm2 s1 and d varied with site source strength per unit length held fixed to simulate the continuum limit; the wave speed in the continuum limit, vc, is achieved when D/d  vc (dotted lines). Other parameters as in Fig. 1.

The classical continuum limit for these simulations involves shrinking the separation between sites (d) to zero while maintaining a fixed release and re-uptake rate per unit length. Thus we have investigated the continuum limit by taking J_site, J_leak, and J_serca in Eq. 10 proportional to d and repeating the simulations in Fig. 1 with successively smaller values of d. A plot of the wave speed versus D/d is given in the inset to Fig. 2. The speed in the continuum limit is indicated there by v_c, and the transition to the continuum limit is seen to occur when D/d  v_c. These results make it clear that site separation significantly alters the mode of propagation of the wave.

	FIRE-DIFFUSE-FIRE MODEL

To investigate why spark-mediated wave propagation differs so much from continuum propagation, we consider a caricature of the spark-mediated wave. In this simplified model release sites are located at the points x = nd (n = 0, ±1, ±2, ···) and instantaneously release a fixed amount, , of Ca²⁺ when c (=[Ca_i²⁺]) at a site exceeds a threshold value, c*. After release the site becomes refractory. However, the released Ca²⁺ diffuses and may trigger another instantaneous release (a spark) at neighboring sites. This is illustrated in Fig. 3 using overbars for the dimensionless variables:  = x/d (distance measured in terms of the site separation),  = tD/d² (time measured in terms of the time required to diffuse between sites), and  = c/c* (concentration measured in terms of the threshold concentration for Ca²⁺ release). We refer to this model as "fire-diffuse-fire" since a wave propagates by sequentially triggering Ca²⁺ sparks to the right (or left) by diffusion. Although the fire-diffuse-fire model is greatly simplified, it includes the essential features of fast adaptation and refractivity in that release does not occur over a sustained period and that once a site has released Ca²⁺, it cannot release Ca²⁺ again.

View larger version (8K): [in this window] [in a new window]   FIGURE 3   Schematic illustration of propagation for the "fire-diffuse-fire" model of spark-triggered waves (see text). Nondimensional variables  = x/d,  = tD/d2, and  = c/c* are used, giving the threshold * = 1 shown by the dashed line. The dimensional release rate for all release sites is localized in space and time using Dirac delta function spikes, i.e., Jsite = (x  nd)(t  tn), with tn the time that the nth site fires. The nondimensional parameter for the amount of Ca2+ release is  = c*d/. The dotted vertical lines schematically show the instantaneous release of Ca2+ from sites  = ±1 at  = 0.5 due to a release event at  = 0 from the site  = 0. The Gaussian curves were calculated using the analytical formula (Murray, 1989)  = exp(2/4)/ with  = 1/. This is the largest value of  for which release from the site at  = 0 can equal the threshold * = 1 and trigger release at its neighbors, which can be seen setting  = 1 in the analytical expression for and setting  = 1.

The partial differential equation governing this model is

(1)

where  is the Dirac delta function and _i is the time that site i fires. The dimensionless parameter,  = c*d/, which governs the dynamics of the fire-diffuse-fire model, is the ratio of the threshold concentration for CICR (c*) to the concentration due to release by a single site (/d). The value of _i can be obtained recursively (see Eq. 6).

Thus the mean speed of a right-going wave front at site n is given by the simple formula v_n = d/_n where _n = t_n  t_n1 is the time interval between the firing of the spark at site n  1 and site n. Or in terms of _n,

(2)

For the release event shown in Fig. 3 to trigger release from its neighboring sites, the value of at  = ±1 must equal 1. The first time that this occurs, if ever, is ₁, which can be obtained from the relationship  = 1, i.e.,  = exp(1/4₁)/ (see legend to Fig. 3). This can occur only when   1/  0.24. This agrees with the intuition that wave initiation is favored by a low threshold (c*), sites that are close together (d), and large releases of Ca²⁺ (); thus not all regular arrays of sites can initiate a wave from the firing of a spontaneous spark. For the initiation step illustrated in Fig. 3, ₁ = 0.5 and the initial speed is v₁ = 2D/d.

The speed of the wave front increases as subsequent release events contribute to the Ca²⁺ profile. The Ca²⁺ released by a site at n contributes to the overall profile of via the formula (Murray, 1989)

(3)

If the wave speed becomes constant (as in Fig. 1), then the interval between the firing of successive sites becomes a constant, i.e., _n =  for n large enough. In this case, by summing up the contributions in Eq. 3 over all sites it can be shown that

(4)

which defines g().

Equation 4 has a single root, (), that gives the steady wave speed

(5)

Since Eq. 4 is independent of the diffusion constant, the wave speed is proportional to D, as found in the simulations in Fig. 2. This result also explains why the speeds of the spark-dependent waves in Figs. 1 and 2 are approximately proportional to D/d. For   1 (i.e., large source), Eq. 4 implies that  = 1/[4 ln(1/)] = 1/[4 ln(/dc*)]. Thus v = 4(D/d)ln(/dc*) and the weak dependence of the wave speed on ln(1/d) would be difficult to detect. The value of () can be obtained graphically from the plot of g() in Fig. 4.

View larger version (14K): [in this window] [in a new window]   FIGURE 4   Graphical calculation of the firing interval, , between adjacent sites at the front of a wave. The value of for which the function in Eq. 4, g(), equals  is the firing interval for that value of . Since g()  1, waves do not propagate for  > 1. This result is easily generalized for a spatially uniform uptake and release rate of the form (  0). For non-zero values of the right-hand side of Eqs. 3 and 4 are multiplied by exp() and exp(), respectively. If 0 > 0, then the left-hand side of Eq. 4 must also be multiplied by 1  (0/*). The resulting sum is plotted for  = 0.1 and 0.5.

Despite the fact that a spontaneous release event cannot initiate a wave if  is larger than ~0.24, the simultaneous release from several sites might initiate a wave with larger values of . Indeed, since Fig. 4 shows that g can get as large as 1, it seems possible that waves could propagate with 0.24  g  1. We have explored this further by applying the formula in Eq. 3 and solving iteratively for _n. Indeed, when all the _i for (n  1)  i  n  1 are known, then _n = _n  n1 can be obtained by solving the following equation for _n:

(6)

The results of this procedure, after waiting long enough for convergence, are shown in Fig. 5. Below  = 0.512 the interval between successive firings, , converges to a constant. This value of , however, is a critical point at which a period doubling bifurcation occurs, i.e., successive firing intervals alternate between a longer and a shorter value. This period doubling continues, leading to an apparently chaotic state for  large enough. For values of   0.535 the chaotic attractor ceases to exist and waves do not propagate. Thus propagation failure occurs via period doubling to chaos, well below the limit of  set by Fig. 4.

View larger version (13K): [in this window] [in a new window]   FIGURE 5   The "firing interval", i.e., the time interval between adjacent sparks at the front of the wave, is obtained by successively calculating n for n = 1, 2, ...  using Eq. 3 and the criterion (n, n) = 1 for a range of values of . The wave was initiated by simultaneously firing all the sites for 15  n  15. A period doubling cascade begins at   0.512, which terminates in a chaotic state at   0.535, beyond which waves do not propagate.

The period doubling cascade produces a rhythmic alteration in the progress of the wave front. For  < 0.512 the time interval between Ca²⁺ release at the frontmost site (n) and the next site (n + 1) is fixed. Using Eq. 5, this implies that the wave front propagates at a fixed speed. This steady propagation changes at the first period doubling bifurcation, where the time interval between release at site n and n + 1 becomes slightly longer than the interval between sites n + 1 and n + 2; thus the speed of the front alternates between a slower and faster value, giving a slightly jerky appearance to the front. At the second period doubling bifurcation the wave propagation is more complex, with four alternating speeds. This continues with propagation becoming increasingly more complex at each period double bifurcation until finally propagation failure occurs at the chaotic state. Although these specific dynamical features are not seen in simulations with a two-dimensional array of sites (not shown), a rich variety of complex dynamics are still observed near the propagation failure limit.

Similar results are obtained if we introduce a linear, spatially uniform uptake and release of Ca²⁺ from the SR with rate constant , as described in the legend to Fig. 4. The plots of g() for  = 0.1 and 0.5, however, now have a maximum. Thus for  > 0, the solution to Eq. 4 has two roots. The larger value of corresponds to the intersection on the declining branch of the curve in Fig. 4 and a wave that is slower and unstable.

	DISCUSSION

Top Abstract Introduction Discussion Appendix References

The saltatory nature of Ca²⁺ wave propagation in cardiac myocytes has been revealed only recently using high-speed, high-resolution line scan images (Cheng et al., 1996). The correlation between the underlying Ca²⁺ sparks and the regular array of t-tubule structures in these images provided the motivation for the one-dimensional simulations that we report here. By using a model of a release site based on the kinetics of ryanodine receptors (Smith, 1996; Keizer and Smith, 1998) we have found that these waves propagate with a speed that is proportional to the diffusion coefficient of Ca²⁺, rather than its square root as prediction by conventional reaction-diffusion equations (Jaffe, 1993; Jafri and Keizer, 1995). Our results in Fig. 2 (inset) illustrate the connection between these two extreme types of wave propagation, saltatory and continuous. In the saltatory case the value of [Ca_i²⁺] at the wave front is dominated by release from a single site (c.f. Fig. 1). For continuous propagation, release is distributed continuously in space and many sites at the front release Ca²⁺ simultaneously.

The "fire-diffuse-fire" model, which we have introduced to help explain saltatory propagation, is a simplified model of CICR by release sites with a refractory state. In this model a site releases Ca²⁺ instantaneously ("fires") when the value of [Ca_i²⁺] at the site exceeds a threshold value. To mimic a long-lasting refractory state, once a site has released Ca²⁺, it can no longer fire again. To mimic the regular array of t-tubules in myocytes, the release sites are located with a fixed separation, d, and Ca²⁺ released at one site diffuses continuously with an "effective" diffusion constant, D, due to the presence of myoplasmic buffers. In this model the speed of the wave front is determined by the time it takes Ca²⁺ released by the site at the front to diffuse to the next active site and raise the value of [Ca_i²⁺] there to the threshold (c.f. Eq. 2). In contrast to the kinetic model, the fire-diffuse-fire model leads to analytical expressions for the wave shape and the wave speed, and therein lies its value. Indeed, the simple result illustrated graphically in Fig. 3 makes it clear that the wave speed is proportional to the diffusion constant since the time to diffuse between release sites (d²/D) is inversely proportional to the diffusion constant. The model also suggests that propagation failure of saltatory waves may be quite complex.

The simplifying assumption of instantaneous release of Ca²⁺ in the fire-diffuse-fire model is not responsible for the saltatory nature of the waves. Indeed, the simulations with the kinetic model in Fig. 1 do not make this assumption and yet exhibit saltatory propagation. We have investigated this further using a generalization of the fire-diffuse-fire model in which release is not instantaneous (J. Pearson, J. Keizer, and S. Ponce-Dawson, unpublished observations). We find that a key dimensionless number is D/d² where  is the mean time that a site is open and d²/D is the intersite diffusion time. When D/d²  1, propagation is saltatory and the wave speed is proportional to D, as we have shown. In the saltatory limit, propagation consists of isolated bursts of Ca²⁺ that occur as each consecutive site fires. When D/d²  1, propagation is continuous, the velocity is proportional to , and many sites are releasing Ca²⁺ simultaneously.

This analysis explains the transition from saltatory to continuous propagation shown in the inset to Fig. 2. According to the Luther equation the continuum wave speed should be given by the expression v_c =  (Jaffe, 1993). In terms of v_c, the dimensionless parameter D/d² can, therefore, be rewritten as (D/dv_c)². Thus the analysis in the previous paragraph predicts that the transition between saltatory and continuous propagation occurs when D/d  v_c. Moreover, due to the second-power dependence of the ratio (D/dv_c)², the transition should be relatively sharp. This agrees with both the location and sharpness of the transition for the simulations plotted in the inset of Fig. 2. It is easy to evaluate the ratio D/d² using experimental data for cardiac myocytes (D = 30 µm² s¹,  = 14 ms, and d = 2.0 µm; Smith et al., 1998). This gives D/d² = 0.1, which is well into the saltatory regime. This prediction is compatible with the punctate images of the wave front that result from enhancement of Ca²⁺ waves in myocytes (Cheng et al., 1996).

The fire-diffuse-fire model suggests experimental tests that distinguish saltatory from continuously propagating Ca²⁺ waves. The linear dependence of v on D should show up in the dependence of wave speed on temperature. Indeed, from the Arrhenius equation, D = D_oexp(E_a/RT), where E_a is the activation energy for diffusion. Generally, this is much smaller than the activation energy for the biochemical processes that determine the activation energy associated with . This implies that the speed of a saltatory wave, which is proportional to D/d, should be less sensitive to changes in temperature than predicted by the Luther equation, which gives v_c = . Using exogenous Ca²⁺ buffers it also should be possible to manipulate the diffusion constant of Ca²⁺ (Wagner and Keizer, 1994), which would alter the wave speed differently for saltatory and continuous propagation.

Several predictions of the fire-diffuse-fire model are in agreement with published observations of waves in cardiac myocytes. Thus the saltatory structure of the waves in the fire-diffuse-fire model, when displayed as in Fig. 1, are similar to the images obtained from line-scan data in myocytes (Cheng et al., 1996). The fact that there is a maximum value for  = c*d/ above which waves do not initiate is compatible with the fact that Ca²⁺ overloading of myocytes is a prerequisite for waves, i.e., the source strength, , must be sufficiently large (Cheng et al., 1996). The dependence of the wave speed on  (c.f. Fig. 5) may help explain irregularities in the speed of the wave front as it moves across a myocyte (Cheng et al., 1996) since both d and store content () may vary somewhat within a cell. Finally, the observation that wave propagation is irregular just before propagation failure may be related to the onset of chaos that is seen in the fire-diffuse-fire model (Cheng et al., 1996).

What do these results tell us about the physiological role for saltatory propagation of Ca²⁺ waves in cells? In fact, the saltatory Ca²⁺ waves observed in cardiac myocytes and the immature Xenopus oocyte (Callamaras et al., 1998) both occur only under extreme physiological conditions. Therefore, it is possible that the punctate arrangement of release sites in these cells actually functions to inhibit Ca²⁺ waves by guaranteeing propagation failure under normal physiological conditions. This would contrast with saltatory propagation of action potentials in myelinated nerve, where the nodes of Ranvier are analogous to release sites (Fitzhugh, 1962). In this case, the saltatory wave speed has been reported to exceed that for the continuum limit. This difference from what we find for the Ca²⁺ wave in myocytes (c.f. Fig. 2, inset) is probably due to the finite size of the nodes and may reflect the different physiological roles of CICR and action potentials in the two cell types. On one hand, in muscle Ca²⁺ release is believed to act locally and, thus, failure of a propagating Ca²⁺ signal would be appropriate physiologically. In myelinated nerve, on the other hand, action potentials are transmitted to distant cells and nodes are utilized to reinforce their propagation.