## Abstract

The amplitude and time course of stimulus-evoked second messenger signals carried by intracellular changes in free calcium ([Ca]_{free}) depend on the total influx of Ca^{2+}, the fraction bound to endogenous buffer and the rate of extrusion. Estimates of the values of these three parameters in proximal dendrites of 15 mouse α retinal ganglion cells were made using the “added buffer” method and found to vary greatly from one experiment to the next. The variations in the measured parameters were strongly correlated across the sample of cells. This reduced the variability in the amplitude and time course of the dendritic Ca^{2+} signal and suggests that the expression of Ca^{2+} channels, binding proteins and extrusion mechanisms is homeostatically coordinated to maintain the amplitude and kinetics of the Ca^{2+} signal within a physiologically appropriate range.

## Introduction

Calcium is a ubiquitous intracellular second messenger that can regulate enzymes, operate switch proteins and modulate ion channels to govern a wide variety of essential cellular actions. As for most molecular signals, the functional consequences of messages carried by calcium (Ca^{2+}) depend on the amplitude and time course of changes in its intracellular concentration. In retinal ganglion cell (RGC) dendrites, light-evoked synaptic input causes a change in membrane potential that opens voltage-gated Ca^{2+} channels (VGCCs) leading to an influx of Ca^{2+} ions and a transient increase in its intracellular concentration (Denk and Detwiler, 1999; Oesch et al., 2005; Margolis et al., 2010). While Ca^{2+} signals have been recorded from RGC dendrites using fluorescent indicators, the endogenous dynamics of the change in Ca^{2+} are distorted by the presence of the indicator. By binding Ca^{2+} the indicator acts as a buffer that reduces the amplitude and slows the recovery kinetics of the cell's normal Ca^{2+} signal. To better understand Ca^{2+} signals and the mechanisms that regulate them, both in the presence and in the absence of indicator, the ”added buffer” method (Neher and Augustine, 1992; Zhou and Neher, 1993; Neher, 1995) was used to measure the three parameters that govern Ca^{2+} dynamics in small neuronal compartments, i.e., the endogenous Ca^{2+} binding ratio (κ_{e}), which is the fraction of added Ca^{2+} that is bound by buffer, the rate of Ca^{2+} extrusion (Γ) and the amplitude of the stimulus-evoked Ca^{2+} influx ([Ca]_{tot}). A similar method has previously been used to document age-dependent variation of these properties in the somas of a population of rat retinal ganglion cells (Mann et al., 2005).

As with many biological parameters (Edelman and Gally, 2001; Marder and Goaillard, 2006) there was significant heterogeneity in the estimated values of the three primary parameters κ_{e}, Γ, and [Ca]_{tot,} with ∼10- to 60-fold variation over the population. Their values for each cell were correlated, however, in that cells with a higher than average Ca^{2+} binding ratio were also likely to extrude Ca^{2+} at a higher rate and to sustain larger Ca^{2+} influxes. We found that the correlations have the effect of reducing cell-to-cell variability in the amplitude and time course of the Ca^{2+} signal that would otherwise be present in an uncorrelated population. This suggests that cells employ a homeostatic mechanism to balance the expression levels and localization of Ca^{2+} binding proteins, pumps and ion channels so that their function is robust to variation in protein expression and that [Ca]_{free} is maintained within a physiologically appropriate range.

## Materials and Methods

##### Animals and tissue preparation.

The Administrative Panel on Laboratory Animal Care at the University of Washington approved all experimental protocols. Light adapted adult (5–8 weeks) male wild-type C57BL/6 mice (Jackson Laboratory) were killed by cervical dislocation. Both eyes were enucleated and placed in room-temperature Ames' medium (Sigma) bubbled with 5% CO_{2}/95% O_{2}. The intact retina was isolated from the retinal pigment epithelium in room light and then adhered photoreceptor-side down onto translucent Anodisc filter paper (Whatman) by wicking away excess solution. The flat-mount retina was put into a recording chamber with a clear glass bottom, placed on the stage of a custom made upright two-photon laser scanning fluorescence microscope (Euler et al., 2009) and perfused with warmed (29−33°C) Ames' medium with tetrodotoxin (TTX) (100 nm, Sigma) at a flow rate of ∼6 ml/min.

##### Cell targeting and classification.

Retinal ganglion cells with large (∼20 μm diameter) somas, which we refer to as α cells (αRGCs) (Peichl, 1991) were targeted for whole-cell recording. Each cell was further classified as either OFF or ON subtype based on the depth of its dendritic stratification and on its responses to full-field flashes of 577 nm light focused on the photoreceptor layer via the imaging objective (Euler et al., 2009). Cells that did not have a strong light response or conclusive dendritic morphology were not classified (NC) and are labeled as such. The experiment and analysis were attempted on ∼90 αRGCs and successfully completed on a total of 15 cells (5 ON-type; 7 OFF-type; 3 NC).

##### Electrophysiology.

Retinal ganglion cells were visualized with a CCD camera (Watec) using infrared (860 nm) illumination and exposed by micro-dissecting the inner limiting membrane with an empty patch pipette. Patch-clamp recordings were made of single cells in whole-cell voltage-clamp mode using pipettes (3–6 MΩ) filled with an internal solution that contained (in mm): 120 K-gluconate, 5 NaCl, 5 KCl, 5 HEPES, 1 MgCl_{2}, 1 ATP, 0.1 GTP, 0.11 Fluo-5F, and 0.06 Alexa 594 hydrazide (Invitrogen) adjusted to pH 7.4 with KOH. Series resistance was monitored closely and data were discarded *post hoc* if resistance increased during a recording by >65% or if the estimated loading time constant was >10 min. Experiments were rejected when there were obvious changes to the integrity or quality of the electrical and/or optical recording caused by either large variation in series-resistance, increased leakage current, hindered indicator loading or optical changes resulting from tissue movement or photodynamic injury. Voltages were corrected for the −11 mV liquid junction potential. Signals were amplified with a Multiclamp 700A (Molecular Devices) and digitized with an ITC-18 (InstrucTech). Data were analyzed offline using custom Matlab (MathWorks) scripts in conjunction with WinBUGS (Thomas et al., 1992), a free software package for implementing hierarchical, random-effects Bayesian inference models. Results are reported as the mean ± the SEM unless stated otherwise.

##### Imaging.

RGCs were filled with Fluo-5F and Alexa 594 and excited by two-photon absorption of ∼ 200 fs pulses of 910 nm light from a 532 nm pumped titanium/sapphire laser (Verdi-V10 and Mira-900; Coherent Inc.). Green and red epifluorescence signals were collected by a 20× water immersion objective (0.95 numerical aperture XLUMPLANFI, Olympus) and measured simultaneously by two independent photomultiplier detectors (Hamamatsu H7422P-40) after passing through either a green (HQ 510/50m-2p-18deg) or red (HQ 622/36m-2p-18deg) dichroic filter (Chroma Technology). The excitation laser (average power at back aperture of the objective ∼5 mW) was scanned in a line across the width of a proximal dendrite (<10 μm from the soma) at 2 ms per line. Line scan data were digitized and stored using custom hardware and CfNT software originally written at Bell Labs (Murray Hill, NJ) by R. Stepnoski and modified extensively by M. Müller at the Max Planck Institute for Medical Research in Heidelberg, Germany (Euler et al., 2009).

##### Fluorescence calibration.

The ratio of the intensities of the green to red epifluorescence signals is used to calculate the concentration of free calcium ions ([Ca^{2+}]) in the cytosol using a ratiometric calibration method detailed in previous studies (Maravall et al., 2000; Yasuda et al., 2004). Briefly, the dendrite is considered to be a single compartment containing Ca^{2+}, endogenous (*B*_{e}) and exogenous Ca^{2+} buffers, which in the latter case is the indicator Fluo-5F (*B*_{ind}), and a mechanism that extrudes Ca^{2+} with a rate constant of Γ (s^{−1}). The binding kinetics of both buffers are assumed to be fast enough that the compartment is always in a state of dynamic equilibrium represented by the chemical equation:
where Ca*B*_{e} and Ca*B*_{ind} represent Ca^{2+} ion bound to endogenous buffer and indicator, respectively. This is a reasonable assumption based on reports of Ca^{2+} buffers with *k*_{on} = 6 × 10^{8} m^{−1} s^{−1} that equilibrate with a physiological change in Ca^{2+} in <50 μs (Klingauf and Neher, 1997; Matveev et al., 2004; Cornelisse et al., 2007). At equilibrium the bound fraction of indicator can be expressed as a function of free Ca^{2+} ([Ca_{free}]) and the affinity of Fluo-5F in cytosol (*k*_{d}^{ind}):
where [*B*_{ind}]_{tot} is the total concentration of bound and unbound indicator present and *k*_{d}^{ind} = 1.3 μm for Fluo-5F in cytosol (Woodruff et al., 2002; Yasuda et al., 2004). The bound fraction of buffer is also related to the strength of the background-subtracted green fluorescence signal (*F*_{green}), which allows fluorescence to be related to Ca^{2+} using the equation:
where *F*_{green}^{max} and *F*_{green}^{min} are the intensities of green fluorescence when all and none of the available indicator are bound to Ca^{2+}, respectively. When the dynamic range of the indicator, defined as the ratio *F*_{green}^{min}/*F*_{green}^{max}, is large, as is the case for Fluo-5F with reported values between 40 and 240, Equation 3 can be simplified (Yasuda et al., 2004) as:
To correct for the effects of variations in the concentration of Fluo-5F and possible variations in optical parameters that could affect the fluorescence measurement (e.g., path length and laser pulse intensity) the ratio of green fluorescence to the red, Ca^{2+}-independent fluorescence is a more robust indicator of free Ca^{2+}. Substituting this ratio into Equation 4 we have:
where (*F*_{green}/*F*_{red})_{max} is the ratio of green to red fluorescence when all of the available Fluo-5F is bound to Ca^{2+}. This saturating value was measured as the peak signal during a 2 s, 100 mV depolarizing step from −91 to + 9 mV during which [Ca]_{free} rose and then plateaued.

##### Analysis of Ca^{2+} transients.

In our experiments a transient increase in [Ca]_{free} is elicited by a 20 ms somatic voltage step from −91 to +9 mV. Since Ca^{2+} influx is much faster than extrusion it is modeled as a delta function whose integral is equal to the total increase in [Ca]_{free} before reaching equilibrium with intracellular buffers. The resulting peak amplitude of the Ca^{2+} transient depends on the fraction of Ca^{2+} that is subsequently bound by endogenous buffer and indicator, which is known as the Ca^{2+} binding ratio of the buffer (κ_{e}) and the indicator (κ_{ind}), respectively. For a specified influx of Ca^{2+} ([Ca]_{tot}) the resulting amplitude (Δ[Ca]) is given by the equation:
where κ_{e} is defined as the endogenous differential binding ratio ∂[Ca*B*_{e}]/∂[Ca] and κ_{ind} is the differential binding ratio of the indicator, ∂[Ca*B*_{ind}]/∂[Ca] (Neher, 1995). For a known concentration of Fluo-5F, its binding ratio can be calculated using the equation:
where [*B*_{ind}]_{tot} is the total concentration of indicator. For simplicity we assume that κ_{ind} does not change during the transient change in [Ca]_{free} and instead use an estimate called the incremental Ca^{2+} binding ratio (Neher and Augustine, 1992; Helmchen et al., 1996) given by the equation:
where [Ca]_{rest} = 100 nm_{,} the concentration of Ca^{2+} before stimulation and Δ[Ca]_{est} is an estimate of the peak change in Ca^{2+} during stimulation. For a single cell, Δ[Ca]_{est} is a constant equal to the mean amplitude of all the transients.

Following a brief influx, Ca^{2+} is extruded from the compartment resulting in a slow, exponential decay in [Ca]_{free} approaching [Ca]_{rest}. The combined action of various pumps and exchangers is assumed to be linear such that the rate of Ca^{2+} efflux is given by the equation:
where Γ is the extrusion rate constant. Integration of Equation 9 combined with an equation for Ca^{2+} influx yields the equation for a Ca^{2+} transient:
where the time constant of decay (τ) is proportional to the total Ca^{2+} binding ratio of the endogenous buffers and the indicator:

##### Experimental determination of κ_{e}, Γ, and [Ca]_{tot}.

The “added buffer” method for determining the endogenous Ca^{2+} binding ratio (κ_{e}), the extrusion rate constant (Γ) and the total Ca^{2+} influx elicited by the stimulus ([Ca]_{tot}) relies on systematic variation of the concentration of indicator used to measure evoked calcium transients. This is achieved by eliciting Ca^{2+} transients as the solution in the patch pipette diffuses into the cell. Following the rupture of the membrane patch the two fluorophores in the intracellular dialysis solution enter the cell and their intracellular concentrations increase as the cell equilibrates with the concentrations in the pipette. As expected from previous studies of the intracellular dialysis of small molecules from a whole-cell patch pipette (Pusch and Neher, 1988; Mathias et al., 1990) the increase in *F*_{red} follows a monoexponential time course with a time constant (τ_{load}) given by:
where *t _{i}* is the time of the

*i*th transient after rupturing the patch and

*F*

_{red}

^{max}is the maximum red fluorescence signal when the concentration of Alexa 594 in the cell is equilibrated with the concentration in the whole-cell dialysis solution. The observed loading time constants were proportional to the access resistance and ranged from 2 to 7 min. Since Fluo-5F and Alexa 594 have similar molecular weights (931 MW and 759 MW, respectively) the increase in Fluo-5F after rupture was expected to follow the same time course and τ

_{load}was used to estimate the time-dependent increase in the intracellular concentration of the Ca

^{2+}indicator during the loading period and the associated increase in the Ca

^{2+}binding ratio. From the relationships in Equation 6 and 11 we expect that the change in Ca

^{2+}binding ratio leads to a proportional increase in both the inverse amplitude of the Ca

^{2+}transient (Δ[Ca]

^{−1}) and the decay time constant (τ). A plot of τ versus the added Ca

^{2+}binding ratio of the indicator during each transient (κ

_{ind}) illustrates this linear relationship (see Fig. 2

*B*). The slope and negative

*x*-intercept of a fitted line provide an estimate of 1/Γ and κ

_{e}, respectively. Similarly, a plot of Δ[Ca]

^{−1}versus κ

_{ind}can be fitted with a line whose slope and negative

*x*-intercept yield 1/[Ca]

_{tot}and κ

_{e}respectively (see Fig. 2

*C*). Analyzing the Ca

^{2+}transients in this way yields two estimates of κ

_{e}. Differences in these estimates are not easily reconciled since there is no straightforward way to assess the error in the fitted parameter values. This is due to the difficulty of propagating error estimates through multiple fitting steps that combine multiple datasets.

##### Random-effects model.

To overcome the limitations of the least-squares (LSQ) approach we used a hierarchical random-effects model of the loading experiment. In such a model, the parameters are thought of as random variables whose probability distributions are related to the data through a specified set of equations. In this case the parameters include κ_{e}, Γ, [Ca]_{tot} and τ_{load} and others, which are related to the fluorescence measurements through Equations 5–12. Each source of variability in the data –biological or otherwise– is accounted for in the model. Since fluorescence measurements are made optically there is instrument noise in the measurements of *F*_{red} during the *i*th transient (*F*_{red,i}) that is modeled as a normal distribution with unknown variance σ_{red}^{2} and time-varying mean μ_{red,i}:
The mean red fluorescence of the *i*th transient (μ_{red,i}) is a function of three additional unknown parameters *F*_{red}^{init}, *F*_{red}^{max}, and τ_{load} that is derived from Equation 12:
where the elapsed time since rupturing the patch is given by *t _{i}*.

The free Ca^{2+} concentration at the *j*th time-point of the *i*th transient ([Ca]* _{ij}*) is also modeled as a normal distribution whose variance (σ

_{Ca,i}

^{2}) is unknown: and whose time-varying mean (μ

_{Ca,ij}) is derived from the equation for a Ca

^{2+}transient (Eq. 10): In this equation [Ca]

_{rest}is a previously defined constant,

*t*is the time at which each measurement is made and Δ[Ca]

_{ij}*and τ*

_{i}*are intermediaries that represent the amplitude and the decay time constant of the*

_{i}*i*th transient, respectively. It is at this stage that biological sources of noise in the experiment need to be accounted for in the model. Due to changes in the dendrite that may occur over the course of an experiment, the stimulus may not elicit a Ca

^{2+}influx with precisely consistent amplitude each time. This source of noise is incorporated by modeling the amplitude of each transient as a normal distribution whose variance (σ

_{Δ[Ca]}

^{2}) is an unknown parameter: The mean of the distribution at the

*i*th transient (μ

_{[Ca],i}) is an intermediary given by Equation 6: where Δ[Ca]

_{tot}and κ

_{e}are unknown parameters and κ

_{ind,i}is an intermediary representing the Ca

^{2+}binding ratio of the indicator during the

*i*th transient. Its value during each transient is given by an equation that comes from combining Equations 8 and 12: which is dependent on the unknown parameter τ

_{load}and the previously defined constants

*k*

_{d}

^{ind}, [Ca]

_{rest}, Δ[Ca]

_{est}and

*t*. The concentration of indicator in the pipette ([

_{i}*B*

_{ind}]

_{pipette}) is also a constant.

Additionally, changes in the dendrite that affect Ca^{2+} extrusion may cause the decay time constant of each transient to vary. This source of noise is incorporated by modeling the decay time constant of each transient (τ* _{i}*) as a normal distribution whose variance (σ

_{τ}

^{2}) is an unknown parameter: The mean of the distribution (μ

_{τ,i}) is an intermediary given by Equation 11: where Γ and κ

_{e}are unknown parameters and κ

_{ind,i}was previously defined in Equation 19.

Though in our presentation of the model we have broken it into many separate pieces, it is important to realize that all of the pieces (Eqs. 13–21) are interdependent. For example, an estimate of τ_{load} will depend not only on the red fluorescence data, but also on the Ca^{2+} data. This statement is true of all the unknown parameters of which there are 9 + *k*, where *k* is the number of transients recorded. Notably, in this random-effects model the endogenous Ca^{2+} binding ratio of the cell (κ_{e}) is a single parameter. By effectively combining the steps of the sequential LSQ method we have eliminated the two independent estimates of κ_{e} from the analysis.

##### Bayesian estimates of the unknown parameters.

The unknown parameters of a random-effects model can be estimated from data using a variety of methods depending on the complexity of the model. We chose to implement the model within a Bayesian framework because of the ease with which the results can be interpreted and because there is freely available software that aids in parameter estimation using empirical Bayes methods (WinBUGS) (Thomas et al., 1992). Each parameter was assigned a minimally informative prior distribution, which was relatively uniform over a broad range of values (see Fig. 3). The prior distribution is representative of the best estimate of each parameter before the experiment. Bayes' theorem is applied iteratively with each piece of data to yield a posterior distribution for each parameter—a conditional likelihood distribution given the observed data.

There is an additional step that is required to visualize the posterior distributions and to make biological inferences. Typically a probability distribution is visualized by evaluating the probability density over a range of parameter values using a fine mesh. In the case of a normal distribution this would yield a “bell curve” shaped plot that could be used to identify the median the mean and other likely parameter values. However, since a posterior distribution has one dimension for every parameter it becomes intractable to fully evaluate a probability distribution that has more than a few parameters. Instead, we used Markov Chain, Monte Carlo methods for drawing unbiased samples from the multidimensional posterior probability distributions. As the number of samples grows it converges on the actual posterior distribution. From these distributions inferences were made about the most likely value of each parameter and the reliability of the estimate given the data. Though the technical details of the sampling methods we used are omitted, they are handled mostly by the WinBUGS software package and are described fully previously (Gelfand and Smith, 1990; Gilks et al., 1996).

## Results

In intact whole-mount mouse retinas, α retinal ganglion cells (αRGCs) were selectively targeted for whole-cell voltage-clamp recording and identified as either ON or OFF-type (see Materials and Methods). The results and analysis of the standard experiment are presented first using a single αRGC as an example, followed by a summary of the results of the identical experiment performed on 14 additional cells.

### Amplitude and recovery time course of a Ca^{2+} transient changes during loading period

The goal was to monitor the changes in the amplitude and recovery kinetics of stimulus-evoked increases in [Ca]_{free} during the “loading period” as the cell filled with Ca^{2+} indicator. Fluorescent Ca^{2+} signals were optically recorded from a proximal (<10 μm from the soma) dendrite (Fig. 1*B*) during a brief influx of Ca^{2+} ([Ca]_{tot}) that was evoked reproducibly by a brief (20 ms) depolarizing step of the membrane potential from −91 to +9 mV. The first stimulus was delivered shortly after patch rupture and at seven additional times during the loading period at intervals ranging from 20 to 120 s (Fig. 1). All experiments were conducted in the presence of TTX (100 nm) to eliminate voltage-gated Na^{+} current and thereby reduce potential variations in the stimulus introduced by fast, unclamped voltage changes. The depolarizing pulse elicited a transient inward current followed by a net outward current that lasted for the duration of the step. The early inward current was not an accurate measure of Ca^{2+} current (i.e., Ca^{2+} influx) because depolarizing stimuli were delivered abruptly after rupturing the patch, which did not leave sufficient time for proper adjustment of capacitance and series resistance compensation. In any given cell, however, the total stimulus-evoked Ca^{2+} influx ([Ca]_{tot}) was unknown but, constant over the loading period, since the depolarizing pulse was large enough to far exceed the threshold of maximum conductance for voltage-gated Ca^{2+} entry (Margolis et al., 2010).

During the loading period, as the concentration of indicator in the dendrite increased, the peak amplitude of the evoked Ca^{2+} transients decreased and their decay time constant grew longer (Fig. 1). In addition, the Ca^{2+} independent fluorescence of Alexa 594 (*F*_{red}) increased due to its increasing concentration in the dendrite as both fluorophores approached their concentrations in the pipette. These observations are consistent with a Ca^{2+} binding ratio of the dendrite that increases as it fills with exogenous Ca^{2+} indicator.

### Analysis of the Ca^{2+} transients reveals endogenous properties of Ca^{2+} dynamics

The changes in the amplitude and kinetics of evoked Ca^{2+} transients were used to estimate three primary parameters of Ca^{2+} dynamics, (1) the endogenous Ca^{2+} binding ratio (κ_{e}), (2) the extrusion rate constant (Γ) and (3) the total Ca^{2+} influx elicited by the stimulus ([Ca]_{tot}). From these parameters we could also estimate the amplitude (Δ[Ca]_{0}) and decay time constant (τ_{0}) of a stimulus-evoked Ca^{2+} transient in the absence of indicator. The set of Ca^{2+} transients and the Ca^{2+} independent fluorescence (*F*_{red}) were analyzed in two ways. The first method, which is referred to as the sequential LSQ method, is the most common and perhaps the most intuitive approach. It begins by converting the change in fluorescence triggered by each voltage pulse into a change in [Ca]_{free} for each stimulus using Equation 5. The falling phase of the response to each stimulus was fitted with a monoexponential to obtain an estimate of the amplitude (Δ[Ca]) and decay time constant (τ) of each transient. The additional Ca^{2+} binding ratio provided by the indicator during each transient was inferred from *F*_{red} by first fitting a plot of *F*_{red} versus time (Fig. 2*A*) with an exponential “loading” curve (Eq. 12). In the exemplar cell the derived loading time constant (τ_{load}) was 2.7 min and this was used to calculate the concentration of indicator ([*B*_{ind}]) at the time of each transient using the equation:
where *t _{i}* is the time of the

*i*th Ca

^{2+}transient and [

*B*

_{ind}]

_{pipette}is the concentration of indicator in the pipette. From [

*B*

_{ind}]

_{i}the Ca

^{2+}binding ratio of the indicator (κ

_{ind,i}) was calculated (Eq. 8). For a single compartment model (see Materials and Methods) the decay time constant (τ) of the Ca

^{2+}transient is expected to increase linearly with κ

_{ind}with a reciprocal slope given by the Ca

^{2+}extrusion rate constant (Γ) and a negative

*x*-intercept that corresponds to the endogenous Ca

^{2+}binding ratio (κ

_{e}) (Eq. 11). The exemplar cell considered here had an extrusion rate constant of 0.19 ms

^{−1}and an endogenous Ca

^{2+}binding ratio (κ

_{e}) of 10.4 (Fig. 2

*B*). The

*y*-intercept (κ

_{ind}= 0) of the extrapolated regression line shows that the decay time constant of a Ca

^{2+}transient in the absence of added indicator (τ

_{0}) would be 69 ms. Analogously, the inverse slope and the negative

*x*-intercept of the regression line fitted to a plot of the inverse amplitude (Δ[Ca]

_{i}

^{−1}) of each transient versus the Ca

^{2+}binding ratio of the indicator (κ

_{ind,i}) (Eq. 6) gave estimates of the total stimulus-evoked Ca

^{2+}influx ([Ca]

_{tot}) and the endogenous Ca

^{2+}binding ratio (κ

_{e}), which in this case were 16.2 μm and 12.0, respectively (Fig. 2

*C*). From the

*y*-intercept of the extrapolated regression line, the amplitude of a Ca

^{2+}transient in the absence of added buffer (Δ[Ca]

_{0}) would be 1.05 μm.

While the sequential LSQ analysis meets the objectives of the experiment by yielding estimates for each of the parameters of interest (i.e., κ_{e}, Γ and [Ca]_{tot}) it has a number of shortcomings. One issue arises from the fact that the eight Ca^{2+} transients are fitted independently of each other; the amplitude of each transient is determined using data exclusively from that transient. The model we have adopted stipulates, however, that the transients are inextricably linked to each other by the parameters κ_{e}, Γ and [Ca]_{tot}, which remain constant during loading and the intracellular concentration of the Ca^{2+} indicator whose concentration increases as the cell equilibrates with the contents of the pipette solution. Fitting each transient independently ignores these links and adds many free parameters to the fitting process. One symptom of this problem is that the analysis yields two estimates of the endogenous Ca^{2+} binding ratio of the cell, κ_{e}. One is based on the amplitudes of the transients and the other is based on their decay time constants, despite that in the single compartment model the Ca^{2+} binding ratio is specified as a single parameter that relates the amplitude and decay of a transient to one another.

Another limitation of the LSQ analysis is that there is no straightforward way to assess the error in the estimates of the parameter values. When raw data are fitted using linear regression, the correlation coefficient can be used to provide confidence intervals. However, the linear regression that yields the estimates κ_{e}, Γ and [Ca]_{tot} uses previously analyzed data (i.e., fitted parameters of Ca^{2+} transients), so the correlation coefficient does not accurately reflect the reliability of the fit. Since imaging measurements and biological systems are inherently noisy, it is important to appropriately propagate estimates of error through the analysis so that accurate measures of certainty can be placed on each of the parameters of interest.

To avoid the flaws of the sequential LSQ method, a hierarchical random-effects model of the experiment was created (see Materials and Methods for details). The result of the analysis is a probability distribution of possible values of each parameter, given the experimental observations. From this distribution, called the Bayesian marginal posterior distribution (Fig. 3), we can infer the most likely value of each parameter from the location of the peak and also a range, which has a 95% chance of encompassing the true value, termed the 95% Bayesian credible interval (CI). For the example cell, the 95% CIs of κ_{e}, [Ca]_{tot} and Γ are 6.1–18.9, 13.3–18.3 μm, and 0.17–0.21 ms^{−1}, respectively. For every parameter these ranges are inclusive of the estimates obtained using the LSQ method (Fig. 3), which is validation of the Bayesian model we used.

### Ca^{2+} dynamics are heterogeneous in the population of RGCs

The results from the Bayesian analysis of the 15 completed experiments are presented in Figure 4. The most likely estimate of the endogenous Ca^{2+} binding ratio in each cell is remarkably variable, ranging from 2 to 123. It is important to note that because we have determined credible intervals for each estimate, we believe that this variability represents true heterogeneity in the population as opposed to noise in our measurements. For example, in one cell the data indicate that there is a 95% chance that the true value of κ_{e} is between 42 and 160 whereas in another cell there is a 95% chance that the κ_{e} is between 1.5 and 6.8. There is also notable heterogeneity in the extrusion rate constant (Γ), which ranges from 0.07 to 0.9 ms^{−1} and in the amplitude of the stimulus-evoked Ca^{2+} influx ([Ca]_{tot}), which ranges from 10 to 82 μm (Fig. 4*B*,*C*). The heterogeneity that we have observed may represent cell-to-cell differences or location-to-location differences along the dendrite, since the recording is made from only a narrow section of dendrite (see Discussion for further consideration).

To determine whether some of this heterogeneity could be explained by differences in αRGC cell types, a statistical hypothesis test was performed on the parameter estimates for ON and OFF αRGC subtypes. Since the most likely estimates of each parameter were not normally distributed across the sample population, a permutation test was used instead of a Student's *t* test to determine whether the mean of the ON cells was significantly different from that of the OFF cells. For all parameters except [Ca]_{tot} the ON cells were statistically indistinguishable from the OFF cells. However, in response to the same voltage stimulus, OFF αRGCs sustain a Ca^{2+} influx that is 24.2 μm greater on average than ON αRGCs (*p* < 0.06). This is consistent with previous findings that OFF-type αRGCs possess low voltage-activated T-type Ca^{2+} channels in addition to high-voltage activated L-type Ca^{2+} channels, which are present in both ON- and OFF-type cells (Margolis and Detwiler, 2007; Margolis et al., 2010). These channels mediate stimulus evoked Ca^{2+} influx, which has not been shown to be augmented by Ca^{2+}-induced Ca^{2+} release from internal stores.

### Correlated variability of κ_{e}, Γ and [Ca]_{tot} reduces variability in the Ca^{2+} transients

While the estimated values of the primary parameters in the absence of exogenous buffer varied greatly across our sample of recorded αRGCs, there was strong positive correlation between κ_{e} and both Γ and [Ca]_{tot} (Fig. 5*A*,*B*). Cells with higher endogenous Ca^{2+} binding ratios also had higher extrusion rate constants and larger Ca^{2+} influxes. Since the distributions of the parameter values were not Gaussian, the correlations were quantified using Spearman's rank correlation test; a nonparametric statistic that provides a measure of the monotonicity of the pairing of two parameters that does not depend on the underlying distribution of either parameter. A rank correlation of +1 indicates that the relationship between the variables is perfectly described by a positive monotonic function. Values between 0 and 1 are comparable to those based on a *R*^{2} statistic if the distributions were in fact Gaussian.

Across the population of αRGCs the rank correlations between κ_{e} and Γ and between κ_{e} and [Ca]_{tot} were 0.96 and 0.57, respectively (Fig. 5*A*,*B*). Both correlations were statistically significant with *p*-values < 0.03, based on a two-tailed, nonparametric hypothesis test (Press et al., 2007). The rank correlations between κ_{e}, and Γ and [Ca]_{tot} based on estimates obtained from the LSQ approach were also statistically significant (*p* < 0.05). This shows that the covariation of the primary parameters is not an artifact of the Bayesian analysis and supports the conclusion that the correlations between the primary parameters are real and biologically relevant.

In the single-compartment model of Ca^{2+} dynamics κ_{e}, Γ and [Ca]_{tot}, determine the amplitude (Δ[Ca]_{0}) and decay time constant (τ_{0}) of the Ca^{2+} signal (Eqs. 6, 11). To evaluate the effect of the correlations between the primary parameters on the distribution of these secondary parameters, we calculated for each cell τ_{0} and Δ[Ca]_{0} using either the observed estimates of κ_{e}, Γ and [Ca]_{tot} or a random pairing of κ_{e}, Γ and [Ca]_{tot} from different cells. For a sample of 15 cells there are 225 possible values of τ_{0} and of [Ca]_{tot} from which 15 values can be drawn randomly for each, to give a simulated sample of secondary parameter values to be compared with the observed sample. This comparison was performed 10^{6} times and collectively they are the basis for what is known as a permutation test (Efron and Tibshirani, 1993). Note that the number of possible combinations of drawing 15 samples randomly from a distribution of 225 without replacement is 225!/15!(225–15)!; an astronomically large number. For each simulated sample of 15 cells the rank correlation of the primary parameters and the SD of τ_{0} and Δ[Ca]_{0} was calculated and compared with those of the observed distribution. Of the 10^{6} simulated samples none had a smaller SD of τ_{0} values than the observed sample (*p* < 10^{−6}) and <5% of the permutated distributions of Δ[Ca]_{0} had a smaller SD than 2.0 μm, the SD of the observed distribution (Fig. 5*C*,*D*). These permutation tests show that the rank correlation of the primary parameters significantly reduces the inter-cell variability of τ_{0} and Δ[Ca]_{0} across the observed sample.

The effect of the correlations between the primary parameters on the distributions of τ_{0} and Δ[Ca]_{0} were tested further by statistically generating a larger population of cells to examine how the correlations in the primary parameters influence the distributions of τ_{0} and Δ[Ca]_{0}. This was done by fitting the observed values of κ_{e}, Γ and [Ca]_{tot} with statistical distributions from which samples could be drawn to create simulated αRGCs with parameters that match the statistics of the observed sample. The observed sample of the primary parameters are not fit well by Gaussian probability density functions; the parameters are biologically constrained to be greater than zero and their mean values are much larger than their median values, indicative of long positive tails in their distributions. The parameter values were, however, well fitted with a lognormal distribution, a common distribution in biology (Koch, 1966). From these fitted distributions of κ_{e}, Γ and [Ca]_{tot}, 10^{6} fake cells were drawn at random and τ_{0} and Δ[Ca]_{0} were calculated for each cell (Fig. 5*E*,*F*). The resulting distributions of τ_{0} and Δ[Ca]_{0} had lower peaks and greater SDs (mode: 8.8 ms, SD = 830 ms and mode: 0.18 μm, SD = 5.6 μm) than the observed distributions of the secondary parameters (mode: 65 ms, SD = 32 ms and mode: 0.49 μm, SD = 2.0 μm). When rank correlations of 0.96 and 0.57, equivalent to the correlations present in the observed distribution of τ_{0} and Δ[Ca]_{0}, respectively, were induced in the simulated distributions of the primary parameters (Iman and Conover, 1982), τ_{0} and Δ[Ca]_{0} were recalculated and the mode and SD of the resulting distributions of τ_{0} (66 ms, SD = 46 ms) and Δ[Ca]_{0} (0.48 μm, SD = 2.17 μm) were similar to those of the observed distributions (Fig. 5*E*,*F*). This shows that inducing the observed rank correlation in a population of cells whose primary parameters are randomly distributed, is sufficient to reduce the SD of the distributions of τ_{0} and Δ[Ca]_{0} in an uncorrelated population to match those that we observed. More generally the analysis demonstrates how homeostatic regulation (i.e., correlation) of the parameters that govern Ca^{2+} dynamics could reduce variability in the Ca^{2+} signal across a heterogeneous population.

Overall the results of the analysis indicate that the observed correlation between the primary parameters reduces variability in the amplitude and time course of dendritic Ca^{2+} signals across the population of αRGCs.

### Numerical simulation of Ca^{2+} transients allows comparison of a global or a localized influx

The parameters that were experimentally derived from Ca^{2+} transients evoked by brief depolarization of the soma (global Ca^{2+} signal) were used in a NEURON-based electrotonic model (Hines and Carnevale, 2001) of an αRGC dendrite to investigate the properties of dendritic Ca^{2+} signals evoked by spatially localized stimuli. The model dendrite consists of a cylinder with a diameter of 1.5 μm and a total length of 5 μm that is divided into 101 cylindrical slices. This length of dendrite was sufficient for modeling the full spread of a Ca^{2+} signal. The simulation calculates the concentrations of free Ca^{2+}, unbound buffer and bound buffer in each slice, over time. Four processes affect these concentrations in the simulation: 1) instantaneous influx, 2) longitudinal diffusion, 3) binding/unbinding of buffer, 4) linear first-order extrusion. In the simulation the buffer is nondiffusible, consistent with previous studies showing that this is commonly the case for a majority of cellular Ca^{2+} binding sites (Allbritton et al., 1992; Zhou and Neher, 1993; Neher, 1995). The values of the three parameters of Ca^{2+} dynamics in the simulation (i.e., κ_{e}, Γ and [Ca]_{tot}) were set using the modes of log-normal distributions fitted to the population of experimentally determined values for each parameter. This resulted in Ca^{2+} signals with an amplitude (Δ[Ca]) of 0.48 μm and a decay time constant (τ) of 65 ms.

The simulation was validated using a global Ca^{2+} influx in which [Ca]_{free} was elevated instantaneously in every slice of the dendrite. The simulation was run several times, each time choosing a different value for κ_{e} from the entire range of values observed in the population (κ_{e} = 2–123) while altering Γ and [Ca]_{tot} proportionally so as to mimic homeostatic compensation. As expected the amplitude and decay time constant of the simulated transient did not change (Fig. 6*A*). The time course of transients evoked by a global Ca^{2+} influx is dominated by the rate of Ca^{2+} extrusion. In this scenario there is no longitudinal diffusion of Ca^{2+} because the global influx of Ca^{2+} is spatially uniform.

To address questions about the Ca^{2+} signal produced by a local influx, the model was used to predict the change in free Ca^{2+} that would result from an instantaneous increase in Ca^{2+} in a small compartment (0.55 μm wide) in the center of the 5 μm length of dendrite, similar to the width of a postsynaptic density that may experience a localized influx of Ca^{2+}. The endogenous Ca^{2+} binding ratio, the extrusion rate constant and the total Ca^{2+} influx were varied using the same values as above. Because a localized Ca^{2+} influx generates a longitudinal concentration gradient of [Ca]_{free}, Ca^{2+} begins to diffuse immediately following its transient influx (*D*_{Ca} = 0.1 μm^{2}/ms). This causes [Ca]_{free} in the small central compartment to decrease more rapidly than in the case of the global Ca^{2+} influx (Fig. 6*A*). Under these conditions the decay kinetics of the local change in Ca^{2+} is strongly dependent on the Ca^{2+} binding ratio even in the presence of homeostatic compensation of the extrusion rate constant. A higher Ca^{2+} binding ratio slows decay of [Ca]_{free} because the nondiffusible buffer binds Ca^{2+}, slows its diffusion and restricts its spatial spread (Fig. 6*B–E*). The precise decay time constant and the spatial extent of [Ca]_{free} in an actual RGC dendrite would also depend on the geometry of the dendrite, the width of the stimulated region and on the diffusion coefficient of free Ca^{2+}.

## Discussion

We set out to measure the values of the three primary parameters that control the amplitude and time course of stimulus-evoked changes in [Ca]_{free}. Since the targeted cells were all members of the α class of large soma RGCs with mono-stratified dendrites, the measurements were expected to provide a precise estimate of the value of each of the primary parameters for the population. The results showed, however, that the cell-to-cell variation in our measurements was much greater than could be accounted for by experimental error. The dynamics of Ca^{2+} signaling in αRGC dendrites could not be summarized by a single set of primary parameter values and could only be fully described by the distribution of observed values. This is demonstrated most plainly by the absolute range over which each parameter varied. In different cells the estimates of the endogenous Ca^{2+} binding ratio, extrusion rate and total Ca^{2+} entry varied between 2 and 123, 0.07 to 1 ms^{−1} and 10 to 80 μm, respectively; not inconsistent with the range of parameter values reported for other types of neurons (Helmchen et al., 1996; Maravall et al., 2000; Kaiser et al., 2001; Sabatini et al., 2002; Mann et al., 2005; Cornelisse et al., 2007).

Where does this variability come from? Some of it must arise from the stochastic nature of molecular biological and biochemical mechanisms, which have been shown to introduce large variations in the expression of mRNA and protein, including ion channels and neurotransmitter receptors, in the same, or even genetically identical, cell types (Edelman and Gally, 2001; Davis, 2006; Marder and Goaillard, 2006). In addition to the biological variations between cells, there are also regional differences in the protein expression levels of subcellular compartments in single cells. These within-cell differences are consistent with compartmentalized Ca^{2+} signaling in neurons, as shown, for example, by differences in the spatial distribution of voltage-gated Ca^{2+} channels and Ca^{2+} sensors in neurites and axon terminals (Magee, 2008). Since our measurements were made using laser line-scanning to probe the Ca^{2+} dynamics of a narrow (<1 μm) slice of dendrite, the results might be expected to vary depending on the local properties of the region of the dendrite that was scanned. While spatial variations in Ca^{2+} signaling have not been observed in the aspiny dendrites of αRGCs, they could arise from the presence of micro-domains of colocalized Ca^{2+} channels, binding-partners and extruders, as has been described for the aspiny dendrites of cortical interneurons (Goldberg et al., 2003).

### Correlated variability

The estimates of the three parameters that control the amplitude and time course of Ca^{2+} signals were strongly correlated such that variation in one parameter was compensated by variation in the others. As a result of this correlation, the Ca^{2+} signals were robust to the heterogeneity in the individual components that control Ca^{2+} dynamics. The positive correlations between κ_{e}, Γ and [Ca]_{tot} across our sample of αRGCs (Fig. 5) are consistent with the expression of the proteins involved in Ca^{2+} influx, buffering and extrusion being coordinated homeostatically. This would ensure that the primary parameters controlling Ca^{2+} dynamics are matched according to the size of the Ca^{2+} load they receive. A region of dendrite where Ca^{2+} influx is large also has higher endogenous buffer capacity and more clearance mechanisms to bring [Ca]_{free} back to its resting level. If the controlling parameters were not coordinated, the properties of the dendritic Ca^{2+} signal would vary along the dendrite in a way that could compromise its ability to serve as a useful signal. Consider for example if the Ca^{2+} binding ratio and extrusion rate constant were low in regions where Ca^{2+} entry was large, the amplitude of change in [Ca]_{free} would exceed the physiological range and last longer than would be compatible with integrating and processing visual signals.

The coordinated expression of the primary parameters would serve to make the dendritic Ca^{2+} signal, i.e., the stimulus-evoked change in free Ca^{2+}, more uniform over the dendrite by compensating for local differences in Ca^{2+} entry. The uniformity of the change in free Ca^{2+} does not preclude, however, the presence of spatial differences in activation of Ca^{2+}-dependent signaling cascades. The correlations between the parameters would tend to associate localized sites of high Ca^{2+} entry with more buffer, i.e., more Ca^{2+} binding proteins for triggering second messenger signaling.

Our results showing correlations between the parameters that control Ca^{2+} signaling are supported by a study that used the “added buffer” method to evaluate the properties of stimulus-evoked Ca^{2+} changes in the soma of newborn (PND 5–10) and adult RGCs (Mann et al., 2005). Their measurements indicate that the time course of the Ca^{2+} signal is essentially the same in newborn and adult RGCs (τ_{0} = 2.1 s vs 1.7 s) despite a large difference in their Ca^{2+} binding ratios (κ_{e} = 568 vs 156, respectively). While not commented on by the authors, the difference in endogenous buffering in the two age groups is compensated for by faster extrusion in newborn than adult RGCs (Γ = 277 vs 52 s^{−1}, respectively), which suggests that the amplitude and time course of the somatic Ca^{2+} signal varies little as the cell progresses to adulthood, despite large changes in buffer expression.

The existence of documented differences in the subcellular properties of different neuronal compartments (Lai and Jan, 2006) makes it seem likely—if not obvious—that coordinated regulation of the proteins controlling Ca^{2+} signaling is a common property of neurons. This has not been well recognized, however. The only studies of coordinated compensation, similar to what we have observed in αRGCs, report that the Ca^{2+} binding ratio is reduced in mutant cells that express VGCCs with lower open probabilities (Dove et al., 2000; Murchison et al., 2002). While the mechanism(s) responsible for the coordinated expression of the proteins that control Ca^{2+} dynamics or participate in other signaling systems is not known (MacLean et al., 2003; Davis, 2006), the presence of Ca^{2+} binding domains in the promoters for several Ca^{2+} binding proteins could provide a way for Ca^{2+} influx to regulate buffering (Arnold and Heintz, 1997).

Homeostatic compensation for variations in the primary parameters would give rise to Ca^{2+} signals that are more consistent across the sample of recorded αRGCs. The mean peak amplitude of the stimulus-evoked change in [Ca]_{free} was 1.6 ± 0.4 μm with a recovery time constant of 97 ± 9 ms (*n* = 15). The amplitude of this signal is more than twice the amplitude of the signal previously recorded from the soma of adult RGCs (Δ[Ca]_{0} = ∼ 600 nm) and it is ∼20 times faster (Mann et al., 2005). This may be the result of inherent differences in the dynamics of somatic and dendritic Ca^{2+} signals. There are, however, notable methodological differences in the present and earlier study that leave this as an open question until such time that Ca^{2+} signals in soma and dendrite are recorded using the same experimental protocols.

### Spatial spread

The change in dendritic Ca^{2+} produced by a local influx of Ca^{2+} is spatially confined by binding to endogenous buffer(s) that are either fixed or slowly mobile and as a result only spreads 1–2 μm from the site of the influx. If the endogenous buffer was substantially mobile with a diffusion coefficient *D*_{B}, a local change in Ca^{2+} would be reduced in amplitude as it spread out over a wider extent given by approximately (2*D*_{B} · τ_{0})^{1/2}, which would be ∼3 μm for *D*_{B} = *D*_{Ca}/2. Note that the intracellular incorporation of a diffusible Ca^{2+} indicator, such as the one used in this study, will act similarly to disperse a local change in Ca^{2+} making it more difficult to detect using fluorescence.

### Function

Retinal ganglion cells express several different members of the EF-hand family of Ca^{2+} binding proteins that include calretinin, calbindin D-28K and parvalbumin (Wässle et al., 1998; Ghosh et al., 2004), all of which are considered to be nondiffusible buffers that bind Ca^{2+} with high affinity ranging from 50 to 1500 nm (Falke et al., 1994; Schäfer and Heizmann, 1996; Schwaller et al., 2002). A number (5 to 8) of calmodulin-like orphan retinal Ca^{2+} binding proteins have also been identified in retina (Haeseleer et al., 2000). While the amplitude of the Ca^{2+} signal in αRGC dendrites is sufficient for it to be sensed by any of these of binding proteins the physiological consequences of such interactions are not known. With a decay time course of ∼90 ms it is possible that Ca^{2+} signals participate in the integration and processing of visually evoked synaptic inputs, especially in the mechanisms that underlie slower adaptational changes in αRGC function.

## Footnotes

Funding was provided by NIH Grants EY02048 (P.B.D.) and 5 T32 GM007108-35 (A.J.G.). A.J.G. submitted the data in this manuscript in partial fulfillment of the requirements for the degree of Doctor of Philosophy at the University of Washington in 2010. We are grateful to Bertil Hille and Fred Rieke for helpful comments on an earlier version of this manuscript, Winfried Denk for designing and aiding the assembly of our two-photon microscope, and Paul Newman for technical assistance.

- Correspondence should be addressed to either Andrew J. Gartland or Peter B. Detwiler, Department of Physiology & Biophysics and Program in Neurobiology & Behavior, Box 357290, University of Washington, Seattle, WA 98195. gartland{at}u.washington.edu or detwiler{at}u.washington.edu