{'fraction_non_alphanumeric': 0.037846970031894404, 'fraction_numerical': 0.019326112356712034, 'mean_word_length': 4.413247863247864, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
SCALPEL: Extracting Neurons from Calcium Imaging Data

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.


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, Dombeck et al., 2007, Huber et al., 2012, Prevedel et al., 2014. 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, voltagegated 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 [Chen et al., 2013, Looger and Griesbeck, 2012, Rochefort et al., 2008. 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 Figure 1: In (a), we display sample frames from the raw calcium imaging video described in the text in Section 4, and analyzed in greater detail in Section 6.1. We wish to construct a spatial map of the neurons, like that shown in (b), and to estimate the calcium trace for each neuron over time, as shown in (c).
A typical calcium imaging video consists of a 500 × 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.

Notation
Let P denote the total number of pixels per image frame, and T the number of frames of the video. We define Y to be a P × T matrix for which the (i, j)th element, y i,j , contains the fluorescence of the ith pixel in the jth frame. We let y i,· = y i,1 y i,2 . . . y i,T represent the fluorescence of the ith pixel at each of the T frames. We let y ·,j = y 1,j y 2,j . . . y P,j represent the fluorescence of all P pixels during the jth 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 Y into a matrix of spatial components, A ∈ R P ×K , and a matrix of temporal components, where K is the total number of estimated neurons. Note that a ·,k specifies which of the P pixels of the image frame are mapped to the kth neuron, and z k,· quantifies the calcium concentration for the kth 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.

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 [Mellen andTuong, 2009, Ozden et al., 2008]. 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 A and Z in (1) by solving In (2), the term a b z b is a rank-one correction for background noise: z b ∈ 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 a b ∈ R P is a spatial representation of the background noise. Elementwise positivity constraints are imposed on A, Z, and a b in (1). While (2) can handle overlapping neurons, there are no constraints on the sparsity of A or Z, or locality constraints on 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 A are sparse and the spatial components 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 Z, and combining a sparsity constraint with intermediate image filtering when estimating the spatial components 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 andHamprecht [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 A and 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 A. We then use our estimate of A in order to obtain an accurate estimate of the temporal components Z, while simultaneously selecting the final set of neurons in 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.

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: 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, Y refers to the calcium imaging data after these three pre-processing steps have been performed.
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 A 0 ∈ R P ×K 0 with a 0 j,k = 1 if the jth pixel is contained in the kth preliminary dictionary element, and a 0 j,k = 0 otherwise. This is discussed further in Section 4.1.
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 A ∈ R P ×K : a j,k = 1 if the jth pixel is contained in the kth dictionary element, and a j,k = 0 otherwise. More details are provided in Section 4.2.
Step 3. Application of the Spatial Component Dictionary: We filter out dictionary elements corresponding to clusters with few members, resulting in A f ∈ R P ×K f , which contains a subset of the columns of A.
We then estimate the temporal components Z corresponding to the filtered elements of the dictionary by solving a sparse group lasso problem with a non-negativity constraint. The kth row ofẐ is the estimated calcium trace corresponding to the kth filtered dictionary element. For some values of the tuning parameters, entire rows ofẐ 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 × 226 pixels and 3000 frames. 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.

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.
1. 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).

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. 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. Figure 2: A summary of the SCALPEL procedure, along with the results of applying each step to an example data set with 205 × 226 pixels and 3000 frames, described in the text in Section 4, and analyzed in greater detail in Section 6.1. Figure 3: In (a), we display a single frame of the example calcium imaging video after performing the preprocessing described in Appendix A.1. In (b), we show the binary image that results after thresholding using the negative of the 0.1% quantile of the video's elements. In (c), we display the seven connected components from the image in (b) that contain at least 25 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 Y . Therefore, we threshold each frame using the negative of the minimum value of Y . We also repeat this procedure using a threshold equal to the negative of the 0.1% quantile of 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 A 0 ∈ R P ×K 0 : the kth column of A 0 is a vector of 1's and 0's, indicating whether each pixel is contained in the kth preliminary dictionary element.

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: 1. 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.

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 A, which will be used in Step 3, discussed in Section 4.3.

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 = (a 0 ·,i ) a 0 ·,j , the number of pixels shared between the ith and jth dictionary elements. When i = j, p i,i is simply the number of pixels in the ith dictionary element. We then define the spatial dissimilarity between the ith and jth dictionary elements to be Thus, d s i,j = 1 if and only if the ith and jth 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 Y B , a thresholded version of the pre-processed data matrix Y (obtained in Step 0), with elements of the form 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 Y above. The temporal dissimilarity between the ith and jth dictionary elements is defined as (Note that the elements of Y B a 0 ·,i ∈ R T represent the thresholded fluorescence of each time frame, summed over all pixels in the ith preliminary dictionary element.) We threshold Y before computing this dissimilarity, because (1) we are interested in the extent to which there is agreement between the peak fluorescences of the ith and jth preliminary dictionary elements; and (2) the sparsity induced by thresholding is computationally advantageous.
Finally, the overall dissimilarity is where ω ∈ [0, 1] controls the relative weightings of the spatial and temporal dissimilarities. We use ω = 0.2 to obtain the results shown throughout this paper. In Figure 4, we illustrate pairs of preliminary dictionary elements with various dissimilarities for ω = 0.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, A ∈ R P ×K , as follows: a j,k = 1 if the jth pixel is contained in the kth 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.

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 Figure 4: Each column displays two pairs of preliminary dictionary elements with overall dissimilarities, as defined in (4), of 0.05, 0.1, 0.15, 0.2, and 0.25. For each preliminary dictionary element, the average thresholded fluorescence over time and the (zoomed-in) spatial map are shown. These results are based on the example calcium imaging video.
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, A f ∈ R P ×K f , which contains the retained columns of 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, where α ∈ [0, 1] and λ > 0 are tuning parameters, andÃ f is defined asã f ·,k = a f ·,k / a f ·,k 2 2 . This scaling of 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 Z [Yuan and Lin, 2006], which leads to selection of dictionary components: for λ(1 − α) 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 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.

Tuning Parameter Selection
SCALPEL involves a number of tuning parameters: • Step 1: Quantile thresholds for image segmentation

•
Step 2: Cut-point for dendrogram Figure 5: In (a), we display the dendrogram that results from applying prototype clustering to the example calcium imaging data set. Three different cut-points are indicated: 0.05 ( ), 0.18 ( ), and 0.4 ( ). In (b), we display the number of clusters that result from these three cut-points. In (c)-(e), we show the refined dictionary elements that result from using these three cut-points. For simplicity, we only display dictionary elements corresponding to clusters with at least 5 members.
• Step 3: λ and α 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 λ 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.

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 Y , the negative of the 0.1% quantile of 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).

Tuning Parameters for Step 2
In Step 2, we must choose a height h ∈ [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 ω in (4), as this guarantees that the dictionary elements within each cluster will have spatial overlap. Note that we do not consider ω to be a tuning parameter, as we recommend keeping it fixed at ω = 0.2.

Tuning Parameters for Step 3
In Step 3, we also must choose values of λ and α in (5). We have had empirical success using a fixed value of α near 1, as we expect each neuron to spike a small number of times. We used α = 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 λ. 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 λ that should be considered in our cross-validation scheme. Contrary to standard practice, we use a thresholded version of 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 Y provides an efficient way to choose λ. We can think of the elements of 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 Y . The closed-form solution for non-overlapping neurons given in (7) implies that a frame ofẑ is zeroed out if its average fluorescence is less than λα. Thus, we choose λ = −quantile 0.1% (Y ) α . Alternatively, we can choose the value of λ by visually inspecting the calcium traces (that is, the rows of Z), as shown for instance in Figure 1(c), at each value of λ.
Any of these approaches can be used to choose the value of λ 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 λ for each group of overlapping neurons. This corresponds to a straightforward modification of the optimization problem (5).

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.

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 × 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 α = 0.9 and λ = 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 σ, which relate to the spatial filtering, and the minimum local correlation c min and the minimum peak-to-noise ratio α 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 σ = 3, c min = 0.85, and α 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 min and α min . We fit CNMF-E for all combinations of c min = 0.5, 0.6, 0.7 and α 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 Figure 6: In (a), we plot the spatial maps for the 29 elements of the final dictionary for the calcium imaging video considered in Section 6.1. In (b), we plot their estimated intracellular calcium concentrations corresponding to λ chosen via a validation set approach. In (c), we compare the outlines of the 29 dictionary elements from (a) to a heat map of the pixel-wise variance of the calcium imaging video. That is, we plot the variance of each pixel over the 3000 frames, with whiter points indicating higher variance.   [Zhou et al., 2016] to the calcium imaging video considered in Section 6.1 for (a) the default parameters c min = 0.85 and α min = 10, (b) the parameters c min = 0.6 and α min = 7, and (c) the parameters c min = 0.5 and α min = 3. The variation in darkness of the neurons estimated by CNMF-E is due to the fact that they take on continuous values, compared to the binary masks produced by SCALPEL. In each plot, the 32 true neurons identified by SCALPEL are outlined in gray. 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.

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 × 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.

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). Figure 9: We present the results for the calcium imaging video analyzed in Section 6.2.1. In (a), we plot the outlines of the neurons identified by the Allen Institute in blue, along with the outlines of the corresponding SCALPEL neurons in orange. In (b), we plot the 25 potential neurons uniquely identified by SCALPEL in color, along with the SCALPEL neurons also identified by the Allen Institute in gray. In (c), we provide evidence for 4 of the 25 unique neurons. Similar plots for all of the potential neurons uniquely identified by SCALPEL are available at www.ajpete.com/software.

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).

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. Figure 10: We present the results for the calcium imaging video analyzed in Section 6.2.2. In (a), we plot the outlines of the neurons identified by the Allen Institute in blue, along with the outlines of the corresponding SCALPEL neurons in orange. Those shown in green are the Allen Institute neurons that appear to actually be a combination of 2 neurons. In (b), we plot the 94 potential neurons uniquely identified by SCALPEL in color, along with the SCALPEL neurons also identified by the Allen Institute in gray. In (c), we provide evidence for 4 of the 94 unique neurons. Similar plots for all of the potential neurons uniquely identified by SCALPEL are available at www.ajpete.com/software.

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 A f in Section 7.3, a result about the tuning parameters α and λ 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.

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 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 , where a f ·,k has binary elements. Therefore, Y ã f ∈ R T is the average fluorescence of pixels in the filtered dictionary element at each of the frames and 1 (ã f ) ã f equals the number of pixels in the dictionary element. When λ = 0,ẑ = , 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 λ for three different settings of α: Figure 11: For a single dictionary element in the example video, we plot the solutionẑ, as given in (7), for a range of λ when (a) only a lasso penalty is used (α = 1), (b) a mixture of penalties is used (α = 0.5), and (c) only a group lasso penalty is used (α = 0).
• α = 1: In this setting,ẑ is the positive part of the soft-thresholded total fluorescence. This encourages elements ofẑ to be exactly zero for frames in which the dictionary element has low fluorescence.
• α = 0: In this setting,ẑ is found by scaling all elements of the total fluorescence by the same amount. Thus, individual elements ofẑ are not encouraged to be 0, thoughẑ = 0 if the dictionary element has a low amount of fluorescence across all frames (i.e., Y ã f 2 small) or λ is very large.
In Figure 11, we illustrate the values ofẑ for the three scenarios described above.
The proof of Lemma 7.2 is in Appendix A.7. Our approach to solving (8) depends on the size of N s . For |N s | = 1, we can simply use the closed-form solution for z Ns,· 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 |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, N 1 , . . . , 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, . . . , S. Details are provided in Algorithm 1.
We typically solve (5) for a sequence of exponentially decreasing λ values. To improve computational performance, Step 2(b) of Algorithm 1 can be implemented using warm starts, in which Z

(0)
Ns,· is initialized as the solution for Z Ns,· for the previous value of λ. 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.

Scaling of A f
In (5), the kth column of the matrixÃ f encodes the spatial mapping of the kth filtered dictionary element, after scaling. To obtainÃ f , we divide the kth column of A f by a f ·,k 2 2 , the number of pixels in the kth 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 a f ·,k 2 2 to be independent of the largest value of λ for whichẑ k,· = 0. The following lemma supports this particular scaling of the columns of A f . Figure 12: We solve (5) for different scalings ofÃ f . For each spatial component, we note the value of λ at which the spatial component enters the model (i.e., the largest λ for whichẑ k,· = 0). We plot the size of each spatial component versus the value of λ at which the spatial component enters for There is a high correlation in the scatterplots in panels (a) and (b), but little correlation in (c). This motivates us to use the scalingã f ·,k = a f ·,k / a f ·,k 2 2 in Step 3, so that dictionary elements receive a fair shot of selection by the sparse group lasso (5), regardless of their size.
Lemma 7.3. Suppose Y = A f Z * where the following conditions hold: . . , K f , and (a f ·,2 ) a f ·,k = 0 for k = 3, . . . , K f and (ii) Z * ∈ R K f ×T with z * 1,· = P z * 2,· for some T × T permutation matrix P . If we solve (5) for Z withÃ f such thatã f ·,k = a f ·,k / a f ·,k 2 2 , thenẑ 1,· = 0 if and only ifẑ 2,· = 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 λ. In Figure 12, we provide empirical evidence for the chosen scaling of A f .

Sparsity of the Solution
We now consider the range of λ for which the solution to (5) is completely sparse (i.e.,Ẑ = 0) for a fixed value of α.
Corollary 7.5. For any α ∈ (0, 1), if , Figure 13: We plot the value of the objective of (5) atẐ(λ), the minimizer of (5) at λ, for a replicate of data as λ varies. We compare two ways of finding a λ large enough such thatẐ(λ) = 0, which results in the objective shown as . We take λ that satisfies Lemma 7.4 ( ) or λ as defined in Corollary 7.5 ( ). The former ( ) gives the smallest λ such thatẐ(λ) = 0. We can see this from the fact that the line ( ) is on the boundary of the shaded box, which indicates the range of λ for which the the objective value at Z(λ) equals the objective value at 0, i.e., the range of λ for whichẐ(λ) = 0.
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 λ.

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, 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 A f ∈ R P ×K f : Then, defineã f ·,k ≡ a f ·,k / a f ·,k 2 , and consider solving (5) for Z with α < 1. Then,ẑ 3,· = 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, a f ·,1 and a f ·,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.

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. Figure 14: For the example calcium imaging video, we plot the median fluorescence of the pixels within each frame, along with the smoothing spline fit ( ) that is used to correct for the bleaching effect.

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 × 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 ∆f /f transformation [Ahrens et al., 2013, Grewe et al., 2010, Grienberger and Konnerth, 2012. For the ith pixel in the jth frame, the standardized fluorescence is equal to , where y 0 i,j is the fluorescence (after smoothing and bleaching correction) of the ith pixel in the jth frame. This differs from the typical ∆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): , and d max 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)   : In (a), we show a sample frame from the raw example calcium imaging video. In (b), we show the same frame after spatial and temporal smoothing has been done and the bleaching effect has been removed. In (c), we show the frame after the ∆f /f transformation has been performed, which is the final result of the pre-processing. Figure 16: We report our spatial dissimilarity metric, as defined in (3), and three alternative metrics, defined in (12), for four pairs of preliminary dictionary elements. We also report the Euclidean distance between each pair of elements. 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, a 0 ·,i − a 0 ·,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. Figure 17: We focus on a single cluster of preliminary dictionary elements from Step 2. The representative dictionary element for this cluster is highlighted in (a). The spatial maps for a random subset of the 36 preliminary dictionary elements in this cluster are shown in (b), along with a spatial map containing all 36 dictionary elements, in which the hue intensity at a given pixel indicates the number of elements containing that pixel. In (c), we plot the average thresholded fluorescence for each of the dictionary elements from (b). Finally, in (d), we show the representative element for that cluster. The gray coloring indicates the union of all preliminary dictionary elements in the cluster.

A.4 Selecting λ for (5)
To choose λ for (5) via a validation set approach, we perform the following steps: 1. ObtainÃ f ∈ R P ×K f by dividing the kth column of A f by a f ·,k 3. Using Algorithm 1, solve (5) on the training set of pixels for a decreasing sequence of 20 λ values, λ 1 , . . . , λ 20 : 4. For each λ i , calculate the validation error, We use a thresholded version of 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 λ as this is the largest value of λ that results in a validation error within 5% of the minimum validation error achieved by any value of λ considered.
5. Solve (5) on all pixels: Figure 18: On the example video considered in Sections 4 and 6.1, we manually inspected each dictionary element resulting from Step 2 of SCALPEL, by examining the frames from which each dictionary element was derived. Based on this manual inspection, we classified each of the 50 dictionary elements as a "neuron" or a "non-neuron". Next, we considered whether simply filtering each dictionary element based on the number of elements in its cluster (as described at the beginning of Step 3 of SCALPEL) would accurately distinguish between "neurons" and "non-neurons". In the figure, the y-axis shows the percentage of "neurons" that would remain after filtering ( ), and the percentage of "non-neurons" that would be eliminated via filtering ( ), as a function of the filtering threshold (shown on the x-axis). We find that in this video, a careful manual analysis of each dictionary element yields very similar results to simply filtering each dictionary element based on the number of elements in its cluster.
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 N 1 , . . . , N S in order to select a different value of λ for each group, or for all groups at once to select a single value of λ. By following steps similar to those described above, λ 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.
Proof. Letβ = argmin β≥0 1 2 y − β 2 2 + λ β 2 and C = {i : y i ≥ 0}. First, we showβ −C = 0. In anticipation of contradiction, assume there exists j such that j / ∈ C andβ j > 0. Defineβ This is a contradiction, so we conclude thatβ i = 0 for all i / ∈ C. It remains to solve By a result in Section 3.1 of Simon et al. [2013], the solution to (13) without the non-negativity constraint on β C is 1 − λ y C 2 + y 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 Note that solving (14) is equivalent to solving (6), as z 1 = 1 z when z ≥ 0. By algebraic manipulation, we can show that where C is a constant that does not depend on z. Therefore, the solution to (14) is the same as the solution to minimize z≥0 1 2 We solve (15)  The last equality follows from the condition of the lemma, which guarantees thatÃ f M(Ns),N s = 0 for all s = s .
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 λ that guarantees thatẐ Ns,· , the solution to (8) Recall that if |N s | > 1, we iteratively solve forẐ Ns,· using Step 2(b) of Algorithm 1. We initialize at the sparse solution Z 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) for l = 1, . . . , T . Thus, when at least one of these two conditions is satisfied for all k = 1, . . . , K f , then the solution to (5) will be sparse.