Learning spatially-correlated temporal dictionaries for calcium imaging

TL;DR

This paper introduces a novel dictionary learning approach for calcium imaging data that emphasizes spatial inference and temporal trace detection, improving neuron activity recovery with hierarchical spatial modeling.

Contribution

The paper proposes a new spatially-correlated temporal dictionary learning method that enhances neuron activity detection and inference in calcium imaging data, with hierarchical spatial filtering.

Findings

Improves initialization robustness in neuron detection.

Implicitly infers the number of neurons.

Detects different neuronal types simultaneously.

Abstract

Calcium imaging has become a fundamental neural imaging technique, aiming to recover the individual activity of hundreds of neurons in a cortical region. Current methods (mostly matrix factorization) are aimed at detecting neurons in the field-of-view and then inferring the corresponding time-traces. In this paper, we reverse the modeling and instead aim to minimize the spatial inference, while focusing on finding the set of temporal traces present in the data. We reframe the problem in a dictionary learning setting, where the dictionary contains the time-traces and the sparse coefficient are spatial maps. We adapt dictionary learning to calcium imaging by introducing constraints on the norms and correlations of the time-traces, and incorporating a hierarchical spatial filtering model that correlates the time-trace usage over the field-of-view. We demonstrate on synthetic and real data that our solution has advantages regarding initialization, implicitly inferring number of neurons and simultaneously detecting different neuronal types.

Figures

Equations

$${\mathbf{Y}}^T=\mathbf{A}{\mathbf{\Phi }}^T+{\mathbf{E}}^T,$$

$$\mathbf{Y}=\mathbf{\Phi }{\mathbf{A}}^T+\mathbf{E},$$

$${\mathbf{y}}_{i,j}=\sum _k{\mathbf{\varphi }}_k{a}_{i,j,k}+{\mathbf{ϵ}}_{i,j},$$

$$\arg \underset{\mathbf{\Phi },a_{i,j}}{\min }\left[\sum _{i,j}{\Vert {\mathbf{y}}_{i,j}-\mathbf{\Phi }{\mathbf{a}}_{i,j}\Vert }_2^2+\lambda {\Vert {\mathbf{a}}_{i,j}\Vert }_1\right],$$

$$\begin{array}{rl}\hat{\mathbf{\Phi }}& =\arg \underset{\mathbf{\Phi }\ge 0}{\min }{\Vert \mathbf{Y}-\mathbf{\Phi }\mathbf{A}\Vert }_F^2+{\kappa }_1\Vert \mathbf{\Phi }{\Vert }_F^2\\ & +{\kappa }_2\Vert \mathbf{\Phi }-\tilde{\mathbf{\Phi }}{\Vert }_F^2+{\kappa }_3||{\Phi }^T\Phi -\text{diag}({\Phi }^T\Phi )|{|}_{sav}\end{array}$$

$${\mathbf{y}}_{i,j}$$

$$a_{i,j,k}$$

$${\lambda }_{i,j,k}$$

$${\mathbf{\varphi }}_k$$

$${\hat{\mathbf{a}}}_{i,j}=\arg \underset{\mathbf{a}\ge 0}{\min }\frac{1}{2{\sigma }_y^2}{\Vert {\mathbf{y}}_{i,j}-\mathbf{\Phi }\mathbf{a}\Vert }_2^2+\sum _k{\lambda }_{i,j,k}|a_{i,j,k}|,$$

$${\lambda }_{i,j,k}=\xi {\left(\beta +|a_{i,j,k}|+{\left[|\mathbf{W}\ast {\hat{\mathbf{A}}}_k|\right]}_{i,j}\right)}^{-1},$$

Full text

Index Terms—  Calcium imaging, Dictionary learning, Sparse coding, Two-photon microscopy, Re-weighted $\ell_{1}$

1 Introduction

Calcium imaging, or the optical recording of calcium concentrations in neural tissue, is an important neural imaging technique due to its ability to simultaneously record large neural populations at single cell resolution in awake behaving animals [1, 2, 3]. This technique enables neuroscientists to investigate the role of large cortical areas, uncovering neural representations that are easily lost in single-unit recordings. In particular, calcium imaging via two-photon microscopy (TPM) leverages nonlinear optical interactions to resolve activity at greater depths, enabling minimally-invasive imaging hundreds of microns beneath the surface. The advances in scan rate, achievable field-of-view (FOV), and fluorescence signal strength that have solidified TPM as a fundamental technique in neuroscience, however, have also created a complex data-analysis problem. The high-dimensionality of the datasets ($10$K-$100$K $512\times 512$ frames) and the complex noise environments must be overcome to accurately extract the individual neurons’ time-traces. Volumetric methods that project 3D volumes to 2D images further complicate demixing by producing highly overlapping neuronal structures [4].

Current TPM analysis methods focus on factoring the data cube into a set of spatial profiles111Spatial profiles are sometimes termed Regions of Interest (ROIs), however we prefer the terminology “profile” as it more accurately reflects that these components are 2D projections of neural processes. (the location within the FOV that a neuron occupies) and a corresponding set of time-traces (the fluorescence activity for each component). While many methods have been devised to isolate these factors, based on deep learning [5], active contours [6], spectral embeddings [7], matrix factorization remains the most popular approach [8, 9, 10, 11, 12]. Given an $N_{x}\times N_{y}$ pixel FOV sampled at $T$ time-points, we have

[TABLE]

where $\bm{Y}\in\mathbb{R}^{T\times N_{x}N_{y}}$ is the fluorescence movie data, $\bm{A}\in\mathbb{R}^{N_{x}Ny\times M}$ are the spatial profiles for the $M$ neural components, $\bm{\Phi}\in\mathbb{R}^{T\times M}$ are their corresponding time-traces and $\bm{E}\in\mathbb{R}^{T\times N_{x}N_{y}}$ is the recording noise, often modeled as $i.i.d.$ Gaussian. Constraints are often added to this model to better represent the data, such as sparse compact spatial support, non-negativity and auto-regressive time-trace models [8, 12, 13]. This matrix factorization equally weights both the spatial and temporal components. Ensuing validation likewise often assesses the accuracy of spatial profiles [14].

The main object of interest, however, for scientific inquiry is the set of time-traces, as they are the variables that are linked to behavior, stimuli, learning, etc. Thus, in this work, we reverse the modeling philosophy of looking for spatial components in order to find temporal traces, and rather directly model temporal traces. Mathematically we transpose the original model as

[TABLE]

such that we reverse the roles of $\bm{\Phi}$ and $\bm{A}$. While a seemingly trivial semantic change, we will show this is a more natural model, and this paradigm reduces sensitivity to initializations and removes the burden of morphological post-processing. Philosophically this reorganization places more importance in finding the temporal activity and treats the spatial statistics as secondary, as opposed to typical analyses [15, 16, 13]. This modeling shift is similar to the difference, for example, in modeling the spatial statistics of hyperspectral imaging data versus the predominant spectral end-member analysis.

To infer $\bm{\Phi}$ under the new model, we adopt a dictionary learning (DL) approach, where the dictionary is composed of temporal components and the sparse coefficients are the spatial profiles (Fig. 1B,C). In the temporal domain, we introduce correlation- and continuity-based regularizers on the learned dictionary. This enables implicit inference of the number of neurons, which remains a challenge for many current methods. We further induce implicit spatial contiguity into the DL framework via Reweighted-$l_{1}$ Spatial Filtering (RWL1-SF) in the coefficient inference [17]. We validate our spatially-filtered DL algorithm on both synthetic and real TPM data.

2 Background

Dictionary learning (DL) is an unsupervised method aimed at finding optimal, parsimonious representations for data given exemplar data [18, 19]. In particular, DL decomposes the dataset $\bm{Y}$ in Equation (2) into the dictionary $\bm{\Phi}$ and coefficients $\bm{A}$, under the assumption that the rows of $\bm{A}$ are sparse (i.e., most elements of $\bm{A}$ are zero). While originally used for image processing tasks such as denoising or inpainting [20], DL has been successful in other domains as well.

Hyperspectral Imagery (HSI) has been one such application [21, 22] that has many parallels to TPM. Both are high-dimensional image volumes, where the non-spatial dimension contains important identity information (time-courses in TPM holding neural identity and, optical reflection spectra in HSI holding material identity). Thus HSI processing has focused in large part on learning spectral signatures, using DL for unsupervised extraction of material spectra from data-cubes [21]. In this case the dictionary contains spectral atoms and the sparse coefficient are spatial abundance profiles. Related to the concept of endmembers (the boundaries of the convex hull containing the HSI data [23]), the spectral dictionary both uncovers valuable semantic material information about the imaged area, and is used for many signal processing tasks (e.g. inpainting, super-resolution [24, 17]). To leverage the coarser spatial information, advanced models correlate the pixel-decompositions of neighboring pixels, for example using joint spatial-spectral dictionaries [24], coefficient regularization [25], or spatial filtering [17].

Current applications of DL for TPM data, however, have focused on learning spatial dictionary atoms [15, 16, 13]. For example Pachitariu et al. [15] developed a spatial generative model based on convolutional sparse block coding with the goal of learning the spatial locations of somas and dendrites from the mean image across time. Diego and Hamprecht [16] extend convolutional sparse coding to video data, extracting the spatial components and their temporal activity while estimating a non-uniform and temporally varying background. SCALPEL [13] similarly build a dictionary of spatial components that is refined through iterative merging and clustering, before the temporal components are inferred.

3 Dictionary Learning for Calcium Imaging

Following the inspiration from HSI data processing, we focus on modeling TPM data in terms of the time-traces of single neurons. In particular, we combine DL methodology with the spatially correlated generative model of [17]. This combination enables us to focus on accurate time-trace estimation, reducing errors that may arise when also trying to accurately estimate the spatial maps (i.e., sensitivity to initialization, errors due to mismatched spatial profiles, etc.).

Given a sequence of motion-corrected TPM frames $\bm{Y}_{l}\in\mathbb{R}^{N_{x}\times N_{y}}$, we consider the temporal activity vector $\bm{y}_{i,j}\in\mathbb{R}^{T}$ at each pixel $\{i,j\}$. In any given video, we model $\bm{y}_{i,j}$ as

[TABLE]

where $\bm{\phi}_{k}\in\mathbb{R}^{T}$ for $k\in[1,\cdots,K]$ are the $K$ neural time-traces, $a_{i,j,k}$ is the strength of each neuron’s fluorescence at each pixel, and $\bm{\epsilon}_{i,j}\in\mathbb{R}^{T}$ represents the sensor noise. We define a cost function of the time-traces $\bm{\phi}_{k}$ and the spatial presence coefficients $a_{i,j,k}$, that captures both the data model and the a-priori information that few neurons overlap at any given pixel ($\bm{a}_{i,j}=[a_{i,j,1},\dots,a_{i,j,K}]^{T}\in\mathbb{R}^{K}$ is sparse) and that both $\bm{\Phi}$ and $\bm{A}$ are non-negative. Mathematically, the basic DL problem is

[TABLE]

where $\bm{\Phi}=[\bm{\phi}_{1},\dots,\bm{\phi}_{K}]\in\mathbb{R}^{T\times K}$ is the time-trace dictionary, and $\lambda$ is a parameter that trades off data fidelity and sparsity in $\bm{a}_{i,j}$.

Standard methods for solving Equation (4) alternate between optimizing $\bm{a}_{i,j}$ given an initial $\bm{\Phi}$, which is in turn updated based on the inferred $\bm{a}_{i,j}$. As there is no inherent penalization on $\bm{\Phi}$, a constraint on the norm or maximum value of each $\bm{\phi}_{k}$ is usually imposed to prevent the trivial solution of $\|\bm{\Phi}\|_{F}^{2}\rightarrow\infty$. We extend this basic model for TPM data in two important ways. First build in additional constraints on $\bm{\Phi}$ to isolate single traces despite modeling nonlinear data with a linear-generative model. Second we expand the sparsity model over $\bm{a}_{i,j}$ to include spatial contiguity.

Temporal dictionary learning: We introduce two penalties over $\bm{\Phi}$ to the cost function to facilitate both automatically adapting the number of neural time-traces, and reducing the impact of subtle nonlinearities in the spatial expression of a given time-trace. A Frobenius cost $\|\bm{\Phi}\|_{2}^{2}$ serves to remove unused dictionary elements, and a penalty over intra-dictionary correlations $\|\bm{\Phi}^{T}\bm{\Phi}-\operatorname{diag}(\bm{\Phi}^{T}\bm{\Phi})\|_{sav}=\sum_{i\neq k}\bm{\phi}_{i}^{T}\bm{\phi}_{k}$ penalizes time-traces with trivial differences. To ensure stable convergence, we also include a continuation term that penalizes the change in $\widehat{\bm{\Phi}}$ from the previous estimate $\widetilde{\bm{\Phi}}$. Overall, the update for $\bm{\Phi}$ given $\bm{a}_{i,j,k}$ is

[TABLE]

Spatially filtered sparse coding: To remedy the lack of spatial cohesion in traditional DL, we adapt the re-weighted $\ell_{1}$ spatial filtering (RWL1-SF) previously used to infer sparse, clustered coefficients [17] (Fig. 1C). RWL1-SF, which is an expansion of re-weighted $\ell_{1}$ [26] and the Laplacian-scale mixture model [27], places a hierarchical layer above the Laplacian prior (Fig. 1C). Specifically the data $\bm{Y}$ at each pixel $\{i,j\}$ given $\bm{\Phi}$ is modeled as

[TABLE]

Inferring $\bm{a}_{i,j}$ (via approximate expectation-maximization) under this model amounts to iteratively solving

[TABLE]

for all $\{i,j\}$ and then updating the weights $\lambda_{i,j,k}$ as

[TABLE]

where $\bm{A}_{k}\in\mathbb{R}^{N_{x}\times N_{y}}$ contains the presence coefficients for a single neuron across the FOV and $\bm{W}\ast\widehat{\bm{A}}_{k}$ indicates a 2D convolution. Here $\xi$ and $\beta$ depend only on model parameters $\alpha$, $\theta$, $\sigma_{y}^{2}$ [17].

The weights $\lambda$ in RWL1-SF incorporate spatial information into per-pixel solutions by sharing second-order statistics. Coefficients “activated” in the initial optimization lower the weights for neighboring coefficients, encouraging them to activate in subsequent iterations, whereas non-active coefficients penalize activation with higher weights. The kernel $\bm{W}$ specifies the coefficient influence radius, dictating the neighborhood where interactions are strongest.

4 Experimental Results

**Simulated data: ** We first validate our method on simulated data composed of 14 time-traces generated by convolving a randomly generated weighted spike train with an autoregressive filter with an single pole at -0.7 (Fig. 2A). Spatial presence maps were created via draws from a 2D Gaussian Process which was windowed by a truncated 2D Gaussian filter. One component was left unwindowed to simulate distributed neuropil signals.

We learned time-traces for this data using Algorithm 1 with 16 components, $\kappa_{1}=0.3$, $\kappa_{2}=0.4$, $\kappa_{3}=0.2$, $\xi=2$, and $\beta=0.1$. While the stopping criteria was $\|\Delta\widehat{\bm{\Phi}}\|_{F}^{2}/\|\widehat{\bm{\Phi}}\|_{F}^{2}\leq 10^{-5}$, the algorithm completed in under 20 iterations. Most time-traces converged quickly to a gross estimate and only took longer to remove small, spurious transients (Fig. 2B). The learned dictionary correctly identified the 14 time-traces, and also correctly reduced the extra time-traces to negligible levels (Fig. 2C). These results encouraged us to apply this method to real two-photon calcium imaging data.

Somatic Imaging We ran our approach on a $61\times 46$ pixels by 3000 frames (at 3Hz) dataset from Neurofinder [14]. We initialize the dictionary with 14 components, as compared to the 7 labeled ground-truth neurons. We compare our results with Suite2p [11] which is a modern, state-of-the-art ROI extraction method.

Seven of the resulting learned time-traces (Fig. 3, green) display a close match with the time-traces estimated from the ground-truth labeled pixels (Fig. 3, black) and contain less neuropil induced dips than the Suite2p estimates (Fig. 3, orange). Of the 7 labeled neurons, we found 7 neurons and a background component while Suite2p found 6. Additionally, we detect 3 apical dendrites and a fourth component composed of basal dendrites not included in the provided ground-truth. Detecting these dendrites is possible because of our modeling emphasizing the temporal decompositions rather than relying on specific initialization of spatial components, or a-priori assumptions of neuron’s sizes and shapes.

One interesting case was an apical dendrite located adjacent to a soma (Fig. 4A). Several example frames from the raw video at peaks of activity for one of the components shows that there is indeed an apical (blue) with separate activity from the soma (green). Averaging the activity of the pixels in the ground-truth labeling (black trace) picks up the activity from the dendrite (blue trace). Suite2p is also effected by this overlap, depicting dips in the baseline estimate every time the apical dendrite fires. A second interesting case from a different region of the same Neurofinder dataset (Fig 4B) has our method inferring one dictionary element (green) that is active in two locations that Suite2p marks as two different components (red and burgundy). Example frames from the video confirm that these two components always fire together (Fig 4B).

5 Conclusions

We propose a new approach for extracting neural time-traces from two-photon calcium imaging data. Our approach is based on expanding a dictionary learning algorithm that focuses on time-trace estimation by leveraging weaker spatial information than other methods that rely heavily on (e.g., morphological classification). We validate our method both in simulation and real calcium imaging data. Our method correctly picks out known components in both cases, and furthermore can uncover more subtle features in the real datasets, such as spatial distant pieces of the same neuron. Furthermore, the lack of overly-specific spatial regularization means that our method is applicable to other imaging modalities, such as wide-filed calcium imaging, where exact spatial structure is difficult to quantify.

A number of practical considerations arise in our method, specifically initialization, selecting the number of neurons, and parameter selection. First, we initialize $\bm{\Phi}$ with random values, demonstrating a reduced sensitivity to initialization than other approaches [12]. It may still be possible to improver performance with targeted initializations, however; a possibility that should be explored in future work. Second, the number of dictionary components should be set to more than the expected number of neurons and a background components. The sparsity and Frobenius norms serve to decay unused components (implicitly estimating the number of neurons), but cannot add new elements. Finally, we hand-tuned a number of parameters, e.g., setting the kernel $W$ as a $7\times 7$ Gaussian kernel with variance 3. Future work should address automatic parameter selection (e.g., using variational methods) and sensitivity analysis.