GLOSSARY

TOP ABSTRACT GLOSSARY INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX A APPENDIX B REFERENCES

Most symbols were defined in the companion paper (Izu et al., 2001). Here we list new symbols and those that occur frequently in this paper.

x, y	spatial coordinates (µm)
t	time coordinate (ms)
C	free Ca2+ concentration (µM)
DCx, DCy	Ca2+ diffusion coefficients along x and y (µm2/ms)
S	stochastic switching function equaling either 0 or 1
Pmax	maximum probability of Ca2+ spark occurrence/calcium release unit/ms
P	probability of Ca2+ spark occurrence/calcium release unit/ms
n	Hill coefficient in definition of P
K	Ca2+ sensitivity parameter in definition of P (µM)
x, y	spatial separation of Ca2+ release units (CRUs) along x and y (µm)
	molar flux of CRU for 2-dimensional diffusion domains (pmol/ms/µm)
ISR	current through the CRU (pA)
Topen	duration of current flow through CRU (ms)
C0	baseline Ca2+ concentration (µM)
	waiting time distribution
p(t < )	probability of a CRU not firing in time t <
P(X  1, t < )	probability of at least one CRU firing in time t <

	INTRODUCTION

TOP ABSTRACT GLOSSARY INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX A APPENDIX B REFERENCES

Slowly (~100 µm/s) propagating waves of elevated Ca²⁺ concentration (Ca²⁺ waves) appear to be a ubiquitous finding and have been observed in a wide diversity of cell types including skeletal (Endo, et al., 1970) and cardiac muscle (Fabiato and Fabiato, 1972), medaka eggs (Ridgway, et al., 1977), and astrocytes (Cornell-Bell and Finkbeiner, 1991). In cardiac atrial cells that lack transverse tubules (T tubules), physiological activation of the myofibrils in the interior of the cell depends on the propagation of a Ca²⁺ wave from the sarcolemma to the cell's center (Berlin, 1995; Hüser et al., 1996). In cardiac ventricle cells, Ca²⁺ waves are not physiological and are believed to be a pathological manifestation of Ca²⁺ overload and might trigger ventricular arrhythmias (Lakatta and Guarnieri, 1993).

Ca²⁺ waves are a natural consequence of regenerative Ca²⁺ release by the sarcoplasmic reticulum or Ca²⁺-induced Ca²⁺-release (CICR) (Endo, et al., 1970; Ford and Podolsky, 1970). All models of Ca²⁺ waves in skeletal and cardiac muscle have CICR at their core and the evolution of Ca²⁺ wave models reflects the growth of knowledge of excitation-contraction coupling in muscle. For example, an early mathematical model of Ca²⁺ waves (Backx et al., 1989) distributed Ca²⁺ release sites uniformly throughout the cytoplasm. However, in normal muscle cells, Ca²⁺ release sites are distributed in discrete bands at the z-lines (Carl, 1995). With the discovery of Ca²⁺ sparks (Cheng et al., 1993) the discrete nature of Ca²⁺ release was resolved to the opening of a few ryanodine receptors (RyRs) on the z-line in cardiac muscle (Shacklock et al., 1995). In addition, Cheng et al. (1993, 1996), speculated that Ca²⁺ waves arise from the collective firing of Ca²⁺ sparks.

Contemporary models of Ca²⁺ sparks and Ca²⁺ waves restrict Ca²⁺ release to the z-lines (Keizer et al., 1998; Keizer and Smith, 1998; Izu et al., 1999; Lukyanenko et al., 1999). However, understanding the relationship between Ca²⁺ sparks and Ca²⁺ waves in living cells has been complicated by asymmetries in cell structure, asymmetries in the Ca²⁺ spark profile, anisotropic diffusion of the Ca²⁺-bound fluo-3, uncertainties regarding the Ca²⁺ sensitivity of the Ca²⁺ release units (CRU), and the Ca²⁺ currents required to generate Ca²⁺ sparks.

In this paper, we present a model that provides a unified framework for studying both Ca²⁺ sparks and Ca²⁺ waves. The following are key elements in this model.

1.	Large CRU currents. By using relatively large CRU currents (5-20 pA; Izu et al., 2001) the model generates realistic Ca2+ sparks and realistic Ca2+ waves while maintaining CRU Ca2+-sensitivities compatible with experimental measurements (Lukyanenko and Györke, 1999).
2.	Stochastic triggering of Ca2+ release. Previous models of Ca2+ waves (except Keizer and Smith, 1998) have used a deterministic rule for triggering Ca2+ release; release occurred when the Ca2+ concentration exceeded some fixed value C*. The use of stochastic instead of deterministic dynamics is important for two reasons. First, sparks occur spontaneously and randomly (Cheng et al., 1993) so a unified framework for studying sparks and waves must provide for the stochastic firing of the CRUs. Second, deterministic and stochastic systems can behave very differently when the Ca2+ sensitivity of the CRUs is high. In a deterministic system, waves can occur even when the sensitivity is arbitrarily high. However, in a stochastic system, even with moderately high CRU Ca2+ sensitivity (~5 µM), the large number of spontaneous sparks would trigger so many Ca2+ waves at once that observing a well-defined wave would be almost impossible.
3.	Asymmetric distribution of discrete CRUs, and
4.	anisotropic diffusion of Ca2+ and of mobile buffers. An important difference between our model and that of Keizer and coworkers (Keizer et al., 1998; Keizer and Smith, 1998) is how the Ca2+ buffers (endogenous and Ca2+ indicator) are handled. The Keizer models do not include buffers, and they compensate for the absence of buffers by reducing the free Ca2+ diffusion coefficient ~10 fold. We will show that buffers endow the system with a property we call "superadditivity" that has profound effects on Ca2+ signaling, which extend much further than simply slowing Ca2+ diffusion.

The effect of three of the four key features in this model on Ca²⁺ wave propagation has been investigated by others. The effect of discrete, asymmetric distribution of CRUs was studied by Bugrim et al. (1997) and anisotropic Ca²⁺ diffusion was addressed by Kargacin and Fay (1991) and Girard et al. (1992). Keizer and Smith (1998) examined the effect of stochastic firing of CRUs on wave propagation. What is novel in our model is the combination of large CRU currents, stochastic triggering of CRU release, asymmetric distribution of CRUs, and anisotropic diffusion.

	METHODS

TOP ABSTRACT GLOSSARY INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX A APPENDIX B REFERENCES

Figure 1 shows the geometry of our model. A is the 3-dimensional (3D) schematic of a ventricular cell. The x-direction is the cell's longitudinal axis. The three vertical planes occur at the z-lines and are spaced a distance _x (=2 µm) apart. The black dots in the y-z plane are the CRUs. The horizontal plane is the 2-dimensional (2D) slice on which we will carry out our simulations. At present, we cannot do a full 3D simulation because of limitations in computational power. Figure 1 B shows CRUs in a y-z plane. Ca²⁺ release by discrete CRUs will generate concentration gradients in all directions. To eliminate gradients along z, the discrete CRUs along z are replaced by line sources that extend from  < z <  and spaced a distance of _y (0.4 or 0.8 µm) along y as shown in Fig. 1 C. These infinite line sources induce symmetry in z making planes at any z equivalent, thereby reducing the problem from three to two dimensions. Figure 1 D shows the plane on which the model is defined. The line sources intersect the x-y plane and at regular intervals of _x along x and _y along y. These intersections are called lattice sites.

View larger version (19K): [in this window] [in a new window]   FIGURE 1   Geometry of model.(A) Three-dimensional schematic of a cardiac ventricular cell. The longitudinal axis of the cell is in the x-direction. The three vertical planes represent the z-lines and are spaced a distance of x (= 2 µm) apart. The y-z plane at the z-line contain the CRUs (dots) that are symmetrically spaced along the y and z directions. The horizontal plane shown in the center is the domain on which the model equations are defined. (B) An individual y-z plane and its CRUs. Ca2+ release by discrete CRUs will generate concentration gradients in all directions. To eliminate gradients along z, the discrete CRUs along z are replaced by line sources that extend from  < z <  as shown in C. (D) The 2D slice of the ventricular cell with CRUs represented by line sources spaced x apart along x and y (= 0.4 or 0.8 µm) along y.

The model equations are the same as Eqs. 2-7 in Izu et al. (2001), that describe the so-called "Smith buffer model." Apart from the restriction to two spatial dimensions, the only difference here is the differential equation for the free Ca²⁺ concentration C(x, y, t) that is now

(1)

The summation term is new to this paper. Each term in the summation represents a point source (the CRU) located at the lattice site (x_i, y_j) that produces a molar flux . The lattice sites are spaced _x apart (2 µm, the sarcomere length) along the longitudinal axis of the cell (x-axis) and _y apart (0.4 or 0.8 µm) in the transverse direction (y-axis). Figure 3, inset, and Fig. 4 show the lattice. CRUs on a column (fixed x) are said to be on a z-line. S is a stochastic function taking values of either 0 or 1, so switches the CRU on (firing) or off. After S becomes 1, it stays at this value for time T_open. The probability that the CRU will fire in time t is P(C(x, y, t), K, n) · t, where P is the probability of firing per unit time. P is a function of the ambient Ca²⁺ concentration C(x_i, y_j, t) and is given by

(2)

P was determined as follows. Let r be the number of sparks/100 µm linescan/ms. Assume that the microscope's lateral (y) and axial resolution (z) is 0.5 and 1 µm, respectively. If CRUs are arranged in a square lattice on the z-line plane with spacing _y = 0.8 µm, then there will be ~2 CRUs (or ~6 for _y = 0.4) in a 0.5 µm × 1 µm confocal sample area at each z-line. Assuming a sarcomere length of _x = 2 µm, then, in a 100-µm confocal linescan, there are N_CRU = 2 CRU/z-line × 50 z-lines/100 µm linescan = 100 CRU/100 µm linescan. Then P = r/N_CRU. Note that P_max has units of sparks/ms/CRU. The spark rate, the Hill coefficient n = 1.6 and the Ca²⁺-sensitivity factor K = 15 µM are taken from Lukyanenko and Györke (1999). Using a spark frequency of 10 sparks/100 µm linescan/s at a Ca²⁺ concentration of 100 nM, for N_CRU = 1 CRU/µm, P_max equals 0.3/CRU/ms (or 0.05/CRU/ms for _y = 0.4). Note that P_max is kept at these values (0.3 or 0.05) in simulations where K is varied.

The source strength  requires some explanation. Because a line source is not equivalent to a linear array of discrete CRUs (shown in Fig. 1 B), we need to somehow adjust the molar flux of the line source to approximate the molar flux of a point source. If Ca²⁺ were being released from a CRU into a 3D volume, then ₃ = I_SR/2F, where I_SR is the current and F is the Faraday. Note that ₃ has units of mole/ms. Because the model is 2D,  = ₂ has units of mole/ms/µm. What value should we use for ₂ in place of ₃? Although no value of ₂ will give the identical space-time Ca²⁺ distribution in 2D as ₃ will in 3D, we define an "equivalent" source strength ₂ as one that gives the same concentration as ₃ at r =  (the Euclidean distance) and t =  in a linear system. Let C(r, t, d) be the concentration at (r, t) in d = 2, 3 dimensions. Then C(r, t, 3) = (₃/(4Dr))erfc(z), C(r, t, 2) = (₂/(4D))E₁(z²) where z = r/ and E₁ is the exponential integral (Appendix A). For  = 0.5 µm,  = 5 ms, and D = 0.2 µm²/ms (the geometric mean of D_Cx and D_Cy), ₂ = 0.64 ₃/µm. We used this numerical conversion in all our simulations. With this conversion factor, the sparks in the 2D model has about the same spatial spread and time course as the 3D spark. For example, for I_SR = 20 pA the spark FWHM along x at the end of 5 ms was 2 µm in 3D (Izu et al., 2001) and 1.9 µm in 2D.

We will vary I_SR, T_open, K, and _y. The remaining parameters use the following (standard) values (see Izu et al., 2001 for references): _x = 2 µm (Shacklock et al., 1995), n = 1.6 (Lukyanenko and Györke, 1999), H_B = 123 µM, k_B⁺ = 100/µM/ms, k_B = 100/ms, H_D = 50 µM, k_D⁺ = 80 µM/s, k_D = 90/s, V_p = 200 µM/s, K_p = 184 nM, n_p = 4, D_Cx = 0.3 µm²/ms, D_Cy = 0.15 µm²/ms, D_Dx = 0.02 µm²/ms, and D_Dy = 0.01 µm²/ms, and C₀ = 0.1 µM.

We used the reaction of Ca²⁺ with fluo-3 described by Smith et al. (1998) instead of the more complex set of reactions proposed by Harkins et al. (1993) because of computational limitations. In our companion paper (Izu et al., 2001), we found that the Harkins buffer model produced sparks with more realistic F/Fo and slightly larger spatial spread than did the Smith model. However, computations took considerably longer with the Harkins model and required a finer spatial discretization than the Smith model, making the computational costs prohibitive. We consider possible ramifications of using the Harkins model in the Discussion.

Numerical methods

The model equations were discretized and numerically solved using Facsimile as described in Izu et al. (2001). The stochastic opening of a CRU was simulated as follows. At regular intervals of time t (typically 1 ms), a uniformly distributed random number u between 0 and 1 is generated. Note that, as C  , P/P_max  1 so the probability that u is less than P/P_max approaches 1. Thus, the CRU fires when u < P/P_max by setting S to 1, and S remains unity for T_open ms. After closing, the channel does not reopen.

The computational grid was a rectangle 20 µm × 20 µm with mesh size of 0.1 µm along both x and y directions. We were constrained to use these small domains because a typical simulation of 150 ms takes ~24 hrs on a 400 MHz Pentium II processor with 256 MB RAM, and computation time grows faster than computational area. Neumann boundary conditions were imposed on all edges.

	RESULTS

TOP ABSTRACT GLOSSARY INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX A APPENDIX B REFERENCES

The basic principle behind wave evolution in this model is simple. Waves originate and are sustained as Ca²⁺ released from one or more CRUs diffuse to neighboring CRUs and raise the probability that these neighbors will release Ca²⁺ (or "fire") and, in turn, induce other neighbors to fire. Our model of Ca²⁺ waves is analogous to piles of gunpowder wherein igniting one pile may ignite the neighboring pile. The probability that the neighbor will ignite depends on how far the piles are from each other (_x, _y), the rate of heat production (I_SR), the total amount of heat produced from each pile (I_SR × T_open), and the thermal stability of the gunpowder (K). The challenge is to get a quantitative understanding of the conditions needed to initiate a wave, of how the wave evolves, the wave velocity, etc.

We start by analyzing a deterministic linear model for waves on lattices even though it turns out to be inadequate for understanding waves in muscle cells for the following reasons.

1.	The linear analysis is a natural first attempt because linear systems are simple and exactly solvable. So if they prove to be accurate descriptions of waves in nonlinear systems, then virtually everything about waves will be known.
2.	Because the linear solution is given in terms of known functions (see Eq. 3), it provides a standard to check the accuracy of the numerical algorithm used to solve the model equations.
3.	Most importantly, the failure of the linear approximation to accurately predict properties of waves in the nonlinear system reveals an essential difference between nonlinear and linear systems: superadditivity. We examine the origin of superadditivity and show its importance in understanding waves in cells.

Next, we examine the propagation of a Ca²⁺ wave and extract the basic steps in the propagation of waves on rectangular lattices. Because the firing of CRUs is stochastic, we need to develop the probabilistic machinery to estimate whether a given set of model parameters will be able to support a wave, predict the wave's speed, and estimate the number of CRUs required to trigger a wave.

Linear systems

We would know virtually everything about waves if we could predict the Ca²⁺ concentration at any point (x, y, t) given the history of firings of CRUs. This prediction is not possible in nonlinear systems but is possible in linear systems. Linear systems, endowed with the property of additivity, allow us to calculate C(x, y, t) knowing the space-time distribution of Ca²⁺ from a single spark. C(x, y, t) is given by

(3)

The addend, derived in Appendix A, gives the space-time Ca²⁺ distribution for a single CRU with constant molar flux q_ij that is open for time T_open. Note that the solution takes account of anisotropic diffusion. The summation is taken over all sites that have fired. t_ij is the time the channel at (x_i, y_j) had fired. This fundamental equation contains all the information about wave propagation on discrete lattices in a linear deterministic system.

The analytic power of linear systems stems from the ability to predict the Ca²⁺ concentration at any point (x, y, t) due to the firing of an arbitrary distribution of CRUs given the Ca²⁺ distribution due to a single CRU. The objective here is to fit the Ca²⁺ distribution from a single spark in the presence of buffers and pump (nonlinear system) to that of a linear system (see Appendix A),

(4)

that defines an equivalent flux q and diffusion coefficient D_x. D_y always equaled D_x/2. Figure 2, A and B, show the concentration profiles from the simulations (circles) and the fitted function (solid curve) at t = 5 ms (left profile) and 10 ms (right profile) after the channel opens. Panel A shows the profile along x, and B, the profile along y. Figure 2, C and D, show q and D_x; in a linear system, these parameters would be time independent. Note that the values of D_x and q determined from the x profile (up triangle) and y profile (down triangle) are almost identical. This equality is important because it means that the Ca²⁺ concentration at an arbitrary point (x, y, t) can be calculated using Eq. 3 using a single pair (q, D_x) for any given moment. Time-dependent changes in q and D_x simply arise from fitting the solution to the nonlinear problem to that of a linear problem. If the buffers were absent, then the system would be linear and q and D_x would be invariant. There is no simple physiological interpretation for the time dependency of q and D_x.

View larger version (16K): [in this window] [in a new window]   FIGURE 2   Determination of "equivalent" q and Dx. Single spark profiles from simulations (circles) were fit to Eq. 4 by adjusting q and Dx. (A) The profile along x at t = 5 ms (r < 0) and t = 10 ms (r > 0). (B) Similar to A, but shows the profile along y. The best fit values of q and Dx as functions of t are shown in C and D. Up-triangles show the fit of the spark profile along x; down-triangles, along y. Simulations carried out using ISR = 20 pA, Topen = 10 ms.

It is important to observe in Fig. 2 that the Ca²⁺ concentration is a very steep function of distance; it drops 1000-fold within 1.5 µm from the source. A CRU at x = 2 µm would "see" a Ca²⁺ concentration of just 179 nM, 10 ms after the channel opens. Even if this concentration could be maintained indefinitely, the mean waiting time for the CRU to fire is 1/P(179 nM, 15 µM, 1.6) = 3.98 s (see Appendix B), so a single spark is unlikely to trigger a spark at an adjacent sarcomere. This calculation is in accord with Parker et al.'s (1996) observation that the probability of a spark triggering another 2 µm away was zero.

Triggering of sparks on adjacent sarcomeres is the sine quo non of wave propagation in our model. Thus multiple channels must fire within a short time interval to raise the Ca²⁺ concentration at an adjacent sarcomere enough to trigger Ca²⁺ release. We calculated C(x, y, t) using Eq. 3 and q, D_x values at t = 10 ms for different configurations of simultaneously firing CRUs shown in Fig. 3. The channel configurations are labeled (a) 1 × 1, (b) 1 × 3, (c) 1 × 7, (d) 2 × 7, and (e) 1 × 9, where the first and second numbers refer to the number of columns and rows of CRUs, respectively. The inset in Fig. 3 shows a schematic of the lattice, the arrangement of CRUs in a 2 × 3 configuration. The CRUs in a column are symmetrically displaced about the origin so, for example, the coordinates of the 3 CRUs in the 1 × 3 configuration are (0, 0), (0, +_y), and (0, l_y), where _y = 0.8 µm. The rationale for choosing these configurations will become evident later. For configuration (d) the 2 columns are at x = 0 and 2 µm. Fig. 3 A shows the concentration along the x-axis (y = 0) at t = 10 ms for configurations (a-e). Because of the steep Ca²⁺ concentration profile, the firing of additional CRUs in b-e do not add substantially to the concentration due to the single CRU at the origin (a, lower curve). At this resolution, the profiles for channel configurations b, c, d, and e are indistinguishable. Thus the linear model predicts that none of these configurations of firing channels could raise the Ca²⁺ concentration at x = 2 µm sufficiently to substantially increase the probability of firing within 10 ms.

View larger version (21K): [in this window] [in a new window]   FIGURE 3   Ca2+ concentration (C(x, y, t)) in (A) additive and (B, C) superadditive systems. (A) and (B) show the Ca2+ spatial distribution along x at y = 0 at t = 10 ms due to firing of CRUs in the following configurations: (a) 1 × 1, (b) 1 × 3, (c) 1 × 7, (d) 2 × 7, and (e) 1 × 9. Inset, a schematic of the lattice with a 2 × 3 configuration of CRUs. The darker circle indicates the CRU at the origin. In an additive system (A) simultaneous firing of multiple CRUs (b-e) do not add significantly to the Ca2+ concentration from a single CRU (a). A buffered system exhibits superadditivity, and the Ca2+ distribution due to firing of multiple CRUs is remarkably different. Note that there is little difference due to firing of 7 or 9 CRUs (c, e). (C) The temporal behavior of the Ca2+ concentration at x = 2, y = 0 in a buffered system. Simulation parameters: ISR = 20 pA, Topen = 10 ms.

Superadditivity in buffered systems

A buffered system behaves differently, however. Fig. 3 B shows Ca²⁺ concentration profiles derived from simulations of the nonlinear system (buffers reactions and pump included) for the same five configurations (the graphs for c and e are virtually identical). In contrast to the linear system, the firing of the flanking CRUs greatly increased the Ca²⁺ concentration along the x-axis.

The concentration profiles in Fig. 3 A are expected if the effect of two or more sparks were additive. However, Fig. 3 B shows that the actual concentration exceeds the expected values. We call this property superadditivity. Superadditivity arises simply because the free Ca²⁺ buffers become depleted. The addition of m moles of Ca²⁺ causes a drop in the available free buffer, so subsequent additions of m moles result in a larger increase in the Ca²⁺ concentration.

A consequence of superadditivity important for understanding Ca²⁺ wave propagation is that, when CRUs fire simultaneously, the range of Ca²⁺ signaling increases tremendously. This is seen in Fig. 3 B and is also illustrated in Fig. 3 C, where the concentration at (x = 2 µm, y = 0) is plotted against time. Note the tremendous amplification in the Ca²⁺ concentration as more CRUs fire. The maximum concentration attained 2 µm away from a single spark is 278 nM (a), whereas, at the same distance from the center of seven sparks (c), the maximum concentration is 11.3 µM, 41 times larger.

Wave propagation in nonlinear systems

For a linear (additive) system, we can make a priori predictions of wave properties, such as the wave speed based on the properties of a single spark. However, because of superadditivity in buffered systems, such predictions are not possible. This does not mean, however, that there are no predictive tools for nonlinear waves. Instead, the probabilistic tools we will develop below require more information than that provided from a single spark as input data and require some knowledge of how waves propagate on lattices. Thus we first examine in detail how a wave propagates on a rectangular lattice. Once we understand the basic steps of wave propagation, we can set up approximations to the Ca²⁺ distribution in a wave, which are used in the probabilistic equations. These equations are used to predict wave properties, such as whether waves can exist for a given set of parameters and the wave velocity.

Before proceeding, let us explain why we need these predictive tools and why we cannot rely solely on numerical solutions of the nonlinear model equations. The primary reason is that these tools give us insight into the factors that shape the evolution of waves. The second reason is pragmatic. A typical simulation takes about 24 hrs of computation time, so we rely on our predictive tools to guide us in making judicious parameter choices. The input data for the probabilistic equations require only about 1 hr to compute.

Basic steps in propagating waves on discrete rectangular lattices

Figure 4 shows how a Ca²⁺ wave propagates on a discrete rectangular lattice. The 11 images are of the Ca²⁺ concentration distribution in the x-y plane at t = 0, 26, 29, 34, 38, 55, 65, 75, 85, 95, and 105 ms (left to right, top to bottom). The white dots are the locations of the CRUs spaced 2 µm apart longitudinally and 0.8 µm apart transversely. The white curves show where C is between 10 µM and 20 µM. The origin is at the lower left corner. At t = 0, C(x, y, 0) = 0.1 µM except on the 4 × 4 µm² square at the origin where C = 30 µM. Note that the buffers are in equilibrium with the Ca²⁺ everywhere. We pick up the progress of the wave at t = 26 ms where we see five firing CRUs at x = 6 µm; the CRU at y = 2.4 µm fired at t = 15 ms and is no longer releasing Ca²⁺. At this time there is a "wall" of high Ca²⁺ concentration (white curve) bearing down on the CRUs at x = 8 µm. The next three images at t = 29, 34, and 38 ms illustrate the rapid activation of six contiguous CRUs along x = 8 µm. This rapid activation is shown also in the 12th image, which is a transverse linescan through x = 8 µm. This linescan image is oriented with the y-coordinate in the same direction as the x-y images, so there is a point-to-point correspondence with the x-y image. The six contiguous CRUs fire in rapid succession, giving rise to the nearly flat apex of the transverse wave.

View larger version (81K): [in this window] [in a new window]   FIGURE 4   Wave propagation on rectangular lattices. These snapshots of Ca2+ distribution in the x-y plane were taken at 0, 26, 29, 34, 38, 55, 65, 75, 85, 95, and 105 ms (left to right, top to bottom). The white dots are the locations of the CRUs spaced 2 µm apart longitudinally and 0.8 µm apart transversely. The white curves show where 10 µM < C < 20 µM. The simulation starts with C set to 30 µM on 0 < x < 4 µm, 0 < y < 4 µm (image 1); buffers (endogenous and fluo-3) are in equilibrium with Ca2+. Elsewhere, all chemical species are set to their baseline values. Images 2-5 illustrate how Ca2+ waves propagate on rectangular lattices. Image 2 shows a "wall" of high Ca2+ concentration (white curve) bearing down on the column of CRUs at x = 8 µm. This wall of Ca2+ causes the nearly simultaneous activation of about 7 CRUs (images 3-5). The sequence of firings along x = 8 is shown in the transverse linescan (image 12, time flows from left to right). The near simultaneous firing of ~6 CRUs shows up as the flat part of the linescan image. Subsequently, CRUs along y fire sequentially at nearly regular intervals, seen as the sloping part of the linescan image. About 5-7 contiguous CRUs are needed to raise the Ca2+ concentration high enough to trigger CRUs 2 µm away at the next column to sustain the wave as was seen in the previous figure. Note spontaneous sparks occur in images 6, 10, and 11. "Jumping" of the wave from one column to the next is not evident in these static images; this solution and others in the form of MPEG movies can be downloaded from our ftp site ftp://ntcv.umaryland.edu/pub/izu/. (The readme.txt file gives a guide to the movies.) The scale bar in panel 12 represents 50 ms.

In the 5th image (t = 38 ms), there is the wall of high Ca²⁺ concentration approaching the CRUs at x = 10 µm. In time, CRUs at x = 10 µm will fire, sending out another wall of Ca²⁺ to x = 12 µm, thus regenerating the cycle of wave propagation.

The sequence of events just describedfiring of 5-7 contiguous CRUs in rapid succession, generating a wall of high Ca²⁺ that, in turn, triggers CRU firing at the next sarcomereis the basic pattern for Ca²⁺ wave propagation on lattices. To generate enough Ca²⁺ to saturate the buffers and to sufficiently raise the Ca²⁺ concentration at the next sarcomere, about 5-7 CRUs must fire in rapid succession as was seen in Fig. 3 C, illustrating superadditivity.

The rapid firing of 5-7 contiguous CRUs on a z-line shows up as the flat apex in a transverse linescan image (image 12, left-hand side). However, after this initial flurry of firings, subsequent firings on that z-line occur in a step-by-step fashion at regular intervals as shown in the sloping edge of the linescan image.

Images 6-11 are taken at 10-ms intervals and illustrate the steady progression of the wave. Two spontaneous sparks are seen in image 6 (t = 55 ms) at (10, 9.6) (near the center) and at (18, 14.4). The elliptical diffusion of Ca²⁺ is clearly seen in the (10, 9.6) spark. These sparks do not trigger firing from adjacent sarcomeres 2 µm away. In the 10th (t = 95 ms) image image a spark occurs at (18, 13.6) and, in the 11th (t = 105 ms) image, sparks occur at (18, 15.2) and (18, 12.8). The sparks at (18, 13.6), (18, 15.2), and (18, 12.8) were probably triggered by the spark at (18, 14.4). Using the probabilistic equations in the next section, we calculate the probability of these CRUs firing to be ~0.5. In contrast, the probability of firing of a CRU far from the (18, 13.6) spark, say at (12, 14.4), during 55 < t < 105 is only 0.014. Animated images of the simulation shown in Fig. 4 and other simulations, in the form of MPEG movies, can be downloaded from our ftp site ftp://ntcv.umaryland.edu/pub/izu/.

Wave statistics

Earlier deterministic models for Ca²⁺ wave propagation (Backx et al., 1989; Keizer et al., 1998; Lukyanenko et al., 1999) assumed that the CRU fires as soon as the local Ca²⁺ concentration exceeds the threshold C*. In a stochastic system, this would be true if n were so large that P(C, K, n) would be nearly a step function with C*  K. Because n ~ 2, however, P(C, K, n) has a shallow slope, so there is a distribution of times before a CRU fires. Thus exactly when a CRU will fire cannot be known, so, instead, we calculate the probability of firing within some time interval.

If the Ca²⁺ concentration at a CRU is fixed at , then the waiting time before it fires follows an exponential distribution and the mean waiting time (MWT) before firing is 1/P(, K, n). This standard result is not applicable to studying waves because the Ca²⁺ concentration at a CRU varies with time. In Appendix B, we derive, assuming time-varying P(C, K, n), expressions for p(x, y, t < ), the probability that a CRU at (x, y) does not fire in time t <  (Eq. B3); (x, y, t), the waiting time distribution (Eqs. B5 and B8); the MWT before the CRU fires (Eq. B6), and P(X  1, t < ), the probability that at least one CRU has fired in time t <  (Eq. B7).

Effect of I_SR and T_open on waves

We apply these formulas to study the effects I_SR and T_open have on waves. First, we approximate the Ca²⁺ distribution in a wave whose front is at column i = 0 by two methods. In method a we simultaneously fire two columns of seven contiguous CRUs (y = _yj, j = 3, ... , 3, x = _x and x = 0); in method b we fire one column of seven CRUs (y = _yj, j = 3, ... , x = 0) and set the initial Ca²⁺ (and the other chemical species) to a high value, comparable to that seen in waves, in the region x < 0, 3_y < y < 3_y. The choice of seven CRUs comes from observing, in Ca²⁺ waves, that ~5-7 CRUs must fire on column i before any CRU on column i + 1 fires. In method b, the choice of the initial Ca²⁺ concentration is arbitrary, but the range is guided by results from some simulations. The distribution in a was chosen to approximate the distribution in a nascent wave, whereas that in b, a well-developed wave. Needless to say, both methods approximate the concentrations in a wave very crudely, but they nevertheless provide useful guides in choosing model parameters for simulations. The time cost of these rough calculations, ~1-2 hrs, is well warranted before carrying out an ~24-36-hr simulation.

Method a underestimates the velocity because the Ca²⁺ concentration gradient on x < 0 is greater than that in a wave, so less Ca²⁺ flows to the CRUs on column i + 1. Method b overestimates the velocity because, in a wave, the CRUs on the wave front fire in rapid succession but not simultaneously as assumed in method b. Thus, the velocity from the simulation is bracketed by the estimates from the two methods.

Each CRU at (x = 2, y = 0.8j), j = 3, 2, ... , 2, 3 sees a changing Ca²⁺ concentration similar to that in curve d in Fig. 3 C. Knowing C(x = 2, y = 0.8j, t), we then calculate the probability that at least one CRU at x = 2 will fire within  ms, P(X  1, t > ). The waiting time distribution (WTD)  before a CRU at x = 2 fires is then used to calculate the MWT for firing.

Figure 5 A shows the probability of at least one firing and 5 B the waiting time distributions for different combinations of (I_SR, T_open). The curves labeled with circles, squares, and down-triangles derive from I_SR = 20 pA and T_open = 8, 7, and 5 ms, respectively. The curve with up-triangles comes from (15, 8). The a and b sublabels signify the method used to generate the concentration distributions. These curves illustrate the striking sensitivity of P(X  1, t < ) and the WTDs to I_SR and T_open.

View larger version (22K): [in this window] [in a new window]   FIGURE 5   Probability of firing (Aa, Ab) and waiting time distributions (Ba, Bb) for different combinations of (ISR, Topen). The sublabels (a,b) refer to the method the Ca2+ concentration in the wave was approximated (see text). Circles, (ISR, Topen) = (20, 8); squares, (20, 7); up-triangle, (15, 8); down-triangle, (20, 5).

First consider the case where I_SR = 20 pA and T_open = 8 ms (circles). Using the distribution of method a, we see (Aa) the probability of firing starts rising starting ~10 ms after the channels at x = 0 and x = 2 have fired, and, by 20 ms, it is almost certain that at least one CRU at x = 2 has fired. The WTD  (Ba) also rises rapidly, peaks at 15 ms, then declines rapidly. The decline of  to almost zero beyond ~25 ms means that the probability of needing to wait more than 25 ms before a CRU fires is virtually zero. The MWT for this distribution is 15.89 ± 4.47 ms, from which we estimate the longitudinal velocity by _x = _x/MWT = 2 µm/15.89 ms = 126 µm/s.

The parallel calculations of P(X  1, t < ) and  using the distribution of method b are given in Ab and Bb. The CRUs at x = 2 see a higher Ca²⁺ concentration than in method a because the shallower concentration gradient on x < 0 causes the large amount of Ca²⁺ in the reservoir (x < 0, 3_y < y <3_y) and Ca²⁺ released from the CRUs to flow mostly rightward. Not surprisingly then, P(X  1, t < ) approaches 1 more quickly, and  is sharper for method b. The MWT is now 10.06 ± 1.91 ms and the estimated velocity is 199 µm/s. The wave velocity measured in a simulated wave (see below) with these parameters was 155 µm/s, which is between the estimates from methods a and b.

Firing properties change dramatically when T_open decreases by just 3 ms from 8 ms to 5 ms (down triangles). Whereas the probability of CRU firing at x = 2 was almost 1 at t = 20 ms for T_open = 8 ms, when T_open = 5 ms, the probability of firing does not reach 0.5 until 80 ms using method a or 37 ms using method b. Note the waiting time distributions are very broad with a barely perceptible peak at 22 ms (Bb) and 47 ms (Ab). The MWTs are 69.03 ± 35.40 ms (a) and 40.63 ± 30.22 ms (b). The broad waiting time distributions mean that a CRU at x = 2 will have an equally low chance of firing at almost any time after ~10 ms after CRUs at x = 2 and x = 0 have fired.

We do not expect waves to propagate with perfect regularity in a stochastic system, but there must be a modicum of predictability. The broadening of the waiting time distribution curves signals loss of predictability. For the parameter pairs (20, 5) and (15, 8), down- and up-triangles, the waiting time distributions are so broad that they preclude a reasonable estimate of the time a wave front will reach the next column of CRUs. Hence these calculations indicate that Ca²⁺ waves are unlikely to occur with these parameter combinations. These predictions were confirmed by simulations started with method a.

Physiological meaning of the waiting time distributions

We get a clearer understanding of the physiological meaning of the WTDs by examining them in conjunction with linescan images of Ca²⁺ waves shown in Fig. 6. The three waves were produced with parameter combinations (I_SR = 20 pA, T_open = 8 ms) (panel Aa), (20, 5) [panel Ba], and (5, 32) (panel Ca). These linescans are oriented with time flowing from left to right, and the linescan was along the longitudinal axis at y = 2.4 µm. In Panel A the wave produces a sharp edge in the linescan image whose slope (measured with respect to the horizontal) equals the wave velocity of 155 µm/s. The sharpness of the edge stems from the almost perfectly regular firing intervals between CRUs shown in the accompanying time plot of the Ca²⁺ concentration at the release sites at x = 6, 8, 10, ... , 18 µm (Fig. 6 Ab). The mean time between firings is 14.1 ms with a standard deviation (SD) of 2.1 ms. In this and remaining time plots, the Ca²⁺ concentration rises rapidly in the first few milliseconds after the CRU fires, then more slowly as diffusion balances influx. The sudden drop occurs when the CRU closes.

View larger version (34K): [in this window] [in a new window]   FIGURE 6   Regularity of wave propagation for different (ISR, Topen). Local variations in the wave velocity are reflected in the longitudinal linescan images Aa (ISR = 20, Topen = 8), Ba (20, 5), and Ca (5, 32). In Aa, the wave front progresses from z-line to z-line (bright spots showing CRU firing) at a regular pace, giving rise to a sharp edge in the linescan image. Ab is the corresponding time plot of log(C(x, y = 2.4, t)) for x = 6, 8, 10, ... , 16 µm showing the wave moving with almost deterministic precision despite the firing of CRUs being a random variable. The mean time between firings is 14.1 ms with standard deviation (SD) of only 2.1 ms. This great regularity is predicted from the very narrow waiting time distribution (WTD), curve A, panel D, calculated using an approximation of the Ca2+ distribution in a wave using method b (see text for explanation). Ba shows a ragged wave indicating very irregular wave propagation caused by reducing Topen by just 3 ms to 5 ms. The time plot of log(C(x, y = 2.4, t)) (Bb) shows not only large variability in the time of CRU firing but also the local Ca2+ concentration at the time of firing (circles). Note that the CRU at x = 12 sees a Ca2+ of 3 µM for ~30 ms but does not fire until C  10 µM. The mean time between firings is 41.6 ms with a large SD of 29.1 ms. The WTD (curve B, panel D) is also very broad; its mean is 40.6 ms and SD is 30.2 ms. The parameters for wave C were chosen so that the total Ca2+ released, ISR × Topen, was the same as in A. Although the parameters for A and C appear vastly different, the waves are more similar to each other than the wave in B. Linescans A and C are scaled between 0.1 and 350 µM and are about equally bright because the total Ca2+ released is identical; linescan B is scaled between 0.1 and 175 µM and is considerably dimmer than A and C. C's wave velocity (81 µm/s) is about half of A's (155 µm/s) reflecting the slower rate of release (ISR), but the wave front is quite sharp. The mean time between firing (from the time plot Cb) is 24.6 ms with SD = 8.4 ms; the WTD has mean of 23.3 ms and SD of 4.6 ms.

In contrast to this wave that propagates with almost deterministic precision, the wave in Fig. 6 Ba propagates unpredictably. The unpredictable timing of firing of successive CRUs is seen in the time plot (Fig. 6 Bb) and gives the linescan image a ragged edge. The stochastic nature of CRU firing is clearly evident in the time plot. Note that the Ca²⁺ concentration at which the CRU fires (circles) varies widely from 1.3 to 15.3 µM. The CRU at x = 12 µm (4th curve from left) sees a Ca²⁺ concentration of ~2 µM for ~30 ms without firing, although CRUs at x = 6, 8, and 10 µm had fired at this concentration. In fact, the CRUs at x = 12 and 16 µm do not fire until the Ca²⁺ concentration is raised by the firing of neighboring CRUs (off of the linescan) to ~10 µM.

The widely differing firing patterns of waves in Aa and Ba are reflected in their respective WTDs shown in Fig. 5 and reproduced in Fig. 6 D. The nearly fixed firing intervals in A can be anticipated from the narrow WTD, which has a SD of 1.91 ms, close to the measured dispersion of firing intervals in the wave. The mean of this WTD (i.e., the MWT) is 10.9 ms, somewhat less than the mean firing interval. As mentioned earlier, the MWT calculated using method b overestimates the wave speed.

The WTD for (20, 5) (curve B) is very broad; its mean is 40.6 ms with a SD of 30.2 ms. The large dispersion in the WTD is the reason for the large variability in firing intervals for the wave in Panel B. The MWT and SD for WTD (B) are close to the mean firing interval measured from time plot Bb of 41.6 ms with SD of 29.1 ms.

Total Ca²⁺ released is a strong determinant of wave properties

The waves in panels A, B, and C of Fig. 6 were generated with (I_SR, T_open) of (20, 8), (20, 5), and (5, 32), respectively. Although the parameters are similar for A and B, the waves are quite different. In contrast, (20, 8) and (5, 32) are dissimilar yet the waves they engender are similar in terms of their regularity. What is common to both (20, 8) and (5, 32) parameter sets is the total amount of Ca²⁺ released per CRU (I_SR × T_open) is the same. The mean and SD for the (5, 32)-WTD, curve C, are 23.33 and 4.55; the mean firing interval measured from the time plots in Fig. 6 Cb is 24.56 ms and the SD is 8.44 ms. We expect that the MWT for (5, 32) would be longer than that for (20, 8) simply because the rate of release is 4 times smaller. Surprisingly, however, the MWT is only about twice as long, not four times longer (24.6 versus 14.1). In contrast, the MWT for (20, 5) is about threefold larger (41.6 versus 14.1) than the MWT for (20, 8).

These results show that there is no simple relationship between the rate of release I_SR, the open time, and the wave velocity. However, the total Ca²⁺ released, I_SR × T_open, is a fairly good predictor of the SD of the WTD, hence a fairly good predictor of the sharpness of the wave front. Two waves having equal I_SR × T_open but different I_SR, would appear about equally sharp in linescan images but would have different velocities. In actual linescan images, we do not know I_SR × T_open, but, because it equals the total Ca²⁺ released, cells with equal I_SR × T_open would be about equally bright (compare time plots in Fig. 6, Ab and Cb). "Equally bright" must be interpreted cautiously, however, when using a nonratiometric dye such as fluo-4.

Initiation of waves

The distinction between wave propagation and wave initiation is essential for understanding the stability of the Ca²⁺ control system in muscle. For (I_SR = 20, T_open = 5) the MWT before activating CRUs 2 µm away is 69 ± 35 ms when the Ca²⁺ distribution is set up using method a (Fig. 5 Ba, down triangles). Because of the large MWT and large SD, we guessed that a wave would not be initiated with these parameters and initial conditions. This guess was, in fact, correct. We carried out a simulation using the same initial conditions as method a, firing two columns of seven CRUs that raised the Ca²⁺ concentration locally to high levels, but not enough Ca²⁺ was released to trigger CRU release on adjacent sarcomeres. However, raising the Ca²⁺ concentration to 30 µM in a 4 × 4 µm region was enough to trigger the wave shown in Fig. 6 B.

These examples illustrate the inherent stability of discretely spaced CRUs to even fairly large increases in Ca²⁺ concentration that might arise from injury or the random firing of CRUs. Large enough perturbations can trigger waves, however.

To assess whether a given perturbation, i.e., nonequilibrium Ca²⁺ distribution, can trigger a wave, we calculate the WTD and P(X  1, t <) for that perturbation. Fig. 6 D shows the WTD for CRU firing at x = 6 µm for the standard initial conditions used to start wave calculations, that is, C(x, y, 0) = 30 µM on 0 < x < 4, 0 < y < 4 (dashed line). Note, in this region, that the buffers are in equilibrium with Ca²⁺, so there is a large reservoir of Ca²⁺. This WTD has a mean of 18.02 ± 5.43 ms. Given the short MWT and low SD, we expect that this IC will initiate firing of CRUs at x = 6. Initiation does occur and, in Fig. 6 B, the CRU at x = 6 (3rd spark from the left) fires on schedule at t = 16 ms. The reason this perturbation (30 µM on a 4 × 4 µm² region) could initiate a wave but the perturbation of method a (two columns of seven firing CRUs) could not is simply related to the amount of Ca²⁺ available to diffuse. In the former case, the large Ca²⁺ reservoir can raise the Ca²⁺ concentration at x = 6 µm rapidly and sustain the high concentration; in the latter case, the Ca²⁺ released is rapidly taken up by buffers, so free Ca²⁺ concentration at x = 6 µm is low.

Initiating CRU firing on an adjacent sarcomere does not guarantee wave propagation, however. The Ca²⁺ distribution set up by method b serves to assess whether a given set of parameters can support a wave. The WTD using method b for (I_SR = 20, T_open = 5), which we have seen before, is shown again in Fig. 6 D (B). The large MWT and SD (40.6 ± 30.2 ms) for this WTD make it ambiguous whether a wave might be supported. This ambiguity is reflected in the wave itself as we see in Fig. 6 B, where there are times when the wave "stalls" and might seem to cease propagating (e.g., at x = 12 µm and 16 µm).

In summary, even if a given set of parameters could support wave propagation, not every disturbance initiates a wave. To initiate a wave, the disturbance must raise the Ca²⁺ concentration sufficiently high over a sufficiently large region. The magnitude and spatial extent of the needed perturbation depends naturally on the parameters and can be estimated from the WTD.

Effect of CRU Ca²⁺ sensitivity

Up to now, we have fixed the CRUs' Ca²⁺ sensitivity K to 15 µM and focused on the effects that I_SR and T_open have on wave behavior. K is affected by the Mg²⁺ concentration (Rousseau and Meissner, 1989), drugs such as caffeine (Rousseau and Meissner, 1989) and tetracaine (Györke et al., 1997), and SR Ca²⁺ loading (Rousseau and Meissner, 1989; Györke and Györke, 1998). Figure 7 shows the predicted (curves) and actual (open squares) velocity dependence of _x on K for I_SR = 20 pA and T_open = 7 ms. The lower and upper curves were calculated from the WTD generated using methods a and b, respectively. The upper and lower boundaries decline almost linearly with K. Similar calculations using I_SR = 10, T_open = 10 and I_SR = 20, T_open = 8 show also the linear decline of the velocity boundaries, but the slopes differ for each pair of (I_SR, T_open).

View larger version (14K): [in this window] [in a new window]   FIGURE 7   Longitudinal velocity dependence on K. Curves give the lower and upper x estimates using methods a and b, respectively. Squares give the actual velocity from wave simulations. Inset shows velocity dependence as K approaches the baseline Ca2+ concentration of 100 nM. Computation parameters: ISR = 20 pA, Topen = 7 ms.

Increasing sensitivity and spontaneous waves

Halving K from 15 µM to 7.5 µM doubles the Ca²⁺ sensitivity of the CRU and triples the baseline spark rate, P. Figure 8 shows snapshots of the Ca²⁺ distribution for this high-sensitivity condition. This simulation was started with all chemical species in equilibrium so all sparks arise spontaneously. The first three panels (left to right) were taken at t = 21, 27, and 35 ms, and the arrows point to CRUs at (2, 4), (2, 3.2), and (2, 2.4), respectively. These three contiguous CRUs in the transverse direction fire sequentially and can be seen more clearly in the transverse linescan along x = 2 shown in the 12th panel. This sequential firing of sparks in the transverse direction is similar to those observed by Parker et al. (1996) in ventricular cells. Note that not all sparks trigger firing of contiguous CRUs. For example, in panel 2, the spark at (12, 9.6) (near the center) does not trigger any other spark.

View larger version (39K): [in this window] [in a new window]   FIGURE 8   High Ca2+ sensitivity (K) causes spontaneous wave. Snapshots of Ca2+ distribution at 21, 27, 35, 45, 61, 65, 72, 80, 90, 110, and 130 ms showing the spontaneous appearance of a Ca2+ wave. K is 7.5 µM, twice the usual sensitivity, giving rise to numerous spontaneous sparks. The first three panels show the sequential firing of three CRUs at x = 2 (arrows). A transverse linescan showing the dye concentration (time flowing from left to right) at x = 2 is in panel 12. Although there are numerous spontaneous sparks, most do not trigger a wave. The birth of the wave is seen in panels 5-7. The wave starts from spontaneous sparks on adjacent sarcomeres (x = 14 and 16) occurring at about the same time at about the same y value (~8). Scale bar 20 ms and 2 µm, color bar ranges from 0.1 to 350 µM for Ca2+ and 4 to 50 µM for the Ca-bound dye. Computation parameters: ISR = 20 pA, Topen = 7 ms, K = 7.5 µM.

Because of the high spontaneous spark rate, it is not surprising to see a wave arise spontaneously. The birth of a wave is seen in panels 4-7. At t = 45 ms, the CRU at (16, 9.6) (panel 4, arrow) has fired and will sequentially trigger firing of contiguous CRUs on x = 16. By t = 61 ms (panel 5), CRUs at (16, 9.6) and (16, 8.8) have fired and turned off and the CRU at (16, 8) is firing. These sparks on x = 16 have raised the Ca²⁺ concentration at (14, 9.6) (arrow) only slightly above the baseline, so the firing of the CRU at (14, 9.6) is pure happenstance. If the CRUs on x = 16 have not fired, the fate of CRUs on x = 14 might be similar to those we have seen: CRU (14, 9.6) might be an isolated spark or it might trigger one or two other sparks.

But the firing of CRUs on x = 16 contribute Ca²⁺ to sites on x = 14 and also limit the diffusion of Ca²⁺ away from site (14, 9.6). Thus the probability of firing of adjacent CRUs at (14, 9.6 ± 0.8) is slightly higher than the CRU at (2, 3.2) following the CRU (2, 4) firing. Each additional firing of CRUs on x = 16 leads to an ever increasing likelihood of CRU firing on x = 14 and vice versa, thus stabilizing the nascent wave. After a certain number of CRUs have fired, the probability of the nascent wave dying is virtually zero. This is the probabilistic equivalent of the "critical curvature" a wave must have before it can propagate (Wussling et al., 1997).

Effect of transverse lattice spacing

What effect would reducing the transverse CRU spacing from _y = 0.8 to 0.4 µm, keeping _x fixed to 2 µm, have on the ratio of _x to _y? When _y = 0.8 µm, _x/_y  1.3; for the wave in Fig. 4, _x/_y = 155 µm/s/124 µm/s = 1.25. When _y is reduced to 0.4 µm, the distance to the adjacent CRU in the y-direction is only of that to the CRU on the next sarcomere. Thus the shorter distance should compensate for the 2:1 diffusion anisotropy, and we might expect _x/_y to be less than unity. This expectation is wrong.

Figure 9 shows snapshots of a Ca²⁺ wave propagating on a lattice with _x = 2 µm and _y = 0.4 µm at 37-ms time intervals starting at t = 0 ms. As in Fig. 4, the white dots mark the CRU locations and the white curve shows where 10 µM < C < 20 µM. Contrary to our naive expectation, the longitudinal velocity is almost twice as large as the transverse velocity, 120 µm/s versus 67 µm/s. In this simulation I_SR = 15, T_open = 6, and K = 15 µM. Simulations using (I_SR, T_open) = (10, 10) or (10, 16) gave _x/_y ratios of 1.7. Thus the change in ellipticity of the Ca²⁺ wave went opposite to our expectation.

View larger version (66K): [in this window] [in a new window]   FIGURE 9   Effect of transverse lattice spacing on wave propagation. The transverse spacing of CRUs has been reduced from 0.8 to 0.4 µm. The naive expectation that the wave will travel faster transversely than longitudinally when y is reduced is not borne out as shown in these Ca2+ concentration images. The longitudinal velocity is about twice as large as the transverse velocity