INTRODUCTION

Top Abstract Introduction Methods Results Discussion References

Calcium (Ca²⁺) "sparks" are brief, spatially localized Ca²⁺ release events resulting from the opening of one or a cluster of sarcoplasmic reticulum (SR) Ca²⁺ release channels. The combination of laser scanning confocal microscopy and the fluorescent Ca²⁺ indicator, fluo-3, has revealed Ca²⁺ sparks in single cardiac cells (Cheng et al., 1993; López-López et al., 1995), cardiac trabeculae (Wier et al., 1997), skeletal muscle cells (Tsugorka et al., 1995; Klein et al., 1996), and smooth muscle cells (Nelson et al., 1995).

The observation of spontaneous Ca²⁺ sparks (Cheng et al., 1993) and the description of the voltage dependence of evoked Ca²⁺ sparks (López-López et al., 1995; Cannell et al., 1995) provided important experimental support for the local control theory of cardiac excitation-contraction coupling (Stern, 1992). It is now generally accepted that Ca²⁺ current through L-type Ca²⁺ channels locally triggers SR Ca²⁺ release. Evidence for this comes indirectly from measurements of Ca²⁺ sparks under whole cell voltage clamp (López-López et al., 1994, 1995; Santana et al., 1996) and more directly from measurements of Ca²⁺ sparks localized to the region under a cell-attached membrane patch (Shorofsky et al., 1998).

While the trigger for SR Ca²⁺ release is understood, how the release is controlled is still unclear. A myriad of factors appear to control SR Ca²⁺ release. For example, 1) SR load affects both the probability of Ca²⁺ spark occurrence and the Ca²⁺ spark amplitudes (Satoh et al., 1997; Györke et al., 1997; Song et al., 1997); 2) the number of SR Ca²⁺ release channels opening to generate a Ca²⁺ spark might be variable (Lipp and Niggli, 1996) and produce Ca²⁺ sparks of different amplitudes (Parker and Wier, 1996); 3) the SR Ca²⁺ release channel appears to have a subconductance that results in two populations of Ca²⁺ sparks with different amplitudes (Cheng et al., 1993; Xiao et al., 1997); 4) Ca²⁺ sparks can trigger other Ca²⁺ sparks (Klein et al., 1996; Parker et al., 1996; Blatter et al., 1997), that might appear as population of Ca²⁺ sparks with different amplitudes (Klein et al., 1996); and 5) the proportion of Ca²⁺ sparks in different amplitude populations appears to be altered in certain disease states (Shorofsky et al., 1996, 1997).

To better understand the mechanisms underlying the control of Ca²⁺ release we need to estimate the current through the SR Ca²⁺ release channel, but this current cannot be measured directly in an intact cell and must be estimated from the amplitude of the Ca²⁺ spark. This estimate is complicated by 1) the kinetics and capacity of Ca²⁺ buffering by endogenous Ca²⁺ buffers and exogenous Ca²⁺ buffers, such as the Ca²⁺ indicator fluo-3 (Balke et al., 1994; Smith et al., 1996); 2) the diffusion properties of free Ca²⁺ and mobile Ca²⁺ buffers (Wagner and Keizer, 1994; Smith et al., 1996); and 3) the uncertainty of the distance between the site of Ca²⁺ release and the position of the confocal linescan (Pratusevich and Balke, 1996; Shirokova and Ríos, 1997).

Pratusevich and Balke (1996) showed that the random placement of the confocal linescan relative to the Ca²⁺ release sites results in a broad distribution of measured Ca²⁺ spark amplitudes even if all Ca²⁺ sparks were generated identically. Thus it becomes difficult to distinguish variations in Ca²⁺ spark amplitude due to intrinsic variations in SR Ca²⁺ release channel current from positional variations arising from varying distance between the confocal linescan and the site of Ca²⁺ release. In this paper we show, in theory, that variations in Ca²⁺ spark amplitude arising from intrinsic and positional changes can be separated. We do this by deriving an integral equation that gives the relationship between the probability distribution of source strengths (SR Ca²⁺ release channel current or channel open time) and the Ca²⁺ spark amplitude histogram. This integral equation takes into account the effect of positional variations on the Ca²⁺ spark amplitude histogram. The integral equation can be solved analytically, allowing us to explicitly solve for the probability distribution of source strengths when given the measured spark amplitude histogram. In this way the variations in Ca²⁺ spark amplitude due to intrinsic and positional changes, which are intertwined in the amplitude histogram, can be separated.

	METHODS

Top Abstract Introduction Methods Results Discussion References

To understand the relationship between the measured Ca²⁺ spark properties and the underlying events we need to simulate the processes that influence the formation of the Ca²⁺ spark. These processes are (1) the release of Ca²⁺ from the SR by the opening of a Ca²⁺ release channel; (2) the diffusion of Ca²⁺ into the cytoplasm; (3) the reaction of Ca²⁺ with endogenous buffers, such as troponin-C, and the fluorescent Ca²⁺ indicator fluo-3; (4) the formation of the optical image of the Ca²⁺-bound fluo-3 dye; (5) the generation of a linescan image from the optical signal; (6) the generation of random fluctuations of the fluorescent signal due to photon and other sources of noise; and (7) the detection of the Ca²⁺ spark.

Reaction-diffusion equations

Processes (1), (2), and (3) are captured in the set of partial differential equations describing the reaction of Ca²⁺ with buffers and the diffusion of Ca²⁺ in the cytoplasm. The solution of these equations gives the 3-dimensional distribution of the Ca²⁺-bound fluo-3 as a function of time.

The chemical species included in the model equations are Ca²⁺ (concentration denoted by C), immobile endogenous buffers both free (F_b) and bound to Ca²⁺ (G_b), mobile fluo-3 (free F_m, bound G_m), and immobilized fluo-3 (free F_i, bound G_i). Mass transport of Ca²⁺ and mobile fluo-3 is assumed to follow Fick's law and the reaction rates (R_j) are governed by mass action kinetics. Thus the reaction rate between Ca²⁺ and any buffer is given by

(1)

where j = i, m, or b and k_j⁺ and k_j are the forward and reverse rate constants, respectively.

We also make the following assumptions: (1) Ca²⁺ released from the SR is approximated by a point source; (2) both reaction and diffusion occur radially symmetrically; (3) the diffusion coefficients of the Ca²⁺-bound mobile fluo-3 and the free mobile fluo-3 are identical; and (4) before the opening of the Ca²⁺ channel, all chemical species are at their steady-state value and there are no spatial gradients.

Under these assumptions the reaction-diffusion equations are

(2)

(3)

(4)

(5)

where

(6)

and

(7)

The Laplacian operator ² for the radially symmetric domains is

(8)

The point source of Ca²⁺ release by the SR is located at the origin and is given by J_SR(r), where (r) is the Dirac delta-function. J_SR is related to the Ca²⁺ release channel current I_SR by J_SR = I_SR/(z) where z = 2 is the Ca²⁺ valence and  is Faraday's constant. A typical value of I_SR is 1.4 pA.

The reason there are no equations for G_j is because under assumption (3) the sum H_j satisfies the linear diffusion equation H_j/t = D_j²H_j, and under assumption (4) the initial condition satisfies H_j(r, 0) = _j = constant, so the diffusion equation has the solution H_j(r, t)  _j. Accordingly, F_j and G_j satisfy the algebraic relationship G_j = H_j  F_j = _j  F_j.

For simplicity, we did not include a Ca²⁺ pump in the model because others (Gómez et al., 1996) found that 80% of the decline of the Ca²⁺ fluorescence signal could be accounted for by diffusion and buffering.

The diffusion coefficient for Ca²⁺ was set to 6 × 10⁶ cm²/s. The apparent diffusion coefficient of fluo-3 is found to be 0.2 × 10⁶ cm²/s in frog skeletal muscle, which is about a factor of 5 times smaller than predicted from its molecular weight (Harkins et al., 1993). To account for this difference Harkins et al. (1993) estimated that 78% of the dye is bound to immobile myoplasmic constituents and only 22% is freely mobile. In our simulations the ratio of concentrations () of immobile to mobile fluo-3,  was set to 5 or in some cases 2. This latter value comes from prior determination of the ratio of immobile to mobile fura-2 in guinea pig heart cells (Blatter and Wier, 1990). D_m was set to 0.9 × 10⁶ cm²/s. The free mobile fluo-3 concentration was fixed to 50 µM in all simulations.

The rate-limiting step of most reactions involving Ca²⁺ is the dehydration of the calcium ion and is ~200-700/µM s (Hague, 1977). We chose a value for the forward rate constant k_m⁺ to be near the middle of the range, 400/µM s. The reverse rate constant, k_m = 160/s, was calculated using the dissociation constant value of 400 nM.

The total concentration of endogenous buffers _b was set to 123 µM (Berlin et al., 1994). The forward rate constant k_b⁺ was chosen to be 100/µM s and the reverse rate constant of k_b = 100/s to give the endogenous buffer dissociation constant of 1 µM, close to the value (0.96 µM) found by Berlin et al. (1994).

The experimental parameters were determined at room temperature (20-25°C) except for Harkins et al.'s estimate of  = 5, which was measured at 16°C. No temperature compensation was made.

Numerical solution

Eqs. 2-5 were solved using Facsimile (AEA Technologies, Harwell, UK), which solves each of the equations in the time variable on a workstation (RS 6000, IBM Corp., Armonk, NY). The spatial domain extending from 0 < r < L was discretized into N compartments of equal lengths h = L/N. The ith compartment is the spherical shell bounded by r_i = ih and r_i+1 = (i + 1)h. The time rate of change of the concentration in the ith compartment is

(9)

where A_i = 4i²h² is the surface area of the sphere of radius r_i and V_i = 4h³(i² + i + 1/3) is the volume between the spheres of radius r_i and r_i+1. R_i is the reaction term. J_i is given by Fick's law J_i = (D/h)(c_i  c_i1). The differential equation for the Nth compartment is derived similarly, but imposing the zero-flux boundary condition that requires the fictitious point c_N+1 = c_N. For i = 0, which contains the point source, material balance yields

(10)

where V₀ = 4h³/3.

We used L = 6 µm and N = 600. This code was tested on the linear problem obtained by setting all reaction terms to zero for which the analytic solution was available for comparison. Except at very early times following channel opening (~10 µs) the relative error was within 5% even at the smallest resolvable distance of r = 0.015 µm. In all simulations the concentration at r = 6 µm did not vary over the short (<200 ms) integration time so any of the usual boundary conditions (Dirichlet, Neumann, or Robin) would give essentially identical results.

Formation of the optical image

The Ca²⁺ spark is the optical image of the distribution of Ca²⁺-bound fluo-3, G_m(r, t) + G_i(r, t). Any optical instrument forms a blurred image of the object and the extent of the blurring is given by the instrument's point spread function (PSF). We used a 3-dimensional Gaussian function as the PSF of the confocal microscope

(11)

where N = (_xy²_z^3/2)¹ normalizes the integral of the PSF over all space to unity. The standard deviation  is related to the confocal full-width at half-maximum (FWHM) by  = FWHM/[2(log 2)^1/2]. A typical value for the lateral FWHM is FWHM_xy = 0.4 µm. Values for the axial FWHM range from 0.41 µm (Parker et al., 1997) to 1.3 (Pratusevich and Balke, 1996).

The intensity contributions to the image of the Ca²⁺-bound fluo-3 at any point (x, y, z) measured from the point source at the origin is proportional to the convolution

(12)

The convolution was carried out by multiplication of the discrete approximation of the Fourier transforms of G_m + G_i and PSF then taking the inverse transform. The size of the volume to carry out the convolution was dictated by the extent of spreading of the bound indicator and the values of the axial and lateral FWHM. A typical computational volume (of 64³ elements) measures 3.6 µm along the x- and y-directions and 4.5 µm in the z-direction.

Generating a linescan image

Generation of the linescan image of a Ca²⁺ spark starts by choosing the linescan position (y*, z*) in the y  z plane, which is perpendicular to the linescan direction along x. Then for each time point t_j values of I(x, y*, z*, t_j) for all x are collected. Stacking these one-dimensional arrays for all the computed times (0 < t_j < 180 ms) produces the linescan image of the Ca²⁺ spark. The length of the linescan image along x is 4 µm and 180 ms in time. This small linescan image is embedded in a larger array whose values are set to the image value of bound fluo-3 (G_m + G_i) at equilibrium. Additionally, multiple Ca²⁺ sparks can be embedded in the large array at random positions, with the constraint that Ca²⁺ sparks do not overlap. The result of this embedding is an image that looks qualitatively like a linescan image from a real confocal microscope. A sample image in which signal fluctuations have been added is shown in Fig. 1 B.

View larger version (115K): [in this window] [in a new window]   FIGURE 1   Comparison of an actual confocal linescan image (A) and a simulated linescan (B). Periodically spaced horizontal lines in (A) are located at the t-tubules and may represent inhomogeneous dye distribution. No dye inhomogeneities exist in the simulations, so the background is uniform in (B). The simulated Ca2+ sparks are qualitatively similar to real Ca2+ sparks. The number of Ca2+ sparks per linescan image was random. The space and time positions of the Ca2+ sparks were also random, but were constrained not to overlap. Ca2+ sparks labeled a and b arise from linescan positions marked in Fig. 3 A. Each pixel is 0.1 µm by 3 ms.

We created realistic linescan images because we wanted observers to identify these simulated Ca²⁺ sparks in order to study the role of subjective factors on Ca²⁺ spark identification.

Random fluctuations of the fluorescent signal

Random fluctuations of the fluorescent signals arise from the intrinsic granularity of photons and from electronic noise (Pawley, 1995). To accurately model noise in simulated confocal linescan images, we measured the noise properties in linescan images made from the homemade confocal microscope (Parker et al., 1997). The mean signal level and the standard deviation were calculated in a 10 × 10-pixel sample area in three regions: background, where the fluorescence intensity was low and uniform; regions containing a narrow band of elevated fluorescence at the site of the t-tubules (Shacklock et al., 1995; Klein et al., 1996); and at Ca²⁺ sparks. Fig. 2 shows that the standard deviation of the values in the sample areas increase linearly (slope = 0.3) with the mean fluorescence, and this linear relationship holds regardless of the sample region. Moreover, the distribution of noise values is approximately Gaussian (data not shown). We therefore added to the value of each point in our simulated linescan image a random number from a Gaussian distribution whose standard deviation was 0.3 times the value at that point.

View larger version (17K): [in this window] [in a new window]   FIGURE 2   Statistical fluctuation measured by the standard deviation of fluorescence signal as a function of mean fluorescence level. Data were collected in a 10 × 10 pixel area in three regions: background (circles), t-tubular region (squares), and Ca2+ sparks (triangles). Measurements were made on linescan images collected using the homemade confocal microscope. The best fit line has a slope of 0.3.

Automatic detection of Ca²⁺ sparks

We developed a program to automatically identify Ca²⁺ sparks in linescan images. This program relieves the tedium of manually identifying Ca²⁺ sparks and ensures a more consistent choice of Ca²⁺ sparks than might be achieved by observers. This program identifies as Ca²⁺ sparks regions that have a sufficiently high density of pixels that exceed some threshold level. Identification of Ca²⁺ sparks starts by creating a binary image in which all pixels in an image whose value is less than the background level + threshold are set to zero, and all other pixels to unity. The threshold equals the standard deviation of the background signal times a factor (typically 1.4) that can be varied by the user.

High density regions of non-zero pixels are identified by using the following procedure iteratively. At every pixel in the image (i, j), the number of non-zero pixels within a square neighborhood of size N_size centered on (i, j) are counted. If this number is less than N_live, then the (i, j) pixel is set to 0 (i.e., the pixel "dies"); otherwise it is set to 1 (i.e., the pixel "survives" or is "born"). This procedure is repeated N_generation times for every pixel. As the notation suggests, this algorithm is based on ideas gleaned from modeling density-dependent population growth using cellular automata. Although the procedure appears slow and tedious, it in fact runs quickly with the array-oriented programming language IDL (Research Systems, Inc., Boulder, CO). The number of "live" neighbors a pixel has is found by doing a boxcar averaging of size N_size × N_size (typically N_size = 7) on the binary image. This smoothed array is thresholded-setting all pixels whose value is less than N_live to 0 and 1 otherwise.

Before processing actual linescan images, the prominent horizontal lines seen in many images (see Fig. 1 A) are removed to avoid being identified as potential Ca²⁺ sparks by the detection program. This is done by setting the zero frequency component (corresponding to time) of the image's Fourier transform to zero. The linescan image without horizontal lines is recovered by inverting the modified transform.

	RESULTS

Top Abstract Introduction Methods Results Discussion References

Confocal images of Ca²⁺ sparks

Fig. 1 A shows a linescan image of a rat ventricular cell obtained by using our homemade confocal microscope system (Parker et al., 1997). This system has a lateral FWHM of 0.31 µm and an axial FWHM of 0.41 µm. The bottom panel shows a simulated linescan in which both axial and lateral FWHM values were set to 0.35 µm. Setting the axial FWHM equal to the lateral FWHM greatly simplifies the analytic calculations without sacrificing much accuracy.

Before going to a quantitative description of the simulated Ca²⁺ sparks we point out two features of the simulated linescan. First, 1 A shows prominent horizontal lines that are spaced ~2 µm apart vertically; they originate on the t-tubules (Shacklock et al., 1995) and may arise from inhomogeneous distribution of dye. These lines are absent in the simulated image as we have assumed that the dyes, both mobile and immobile, are initially homogeneously distributed. Second, apart from the absence of the streaks, the simulated linescan image looks qualitatively like the linescan from an actual experiment.

Fig. 3 A illustrates how the simulated Ca²⁺ sparks are generated. The point source of Ca²⁺ release is assumed to be at the origin in panels A-C. The circles in the figure show the randomly chosen positions of the linescan in the y-z plane. Each circle corresponds to at least one Ca²⁺ spark (some positions may be chosen more than once) in the 200 linescan images made for this particular simulation. The triangles indicate the position of those Ca²⁺ sparks that were detected by the Ca²⁺ spark detection program. The arc encloses the region where 90% of the Ca²⁺ sparks were detected; the arc radius R₉₀ is 0.56 µm. The linescan positions of Ca²⁺ sparks labeled a and b in Fig. 1 are correspondingly labeled in Fig. 3 A. As expected, the bright Ca²⁺ spark (a) arose when the linescan was near the point source and the dim Ca²⁺ spark (b) arose when the linescan was far away.

View larger version (49K): [in this window] [in a new window]   FIGURE 3   Spatial distribution of detected Ca2+ sparks as the microscope FWHM changes. In (A)-(C) the small circles mark the locations of the linescan in the y-z plane with respect to the source located at the origin. The triangles show the linescan positions at which the Ca2+ spark could be detected. The arc encloses the region where 90% of the detected Ca2+ sparks were found. In (A) the lateral and axial FWHM were FWHMx,y = FWHMz = 0.35 µm. The detected Ca2+ sparks are symmetrically disposed about the origin as expected. The points marked a and b are the linescan positions that generated the similarly labeled Ca2+ sparks in Fig. 1 B. In (B), FWHMx,y = 0.2 µm and FWHMz = 0.6 µm. Although FWHMz is three times as large as FWHMx,y the linescan positions of detected Ca2+ sparks are still symmetrically arrayed about the origin. In (C), FWHMx,y = 0.2 µm and FWHMz is six times as large. The linescan positions of detected Ca2+ sparks are no longer symmetric about the origin, but show an elliptical pattern. In the absence of noise, all Ca2+ sparks are detected by the spark detection program. (D) shows the contrast increase of a Ca2+ spark as the axial FWHM decreases from 1.2 to 0.6 to 0.2 µm while keeping the lateral FWHM fixed to 0.2 µm.

Performance of the Ca²⁺ spark detection program

As shown in Fig. 3 A a large number of Ca²⁺ sparks in the linescan images go undetected. The sensitivity of the program to pick out dim Ca²⁺ sparks can be altered by changing the detection parameters N_size, N_live, and N_generation. Decreasing N_size or N_live increases the sensitivity of the program allowing detection of dimmer Ca²⁺ sparks, but at the expense of making more false identifications. Increasing N_generation has only a small effect on the sensitivity but reduces the number of false identifications. We could check the false identification rate because the positions of all Ca²⁺ sparks in the linescans were known. Note that falsely identified Ca²⁺ sparks were excluded from our measurements. The program parameters were adjusted empirically to achieve a balance between sensitivity and low false identification rate. We found that by using N_size = 7, N_live = 12, and N_generation = 3 the program identified all Ca²⁺ sparks correctly identified by observers and correctly identified dim Ca²⁺ sparks not identified by observers, while maintaining a false identification rate of ~2-5%. The number of dim Ca²⁺ sparks found by the program that was not detected by observers varied between observers but the program typically found ~50% more of the dimmest detectable Ca²⁺ sparks. The processing time for 100 linescan images, 166 pixels (500 ms) by 256 pixels (25.6 µm) in size, is ~2 min on an IBM RS-6000 workstation.

Properties of simulated Ca²⁺ sparks

Ca²⁺ sparks shown in Fig. 1 B were generated using a channel current of 1.4 pA, a channel open time of 10 ms (Rousseau and Meissner, 1989; Lukyanenko et al., 1996) with  = 5. We also ran an identical simulation except with  = 2. Ca²⁺ spark characteristics from both simulations are shown in Table 1.

The time for the fluorescence (that is, G_i + G_m) to decrease from peak value (measured at the brightest point of the Ca²⁺ spark) to half its value to the baseline is t_1/2. The peak ratio, or Ca²⁺ spark amplitude a, equals F/Fo where F is the peak fluorescence value and Fo is the baseline fluorescence value. The mean Ca²⁺ spark amplitude is given by F/Fo and the maximum ratio [occurring when (y, z) = (0, 0)] is F/Fo(max). The spatial spread of the Ca²⁺ spark at the time of the peak fluorescence is characterized by the FWHM. Because of the large variation in t_1/2 and FWHM, we also calculated these values [t_1/2(bright), FWHM(bright)] using only the 10 brightest Ca²⁺ sparks.

Typical t_1/2 values for Ca²⁺ sparks from heart cells is ~20 ms (Cheng et al., 1993), which is close to that found when  = 5 but not when  = 2. Note that the standard deviations are quite large, about half the mean t_1/2 value. The reason for this large variation is shown in Fig. 4 A where t_1/2 is plotted against the Ca²⁺ spark amplitude. The variation in t_1/2 is fairly small for the large amplitude Ca²⁺ sparks but is large for the low amplitude Ca²⁺ sparks because of noise.

View larger version (10K): [in this window] [in a new window]   FIGURE 4   Plots of Ca2+ spark decay time (t1/2) as a function of Ca2+ spark amplitude. Ca2+ sparks in simulated linescan images were identified with the Ca2+ spark detection program. All Ca2+ sparks were generated identically, so Ca2+ spark amplitude simply reflects distance between linescan and source. Thus decay time should reflect diffusion time within the sample volume and rise monotonically as the Ca2+ spark amplitude decreases, as shown in (B) obtained from noise-free linescans. The upward trend of decay time for decreasing amplitude is still evident when noise is present (A) but there is tremendous variability in the decay times when the Ca2+ spark amplitude is small.

Fig. 4 B shows a plot of t_1/2 against amplitude for the same set of simulations in Fig. 4 A but in the absence of noise. Since the Ca²⁺ sparks were generated identically, amplitude variations are due solely to variations in distance between linescan and Ca²⁺ spark origin. The decay time of identically generated Ca²⁺ sparks is controlled by the diffusion of Ca²⁺ into the scanned volume, so it increases with distance and, equivalently, decreases with Ca²⁺ spark amplitude. Viewed in isolation, Fig. 4 B suggests that the decay time could be used to distinguish whether a Ca²⁺ spark has a small amplitude because the linescan was far from the source or because the source strength was small. The results in Fig. 4 A cautions against such a method as virtually any decay time may be obtained for small amplitude Ca²⁺ sparks.

The mean Ca²⁺ spark amplitude is almost identical for  = 2 and 5 and is typical for experimentally measured Ca²⁺ sparks. F/Fo(max) values are also similar for the two values of , indicating that despite the larger amount of dye available when  = 5 the amount of Ca²⁺ released is sufficient to saturate the dye.

The FWHM values for the 10 brightest Ca²⁺ sparks is ~2 µm, which is about half the value reported by Gómez et al. (1996) for rat ventricular cells. Simulations carried out with longer open times or larger channel currents did not greatly alter the FWHM values.

Effect of changing microscope's FWHM

The triangles in Fig. 3 A showing the linescan positions at which the Ca²⁺ spark could be detected are symmetrically distributed around the origin, as expected since the axial and lateral FWHM values are equal. To study how this distribution changes when the blurring kernel is asymmetric, we increased the _z/_xy ratio to 3 (Fig. 3 B) and 6 (Fig. 3 C) where _xy was fixed to 0.2 µm. Note that the confocal parameters are different from those used to generate Fig. 3 A. The case where _z/_xy = 1 is not shown since the distribution of detected Ca²⁺ sparks is symmetric, as in Fig. 3 A. (Fig. 3, A-C may be interpreted in two equivalent ways: the point source is at the origin and the circles represent the linescan positions, or the linescan is fixed at the origin and the circles mark the point source locations. We take the latter viewpoint now so we can talk about the distribution of detected Ca²⁺ sparks instead of the more unwieldy distribution of linescan positions at which the Ca²⁺ spark was detected.) Fig. 3 B shows, surprisingly, that the distribution of detected Ca²⁺ sparks is still symmetric about the origin despite the axial FWHM being 3 times larger than the lateral FWHM. The distribution of detected Ca²⁺ sparks becomes asymmetric, however, when _z/_xy = 6, as shown in Fig. 3 C.

We were initially surprised to see the distribution in Fig. 3 B because we had expected to see an ellipsoidal distribution that parallels the elongation of the PSF along the z-axis. With increases in the depth of field (increasing _z) comes a loss in contrast of the Ca²⁺ spark, making it more difficult to detect the Ca²⁺ spark. This decrease in contrast with increases in _z is shown in Fig. 3 D (right to left) where the same Ca²⁺ spark is imaged with _z equaling 0.2, 0.6, and 1.2 µm, respectively. (No microscope to date has achieved an axial resolution of ~0.2 µm, but we have used this value for illustration.) The F/Fo values for these three cases are 2.19, 2.64, and 3.12 (left to right) yielding Ca²⁺ concentration values of 312, 447, and 664 nM (Cheng et al., 1993). Thus the simple act of opening the confocal pinhole, which increases both axial and lateral FWHM, can reduce the Ca²⁺ concentration estimates.

The physical reason for the decrease in contrast as _z increases is that because the light energy is spread over a larger volume, the intensity must be lower to maintain energy conservation. Mathematically, this constraint is expressed in the larger denominator (~_z) in the normalization factor of the Gaussian kernel in Eq. 11.

Ca²⁺ spark amplitude distribution

Ca²⁺ spark amplitude distributions obtained from simulated linescan images are shown in Fig. 5. All Ca²⁺ sparks were generated with a channel current of 1.4 pA and channel open time of 10 ms; only the linescan positions were varied randomly. In the absence of noise, panel A, the Ca²⁺ spark amplitude distribution decreases monotonically except for statistical sampling variations. [Nonmonotonicity in the Ca²⁺ spark amplitude distribution due to sampling variation can be distinguished from intrinsically multimodal distributions (for example when the SR Ca²⁺ release channels are arranged on a lattice, see below) by increasing the sample size or by changing the seed value of the random number generator. Intrinsically multimodal Ca²⁺ spark amplitude distributions are unaffected by these changes.]

View larger version (20K): [in this window] [in a new window]   FIGURE 5   Amplitude histograms and 1/fa(a) obtained using Ca2+ sparks in simulated linescan images identified by the Ca2+ spark detection program. (A) shows the Ca2+ spark amplitude distribution obtained from noise-free linescan images. All 176 Ca2+ sparks were detected by the Ca2+ spark detection program without false positives. (C) shows the plot of 1/fa(a) (squares) calculated using the bin values in (A). The theoretical line (solid) is virtually identical to the best-fit line (dashed). In (B) noise was added to the linescan images so the very dim Ca2+ sparks were not detected. (D) shows the 1/fa(a) values calculated from the bin values in (B) (squares); solid and dashed lines are as in (C). The maximum amplitude occurs when the confocal linescan goes through the spark origin. In noise-free images (A), amax = 2.62 while in noisy images (B), amax = 2.85.

These graphs illustrate the inherent difficulty in assessing the source strength distribution. Although all Ca²⁺ sparks in the linescan images were generated identically, because of the arbitrary placement of the linescan relative to the source, there is a broad distribution of measured Ca²⁺ spark amplitudes instead of a single narrow bin or narrow Gaussian distribution. A Gaussian distribution has been interpreted to indicate that Ca²⁺ sparks have stereotypic origins. However, Fig. 5 shows that, in our model, Ca²⁺ sparks generated identically do not generate a narrow Ca²⁺ spark amplitude distribution. This result is similar to that obtained by Pratusevich and Balke (1996).

One way that a monotonically decreasing Ca²⁺ spark amplitude distribution, Fig. 5 A, might be transformed into a Gaussian-like distribution is suggested by Fig. 5 B. In the presence of noise, Ca²⁺ sparks whose amplitude was <~1.2 were not detected by the detection program. Moreover, more Ca²⁺ sparks whose amplitudes were in the range 1.3-1.4 were detected than those Ca²⁺ sparks with amplitudes of 1.2-1.3. When noise is present the Ca²⁺ sparks of low amplitude are not detected with the same reliability as the large amplitude Ca²⁺ sparks. Thus, although there were actually more low amplitude Ca²⁺ sparks in the linescan images (Fig. 1 B) these Ca²⁺ sparks are masked by noise and appear to occur less frequently. The difference in reliability is quantified by a visibility function proposed by Pratusevich and Balke (1996). The sigmoidal visibility function gives the probability of detecting a Ca²⁺ spark of a given amplitude and ranges from 0 for amplitudes near 1 and rises to unity as the Ca²⁺ spark amplitude increases. The Ca²⁺ spark amplitude distribution that is measured is then the product of the "ideal" amplitude distribution, obtained by a perfect detector in the absence of noise (Fig. 5 A), and the visibility function. Multiplying an appropriately shaped visibility function with an amplitude distribution such as in Fig. 5 A can give a Ca²⁺ spark amplitude distribution that is Gaussian-like and similar to those reported in the literature (Klein et al., 1996; Shorofsky et al., 1996, 1997; Shirokova and Ríos, 1997; Xiao et al., 1997; Wier et al., 1997).

The key question is whether the intrinsic properties of the SR Ca²⁺ release channel, not detector characteristics, produce these experimentally measured Ca²⁺ spark amplitude distributions. To answer this question we need to establish the relationship between the Ca²⁺ spark amplitude distribution and the underlying source strength distribution.

Relationship between the Ca²⁺ spark amplitude distribution and source strength distribution

Let f_a(a) be the probability density function (pdf) of Ca²⁺ spark amplitudes. That is, the probability of finding a Ca²⁺ spark whose amplitude is between a  a/2 and a + a/2 is f_a(a)a. Likewise, let f() be the pdf of the source strength. The source strength  may refer to the SR Ca²⁺ release channel current for a fixed channel open time or to the open time for a fixed channel current.

To establish the link between f_a and f consider a simple and intuitive example. Suppose that a light bulb located at the origin flashes with intensity ' with probability p(') and flashes with intensity " with probability p(") = 1  p('). The light intensity a that an observer measures depends on his/her distance r from the lamp and the lamp intensity , and is given by the observation function, g(, r)

(13)

Suppose when the observer is at r' the lamp flashes with intensity ' and the observer measures intensity a₁ = g(r', '). If the observer moves randomly then the mean number of times that he/she measures an intensity a₁ is proportional to the probability of being at a distance r', p(r'), times the probability that the lamp flashed with intensity ', that is

(14)

The observer will also measure intensity a₁ when the lamp flashes with intensity " and his/her distance r" is adjusted accordingly to give a₁ = g(r", "). The appropriate distance is given by The function (a, ) can always be found provided the observation function g(, r) is a strictly monotonic function of r. Thus the probability p_a of measuring intensity a₁ becomes

(15)

To extend the argument to a continuum of source strengths let F_a(a) be the probability that the measured Ca²⁺ spark amplitude g(, r) < a. F_a(a) is the cumulative distribution function

(16)

where f_r(r) is the pdf of being at a distance r from the origin.

Although  and r are independent random variables, the values of  and a constrain the lower limit of integration of r. In order to satisfy g(, r) < a, the lower bound of r must be (, a). Thus,

(17)

The largest  compatible with a given a is given by _max = g¹(a, r). Thus Eq. 17 becomes

(18)

Differentiating F_a(a) yields the probability density function f_a(a)

(19)

This integral equation relating f_a(a) to f() is the main result. We now need to find specific forms of f, f_r, and . If the linescan can be at any position between 0  r  R with equal probability, then f_r(r) = 2r/R². Note that the use of r and not (y, z) comes from the implicit assumption that the blurring along the lateral dimensions x and y is the same as along the z-axis. Another assumption implicit in the use of r is that the diffusion is radially symmetric.

Explicit form for the observation function g(, r)

Because of the nonlinear buffer reactions, the observation function cannot be found analytically. We determined g(, r) empirically using the following procedure. Linescan images (100-200) containing a total of ~150-300 Ca²⁺ sparks were generated with a set of parameters for the reaction-diffusion simulations and a channel current of , say 1.4 pA, and a fixed channel open time (10 ms). Ca²⁺ sparks were found using the Ca²⁺ spark detection program and their amplitudes (a = F/Fo) calculated. Since the (y, z) coordinates of each linescan were known, the amplitude at the distance r = (y² + z²)^1/2 could be calculated. The pairs of (a, r) were fit to the function

(20)

This procedure was repeated for different channel currents to determine A(), B(), and C().

The observation functions for four different channel currents are shown in Fig. 6 A. The solid curve shows the best fit to the data and, for clarity, data points are only shown for  = 0.7 pA and  = 2.8 pA. A() was fit to the hyperbolic function

(21)

and C() to the line

(22)

shown in Figs. 6 B and 6 C. No theoretical significance is attached to the specific forms of A() and C(); they were simply chosen for simplicity. B() was essentially independent of  varying between 1.20 and 1.25. This is expected since B should only reflect the sensitivity of the detection program and the amount of added noise. Moderate (3-5-fold) reduction of k_m⁺, k_m, k_i⁺, and k_i values did not affect the functional form of the observation function or of A() and C().

View larger version (17K): [in this window] [in a new window]   FIGURE 6   Determining the observation function for variable channel current. Ca2+ sparks from simulated linescan images were identified with the Ca2+ spark detection program. (A) shows the measured amplitude as a function of distance between the linescan and source for each identified Ca2+ spark (symbols). The data points were fit to the function A() exp{[C()r]2} + B(), where  is the channel current. For clarity only the data points for  = 0.7 (circles) and 2.8 (triangles) pA are shown. Intermediate curves are for  = 1.4 and 2.1 pA. (B) and (C) show the fit parameters A() and C() from (A) as functions of .

In the next section we will derive a specific relationship between f and f_a that will allow us to examine the effects that different source strength distributions have on the Ca²⁺ spark amplitude distribution. The specific relationship between f and f_a depends, of course, on our assumption that the observation function is Gaussian. Different observation functions yield different relationships between f and f_a. Thus, it is worthwhile to examine the range of conditions under which the observation function is likely to be Gaussian. When Ca²⁺ release comes from a point source and the source strength is sufficiently weak so that the dye does not saturate, then the Ca²⁺ bound fluo-3 distribution is approximately Gaussian. We assumed that the PSF is Gaussian, which well approximates the actual PSF for a correctly aligned confocal microscope with a fairly small pinhole. The convolution of the Gaussian Ca²⁺ bound fluo-3 concentration profile with the Gaussian PSF gives a Gaussian image; the observation function is the profile of this convolution. We note that the spatial profiles of many Ca²⁺ sparks are approximately Gaussian (Parker et al., 1996; Gómez et al., 1996).

The observation function will deviate from a Gaussian when the source is extended [see Smith et al. (1998) for a discussion of extended sources], when the dye is saturated, or the confocal microscope is poorly aligned. Under these conditions the observation function must be amended. Later we will see the effect of dye saturation on f_a.

Explicit relationship between f and f_a

From Eq. 20 it follows that (a, r) is given by

(23)

and

(24)

Now suppose that all Ca²⁺ sparks are generated identically; that is, there are no variations in the source strength then the source strength pdf is f() = (  _o), where  is the Dirac delta-function. In this case Eq. 19 becomes

(25)

The smallest amplitude that can be attained with this  is a_min(_o) = g(_o, R) and the largest is a_max(_o) = g(_o, 0). Using Eqs. 23, 24, and f_r() = 2/R², the explicit expression for f_a(a) in Eq. 25 is

(26)

The difference of the Heaviside functions, H, limits f_a(a) to a_min < a < a_max. [The Heaviside function H(a  x) is a step function that equals 1 for x  a and 0 otherwise.] Between these limits f_a(a) ~ (a  B)¹.

Equation 26 is one of the key results of this paper. It implies that if all Ca²⁺ sparks were generated identically and if the observation function were Gaussian (Eq. 20), then the resulting Ca²⁺ spark amplitude histogram as measured by confocal microscopy should be hyperbolic, not Gaussian. Accordingly, a plot of 1/f_a(a) against a yields a straight line.

Relationship between f_a and the Ca²⁺ spark amplitude histogram N(a)

Let N(a) be the number of Ca²⁺ sparks having amplitudes between a  /2  a  a + /2, where  is the binwidth. Then

(27)

where N_total is the total number of Ca²⁺ sparks. Equation 27 can be turned around to get an estimate of f_a, f_a^est,

(28)

We can now compare the theoretical curve f_a(a) given by Eq. 26 to that given by Eq. 28. In Fig. 5 C we have plotted 1/f_a(a) = N_total/N(a) (squares), where N(a) is the data from Fig. 5 A, N_total = 176 Ca²⁺ sparks, and  = 0.1. The solid line is the theoretical f_a calculated using Eq. 26 with C( = 1.4) = 2.86 and R = R₉₀ = 0.80 µm; this line is the best descriptor of the data points as it is virtually coincident with the best fit line (dashed line). This agreement between simulation and theoretical results is important because it provides a check on the derivation of the relationship between f_a and f. Thus we can simulate the distribution of Ca²⁺ spark amplitudes in a new way. Instead of making linescan images, detecting the Ca²⁺ sparks, and then calculating their amplitudes, we used the following method. The confocal linescan position was chosen randomly in the y-z plane and its distance r from the Ca²⁺ spark at the origin was calculated. The amplitude was then calculated using the observation function. With this new method we could simulate conditions that would be extremely tedious or impossible by the old method.

Estimating  from the Ca²⁺ spark amplitude histogram

In this instance  was known so the theoretical line could be calculated. In practice  is unknown but can be calculated from the information available in the Ca²⁺ spark amplitude histogram as follows. If the plot 1/N(a) against a falls on a single straight line then the data are consistent with a delta-function source strength pdf, (  _o). (See below for f_a when f is more complicated than a single delta-function.) _o is calculated using the largest measured Ca²⁺ spark amplitude using Eqs. 20 (with r = 0) and 21. In this case a_max = 2.85, which gives _o = 2.0, precisely the value used in the simulations. Having calculated _o, R can be calculated for each a using Eqs. 22 and 20. The calculated values will naturally depend on the simulation parameters such as the amount of buffer available and their kinetics of reaction with Ca²⁺.

f_a of more complicated f

Suppose instead of f being a single Dirac -function, f is the weighted sum of -functions

(29)

where _i gives the probability of the source strength being _i so the  values satisfy _ii = 1. Since Eq. 26 holds for all  it follows that

(30)

Since a_min() and a_max() are increasing functions of , f_a is the sum of terms (a  B)