Rank-penalized nonnegative spatiotemporal deconvolution and demixing - PowerPoint PPT Presentation

Rank-penalized nonnegative spatiotemporal deconvolution and demixing of calcium imaging data Eftychios A. Pnevmatikakis 1 Joint work with Tim Machado 1 , Logan Grosenick 2 , Ben Poole 2 , Joshua Vogelstein 3 , and Liam Paninski 1 1 Columbia University 2 Stanford University 3 Duke University

Calcium imaging: Opportunities and challenges ◮ Enables imaging of large neuronal ensembles or detailed dendritic structures. ◮ Numerous applications: Neural connectivity, dendritic computation. ◮ Challenges: noisy slow signal, limited temporal resolution, spatial mixing of active components. Key problem Denoise, unmix, and optimally extract spiking signals from spatiotemporal calcium imaging data.

Spike inference from 1-d fluorescence traces Nonnegative deconvolution (FOOPSI) Assume simple autoregressive dynamics: ∆ C ( t ) = − ∆ t /τ C ( t − 1 ) + As ( t ) Y ( t ) = C ( t ) + ε t Approximate maximum-a-posteriori (MAP) estimate (Vogelstein et. al. 2010): 1 + 1 � � ( Y ( t ) − C ( t )) 2 min s ( t ) 2 σ 2 A λ C , s ≥ 0 t t = 1 � �� � � �� � Loss function Spike sparsity penalty ◮ Outperforms optimal linear (Wiener) filter (nonnegativity is important). ◮ Other approaches: "Peeling" (Grewe et. al. 2010).

Spike inference from 1-d fluorescence traces Nonnegative deconvolution (FOOPSI) Assume simple autoregressive dynamics: ∆ C ( t ) = − ∆ t /τ C ( t − 1 ) + As ( t ) Y ( t ) = C ( t ) + ε t Approximate maximum-a-posteriori (MAP) estimate (Vogelstein et. al. 2010): 1 + 1 � � ( Y ( t ) − C ( t )) 2 min s ( t ) 2 σ 2 A λ C , s ≥ 0 t t = 1 � �� � � �� � Loss function Spike sparsity penalty ◮ Outperforms optimal linear (Wiener) filter (nonnegativity is important). ◮ Other approaches: "Peeling" (Grewe et. al. 2010). Our goal Extend to the high-d setup with multiple (possible overlapping) neurons by sharing information across different pixels optimally.

Single neuron spatiotemporal fluorescence contribution Neuron 1 Spikes Spikes from neuron 1: s 1 ( t )

Single neuron spatiotemporal fluorescence contribution Neuron 1 Spikes Spikes from neuron 1: s 1 ( t ) Neuron 1 Trace Calcium trace from neuron 1: p � C 1 ( t ) = γ j C 1 ( t − j ) + s 1 ( t ) j = 1

Single neuron spatiotemporal fluorescence contribution Neuron 1 Spikes Spikes from neuron 1: s 1 ( t ) Neuron 1 Trace Calcium trace from neuron 1: p � Neuron 1 Location C 1 ( t ) = γ j C 1 ( t − j ) + s 1 ( t ) j = 1 Neuron 1 location: a 1

Single neuron spatiotemporal fluorescence contribution Neuron 1 Spikes Spikes from neuron 1: s 1 ( t ) Neuron 1 Trace Calcium trace from neuron 1: p � Neuron 1 Location Spatiotemporal fluorescence contribution of neuron 1 C 1 ( t ) = γ j C 1 ( t − j ) + s 1 ( t ) j = 1 Pixel Neuron 1 location: a 1 Time Neuron 1 spatiotemporal fluorescence: F 1 ( t ) = a 1 C 1 ( t )

Single neuron spatiotemporal fluorescence contribution Neuron 2 Spikes Spikes from neuron 2: s 2 ( t ) Neuron 2 Trace Calcium trace from neuron 2: p � Neuron 2 Location Spatiotemporal fluorescence contribution of neuron 2 C 2 ( t ) = γ j C 2 ( t − j ) + s 2 ( t ) j = 1 Pixel Neuron 2 location: a 2 Time Neuron 2 spatiotemporal fluorescence: F 2 ( t ) = a 2 C 2 ( t )

Multi-d Case Spatiotemporal dynamics are of low rank Key observation Each neuron contributes a rank-1 term to the spatiotemporal matrix. F = A X C N d d T N T d : # of pixels. T : # of timesteps. N : # of neurons. ◮ rank ( F ) ≤ N ≪ d , T (degrees of freedom ( d + T ) N ≪ dT ) p � CG T = S C i ( t ) = γ j C i ( t − j ) + s i ( t ) Neuron spiking j = 1 N � F ( t ) = a i C i ( t ) F = AC Noiseless image i = 1 Y ( t ) = F ( t ) + ε , Y = F + E Noisy observations

Multi-d Case Rank-penalized denoising ◮ Plain SVD/PCA methods to exploit the low rank structure do not retain nonnegativity of the spikes.

Multi-d Case Rank-penalized denoising ◮ Plain SVD/PCA methods to exploit the low rank structure do not retain nonnegativity of the spikes. ◮ Instead we want to penalize the rank in an structured way: λ 1 � FG T � 1 minimize L ( Y , F ) + + λ NN � F � ∗ F : FG T ≥ 0 � �� � � �� � � �� � Quadratic loss function Spike sparsity penalty Rank penalty � �� � Nonnegative spikes ◮ The nuclear norm (NN) � · � ∗ (sum of the singular values) is the convex relaxation of the rank ( · ) function.

Multi-d Case Rank-penalized denoising ◮ Plain SVD/PCA methods to exploit the low rank structure do not retain nonnegativity of the spikes. ◮ Instead we want to penalize the rank in an structured way: λ 1 � FG T � 1 minimize L ( Y , F ) + + λ NN � F � ∗ F : FG T ≥ 0 � �� � � �� � � �� � Quadratic loss function Spike sparsity penalty Rank penalty � �� � Nonnegative spikes ◮ The nuclear norm (NN) � · � ∗ (sum of the singular values) is the convex relaxation of the rank ( · ) function. ◮ The NN penalty (ideally) shrinks the components due to noise to zero, and provides a robust estimate on the underlying number of neurons. ◮ The NN convex relaxation method has been used successfully in many machine learning applications (e.g. low rank matrix completion (Candès and Recht, 2009)).

Spatiotemporal deconvolution and demixing Nonnegative Matrix Factorization (NMF) We can estimate the spatial and temporal components sequentially using alternating matrix regression: C k + 1 = arg min + λ � CG T � 1 L ( Y , A k C ) Spike deconvolution C : CG T ≥ 0 � �� � � �� � Loss function Spike sparsity A k + 1 = arg min L ( Y , AC k + 1 ) + λ 1 � A � 1 Component demixing A : A ≥ 0 � �� � � �� � Loss function Component sparsity ◮ The spatial component is initialized from the nuclear norm penalized solution using clustering methods. ◮ The sparsity penalties can be consistently estimated using AIC like criteria.

Results: Single-neuron patch ◮ GCaMP6S data, spinal cord neuron, in vitro. ◮ Antidromic stimulation (spike times are known).

Results: Spinal cord synchronized neurons ◮ GCaMP3 data, spinal cord neurons, in vitro.

Results: In vivo cortical dendritic data (rodent barrel cortex, data from Clay Lacefield and Randy Bruno)

Additional details and extensions Bayesian methods ◮ Use sampling methods to quantify uncertainty and estimate parameters. ◮ Can be done efficiently (in just O ( T ) time) in the low-SNR regime, using an augmented block-Gibbs sampler (Martens and Sutskever 2010). Compressive fluorescence imaging ◮ Our methods are applicable to “scanless" imaging approaches that can achieve higher imaging rates (Nikolenko et. al. 2008). ◮ Using compressed sensing ideas, only a few measurements are sufficient. Parameter estimation The parameters are estimated using a method-of-moments approach on the autoregressive dynamics.

Conclusions Structured approaches to calcium imaging denoising problems ◮ Modern statistical approaches provide flexible and powerful methods for dealing with spatiotemporal Ca 2 + imaging denoising problems in a structured way. ◮ Convex optimization leads to efficient MAP inference and can be used to initialize computationally more intensive Bayesian methods (Gibbs sampling, particle filtering). ◮ More broadly, the rapid development of novel experimental methods provides many new challenges and opportunities for breakthroughs based on statistical ideas.

Thank you - Questions?

Extensions: Fully Bayesian approaches Beyond MAP inference Calcium trace estimate 0.5 Raw Data Standard Error 0.4 Mean Gibbs 0.3 MAP 0.2 0.1 0 MAP estimate of Spikes ◮ Use sampling methods to quantify uncertainty and estimate parameters. ◮ Can be done efficiently (in just O ( T ) time) in the low-SNR regime, using an augmented block-Gibbs sampler (Martens Gibbs probabilities of spikes and Sutskever 2010). ◮ Can introduce intermediate timebins to achieve spike superresolution. ◮ Challenge: Perform tractable fully 0 200 400 600 800 1000 1200 1400 Frame Bayesian inference for the general Spike amplitude Baseline Mean spike probability spatiotemporal case. 0.1 0.07 0.2 0.06 0.08 0.05 0.15 0.06 0.04 0.04 0.03 0.1 0.02 0.02 Repetition # Figure : GCaMP3 data, spinal cord neuron in vitro, antidromic stimulation.

Extensions: Compressive fluorescence imaging ◮ Compressive fluorescence imaging, e.g. spatial light modulator (SLM) microscopy (Nikolenko et. al. 2008), takes generalized linear measurements at any frame, allowing for faster imaging rates. ◮ Exploit spike sparsity. # of measurements ≪ # of neurons In high-SNR: n t ∼ O ( pN log ( 1 / p )) ◮ Potential to dramatically increase the number of neurons we can record from given a fixed number of measurements. ◮ Challenge: Accurately estimate the model parameters and apply to real data. Simulated data, 50 neurons.

Parameter estimation γ j and σ 2 are estimated from the autoregressive statistics ◮ For τ > p it holds for the autocovariance C Y : p � C Y ( τ ) = γ j C Y ( τ − j ) . j = 1 ◮ For 0 < τ ≤ p , becomes p � γ j C Y ( τ − j ) − σ 2 γ τ . C Y ( τ ) = j = 1 To estimate b and A λ use A λ E ( y ) = b + 1 − � γ j , and the approximate relation σ ∝ b .