This Article Abstract Full Text (PDF) All Versions of this Article: biophysj.104.045260v1 88/4/2403    most recent 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 Ventura, A. C. Articles by Ponce Dawson, S. Search for Related Content PubMed PubMed Citation Articles by Ventura, A. C. Articles by Ponce Dawson, S.

A Model-Independent Algorithm to Derive Ca²⁺ Fluxes Underlying Local Cytosolic Ca²⁺ Transients

Alejandra C. Ventura ^*, Luciana Bruno ^*, Angelo Demuro , Ian Parker  and Silvina Ponce Dawson ^* 

^* Departamento de Física, Facultad de Ciencias Exactas y Naturales, Universidad de Buenos Aires, Buenos Aires, Argentina; Department of Neurobiology and Behavior, University of California, Irvine, California; and T10-Theoretical Biology and Biophysics Group, Los Alamos National Laboratory, Los Alamos, New Mexico

Correspondence: Address reprint requests to Silvina M. Ponce Dawson, PhD, Tel.: 54-11-4576-3353; E-mail: silvina{at}df.uba.ar.

	ABSTRACT

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS APPENDIX 1: MOTIVATION FOR... APPENDIX 2: DESCRIPTION OF... ACKNOWLEDGEMENTS REFERENCES

 
Local intracellular Ca²⁺ signals result from Ca²⁺ flux into the cytosol through individual channels or clusters of channels. To gain a mechanistic understanding of these events we need to know the magnitude and spatial distribution of the underlying Ca²⁺ flux. However, this is difficult to infer from fluorescence Ca²⁺ images because the distribution of Ca²⁺-bound dye is affected by poorly characterized processes including diffusion of Ca²⁺ ions, their binding to mobile and immobile buffers, and sequestration by Ca²⁺ pumps. Several methods have previously been proposed to derive Ca²⁺ flux from fluorescence images, but all require explicit knowledge or assumptions regarding these processes. We now present a novel algorithm that requires few assumptions and is largely model-independent. By testing the algorithm with both numerically generated image data and experimental images of sparklets resulting from Ca²⁺ flux through individual voltage-gated channels, we show that it satisfactorily reconstructs the magnitude and time course of the underlying Ca²⁺ currents.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS APPENDIX 1: MOTIVATION FOR... APPENDIX 2: DESCRIPTION OF... ACKNOWLEDGEMENTS REFERENCES

 
The liberation of Ca²⁺ ions from intracellular stores into the cytosol is used as a signaling mechanism by virtually all cell types to regulate different functions. Recent advances in imaging techniques—including enhanced fluorescent Ca²⁺ indicators and confocal imaging systems with improved temporal (1 ms) and spatial (0.5 µm) resolution (Pawley, 1995)—provide a powerful tool to study these signals. They have revealed the existence of a wide range of events involving the release from single channels, clusters of channels, and global waves (Cheng et al., 1993; Yao et al., 1995). A complete description of these signals requires a detailed knowledge of the magnitude and kinetics of the underlying Ca²⁺ flux or current. These cannot readily be inferred from fluorescence Ca²⁺ images. In the first place, the images reflect Ca²⁺ ions bound to the fluorescent indicator. Secondly, but more difficult to solve, Ca²⁺ ions do not diffuse freely upon release: endogenous buffers and pumps affect their spatiotemporal dynamics in ways that are hard to determine.

Various methods for estimating Ca²⁺ release fluxes have been presented in the literature, especially in the case of sparks in cardiac, skeletal, and smooth muscle. As described in Rios and Brum (2002), they can be classified in two categories:

1. Forward calculations, in which a model of intracellular Ca²⁺ dynamics is numerically simulated, and the fluorescent image that such a model produces is compared to the experimental observation (see, e.g., Smith et al., 1998; Baylor et al., 2002; Izu et al., 2001).
   
2. Backward calculations, in which, starting from a given image, the underlying Ca²⁺ current is computed (see, e.g., Blatter et al., 1997; Rios et al., 1999; Rios and Brum, 2002; Soeller and Cannell, 2002).
   

Regardless of the approach, all of these methods require a working model of cytosolic Ca²⁺ dynamics. Moreover, quantitative data including the kinetics and concentrations of endogenous buffers and pumps are needed (Smith et al., 1998; Izu et al., 2001; Rios and Brum, 2002; Soeller and Cannell, 2002; Blatter et al., 1997; Lukyanenko et al., 1998; Rios et al., 1999), and have been obtained either from experimental studies or by adjustment of parameters within the model to provide the best agreement with the observations (Blatter et al., 1997; Timmer et al., 1998). In most cell types, however, knowledge of their intracellular Ca²⁺ buffers is much less complete than in well-characterized muscle cells, so that these approaches are not generally applicable.

We propose here a novel, model-independent algorithm to obtain Ca²⁺ release fluxes from intracellular fluorescence images, which is particularly suitable for images obtained by optical single ion channel experiments. The method is applicable when Ca²⁺ release is clearly localized in space and time. It only needs a priori information on the kinetics and diffusion parameters of the dye and relies on geometric assumptions that allow us to extrapolate line-scan information into information on the whole space. It does not need, however, any a priori information on the endogenous buffers, pumps, or leaks present in the experimental system. The underlying concept is that the spatiotemporal series of the Ca²⁺ concentration inherently carries information of how the various processes at work affect its evolution, an idea that is drawn from the theory of time-series analysis (Packard et al., 1980; Mindlin et al., 1990, 1998; Gouesbet, 1992). In this article we present a detailed description of the algorithm, and verify its accuracy by testing it with both numerically generated images and experimental confocal fluorescence images of Ca²⁺ flux through single voltage-gated channels expressed in Xenopus oocytes.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS APPENDIX 1: MOTIVATION FOR... APPENDIX 2: DESCRIPTION OF... ACKNOWLEDGEMENTS REFERENCES

 
The algorithm: basic assumptions
Fluorescence images provide information on the distribution of Ca²⁺-bound dye, [CaB]. Given that we are interested in images obtained in oocytes, we assume that the dye only reacts with Ca²⁺ according to the kinetic scheme

(1)

where k_on and k_off are the kinetic rate constants from which the dissociation constant, K_d k_off/k_on, is obtained. We assume that the dye, B, is initially distributed uniformly in space, that there are no dye sources or sinks during the time course of the experiment, and that its diffusion coefficient is not changed when Ca²⁺ binds to it. Then, the total dye concentration, [B]_T = [B] + [CaB], remains equal to its uniform (known) initial value at any given time and spatial point. Thus, the evolution of [CaB] is given by

(2)

where D_B is the dye diffusion coefficient and R is the reaction term,

(3)

We can compute the time and spatial derivatives of [CaB] from the experimental records. In this section we will assume that they can be computed with perfect accuracy and discuss their errors later. Now, line-scan images only contain data along one spatial dimension, x. We then assume that the release of Ca²⁺ has spherical symmetry (this is satisfied, in particular, if the release occurs over a very small region that can be treated as a point), and that there is nothing else around that can break this symmetry. With these assumptions, the solution always has spherical symmetry (i.e., it only depends on time and the distance, r, from the site of release). If the line-scan goes through the center of the release region, we can compute ²[CaB] using only the information along the direction, x. We analyze later what happens when the line-scan does not go exactly through the center of the release site. Calculating ²[CaB] and [CaB]/t we can obtain from Eq. 2 the value of R at every spatial point and time. Using Eq. 3, we can then obtain a spatiotemporal series of the variable of interest, [Ca²⁺]. Up to this point, our method is similar to previous ones (Blatter et al., 1997). It is also similar to the method proposed recently in Soeller and Cannell (2002), although it differs in the fact that here we assume that the dye only binds to Ca²⁺.

Upon its release Ca²⁺ diffuses, reacts with buffers, and is recaptured by different pumps. Of all these processes, we want to single out three—two for which we know some details, the reaction with the dye and diffusion; and the other that we want to compute, the Ca²⁺ source, Q_Ca. In this way we write the evolution equation for Ca²⁺ as

(4)

where D_Ca is the free diffusion coefficient of Ca²⁺ (220 µm²/s, Allbritton et al., 1992) and the term M represents the sum of all the processes whose details are poorly known in each experiment, which usually depend on the peculiarities of the cell that is being observed. If we had a complete model of the intracellular Ca²⁺ dynamics, then we would be able to write an expression for M, which would depend on other variables, such as the concentrations of the endogenous Ca²⁺ buffers and of the pumps (this is what is done in Blatter et al., 1997; Rios and Brum, 2002; Rios et al., 1999; Soeller and Cannell, 2002). Therefore, Eq. 4 would have to be solved coupled to evolution equations for these other variables (or, equivalently, as done in Blatter et al., 1997; Rios and Brum, 2002; Rios et al., 1999; Soeller and Cannell, 2002, once [Ca²⁺] is known as a function of space and time, the evolution equations for these other variables could be integrated, provided that we knew all the kinetic parameters and initial concentrations that, together with [Ca²⁺], determine their evolution). Unfortunately, when the information about the buffers, pumps, or leaks in the experimental condition is insufficient, it is not possible to write an accurate expression for M without introducing ad hoc assumptions. Even when this can be done, the model validation itself can be a difficult (if not impossible) task and usually, the model lacks universality.

As shown in Appendix 1, all the dependence of M(r, t) on variables other than [Ca²⁺] can be recast in terms of the values of [Ca²⁺] at (r, t) and at other spatial points and previous times, (r', t'). Due to the dissipative and diffusive nature of the processes involved in M, the dependence with [Ca²⁺](r', t') becomes weaker as (r', t') gets further away from (r, t). Performing a series expansion around (r, t), we can rewrite M in terms of [Ca²⁺](r, t) and space and time derivatives of [Ca²⁺] of relatively low order. The distinguishing feature of our method is the assumption that M can be approximated by this type of ansatz and that its functional dependence can readily be determined from the data. Namely, we assume that M can be written as

(5)

where g, f, h, and k are unknown functions of the free Ca²⁺ concentration, [Ca²⁺], to be determined by the algorithm. We see from Eq. 5 that there is no term of the form v · [Ca²⁺] (a term proportional to a first derivative with respect to the spatial coordinates). This is so because there is no privileged direction in the problem due to the assumption of the spherical symmetry. Thus, we have gone to second order in space, keeping those terms that respect the spherical symmetry.

Even if we insert Eq. 5 in Eq. 4, Eq. 4 still has the unknown that we want to determine: the Ca²⁺ source, Q_Ca. The fact that Ca²⁺ release is localized in space and time implies that Q_Ca vanishes away from the source and everywhere in space when Ca²⁺ release ceases. Thus, there are points of the experimental record where Q_Ca = 0 and M is then the only unknown term of Eq. 4. Using the ansatz (Eq. 5), we can obtain numerical approximations to the unknown functions g, f, h, and k from the points of the record where Q_Ca = 0. Once we have these functions, we then know M as a function of [Ca²⁺] and its derivatives. Therefore, the only unknown of Eq. 4 is now Q_ca. So, computing the various space and time derivatives that are needed at the points where Q_Ca 0, we can obtain Q_ca from Eq. 4 simply by doing Q_Ca = – R – D_Ca²[Ca²⁺] – M + [Ca²⁺]/t. Implicit in this method is the assumption that the way in which the other processes affect the spatiotemporal behavior of [Ca²⁺] can be written uniformly throughout the cytosol as a function of [Ca²⁺] and its derivatives. This is equivalent to say that the distribution of pumps and the total buffer concentrations are uniform in the cytosol. These assumptions are always present in the backward methods used to obtain the Ca²⁺ current from fluorescent images.

The algorithm: numerical implementation
We now describe the algorithm we have developed following the ideas outlined in the previous section. It has six main steps. Fig. 1 shows an overview of the procedure. The images of Fig. 1 correspond to numerically generated experiments obtained by simulating Eqs. 15–17 with a source of the form (Eq. 19), with R_f = 0.15 µm, t_s = 3 ms, and t_d = 10 ms. The experiment on the top has I = 0.1 pA and the one on the bottom I = 0.3 pA.

View larger version (57K): [in this window] [in a new window]   FIGURE 1  Outline of the main steps involved in the flux reconstruction algorithm. (A) Space-time plot of the fluorescence distribution of the two (numerically generated) experiments with which we are going to work. Both experiments are needed for the application of the algorithm but the calcium release flux will only be determined for the one on the top. (B) [CaB] is shown on the same space-time plot as in A. [CaB] is obtained from the fluorescence data using Eq. 6. (C) [Ca2+] is shown on the space-time plot of A. [Ca2+] is obtained from [CaB] using Eq. 8. (D) Same figure as in C but with a grid superimposed onto the regions where we assume that QCa = 0. For illustration purposes, the grid that is drawn is much more sparse than the one that we actually use for our calculations (which would be undistinguishable on the scale of the figure). We compute M at each point in the gridded regions using Eq. 9. (E) In Step 4, the values of [Ca2+] are binned according to Eq. 10 (in the example we use N = 50). We show in the figure the number of data points, (ri, ti), in the gridded regions of D as a function of Ci. At least four data points are necessary for each Ci to apply the method. (F) The functions g (dimensionless), f (in units of µm2/s), h (in units of µm2/(µM/s)), and k (in units of µM/s), as a function of [Ca2+], obtained by interpolating the solution of the overdetermined system of the expressions in Eq. 11. Circles correspond to the solutions that minimize the merit function (Eq. 35) at each binned [Ca2+] value (see Appendix 2). The smooth curves correspond to fits to the circles using fifth-degree polynomials. Knowing the smooth functions, we have an expression for M as a function of [Ca2+] and its derivatives via Eq. 5, for all the values of [Ca2+] that the experiment on the top of A attains. Inserting this expression in Eq. 4, we obtain the source distribution, QCa, in space and time for this experiment as shown in G.

(6)

with F(x_j, t_k) the fluorescence at pixel (j, k), of coordinates x_j = jx and t_k = kt, where x and t are the spatial and temporal resolutions of the experiment, respectively, and the spatial coordinate x is the distance along the line-scan. F_min and F_max are the minimum and maximum values of the fluorescence corresponding to Ca²⁺-free and Ca²⁺-saturated dye. The ratio F_max/F_min can be determined, a priori, for each experimental condition. The meaningful information is contained in all the data points with either x > 0 or x < 0, given that the point x = 0 coincides with the center of the source. Thus, there is redundant information if the spherical symmetry assumption is correct. We then average the information at each point x > 0 with the information from the corresponding negative value (x < 0) and work with a half-signal image.

Fig. 1 A shows the spatiotemporal evolution of the fluorescence, as typically obtained in line-scan images, for the two numerical experiments. Fig.1 B displays the distribution, [CaB], computed from Eq. 6 using the fluorescent distribution of Fig. 1 A after taking the average of F at plus and minus x, for each x. We use a pseudo-color scale to quantify concentration values: red corresponds to the highest concentration whereas blue corresponds to the basal one.

Step 2: obtaining the free calcium concentration
Once having [CaB], we obtain the reaction term R from Eq. 2,

(7)

[Ca²⁺](r, t) is computed using Eqs. 3 and 7,

(8)

Notice that Eq. 8 contains a term of the form [B]_T – [CaB] in the denominator. Therefore, the method will work provided that this term is away from zero, a condition that is met when the dye is not saturated. That is exactly the condition under which imaging experiments are done: dye saturation is avoided to follow [Ca²⁺] changes as closely as possible (e.g., compare the maximum value of [CaB] in Figs. 1 and 7 with the total dye concentration, [B]_T = 40 µM). Also notice that obtaining R and [Ca²⁺] involves taking numerical derivatives with respect to time and space. As numerical differentiation is very sensitive to noise in the input data, filters might be used to minimize noise amplification. As explained before, for the calculation of the Laplacian we assume spherical symmetry around the center of the source.

View larger version (46K): [in this window] [in a new window]   FIGURE 7  Experimental sparklets. (A) Line-scan image illustrating single channel events (sparklets) that arise in response to a 3-s-long depolarizing pulse. (B) Spatiotemporal spreading of the calcium-bound dye for the events in the box indicated in A. The scale bars in B are also valid for C, E, and G. (C) Spatiotemporal, M. (D) The curve M = k([Ca2+]) and its smooth fitting. (E) Source obtained by the algorithm. (F) Temporal profiles of the reconstructed current that the algorithm predicts obtained by volume integration (Eq. 13) over a semispherical region of radius r 0.225 µm (solid line) and of the one derived from the simplified, piecewise constant source, QCa, used in the forward integration of the evolution equations (dashed line). (G) Spatiotemporal distribution of the Ca2+-bound dye complex, [CaB]s, obtained by numerical integration of Eq. 4 with M(Ca2+) given by the smooth curve shown in D and the simplified source from which the dashed curve of F is obtained. (H) Temporal profile of the reconstructed current (solid line) and of the measured fluorescence averaged over a 0.45-µm region around the center of release (dashed line). For clarity, the error of the reconstructed current (27%) is not shown.

(9)

In the case of Fig. 1 we know that the Ca²⁺ source spreads over a region of radius 0.15 µm; it turns on at t = 3 ms and lasts for 10 ms. However, for experimental data, mild a priori assumptions about the spatiotemporal distribution of the release source must be made at this point, since the precise morphology of the source is not generally known. In these cases, we discard a region of the spatiotemporal record that overestimates the typical source size and duration, trying to keep enough data points in the gridded regions with concentrations different from basal. If there are N_tot data points within the gridded regions, we have N_tot equalities of the form (Eq. 9) from which we can obtain M as a function of (r, t) for values of r and t such that Q_Ca(r, t) = 0. As explained in the next step, N_tot needs to be relatively large for the algorithm to work.

Step 4: binning the values of [Ca²⁺] attained during a series of experiments
The previous step provides M as a function of (r, t) for values of r and t such that Q_Ca(r, t) = 0. However, we want M as a function of [Ca²⁺] and its derivatives so that we can use Eq. 5, the ansatz, to determine the four unknown functions of [Ca²⁺], f, g, h, and k. To this end, we compute the minimum, [Ca²⁺]_min, and maximum, [Ca²⁺]_max, values of [Ca²⁺], attained in regions with Q_Ca = 0 by the whole set of experiments that we are going to use (all data points in the gridded regions of Fig. 1 D). Then, we bin the interval [[Ca²⁺]_min, [Ca²⁺]_max] in a finite set of values. Thus, we define a vector of free Ca²⁺ concentration values as

(10)

where C ([Ca²⁺]_max – [Ca²⁺]_min)/n and n is the number of bins of [Ca²⁺] values that we will consider. We then determine all the points, _{j, k}(C_i) (r_j, t_k)(C_i), in the gridded regions of the images such that C_i [Ca²⁺](r_j, t_k)(C_i) < C_i+1. Clearly, there is more than one data point, _{j, k}(C_i), for each value of C_i. This is illustrated in Fig. 1 E where we plot the total number of grid points, N (coming from all the experiments that are being used), whose [Ca²⁺] value is contained in the bin of concentration, C. As C increases, N decreases significantly. In particular, N(C) 4, to be able to apply the algorithm. Therefore, all values of C such that N(C) < 4 are discarded. The arrow in Fig. 1 E indicates this cutoff and the inset shows a zoom of the figure around the cutoff value in this example. As described later in further detail, the value of C must be chosen so as to guarantee the applicability of the algorithm for a sensible range of [Ca²⁺] values.

Step 5: obtaining the functions g, f, h, and k
The binning of [Ca²⁺] reduces the functions f, g, h, and k to vectors defined as f_i = f(C_i), g_i = g(C_i), h_i = h(C_i), and k_i = k(C_i), 1 i n. Therefore, we now seek approximations to these vectors. For each value of C_i, we have as many binned versions of Eq. 5 of the form

(11)

as data points _{j, k}(C_i). The only unknowns in Eq. 11 are f_i, g_i, h_i, and k_i, 1 i n. Thus, to determine their values, we need at least four data points, _{j, k}(C_i) for each i, 1 i n. This imposes some constraints on the allowed values of C. We choose C so that the system (Eq. 11) is overdetermined (the number of data points, _{j, k}(C_i) 4, for every i). In this way, f_i, g_i, h_i, and k_i may be found using standard minimization or singular value decomposition techniques (see Appendix 2). Once g_i, f_i, h_i, and k_i are found, a global fitting of each function dependence on [Ca²⁺] is performed (Fig. 1 F). To this end, we use a feedforward neural network which provides a nonlinear curve fitting in terms of hyperbolic tangents. This allows us to interpolate the value of M at any Ca²⁺ concentration within the range [C_min, C_max]. As explained in Results, below, when working with real experimental records, only a small proportion of the records has [Ca²⁺] values that lie beyond C_max and are used only for the estimation of the functions f, g, h, and k.

Step 6: calcium release reconstruction
Finally, once the functions g, f, h, and k are known, M can be computed at any spatial point and time using Eq. 5, even where and while there is release, since both [Ca²⁺] and its space and time derivatives are known at every data point, (r, t). Now, in a given experiment, [Ca²⁺] will be larger in the region with Q_Ca 0 than in the region where there is no source. However, the collection of M values as a function of [Ca²⁺] and its derivatives is obtained in the region with Q_Ca = 0. If we use the data from one experiment with only one release event, we will not have direct information on the values of M for the larger [Ca²⁺] values that occur in the region with Q_Ca 0. One possibility to overcome this difficulty is to extrapolate the functions, f, g, h, and k, to obtain estimates of their values for these larger [Ca²⁺], using the expressions obtained for the lower concentrations. The other possibility is to determine these functions using various spatiotemporal regions, coming either from around different release sites or from different experiments done on the same cell (as illustrated in Fig. 1 A, where the upper image corresponds to a low influx Ca²⁺ current and the bottom image to a higher current one). In this way, the information that comes from the region around the site with the largest release will be used only to determine g, f, h, and k. The algorithm will then be applied to obtain the release flux from the regions with smaller release. We mainly follow this last approach. In this way, we cannot obtain an estimate of what the underlying Ca²⁺ current is at the site of largest release. As discussed in Appendix 2, we obtained better current estimates setting f = g = h = 0 from the beginning and using only k 0. All figures shown in the article, with the exception of Fig. 1, were done in this way.

Once we have M at every (r, t) for a particular experiment we can go back to Eq. 4 and solve for Q_Ca,

(12)

Eq. 12 gives the source distribution in space and time in units of concentration per unit time, e.g., µM/s (Fig. 1 G). Integrating Q_Ca over the volume where it is non-negligible we can also obtain the number of Ca²⁺ ions that enter the region per unit time which, in turn, can be expressed as an electrical current, I, using the known information on the ion's charge. If Q_Ca(r, t) is in µM/s, then the current could be obtained, in pA, as

(13)

where = 1.92 x 10^–4 pC/µm³ and we choose r_source(t) so that

(14)

with = 0.5, and = 0.3.

Numerical simulations
Most of the numerical simulations that we use to test the accuracy of the reconstruction method correspond to a model with Ca²⁺, dye, and one extra buffer. The dye, B, and the extra buffer, E, react with Ca²⁺ according to a kinetic scheme like the one in Eq. 27 with reaction rates of and respectively. Assuming that both the dye and the buffer are initially uniformly distributed in space and that the Ca²⁺ bound and the free corresponding buffer diffuse at equal rates, then the evolution equations are

(15)

(16)

(17)

where [B]_T and [E]_T are the dye and buffer total concentrations, D_Ca, D_CaB, and D_CaE are the diffusion coefficients of free Ca²⁺, dye, and extra buffer, respectively, and q_uptake represents the Ca²⁺ uptake, which is given by (Soeller and Cannell, 2002)

(18)

All parameters of the model are listed in Table 1 (the parameters of the dye B correspond to those of fluo-4 dextran, and the buffer E to those of EGTA). For some numerical experiments we take k₁ = 0. The set of expressions in Eqs. 15–17 is not a realistic model of Ca²⁺ dynamics inside a cell, but is useful to test the performance of the algorithm. For the Ca²⁺ release term, Q_Ca, we assume that it is different from zero on a sphere of radius R_f (0.15 µm R_f 1.0 µm). The Ca²⁺ release time course is either a square step

(19)

or a step with an exponential tail,

(20)

where so that Q_ca is in µM/s and I in pA. All these simulations were started from a uniform equilibrium distribution of all species at a resting Ca²⁺ concentration of 50 nM. The model equations were solved numerically using a finite difference scheme on a spherically symmetric grid, with spatial and temporal steps of 0.01 µm and 10^–7 s, respectively.

(21)

where (w = xy, z). In some cases, to analyze the effects of the blurring process by itself, we applied Eq. 21 to each data point of the simulation (i.e., using d equal to the spatial resolution of the simulation instead of 150 nm). In certain cases, we also applied a deblurring procedure to the blurred data using the inverse of Eq. 21. We also analyzed the effect of noise of realistic amplitude on the application of the algorithm to numerical simulations. Namely, numerical simulations were performed as described before. The resulting data were subsequently blurred according to Eq. 21, keeping points with the space (d = 0.15 µm) and time (t = 8 ms) resolutions of the experiments. Noise with normal distribution of standard deviation 0.12 in units of the rescaled fluorescence (see Eq. 25) was then added to the blurred data. This amount of noise is similar to the one observed in the experimental records that we analyze in this article.

The assumption that the signal has spherical symmetry and that the scanning line, x, corresponds to a radial direction, holds only if the scanning line goes through the center of the release site. Although experimentally the focal plane is adjusted so that the release site belongs to it, there are always small errors and a certain offset, d_o, can be expected. To study how the algorithm performs if the scanning line does not go through the center of the release site, we took the results of the (spherically symmetric) simulations at a chosen distance or offset, d_o, from the center of the source to represent out-of-focus data. We then applied the algorithm to these out-of-focus numerically generated data.

Finally, we also performed some tests on a more realistic model of intracellular Ca²⁺ dynamics, more specifically, on a modified version of the model presented in Baylor et al. (2002) to describe Ca²⁺ sparks in frog skeletal muscle. The modifications we introduced are related to the way we deal with the dye and with the Ca²⁺ pump. First, we assume that the dye only binds Ca²⁺. As it currently stands, our method cannot handle situations in which significant amounts of dye bind to proteins and fluoresce. Second, we use Eq. 18 for the Ca²⁺ uptake into internal stores (Soeller and Cannell, 2002). Other than that, we keep all the reactions considered in Baylor et al. (2002), as described in Table 2, and use a source of the form (Eq. 19) with R_f = 0.15 µm, t_d = 10 ms, and different values of I. The main goal of these simulations is to test the performance of the algorithm when several buffers of various kinetics are present and this modified model serves for this purpose.

(22)

using the source predicted by the algorithm and the one used in the simulation. In the case of sources of the form (Eq. 19), we also compute the time average

(23)

using the source obtained with the algorithm and the one used in the simulation. This provides information on how well the algorithm determines the spatial spread of the source.

Experiments
Experiments were performed on defolliculated Xenopus laevis oocytes expressing N-type voltage-activated Ca²⁺ channel from _1B–d and ß₃-subunits as previously described in Demuro and Parker (2003). Briefly, 2–4 days after mRNA injection, oocytes were injected with fluo-4 dextran (final intracellular concentration 40 µM) and 1 h later were placed in a chamber bathed in Ringer's solution including 6 mM calcium. Fluorescence signal was detected through a confocal pinhole providing lateral and axial resolutions of 300 and 500 nm, respectively. Experiments were performed in line-scan mode (0.15 µm/pixel; 8 ms/line) and near-membrane images were constructed and analyzed using routines written in the IDL programming language. Membrane voltage was normally clamped at –60 mV using a conventional two-electrode voltage-clamp technique, and transient Ca²⁺ signals were evoked by stepping the oocyte membrane to more positive voltage than –25 mV. Strong polarization up to +30 mV evoked spatially homogeneous fluorescence signals that progressively become more localized and stepwise during weaker voltage steps approaching –25 mV. On the evidences presented in Demuro and Parker (2003) we consider these pulsatile signals to rise from Ca²⁺ influx through a single N-type Ca²⁺ channel and refer to them as sparklets.

All data presented here were obtained using fluo-4 dextran as the Ca²⁺ indicator. This low affinity dye showed a low fluorescence signal at basal cytosolic free [Ca²⁺] and a large increase in fluorescence signal in the presence of high Ca²⁺ concentration, such that local fluorescence signals (F/F) of fivefold or greater could be produced. The parameters we considered for fluo-4 dextran (low affinity, K_d 4 µM) are: F_max/F_min 20 in the oocyte, k_on = 100 µM^–1 s^–1, and k_off = 400 s^–1, whereas D = 50 µm² s^–1. Dye concentration in the oocyte was 40 µM and resting Ca²⁺ concentration was considered to be [Ca²⁺]_rest 50 nM.

Processing of experimental signals
In the case of the experimental signals, the starting point (unless otherwise noted) is the smoothed-out fluorescent distribution that is obtained by averaging the raw one at each point of the space-time data file with its eight adjacent neighbors. We then compute the (theoretical) average basal fluorescence, F_b, inverting Eq. 6, using the basal concentration [CaB]_b = [Ca²⁺]_b[B_T]/(K_d+[Ca²⁺]_b) that corresponds to [Ca²⁺]_b = 50 nM, K_d = 4 µM, and [B_T] = 40 µM. This gives

(24)

We then compute the mean resting fluorescence at position x_j on the scan line, as the time average of the fluorescence value at spatial point, x_j, before there is Ca²⁺ release. In the absence of inhomogeneities, and assuming that [Ca²⁺]_b = 50 nM before there is Ca²⁺ release, it should be and this should be independent of (x_j); i.e., In any real experiment, this equality is only approximate. Thus, inserting in Eq. 24, we obtain

(25)

We use this approximation to compute [CaB] from Eq. 6. The ratio in Eq. 25 is the rescaled fluorescence, F/F₀, which we display in the plots of the experimental signals.

To test to what extent the results of the algorithm are affected by the presence of noise, we smoothed out the experimental images even further using a nonlinear function fitting obtained with a feedforward neural network of the MATLAB Neural Network Toolbox (The MathWorks, Natick, MA). The parameters of the fitting, which can be written in terms of hyperbolic tangents, minimize the mean-square error between the smooth function and the data points. We compared the results obtained starting from these smoother images with those obtained starting from an image in which only the eight-neighboring preprocessing had been done. We show that comparison in Results, below.

Error estimation
To estimate the error involved in Eq. 25, we compute the average resting fluorescence over all the pixels of the image, N_pixels, before Ca²⁺ release, and its standard deviation, We assume that F_rest=F_b. We then assume that the standard deviation of the basal fluorescence level, at any given point (before Ca²⁺ release), is given by Therefore, we assume that the error of approximating by Eq. 25 is Furthermore, the measurement of the fluorescence itself always carries an error, F. Taking these two sources of error into account and using Eqs. 24, 6, and F_b = F_rest, these errors translate into an error in [CaB] given by

(26)

For the experimental signals, we assume that the photon counting process is a Poisson process, so that the square of the standard deviation and the mean number of photons counted for an event of a given size are the same. Thus, we assume that Equation 26 with is the starting point for the error calculations in the case of the experimental signals. In the case of the noisy numerical simulations in which the normally distributed noise is added directly to the (blurred) [CaB] distribution, we use Eq. 26 with and _F = 0. In the case of the noiseless numerical simulations we assume that [CaB] has no errors.

Once [CaB] is known, the next step of the algorithm requires the calculation of time and space derivatives. To calculate the first derivative of a function, , at time t_i, we take a total of three data points, (t_i–1), (t_i), and (t_i+1) and compute the slope of the straight line that gives the best fit to those three data points using least squares. In this way we also obtain an error associated to the derivative. Space and higher derivatives are computed similarly. These errors are then propagated so as to obtain the errors of all the quantities of interest. The fitting procedure in Step 5 of the algorithm introduces new errors. As explained in Appendix 2, this fitting is done by minimizing, for every binned Ca²⁺ concentration value, C_i—the sum of the square of the difference between the left- and right-hand sides of Eq. 11. As a result of this process, the errors in g_i, f_i, h_i, and k_i are determined in a standard way (see Appendix 2). These errors are then propagated so as to obtain an error in M which then propagates into an error for Q_Ca and, consequently, for the current that the algorithm gives.

Kinetic data obtained from the experimental records
The experimental records provide a visualization of Ca²⁺ release upon individual channel openings. Thus, given the time series of the reconstructed current, I, it is possible to obtain some kinetic information, such as the mean rising, open, and closing times of the channel. In this article we focus on the open time, t_o, which, for each event, we compute as t_o = t₂ – t₁ where I(t₂) = I(t₁) = max I/2, with max I as the maximum value attained by the current during the event of interest. For comparison, we also compute it using the time series of the fluorescence averaged over a 0.45-µm region around the center of the release, as done in Demuro and Parker (2003).

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS APPENDIX 1: MOTIVATION FOR... APPENDIX 2: DESCRIPTION OF... ACKNOWLEDGEMENTS REFERENCES

 
Numerical tests of the reconstruction algorithm
We show here the results of applying the algorithm to numerically simulated images for which we know all the details, so that we can compare the estimates that the algorithm gives with the actual values used in the simulations. In particular, we investigate how well the algorithm can determine the source size and geometry, the time course of Ca²⁺ release, and the total current. The numerical simulations were performed with the resolutions listed in Table 1. However, only points separated by 10^–4 s in time were kept for the images and the reconstruction, unless otherwise noted.

Steplike single source
The first set of tests corresponds to simulations of Eqs. 15–17 with Ca²⁺ release of the form in Eq. 19. One such example is the one displayed in Fig. 1, for which R_f = 0.15 µm, t_s = 3 ms, and t_d = 10 ms. We can observe in Fig. 1 G that the algorithm determines the spatial distribution of the Ca²⁺ source and its temporal time course quite well: the source fades away for radii greater than 0.15 µm and times larger than 13 ms. We can also observe that the fitted functions f, g, and h are relatively small for large enough [Ca²⁺], and that the individual values, f_i, h_i and g_i, obtained for each binned [Ca²⁺] value are widely dispersed around the mean that the global fitting provides. For this reason, we tested the performance of the algorithm when we set f, g, and h exactly equal to zero obtaining current estimates that were 20% better (see Appendix 2). All the figures that we subsequently show in this article have been obtained in this way.

We show in Fig. 2 a plot of the magnitude of the reconstructed current obtained via Eqs. 13 and 22 as a function of the one used in the simulations, I, for simulations done with R_f = 0.15 µm, t_s = 3 ms, and t_d = 10 ms. The slope of the linear fitting is 0.96 ± 0.04. This figure was done using simulations with I up to 3.9 pA.

View larger version (14K): [in this window] [in a new window]   FIGURE 2  Application of the algorithm to data generated by numerical integration of Eqs. 15–17 with Eq. 19 as a source. The currents obtained with the algorithm are shown as a function of the ones used in the simulations. The best linear fit to these data of slope 0.96 ± 0.04 (dotted line) and the identity function are also plotted (solid line).

View larger version (22K): [in this window] [in a new window]   FIGURE 3  Algorithm's ability to distinguish different source parameters. The simulated data correspond to Eq. 19 (A and B), or to Eq. 20, with Rf = 0.15 µm, td = 5 ms, = 2 ms, and k1 = 0 (C and D). (A) Currents as a function of time, obtained with the algorithm (thick dashed, dash-dotted, and solid lines) and used in the simulations (thin dashed lines). Simulations done with I = 1 pA, Rf = 0.15 µm, and ts = 1 ms (indicated by an up-arrow), and td = 5, 10, and 15 ms, respectively (indicated by down-arrows). (B) Ratio between the time average of the source obtained with the algorithm, and the source amplitude used in the simulation, as a function of space. Simulations done with td = 10 ms and I = 0.1 pA and Rf = 0.15 µm; I = 3.7 pA and Rf = 0.5 µm; and I = 29.6 pA and Rf = 1 µm, respectively (source radii indicated by thin dotted lines). (C) Current used in the simulation (solid line) and obtained with the algorithm (dotted line), as a function of time, for a case with exponential decay and I = 1 pA. (D) Ratio between the decay time constants obtained with the algorithm and those used in the simulations as a function of I.

Steplike multi-opening source
We study now whether the algorithm can distinguish successive openings and closings of a single spherical source. Fig. 4 A shows the Ca²⁺ bound dye distribution from a simulation of Eqs. 15– 17 with Ca²⁺ release of Eq. 19, with R_f = 0.15 µm, I = 1 pA but where the source opens and closes three times with pulse durations of 1, 2.5, and 1 ms, respectively. Between a closing and the successive opening no activity is present. In this case, only data points with r > 0.15 µm (at all times) were used to implement Steps 4 and 5 of the algorithm. Fig. 4 B shows the reconstructed source where it can be seen how well the algorithm distinguishes between active and inactive regions. The algorithm also recognizes the three successive openings as well as the pulse durations, as can be observed in Fig. 4 C, where the current obtained with the algorithm (dotted line) is compared with the one used in the simulation (solid line). As is apparent from the figure, the current magnitude is estimated with the same accuracy regardless of the duration of the release event.

View larger version (24K): [in this window] [in a new window]   FIGURE 4  Algorithm's ability to recognize successive openings and closings of the source. The simulated data correspond to Eq. 19, with Rf = 0.15 µm, I = 1 pA, and ts = 0, td = 1 ms; ts = 1.5 ms, td = 2.5 ms; ts = 5 ms; and td = 1 ms, respectively. (A) Distribution of calcium bound dye from the simulation. (B) Reconstructed source. (C) Currents used in the simulations (solid line) and reconstructed ones (dotted line), as a function of time.

View larger version (22K): [in this window] [in a new window]   FIGURE 5  Effects of defocusing, blurring, finite space and time resolution, and noise on the reconstructed current. The simulations correspond to Eq. 19, with Rf = 0.15 µm, td = 10 ms (except for G and H, where td = 100 ms), and different values of I. The optical blurring was simulated using Eq. 21 with different values of xy and z. Only the simulations in G and H have added noise. (A) Effects of defocusing without blurring. The offsets are do = 0 (circles), do = 0.0707 µm (squares), do = 0.1118 µm (diamonds), and do = 0.1414 µm (asterisks). The solid line is the expected value for the current. (B) Effects of blurring and defocusing combined. Results with no defocusing or blurring (circles), with realistic blurring only (squares), with blurring and defocusing of do = 0.1414 µm (triangles), and with only defocusing of do = 0.1414 µm (asterisks). (C and D) Effects of blurring and spatial resolution without defocusing. Results with zero blurring (circles), with xy = 0.3 µm, z = 0.7 µm (squares), and with xy = 0.5 µm, z = 0.9 µm (diamonds). In C all the data points of the simulation were used, whereas only points separated by 0.15 µm were used in D. The slopes of the best linear fits are 0.96 (circles), 0.85 (squares), 0.78 (diamonds) in C and 0.76 (circles), 0.56 (squares), 0.48 (diamonds) in D. (E) Time-average of the source used in a simulation as a function of the line-scan coordinate (thin solid line), and corresponding ones obtained by applying the algorithm to the simulated data after it had been blurred using Eq. 21 with d = 0.15 µm (thick solid line), and after it had been subsequently deblurred by application of the inverse of Eq. 21 (dashed line). (F) Similar to E, but for QCa(x = 0) as a function of time. (G) Effects of noise and finite temporal resolution on the magnitude of the reconstructed current. Comparison of results with realistic blurring, finite space and time resolutions without noise (circles) and with noise (diamonds). The solid line is the expected value for the current. (H) Effects of noise and finite temporal resolution on the time course of the reconstructed current. Reconstructed currents (solid lines) and profiles used in the simulations (dashed lines) for an example with 0.1 pA (upper panel) and one with 0.3 pA (lower panel).

We show in Fig. 5 E the time average of the source used in the simulation (dotted line), and the ones obtained by applying the algorithm to the simulated clean data with d = 0.15 µm (thin solid line), to the simulated data after it had been blurred using Eq. 21 with d = 0.15 µm (thick solid line), and to the same data but after it had been deblurred by application of the inverse of Eq. 21 (dashed line). We also applied the algorithm to previously deblurred data, keeping all the data points of the simulation, obtaining a source that is indistinguishable from the one obtained using clean data (data not shown). We then conclude that the difference between the source obtained from clean data (thin solid line) and the one obtained from blurred-deblurred data (dashed line) is a consequence of the 0.15 µm space resolution. Fig. 5 F is similar to E, except for Q_Ca being at the center of the release region (x = 0) as a function of time. The dotted line corresponds to the simulated source and the thin solid line to the one obtained by applying the algorithm to the clean simulated data. We observe in Fig. 5 F that the algorithm slightly overestimates the value of Q_Ca(x = 0), but that this is compensated by an underestimation in other spatial points. We can also observe that the spatial distribution of the source is the feature that differs the most between the reconstructed source and the one used in the simulation. The temporal profile, on the other hand, including both the times at which the source is turned on and off, are very well captured by the reconstruction algorithm in all cases.

Effects of noise and finite time resolution
Fig. 5, G–H, shows the effects of noise and finite time resolution when there is realistic blurring and no offset. The simulations are done as in Fig. 5, C–F, blurring the simulated data with _xy = 0.3 µm, _z = 0.7 µm, and keeping points every 0.15 µm in space and every 8 ms in time. The noisy data correspond to the addition of normally distributed noise to the blurred fluorescence, as explained in Methods, above. We show in Fig. 5 G the magnitude of the reconstructed currents as a function of the ones used in the simulations when the algorithm is applied to noiseless (circles) and to noisy (diamonds) data. The slopes are 1.05 ± 0.08 and 1.02 ± 0.18, respectively. The solid line is the identity function. As can be seen, there are no significant differences in the magnitude of the reconstructed currents for the noiseless and noisy cases. It must be noted that the current estimation is improved with respect to that of Fig. 5 D because the criterion given by Eq. 14 prescribes a larger value of r_source(t) for the data with 8-ms time-resolution. Fig. 5 H shows two examples in which the magnitude and time course of the currents obtained with the algorithm when applied to the noisy data (solid lines) can be compared with the ones used in the simulations (dashed lines). The release events can be distinguished in both cases, although the uniform noise level that we have added, which is relatively large in the regions with low [Ca²⁺], smears some of the release events in the case with 0.1 pA.

Effects of the presence of endogenous buffers with different kinetics
We study now the performance of the algorithm when applied to a model with several endogenous buffers of different kinetics. We use the model of Baylor et al. (2002) with the modifications described in Methods, above, and the parameters of Table 2, with [Ca²⁺]_basal = 0.05 µM always, regardless of the total amount of dye used. This is a model with very large amounts of ATP and Parvalbumin, which results in a fast capture of Ca²⁺. In all cases, the algorithm is able to determine the duration and spatial spread of the source as accurately as in the cases illustrated in Fig. 3 (data not shown). However, due to the effective Ca²⁺ buffering by ATP and Parvalbumin, the Ca²⁺ current can be more or less underestimated, depending on the dye concentration. We show in Fig. 6 a plot of the currents that the algorithm gives as a function of the ones used in the simulations for different values of the total dye concentration. We observe that, for a total dye concentration, [B]_T = 40 µM, the estimated current is 67% of the simulated one (solid squares) and for [B]_T = 1 mM (asterisks), the current is 74% of the simulated one. If [B]_T is not large enough, the dye cannot compete with all the buffers and it is unable to sense the fast [Ca²⁺] changes that follow release. Since the algorithm relies on information that can be extracted from [CaB], the released current is underestimated. However, if the [Ca²⁺] concentration dynamics can be described correctly by Eq. 4 with M given by the ansatz (Eq. 5), then the form of the functions f, g, h, and k should be independent of the properties of the dye (including its total concentration) and of the release, Q_Ca. However, when the dye cannot sense the fast variations of [Ca²⁺], the functions f, g, h, and k obtained from its corresponding [CaB] are not determined correctly and do change with these properties. This results in a poorer determination of the [Ca²⁺] current, for which a correction factor could be estimated by analyzing the change of the functions f, g, h, and k with properties from which they should be independent. If the functions are correctly determined, we expect the algorithm to give a relatively good estimate of the [Ca²⁺] current. To check this, we used the functions obtained for [B]_T = 1 mM to analyze the simulations done with other values of [B]_T. Namely, we applied Steps 1–5 of the algorithm to the set of numerically simulated images generated with [B]_T = 1 mM and determined the function k (we set f = g = h = 0). We then used this function to approximate M (Eq. 5) and, using this expression for M, we applied Step 6 of the algorithm to the numerically simulated images generated with [B]_T = 40 µM. We show in Fig. 6 the currents we obtain in this way as a function of the ones used in the ([B]_T = 40 µM) simulations (open squares), which can be fitted by a linear relation of slope 0.74. The improvement in the current value determination (from 67 to 74% of the actual one) means that the ansatz (Eq. 5) may capture the dynamics of the poorly known processes, even in the presence of various buffers of different kinetics, but that the functions that define it can only be determined with the necessary accuracy if the dye concentration is large enough to sense the fastest [Ca²⁺] changes.

View larger version (14K): [in this window]