DEFINITION OF SYMBOLS

TOP ABSTRACT DEFINITION OF SYMBOLS INTRODUCTION METHODS NUMERICAL METHODS AND COMPUTER... RESULTS DISCUSSION AND CONCLUSION APPENDIX A REFERENCES

Symbols in the text surrounded by brackets [ ] represent concentrations of the indicated species.

Ca2+	=	Calcium
CaATP	=	Calcium ATP
CaD	=	Calcium dye
CaP	=	Calcium pump
CaPARV	=	Calcium parvalbumin
CaTN	=	Calcium troponin
d	=	Depth of the computational domain
D	=	Dye
DATP	=	Diffusion coefficient of ATP
DC	=	Diffusion coefficient of Ca2+
DCaATP	=	Diffusion coefficient of CaATP
DCaD	=	Diffusion coefficient of CaD
DD	=	Diffusion coefficient of D
DMgATP	=	Diffusion coefficient of MgATP
F	=	Fluorescence signal
F0	=	Resting fluorescence signal
FDHM	=	Full duration at half maximum
FWHM	=	Full width at half maximum
h	=	Height of the computational domain
I	=	Magnitude of Ca2+ current
koff,CaATP	=	Backward rate constant for Ca2+ and ATP
koff,CaD	=	Backward rate constant for Ca2+ and D
koff,CaP	=	Backward rate constant for Ca2+ and pump
koff,CaPARV	=	Backward rate constant for Ca2+ and parvalbumin
koff,CaTN	=	Backward rate constant for Ca2+ and troponin
koff,MgATP	=	Backward rate constant for Mg2+ and ATP
koff,MgPARV	=	Backward rate constant for Mg2+ and parvalbumin
kon,CaATP	=	Forward rate constant for Ca2+ and ATP
kon,CaD	=	Forward rate constant for Ca2+ and D
kon,CaP	=	Forward rate constant for Ca2+ and pump
kon,CaPARV	=	Forward rate constant for Ca2+ and parvalbumin
kon,CaTN	=	Forward rate constant for Ca2+ and troponin
kon,MgATP	=	Forward rate constant for Mg2+ and ATP
kon,MgPARV	=	Forward rate constant for Mg2+ and parvalbumin
kP	=	Rate constant for SR Ca2+ pump
l	=	Length of the computational domain
L	=	Ca2+ leak from SR
LSCI	=	Laser scanning confocal imaging
LSCM	=	Laser scanning confocal microscope
Mg2+	=	Magnesium
MgATP	=	Magnesium ATP
MgPARV	=	Magnesium parvalbumin
MPE	=	Maximum percentage error
MSPE	=	Mean square percentage error
P	=	SR Ca2+ pump
PARV	=	Parvalbumin
PSF	=	Point spread function
RF	=	Reducing factor of diffusion coefficients at SR
Rr	=	Calcium release rate at the source
RMSPE	=	Square root of MSPE
RTime	=	Rise time, time to go from 10 to 90% of peak value
SR	=	Sarcoplasmic reticulum
Str	=	Starting values
TN	=	Troponin
V	=	SR Ca2+ pump transport
V0	=	Resting SR Ca2+ pump transport
t	=	Time step size
x	=	Spatial step size in x direction
y	=	Spatial step size in y direction
z	=	Spatial step size in z direction

	INTRODUCTION

TOP ABSTRACT DEFINITION OF SYMBOLS INTRODUCTION METHODS NUMERICAL METHODS AND COMPUTER... RESULTS DISCUSSION AND CONCLUSION APPENDIX A REFERENCES

Since the detection and report of Ca²⁺ "sparks" in cardiac (Cheng et al., 1993) and skeletal muscle (Tsugorka et al., 1995; Klein et al., 1996), considerable effort has been devoted to the investigation of the mechanism underlying the observed sparks. However, it is not obvious which aspects of Ca²⁺ sparks are determined by the activity pattern of the SR Ca²⁺ channel or channels which underly the generation of sparks, and which properties are determined by the diffusion and Ca²⁺ binding kinetics of both the intrinsic binding sites in the fiber and the Ca²⁺ indicator dye. Numerical simulation is a feasible approach to evaluate these factors in studying the Ca²⁺ movement underlying the sparks.

The complete simulation of Ca²⁺ sparks, as measured experimentally using laser scanning confocal microscopy in a skeletal muscle fiber, consists of two main parts. The first part is the simulation of the Ca²⁺-dye formation process and its spatiotemporal distribution by a diffusion and binding model. The second part is the simulation of the imaging process based on the results from the first part and the optical properties of the confocal microscope. This entire process is actually in close analogy to the experimental observation of the sparks, in which, first, the sparks are produced by the release and spread of Ca²⁺ from the Ca²⁺ release channel source in the junctional SR and the reactions of Ca²⁺ with fluo-3 and other myoplasmic Ca²⁺ binding sites. The change in Ca²⁺ dye is monitored as a change in fluorescence by the LSCM. In the present article, we first simulated the spread of Ca²⁺ dye due to the diffusion and binding processes, then simulated the image acquisition by the LSCM by convolving an approximation of the theoretical PSF of the LSCM with the simulated data from the diffusion and binding model.

Our model of diffusion and binding is a combination of structural domain configuration and functional equation formulation. The model was based on the physical structure of a sarcomere within a myofibril of a skeletal muscle fiber. The mathematical equation system was formulated on mass diffusion theory and chemical reaction kinetics. The following aspects represent the most important developments of the present three-dimensional (3D) simulation model from other published models of Ca²⁺ movements underlying Ca²⁺ sparks: 1) the diffusion of Ca²⁺ dye (D) and Ca dye (CaD) were included in the starting simulation; 2) several main Ca²⁺ binding sites were treated separately instead of being treated as a single lumped ligand (Blatter et al., 1997; Pratusevich and Balke, 1996), and the reaction rate constants were chosen from literature or experimental results; 3) the distributions of the binding sites were modeled in close analogy to the biological structure; 4) the computational domain was discretized into a graded network for the finite difference solution of the differential equation system, which ensured the necessary accuracy near the Ca²⁺ release source and reduced the computation load far from the source; 5) the numerical solutions of the diffusion processes of both an instantaneous source and a continuous source were verified by the comparison with the analytical solutions of the same processes under the condition of no binding reactions.

The main applications of the model have been the simulation of experimentally measured sparks; the comparison of sparks viewed at different positions relative to the Ca²⁺ source; and the investigation of the effects of changes in the Ca²⁺ release current source and of changes in the parameters for the Ca²⁺ binding sites, in the diffusion coefficients of D and CaD and in Ca²⁺ dye kinetics and concentration on the properties of the sparks. With appropriate sets of parameter values, the simulated Ca²⁺ sparks presented here reproduce many features of the amplitude and temporal properties of Ca²⁺ sparks recorded experimentally from frog skeletal muscle fibers (e.g., Lacampagne et al., 1996, 1998). However, the spatial spread of the simulated sparks, calculated as the spatial FWHM, was consistently less than experimentally measured values.

	METHODS

TOP ABSTRACT DEFINITION OF SYMBOLS INTRODUCTION METHODS NUMERICAL METHODS AND COMPUTER... RESULTS DISCUSSION AND CONCLUSION APPENDIX A REFERENCES

The diffusion and binding model

The central aspect of our simulation was the mass diffusion theory. The general equation for diffusion of substance A in a homogeneous and isotropic medium is (Crank, 1975)

(1)

where ² is the Laplacian operator. D_A is the diffusion coefficient of substance A and [A] is its concentration. This is a parabolic partial differential equation, the numerical solution of which is to track the time evolution of variable [A]. [A] is implicitly understood to be a function of time t and spatial coordinates x, y, and z, i.e., [A](t, x, y, z), but the time and space coordinates are not explicitly stated to simplify the notation. In Cartesian coordinates, ² = (²/x² + ²/y² + ²/z²).

When the diffusion process is involved with chemical reactions, the general Eq. 1 for diffusion, with the diffusing substance A given by free Ca²⁺ as in this study, is (cf., Cannel and Allen, 1984)

(2)

where

(3)

represents the effect of chemical reactions on the rate of change of [Ca²⁺]. B_i represent different binding sites and buffers. In the present study, we use the term "site" to indicate non- or slowly diffusible reactants binding Ca²⁺ and the term "buffer" to indicate diffusible reactants binding Ca²⁺. Eq. 2 is also the general equation used in the present study from which the structural domain configuration and the functional equation systems were developed.

Structural domain configuration

View larger version (38K): [in this window] [in a new window]   FIGURE 1   Schematic representation of the model structure. (A) The central cross section of the simulation domain, which cuts the entire simulation domain into two one-sarcomere sections. The fiber and myofibril axis is perpendicular to the plane shown. The whole simulation domain composed of 12 neighboring elementary sarcomere blocks with four-sarcomere deep (vertical), three-sarcomere wide (horizontal), and two-sarcomere long (not shown), which has the sizes of 4.0 × 3.0 × 8.0 µm. The cross section and the space of the computational domain are 1/4 and of those of the simulation domain, respectively. (B) The 3D view of the computational domain, which consists of a fine network block F and three coarse network blocks A, B, and C. Block F has 16,337 cubic voxels of the size of 0.03333 µm3 each. Blocks A, B, and C have 6480, 8784, and 31,232 cubic voxels, respectively, and the size of their voxels is 0.06663 µm3 each. The simulation results in coarse blocks A, B, and C will be interpolated and the total numbers of voxels in the entire computational domain are 62,833 and 373,527 before and after the interpolation, respectively. The x axis is along the long axis of the computational domain, and the y and z axes are the horizontal and vertical short axes of the computational domain. (C) and (D) The first Y-Z and X-Z planes of the computational domain.

Formulation of the functional equation system

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

(24)

(25)

View larger version (13K): [in this window] [in a new window]   Scheme 1

Simulation of laser scanning confocal imaging

The fluo-3 fluorescence signal recorded experimentally from a muscle fiber is considered to be proportional to [CaD] (Klein et al., 1996). However, because of the limitation in the spatial resolution of LSCM, the spark images obtained under the conditions of experimental observation contain distortion and blurring of the actual distribution of [CaD] in the muscle. Therefore, no individual line or plane of the domain of distribution of [CaD] produced by the diffusion and binding simulation can be directly viewed as equivalent to the spark image recorded by the LSCM. The optical distortion and blurring of the LSCM system must also be simulated to produce simulated Ca²⁺ sparks comparable to the experimentally detected sparks. The way in which a microscope distorts and blurs an individual point in the specimen is described by its PSF.

A Gaussian kernel was used to approximate the microscopic PSF in the present simulation having spatial FWHM dimensions of 0.5 µm, 0.5 µm, and 1.0 µm, in length (x), depth (y), and height (z), respectively, corresponding to the approximate resolution of our confocal microscope. These values were selected to be slightly larger than those determined from examination of subresolution fluorescent beads in air on our confocal system (0.4 µm, 0.4 µm, 0.8 µm; Schneider and Klein, 1996). The slightly larger half-width values were used in the present simulations in an attempt to account for the likely greater optical distortion within a fiber.

The image of the simulated spark was generated as done from the experimental recordings (Klein et al., 1997) by first subtracting the value before the Ca²⁺ release, [CaD]₀ or F₀, from the [CaD] or F values after the Ca²⁺ release starts, and then normalizing these [CaD] or F differences to [CaD]₀ or F₀ to give

(26)

	NUMERICAL METHODS AND COMPUTER PROGRAMMING

TOP ABSTRACT DEFINITION OF SYMBOLS INTRODUCTION METHODS NUMERICAL METHODS AND COMPUTER... RESULTS DISCUSSION AND CONCLUSION APPENDIX A REFERENCES

Numerical methods for the simulation of diffusion and binding

In the present simulation, the finite difference method was chosen as the numerical method for the solution of the partial differential equations. However, the explicit implementation of the finite difference method is only conditionally stable. In the present study, the stability criterion of the explicit method was satisfied by reducing the time step used in the simulation. The explicit Euler method (Atkinson, 1989) was used to solve the ordinary differential equations in the equation system, Eqs. 10-25. The complete set of finite difference equations for Eqs. 10-19 is listed in Appendix A.

The first step of the implementation of the finite difference method in the present simulation was the discretization of the computational domain. In this study, we adopted the definition that a voxel is the discrete volume element of the continuous specimen and the 3D equivalent of a pixel, which is the basic element of a 2D image (Pawley, 1995). Therefore the terms voxel, pixel, and grid point are often used interchangeably in the text. We discretized the computational domain into even-sided cubic voxels much smaller than the image pixels of the confocal microscope (Fig. 1, B-D). Using even-sided cubic voxels ensures the same accuracy of the finite differentiation in all three dimensions. A finer network near the Ca²⁺ source is necessary to ensure adequate computing accuracy of diffusion near the source, whereas a coarser network far away from the calcium source is more computationally efficient in significantly reducing the computing time of the simulation and required storage space of the results. Therefore, the entire computational domain was divided into a total of four network blocks, one fine network block (F) and three coarse network blocks (A, B, and C). A schematic diagram of the 3D view of the whole computational domain with the arrangement of the four blocks is shown in Fig. 1 B. The Ca²⁺ spark is considered to originate from a point source located at the origin of this coordinate system. The spatial distance between the neighboring points of the coarse network was set as two times that of the fine network, so that the points of the coarse networks at the interface coincide with half of the points of the fine network at the interface. Figure 1, C and D show the two side views of the computational domain. The simulation results within the three coarse blocks were linearly interpolated to match the spacing of the fine grid.

The main mathematical operation involved with the simulation of LSCI is the 3D discrete convolution of the spatial-temporal distribution of [CaD], as obtained from the solution of the diffusion and binding simulation, with the Gaussian kernel for the confocal microscope.

The code for the diffusion and binding simulation was written in FORTRAN 77 language under IRIX 6.2 on an SGI workstation (Mountain View, CA). Double precision variables were used in the computation to reduce rounding errors. All code was optimized manually as well as by setting the appropriate options for the FORTRAN compiler for efficient parallel processing. Diffusion and binding simulations were performed on a Power Challenge 10000 XL Series of SiliconGraphics Computer System, which has twenty 194-MHz IP25 processors.

In experimental recording of Ca²⁺ sparks, the line-scan fluorescence image is constructed by repeatedly scanning the same line along the fiber at 2-ms intervals (Schneider and Klein, 1996). The diffusion and binding simulation results were output at 1-ms intervals, for a total of 30 ms, which produces a higher temporal resolution in the simulated images of LSCM than that in the experimentally recorded images.

All the simulation outputs were saved as single precision data in binary file to reduce disk storage space. In addition, a text log file containing all the parameters used in the current simulation was produced at the beginning of each simulation for the reference of the analysis of the results.

Algorithm development for LSCI simulation

All the results from the diffusion and binding simulation need to be further processed for visualization and analysis. A program written in MATLAB5 was used to carry out the 3D discrete convolution for the simulation of the LSCI in this study. Because the main simulated scan line coincides with the long axis of the computational domain passing through the Ca²⁺ release source (the x axis of the coordinate system in Fig. 1 B), the consolidated domain was extended (flipped from the positive axis directions) to the x, y, and z directions to such an extent that there would be enough points around the x axis when the 3D discrete convolution of those points with the Gaussian kernel was performed along the x axis. Finally, the 3D discrete convolution was performed on this consolidated and extended domain at specified locations to produce the simulated sparks viewed from different positions and defocuses relative to the Ca²⁺ releasing point. Because of the symmetries in the simulation domain, any simulated spark produced in the computational domain is just half of the entire spark in the positive x direction. Each spark image was converted to a ratio image using Eq. 26 and "flipped" to the other half of the simulation domain in the negative x direction to form a complete spark.

Initialization of the simulation

We assumed that, at the far-end boundaries of the computational domain,

(27)

where h, d, and l are the height, depth, and length of the computational domain respectively (Fig. 1 B), and [Ca²⁺]₀ is the resting value of the [Ca²⁺]. Because the three near-end boundaries of the computational domain are the symmetrical planes of the simulation domain, no crossing flux exists. Therefore, the near-end boundary conditions are (Fig. 1 B)

(28)

The boundary conditions for D and CaD were treated by equations analogous to Eqs. 27 and 28 for Ca²⁺.

Before the start of the simulation, [Ca²⁺] and [D] were assumed to be uniformly distributed in the sarcomere, and all Ca²⁺ binding sites and buffers were assumed to be at their equilibrium levels corresponding to a resting [Ca²⁺]₀ of 80 nM (Klein et al., 1996). We use [CaTN]₀, [CaPARV]₀, [MgPARV]₀, [CaP]₀, and [CaD]₀ to denote the initial values of the concentration of the Ca²⁺ binding sites and buffers. Note that, even though spatial point indices have been ignored, the total concentration of troponin ([TN]_T) and SR Ca²⁺ pump sites ([P]_T) are different in different spatial locations, so the [CaTN]₀ and [CaP]₀ will also vary with locations. The expressions for the initial values of all the reactants are given below.

(29)

(30)

(31)

(32)

(33)

(34)

(35)

(36)

(37)

(38)

(39)

(40)

where


All diffusion and binding simulations started with a release of Ca²⁺ from a source at the origin. The rate, R_r, at which free [Ca²⁺] would increase in the origin voxel due to this release in the absence of diffusion and binding was given by

(41)

where I is the Ca²⁺ current,  is the Faraday constant, z is the valence of Ca²⁺ and v is the volume of the voxel of Ca²⁺ release. For our "starting" simulation, I is a square pulse with a magnitude of 1.4 pA (Pratusevich and Balke, 1996),  = 9.65 × 10⁴ Cmol¹, z = 2, and v = 0.0333³ = 3.6926 × 10⁵ µm³,

(42)

The start of this constant Ca²⁺ release was set at 1 ms after the start of the simulation, so that the first simulation output file at the end of the first millisecond contained the resting value under the simulation without release. This resting value was used in Eq. 26. The Ca²⁺ release time was chosen as a fixed 8.0 ms (Table 2) to match the mean rise-time (10-90%) of experimentally measured sparks (Klein et al., 1997). The total simulation time of 30.0 ms (Table 2) is long enough because sparks from experimental recordings have been identified with a 3.3-µm × 20-ms box (Klein et al., 1997).

Program testing

The FORTRAN programs for the diffusion and binding simulation were tested in several limiting cases. The computational and kinetic stabilities were tested under the condition of no Ca²⁺ release for the entire simulation period. After 1.35 × 10⁵ time increments at 2.2 × 10⁴ ms per increment, with repeated computations at each of over 60,000 grid points for each time increment, the MPE in [Ca²⁺], [D], or [CaD] among almost 400,000 grid points in the consolidated computational domain was less than 0.38%, demonstrating stability of the model in the absence of Ca²⁺ release. Similar results were obtained when simulation was performed with D_C = D_D = D_CaD = 0.0 and [Ca²⁺] = 100 µM at the release source, indicating that the diffusion coefficients were properly coded in the program. Another important test was to compare the numerical and analytical solutions for the case of diffusion from both an instantaneous and a continuous point source of Ca²⁺ under the condition of no binding reactions. The numerical simulation of pure diffusion was generated by setting all the kinetic rate constants and the initial [Ca²⁺], [D] and [CaD] to zero except the [Ca²⁺] at the origin and setting the reducing factor of the diffusion coefficients to 1.0 before the start of the simulations. The error in the numerical solution compared to the analytical solution was estimated by examining the MPE and the MSPE or the RMSPE. The MPE between the numerical and analytical solutions for diffusion along a specified line is the maximum of all the percentage errors at all points along the line. The percentage error at any point was calculated relative to the peak of the corresponding analytical solution along a line on which that point resides. The MSPE or RMSPE was used to describe quantitatively the general agreement of the numerical and analytical solutions. The MSPE for a comparable line was calculated by summing the squares of the percentage errors at all the points on the line and dividing by the total number of the points on the line. For the case of pure diffusion from an instantaneous point source, the starting [Ca²⁺] at the origin was set to 100 mM and then allowed to diffuse throughout the surrounding medium. The numerical and analytical solutions were compared at 0.133 ms after the diffusion started, before the diffusion reached the nearest boundary of the computational domain. The MPE for the comparisons along x, y, and z axes of the computational domain were 0.26, 1.14, and 0.26%, respectively. The MSPE (RMSPE) for these comparisons were 0.02 (0.13), 0.24 (0.49), and 0.01 (0.12), respectively. Pure diffusion from a continuous point source of constant Ca²⁺ release rate, corresponding to a 1.4 pA Ca²⁺ current as in the spark simulation, was used to generate both numerical and analytical solutions. Because the analytical solution used in this comparison represented a continuous point source at a constant release rate, we set the Ca²⁺ efflux time to the entire calculation period instead of only a fraction of the entire calculation period as in spark simulations, to obtain a pseudosteady-state distribution of Ca²⁺ for comparison with the analytical solution for diffusion from a point source (Crank, 1975). Two comparisons were made on the numerical and analytical solutions of the diffusion from a continuous source. One was the comparison of the numerical and analytical time courses of the [Ca²⁺] at a point about 0.2 µm from the source, and the other a comparison of the numerical and analytical spatial spreads of [Ca²⁺] along three coordinate axes at the end of 0.5 ms. For the time course of about 3.0-ms duration the MPE was 1.42% and the MSPE (RMSPE) was 0.52 (0.72). For the spatial spread comparisons in x, y, and z directions, the MPE and MSPE (RMSPE) were 0.83%, 0.03 (0.17); 0.75%, 0.14 (0.37); and 0.84%, 0.03 (0.16), respectively. These comparison results indicated that the implementation of the numerical method in the FORTRAN programs was satisfactory.

	RESULTS

TOP ABSTRACT DEFINITION OF SYMBOLS INTRODUCTION METHODS NUMERICAL METHODS AND COMPUTER... RESULTS DISCUSSION AND CONCLUSION APPENDIX A REFERENCES

Selection of starting values of model parameters

The simulation of Ca²⁺ sparks began with the parameter values listed in Tables 1 and 2, which serve as the "starting" parameter values for our simulation. These values were selected to agree with parameter values used in various earlier simulations of average or local changes in free and bound Ca²⁺ during Ca²⁺ release in muscle fibers. Our choice of these starting values represents a somewhat arbitrary starting point in parameter space, from which we can explore the effects of changing various values on the simulated sparks.

Adjustments were made to the values of [TN]_T and [P]_T given in Table 2 based on the nonuniform distribution of troponin and the SR Ca²⁺ pump within the sarcomere. [TN]_T = 240 µM is the average concentration of troponin C Ca²⁺-specific binding sites over the entire fiber volume. Since troponin C is located only in the I band, which constitutes only half of the sarcomere when stretched to 4 µm per sarcomere, we doubled the troponin concentration to [TN]_T,I = 480 µM for the total troponin concentration in the I band region of the sarcomere in our simulation within a 1-µm longitudinal distance from the Z line and set it to zero in the remainder of the sarcomere. The Ca²⁺ pump was considered as being distributed only on the surface of the myofibril where the SR is located instead of distributed uniformly throughout the whole volume of the sarcomere. Therefore, an adjustment was made to compensate for the difference between the [P]_T on the surface of sarcomere and that averaged over the entire volume of the sarcomere. In our fine-mesh network, the number of voxels on the surface of the myofibril is the number of voxels in the volume of the myofibril. The [P]_T value on the surface, [P]_T,S, was thus assumed to be eight times of that of the average concentration over the entire volume (200 µM), and therefore, a value of [P]_T,S = 1600 µM was used in the surface voxels in the fine mesh region in the starting simulation.

Spark simulation using starting parameter values

The characteristics of Ca²⁺ sparks are usually described by the peak amplitude of the spark, the rise time (here defined as the time to go from 10 to 90% of peak amplitude), the temporal FDHM and the spatial FWHM. We therefore examined the simulated Ca²⁺ sparks by plotting the values along two orthogonal lines passing through the location of the peak of the spark, giving the time course and the spatial spread of the simulated spark as shown in Fig. 2, A and B, respectively. This spark was simulated with the parameters from Tables 1 and 2 for a 1.4-pA point source of Ca²⁺ release located in the terminal SR and active for 8 ms. The spark time course at the center of the simulated spark (Fig. 2 A) rises monotonically throughout the 8-ms period of simulated Ca²⁺ release, and then abruptly begins to decline at the cessation of Ca²⁺ release. The spatial profile (Fig. 2 B) at the time of the peak of the simulated spark exhibits a bell-shaped distribution. The Peak amplitude of the simulated spark is 1.12 in the dimensionless units of [CaD]/[CaD]₀, which correspond to the dimensionless units of F/F₀ of experimentally recorded Ca²⁺ sparks, and its rise time, FDHM and FWHM were 6.1 ms, 7.7 ms, and 1.0 µm, respectively. Because the simulated spark (Fig. 2, A and B) was generated from a conservative set of parameters (Tables 1 and 2), it was considered as our starting spark for comparison with sparks simulated using other alternative sets of parameters (below).

Comparison with experimentally recorded sparks

An important issue concerning the Ca²⁺ sparks simulated with any model is the extent to which the simulated sparks agree or disagree with experimentally observed sparks. Ultimately, the goal of modelling a Ca²⁺ spark is to reproduce as closely as possible any experimentally recorded Ca²⁺ spark. This long-range goal is, however, beyond the scope of the present paper, which primarily introduces our present model and examines how the values of various parameters in the model influence the simulated sparks. None the less, it is still important that the results of our present simulations exhibit characteristics that are generally consistent with the average properties of experimentally observed Ca²⁺ sparks.

Figure 2, C and D present the time course and spatial distribution at time of peak obtained by spatiotemporally superimposing and averaging six individual experimentally recorded Ca²⁺ sparks in a train of repetitive mode sparks (Klein et al., 1999). Repetitive mode events were selected for this display because they arise at a single triad, apparently from the repetitive opening of either the same group of Ca²⁺ release channels or of the same single channel (Klein et al., 1999). These repetitive events are thus not influenced by the effects of variations in the distance of the scan line from the release channels (below). The spark simulated with the starting parameter values (Fig. 2 A) has only about half the amplitude of the experimental spark in Fig. 2 C, indicating that the release current magnitude may be too low to account for this experimental record (see below), but the rise time (6.1 ms) and FDHM (7.7 ms) of the simulated spark are roughly similar to those of the experimental record (4.7 and 8.4 ms, respectively). The spatial distribution (Fig. 2 B) of the spark simulated with the starting parameter values is considerably more narrow (FWHM 1.0 µm) than that of the experimental spark (Fig. 2 D; FWHM 2.2 µm). A more extensive comparison of the properties of simulated and experimentally measured sparks will be presented in the Discussion.

View larger version (23K): [in this window] [in a new window]   FIGURE 2   Ca2+ sparks simulated using the starting set of parameter values in the sarcomeric diffusion and binding model. (A) Time course of the peak of the simulated Ca2+ spark. In most cases this time course was used to characterize Ca2+ sparks simulated with different parameters. (B) Spatial spread passing through the peak of the simulated starting spark. This simulated spark was considered to be the starting spark, and the parameters (Tables 1 and 2) used in its simulation were defined as our starting parameters. (C) and (D) Time course and spatial spread of experimentally observed spark. These records are the average of six spontaneously activated repetitive events arising at the same location (Klein et al., 1999). (E) and (F) Comparisons of normalized, with corresponding peaks of each line, time courses, and spatial spreads of the simulated [CaD]/[CaD]0 before () and after () its convolution with the Gaussian kernel. The peak values of [CaD]/[CaD]0 were 4.88 before convolution and 1.12 after convolution.

Effects of the confocal imaging system

The time courses and lateral distributions of both the theoretical (Fig. 2, A and B) and experimentally observed (Fig. 2, C and D) confocal fluorescence signals were each calculated or recorded with finite optical and spatial resolution through the confocal microscope system. The theoretical records were obtained by first calculating the actual value of [CaD]/[CaD]₀ for each voxel in the solution space and then simulating the confocal line-scan image of this calculated spatial distribution (see Methods). In this case, it is thus possible to directly compare the theoretical confocal records with the underlying calculated change occurring within any individual voxel. Figure 2 E presents a comparison of the normalized theoretical time courses calculated within the origin voxel into which SR Ca²⁺ release occurs, with the theoretical time course for a confocal recording centered at the same voxel. In contrast to the rather smoothly declining rate of rise and fall of the simulated confocal time course (points), the local time course within the origin voxel (line) exhibits a more biphasic rise and fall, attaining more than 70% of its final level within the first millisecond and then gradually approaching its final level. This difference arises from the finite spatial resolution of the confocal system, which samples a range of voxels within the Gaussian kernel of the origin voxel, representing the spatial range of sampling in confocal recording. Such extended spatial sampling tends to smooth out the abrupt rise and fall of fluorescence calculated for the origin voxel itself. As a result of the much broader sampling of the confocal simulation compared to the individual origin voxel, the actual peak changes in the theoretical [CaD]/[CaD]<