SCALPEL: Extracting Neurons from Calcium Imaging Data

TL;DR

This paper introduces SCALPEL, a novel method utilizing dictionary learning and sparse group lasso to automatically identify and analyze neurons in calcium imaging videos, facilitating large-scale neural data analysis.

Contribution

The paper presents a new dictionary learning-based approach for neuron extraction from calcium imaging data, implemented in an accessible R package.

Findings

Successfully applied to three calcium imaging datasets.

Provides accurate neuron identification and activity estimation.

Enables scalable analysis of large neural populations.

Abstract

In the past few years, new technologies in the field of neuroscience have made it possible to simultaneously image activity in large populations of neurons at cellular resolution in behaving animals. In mid-2016, a huge repository of this so-called "calcium imaging" data was made publicly-available. The availability of this large-scale data resource opens the door to a host of scientific questions, for which new statistical methods must be developed. In this paper, we consider the first step in the analysis of calcium imaging data: namely, identifying the neurons in a calcium imaging video. We propose a dictionary learning approach for this task. First, we perform image segmentation to develop a dictionary containing a huge number of candidate neurons. Next, we refine the dictionary using clustering. Finally, we apply the dictionary in order to select neurons and estimate their corresponding activity over time, using a sparse group lasso optimization problem. We apply our proposal to three calcium imaging data sets. Our proposed approach is implemented in the R package scalpel, which is available on CRAN.

Equations

$$\mathbf{Y}\approx \mathbf{A}\mathbf{Z},$$

$$\underset{\mathbf{A}\ge 0,\mathbf{Z}\ge 0,{\mathbf{a}}_b\ge 0}{\mathrm{minimize}}\frac{1}{2}\Vert \mathbf{Y}-\mathbf{A}\mathbf{Z}-{\mathbf{a}}_b{\mathbf{z}}_b^{\top }{\Vert }_F^2.$$

$$d_{i,j}^s=1-\frac{p_{i,j}}{\sqrt{p_{i,i}{p}_{j,j}}}.$$

$${\left[{\mathbf{Y}}^B\right]}_{j,k}=\{\begin{array}{ll}{\left[\mathbf{Y}\right]}_{j,k}& \text{if }{\left[\mathbf{Y}\right]}_{j,k}>-{\text{quantile}}_{0.1\%}(\mathbf{Y})\\ 0& \text{otherwise}\end{array}.$$

$$d_{i,j}^t=1-\frac{{\left({\mathbf{a}}_{\cdot ,i}^0\right)}^{\top }{\mathbf{Y}}^B{\left({\mathbf{Y}}^B\right)}^{\top }{\mathbf{a}}_{\cdot ,j}^0}{{\Vert {\left({\mathbf{Y}}^B\right)}^{\top }{\mathbf{a}}_{\cdot ,i}^0\Vert }_2{\Vert {\left({\mathbf{Y}}^B\right)}^{\top }{\mathbf{a}}_{\cdot ,j}^0\Vert }_2}.$$

$$d_{i,j}=\omega d_{i,j}^s+(1-\omega )d_{i,j}^t,$$

$$\underset{\mathbf{Z}\in {\mathbb{R}}^{K_f\times T},\mathbf{Z}\ge 0}{\mathrm{minimize}}\frac{1}{2}{\Vert \mathbf{Y}-{\tilde{\mathbf{A}}}^f\mathbf{Z}\Vert }_F^2+\lambda \alpha \sum _{k=1}^{K_f}{\Vert {\mathbf{z}}_{k,\cdot }\Vert }_1+\lambda (1-\alpha )\sum _{k=1}^{K_f}{\Vert {\mathbf{z}}_{k,\cdot }\Vert }_2,$$

$$\underset{\mathbf{z}\in {\mathbb{R}}^T,\mathbf{z}\ge 0}{\mathrm{minimize}}\frac{1}{2}{\Vert \mathbf{Y}-{\tilde{\mathbf{a}}}^f{\mathbf{z}}^{\top }\Vert }_F^2+\lambda \alpha {\Vert \mathbf{z}\Vert }_1+\lambda (1-\alpha ){\Vert \mathbf{z}\Vert }_2$$

$$\hat{\mathbf{z}}={\left(1-\frac{\lambda (1-\alpha )}{{\Vert {\left({\mathbf{Y}}^{\top }{\tilde{\mathbf{a}}}^f-\lambda \alpha 1\right)}_{+}\Vert }_2}\right)}_{+}{\left(\frac{{\mathbf{Y}}^{\top }{\tilde{\mathbf{a}}}^f-\lambda \alpha 1}{({\tilde{\mathbf{a}}}^f{)}^{\top }{\tilde{\mathbf{a}}}^f}\right)}_{+},$$

$$\underset{{\mathbf{Z}}_{{\mathcal{N}}_s,\cdot }\in {\mathbb{R}}_{+}^{|{\mathcal{N}}_s|\times T}}{\mathrm{minimize}}\frac{1}{2}{\Vert {\mathbf{Y}}_{\mathcal{M}({\mathcal{N}}_s),\cdot }-{\tilde{\mathbf{A}}}_{\mathcal{M}({\mathcal{N}}_s),{\mathcal{N}}_s}^f{\mathbf{Z}}_{{\mathcal{N}}_s,\cdot }\Vert }_F^2+\lambda \alpha \sum _{n\in {\mathcal{N}}_s}{\Vert {\mathbf{z}}_{n,\cdot }\Vert }_1+\lambda (1-\alpha )\sum _{n\in {\mathcal{N}}_s}{\Vert {\mathbf{z}}_{n,\cdot }\Vert }_2,$$

$$\underset{{\mathbf{Z}}_{{\mathcal{N}}_s,\cdot }\ge 0}{\mathrm{minimize}}\frac{1}{2}{\Vert {\mathbf{Y}}_{\mathcal{M}({\mathcal{N}}_s),\cdot }-{\tilde{\mathbf{A}}}_{\mathcal{M}({\mathcal{N}}_s),{\mathcal{N}}_s}^f{\mathbf{Z}}_{{\mathcal{N}}_s,\cdot }\Vert }_F^2+\lambda \alpha \sum _{n\in {\mathcal{N}}_s}{\Vert {\mathbf{z}}_{n,\cdot }\Vert }_1+\lambda (1-\alpha )\sum _{n\in {\mathcal{N}}_s}{\Vert {\mathbf{z}}_{n,\cdot }\Vert }_2,$$

$${\hat{\mathbf{z}}}_{{\mathcal{N}}_s,\cdot }={\left(1-\frac{\lambda (1-\alpha )}{{\Vert {\left(({\mathbf{Y}}_{\mathcal{M}({\mathcal{N}}_s),\cdot }{)}^{\top }{\tilde{\mathbf{a}}}_{\mathcal{M}({\mathcal{N}}_s),{\mathcal{N}}_s}^f-\lambda \alpha 1\right)}_{+}\Vert }_2}\right)}_{+}{\left(\frac{({\mathbf{Y}}_{\mathcal{M}({\mathcal{N}}_s),\cdot }{)}^{\top }{\tilde{\mathbf{a}}}_{\mathcal{M}({\mathcal{N}}_s),{\mathcal{N}}_s}^f-\lambda \alpha 1}{({\tilde{\mathbf{a}}}_{\mathcal{M}({\mathcal{N}}_s),{\mathcal{N}}_s}^f{)}^{\top }{\tilde{\mathbf{a}}}_{\mathcal{M}({\mathcal{N}}_s),{\mathcal{N}}_s}^f}\right)}_{+}.$$

$$\lambda (1-\alpha )\ge {\Vert {\left({\left[({\tilde{\mathbf{A}}}^f{)}^{\top }\mathbf{Y}\right]}_{k,\cdot }-\lambda \alpha 1\right)}_{+}\Vert }_2$$

$$\lambda \ge \underset{k=1,\dots ,K_f}{\max }\left[\min \left(\underset{l=1,\dots ,T}{\max }\frac{{\left({\left[({\tilde{\mathbf{A}}}^f{)}^{\top }\mathbf{Y}\right]}_{k,l}\right)}_{+}}{\alpha },\frac{{\Vert {\left({\left[({\tilde{\mathbf{A}}}^f{)}^{\top }\mathbf{Y}\right]}_{k,\cdot }\right)}_{+}\Vert }_2}{1-\alpha }\right)\right],$$

$$y_{i,j}\equiv \frac{y_{i,j}^0-{\text{median}}_{t=1,...,T}(y_{i,t}^0)}{{\text{median}}_{t=1,...,T}(y_{i,t}^0)+{\text{quantile}}_{10\%}({\mathbf{Y}}^0)},$$

$$d_{i,j}^{union}=1-\frac{p_{i,j}}{p_{i,i}+p_{j,j}-p_{i,j}},d_{i,j}^{min}=1-\frac{p_{i,j}}{\min (p_{i,i},p_{j,j})},\text{ and }{d}_{i,j}^{max}=1-\frac{p_{i,j}}{\max (p_{i,i},p_{j,j})}.$$

$$\hat{\mathbf{Z}}({\lambda }_i)=\underset{\mathbf{Z}\in {\mathbb{R}}^{K_f\times T},\mathbf{Z}\ge 0}{\mathrm{argmin}}\frac{1}{2}{\Vert {\mathbf{Y}}_{\mathcal{T},\cdot }-{\tilde{\mathbf{A}}}_{\mathcal{T},\cdot }^f\mathbf{Z}\Vert }_F^2+{\lambda }_i\alpha \sum _{k=1}^{K_f}{\Vert {\mathbf{z}}_{k,\cdot }\Vert }_1+{\lambda }_i(1-\alpha )\sum _{k=1}^{K_f}{\Vert {\mathbf{z}}_{k,\cdot }\Vert }_2.$$

$$er{r}_{\mathcal{V}}({\lambda }_i)=\frac{1}{|\mathcal{V}|}{\Vert {\mathbf{Y}}_{\mathcal{V},\cdot }^B-{\tilde{\mathbf{A}}}_{\mathcal{V},\cdot }^f\hat{\mathbf{Z}}({\lambda }_i)\Vert }_F^2,$$

$${\lambda }^{\ast }=\underset{{\lambda }_i}{\mathrm{argmax}}\left\{{\lambda }_i:\frac{er{r}_{\mathcal{V}}({\lambda }_i)-{\min }_{{\lambda }_j}er{r}_{\mathcal{V}}({\lambda }_j)}{{\min }_{{\lambda }_j}er{r}_{\mathcal{V}}({\lambda }_j)}\le 0.05\right\};$$

$$\underset{\mathbf{Z}\in {\mathbb{R}}^{K_f\times T},\mathbf{Z}\ge 0}{\mathrm{minimize}}\frac{1}{2}{\Vert \mathbf{Y}-{\tilde{\mathbf{A}}}^f\mathbf{Z}\Vert }_F^2+\frac{{\lambda }^{\ast }}{|\mathcal{T}|/P}\alpha \sum _{k=1}^{K_f}{\Vert {\mathbf{z}}_{k,\cdot }\Vert }_1+\frac{{\lambda }^{\ast }}{|\mathcal{T}|/P}(1-\alpha )\sum _{k=1}^{K_f}{\Vert {\mathbf{z}}_{k,\cdot }\Vert }_2,$$

$$\frac{1}{2}{\Vert \mathbf{y}-\tilde{\mathbf{\beta }}\Vert }_2^2+\lambda {\Vert \tilde{\mathbf{\beta }}\Vert }_2<\frac{1}{2}{\Vert \mathbf{y}-\hat{\mathbf{\beta }}\Vert }_2^2+\lambda {\Vert \hat{\mathbf{\beta }}\Vert }_2.$$

$$\underset{{\mathbf{\beta }}_{\mathcal{C}}\ge 0}{\mathrm{minimize}}\frac{1}{2}{\Vert {\mathbf{y}}_{\mathcal{C}}-{\mathbf{\beta }}_{\mathcal{C}}\Vert }_2^2+\lambda {\Vert {\mathbf{\beta }}_{\mathcal{C}}\Vert }_2.$$

$$\underset{\mathbf{z}\in {\mathbb{R}}^T,\mathbf{z}\ge 0}{\mathrm{minimize}}\frac{1}{2}{\Vert \mathbf{Y}-{\tilde{\mathbf{a}}}^f{\mathbf{z}}^{\top }\Vert }_F^2+\lambda \alpha 1^{\top }\mathbf{z}+\lambda (1-\alpha ){\Vert \mathbf{z}\Vert }_2.$$

$${\Vert \mathbf{Y}-{\tilde{\mathbf{a}}}^f{\mathbf{z}}^{\top }\Vert }_F^2+\lambda \alpha 1^{\top }\mathbf{z}={\Vert \frac{{\mathbf{Y}}^{\top }{\tilde{\mathbf{a}}}^f-\lambda \alpha 1}{\sqrt{({\tilde{\mathbf{a}}}^f{)}^{\top }{\tilde{\mathbf{a}}}^f}}-\sqrt{({\tilde{\mathbf{a}}}^f{)}^{\top }{\tilde{\mathbf{a}}}^f}\mathbf{z}\Vert }_2^2+C,$$

$$\underset{\mathbf{z}\ge 0}{\mathrm{minimize}}\frac{1}{2}{\Vert \frac{{\mathbf{Y}}^{\top }{\tilde{\mathbf{a}}}^f-\lambda \alpha 1}{({\tilde{\mathbf{a}}}^f{)}^{\top }{\tilde{\mathbf{a}}}^f}-\mathbf{z}\Vert }_2^2+\frac{\lambda (1-\alpha )}{({\tilde{\mathbf{a}}}^f{)}^{\top }{\tilde{\mathbf{a}}}^f}{\Vert \mathbf{z}\Vert }_2.$$

$${\Vert \mathbf{Y}-\widetilde{{\mathbf{A}}^f}\mathbf{Z}\Vert }_F^2$$

$$=\sum _{s=1}^S{\Vert {\mathbf{Y}}_{\mathcal{M}({\mathcal{N}}_s),\cdot }-\sum _{s^{\prime }=1}^S{\widetilde{{\mathbf{A}}^f}}_{\mathcal{M}({\mathcal{N}}_s),{\mathcal{N}}_{s^{\prime }}}{\mathbf{Z}}_{{\mathcal{N}}_{s^{\prime }},\cdot }\Vert }_F^2$$

$$=\sum _{s=1}^S{\Vert {\mathbf{Y}}_{\mathcal{M}({\mathcal{N}}_s),\cdot }-{\widetilde{{\mathbf{A}}^f}}_{\mathcal{M}({\mathcal{N}}_s),{\mathcal{N}}_s}{\mathbf{Z}}_{{\mathcal{N}}_s,\cdot }\Vert }_F^2.$$

$$\frac{1}{2}{\Vert {\mathbf{Y}}_{\mathcal{M}({\mathcal{N}}_s),\cdot }-{\tilde{\mathbf{A}}}_{\mathcal{M}({\mathcal{N}}_s),{\mathcal{N}}_s}^f{\mathbf{Z}}_{{\mathcal{N}}_s,\cdot }\Vert }_F^2+\lambda \alpha 1^{\top }{\mathbf{Z}}_{{\mathcal{N}}_s,\cdot }1+\lambda (1-\alpha )\sum _{n\in {\mathcal{N}}_s}{\Vert {\mathbf{z}}_{n,\cdot }\Vert }_2$$

$$f\left({\mathbf{Z}}_{{\mathcal{N}}_s,\cdot }\right)\le f\left({\mathbf{Z}}_{{\mathcal{N}}_s,\cdot }^0\right)+\text{mbox}Tr\left[{\left({\mathbf{Z}}_{{\mathcal{N}}_s,\cdot }-{\mathbf{Z}}_{{\mathcal{N}}_s,\cdot }^0\right)}^{\top }\nabla f\left({\mathbf{Z}}_{{\mathcal{N}}_s,\cdot }^0\right)\right]+\frac{1}{2t}{\Vert {\mathbf{Z}}_{{\mathcal{N}}_s,\cdot }-{\mathbf{Z}}_{{\mathcal{N}}_s,\cdot }^0\Vert }_F^2,$$

$$\underset{{\mathbf{Z}}_{{\mathcal{N}}_s,\cdot }}{\mathrm{minimize}}\frac{1}{2t}{\Vert {\mathbf{Z}}_{{\mathcal{N}}_s,\cdot }-\left({\mathbf{Z}}_{{\mathcal{N}}_s,\cdot }^0-t\nabla f\left({\mathbf{Z}}_{{\mathcal{N}}_s,\cdot }^0\right)\right)\Vert }_F^2.$$

Full text

In this paper, we consider the first step in the analysis of calcium imaging data: namely, identifying the neurons in a calcium imaging video. We propose a dictionary learning approach for this task. First, we perform image segmentation to develop a dictionary containing a huge number of candidate neurons. Next, we refine the dictionary using clustering. Finally, we apply the dictionary in order to select neurons and estimate their corresponding activity over time, using a sparse group lasso optimization problem. We apply our proposal to three calcium imaging data sets.

Our proposed approach is implemented in the R package scalpel, which is available on CRAN.

Keywords: calcium imaging, cell sorting, dictionary learning, neuron identification, segmentation, clustering, sparse group lasso

1 Introduction

The field of neuroscience is undergoing a rapid transformation: new technologies are making it possible to image activity in large populations of neurons at cellular resolution in behaving animals (Ahrens et al., 2013; Prevedel et al., 2014; Huber et al., 2012; Dombeck et al., 2007). The resulting calcium imaging data sets promise to provide unprecedented insight into neural activity. However, they bring with them both statistical and computational challenges.

While calcium imaging data sets have been collected by individual labs for the past several years, up until quite recently large-scale calcium imaging data sets were not publicly-available. Thus, attempts by statisticians to develop methods for the analysis of these data have been hampered by limited data access. However, in July 2016, the Allen Institute for Brain Science released the Allen Brain Observatory, which contains 30 terabytes of raw data cataloguing 25 mice over 360 different experimental sessions (Shen, 2016). This massive data repository is ripe for the development of statistical methods, which can be applied not only to the data from the Allen Institute, but also to calcium imaging data sets collected by individual labs world-wide.

We now briefly describe the science underlying calcium imaging data. When a neuron fires, voltage-gated calcium channels in the axon terminal open, and calcium floods the cell. Therefore, intracellular calcium concentration is a surrogate marker for the spiking activity of neurons (Grienberger and Konnerth, 2012). In recent years, genetically encoded calcium indicators have been developed (Rochefort et al., 2008; Looger and Griesbeck, 2012; Chen et al., 2013). These indicators bind to intracellular calcium molecules and fluoresce. Thus, the location and timing of neurons firing can be seen through a sequence of two-dimensional images taken over time, typically using two-photon microscopy (Svoboda and Yasuda, 2006; Helmchen and Denk, 2005).

A typical calcium imaging video consists of a $500\times 500$ pixels frame over 1 hour, sampled at 15-30 Hz. A given pixel in a given frame is continuous-valued, with larger values representing higher fluorescent intensities due to greater calcium concentrations. An example frame from a calcium imaging video is shown in Figure 1(a). We have posted snippets of the three calcium imaging videos we analyze at www.ajpete.com/software.

On the basis of a calcium imaging video, two goals are typically of interest:

• •

Neuron identification: The goal is to assign pixels of the image frame to neurons. Due to the thickness of the brain slice captured by the imaging technology, neurons can overlap in the two-dimensional image. This means that a single pixel can be assigned to more than one neuron. This step is sometimes referred to as region of interest identification or cell sorting. An example of neurons identified from a calcium imaging video is shown in Figure 1(b).

• •

Calcium quantification: The goal is to estimate the intracellular calcium concentration for each neuron during each frame of the movie. An example of these estimated calcium traces is shown in Figure 1(c).

To a certain extent, these two goals can be accomplished by visual inspection. However, visual inspection suffers from several shortcomings:

• •

It is subjective and it is not reproducible. Two people who view the same video may identify a different set of neurons or different firing times.

• •

It does not yield numerical information regarding neuron firing times, which may be needed for downstream analyses.

• •

It may be inaccurate: for instance, a neuron that is very dim or that fires infrequently may not be identified by visual inspection.

• •

It is not feasible on videos with very large neuronal populations or very long durations. In fact, a typical calcium imaging video contains 250,000 pixels and more than 50,000 frames, making visual inspection essentially impossible.

In this paper, we propose a method that identifies the locations of neurons, and estimates their calcium concentrations over time. Previous proposals to automatically accomplish these tasks have been proposed in the literature, and are reviewed in Section 3. However, our method has several advantages over competing approaches. Unlike many existing approaches, it:

• •

Involves few tuning parameters, which are themselves interpretable to the user, can for the most part be set to default values, and can be varied independently;

• •

Yields results that are stable across a range of tuning parameters;

• •

Is computationally feasible even on very large data sets; and

• •

Uses spatial and temporal information to resolve individual neurons, from sets of overlapping neurons, without post-processing.

The methods proposed in this paper can be seen as a necessary step that precedes downstream modeling of calcium imaging data. For instance, there is substantial interest in modeling functional connectivity among populations of neurons, or using neural activity to decode stimuli (see, e.g., Ko et al. (2011); Mishchencko et al. (2011); Paninski et al. (2007)). However, before either of those tasks can be carried out, it is necessary to first identify the neurons and determine their calcium concentrations; those are the tasks that we consider in this paper.

The remainder of this paper is organized as follows. We introduce notation in Section 2. In Section 3, we review related work. We present our proposal in Section 4, and discuss the selection of tuning parameters in Section 5. We apply our method to three calcium imaging videos in Section 6. We discuss the technical details of our proposal in Section 7. We close with a discussion in Section 8. Proofs are in the appendix.

2 Notation

Let $P$ denote the total number of pixels per image frame, and $T$ the number of frames of the video. We define ${\bm{Y}}$ to be a $P\times T$ matrix for which the $(i,j)$th element, $y_{i,j}$, contains the fluorescence of the $i$th pixel in the $j$th frame. We let ${\bm{y}}_{i,\cdot}=\begin{pmatrix}y_{i,1}&y_{i,2}&\ldots&y_{i,T}\end{pmatrix}^{\top}$ represent the fluorescence of the $i$th pixel at each of the $T$ frames. We let ${\bm{y}}_{\cdot,j}=\begin{pmatrix}y_{1,j}&y_{2,j}&\ldots&y_{P,j}\end{pmatrix}^{\top}$ represent the fluorescence of all $P$ pixels during the $j$th frame. We use the same subscript conventions in order to reference the elements, rows, and columns of other matrices.

The goal of our work is to (1) identify the locations of the neurons; and (2) quantify calcium concentrations for these neurons over time. We view these tasks in the framework of a matrix factorization problem: we decompose ${\bm{Y}}$ into a matrix of spatial components, ${\bm{A}}\in\mathbb{R}^{P\times K}$, and a matrix of temporal components, ${\bm{Z}}\in\mathbb{R}^{K\times T}$, such that

[TABLE]

where $K$ is the total number of estimated neurons. Note that ${\bm{a}}_{\cdot,k}$ specifies which of the $P$ pixels of the image frame are mapped to the $k$th neuron, and ${\bm{z}}_{k,\cdot}$ quantifies the calcium concentration for the $k$th neuron at each of the $T$ video frames. We note that the true number of neurons is unknown, and must be determined as part of the analysis.

3 Related Work

There are two distinct lines of work in this area. The first focuses on simply identifying the regions of interest in the video, and then subsequently estimating the calcium traces. The second aims to simultaneously identify neurons and quantify their calcium concentrations.

Methods that focus solely on region of interest identification typically construct a summary image for the calcium imaging video and then segment this image using various approaches. For example, Pachitariu et al. (2013) calculate the mean image of the video and then apply convolutional sparse block coding to identify regions of interest. Alternatively, Smith and Häusser (2010) calculate a local cross-correlation image, which is then thresholded using a locally adaptive filter to extract the regions of interest. Similar approaches have been used by others (Ozden et al., 2008; Mellen and Tuong, 2009). These region of interest approaches often do not fully exploit temporal information, nor do they handle overlapping neurons well due to using summaries aggregated over time.

We now focus on methods proposed to accomplish both goals, neuron identification and calcium quantification, simultaneously. One of the first automatic methods in this area was proposed by Mukamel et al. (2009). This method first applies principal component analysis to reduce the dimensionality of the data, followed by spatio-temporal independent component analysis to produce spatial and temporal components that are statistically independent of one another. Though this method is widely used, it often requires heuristic post-processing of the spatial components, and typically fails to distinguish between spatially overlapping neurons (Pnevmatikakis et al., 2016).

To better handle overlapping neurons, Maruyama et al. (2014) proposed a non-negative matrix factorization approach, which estimates ${\bm{A}}$ and ${\bm{Z}}$ in (1) by solving

[TABLE]

In (2), the term ${\bm{a}}_{b}{\bm{z}}_{b}^{\top}$ is a rank-one correction for background noise: ${\bm{z}}_{b}\in\mathbb{R}^{T}$ is a temporal representation of the background noise (known as the bleaching line, and estimated using a linear fit to average fluorescence over time of a background region), and ${\bm{a}}_{b}\in\mathbb{R}^{P}$ is a spatial representation of the background noise. Element-wise positivity constraints are imposed on ${\bm{A}}$, ${\bm{Z}}$, and ${\bm{a}}_{b}$ in (1). While (2) can handle overlapping neurons, there are no constraints on the sparsity of ${\bm{A}}$ or ${\bm{Z}}$, or locality constraints on ${\bm{A}}$. Thus, the estimated temporal components are very noisy, and the estimated spatial components are often not localized and must be post-processed heuristically.

To overcome these shortcomings, Haeffele et al. (2014) modify (2) so that the temporal components ${\bm{A}}$ are sparse and the spatial components ${\bm{Z}}$ are sparse and have low total variation. Recently, Pnevmatikakis et al. (2016) further refine (2) by explicitly modeling the dynamics of the calcium when estimating the temporal components ${\bm{Z}}$, and combining a sparsity constraint with intermediate image filtering when estimating the spatial components ${\bm{A}}$. Zhou et al. (2016) extend the work of Pnevmatikakis et al. (2016) to better handle one-photon imaging data by (1) modeling the background in a more flexible way and (2) introducing a greedy initialization procedure for the neurons that is more robust to background noise. Related approaches are taken by Diego et al. (2013); Diego Andilla and Hamprecht (2014); Friedrich et al. (2015). Other recent approaches consider using convolutional networks trained on manual annotation (Apthorpe et al., 2016) and multi-level matrix factorization (Diego Andilla and Hamprecht, 2013).

While these existing approaches show substantial promise, and are a marked improvement over visual identification of neurons from the calcium imaging videos, they also suffer from some shortcomings:

• •

The optimization problems (see, e.g., (2)) are biconvex. Thus, algorithms typically get trapped in unattractive local optima. Furthermore, the results strongly depend on the choice of initialization.

• •

Each method involves several user-selected tuning parameters. There is no natural interpretation to these tuning parameters, which leads to challenges in selection. Furthermore, changing one tuning parameter may necessitate updating all of them. Moreover, there is no natural nesting with respect to the tuning parameters: a slight increase or decrease in one tuning parameter can lead to a completely different set of identified neurons.

• •

The number of neurons $K$ must be specified in advance, and the estimates obtained for different values of $K$ will not be nested: two different values of $K$ can yield completely different answers.

• •

Post-processing of the identified neurons is often necessary.

• •

Implementation on very large data sets can be computationally burdensome.

To overcome these challenges, instead of simultaneously estimating ${\bm{A}}$ and ${\bm{Z}}$ in the model (1), we take a dictionary learning approach. We first leverage spatial information in order to build a preliminary dictionary of spatial components, which is then refined using a clustering approach to give an estimate of ${\bm{A}}$. We then use our estimate of ${\bm{A}}$ in order to obtain an accurate estimate of the temporal components ${\bm{Z}}$, while simultaneously selecting the final set of neurons in ${\bm{A}}$. This dictionary learning approach allows us to re-cast (1), a very challenging unsupervised learning problem, into a much easier supervised learning problem. Compared to existing approaches, our proposal is much faster to solve computationally, involves more interpretable tuning parameters, and yields substantially more accurate results.

4 Proposed Approach

Our proposed approach is based on dictionary learning. In Figure 2, we summarize our Segmentation, Clustering, and Lasso Penalties (SCALPEL) proposal, which consists of four steps:

1. Step 0.

Data Pre-Processing: We apply standard pre-processing techniques in order to smooth the data both temporally and spatially, remove the bleaching effect, and calculate a standardized fluorescence. Details are provided in Appendix A.1. In what follows, ${\bm{Y}}$ refers to the calcium imaging data after these three pre-processing steps have been performed. 2. Step 1.

Construction of a Preliminary Spatial Component Dictionary: We apply a simple image segmentation procedure to each frame of the video in order to derive a spatial component dictionary, which is used to construct the matrix ${\bm{A}}^{0}\in\mathbb{R}^{P\times K^{0}}$ with $a^{0}_{j,k}=1$ if the $j$th pixel is contained in the $k$th preliminary dictionary element, and $a^{0}_{j,k}=0$ otherwise. This is discussed further in Section 4.1. 3. Step 2.

Refinement of the Spatial Component Dictionary: To eliminate redundancy in the preliminary spatial component dictionary, we cluster together dictionary elements that co-localize in time and space. This results in a matrix ${\bm{A}}\in\mathbb{R}^{P\times K}$: $a_{j,k}=1$ if the $j$th pixel is contained in the $k$th dictionary element, and $a_{j,k}=0$ otherwise. More details are provided in Section 4.2. 4. Step 3.

Application of the Spatial Component Dictionary: We filter out dictionary elements corresponding to clusters with few members, resulting in ${\bm{A}}^{f}\in\mathbb{R}^{P\times{K_{f}}}$, which contains a subset of the columns of ${\bm{A}}$. We then estimate the temporal components ${\bm{Z}}$ corresponding to the filtered elements of the dictionary by solving a sparse group lasso problem with a non-negativity constraint. The $k$th row of $\hat{\bm{Z}}$ is the estimated calcium trace corresponding to the $k$th filtered dictionary element. For some values of the tuning parameters, entire rows of $\hat{\bm{Z}}$ may equal zero. Thus, we are simultaneously selecting the spatial components and estimating the temporal components. Additional details are in Sections 4.3 and 7.

Step 1 is applied to each frame separately, and thus can be efficiently performed in parallel across the frames of the video. Similarly, parts of Step 0 can be parallelized across frames and across pixels.

Throughout this section, we illustrate SCALPEL on an example calcium imaging data set that has $205\times 226$ pixels and 3000 frames. Figures 1-8 and Figures 11-13, as well as Figures 14, 15, 17, and 18 in the appendix, involve this data set. In Section 6, we present a more complete analysis of this data set, along with analyses of additional data sets.

4.1 Step 1: Construction of a Preliminary Spatial Component Dictionary

In this step, we identify a large set of preliminary dictionary elements, by applying a simple image segmentation procedure to each frame separately.

Threshold Image: We create a binary image by thresholding the image frame. Figure 3(b) displays the binary image that results from thresholding the frame shown in Figure 3(a). 2. 2.

Identify Connected Components: We identify the connected components of the thresholded image, using the notion of 4-connectivity: connected pixels are pairs of white pixels that are immediately to the left, right, above, or below one another (Sonka et al., 2014). Some of these connected components may represent neurons, whereas others are likely to be noise artifacts or snapshots of multiple nearby neurons. 3. 3.

Filter Components: To eliminate noise, we filter components based on their overall size, width, and height. In the examples in this paper, we discard connected components of fewer than 25 or more than 500 pixels, as well as those with a width or height larger than 30 pixels.

We now discuss the choice of threshold used above. After performing Step 0, we expect that the intensities of “noise pixels” (i.e., pixels that are not part of a firing neuron in that frame) will have a distribution that is approximately symmetric and approximately centered at zero. In contrast, non-noise pixels will have larger values. This implies that the noise pixels should have a value no larger than the negative of the minimum value of ${\bm{Y}}$. Therefore, we threshold each frame using the negative of the minimum value of ${\bm{Y}}$. We also repeat this procedure using a threshold equal to the negative of the 0.1% quantile of ${\bm{Y}}$, as well as with the average of these two threshold values. In Section 5.1, we discuss alternative approaches to choosing this threshold.

The $K^{0}$ connected components that arise from performing Step 1 on each frame, for each of the three threshold values, form a preliminary spatial component dictionary. We use them to construct the matrix ${\bm{A}}^{0}\in\mathbb{R}^{P\times K^{0}}$: the $k$th column of ${\bm{A}}^{0}$ is a vector of $1$’s and [math]’s, indicating whether each pixel is contained in the $k$th preliminary dictionary element.

4.2 Step 2: Refinement of the Spatial Component Dictionary

We will now refine the preliminary spatial component dictionary obtained in Step 1, by combining dictionary elements that are very similar to each other, as these likely represent multiple appearances of a single neuron. We proceed as follows:

Calculate Dissimilarity Matrix: We use a novel dissimilarity metric, which incorporates both spatial and temporal information, in order to calculate the dissimilarity between every pair of dictionary elements. More details are given in Section 4.2.1. 2. 2.

Perform Prototype Clustering: We use the aforementioned pair-wise dissimilarities in order to perform prototype clustering of dictionary elements (Bien and Tibshirani, 2011). We also identify a representative dictionary element for each cluster. More details are given in Section 4.2.2.

These elements of this refined dictionary make up the columns of the matrix ${\bm{A}}$, which will be used in Step 3, discussed in Section 4.3.

4.2.1 Choice of Dissimilarity Metric

Before performing clustering, we must decide how to quantify similarity between the $K^{0}$ elements of the preliminary dictionary obtained in Step 1. Dictionary elements that correspond to the same neuron are likely to have (1) similar spatial maps and (2) similar average fluorescence over time. To this end, we construct a dissimilarity metric that leverages both spatial and temporal information.

We define $p_{i,j}=({\bm{a}}^{0}_{\cdot,i})^{\top}{\bm{a}}^{0}_{\cdot,j}$, the number of pixels shared between the $i$th and $j$th dictionary elements. When $i=j$, $p_{i,i}$ is simply the number of pixels in the $i$th dictionary element. We then define the spatial dissimilarity between the $i$th and $j$th dictionary elements to be

[TABLE]

Thus, $d^{s}_{i,j}=1$ if and only if the $i$th and $j$th elements are non-overlapping in space, and $d^{s}_{i,j}=0$ if and only if they are identical. Note that $d^{s}_{i,j}$ is known as the cosine dissimilarity or Ochiai coefficient (Gower, 2006). Alternatives to (3) are discussed in Appendix A.2.

We now define the matrix ${\bm{Y}}^{B}$, a thresholded version of the pre-processed data matrix ${\bm{Y}}$ (obtained in Step 0), with elements of the form

[TABLE]

Note that when a value other than the negative of the 0.1% quantile is used for image segmentation in Step 1, this value can also be used to threshold ${\bm{Y}}$ above. The temporal dissimilarity between the $i$th and $j$th dictionary elements is defined as

[TABLE]

(Note that the elements of $\left({\bm{Y}}^{B}\right)^{\top}{\bm{a}}^{0}_{\cdot,i}\in\mathbb{R}^{T}$ represent the thresholded fluorescence of each time frame, summed over all pixels in the $i$th preliminary dictionary element.) We threshold ${\bm{Y}}$ before computing this dissimilarity, because (1) we are interested in the extent to which there is agreement between the peak fluorescences of the $i$th and $j$th preliminary dictionary elements; and (2) the sparsity induced by thresholding is computationally advantageous.

Finally, the overall dissimilarity is

[TABLE]

where $\omega\in[0,1]$ controls the relative weightings of the spatial and temporal dissimilarities. We use $\omega=0.2$ to obtain the results shown throughout this paper. In Figure 4, we illustrate pairs of preliminary dictionary elements with various dissimilarities for $\omega=0.2$.

4.2.2 Prototype Clustering

We now consider the task of clustering the elements of the preliminary dictionary. To avoid pre-specifying the number of clusters, and to obtain solutions that are nested as the number of clusters is varied, we opt to use hierarchical clustering (Hastie et al., 2009).

In particular, we use prototype clustering, proposed in Bien and Tibshirani (2011), with the dissimilarity given in (4). Prototype clustering guarantees that at least one member of each cluster has a small dissimilarity with all other members of the cluster. To represent each cluster using a single dictionary element, we choose the member with the smallest median dissimilarity to all of the other members. Then we combine the representatives of the $K$ clusters in order to obtain a refined spatial component dictionary. We can represent this refined dictionary with the matrix, ${\bm{A}}\in\mathbb{R}^{P\times K}$, as follows: $a_{j,k}=1$ if the $j$th pixel is contained in the $k$th cluster’s representative, and $a_{j,k}=0$ otherwise.

We apply prototype clustering to the example calcium imaging data set, using the R package protoclust (Bien and Tibshirani, 2015). The resulting dendrogram is in Figure 5(a). In Section 5.2, we discuss choosing the cut-point, or height, at which to cut the dendrogram. Results for different cut-points are displayed in Figures 5(b)–(e). An additional example is provided in Appendix A.3.

4.3 Step 3: Application of the Spatial Component Dictionary

In this final step, we optionally filter the $K$ refined dictionary elements, and then estimate the temporal components associated with this filtered dictionary. There are a couple of ways to perform the optional filtering of the dictionary elements: (1) we can filter based on the number of members in the cluster, as clusters with a larger number of members are more likely to be true neurons, or (2) we can manually filter the elements. This filtering process is discussed more in Appendix A.5. After this filtering, we construct the filtered dictionary, ${\bm{A}}^{f}\in\mathbb{R}^{P\times{K_{f}}}$, which contains the retained columns of ${\bm{A}}$.

We estimate the temporal components associated with the $K_{f}$ elements of the final dictionary by solving a sparse group lasso problem with a non-negativity constraint,

[TABLE]

where $\alpha\in[0,1]$ and $\lambda>0$ are tuning parameters, and $\tilde{\bm{A}}^{f}$ is defined as $\tilde{\bm{a}}^{f}_{\cdot,k}={\bm{a}}^{f}_{\cdot,k}/\|{\bm{a}}^{f}_{\cdot,k}\|_{2}^{2}$. This scaling of ${\bm{A}}^{f}$ is justified in Section 7.3.

The first term of the objective in (5) encourages the spatiotemporal factorization (1) to fit the data closely. The second term in (5) encourages the temporal components to be sparse, so that each neuron is estimated to be active in a small number of frames (Tibshirani, 1996). The third term in (5) encourages group-wise sparsity on the rows of ${\bm{Z}}$ (Yuan and Lin, 2006), which leads to selection of dictionary components: for $\lambda(1-\alpha)$ sufficiently large, only a subset of the ${K_{f}}$ elements in the filtered dictionary will have a non-zero temporal component.

In Section 7, we discuss our algorithm for solving (5), justify the scaling of ${\bm{A}}^{f}$ used in (5), and present conditions under which the solution to (5) is sparse. Readers simply interested in the practical use of our method may disregard this section.

5 Tuning Parameter Selection

SCALPEL involves a number of tuning parameters:

• •

Step 1: Quantile thresholds for image segmentation

• •

Step 2: Cut-point for dendrogram

• •

Step 3: $\lambda$ and $\alpha$ for Equation 5

However, in marked contrast to competing methods, the tuning parameters in SCALPEL can be chosen independently at each step and are very interpretable to the user. Furthermore, we strongly recommend the use of default values for all but the choice of $\lambda$ in Step 3. Therefore, in practice, SCALPEL only involves one tuning parameter that must be selected by the user. We discuss this at greater length below.

5.1 Tuning Parameters for Step 1

The tuning parameters in Step 1 are the quantile thresholds used in image segmentation. In analyzing the example video considered thus far in this paper, we used three different threshold values to segment the video: the negative of the minimum value of ${\bm{Y}}$, the negative of the 0.1% quantile of ${\bm{Y}}$, and the average of these two values. In principle, these quantile thresholds may need to be adjusted; however, it is straightforward to select reasonable thresholds by visually examining a few frames and their corresponding binary images, as in Figures 3(a) and (b).

5.2 Tuning Parameters for Step 2

In Step 2, we must choose a height $h\in[0,1]$ at which to cut the hierarchical clustering dendrogram, as shown in Figure 5(a). Fortunately, the cut-point has an intuitive interpretation, which can help guide our choice. If we cut the dendrogram at a height of $h$, then each cluster will contain an element of the preliminary dictionary that has dissimilarity no more than $h$ with each of the other members of that cluster. Figure 4 displays pairs of preliminary dictionary elements with a given dissimilarity. We have used a fixed cut-point of $0.18$ in order to obtain all of the results shown in this paper. An investigator can either choose a cut-point by visual inspection of the resulting refined dictionary elements, or can simply use a fixed cut-point, such as $0.18$. We recommend choosing a cut-point less than the dissimilarity weight $\omega$ in (4), as this guarantees that the dictionary elements within each cluster will have spatial overlap. Note that we do not consider $\omega$ to be a tuning parameter, as we recommend keeping it fixed at $\omega=0.2$.

5.3 Tuning Parameters for Step 3

In Step 3, we also must choose values of $\lambda$ and $\alpha$ in (5). We have had empirical success using a fixed value of $\alpha$ near 1, as we expect each neuron to spike a small number of times. We used $\alpha=0.9$ to obtain all of the results in this paper.

Since (5) is a supervised learning problem, cross-validation (or a closely-related validation set approach) is a natural approach to select the value of $\lambda$. It is computationally feasible, because solving (5) is actually quite efficient, as detailed in Section 7.2. Furthermore, Corollary 7.5 makes it easy to determine the largest value of $\lambda$ that should be considered in our cross-validation scheme. Contrary to standard practice, we use a thresholded version of ${\bm{Y}}$ when calculating the validation error, as we wish to achieve small reconstruction error on the brightest parts of the video, which align with times when neurons are active. We illustrate this approach in Section 6, and provide further details in Appendix A.4.

As an alternative to cross-validation, the distribution of ${\bm{Y}}$ provides an efficient way to choose $\lambda$. We can think of the elements of ${\bm{Y}}$ as corresponding to noise pixels (i.e., non-neuronal pixels or neuronal pixels in non-firing frames) as well as signal pixels. Due to the pre-processing performed in Step 0, the distribution of the noise pixels is approximately symmetric and approximately centered around 0. The distribution of the signal pixels has a much larger center. Therefore, it is likely that many signal pixels, and very few noise pixels, will exceed the negative of (for instance) the 0.1% quantile of ${\bm{Y}}$. The closed-form solution for non-overlapping neurons given in (7) implies that a frame of $\hat{\bm{z}}$ is zeroed out if its average fluorescence is less than $\lambda\alpha$. Thus, we choose $\lambda=\frac{-\text{quantile}_{0.1\%}({\bm{Y}})}{\alpha}$.

Alternatively, we can choose the value of $\lambda$ by visually inspecting the calcium traces (that is, the rows of $\hat{\bm{Z}}$), as shown for instance in Figure 1(c), at each value of $\lambda$.

Any of these approaches can be used to choose the value of $\lambda$ in (5). However, in light of the fact that (5) decomposes into groups of non-overlapping neurons as detailed in Section 7.2, it is also possible to choose a different value of $\lambda$ for each group of overlapping neurons. This corresponds to a straightforward modification of the optimization problem (5).

6 Results for Calcium Imaging Data

In this section, we compare SCALPEL’s performance to that of competitor methods on three calcium imaging data sets. In particular, we consider one-photon and two-photon calcium imaging data sets in Sections 6.1 and 6.2, respectively. In Section 6.3, we compare the run times for the various methods. Additional plots regarding the analyses of these three calcium imaging videos are available at www.ajpete.com/software.

6.1 Application to One-Photon Calcium Imaging Data

We now present the results from applying SCALPEL to the calcium imaging video used as an example in Section 4. This one-photon video, collected by the lab of Ilana Witten at the Princeton Neuroscience Institute, has 3000 frames of size $205\times 226$ pixels sampled at 10 Hz. We used the default tuning parameters to analyze this video. Using the default quantile thresholds that corresponded to thresholds of 0.0544, 0.0743, and 0.0942, Step 1 of SCALPEL resulted in a preliminary dictionary with 2943 elements, which came from 997 different frames of the video. Using a cut-point of 0.18 in Step 2 resulted in a refined dictionary that contained 50 elements. In Step 3, we discarded the twenty-one components corresponding to clusters with fewer than 5 preliminary dictionary elements assigned to them, and then fit the sparse group lasso model with $\alpha=0.9$ and $\lambda=0.0416$, which was chosen using the validation set approach described in Appendix A.4. The results are shown in Figures 6(a) and (b). In Figure 6(c), we compare the estimated neurons to a pixel-wise variance plot of the calcium imaging video. We expect pixels that are part of true neurons to have higher variance than pixels not associated with any neurons. Indeed, we see that many of the estimated neurons coincide with regions of high pixel-wise variance. However, some estimated neurons are in regions with low variance. Examining the frames from which the dictionary elements were derived can provide further evidence as to whether an estimated neuron is truly a neuron. For example, in Figure 7(a), we show that one of the estimated neurons in a low-variance region does indeed appear to be a true neuron, while Figure 7(b) shows evidence that one of the estimated neurons is not truly a neuron.

Instead of filtering on the basis of cluster size, we can manually classify the elements of the refined dictionary from Step 2 before proceeding to Step 3. In order to decide whether to keep or discard an estimated neuron, we can examine the frames from which the corresponding preliminary dictionary elements were extracted. We performed this process for the 50 dictionary elements from Step 2, and found evidence that 32 of the 50 elements correspond to true neurons. Compared to requiring at least 5 members in a cluster, this manual process produced an additional 4 true neurons, while excluding an element with 11 members that did not appear to correspond to a true neuron (element 22 in Figure 6). The spatial and temporal results based on this manual filtering are shown in Figures 1(b) and (c).

We compare the performance of SCALPEL to that of CNMF-E (Zhou et al., 2016), which is a proposal for the analysis of one-photon data that takes a matrix factorization approach, as described in Section 3. The tuning parameters we consider are those noted in Algorithm 1 of Zhou et al. (2016): the average neuron size $r$ and the width of the 2D Gaussian kernel $\sigma$, which relate to the spatial filtering, and the minimum local correlation $c_{\text{min}}$ and the minimum peak-to-noise ratio $\alpha_{\text{min}}$, which relate to initializing neurons. We choose $r=11$ in accordance with the average diameter of the neurons identified using SCALPEL. The default values suggested for the other tuning parameters are $\sigma=3$, $c_{\text{min}}=0.85$, and $\alpha_{\text{min}}=10$. We present the results for these default values in Figure 8(a). Only 14 of the 32 neurons identified using SCALPEL were found by CNMF-E using these default parameters. In order to increase the number of neurons found, we consider lower values for $c_{\text{min}}$ and $\alpha_{\text{min}}$. We fit CNMF-E for all combinations of $c_{\text{min}}=0.5,0.6,0.7$ and $\alpha_{\text{min}}=3,5,7$. To assess the performance of these 9 combinations of tuning parameters, we reviewed each estimated neuron for evidence of whether or not it appeared to be a true neuron, based on reviewing the frames in which the neuron was estimated to be most active. We also determined instances where multiple estimated neurons corresponded to the same true neuron. In Figure 8(b), we present the fit, chosen from the 9 fits considered, that has the smallest number of false positive neurons (i.e., estimated neurons that are noise or duplicates of other estimated neurons). This fit estimated 24 neurons: 21 elements correspond to neurons identified using SCALPEL, one element (element 23 in Figure 8(b)) corresponds to a neuron not identified using SCALPEL, and two elements (elements 22 and 24 in Figure 8(b)) appear to be duplicates of other estimated neurons. In Figure 8(c), we present the fit, chosen from the 9 fits considered, with the highest number of true positive neurons (i.e., neurons that were identified by CNMF-E that appear to be real). This fit estimated 41 neurons: 25 elements correspond to neurons identified using SCALPEL, two elements (elements 27 and 34 in Figure 8(c)) correspond to neurons not identified using SCALPEL, 11 elements appear to be duplicates of other estimated neurons, and three elements (elements 39, 40, and 41 in Figure 8(c)) appear to be noise. So while this pair of tuning parameter values resulted in the identification of most of the neurons, it also resulted in a number of false positives. Some of the estimated neurons in Figure 8(c) are large and diffuse making them difficult to interpret.

6.2 Application to Two-Photon Calcium Imaging Data

We now illustrate SCALPEL on two calcium imaging videos released by the Allen Institute as part of their Allen Brain Observatory. In addition to releasing the data, the Allen Institute also made available the spatial masks for the neurons they identified in each of the videos. Thus we compare the estimated neurons from SCALPEL to those from the Allen Institute analysis. The two-photon videos we consider are those from experiments 496934409 and 502634578. The videos contain 105,698 and 105,710 frames, respectively, of size $512\times 512$ pixels. In their analyses, the Allen Institute down-sampled the number of frames in each video by 8, which we also did for comparability. For these videos, we found using a threshold value slightly smaller than the negative of the 0.1% quantile was effective. Otherwise, default values were used for all of the tuning parameters.

6.2.1 Allen Brain Observatory Experiment 496934409

Using thresholds of 0.250, 0.423, and 0.596, Step 1 of SCALPEL resulted in a preliminary dictionary with 68,630 elements, which came from 11,739 different frames of the video. After refining the dictionary in Step 2, we were left with 544 elements. In the analysis by the Allen Institute, neurons near the boundary of the field of view were eliminated from consideration. Thus we filtered out 259 elements that contained pixels outside of the region considered by the Allen Institute. Of the remaining 285 elements, 32 of these were determined to be dendrites, 131 were small elements not of primary interest, and 10 were duplicates of other neurons found. Thus in the end, we identified the same 87 neurons that the Allen Institute did, in addition to 25 potential neurons not identified by the Allen Institute. In Figure 9(a), we show the neurons jointly identified by SCALPEL and the Allen Institute. In Figure 9(b), we show the potential neurons uniquely identified by SCALPEL, along with evidence supporting them being neurons in Figure 9(c).

6.2.2 Allen Brain Observatory Experiment 502634578

Using thresholds of 0.250, 0.481, and 0.712, Step 1 of SCALPEL resulted in a preliminary dictionary with 84,996 elements, which came from 12,272 different frames of the video. After refining the dictionary in Step 2, we were left with 1297 elements. Once again, we filtered out the 390 elements that contained pixels outside of the region considered by the Allen Institute. Of the remaining 907 elements, 22 of these were determined to be dendrites, 382 were small elements not of primary interest, and 39 were duplicates of other neurons found. Thus in the end, we identified 370 of the 375 neurons that the Allen Institute did, in addition to 94 potential neurons not identified by the Allen Institute. Note that the 5 neurons identified by the Allen Institute, but not SCALPEL, each appear to be combinations of two neurons. SCALPEL did identify the 10 individual neurons of which these 5 Allen Institute neurons were a combination. In Figure 10(a), we show the neurons jointly identified by SCALPEL and the Allen Institute. In Figure 10(b), we show the potential neurons uniquely identified by SCALPEL, along with evidence supporting them being neurons in Figure 10(c).

6.3 Timing Results

All analyses were run on a Macbook Pro with a 2.0 GHz Intel Sandy Bridge Core i7 processor. Running the SCALPEL pipeline on the one-photon data presented in Section 6.1 took 6 minutes for Step 0 and 2 minutes for Steps 1-3. Running CNMF-E on the one-photon data presented in Section 6.1 took 5, 5, and 7 minutes for the analyses presented in Figures 8(a), 8(b), and 8(c), respectively. Running the SCALPEL pipeline on the two-photon data presented in Section 6.2.1 took 12.85 hours for Step 0, 2.33 hours for Step 1, and 0.42 hours for Step 2. Running the SCALPEL pipeline on the two-photon data presented in Section 6.2.2 took 12.50 hours for Step 0, 2.55 hours for Step 1, and 0.43 hours for Step 2.

Further computational gains could be made by parallelizing the implementation of SCALPEL Steps 0 and 1. In the instance of the analysis of the Allen Institute data, this could provide an almost 10-fold decrease in the computational time required. Also, recall that SCALPEL’s most time-intensive step, Step 0, is only ever run a single time for each data set, regardless of whether the user wishes to fit SCALPEL for different tuning parameters.

7 Further Discussion of Step 3

In this section, we elaborate on issues related to solving the sparse group lasso problem (5) in Step 3. The discussion in this section is somewhat technical, and can be skipped by readers only interested in the practical use of SCALPEL. We discuss the solution to (5) when ${K_{f}}=1$ in Section 7.1, our algorithm for solving (5) for any value of ${K_{f}}$ in Section 7.2, the justification for the scaling of ${\bm{A}}^{f}$ in Section 7.3, a result about the tuning parameters $\alpha$ and $\lambda$ that lead to a sparse solution in Section 7.4, and the ability of the group lasso penalty in (5) to zero out unwanted dictionary elements in Section 7.5.

7.1 Single Component Problem

We first consider solving (5) in the setting in which there is a single spatial component (${K_{f}}=1$). While calcium imaging data will not have only a single neuron, this setting provides intuition, and will prove useful when we later solve (5) for ${K_{f}}>1$ in Section 7.2.

Lemma 7.1**.**

The solution to

[TABLE]

is

[TABLE]

where the positive part operator is applied element-wise.

The proof of Lemma 7.1 is in Appendix A.6. We can inspect the solution (7) to gain intuition. Recall that $\tilde{\bm{a}}^{f}_{\cdot,k}\equiv{\bm{a}}^{f}_{\cdot,k}/\|{\bm{a}}^{f}_{\cdot,k}\|_{2}^{2}$, where ${\bm{a}}^{f}_{\cdot,k}$ has binary elements. Therefore, ${\bm{Y}}^{\top}\tilde{\bm{a}}^{f}\in\mathbb{R}^{T}$ is the average fluorescence of pixels in the filtered dictionary element at each of the frames and $\frac{1}{(\tilde{\bm{a}}^{f})^{\top}\tilde{\bm{a}}^{f}}$ equals the number of pixels in the dictionary element. When $\lambda=0$, $\hat{\bm{z}}=\left(\frac{{\bm{Y}}^{\top}\tilde{\bm{a}}^{f}}{(\tilde{\bm{a}}^{f})^{\top}\tilde{\bm{a}}^{f}}\right)_{+}$, which is simply the positive part of the total fluorescence at all pixels in the dictionary element over time. We now consider the impact of $\lambda$ for three different settings of $\alpha$:

• •

$\alpha=1$: In this setting, $\hat{\bm{z}}$ is the positive part of the soft-thresholded total fluorescence. This encourages elements of $\hat{\bm{z}}$ to be exactly zero for frames in which the dictionary element has low fluorescence.

• •

$\alpha=0$: In this setting, $\hat{\bm{z}}$ is found by scaling all elements of the total fluorescence by the same amount. Thus, individual elements of $\hat{\bm{z}}$ are not encouraged to be 0, though $\hat{\bm{z}}={\bm{0}}$ if the dictionary element has a low amount of fluorescence across all frames (i.e., $\|{\bm{Y}}^{\top}\tilde{\bm{a}}^{f}\|_{2}$ small) or $\lambda$ is very large.

• •

$\alpha\in(0,1)$: Both soft-thresholding and soft-scaling are performed, which encourages sparsity of individual elements of $\hat{\bm{z}}$ and the entire vector $\hat{\bm{z}}$, respectively.

In Figure 11, we illustrate the values of $\hat{\bm{z}}$ for the three scenarios described above.

7.2 Algorithm

We now consider how to solve (5) for ${K_{f}}>1$. While generalized gradient descent (Beck and Teboulle, 2009) can be used to solve concurrently for ${\bm{z}}_{1,\cdot},\ldots,{\bm{z}}_{{K_{f}},\cdot}$ in (5), the problem is solved more efficiently by noting that (5) is decomposable into groups of overlapping spatial components.

Let $\mathcal{N}_{1},\ldots,\mathcal{N}_{S}$ denote a partition of the ${K_{f}}$ elements of the filtered dictionary, such that $\mathcal{N}_{s}\cap\mathcal{N}_{s^{\prime}}=\emptyset$ for $s\neq s^{\prime}$, and $\cup_{s=1}^{S}\mathcal{N}_{s}=\{1,\ldots,{K_{f}}\}$. Define the mapping $\mathcal{M}(\mathcal{N}_{s})=\{p\in(1,\ldots,P):(\tilde{{\bm{A}}^{f}}_{p,\mathcal{N}_{s}})^{\top}{\bm{1}}>0\}$. That is, $\mathcal{M}(\mathcal{N}_{s})$ indexes the set of pixels that are active in that subset of neurons.

Lemma 7.2**.**

Suppose that $\mathcal{M}(\mathcal{N}_{s})\cap\mathcal{M}(\mathcal{N}_{s^{\prime}})=\emptyset$ for all $s\neq s^{\prime}$, so that there is no spatial overlap between the sets of filtered dictionary elements $\mathcal{N}_{1},\ldots,\mathcal{N}_{S}$. Then solving (5) gives the same solution as solving

[TABLE]

for $s=1,2,\ldots,S$.

The proof of Lemma 7.2 is in Appendix A.7.

Our approach to solving (8) depends on the size of $\mathcal{N}_{s}$. For $|\mathcal{N}_{s}|=1$, we can simply use the closed-form solution for ${\bm{z}}_{\mathcal{N}_{s},\cdot}$ given by Lemma 7.1. This is advantageous as the calcium imaging data sets that we have analyzed often have some dictionary elements that do not overlap with any others. For $|\mathcal{N}_{s}|>1$, we use generalized gradient descent to solve for the global optimum of (8) (Beck and Teboulle, 2009).

In light of Lemma 7.2, in order to solve (5), we first partition the filtered dictionary elements into $S$ sets, $\mathcal{N}_{1},\ldots,\mathcal{N}_{S}$, such that there is no overlap between the pixels in the $S$ sets, and so that no set can be partitioned further. This can be done quickly, as outlined in Step 1 of Algorithm 1. Then, we solve (8) for $s=1,\ldots,S$. Details are provided in Algorithm 1.

We typically solve (5) for a sequence of exponentially decreasing $\lambda$ values. To improve computational performance, Step 2(b) of Algorithm 1 can be implemented using warm starts, in which ${\bm{Z}}_{\mathcal{N}_{s},\cdot}^{(0)}$ is initialized as the solution for ${\bm{Z}}_{\mathcal{N}_{s},\cdot}$ for the previous value of $\lambda$. Additional details regarding the derivation of the generalized gradient descent algorithm used in Step 2(b) of Algorithm 1 are given in Appendix A.8.

7.3 Scaling of ${\bm{A}}^{f}$

In (5), the $k$th column of the matrix $\tilde{\bm{A}}^{f}$ encodes the spatial mapping of the $k$th filtered dictionary element, after scaling. To obtain $\tilde{\bm{A}}^{f}$, we divide the $k$th column of ${\bm{A}}^{f}$ by $\|{\bm{a}}^{f}_{\cdot,k}\|_{2}^{2}$, the number of pixels in the $k$th filtered dictionary element. This scaling is performed so that the sizes of the dictionary elements do not impact when the components enter the model. That is, we would like $\|{\bm{a}}^{f}_{\cdot,k}\|_{2}^{2}$ to be independent of the largest value of $\lambda$ for which $\hat{\bm{z}}_{k,\cdot}\neq{\bm{0}}$. The following lemma supports this particular scaling of the columns of ${\bm{A}}^{f}$.

Lemma 7.3**.**

Suppose ${\bm{Y}}={\bm{A}}^{f}{\bm{Z}}^{*}$ where the following conditions hold:

1. (i)

${\bm{A}}^{f}\in\mathbb{R}^{P\times{K_{f}}}$* with $({\bm{a}}^{f}_{\cdot,1})^{\top}{\bm{a}}^{f}_{\cdot,2}=0$, $({\bm{a}}^{f}_{\cdot,1})^{\top}{\bm{a}}^{f}_{\cdot,k}=0$ for $k=3,\ldots,{K_{f}}$, and $({\bm{a}}^{f}_{\cdot,2})^{\top}{\bm{a}}^{f}_{\cdot,k}=0$ for $k=3,\ldots,{K_{f}}$ and* 2. (ii)

${\bm{Z}}^{*}\in\mathbb{R}^{{K_{f}}\times T}$* with ${\bm{z}}^{*}_{1,\cdot}={\bm{P}}{\bm{z}}^{*}_{2,\cdot}$ for some $T\times T$ permutation matrix ${\bm{P}}$.*

If we solve (5) for ${\bm{Z}}$ with $\tilde{\bm{A}}^{f}$ such that $\tilde{\bm{a}}^{f}_{\cdot,k}={\bm{a}}^{f}_{\cdot,k}/\|{\bm{a}}^{f}_{\cdot,k}\|_{2}^{2}$, then $\hat{\bm{z}}_{1,\cdot}={\bm{0}}$ if and only if $\hat{\bm{z}}_{2,\cdot}={\bm{0}}$.

The proof of Lemma 7.3 is in Appendix A.9. Lemma 7.3 indicates that two non-overlapping spatial components, possibly of different sizes, whose temporal components are identical up to a permutation, will enter the model at the same value of $\lambda$. In Figure 12, we provide empirical evidence for the chosen scaling of ${\bm{A}}^{f}$.

7.4 Sparsity of the Solution

We now consider the range of $\lambda$ for which the solution to (5) is completely sparse (i.e., $\hat{\bm{Z}}={\bm{0}}$) for a fixed value of $\alpha$.

Lemma 7.4**.**

For any $\alpha\in[0,1]$, the solution to (5) is completely sparse if and only if

[TABLE]

for $k=1,\ldots,{K_{f}}$.

Unfortunately, when $\alpha\in(0,1)$, $\lambda$ is on both sides of the inequality in (11). Though we can solve for $\lambda$ in (11) using a root finder when $\alpha\in(0,1)$, the following corollary provides a simple alternative.

Corollary 7.5**.**

For any $\alpha\in(0,1)$, if

[TABLE]

then the solution to (5) is completely sparse.

The condition in Corollary 7.5 is sufficient, but not necessary. Proofs of Lemma 7.4 and Corollary 7.5 can be found in Appendices A.10 and A.11, respectively. An illustration of Lemma 7.4 and Corollary 7.5 is provided in Figure 13. In Section 5.3, we discuss how the results from Lemma 7.4 and Corollary 7.5 can assist in selecting the value of the tuning parameter $\lambda$.

7.5 Zeroing Out of Double Neurons

Some of the elements in the preliminary dictionary obtained in Step 1 may be double neurons, i.e., elements that are a combination of two separate neurons. This occurs when two neighboring neurons are active during the same frame. In Step 2 of SCALPEL, these double neurons are unlikely to cluster with elements representing either of the individual neurons they combine, and thus there may be double neurons that remain in the filtered set of dictionary elements, ${\bm{A}}^{f}$, used in Step 3 of SCALPEL. Fortunately, as detailed in the following lemma, the group lasso penalty in (5) filters out these double neurons by estimating their temporal components to be the zero vector.

Lemma 7.6**.**

Suppose that the following conditions hold on ${\bm{A}}^{f}\in\mathbb{R}^{P\times{K_{f}}}$:

1. (i)

${\bm{a}}^{f}_{\cdot,3}={\bm{a}}^{f}_{\cdot,1}+{\bm{a}}^{f}_{\cdot,2}$, 2. (ii)

$({\bm{a}}^{f}_{\cdot,1})^{\top}{\bm{a}}^{f}_{\cdot,2}=0$, 3. (iii)

$({\bm{a}}^{f}_{\cdot,1})^{\top}{\bm{a}}^{f}_{\cdot,k}=0$* for $k=4,\ldots,{K_{f}}$, and* 4. (iv)

$({\bm{a}}^{f}_{\cdot,2})^{\top}{\bm{a}}^{f}_{\cdot,k}=0$* for $k=4,\ldots,{K_{f}}$.*

Then, define $\tilde{\bm{a}}^{f}_{\cdot,k}\equiv{\bm{a}}^{f}_{\cdot,k}/\|{\bm{a}}^{f}_{\cdot,k}\|_{2}$, and consider solving (5) for ${\bm{Z}}$ with $\alpha<1$. Then, $\hat{\bm{z}}_{3,\cdot}={\bm{0}}$.

The proof of Lemma 7.6 is in Appendix A.12. Note that Lemma 7.6 assumes that the individual elements for the neighboring neurons, ${\bm{a}}^{f}_{\cdot,1}$ and ${\bm{a}}^{f}_{\cdot,2}$, do not overlap at all. The group lasso penalty can also be effective at zeroing out double neurons resulting from overlapping neurons, though this depends on the amount of overlap, among other factors.

8 Discussion

We have presented SCALPEL, a method for simultaneously identifying neurons from calcium imaging data and estimating their intracellular calcium concentrations. SCALPEL takes a dictionary learning approach. We segment the frames of the calcium imaging video to construct a large preliminary dictionary of potential neurons, which is then refined through the use of clustering using a novel dissimilarity metric that leverages both spatial and temporal information. The calcium concentrations of the elements of the refined dictionary are then estimated by solving a sparse group lasso with a non-negativity constraint.

Future work could consider alternative ways of deriving a preliminary dictionary in Step 1. Currently, we perform image segmentation via thresholding with multiple quantiles. This approach assumes that active neurons will have brightness, relative to their baseline fluorescence levels, that is within the range of our image segmentation threshold values. In practice, there is evidence that some neurons have comparatively lower fluorescence following spiking, which presents a challenge for optimal identification. Though SCALPEL performed well on the one-photon calcium imaging video we considered in Section 6.1, other one-photon videos may have more severe background effects. If this is the case, it may be desirable to incorporate more sophisticated modeling of the background noise, like that employed in Zhou et al. (2016). Additionally, in future work, we could modify Step 3 of SCALPEL to make use of a more refined model for neuron spiking, as in Friedrich et al. (2015); Vogelstein et al. (2010).

Our SCALPEL proposal is implemented in the R package scalpel, which is available on CRAN. A vignette illustrating how to use the package and code to reproduce all results presented in this paper are available at ajpete.com/software.

Acknowledgements

We thank Ilana Witten and Malavika Murugan at the Princeton Neuroscience Institute for providing guidance and access to the calcium imaging data set analyzed in Section 6.1. We also thank Michael Buice at the Allen Institute for Brain Science for his helpful comments and providing access to the calcium imaging data sets analyzed in Section 6.2. We also thank Liam Paninski and Pengcheng Zhou for helpful suggestions and for providing software implementing CNMF-E.

Appendix A Supplementary Material

A.1 Data Pre-Processing

To begin, we perform three pre-processing steps on the raw data. These are briefly described in Step 0 of Section 4. First, we smooth the raw $P\times T$ data matrix spatially and temporally using a Gaussian kernel smoother with a bandwidth of 1 pixel (Hastie et al., 2009). Second, we adjust for any bleaching effect over time. Specifically, we fit a smoothing spline with 10 degrees of freedom to the median fluorescence for each frame over time, and subtract the frame-specific smoothed median from the corresponding frame. An example of a smoothing spline fit for the example calcium imaging video is shown in Figure 14.

Finally, we apply a slight variation of the often-used $\Delta f/f$ transformation (Grewe et al., 2010; Grienberger and Konnerth, 2012; Ahrens et al., 2013). For the $i$th pixel in the $j$th frame, the standardized fluorescence is equal to

[TABLE]

where $y^{0}_{i,j}$ is the fluorescence (after smoothing and bleaching correction) of the $i$th pixel in the $j$th frame. This differs from the typical $\Delta f/f$ transformation in that (1) we standardize using the median across image frames instead of the mean; and (2) we add a small number to the denominator. This adjustment in the denominator prevents small fluctuations in the amount of fluorescence at pixels with very little overall fluorescence from resulting in extremely high standardized fluorescences. In Figure 15, we show the resulting images after each stage of pre-processing for a sample frame.

A.2 Alternative Spatial Dissimilarity Metrics

We consider three alternatives to the spatial dissimilarity proposed in (3):

[TABLE]

We refer to the dissimilarities in (12) as the union, min, and max dissimilarities, respectively. Note that all of the dissimilarities in (3) and (12) take on values in $[0,1]$.

To see why we prefer the dissimilarity measure in (3) over those in (12), consider Figure 16(d), in which one component is located entirely within the other, and the two are of substantially different sizes. Min dissimilarity assigns a value of 0 in this case, despite the fact that one component is much smaller than the other. Union dissimilarity and max dissimilarity yield a very large value of 0.62, despite the fact that there is substantial overlap between the two components. In contrast, the dissimilarity in (3) yields a modest but non-zero dissimilarity, which seems reasonable given that the two components in Figure 16(d) are similar but not identical. More generally, the dissimilarity in (3) tends to yield reasonable values across a range of settings (Figures 16(a)-(c)).

Often Euclidean distance, $\|{\bm{a}}^{0}_{\cdot,i}-{\bm{a}}^{0}_{\cdot,j}\|_{2}$, is used as the distance metric in clustering. Here, the squared Euclidean distance equals the number of differing pixels between components. This is problematic as a distance metric for our purposes — it treats $a^{0}_{n,i}=a^{0}_{n,j}=0$ the same as $a^{0}_{n,i}=a^{0}_{n,j}=1$, though the latter is much more informative in terms of the components’ similarity. For example, we see that the pairs of components in Figures 16(b) and (c) have similar Euclidean distances, though the components in Figure 16(b) are non-overlapping and could be arbitrarily far apart without affecting the Euclidean distance.

A.3 Example of a Cluster in Step 2

To see how the preliminary dictionary element in a cluster relates to the other preliminary dictionary elements assigned to that cluster, we give an example in Figure 17.

A.4 Selecting $\lambda$ for (5)

To choose $\lambda$ for (5) via a validation set approach, we perform the following steps:

Obtain $\tilde{\bm{A}}^{f}\in\mathbb{R}^{P\times{K_{f}}}$ by dividing the $k$th column of ${\bm{A}}^{f}$ by $\|{\bm{a}}^{f}_{\cdot,k}\|_{2}^{2}$. 2. 2.

Construct a training set $\mathcal{T}$ by sampling 60% of the pixels in each overlapping group of neurons. That is, we sample 60% of the elements in $\mathcal{M}(\mathcal{N}_{1}),\mathcal{M}(\mathcal{N}_{2}),\ldots,\mathcal{M}(\mathcal{N}_{S})$. Assign the remaining pixels to the validation set, $\mathcal{V}=\{v\in(1,\ldots,P):v\notin\mathcal{T}\}$. 3. 3.

Using Algorithm 1, solve (5) on the training set of pixels for a decreasing sequence of 20 $\lambda$ values, $\lambda_{1},\ldots,\lambda_{20}$:

[TABLE] 4. 4.

For each $\lambda_{i}$, calculate the validation error,

[TABLE]

where $\left[{\bm{Y}}_{\mathcal{V},\cdot}^{B}\right]_{j,k}=\begin{cases}\left[{\bm{Y}}_{\mathcal{V},\cdot}\right]_{j,k}&\text{if }\left[{\bm{Y}}_{\mathcal{V},\cdot}\right]_{j,k}>-\text{quantile}_{0.1\%}({\bm{Y}})\\ 0&\text{otherwise}\end{cases}$. We use a thresholded version of ${\bm{Y}}$ when calculating the validation error, as we only wish to have small reconstruction error on the brightest parts of the video. Select the optimal value of $\lambda$ as

[TABLE]

this is the largest value of $\lambda$ that results in a validation error within 5% of the minimum validation error achieved by any value of $\lambda$ considered. 5. 5.

Solve (5) on all pixels:

[TABLE]

where we have scaled the tuning parameter by the percent of pixels in the training set, in order to account for the fact that the sum of squared errors in the loss function is not scaled by the number of pixels.

This process can be done separately for each group of overlapping neurons $\mathcal{N}_{1},\ldots,\mathcal{N}_{S}$ in order to select a different value of $\lambda$ for each group, or for all groups at once to select a single value of $\lambda$. By following steps similar to those described above, $\lambda$ can alternatively be selected via cross-validation.

A.5 Filtering Dictionary Elements Prior to Fitting the Sparse Group Lasso

Optionally, at the beginning of Step 3, the elements in the refined dictionary from Step 2 may be filtered prior to fitting the sparse group lasso in (5). One way to filter the dictionary elements is on the basis of the number of members in the clusters. That is, we can discard any dictionary elements representing clusters containing fewer than some minimum number of members. This number should be chosen based on the goals of the analysis. If we retain all clusters, regardless of size, then we may include some non-neuronal dictionary elements in the sparse group lasso problem. In contrast, if we discard small clusters (e.g., those with fewer than 5 members), then we may erroneously filter out some true neurons.

As an alternative to filtering dictionary elements on the basis of cluster size, we can instead manually inspect each dictionary element, by examining the frames of the video from which it was derived. In Figure 18, we compare the results obtained by manually filtering each dictionary element, versus simply filtering dictionary elements on the basis of cluster size, on the example video considered in Sections 4 and 6.1. We find that the two approaches give similar results.

A.6 Proof of Lemma 7.1

We first prove a result that we will use later.

Lemma A.1**.**

The solution to $\underset{{\bm{\beta}}\geq{\bm{0}}}{\mathrm{minimize}}\quad\frac{1}{2}\left\|{\bm{y}}-{\bm{\beta}}\right\|_{2}^{2}+\lambda\left\|{\bm{\beta}}\right\|_{2}$ is $\hat{\bm{\beta}}=\left(1-\frac{\lambda}{\left\|\left({\bm{y}}\right)_{+}\right\|_{2}}\right)_{+}\left({\bm{y}}\right)_{+}$.

Proof.

Let $\hat{\bm{\beta}}=\underset{{\bm{\beta}}\geq{\bm{0}}}{\mathrm{argmin}}\;\frac{1}{2}\left\|{\bm{y}}-{\bm{\beta}}\right\|_{2}^{2}+\lambda\left\|{\bm{\beta}}\right\|_{2}$ and $\mathcal{C}=\{i:y_{i}\geq 0\}$. First, we show $\hat{\bm{\beta}}_{-\mathcal{C}}={\bm{0}}$. In anticipation of contradiction, assume there exists $j$ such that $j\notin\mathcal{C}$ and $\hat{\beta}_{j}>0$. Define $\tilde{\bm{\beta}}$ as $\tilde{\beta}_{i}=\begin{cases}\hat{\beta}_{i}&\text{if }i\neq j\\ 0&\text{if }i=j\end{cases}.$ Then

[TABLE]

This is a contradiction, so we conclude that $\hat{\beta}_{i}=0$ for all $i\notin\mathcal{C}$. It remains to solve

[TABLE]

By a result in Section 3.1 of Simon et al. (2013), the solution to (13) without the non-negativity constraint on ${\bm{\beta}}_{\mathcal{C}}$ is $\left(1-\frac{\lambda}{\left\|{\bm{y}}_{\mathcal{C}}\right\|_{2}}\right)_{+}{\bm{y}}_{\mathcal{C}}$, which has all non-negative elements. Therefore, it is also the solution to (13). ∎

We now proceed to prove Lemma 7.1.

Proof.

Our goal is to solve

[TABLE]

Note that solving (14) is equivalent to solving (6), as $\|{\bm{z}}\|_{1}={\bm{1}}^{\top}{\bm{z}}$ when ${\bm{z}}\geq{\bm{0}}$. By algebraic manipulation, we can show that

[TABLE]

where $C$ is a constant that does not depend on ${\bm{z}}$. Therefore, the solution to (14) is the same as the solution to

[TABLE]

We solve (15) by applying Lemma A.1. ∎

A.7 Proof of Lemma 7.2

Proof.

The result follows simply from observing that

[TABLE]

The last equality follows from the condition of the lemma, which guarantees that $\tilde{{\bm{A}}^{f}}_{\mathcal{M}(\mathcal{N}_{s}),\mathcal{N}_{s^{\prime}}}=0$ for all $s\neq s^{\prime}$. ∎

A.8 Details of Step 2(b) of Algorithm 1

Note that minimizing the objective in (9) subject to ${\bm{Z}}_{\mathcal{N}_{s},\cdot}\geq{\bm{0}}$ is equivalent to minimizing

[TABLE]

subject to ${\bm{Z}}_{\mathcal{N}_{s},\cdot}\geq{\bm{0}}$, since ${\bm{1}}^{\top}{\bm{Z}}_{\mathcal{N}_{s},\cdot}{\bm{1}}=\sum_{n\in\mathcal{N}_{s}}\left\|{\bm{z}}_{n,\cdot}\right\|_{1}$ when ${\bm{Z}}_{\mathcal{N}_{s},\cdot}\geq{\bm{0}}$.

Let $f\left({\bm{Z}}_{\mathcal{N}_{s},\cdot}\right)=\frac{1}{2}\left\|{\bm{Y}}_{\mathcal{M}(\mathcal{N}_{s}),\cdot}-\tilde{\bm{A}}^{f}_{\mathcal{M}(\mathcal{N}_{s}),\mathcal{N}_{s}}{\bm{Z}}_{\mathcal{N}_{s},\cdot}\right\|_{F}^{2}+\lambda\alpha{\bm{1}}^{\top}{\bm{Z}}_{\mathcal{N}_{s},\cdot}{\bm{1}}$, which is the differentiable part of (16), and let $g\left({\bm{Z}}_{\mathcal{N}_{s},\cdot}\right)=\lambda(1-\alpha)\sum_{n\in\mathcal{N}_{s}}\left\|{\bm{z}}_{n,\cdot}\right\|_{2}$, the non-differentiable part.

Generalized gradient descent (Beck and Teboulle, 2009; Parikh and Boyd, 2014) is a majorization-minimization scheme. First, we find a quadratic approximation to $f\left({\bm{Z}}_{\mathcal{N}_{s},\cdot}\right)$ centered at our previous estimate for ${\bm{Z}}_{\mathcal{N}_{s},\cdot}$, ${\bm{Z}}_{\mathcal{N}_{s},\cdot}^{0}$, that majorizes $f\left({\bm{Z}}_{\mathcal{N}_{s},\cdot}\right)$. That is,

[TABLE]

where $t$ is the step size such that $\nabla^{2}f(\cdot)\preceq\frac{1}{t}{\bm{I}}$. After completing the square, we can see that minimizing the quadratic approximation to $f\left({\bm{Z}}_{\mathcal{N}_{s},\cdot}\right)$ gives the same solution as solving

[TABLE]

Thus we perform this minimization with $g\left({\bm{Z}}_{\mathcal{N}_{s},\cdot}\right)$ added to the objective function, which gives the proximal problem

[TABLE]

where $\tilde{\bm{Y}}_{\mathcal{N}_{s},\cdot}={\bm{Z}}_{\mathcal{N}_{s},\cdot}^{0}-t\left(-(\tilde{\bm{A}}^{f}_{\mathcal{M}(\mathcal{N}_{s}),\mathcal{N}_{s}})^{\top}\left({\bm{Y}}_{\mathcal{M}(\mathcal{N}_{s}),\cdot}-\tilde{\bm{A}}^{f}_{\mathcal{M}(\mathcal{N}_{s}),\mathcal{N}_{s}}{\bm{Z}}_{\mathcal{N}_{s},\cdot}^{0}\right)+\lambda\alpha{\bm{1}}{\bm{1}}^{\top}\right)$. The minimization in (17) is separable in ${\bm{z}}_{n,\cdot}$, so for $n\in\mathcal{N}_{s}$, we solve

[TABLE]

By Lemma A.1 in Appendix A.6, the solution to (18) is $\hat{\bm{z}}_{n,\cdot}=\left(1-\frac{\lambda(1-\alpha)t}{\left\|\left(\tilde{\bm{y}}_{n,\cdot}\right)_{+}\right\|_{2}}\right)_{+}\left(\tilde{\bm{y}}_{n,\cdot}\right)_{+}$.

It only remains to derive a suitable step size $t$ so that $\nabla^{2}f(\cdot)=(\tilde{\bm{A}}^{f}_{\mathcal{M}(\mathcal{N}_{s}),\mathcal{N}_{s}})^{\top}\tilde{\bm{A}}^{f}_{\mathcal{M}(\mathcal{N}_{s}),\mathcal{N}_{s}}\preceq\frac{1}{t}{\bm{I}}$. A sufficient condition for $\frac{1}{t}{\bm{I}}-(\tilde{\bm{A}}^{f}_{\mathcal{M}(\mathcal{N}_{s}),\mathcal{N}_{s}})^{\top}\tilde{\bm{A}}^{f}_{\mathcal{M}(\mathcal{N}_{s}),\mathcal{N}_{s}}$ to be positive semi-definite is that $\frac{1}{t}{\bm{I}}-(\tilde{\bm{A}}^{f}_{\mathcal{M}(\mathcal{N}_{s}),\mathcal{N}_{s}})^{\top}\tilde{\bm{A}}^{f}_{\mathcal{M}(\mathcal{N}_{s}),\mathcal{N}_{s}}$ be diagonally dominant. That is,

[TABLE]

for all $n\in\mathcal{N}_{s}$. Thus we choose $t=(\max_{n\in\mathcal{N}_{s}}\sum_{j\in\mathcal{N}_{s}}(\tilde{\bm{a}}^{f}_{\mathcal{M}(\mathcal{N}_{s}),j})^{\top}\tilde{\bm{a}}^{f}_{\mathcal{M}(\mathcal{N}_{s}),n})^{-1}$.

A.9 Proof of Lemma 7.3

Proof.

Since ${\bm{a}}^{f}_{\cdot,1}$ does not overlap any other spatial components (i.e., $({\bm{a}}^{f}_{\cdot,1})^{\top}{\bm{a}}^{f}_{\cdot,k}=0$ for $k=2,\ldots,{K_{f}}$), we know by the results in Lemmas 7.1 and 7.2 that

[TABLE]

Note that ${\bm{Y}}^{\top}\tilde{\bm{a}}^{f}_{\cdot,1}=({\bm{Z}}^{*})^{\top}({\bm{A}}^{f})^{\top}{\bm{a}}^{f}_{\cdot,1}/\|{\bm{a}}^{f}_{\cdot,1}\|_{2}^{2}$, so ${\bm{Y}}^{\top}\tilde{\bm{a}}^{f}_{\cdot,1}={\bm{z}}_{1,\cdot}^{*}$ since $({\bm{a}}^{f}_{\cdot,1})^{\top}{\bm{a}}^{f}_{\cdot,k}=0$ for $k=2,\ldots,{K_{f}}$. Thus $\hat{\bm{z}}_{1,\cdot}={\bm{0}}$ if and only if $\lambda(1-\alpha)\geq\|({\bm{z}}_{1,\cdot}^{*}-\lambda\alpha{\bm{1}})_{+}\|_{2}$. Similarly, $\hat{\bm{z}}_{2,\cdot}={\bm{0}}$ if and only if $\lambda(1-\alpha)\geq\|({\bm{z}}_{2,\cdot}^{*}-\lambda\alpha{\bm{1}})_{+}\|_{2}$. Thus, it remains to show that $\|({\bm{z}}_{1,\cdot}^{*}-\lambda\alpha{\bm{1}})_{+}\|_{2}=\|({\bm{z}}_{2,\cdot}^{*}-\lambda\alpha{\bm{1}})_{+}\|_{2}$. Using the properties of permutation matrices that ${\bm{1}}={\bm{P}}{\bm{1}}$ and ${\bm{P}}^{\top}{\bm{P}}={\bm{I}}$, we see

[TABLE]

and therefore, $\hat{\bm{z}}_{1,\cdot}={\bm{0}}$ if and only if $\hat{\bm{z}}_{2,\cdot}={\bm{0}}$ . ∎

A.10 Proof of Lemma 7.4

Proof.

Recall that solving (5) gives the same solution as solving (8). Thus we focus on deriving a condition on $\lambda$ that guarantees that $\hat{\bm{Z}}_{\mathcal{N}_{s},\cdot}$, the solution to (8), equals zero for $s=1,\ldots,S$. If $|\mathcal{N}_{s}|=1$, we see from (10) that $\hat{\bm{z}}_{\mathcal{N}_{s},\cdot}={\bm{0}}$ if and only if

[TABLE]

Recall that if $|\mathcal{N}_{s}|>1$, we iteratively solve for $\hat{\bm{Z}}_{\mathcal{N}_{s},\cdot}$ using Step 2(b) of Algorithm 1. We initialize at the sparse solution ${\bm{Z}}_{\mathcal{N}_{s},\cdot}^{(0)}={\bm{0}}$ and thus for $n\in\mathcal{N}_{s}$

[TABLE]

where $\tilde{\bm{Y}}_{\mathcal{N}_{s},\cdot}=t(\tilde{\bm{A}}^{f}_{\mathcal{M}(\mathcal{N}_{s}),\mathcal{N}_{s}})^{\top}{\bm{Y}}_{\mathcal{M}(\mathcal{N}_{s}),\cdot}-t\lambda\alpha{\bm{1}}{\bm{1}}^{\top}$. We will have $\hat{\bm{Z}}_{\mathcal{N}_{s},\cdot}={\bm{0}}$ if ${\bm{z}}_{n,\cdot}^{(1)}={\bm{0}}$ for all $n\in\mathcal{N}_{s}$. Note that ${\bm{z}}^{(1)}_{n,\cdot}={\bm{0}}$ if

[TABLE]

By algebraic manipulation, the sparsity conditions given in (19) and (20) can be shown to be equivalent to the condition given in Lemma 7.4. Alternatively, this lemma’s result also follows from inspection of the optimality condition for (5). ∎

A.11 Proof of Corollary 7.5

Proof.

The sufficient condition given in Corollary 7.5 follows from noting that (20) is satisfied if $\lambda(1-\alpha)\geq\left\|\left(\left[(\tilde{\bm{A}}^{f})^{\top}{\bm{Y}}\right]_{k,\cdot}\right)_{+}\right\|_{2}$ or if $\lambda\alpha\geq\left(\left[(\tilde{\bm{A}}^{f})^{\top}{\bm{Y}}\right]_{k,l}\right)_{+}$ for $l=1,\ldots,T$. Thus, when at least one of these two conditions is satisfied for all $k=1,\ldots,{K_{f}}$, then the solution to (5) will be sparse. ∎

A.12 Proof of Lemma 7.6

Proof.

Let $\hat{\bm{Z}}$ be the solution to (5). In anticipation of contradiction, assume there exists $j\in\{1,\ldots,T\}$ such that $\hat{\bm{z}}_{3,j}>0$. Define $\tilde{\bm{Z}}$ as $\tilde{\bm{z}}_{1,\cdot}=\hat{\bm{z}}_{1,\cdot}+\left(\frac{{\bm{1}}^{\top}{\bm{a}}^{f}_{\cdot,1}}{{\bm{1}}^{\top}({\bm{a}}^{f}_{\cdot,1}+{\bm{a}}^{f}_{\cdot,2})}\right)\hat{\bm{z}}_{3,\cdot}$, $\tilde{\bm{z}}_{2,\cdot}=\hat{\bm{z}}_{2,\cdot}+\left(\frac{{\bm{1}}^{\top}{\bm{a}}^{f}_{\cdot,2}}{{\bm{1}}^{\top}({\bm{a}}^{f}_{\cdot,1}+{\bm{a}}^{f}_{\cdot,2})}\right)\hat{\bm{z}}_{3,\cdot}$, $\tilde{\bm{z}}_{3,\cdot}={\bm{0}}$, and $\tilde{\bm{z}}_{k,\cdot}=\hat{\bm{z}}_{k,\cdot}$ for $k=4,\ldots,K^{f}$. Let $\text{obj}({\bm{Z}})$ be the value of the objective function of (5) at ${\bm{Z}}$ for some fixed $\lambda$ and $\alpha$. We have

[TABLE]

This is a contradiction, so we conclude $\hat{\bm{z}}_{3,\cdot}={\bm{0}}$. ∎