PaperlayerView exact revision on arXiv

Source-linked AI summary

Fast Online Deconvolution of Calcium Imaging Data

Johannes Friedrich, Pengcheng Zhou, Liam Paninski

arXiv:1609.00639v3q-bio.NCq-bio.QMstat.AP

TL;DR

Extracting spikes from fluorescence traces is difficult, while real-time and large-scale imaging require faster online processing. The paper develops OASIS, an online active-set method based on PAVA and extended to AR(p) models. OASIS achieves large speed-ups, supports warm starts, and enables real-time spike inference for large imaging datasets, while its model assumes linear-Gaussian fluorescence dynamics.

  • Problem

    Real-time and large-scale calcium imaging require scalable spike inference beyond interior-point methods that cannot handle whole-brain zebrafish data online.

  • Method

    OASIS generalizes PAVA into an online active-set method for sparse non-negative deconvolution and extends it to AR(p) calcium dynamics with warm-started ONNLS refinement.

  • Results

    More than one order of magnitude speed-up over interior-point methods enables real-time online spike inference, with linear scaling of 1 µs/frame and maintained 1–2-order computational-time advantages on long traces.

  • Takeaways & Limitations

    OASIS supports real-time neural activity estimation during imaging and warm-started optimization useful for CNMF and model-hyperparameter adjustment.

  • Takeaways & Limitations

    The method assumes a linear-Gaussian fluorescence forward model, while nonlinear saturation and mean-dependent photon-count noise remain future extensions.

Abstract

from arXiv · show

Fluorescent calcium indicators are a popular means for observing the spiking activity of large neuronal populations, but extracting the activity of each neuron from raw fluorescence calcium imaging data is a nontrivial problem. We present a fast online active set method to solve this sparse non-negative deconvolution problem. Importantly, the algorithm progresses through each time series sequentially from beginning to end, thus enabling real-time online estimation of neural activity during the imaging session. Our algorithm is a generalization of the pool adjacent violators algorithm (PAVA) for isotonic regression and inherits its linear-time computational complexity. We gain remarkable increases in processing speed: more than one order of magnitude compared to currently employed state of the art convex solvers relying on interior point methods. Unlike these approaches, our method can exploit warm starts; therefore optimizing model hyperparameters only requires a handful of passes through the data. A minor modification can further improve the quality of activity inference by imposing a constraint on the minimum spike size. The algorithm enables real-time simultaneous deconvolution of $O(10^5)$ traces of whole-brain larval zebrafish imaging data on a laptop.

1 Introduction

Calcium imaging provides optical access to neural population activity, but recovering spike trains from fluorescence traces remains difficult. The paper motivates scalable online inference for real-time experiments and large whole-brain datasets, and introduces OASIS to address this need.

  • 1 Introduction: Calcium imaging records neural action potentials through optically measurable calcium transients, but extracting each neuron’s activity from raw fluorescence is nontrivial.Existing approaches include template matching, linear deconvolution, and sparse non-negative deconvolution.
  • 1 Introduction: Real-time processing is needed for closed-loop experiments that feed back circuit interventions on physiologically relevant timescales.Such processing can support testing inferred models of neural dynamics, connectivity, and causation in vivo.
  • 1 Introduction: Whole-brain larval zebrafish imaging and CNMF require methods that scale to large datasets because deconvolution occurs inside the CNMF inner loop.CNMF simultaneously denoises, deconvolves, and demixes spatiotemporal calcium imaging data.
  • 1 Introduction: Existing interior-point algorithms process hundreds of neurons in recording time but cannot handle whole-brain zebrafish imaging in real time or exploit warm starts.The paper frames the task as sparse non-negative deconvolution and identifies scalability and hyperparameter optimization as limitations.
  • 1 Introduction: OASIS generalizes PAVA to provide fast, sequential spike inference with at least one order-of-magnitude speed-ups over interior-point methods on simulated and real data.Its online progression enables real-time inference during imaging once neuronal spatial shapes have been identified.

2 Methods

The methods model calcium fluorescence with autoregressive dynamics and solve sparse non-negative deconvolution using online active-set algorithms generalized from isotonic regression. Pooling, warm starts, and extensions to higher-order AR models improve scalability while preserving exactness where specified.

  • 2.4 General AR(p) models: For AR(p) models, the dual active-set method is greedy when p > 1, so ONNLS warm-starts a primal active-set method to obtain the exact result.The higher-order implementation tracks multiple elements per pool because the dynamics are not first-order Markov.
  • 2.1 Model for calcium dynamics: Calcium dynamics are modeled as a stable AR(p) process, with fluorescence generated from calcium concentration plus Gaussian noise.The model typically uses p = 1 or 2 and assumes fixed parameters, with a non-negative scale and offset.
  • 2.3 Online Active Set method to Infer Spikes (OASIS): Pooling updates adjacent variables together, making merges O(1) in pool length and reducing the cost of repeated active-set updates.The pool representation maintains sufficient summary elements rather than explicitly recalculating long pools.
  • 2.3 Online Active Set method to Infer Spikes (OASIS): OASIS generalizes PAVA by replacing monotonicity with AR(1) decay constraints and solves the resulting convex problem exactly in O(T).For γ = 1 and λ = 0, the algorithm reduces to PAVA.
  • 2.3.1 Dual formulation with hard noise constraint: Warm starts reduce hyperparameter-update costs from k1T to k2z, where z is the current number of pools and is typically much smaller than T.After updating λ, a warm-started pass may only verify constraints across the existing pools.
  • 2.4.2 Online non-negative least squares (ONNLS): The ONNLS algorithm runs in O(T) time, with computational cost additionally depending on the number of positive spikes within each considered window.For AR(p), the required matrix solves can alternatively use linear-time Kalman filter-smoother methods.

3 Results

OASIS matched convex solvers on spike extraction while substantially reducing computation time, including for real calcium-imaging data. Its online and warm-start capabilities also improved scalability and enabled efficient hyperparameter optimization, while thresholding introduced an accuracy–timing trade-off.

  • 3.1 Benchmarking OASIS: OASIS extracted spikes well and produced results similar to convex solvers for simulated AR(1), simulated AR(2), and real data.For AR(2), denoised fluorescence traces were virtually identical to those from convex solvers, although OASIS was approximate.
  • 3.1 Benchmarking OASIS: Two magnitudes faster for AR(1) and more than one magnitude faster for AR(2) with known λ than competing convex solvers.For the noise-constrained problem, OASIS still outperformed other algorithms by a huge margin despite taking about three times longer to find λ.
  • 3.1 Benchmarking OASIS: 1 µs/frame proportionality constant confirmed linear scaling up to T = 300,000 frames.Estimating λ from an initial 1,000-frame batch sped inference by a factor of three without compromising quality: correlation was 0.882 ± 0.001 versus 0.881 ± 0.002.
  • 3.1 Benchmarking OASIS: 0.101±0.005 s per trace versus 2.37±0.12 s for the fastest other methods on a 14,400-frame real GCaMP6s dataset.The dataset was recorded at 60 Hz and evaluated with the hard noise constraint and p = 2.
  • 3.1 Benchmarking OASIS: 745 s for all neurons was below the 1,500 s recording duration, whereas ECOS and other candidates required over 25,780 s.This estimate was based on the reported per-trace computation times for the benchmark.
  • 3.2 Hyperparameter optimization: Optimizing the baseline b and AR(1) coefficient γ improved inference quality but increased computational cost.Estimating parameters from a limited number of pools or tenfold-decimated traces reduced cost without harming inference quality.
  • 3.3 Hard thresholding: Exact AR(2) solutions consistently improved spike-train correlation over faster greedy methods.After Gaussian smoothing by one bin, correlations were 0.731 ± 0.008 for the greedy thresholded solution and 0.800 ± 0.007 after ONNLS refinement.

4 Conclusion

The paper presents OASIS as an online active set method for sparse non-negative calcium-imaging deconvolution, while noting current scope limitations and implementation availability.

  • The method assumes a linear-Gaussian forward model for calcium imaging data.
  • A minor algorithm modification imposes a minimum positive spike size, while future work is needed to incorporate nonlinear observation models and refractoriness.
  • The formulation uses an ℓ1 penalty, preserving convexity while achieving an order-of-magnitude speedup over interior point methods and supporting online applications.
  • Python and MATLAB implementations are provided with scripts reproducing the paper’s figures and Table 1.
  • OASIS was added as a deconvolution subroutine to CaImAn for temporal calcium-imaging data extracted and demixed from raw pixels.

S.1 Algorithm for isotonic regression without pooling

The unpooled isotonic-regression implementation advances through the data and backtracks when monotonicity is violated, but its repeated updates are inefficient.

  • The naïve implementation repeatedly updates all values in the backtracked interval using Eq. (7).
  • O(T^2) complexity results because pooling and reuse of previously computed sums are not exploited.Introducing pools addresses both inefficiencies.
  • The algorithm initializes x with y and processes time points sequentially from 2 through T.
  • When x_t′ falls below x_t′−1, the procedure backtracks to merge the violating region.

S.2 Weighted regression

The weighted-regression extension introduces time-point weights to the isotonic-regression formulation and adapts pool initialization and updates accordingly.

  • Weighted regression assigns a weight θ to each time point.
  • Weights can downweight high-variance observations, including fluorescence data with Poissonian photon-count statistics.
  • A nonlinear observation model with saturation effects can be handled by Newton steps whose weighted constrained regressions are solved efficiently with OASIS.
  • For an AR(1) process, weighting modifies Eq. (10) and adds weights to the solution of Eq. (11).
  • Each pool is initialized using the weighted observation, weight, time index, and length, with updates preserving the weighted regression identity.
  • The same weighted-regression modification is also formulated for an AR(p) process.

S.3 Validity of updates according to equations (12-14)

The validity proof establishes by induction that pool updates preserve the weighted and unweighted identities required by the algorithm.

  • Theorem 1 states that updates according to Eqs. (12–14) preserve Eqs. (11, S4) for all pool values.
  • The proof proceeds by induction, beginning with an assumption that the numerator and denominator relations hold.
  • The base case follows from initializing each pool so that Eqs. (S6, S7) hold.
  • The induction step merges two adjacent pools satisfying the relations and shows that the merged pool preserves them.

S.4 Initial calcium fluorescence

The method models elevated initial fluorescence as residual calcium from prior spiking rather than automatically assigning it to a spike at the first time point. For AR(2), this is represented by an exponentially decaying initial pool that can merge with the next pool when its constraint is violated.

  • Positive initial fluorescence in the p = 1 case can produce a positive s1, but the method sets s1 = 0 to represent previous spiking activity.
  • For p > 1, the first pool models exponential decay from prior spiking instead of treating the noisy c1 value as ground truth.
  • The initial pool begins with c1 = v1 and assigns fluorescence values dt−1c1 for t = 1, ..., l1.
  • The first pool merges with the second whenever the constraint v2 ≥ dl1v1 is violated.

S.5 Explicit expressions of the hyperparameter updates for AR(2)

The AR(2) supplement gives explicit formulas for updating the sparsity and baseline parameters while maintaining the noise constraint. It expresses pool values and their updates through recurrence relations, residual equations, and a linear dependence on the hyperparameter increment.

  • The noise-constrained solution increases λ from a small starting value until the noise constraint becomes tight, beginning with an initial pool partition.
  • For AR(2), the penalty components are parameterized as µt = µκt, with µ = λ(1 − γ1 − γ2) except for the final two components.
  • The initial calcium pool uses the decay factor d, while subsequent pool values satisfy the AR(2) recurrence ˆcti+t = γ1ˆcti+t−1 + γ2ˆcti+t−2.
  • All equations define ˆc(µ), and updating µ by ∆µ changes the fitted calcium trace linearly in ∆µ.
  • The update increment is evaluated by inserting the shifted trace into the residual noise constraint and solving the resulting quadratic relation.
  • When the baseline is optimized jointly, the total shift is φ = b + µ, and the optimal baseline is ˆb = ⟨y − ˆc⟩.

Video S1

Video S1 illustrates PAVA as it sweeps through a time series and enforces increasing order constraints across the sequence.

  • The video marks the current pool with blue crosses while the algorithm enforces x1 ≤ ... ≤ xT during its time-series sweep.

Video S2

Video S2 illustrates OASIS for an AR(1) process, showing pool formation, fluorescence updates, and the resulting deconvolved trace during the sweep.

  • The video shows true spike times in red, pools in the shaded background, and the current pool with blue crosses.
  • It displays the evolving calcium fluorescence trace c′ above and the deconvolved trace s′ = Gc′ below as the algorithm runs.
Loading 1609.00639v3…