Fast non-negative deconvolution for spike train inference from population calcium imaging

Joshua T. Vogelstein; Adam M. Packer; Tim A. Machado; Tanya Sippy; Baktash Babadi; Rafael Yuste; Liam Paninski

Source-linked AI summary

Fast non-negative deconvolution for spike train inference from population calcium imaging

Joshua T. Vogelstein, Adam M. Packer, Tim A. Machado, Tanya Sippy, Baktash Babadi, Rafael Yuste, Liam Paninski

arXiv:0912.1637v1Dec 9, 2009q-bio.QMq-bio.NC

TL;DR

Calcium imaging produces fluorescence movies, but recovering the underlying spike trains is difficult because of limited signal quality, temporal resolution, unknown parameters, and computational intractability. The paper develops a fast non-negative approximate MAP filter using a log-concave relaxation and efficient banded optimization, yielding more accurate estimates than Wiener filtering while supporting real-time population inference. Spatial filtering and broader model extensions can further refine estimates, although nonlinear generalizations are less robust.

Problem
Recovering the most likely spike trains from calcium fluorescence is challenging because of limited SNR and temporal resolution, unknown parameters, and computational intractability.
Method
The method replaces the Poisson spike prior with an exponential relaxation, enforces non-negative spike estimates, and solves the resulting log-concave problem using efficient optimization.
Results
The fast filter yields more accurate spike-train estimates than optimal linear deconvolution and processes approximately 100-neuron populations faster than real-time.
Takeaways & Limitations
Spatial filtering can be incorporated to further improve estimates, and the filter supports online inference from fluorescence observations.
Takeaways & Limitations
Exact MAP inference remains computationally intractable, and nonlinear fluorescence generalizations are less robust than the linear model.

Abstract

from arXiv · show

Calcium imaging for observing spiking activity from large populations of neurons are quickly gaining popularity. While the raw data are fluorescence movies, the underlying spike trains are of interest. This work presents a fast non-negative deconvolution filter to infer the approximately most likely spike train for each neuron, given the fluorescence observations. This algorithm outperforms optimal linear deconvolution (Wiener filtering) on both simulated and biological data. The performance gains come from restricting the inferred spike trains to be positive (using an interior-point method), unlike the Wiener filter. The algorithm is fast enough that even when imaging over 100 neurons, inference can be performed on the set of all observed traces faster than real-time. Performing optimal spatial filtering on the images further refines the estimates. Importantly, all the parameters required to perform the inference can be estimated using only the fluorescence data, obviating the need to perform joint electrophysiological and imaging calibration experiments.

1 Introduction

Calcium imaging provides fluorescence movies while researchers seek underlying spike trains, making inference difficult under limited SNR, temporal resolution, unknown parameters, and computational constraints. This work introduces a fast approximately maximum-likelihood filter that enforces non-negative spikes and supports population-scale, online inference.

Calcium imaging records fluorescence traces, but the desired neural signal is the spike train of each observable neuron.
Existing approaches include thresholding, template matching, and Wiener filtering, but optimal linear deconvolution does not use the fact that spikes are always positive.
The proposed filter approximates the most likely spike train while restricting inferred spikes to be positive.
50,000 images can be processed in about one second, and filtering roughly 100 neurons runs faster than real-time.
Spatial filtering of the raw movie can be incorporated, and biological datasets demonstrate the filter’s efficacy.

2 Methods

The method models fluorescence as a noisy observation of calcium dynamics driven by spikes, with calcium jumping after spikes and decaying toward baseline. The generative model uses fluorescence, calcium, and Poisson spike-train components while accounting for parameter non-identifiability.

Figure 1 motivates a first-order model in which spikes are convolved with an exponential and combined with Gaussian noise to explain fluorescence.The figure compares measured patch-electrode spikes, the convolved calcium-related signal, and fluorescence residuals.
Fluorescence is modeled as a linear-Gaussian function of intracellular calcium concentration.
Each spike increases calcium by A µM, after which calcium decays toward Cb µM with time constant τ.
The calcium state evolves in discrete frames, with nt denoting the number of spikes in frame t and ∆ equal to the frame duration.
The fluorescence and calcium scales, α and A, and offsets, β and Cb, are pairwise non-identifiable and can be normalized without loss of generality.
Spikes are initially sampled from a Poisson distribution whose expected firing rate per bin is λ∆.

2.2 Goal

The inference target is the MAP spike train given fluorescence, but exact optimization is computationally intractable because it searches over infinitely many integer-valued spike trains. The paper therefore motivates an approximate formulation.

The goal is to find the maximum a posteriori spike train given the observed fluorescence trace.
The posterior combines the fluorescence likelihood and spike-train prior through Bayes’ rule.
The calcium trajectory is deterministic given the spike train, allowing the posterior likelihood to factor through the calcium state.
Exact solution is computationally intractable because it requires nonlinear search over infinitely many spike trains.
Even with at most k spikes per frame, computational complexity scales exponentially with the number of frames as k^T.

2.3 Inferring the most likely spike train, given a fluorescence trace

The inference problem is made tractable by replacing Poisson spikes with an exponential relaxation and optimizing a concave, non-negative objective with a barrier method. The resulting Newton–Raphson algorithm exploits a tridiagonal Hessian to reduce computation to linear time.

2.3 Inferring the most likely spike train, given a fluorescence trace: Replacing the Poisson distribution with an exponential distribution relaxes integer spike counts to non-negative values and makes the optimization concave.The relaxation permits global optimization but can produce partial spikes, which may be thresholded or interpreted probabilistically.
2.3 Inferring the most likely spike train, given a fluorescence trace: The logarithmic barrier term enables optimization under the non-negativity constraint while preserving convergence to the constrained solution as its weight approaches zero.The objective remains concave and twice differentiable, allowing Newton–Raphson updates.
2.3 Inferring the most likely spike train, given a fluorescence trace: The matrix formulation represents the spike-transition constraint with a bidiagonal matrix and enforces non-negative transitions through MC ≥0.The resulting objective combines fluorescence reconstruction, spike-rate regularization, and an element-wise logarithmic barrier.
2.3 Inferring the most likely spike train, given a fluorescence trace: Each Newton–Raphson step updates the calcium estimate using the gradient, Hessian, step direction, and a backtracking line search.The line search selects a step between zero and one that increases the posterior.
2.3 Inferring the most likely spike train, given a fluorescence trace: Linear-time inference replaces cubic-time Hessian inversion because the Hessian is tridiagonal and can be solved with banded Gaussian elimination.The approximation reduces computation from exponential time to linear time in the trace length T.

2.4 Learning the parameters

Unknown model parameters are estimated iteratively from the fluorescence data and inferred spike train. Robust initialization and MAP-based parameter updates allow the procedure to converge after several iterations.

2.4 Learning the parameters: The parameter-learning procedure alternates between estimating the spike train with current parameters and updating parameters from the new estimate.Iterations continue until a convergence criterion is met.
2.4 Learning the parameters: Fluorescence is normalized to [0,1], α is set to one, β is initialized by the median, and σ by the median absolute deviation.γ is initialized from a one-second time constant and λ at 1 Hz.
2.4 Learning the parameters: The intractable integrals for marginal parameter estimation are approximated by maximizing parameters at the MAP estimate of the hidden variables.This approximation is described as a first-order Laplace approximation and is accurate when posterior mass concentrates near the MAP sequence.
2.4 Learning the parameters: The observation parameters α, β, and σ are estimated by separately optimizing the observation-model terms.The scale α can be fixed at one because the scale of calcium is otherwise arbitrary.
2.4 Learning the parameters: The estimated offset is obtained by averaging fluorescence residuals, while σ is the root-mean-square residual error.The rate parameter λ is estimated from the inferred average firing rate.
2.4 Learning the parameters: Parameter iterations stop at an iteration limit or when relative likelihood change falls below a threshold, typically converging after several iterations.

2.5 Spatial filtering

Spatial filtering extends the scalar fluorescence model to the raw pixels in a region of interest. The movie is projected onto a learned spatial filter, reducing inference to a one-dimensional time series while retaining a general pixel-level model.

2.5 Spatial filtering: The preprocessing pipeline segments the movie into regions of interest and projects each ROI’s pixel vector at time t into a scalar fluorescence input.The scalar trace is the assumed input to the filter.
2.5 Spatial filtering: The pixel-level model assigns each pixel its own scale α_x while sharing calcium dynamics and baseline structure across the ROI.Pixel noise is assumed spatially and temporally white with standard deviation σ, though the covariance can be modified.
2.5 Spatial filtering: The generalized model represents the movie as an N_p × T matrix with a spatial filter vector, baseline vector, identity covariance, and Frobenius-norm residual.
2.5 Spatial filtering: Projecting the N_c × T movie onto the spatial filter reduces computation to evaluating a T × 1 vector norm.
2.5 Spatial filtering: The spatial filter is initialized from the median image frame and the offset from the total movie median.Each pixel scale is then estimated by regressing the inferred calcium-plus-offset signal onto that pixel’s fluorescence trace.
2.5 Spatial filtering: Coordinate ascent alternates spatial-filter and offset updates with estimation of the remaining model parameters and spike train.Convergence is defined by stabilization of the spike-train shape because the absolute spatial-filter scale is not identifiable.

2.6 Overlapping spatial filters

The model generalizes to ROIs containing multiple neurons by using vector-valued calcium signals and block-tridiagonal structure. This preserves the linear-time efficiency of the scalar inference algorithm.

2.6 Overlapping spatial filters: Multiple-neuron ROIs are modeled by indexing the N_c neurons within an ROI rather than assuming each ROI contains one neuron.
2.6 Overlapping spatial filters: The generalized model assumes neurons are independent and pixels are conditionally independent and identically distributed given the underlying calcium signals.Each pixel shares standard deviation σ under the stated noise model.
2.6 Overlapping spatial filters: The multi-neuron formulation uses vector-valued calcium and spike variables for the N_c neurons.
2.6 Overlapping spatial filters: The scalar tridiagonal matrix becomes block tridiagonal, but the Thomas algorithm still solves the resulting linear equations in linear time.Thus the efficiency gained from the original tridiagonal structure is maintained for overlapping spatial filters.
2.6 Overlapping spatial filters: Spatial-filter parameters are initialized using a low-rank approximation, with the first N_c principal components providing the initial spatial filters.PCA is selected because fast algorithms are available and it is sufficiently effective and efficient.
2.6 Overlapping spatial filters: The spatial-filter update is solved efficiently in Matlab by applying the estimated baseline parameters across time.

2.7 Experimental Methods

Experiments used calcium imaging of mouse neurons with OGB-1 or Fura-2, while electrophysiological recordings supplied controlled stimulation and spike measurements.

Mouse somatosensory slices 350–400 µm thick were prepared from C57BL/6 mice at postnatal day 14.
Neurons were loaded with 50 µM OGB-1 through recording pipettes or bulk-loaded with Fura-2 AM.
Recordings used a Multiclamp 700B amplifier and National Instruments digitization hardware with custom LabView software.
Current commands delivered square pulses with sufficient amplitude to produce the desired number of action potentials.
Fluorescence traces were extracted by segmenting mean images into ROIs, averaging ROI pixels, removing frequencies below 0.5 Hz, and normalizing values between zero and one.

3 Results

The fast filter produced more accurate non-negative spike-train estimates than the Wiener filter across simulations and biological data, while remaining computationally efficient and extensible.

The fast filter inferred spike trains more accurately than optimal linear deconvolution in simulations and biological recordings.
Non-negative inference removes the negative spikes and ringing permitted by the Wiener filter, producing cleaner estimates.
Parameters estimated from fluorescence alone still yielded fast-filter performance that far outperformed the Wiener filter.
The fast filter could resolve the correct number of spikes within events in OGB-1 recordings, unlike the Wiener-filter output.
A 5,000-frame movie of 100 neurons was analyzed in about ten seconds, versus 100 seconds of acquisition at 50 Hz.
Poisson and time-varying-prior model extensions did not significantly improve filtering results in simulations or data.
The nonlinear Hill-equation model removes the guarantee of posterior log-concavity and global convergence, although linear-model initialization often remained nearly as accurate.
Spatial filtering produced cleaner fluorescence traces and spike inferences when using the true spatial filter.

4 Discussion

The work develops a fast approximate MAP spike-train inference method and extends it toward robust, practical analysis of calcium-imaging data.

4 Discussion: The method approximates the MAP spike train because exact inference is computationally intractable.It uses analytic approximations to make inference computationally feasible.
4 Discussion: The algorithm’s parameters can be estimated from fluorescence observations alone, avoiding joint electrophysiology and imaging calibration.
4 Discussion: Better spatial filtering can produce cleaner fluorescence traces and spike inferences than a typical spatial filter.The simulation compares boxcar filtering with a neuron’s true spatial filter.
4 Discussion: Greedy alternatives yield binary spike trains, whereas this method quickly finds the most likely approximate train but allows partial spikes.
4 Discussion: The fast filter can initialize a more general SMC filter to further improve inference quality.
4 Discussion: A future extension would combine automatic ROI segmentation with the filter to convert raw population movies into populations of spike trains.

B Wiener Filter

This section contrasts Gaussian-based Wiener filtering with the paper’s non-negative inference procedure and describes the efficient optimization implementation.

B Wiener Filter: The Gaussian replacement approximates Poisson spike rates and permits real-valued, including negative, inferred spike counts.
B Wiener Filter: The inference algorithm uses an interior-point method that gradually reduces a barrier parameter while maintaining the non-negative constraint.
B Wiener Filter: Newton-Raphson updates compute the gradient and Hessian, with backtracking line searches used to improve the posterior.
B Wiener Filter: The algorithm updates spatial and temporal model parameters when spatial filtering is enabled.
B Wiener Filter: The Wiener filter is the optimal linear solution under a Gaussian spike-distribution assumption.