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

Comparison between the Predictions of Diffusion-Reaction Models and Localized Ca²⁺ Transients in Amphibian Skeletal Muscle Fibers

Department of Physiology, UCLA School of Medicine, Los Angeles, California

Correspondence: Address reprint requests to Julio Vergara, Dept. of Physiology, UCLA School of Medicine, 10833 Le Conte Avenue 53-263 CHS, Los Angeles, CA 90095-1751. Tel.: 310-825-9307; Fax: 310-206-3788; E-mail: jvergara{at}mednet.ucla.edu.

	ABSTRACT

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

 
We developed a three-dimensional cylindrical diffusion-reaction model of a single amphibian myofibril in which Ca²⁺ release occurred only at the Z-line. The model incorporated diffusion of Ca²⁺, Mg²⁺, and all relevant buffer species, as well as the kinetic binding reactions between the buffers and appropriate ions. Model data was blurred according to a Gaussian approximation of the point spread function of the microscope and directly compared with experimental data obtained using the confocal spot methodology. The flux parameters were adjusted until the simulated Z-line transient matched the experimental one. This model could not simultaneously predict key parameters of the experimental M- and Z-line transients, even when model parameters were adjusted to unreasonably extreme values. Even though the model was accurate in predicting the Z-line transient under conditions of high [EGTA], it predicted a significantly narrower Ca²⁺ domain than observed experimentally. We modified the model to incorporate a broader band of release centered at the Z-line. This extended release model was superior both in simultaneously predicting critical features of the Z- and M-line transients as well as the domain profile under conditions of high [EGTA]. We conclude that a model of release occurring exclusively at the Z-line cannot explain our experimental data and suggest that Ca²⁺ may be released from a broader region of the sarcoplasmic reticulum than just the T-tubule-sarcoplasmic reticulum junction.

	INTRODUCTION

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

 
The development of tension in skeletal muscle depends on the binding of Ca²⁺ to the troponin complex which allows the thick and thin filaments to interact. Since this complex is located in a specific subregion of the sarcomere, understanding the spatiotemporal profile of the calcium concentration ([Ca²⁺]) changes within the muscle fiber is crucial to understanding the excitation-contraction (EC) coupling process.

Traditionally, [Ca²⁺] changes in response to physiological action potential (AP) stimulation have been measured in muscle fibers using techniques that report Ca²⁺-dependent fluorescence signals from regions that span many sarcomeres (global Ca²⁺ transients; Kim and Vergara, 1998; Vergara et al., 1991); see Rios and Pizarro (1991) for a review. However, these Ca²⁺ transients report a mix of fluorescence contributions from regions of the sarcomere that have marked differences in [Ca²⁺], resulting in a spatially averaged signal that provides no information about the spatial distribution of the [Ca²⁺] changes (Escobar et al., 1994). The most relevant information about the Ca²⁺ signaling during EC coupling can only be obtained from measurements of intrasarcomeric Ca²⁺ gradients, which require localized (confocal) detection of fluorescence signals from small regions within the sarcomere. To obtain this information, we developed the confocal spot technique (Escobar et al., 1994; Vergara et al., 2001; DiFranco et al., 2002), which allows us to measure the [Ca²⁺] in a 0.3-µm confocal spot in a single myofibril. An AP is elicited and a localized fluorescence transient is recorded from a single spot position. A complete high resolution spatiotemporal profile of the [Ca²⁺] distribution within a single sarcomere is obtained by recording several transients at different spot positions located at 200-nm intervals along the fiber axis (DiFranco et al., 1999; Vergara et al., 2001).

The magnitude and time-course of the [Ca²⁺] gradients depend on many factors:

1. The magnitude and time-course of the Ca²⁺ release flux from the sarcoplasmic reticulum (SR),
   
2. The geometry and extent of the Ca²⁺ release sites,
   
3. The interaction of Ca²⁺ with Ca²⁺ buffers, and
   
4. The magnitude and time-course of the removal of Ca²⁺ by the SR Ca²⁺ ATPase.
   

To evaluate the role of buffers and the magnitude and location of the SR Ca²⁺ release in light of data obtained using confocal detection techniques, we present a three-dimensional cylindrical diffusion-reaction model that is an extension of previous myofibrillar models (Cannell and Allen, 1984; Baylor and Hollingworth, 1998). In our model, Ca²⁺ ions were released into the myofibril from its periphery (outermost 16 nm) either only at the Z-line, which represents the T-SR junction at this location, or from an extended release site. Ca²⁺ and Mg²⁺ ions were allowed to diffuse throughout the myofibril. All relevant mobile buffers—EGTA, Oregon Green 488 BAPTA-5N (OGB-5N), parvalbumin, and ATP—interacted with Ca²⁺ (and Mg²⁺ when relevant) according to their kinetic binding characteristics and were also allowed to diffuse. The [Ca²⁺]-dependent fluorescence of the indicator OGB-5N was calculated according to its properties and the resulting fluorescence was blurred by convolving it with a triple Gaussian approximation of the point spread function (PSF) of the microscope objective. Simulated fluorescence domains were then constructed by assuming that there was a pinholed detector in the image plane of the microscope that measured a localized transient every 200 nm. These simulated domains were then directly compared with experimental ones obtained from fibers stretched to 4-µm sarcomere spacing to determine how well the model was able to predict the experimental observations.

We have found that the currently accepted view of Ca²⁺ release occurring only at the Z-line and dispersing throughout the myofibril solely by diffusion is not compatible with measurements obtained using the confocal spot technique. Several key parameters could not be predicted, including the time-course and magnitude of the M-line transient, the delay before fluorescence that is first observed at the M-line, and the width of the domain observed under conditions of high [EGTA]. By extending the release site most of these discrepancies were minimized.

Parts of the model comparisons presented in this article have previously been published in abstract form (Novo et al., 2001, 2002).

	METHODS

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

 
Confocal spot detection
Experimental traces used to compare to model simulations were obtained using a confocal spot detection methodology, described in detail previously (Escobar et al., 1994; Vergara et al., 2001; DiFranco et al., 2002). Briefly, a single muscle fiber from the dorsal head of the semitendinosus muscle of the frog Rana catesbeiana was mounted in a double Vaseline gap chamber situated on a custom made stage scanning inverted confocal microscope (DiFranco et al., 1999; Vergara et al., 2001; DiFranco et al., 2002). Both ends of the muscle fiber were soaked in an internal solution containing the Ca²⁺ indicator dye OGB-5N (Molecular Probes, Eugene, OR). In experiments with high (either 10 or 20 mM) [EGTA], the fiber was first loaded with 0.5 mM EGTA and transients were acquired as described in DiFranco et al. (1999, 2002) and Vergara et al. (2001). The high [EGTA] was then added to the cut ends along with Ca²⁺ such that the free [Ca²⁺] was 0.1 µM, and allowed to equilibrate for 30 min before transients were recorded.

Given that the dimensions of the confocal spot (see below) are smaller than the diameter of a frog myofibril, it is reasonable to assume that our fluorescence records arise from a single myofibril (DiFranco et al., 2002). Thus, our model, which is designed to generate data that can be quantitatively compared to the experimental data, is an extension of the single myofibril cylindrical scheme proposed by Cannell and Allen in 1984. In our model, Ca²⁺ ions are either released only at the Z-line or from an extended release band and then diffuse throughout the myofibril, all the while interacting with the various Ca²⁺ binding species. Based upon the formation of Ca²⁺-dye complex at different radial and longitudinal positions, localized fluorescence transients are calculated after considering the blurring introduced by the point-spread function (PSF) of the optical system (see below). When these signals are expressed in terms of F/F (see below), they can be directly compared with the experimental data.

Model geometry
The half-sarcomere of a frog single myofibril is represented as a radially symmetrical cylinder with the positional variables being r, the radial coordinate, ranging from 0 to a from the center to the edge of the myofibril, and x, the longitudinal coordinate, ranging from 0 to l from the Z-line to the M-line (Fig. 1 A). Because of the symmetry within a single myofibril it is assumed that there is no net longitudinal flux of Ca²⁺ at the M-line (x = l). In addition, because of symmetry with other myofibrils, it is assumed that there is no net flux of Ca²⁺ in or out of the myofibril except at the Ca²⁺ release sites and SR Ca²⁺ pump (Fig. 1 B). It is also assumed that the myofibril is an isotropic volume in which all molecules can freely diffuse, with a mobility determined by their diffusion coefficients, in any direction (subject to the boundaries of the model).

View larger version (72K): [in this window] [in a new window]   FIGURE 1  Schematic diagrams of simulated myofibril. (A) Model geometry of the cylindrical representation of a half-sarcomere. Ca2+ is released at the periphery of the myofibril either at the Z-line (dark ring at x = 0, r = a), or in an extended release region (darker shaded area). Dotted arrow denotes that the extended release region can be of a variable length. A typical simulation used 32 radial (r) and 128 longitudinal (l) segments. (B) Voxel types present in the myofibrillar model. The numbers indicate which diffusion-reaction equation from Table 3 applies to a given voxel position. In the Z-line release model, Ca2+ enters the myofibril through channels located on the boundaries between voxels 1 and 2 and between 1 and 4 (solid arrows). In the extended release model, Ca2+ enters the myofibril from voxel 1 and a variable number of voxels 2 (hatched arrows). Ca2+ exits via the SR ATPase located in voxels 2 and 3 (unfilled arrows). The dots represent repeats of the voxel type in the longitudinal (horizontal dots) or radial (vertical dots) directions. Typical values for m and n are 128 and 32, respectively. a represents the radius and l represents the length of the half-sarcomere, 0.5 and 2 µm, respectively.

Although either model could implement any Ca²⁺-release driving function, to simulate the release process in response to an AP, we used an exponential rising and falling function described by the equation

(1.1)

The specific values for J_max, _on, and _off were obtained by varying them until the time course of the simulated fluorescence transient at the Z-line matched a typical experimental Z-line transient. J(t) is in units of pmoles/(cm² x ms), in which the cm² refers to the cross-sectional area across which the flux is occurring. The entire magnitude of the flux applies to the half-sarcomere. In the extended release model, the flux (Eq. 1.1) was identical in all release voxels. To attain this function, the [Ca²⁺] increment at each time-step was calculated by considering the flux, outer surface area, and volume of each release voxel. When the extended flux simulations were constrained to release Ca²⁺ only from voxel 1 in Fig. 1 B, the flux required to simulate the experimental Z-line transients was identical to the Z-line release model (data not shown).

Diffusion reaction scheme
The general diffusion-reaction equation that describes the rate of change of reactant concentration at any point is

(1.2)

where c_p is the concentration of a reactant p at a given coordinate (r, x), represents the diffusional component of the rate of change of this reactant's concentration, and W_p is the binding reaction with other molecular species and the transport via the SR (if applicable). The model involves the simultaneous integration of multiple equations of the type of Eq. 1.2 for Ca²⁺, Mg²⁺, and the bound and free forms of each buffer in the system.

Diffusion
In the first term of Eq. 1.2, D_p represents the diffusion coefficient of the p^th molecular species and ² is the Laplace operator (Crank, 1975). After expressing the ² operator in cylindrical coordinates, Eq. 1.2 becomes

(1.3)

where r is the radial variable and x is the longitudinal variable.

To simplify the calculations, one can make the following transformations into dimensionless variables:

(1.4)

(1.5)

(1.6)

Substituting Eqs. 1.4– 1.6 in Eq. 1.3 yields the simplified equation (Crank, 1975),

(1.7)

Henceforth, c₁ and c₂ will represent the free [Ca²⁺] and [Mg²⁺], respectively.

Reaction term for Ca²⁺ buffers
The model supports an arbitrary number of Ca²⁺ buffers. For simplicity, we have considered that each buffer only has a single binding site for which Ca²⁺ and Mg²⁺ can compete. Thus, there are up to three distinct forms for each buffer (free, Ca²⁺-bound, and Mg²⁺-bound) whose concentrations are considered a distinct c_p obeying Eq. 1.7. Therefore a buffer is completely defined by its diffusion coefficient, kinetic rate constants for binding Ca²⁺ and Mg²⁺, initial concentration, and initial distribution (i.e., the buffer could be nonhomogenously distributed within the myofibril). The Ca²⁺-bound, Mg²⁺-bound, and free form of the buffer all diffused independently of each other. At the start of the simulation the buffer was assumed to be at kinetic equilibrium with the free Ca²⁺ and Mg²⁺. If the buffer only bound Ca²⁺, the initial concentration of Ca²⁺-bound buffer (e.g., c_3c for OGB-5N) is given by

(1.8)

where c_3tot is the total dye concentration, is the equilibrium dissociation constant for the reaction of the dye with Ca²⁺, and c₁(0) is the [Ca²⁺] at the start of the simulation (0.1 µM).

If the buffer binds both Ca²⁺ and Mg²⁺, the initial concentration of the buffer species can be calculated, considering that

and are given by

(1.9)

(1.10)

The notation for Eqs. 1.9 and 1.10 is presented in Appendix A.

The sole component of the W_p term in Eq. 1.7 for buffers is their binding to Ca²⁺ and Mg²⁺. The rate of change of free Ca²⁺ as a result of interactions with a Ca²⁺ buffer (e.g., from OGB-5N), is described by

(1.11)

where and are the kinetic association and dissociation rate constants, respectively, for the interaction of OGB-5N with Ca²⁺. Conversely, the reaction term for the bound form of the buffer (e.g., c_3c for OGB-5N) is equal to

Reaction term for Mg²⁺
The reaction term for Mg²⁺ represents the change in [Mg²⁺] due to its interaction with Mg²⁺-binding buffers. The only two buffers for which this was considered were parvalbumin and ATP. The change in free [Mg²⁺] (e.g., from ATP) is described by:

(1.12)

In the case of a buffer that bound both Ca²⁺ and Mg²⁺, the change in free buffer was equal to the sum of Eqs. 1.11 and 1.12.

Reaction term for Ca²⁺
Since Ca²⁺ interacts with all the buffers in the system as well as with the SR Ca²⁺-ATPase, the W₁ term is the sum of two components. The first, Q₁, represents the change in free [Ca²⁺] due to binding to buffers. This term is equal to the sum of all the Ca²⁺ components of the reaction terms for the free buffers of the form of Eq. 1.11. Thus,

(1.13)

The second component of W₁ is the Ca²⁺-ATPase (P₁), located at the SR membrane, that pumps Ca²⁺ against its concentration gradient from the cytoplasm to the lumen of the SR (Ebashi and Endo, 1968). In the model, the SR-ATPase was located in voxels denoted by 2 and 3 in Fig. 1 B. Although it is thought that the SR-ATPase may have complex kinetics (Burmeister-Getz and Lehman, 1997), its operation was approximated by the equation:

(1.14)

where P₁ is the rate of change of [Ca²⁺] due to pumping, V_max is the maximal reaction velocity, and is the Michaelis-Menten constant. The calculation of V_max involves estimating the SR surface area in a typical myofibril (Table 1), which we approximated from the literature ratio of SR surface area to fiber volume (Eisenberg, 1983). Thus,

(1.15)

Then, multiplying Eq. 1.15 by the estimated surface density of pump molecules, 30,000 pumps/µm² of SR, (Franzini-Armstrong and Ferguson, 1985) and the pumping rate per molecule, (0.0076 Ca²⁺ ions/ms) (Vilsen and Andersen, 1992), then dividing by the volume of the peripheral ring in which the SR-ATPase resides (9.17 x 10^-17 L), yields a V_max of 9.74 µM/ms.

Boundary and initial conditions
The diffusion reaction equation from Eq. 1.7 was integrated to conform to the boundary conditions listed in Eq. 1.16. The model assumed that there was no net Ca²⁺ flux across the boundaries, except at the Ca²⁺ release sites. Thus,

(1.16)

For the Z-line release model an additional boundary condition was

For the extended release model,

but Ca²⁺ ions were released into the model at each time-step resulting in a flux described by J(t) at r = a and 0 < x < = , where represents the length of the release band.

The fact that the only flux across the r = a boundary occurs at the release sites implicitly assumes that the surrounding myofibrils are perfectly in register with the simulated one. This assumption is supported by experimental evidence showing that there is no di-8-ANEPPS fluorescence from the T-tubules anywhere in the sarcomere except at the Z-line (DiFranco et al., 2002; Vergara et al., 2001). The boundary conditions for Mg²⁺ and all the other buffers were identical, except that they did not have the source. Thus,

(1.17)

The initial conditions for Ca²⁺ and Mg²⁺ were set as follows:

(1.18)

The initial conditions for the buffers were calculated according to Eqs. 1.8–1.10.

Integration method
Since the Crank-Nicolson integration method for parabolic partial differential equations has excellent stability and convergence and allows for relatively large time-steps (Crank, 1975), we attempted its use to solve the current model equations. However, its advantage was lost when integrating in two dimensions due to the necessity of inverting large matrices for which no efficient algorithms exist. The alternating direction implicit (ADI) (Peaceman and Rachford, 1955) method allows for large time increments without the need for inverting large two-dimensional matrices and has been used to solve diffusion-reaction models in other contexts (Chiu and Hoppensteadt, 2000; Chiu and Walkington, 1997). A complete description of the implementation of the ADI method in our specific case is given in Appendix B.

Fluorescence
As stated above, the Ca²⁺ indicator OGB-5N diffused and interacted with Ca²⁺ identically to any other buffer. Fluorescence changes were calculated in terms of F/F as follows:

(1.19)

where F_max/min is the fluorescence ratio between saturated and free dye, and c_3c(0) is the Ca²⁺-dye complex at the start of the simulation, calculated as per Eq. 1.8.

Comparison with of model predictions with experimental data
To compare F/F values generated by the model with experimental data we considered that the experimental data was blurred due to the microscope objective (Agard et al., 1989; Castleman, 1979; Wilson, 1990). Thus, model fluorescence data was blurred by convolving it with the triple Gaussian approximation of the PSF of the objective. The PSF could be calculated from confocal theory (Castleman, 1979; Wilson, 1990); however, we chose to measure it directly to take into consideration the optical properties of the specific objective used (Nikon, Fluor-100 1.3 NA).

The lateral (x-y) and axial (z) resolution were determined in vitro using 0.1- to 2.0-µm diameter fluorescently labeled latex beads (Fluorospheres, size kit #2, Molecular Probes, Eugene, OR). The beads were attached to a glass coverslip and bathed in water or embedded in 40% glucose and 0.5% agarose to mimic the refractive index of the fiber (Huxley and Niedergerke, 1958). Beads were scanned in the x- and z-directions in 100- and 500-nm steps, respectively. Fluorescence intensity as a function of position was fit to the Gaussian function,

(1.20)

where x_c is the center position, A is the amplitude and full-width at half-maximum (FWHM) is the width at A/2. For comparison between beads of different diameters, fluorescence values for each bead were normalized to its peak. The PSF was considered to be the normalized Cartesian triple Gaussian:

(1.21)

in which _xy and _z are the standard deviations of the Gaussian fits to the smallest resolvable bead. The FWHM of Gaussian fits to the fluorescence profiles of beads smaller than the diffraction limit were 0.35 and 0.75 µm in the x- and z-directions, respectively. This corresponded to a _xy of 0.21 and a _z of 0.45.

To convolve the Cartesian PSF with the model-generated F/F profile at a given time (F/F snapshot), a three-dimensional Cartesian snapshot of the fiber was generated from the cylindrical model. The convolution was performed either in the Fourier domain, using a custom-made program implementing a three-dimensional FFT algorithm (Press et al., 1996) or in the space domain using an algorithm that directly implemented a convolution (Novo and Vergara, unpublished results). Both methods yielded identical results, although the Fourier method was significantly faster. Experimentally measured bead fluorescence profiles closely approximated blurred model bead simulations (data not shown).

To compare model data with the experimentally obtained localized fluorescence transients, we took into account that in our experimental setup the detection pinhole would only accept light from a 0.5-µm disc region of the myofibril (Vergara et al., 2001; DiFranco et al., 2002). As an approximation, the simulated pinhole was considered to be a 0.5-µm square located at the (x, y, 0) plane of the myofibril. To obtain the fluorescence at a particular x-position from 0–l, the square was centered at the edge of the myofibril (y = a), and at that x-coordinate. Blurred fluorescence values from voxels that were within the pinhole were summed to yield a fluorescence value for that x-position.

Transients were characterized by several kinetic properties (DiFranco et al., 2002):

1. The delay time (t_d, in ms) corresponds to the period between the time of stimulus delivery (t = 0) and the initiation of the rising phase of the transient, defined as the first moment when the fluorescence is sustained for 0.2 ms at two standard deviations (SDs) above F_rest. For simulations, which were noiseless, the initiation of the transient was considered to occur when the simulated F/F rose two SDs above the noise of the corresponding experimental transient.
   
2. The rising phase of the transient was characterized by the rise time (rt, in ms), defined as the interval between the initiation and peak of the transient.
   
3. The derivative at any given time point was calculated as the average slope of the lines joining the given point with the preceding and following point. The maximal derivative will be denoted (d(F/F)/dt)_max (in ms^-1).
   
4. The peak F/F (F/F_peak) is the maximum F/F attained during the transient.
   
5. The full-duration at half-maximum (FDHM) is the length of time during which the fluorescence was greater than one-half of the F/F_peak.
   

	RESULTS

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

 
Before we describe realistic model simulations that can be directly compared with experimental data, it is useful to examine simple Z-line release simulations to illustrate several aspects of intrasarcomeric [Ca²⁺] gradients. Fig. 2 depicts the spatial distribution of the [Ca²⁺] at the peak of an arbitrary release flux in a fiber with no buffers (Fig. 2 A) or containing 0.5 mM OGB-5N (Fig. 2, B–D). The Z-line is at longitudinal position 0 and the two M-lines symmetrically flank the Z-line at longitudinal positions ±2 µm. In this model, Ca²⁺ release occurs at the outer boundary of the myofibril at the Z-line (radial position = 0.5 µm). Fig. 2 A shows that at the peak of the release, the [Ca²⁺] at the release site has increased to 155 µM with a steep gradient that falls toward the M-line, where the [Ca²⁺] is 10 µM. An important observation is that the gradient in the radial direction is not nearly as steep as the gradient in the longitudinal direction. The [Ca²⁺] at the center of the myofibril, 0.5 µm from the source radially, is 90 µM, whereas the concentration 0.5 µm longitudinally from the source is 68 µM. This phenomenon is seen irrespective of the simulation conditions and reflects a general property of the cylindrical geometry due to the fact that Ca²⁺ is being released in a ring around the myofibril and converges toward its central axis. Analytical solutions of the diffusion equation for a point source Ca²⁺ release predict that the [Ca²⁺] in the immediate vicinity of the source is inversely proportional to the diffusion coefficient for Ca²⁺ ions (D₁) (Stern, 1992; Pape et al., 1995). In our model, the peak [Ca²⁺] in the voxel immediately adjacent to the Z-line changed from 102 to 242 µM when D₁ was changed from 4 x 10^-6 to 1 x 10^-6 cm²/s, illustrating that the voxel dimension is small enough to give results that are compatible with analytical predictions.

View larger version (65K): [in this window] [in a new window]   FIGURE 2  Effects of OGB-5N and microscope blurring on the fluorescence gradients. A Z-line release simulation was conducted in a standard myofibril (Table 2) with an arbitrary flux on = 4 ms, off = 3 ms). (A) Spatial profile of free [Ca2+] at the peak of the release (2.2 ms after the initiation of Ca2+ release) in a myofibril devoid of buffers. (B) Spatial free [Ca2+] profile at the peak of the release in an identical simulation as A except that 500 µM OGB-5N was included. (C) Spatial profile of OGB-5N bound to Ca2+ at the peak of the release. (D) Spatial F/F profile at the peak of the release blurred according to the PSF of the microscope objective.

Comparison of model predictions with experimental data
Whereas model simulations allow us to predict the spatiotemporal characteristics of [Ca²⁺] profiles in both the radial and longitudinal directions of the myofibril, the spot detection methodology allows only for the acquisition of fluorescence signals at various longitudinal positions of the myofibril that are blurred by the PSF of the objective. As explained above (Fig. 2 D), the major effect of blurring is that it collapses the radial dependence of the fluorescence gradients. Thus, by recording fluorescence signals from several longitudinal positions and plotting them against the distance from the Z-line at which they were acquired, we obtain the spatiotemporal profile of the "radially-averaged" F/F changes in a myofibril (Ca²⁺ domain) (DiFranco et al., 2002; Vergara et al., 2001). A typical experimental domain is shown in Fig. 3 A. To allow comparison with the experimental data, model domains were constructed from simulated data as described in Methods. A F/F domain from a Z-line release simulation using identical parameters as in Fig. 2, B–D, is shown in Fig. 3 B. It is possible to see that F/F gradients develop early in time and that by 20 ms the gradients dissipate with an elevated fluorescence throughout the myofibril. Fig. 3 C is a domain from a model simulation using more realistic conditions. Although the release flux was identical to that used in Fig. 3 B, endogenous Ca²⁺ buffers (Table 2) and the SR-ATPase (Table 1) were included in the simulation that generated Fig. 3 C. These additions resulted in a domain that looks qualitatively similar to the experimental domain (Fig. 3 A). The fact that the F/F decreased almost back to baseline values within 30 ms, and was constrained to a narrower region around the release site, is primarily due to the presence of Ca²⁺ buffers rather than the SR ATPase since removing the latter did very little to the F/F domain at this timescale (data not shown). Also, once the other buffers were included, the dependence of the Ca²⁺ domain on D₁ was almost eliminated. Unlike the case with no indicator present, changing D₁ from 1 x 10^-6 to 4 x 10^-6 cm²/s had a negligible effect on the peak magnitude of the Ca²⁺ domain (data not shown). This is because, in our model, most of the Ca²⁺ is rapidly bound and diffusion of the various Ca²⁺-dye complexes becomes the primary mechanism by which Ca²⁺ redistributes within the myofibril.

View larger version (39K): [in this window] [in a new window]   FIGURE 3  Comparison between model and experimental AP-evoked Ca2+ domains. (A) Experimental domain measured in a fiber equilibrated with 0.5 mM EGTA, 0.5 mM OGB-5N, 5 mM ATP, and 2 mM Mg. (B) Z-line model-simulated domain with only 0.5 mM OGB-5N and identical release flux parameters as Fig. 2. (C) Z-line model-simulated domain using the fiber conditions described in A and a release flux identical to Fig. 2. Simulation also contained parvalbumin and troponin as described in Table 2.

View larger version (11K): [in this window] [in a new window]   FIGURE 4  Comparison between average experimental Z- and M-line transients and Z-line model predictions. (A) The flux parameters (Jmax = 875 on = 5 ms, and off = 1.6 ms) were adjusted such that the simulated Z-line transient (dashed lines) matched the average Z-line transient (solid line) from 32 different sarcomeres (four fibers). The fiber was stimulated at 5 ms (arrows). F/Fpeak = 1.88, rt = 2.8 ms, FDHM = 5.94 ms, and All simulated and experimental transients were obtained under conditions identical to Fig. 3. (B) The M-line transient (dashed line) from the same simulation as A is compared to the averaged experimental M-line transient (solid line). Properties of experimental trace: F/Fpeak = 0.5, rt = 6.3 ms, FDHM = 16.2 ms, and Properties of simulation trace: F/Fpeak = 0.32, rt = 8.53 ms, FDHM = 31.37 ms, and

View larger version (13K): [in this window] [in a new window]   FIGURE 5  Effects of altering critical model parameters on the simulated M-line transient. The experimental AP was initiated at 5 ms (arrow). In all simulations, the flux was adjusted such that the modeled Z-line transient matched the experimental one. (A) The experimental M-line transient (thick solid line) is plotted against simulated M-line transients obtained from the Z-line release model in which D3 was 1.4 x 10-6 (dotted line), 3.8 x 10-6 (dashed line), or 9 x 10-6 (dash-dot line) cm2/s. Modeled transient properties are in the order 1.4 x 10-6, 3.8 x 10-6, and 9 x 10-6 cm2/s: rt = 12.1, 8.32, and 6.53 ms; F/Fpeak = 0.32, 0.46, and 0.70; ; and td = 1.81, 1.12, and 0.67 ms. (B) The experimental M-line transient (thick solid line) is plotted together with simulated M-line transients where the release site was extended to 0.6 (dotted line), 0.8 (dashed line), or 0.9 (dash-dot line) µm from the Z-line. The release flux was identical across the entire release band. Modeled release fluxes and transient properties are in the order 0.6, 0.8, and 0.9 µm. Jmax = 77, 87, and 80 on = 5, 5, and 5 ms; off = 1.3, 1.1, and 1.1 ms; rt = 9.76, 9.2, and 8.6 ms; F/Fpeak = 0.41, 0.49, and 0.51; ; and td = 1, 0.6, and 0.6 ms.

ATP is the most abundant Ca²⁺ buffer in our simulations and has been shown to have the potential to condition Ca²⁺ signaling within a muscle fiber (Baylor and Hollingworth, 1998). The ability of ATP to influence the simulations strongly depends on the steady-state free [ATP] (i.e., before the stimulus), which is assumed to be identical to that in the cut end pools (3 mM) (DiFranco et al., 2002). Even with such a high free [ATP], the D_ATP that was required in Z-line release model simulations to come close to predicting the M-line transients was 3.8 x 10^-6 cm²/s (data not shown), which is much higher than commonly assumed in the literature (Baylor and Hollingworth, 1998).

A critical parameter that can be used to compare experimental results with model simulations is the Z–M delay because it has the potential to distinguish between localized and extended Ca²⁺ release sites within the sarcomere (Hollingworth et al., 2000; Escobar et al., 1994). We have recently addressed this issue by carefully determining, under conditions that are expected to have high free [ATP], that the average Z–M delay in amphibian muscle fibers is 0.64 ± 0.1 ms (mean ± SD, n = 30) (DiFranco et al., 2002). In contrast, quantitative evaluation of the Z–M delay from Z-line model simulations, using the criteria described in Methods, predicts a significantly larger value of 1.5 ms. Only improbable values for diffusion coefficients were able to approximate the kinetic features of the M- and Z-line transients and the Z–M delay. Thus, we explored the possibility of increasing the extent of the release site while maintaining realistic diffusion coefficients. Fig. 5 B shows that M-line transients simulated using a D₃ value of 1.4 x 10^-6 cm²/s, but with an extended release site, can reproduce some of the M-line features (e.g., rise time, (d(F/F)/dt)_max, Z–M delay). It can be seen that, as the extent of the release site is increased from 0.6 µm (dotted line) to 0.9 µm (dash-dot line), the rising phase gets faster, the Z–M delay is significantly decreased, and the F/F_peak increases, yielding a much closer approximation to the rising phase of the experimental M-line transient. As expected, the Ca²⁺ flux required to predict the Z-line transient in the extended release site model decreased as the width of the release band was increased. Some discrepancies between predicted and experimental M-line transients still remain; particularly the falling phase of the experimental M-line transient which is significantly faster than predicted by the model. Additional simulations showed that the falling phase of the M-line transient is far more dependent on the density and kinetics of the SR ATPase than the specific properties of the release sites (data not shown).

So far, we have shown that the Z-line release model cannot predict the experimentally observed localized transients with reasonable values for the critical model parameters, but that the extended release model affords much better agreement. To more directly assess the kinetics and localization of the Ca²⁺ release from the SR, we chose to record and model AP-induced Ca²⁺ transients under conditions where the myofibril was loaded with 10 or 20 mM EGTA, as described in Methods. It has been suggested theoretically that these concentrations are sufficient to constrain the free [Ca²⁺] change to the immediate vicinity of the release site (Neher, 1998; Stern, 1992; Pape et al., 1995) and effectively accelerate the kinetics of the OGB-5N transients (DiGregorio et al., 1999). In Z-line model simulations the flux parameters were adjusted such that the simulated Z-line transient reproduced the experimentally observed one in 0.5 mM EGTA (Fig. 6 A). The Z-line transient from the same fiber equilibrated with 20 mM EGTA is shown in Fig. 6 B. It should be noticed that increasing the [EGTA] caused the F/F_peak to decrease from 1.7 to 0.65 and the full-duration at half-maximum (FDHM) to decrease from 7.4 to 4.1 ms. These effects were well-predicted by the Z-line model by simply changing the [EGTA] from 0.5 mM (Fig. 6 A) to 20 mM (Fig. 6 B). Very similar results were obtained in extended model simulations.

View larger version (21K): [in this window] [in a new window]   FIGURE 6  The effects of high [EGTA] on properties of the Ca2+signals. (A) The experimental Z-line transient (solid line) was recorded from a fiber equilibrated with 0.5 mM EGTA. The simulated transient (dotted line) was obtained using a Z-line release model. The simulated flux parameters were adjusted (Jmax = 850 on = 5.5 ms, and off = 5.8 ms) until the simulated Z-line transient matched the experimental one. (B) The experimental Z-line transient (solid line) was obtained from the same fiber as A, after equilibration with 20 mM EGTA. The simulated transient (dotted line) was obtained using a Z-line release model. (C) Averaged normalized F/Fpeak values from transients recorded at different distances from the Z-line in four different fibers, two with 10 mM EGTA and two with 20 mM EGTA (solid line and circles). The simulated profile (dotted line and open circles) was obtained from transients generated using the Z-line release model. The data was normalized to the F/Fpeak at the Z-line. (D) Normalized F/Fpeak values from transients recorded at different distances from the Z-line in extended model simulations compared to experimental data. The release site was extended to 0.6 (dotted line), 0.8 (dashed line), and 0.9 (dash-dot line) µm.

	DISCUSSION

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

 
In this article, we describe a three-dimensional diffusion-reaction model of a skeletal muscle myofibril that is based on the one described by Cannell and Allen (1984). Similar to Baylor and Hollingworth (1998), we have incorporated the diffusion and reaction of all relevant buffer species within the myofibril. However, we have extended their model in several ways to both improve the accuracy and convergence of the integration as well as to allow for the direct comparison of model predictions with experimental results obtained using the confocal spot detection system (Escobar et al., 1994; DiFranco et al., 2002; Vergara et al., 2001). Instead of the Euler explicit integration method used previously (Baylor and Hollingworth, 1998) we have utilized an implicit ADI method (Appendix B) that eliminates the necessity of using very small time-steps (Peaceman and Rachford, 1955; Smith, 1985), and has been successfully used to model systems of diffusion-reaction equations in other biological problems (Chiu and Walkington, 1997; Chiu and Hoppensteadt, 2000). The typical time-step required to obtain convergence in the Z-line release simulations was 4 µs, although a 2-µs time-step was required in simulations with large concentrations (>5 mM) of any free buffer. Due to their complexity, the extended release simulations required a shorter time-step of 0.2 µs to guarantee stability. These relatively large integration time-steps allowed us to divide the simulated myofibril into many more voxels than previously used (Cannell and Allen, 1984; Baylor and Hollingworth, 1998), thus increasing the accuracy of the model, although still allowing the simulation to complete in a reasonable amount of time.

Although a previous myofibril model (Baylor and Hollingworth, 1998) considered that Ca²⁺ and Mg²⁺ compete for binding to ATP, it incorporated a reaction scheme approximation that avoided an explicit calculation of the competition. However, comparison of our model predictions with and without this approximation showed that it induces an 10% error in the resulting F/F profiles; hence we chose to calculate the competition directly in every simulation.

In addition, several enhancements were required to compare model simulations with experimental data. Firstly, the simulated fluorescence distribution was blurred using a triple Gaussian approximation of the PSF, as determined directly for our experimental setup. Secondly, a localized transient at a single longitudinal position of the sarcomere was assumed to arise from a simulated detection pinhole that measures the blurred fluorescence from a square area centered in the middle z-axis plane of the myofibril. Thirdly, the simulated pinhole was moved longitudinally in 200-nm increments to construct a simulated fluorescence domain, which was compared with the experimental domains obtained from scanned data.

Since the simulated dye fluorescence was compared directly with the experimentally observed fluorescence, it was necessary to understand how the inclusion of dye affected the model predictions. A comparison of Fig. 2, A and B, illustrates how even the low affinity dye OGB-5N can distort the simulated [Ca²⁺] profile. Increasing the affinity of the modeled indicator to that of Fluo-3 (K_d = 0.7 µM, k_off = 0.175 ms^-1; Escobar et al., 1997), while maintaining the same F_max/min of OGB-5N, caused the F/F_peak to decrease by 69% and resulted in much shallower gradients, both radially and longitudinally (data not shown), than those illustrated in Fig. 2 D. This reinforces the idea that using a low affinity indicator affords the minimum distortion, both spatially and temporally, of the underlying Ca²⁺ gradients. Even with a low affinity indicator, the bound dye profile (Fig. 2 C) is a distorted representation of the free [Ca²⁺] in the myofibril, which highlights the fact that equilibrium assumptions are invalid for rapidly changing Ca²⁺ fluxes. However, closer inspection reveals that the majority of the distortion occurs within 100 nm of the release site. The shallower gradients that exist farther from the release site are relatively undistorted by OGB-5N (e.g., compare the longitudinal gradients at the ‘0 radial position in Fig. 2, A and C).

The fact that the [OGB-5N-Ca²⁺] profile in Fig. 2 C looks similar to the real [Ca²⁺] profile in Fig. 2 A is coincidental, due to the fact that there are no other buffers present. There is generally no simple relationship between the fluorescence profile (reporting the concentration of OGB-5N-Ca²⁺ complex) and the [Ca²⁺] that would have occurred in the absence of the dye. Thus, to fully understand the rapid evolution of the Ca²⁺ profile, it is necessary to combine high resolution measurements with computer modeling.

The blurring due to the PSF of the microscope objective introduces a further distortion of the Ca²⁺ signal. However, this effect is mainly noticeable in the radial, as opposed to the longitudinal, direction of the myofibril (Fig. 2 D). Thus, the longitudinal F/F gradients observed after blurring are relatively accurate representations of the true F/F gradients. However, we did introduce blurring in all of our simulations because it does cause distortion at early time points when the gradients are steep even in the longitudinal direction and blurred transients are more representative of the data obtained with the confocal spot detection system.

A convenient way to visualize the spatiotemporal profile of the Ca²⁺ signal is by constructing a three-dimensional plot representing the Ca²⁺ domain (Chad and Eckert, 1984; DiFranco et al., 2002; Neher, 1998). It is interesting that a domain that is qualitatively similar to an experimental one (see Fig. 3 A; see also Vergara et al., 2001; DiFranco et al., 2002) can be approximately predicted from a Z-line release model simulation that simply includes the appropriate Ca²⁺ buffers (Table 2, Fig. 3 C) and does not require the SR Ca²⁺-ATPase. A model of the Ca²⁺ signaling in presynaptic nerve terminals, which does not contain Ca²⁺ extrusion mechanisms, also predicts that the formation of domains around Ca²⁺ entry sites only requires the presence of Ca²⁺ buffers (DiGregorio et al., 1999). The minimal influence of the SR pump in our myofibrillar cylindrical model may be thought of as surprising, considering that we used higher values for V_max than previously suggested (Cannell and Allen, 1984), and results from its localization to the outermost of the 32 radial shells of a 0.5-µm radius myofibril. The outer shell's volume is 6.2% of the total myofibrillar volume, which is comparable to the 5% value calculated based upon electron microscopy stereological measurements (Mobley and Eisenberg, 1975). Models in which the SR pump is present in a larger proportion of the myofibrillar volume (Hollingworth et al., 2000; Baylor and Hollingworth, 1998) probably overestimate the volume from which the SR pump removes Ca²⁺ and therefore its contribution in shaping the Ca²⁺ transients.

Our model provides the first estimate of the Ca²⁺ release flux that is required to generate the F/F_peak values and kinetic properties observed in localized AP-induced experimental transients. As demonstrated in Fig. 4 A, a J_max of 875 pmoles/(cm² x ms) was required to predict the average experimental Z-line transient. Given the kinetics of the release function, the peak flux actually attained was 89 pmoles/(cm² x ms). This value corresponds to a total Ca²⁺ current of 16.4 pA per AP and results in a rate of [Ca²⁺] change of 55.5 µM/ms for a myofibril with a 0.5-µm radius and a 4-µm length. This value is 2.6-fold smaller than a previous estimate for amphibian muscle fibers based upon global Ca²⁺ measurements (Pape et al., 1995). These values are surprisingly similar, given the widely different experimental protocols and Ca²⁺ detection techniques methodologies. Based on comparisons of experimentally measured Ca²⁺ sparks (obtained with techniques more similar to those described in this article) with model predictions, Rios et al. (1999) calculated that sparks result from an 8 pA current. This would suggest that Ca²⁺ sparks represent the spontaneous activation of 50% of the total Ca²⁺ release flux that occurs in response to an AP during skeletal muscle EC coupling, seemingly too large for their proposed role as the elemental process underlying EC coupling in skeletal muscle (Klein et al., 1999; Rios et al., 1999).

Our model is sensitive to the free [ATP] since it is present in high concentrations and acts as a rapidly diffusing, low affinity Ca²⁺ buffer that carries Ca²⁺ throughout the myofibril (Baylor and Hollingworth, 1998). It should be noted that none of the parameters pertaining to ATP-Ca²⁺ interaction (i.e., ATP-k_on, ATP-k_off, D_ATP) have been measured directly inside the muscle fiber and the literature values that we chose are the ones that would generate the most rapid spread of Ca²⁺ throughout the myofibril. Moreover, since we assumed that the [ATP] and [Mg²⁺] inside the muscle fiber were initially equal to their concentrations in the cut end pools, the free [ATP] inside the myofibril was calculated to be 3 mM, which is much higher than typically assumed in other models (Hollingworth et al., 2000; Baylor and Hollingworth, 1998). Also, the value for D₃ of 1.4 x 10^-6 cm²/s used in the simulations is also fast compared to the value used for other indicators in the literature (Baylor and Hollingworth, 1998; Timmerman and Ashley, 1986; Harkins et al., 1993). Despite using the fastest reasonable diffusion parameter values, the Z-line release model could not simultaneously predict the kinetic features of the experimental M- and Z-line transients, including the Z–M delay (Fig. 4). We believe that the true values for the ATP and OGB-5N properties, as determined in vivo, will be slower than the values used in our simulations, and will only increase the discrepancies between the model predictions and the experimental data.

The Z–M delay provides a very sensitive measure of the diffusional distance from the closest Ca²⁺ release site to the M-line with minimum contamination from the kinetics of the Ca²⁺ release flux itself. Thus, the fact that the Z-line release model predicts a Z–M delay of 1.5 ms, as compared to the experimental one of 0.64 ms, is highly significant. Interestingly, simulations of extended release sites >0.6 µm per half-sarcomere have a Z–M delay comparable to the experimental data (Fig. 5 B). However, simply extending the release site does not allow us to predict the entire M-line transient, and discrepancies remain in the falling phase. We tapered the amplitude of the release flux as it approached the M-line in an attempt to simultaneously predict the rising and falling phases of the experimental M-line transient