Photon shot noise limits on optical detection of neuronal spikes and estimation of spike timing

Abstract

Optical approaches for tracking neural dynamics are of widespread interest, but a theoretical framework quantifying the physical limits of these techniques has been lacking. We formulate such a framework by using signal detection and estimation theory to obtain physical bounds on the detection of neural spikes and the estimation of their occurrence times as set by photon counting statistics (shot noise). These bounds are succinctly expressed via a discriminability index that depends on the kinetics of the optical indicator and the relative fluxes of signal and background photons. This approach facilitates quantitative evaluations of different indicators, detector technologies, and data analyses. Our treatment also provides optimal filtering techniques for optical detection of spikes. We compare various types of Ca(2+) indicators and show that background photons are a chief impediment to voltage sensing. Thus, voltage indicators that change color in response to membrane depolarization may offer a key advantage over those that change intensity. We also examine fluorescence resonance energy transfer indicators and identify the regimes in which the widely used ratiometric analysis of signals is substantially suboptimal. Overall, by showing how different optical factors interact to affect signal quality, our treatment offers a valuable guide to experimental design and provides measures of confidence to assess optically extracted traces of neural activity.

Figure 1

Signal detection theory. (A) Distribution of log-likelihood ratios, in cases with a spike (green curve), p_L(f|H⁽¹⁾), and without a spike (red curve), p_L(f|H⁽⁰⁾), for d′ = 1 and d′ = 3 and two different detection thresholds. (B and C) The receiver operating characteristic (ROC) curve describes the trade off between the true positive rate versus the false positive rate for different values of the detection threshold. (B) In the approximation these distributions are Gaussians of equal variance, d′ completely characterizes the ROC curve. For d′ ≥ 5, the ROC curve is visually indistinguishable from the panel axes. (C) The area under the ROC curve is a metric of spike detectability that is independent of the user’s choice of detection threshold. This area equals the probability of an ideal observer correctly identifying the spike in a choice between a pair of instantiations of each hypothesis. (D) Contour plot showing how d′ depends on the transient amplitude, ΔF/F, and $F_0\tau$, the mean number of background photons collected over an interval of τ in duration. Contours are plotted for d′ = 1–5, 10, 15, 20, and 25. Discriminability improves for increased F₀ at fixed ΔF/F, due to the rise in signal photons.

Figure 2

Simulations of spike detection. (A) Probability distributions of log-likelihood ratios, Eq. 6, taken over all possible photon measurements in the presence (green) and absence (red) of a neural spike, based on simulations of neural spiking with shot-noise-limited optical detection at d′ = 1 or d′ = 3. (Solid lines) Theoretically predicted Gaussian distributions. (Histograms) measured distributions from the simulation. The slight disagreements between the two are due to the use of leading-order approximations in the theory and can be remedied by including higher-order terms (see Section S1 in the Supporting Material). (B) Simulated optical traces and detected spikes for several values of d′. (Blue lines) Optical traces shown in units of the standard deviation from the mean photon count. (Green spikes) The true spike train. (Orange spikes) correctly estimated spikes. (Spikes in non orange hues) Spikes estimated with errors in frame timing. (Gray spikes) False positives. (Gray trace) For a moving window of nine time bins, the log-likelihood ratio for the hypotheses that there is (H⁽¹⁾) and is not (H⁽⁰⁾) a spike. (Dashed black line) Spike detection threshold, log C, given equal costs for false positives and false negatives (see Eq. 8). (Purple) Threshold crossings. Spikes were detected using an iterative, greedy algorithm that assigned a spike to the instance of the log-likelihood ratio’s maximum in each iteration. At low d′, few spikes are detected with this choice of threshold (see panel C). Simulation parameters: ΔF/F = 0.05, τ = 0.15 s, ν = 20 Hz, spike rate λ = 0.5 Hz. The cutoff is given by log C = log ((ν/λ) − 1). For d′ = 1, 3, 5, and 7, the probability of detection is 7.8 × 10⁻⁴, 0.61, 0.96, and 0.999, respectively, whereas the expected number of false positives is 0.0092, 1.9, 0.36, and 0.017, respectively. (C) Detection probability, P_D, and false positive probability, P_F, as a function of d′ calculated using Eqs. S17 and S18 in the Supporting Material and the same equal-cost condition as in panel B. For low values of d′, the false positive probability rises with d′, because at d′ = 0 there are very few spikes detected and thus very few false positives. At large d′ values nearly all spikes are detected, and the false positive rate nearly vanishes.

Figure 3

Spike detection with modest photon counts. When the number of photons collected from signal and background sources is insufficient to justify Gaussian approximations, the treatment involving d′ is not guaranteed to be valid. To explore this regime, we randomly sampled the log-likelihood ratio distributions of Eq. 6, computed the ROC curves, and calculated the underlying areas. The area under the ROC curve is plotted as a joint function of background and signal photons, (A), or of background photons and ΔF/F (B). (Dashed lines) Areas under ROC curves determined using the d′ calculation, demonstrating the robustness of the calculation result despite the invalidity of the Gaussian approximation. The contours in both plots are spaced in 0.05 increments of area under the ROC curve. (Note that unlike in Fig. 1D, the x-axes are on a linear scale.) As expected, the agreement between the direct calculations and the d′ approximation degrades as the condition A/F₀ ∼ ΔF/F ≪ 1 is relaxed. This leads to a maximum disagreement of ∼0.025 in the area under the ROC curve. (C) Simulated trace plotted using the same conventions as in Fig. 2. Simulation parameters: ΔF/F = 0.75, τ = 0.15 s, ν = 20 Hz, spike rate λ = 0.5 Hz. Probability of detection is 0.96 and the expected number of false positives is 0.36.

Figure 4

Spike detection using FRET indicators. In the shot noise limit, ratiometric treatment of FRET indicators generally makes suboptimal use of the photon statistics for spike detection (A and B). (A) The decline in spike detectability is plotted as a joint function of the ratios of signal amplitude A⁽¹⁾/A⁽²⁾ and the background brightness F₀⁽¹⁾/F₀⁽²⁾. When the fluorescence transients are equal in amplitudes, a direct statistical analysis does not yield any benefit. However, when the ratio of signal amplitudes differs from unity, a ratiometric analysis is inferior. (B) The decline in spike detectability is plotted as a joint function of the ratios of the background-normalized signal amplitudes (A⁽¹⁾/F₀⁽¹⁾)/(A⁽²⁾/F₀⁽²⁾) and the background brightness. (C) Simulated pairs of traces (blue and green traces) for the donor and acceptor channels of a FRET indicator. Here, we calculate the log-likelihood ratio of a spike with a moving window of 60 observations using Eq. 6 (the sum of the log-likelihood ratios for the donor and acceptor traces) and a Gaussian approximation (Eqs. S49–S51 in the Supporting Material) to the ratio of the donor and acceptor traces (purple trace). The true spike train is shown (green, top spike trains) above the spike trains estimated via direct analysis of the channels (middle spike trains) or a ratiometric analysis (bottom spike trains). (Spikes in non orange hues) Spikes estimated with errors in frame timing. (Gray spikes) False positive spikes. (Light- and dark-gray traces) Log-likelihood ratio of a spike for the ratiometric and direct channel analyses, respectively. Threshold crossing events for these two traces (highlighted in purple and blue, respectively) were analyzed using the greedy algorithm approach to detect spikes, applying the same equal-cost condition to set the detection threshold for both methods. Note the superior log-likelihood ratio values attained using the direct analysis. The strong correlations between the ratiometric and direct channel calculations occur because the same simulation trial was used for both calculations. In these simulations, ΔF/F = 0.05 for both channels, τ = 1.0 s, ν = 20 Hz, spike rate λ = 0.5 Hz.

Figure 5

Chapman-Robbins lower bounds on the estimation of spike times. The ability to accurately estimate a spike’s time is governed by the discriminability, d′, the sampling rate, ν, and the duration of the signal transient, τ. For fixed d′, spike timing precision improves with shorter optical transients and higher sampling rates. The Chapman-Robbins lower bound can be either lower or higher than the sampling rate and exhibits strong dependence on discriminability d′ and transient duration. (A) Comparison of three optical indicators with different decay kinetics. A fast optical indicator (τ ∼ 0.15 s) such as Oregon Green BAPTA-1 (OGB-1) can localize spikes to ∼30 ms, even in regimes with modest SNR regimes (d′ ∼ 3). Slower indicators with time constants resembling those of GCaMP3 (∼0.5 s) exhibit poorer spike timing characteristics. In regimes with modest SNR, even slower indicators with kinetics resembling those of d3cpv or TN-XXL (1.5 s) allow spike localization with ∼100 ms resolution. (B) Comparison of estimation bounds at three different sampling rates. Increasing the sampling rate modestly improves spike timing resolution. (C) Increasing discriminability also significantly improves spike timing resolution. (A–C, white dashed line) Boundary of the regime of temporal super-resolution. (D) Simulations of spike timing resolution. Using a brute-force maximum likelihood method for determining the spike time, we obtained histograms of the spike time error for two indicators with distinct temporal dynamics. Note the different time scales on the two panels. For visual clarity, histograms are shown normalized to a common peak value. (E) Plots of simulated spike timing resolution and the theoretically calculated Chapman-Robbins lower bound. The simulations generally do not attain the Chapman-Robbins lower bound, especially for situations with low SNR and slow temporal dynamics. Simulations in panels D and E were done using ν = 20 Hz.