Andrea Giovannucci
email: agiovannucci@atironinstitute.org

Johannes Friedrich

Pat Gunn

Jérémie Kalfon

Ann Koay

Jiannis Taxidis

Farzaneh Naja

Jeffey L Gauthier

David W Tank

Dmitri Chklovskii

Sue Ann Koay

Je Rey

L Gauthier

Pengcheng Zhou

Eftychios A Pnevmatikakis
email: epnevmatikakis@atironinstitute.org

CaImAn: An open source tool for scalable Calcium Imaging data Analysis

come   L'archive ouverte pluridisciplinaire

Introduction

Understanding the function of neural circuits is contingent on the ability to accurately record and modulate the activity of large neural populations. Optical methods based on the uorescence activity of genetically encoded calcium binding indicators [START_REF] Chen | Ultrasensitive uorescent proteins for imaging neuronal activity[END_REF] have become a standard tool for this task, due to their ability to monitor in vivo targeted neural populations from many di erent brain areas over extended periods of time (weeks or months). Advances in microscopy techniques facilitate imaging larger brain areas with ner time resolution, producing an everincreasing amount of data. A typical resonant scanning two-photon microscope produces data at a rate greater than 50GB/Hour 1 , a number that can be signi cantly higher (up to more than 1TB/Hour) with other custom recording technologies [START_REF] Sofroniew | A large eld of view two-photon mesoscope with subcellular resolution for in vivo imaging[END_REF] This increasing availability and volume of calcium imaging data calls for automated analysis methods and reproducible pipelines to extract the relevant information from the recorded movies, uorescence and/or neural activity (spikes). The typical steps arising in the processing pipelines are the following (Fig. 1a): i) Motion correction, where the FOV at each data frame (image or volume) is registered against a template to correct for motion artifacts due to the nite scanning rate and existing brain motion, ii) source extraction where the di erent active and possibly overlapping sources are extracted and their signals are demixed from each other and from the background neuropil signals (Fig. 1b), and iii) activity deconvolution, where the neural activity of each identi ed source is deconvolved from the dynamics of the calcium indicator.

Related work

Source extraction

Some source extraction methods attempt the detection of neurons in static images using supervised or unsupervised learning methods. Examples of unsupervised methods on summary images include graph-cut approaches applied to the correlation image [START_REF] Kaifosh | SIMA: Python software for analysis of dynamic uorescence imaging data[END_REF]Spaen et al., 2017), and dictionary learning [START_REF] Pachitariu | Extracting regions of interest from biological images with convolutional sparse block coding[END_REF]. Supervised learning methods based on deep neural networks have also been applied to the problem of neuron detection [START_REF] Apthorpe | Automatic Neuron Detection in Calcium Imaging Data Using Convolutional Networks[END_REF][START_REF] Klibisz | Fast, Simple Calcium Imaging Segmentation with Fully Convolutional Networks[END_REF]. While these methods can be e cient in detecting the locations of neurons, they cannot infer the underlying activity nor do they readily o er ways to deal with the spatial overlap of di erent components.

To extract temporal traces together with the spatial footprints of the components one can use methods that directly represent the full spatio-temporal data in a matrix factorization setup e.g., independent component analysis (ICA) [START_REF] Mukamel | Automated analysis of cellular signals from large-scale calcium imaging data[END_REF], constrained nonnegative matrix factorization (CNMF) [START_REF] Pnevmatikakis | Simultaneous Denoising, Deconvolution, and Demixing of Calcium Imaging Data[END_REF] (and its adaptation to one-photon data [START_REF] Zhou | E cient and accurate extraction of in vivo calcium signals from microendoscopic video data[END_REF]), clustering based approaches [START_REF] Pachitariu | Suite2p: beyond 10,000 neurons with standard two-photon microscopy[END_REF], dictionary learning (Petersen et al., 2017), or active contour models [START_REF] Reynolds | ABLE: An Activity-Based Level Set Segmentation Algorithm for Two-Photon Calcium Imaging Data[END_REF]. Such spatio-temporal methods are unsupervised, and focus on detecting active neurons by considering the spatio-temporal activity of a component as a contiguous set of pixels within the FOV that are correlated in time. While such methods tend to o er a direct decomposition of the data in a set of sources with activity traces in an unsupervised way, in principle they require processing of the full dataset, and thus can be rendered intractable very quickly. Possible approaches to deal with the data size include distributed processing in High Performance Computing (HPC) clusters [START_REF] Freeman | Mapping brain activity at scale with cluster computing[END_REF], spatio-temporal decimation (Friedrich et al., 2017a), and dimensionality reduction [START_REF] Pachitariu | Suite2p: beyond 10,000 neurons with standard two-photon microscopy[END_REF]. Recently, Giovannucci et al. (2017) prototyped an online algorithm (ONACID), by adapting matrix factorization setups [START_REF] Pnevmatikakis | Simultaneous Denoising, Deconvolution, and Demixing of Calcium Imaging Data[END_REF][START_REF] Mairal | Online learning for matrix factorization and sparse coding[END_REF], to operate on calcium imaging streaming data and thus natively deal with large data rates.

Deconvolution

For the problem of predicting spikes from uorescence traces, both supervised and unsupervised methods have been explored. Supervised methods rely on the use of ground truth data to train or t biophysical or neural network models [START_REF] Theis | Benchmarking spike rate inference in population calcium imaging[END_REF][START_REF] Speiser | Fast amortized inference of neural activity from calcium imaging data with variational autoencoders[END_REF]. Unsupervised methods can be either deterministic, such as sparse non-negative deconvolution [START_REF] Vogelstein | Fast nonnegative deconvolution for spike train inference from population calcium imaging[END_REF][START_REF] Pnevmatikakis | Simultaneous Denoising, Deconvolution, and Demixing of Calcium Imaging Data[END_REF] that give a single estimate of the deconvolved neural activity, or probabilistic, that aim to also characterize the uncertainty around these estimates (e.g., [START_REF] Pnevmatikakis | Bayesian spike inference from calcium imaging data[END_REF][START_REF] Deneux | Accurate spike estimation from noisy calcium signals for ultrafast three-dimensional imaging of large neuronal populations in vivo[END_REF]). A recent community benchmarking e ort [START_REF] Berens | Community-based benchmarking improves spike inference from two-photon calcium imaging data[END_REF] characterizes the similarities and di erences of various available methods.

CAIMAN

Here we present CAIMAN, an open source suite for the analysis pipeline of both two-photon and onephoton calcium imaging data. CAIMAN includes frameworks for both o ine analysis (CAIMAN BATCH) where all the data is processed at once at the end of experiment, and online analysis on streaming data (CAIMAN ONLINE). Moreover, CAIMAN requires very moderate computing infrastructure (e.g., a 2 of 40

Overview of analysis pipeline

The standard analysis pipeline for calcium imaging data used in CAIMAN is depicted in Fig. 1a.

The data in movie format is rst processed to remove motion artifacts. Subsequently the active components (neurons and background) are extracted as individual pairs of a spatial footprint that describes the shape of each component projected to the imaged FOV, and a temporal trace that captures its uorescence activity (Fig. 1b-d). Finally, the neural activity of each uorescence trace is deconvolved from the dynamics of the calcium indicator. These operations can be challenging because of limited axial resolution of 2-photon microscopy (or the much larger integration volume in one-photon imaging). This results in spatially overlapping uorescence from di erent sources and neuropil activity. Before presenting the new features of CAIMAN in more detail, we brie y review how it incorporates existing tools in the pipeline.

Motion Correction

CAIMAN uses the NORMCORRE algorithm [START_REF] Pnevmatikakis | NoRMCorre: An online algorithm for piecewise rigid motion correction of calcium imaging data[END_REF] that corrects nonrigid motion artifacts by estimating motion vectors with subpixel resolution over a set of overlapping patches within the FOV. These estimates are used to infer a smooth motion eld within the FOV for each frame. For two-photon imaging data this approach is directly applicable, whereas for one-photon micro-endoscopic data the motion is estimated on high pass spatially ltered data, a necessary operation to remove the smooth background signal and create enhanced spatial landmarks. The inferred motion elds are then applied to the original data frames.

Source Extraction

Source extraction is performed using the constrained non-negative matrix factorization (CNMF) framework of Pnevmatikakis et al. (2016) which can extract components with spatial overlapping projections (Fig. 1b). After motion correction the spatio-temporal activity of each source can be expressed as a rank one matrix given by the outer product of two components: a component in space that describes the spatial footprint (location and shape) of each source, and a component in time that describes the activity trace of the source (Fig. 1c). The data can be described by the sum of all the resulting rank one matrices together with an appropriate term for the background and neuropil signal and a noise term (Fig. 1d). For two-photon data the neuropil signal can be modeled as a low rank matrix [START_REF] Pnevmatikakis | Simultaneous Denoising, Deconvolution, and Demixing of Calcium Imaging Data[END_REF]. For microendoscopic data the larger integration volume leads to more complex background contamination [START_REF] Zhou | E cient and accurate extraction of in vivo calcium signals from microendoscopic video data[END_REF]. Therefore, a more descriptive model is required (see Methods and Materials (Mathematical model of the CNMF framework) for a mathematical description). CAIMAN BATCH embeds these approaches into a general algorithmic framework that enables scalable automated processing with improved results in terms of quality and processing speed.

Deconvolution

Neural activity deconvolution is performed using sparse non-negative deconvolution [START_REF] Vogelstein | Fast nonnegative deconvolution for spike train inference from population calcium imaging[END_REF][START_REF] Pnevmatikakis | Simultaneous Denoising, Deconvolution, and Demixing of Calcium Imaging Data[END_REF] and implemented with both the near-online OASIS algorithm (Friedrich et al., 2017b) and an e cient convex optimization framework [START_REF] Pnevmatikakis | Simultaneous Denoising, Deconvolution, and Demixing of Calcium Imaging Data[END_REF]. The algorithm is competitive to the state of the art according to recent benchmarking studies [START_REF] Berens | Community-based benchmarking improves spike inference from two-photon calcium imaging data[END_REF]. Prior to deconvolution, the traces are detrended to remove non-stationary e ects, e.g., photo-bleaching.

Online Processing

The three processing steps described above can be implemented in an online fashion on streaming data using the ONACID algorithm [START_REF] Giovannucci | OnACID: Online analysis of calcium imaging data in real time[END_REF]. The method builds upon the online dictionary learning framework presented in Mairal et al. (2010) for source extraction, by adding the capability of nding new components as they appear and also incorporating the steps of motion correction and deconvolution (Fig. 1e). CAIMAN ONLINE extends and improves the ONACID prototype 4 of 40 algorithm by introducing a number of algorithmic features and a CNN based component detection approach, leading to a major performance improvement.

We now present the new methods introduced by CAIMAN. More details are given in Methods and Materials and pseudocode descriptions of the main routines are given in the Appendix.

Batch processing of large scale datasets on standalone machines

The batch processing pipeline mentioned above can become a computational bottleneck when tackled without customized solutions. For instance, a naive approach to the problem might have as a rst step to load in-memory the full dataset; this approach is non-scalable as datasets typically exceed available RAM (and extra memory is required by any analysis pipeline). To limit memory usage, as well as computation time, CAIMAN BATCH relies on a MapReduce approach [START_REF] Dean | MapReduce: simpli ed data processing on large clusters[END_REF]. Unlike previous work [START_REF] Freeman | Mapping brain activity at scale with cluster computing[END_REF], CAIMAN BATCH assumes minimal computational infrastructure (up to a standard laptop computer), is not tied to a particular parallel computation framework, and is compatible with HPC scheduling systems like SLURM [START_REF] Yoo | Slurm: Simple linux utility for resource management[END_REF].

Naive implementations of motion correction algorithms need to either load in memory the full dataset or are constrained to process one frame at a time, therefore preventing parallelization. Motion correction is parallelized in CAIMAN BATCH without signi cant memory overhead by processing several temporal chunks of a video data on di erent CPUs. CAIMAN BATCH broadcasts to each CPU a meta-template, which is used to align all the frames in the chunk. Each process writes in parallel to the target le containing motion-corrected data, which is stored in as a memory mapped array. This allows arithmetic operations to be performed against data stored on the hard drive with minimal memory use, and slices of data to be indexed and accessed without loading the full le in memory.

More details are given in Methods and Materials (Memory mapping).

Similarly, the source extraction problem, especially in the case of detecting cell bodies, is inherently local with a neuron typically appearing in a neighborhood within a small radius from its center of mass (Fig. 2a). Exploiting this locality, CAIMAN BATCH splits the FOV into a set of spatially overlapping patches which enables the parallelization of the CNMF (or any other) algorithm to extract the corresponding set of local spatial and temporal components. The user speci es the size of the patch, the amount of overlap between neighboring patches and the initialization parameters for each patch (number of components and rank background for CNMF, stopping criteria for CNMF-E). Subsequently the patches are processed in parallel by the CNMF/CNMF-E algorithm to extract the components and neuropil signals from each patch.

Apart from harnessing memory and computational bene ts due to parallelization, processing in patches acts indirectly as a dynamic range equalizer and enables CAIMAN BATCH to detect neurons across the whole FOV, a feature absent in the original CNMF, where areas with high absolute uorescence variation tend to be favored. This results in better source extraction performance.

After all the patches have been processed, the results are embedded within the FOV (Fig. 2a), and the overlapping regions between neighboring patches are processed so that components corresponding to the same neuron are merged. The process is summarized in algorithmic format in Alg. 1 and more details are given in Methods and Materials (Combining results from di erent patches).

Initialization Methods

Initialization methods for matrix factorization problems can impact results due to the non-convex nature of their objective function. CAIMAN BATCH provides an extension of the GREEDYROI method used in Pnevmatikakis et al. (2016), that detects neurons based on localized spatiotemporal activity.

CAIMAN BATCH can also be seeded with binary masks that are obtained from di erent sources, e.g., through manual annotation or segmentation of structural channel (SEEDEDINITIALIZATION, Alg. 2).

More details are given in Methods and Materials (Initialization strategies).

of 40

Automated component evaluation and classi cation

A common limitation of matrix factorization algorithms is that the number of components that the algorithm seeks during its initialization must be pre-determined by the user. For example, Pnevmatikakis et al. (2016) suggest a large number of components which are then heuristically ordered according to their size and activity pattern. When processing large datasets in patches the target number of components is passed on to every patch implicitly assuming a uniform density of (active) neurons within the entire FOV. In general this assumption does not hold and can generate a large number of spurious components. CAIMAN introduces tests to assess the quality of the detected components and eliminate possible false positives. These tests are based on the observation that active components are bound to have a distinct localized spatio-temporal signature within the FOV. We present below unsupervised and supervised tests employed by CAIMAN for component classi cation. In CAIMAN BATCH, they are initially applied after the processing of each patch is completed, and additionally as a post-processing step after the results from the patches have been merged and re ned, whereas in CAIMAN ONLINE they are used to screen new candidate components. We brie y present these tests below and refer to Methods and Materials (Details of quality assessment tests) for more details: Spatial footprint consistency: To test whether a detected component is spurious, we correlate the spatial footprint of this component with the average frame of the data, taken over the interval when the component, with no other overlapping component, was active (Fig. 2c).

The component is rejected if the correlation coe cient is below a certain threshold ✓ sp (e.g.,

✓ sp < 0.5).
Trace SNR: Similarly, for each component we computed the peak SNR of its temporal trace averaged over the duration of a typical transient (Fig. 2d). The component is rejected if the computed SNR is below a certain threshold ✓ SNR (e.g., ✓ SNR = 2).

CNN based classi cation:

We also trained a 4-layer convolutional neural network (CNN) to classify spatial footprints into true or false components (Fig. 2e), where a true component here corresponds to a spatial footprint that resembles the soma of a neuron. The classi er, named batch classi er, was trained on a small corpus of manually annotated datasets (full description given in section Benchmarking against ground truth) and exhibited similar high classi cation performance on test samples from di erent datasets.

While CAIMAN uses the CNMF algorithm, the tests described above can be applied to results obtained from any source extraction algorithm, highlighting the modularity of our tools.

Online analysis with CAIMAN ONLINE

CAIMAN supports online analysis on streaming data building on the core of the prototype algorithm of Giovannucci et al. (2017), and extending it in terms of qualitative performance and computational e ciency:

Initialization: Apart from initializing CAIMAN ONLINE with CAIMAN BATCH on a small time interval, CAIMAN ONLINE can also be initialized in a bare form over an even smaller time interval, where only the background components are estimated and all the components are determined during the online analysis. This process, named BAREINITIALIZATION, can be achieved by running the CNMF algorithm [START_REF] Pnevmatikakis | Simultaneous Denoising, Deconvolution, and Demixing of Calcium Imaging Data[END_REF] over the small interval to estimate the background components and possibly a small number of components. The SEEDEDINITIALIZATION of Alg. 2 can also be used.

Deconvolution: Instead of a separate step after demixing as in Giovannucci et al. (2017), deconvolution here can be performed simultaneously with the demixing online, leading to more stable traces especially in cases of low-SNR, as also observed in Pnevmatikakis et al. (2016).

Online deconvolution can also be performed for models that assume second order calcium dynamics, bringing the full power of Friedrich et al. (2017b) to processing of streaming data.

Epochs: CAIMAN ONLINE supports multiple passes over the data, a process that can detect early activity of neurons that were not picked up during the initial pass, as well as smooth the activity of components that were detected at late stages during the rst epoch.

New component detection using a CNN:

To search for new components in a streaming setup, ONACID keeps a bu er of the residual frames, computed by subtracting the activity of already found components and background signals. Candidate components are determined by looking for points of maximum energy in this residual signal, after some smoothing and dynamic range equalization. For each such point identi ed, a candidate shape and trace are constructed using a rank-1 NMF in a local neighborhood around this point. In its original formulation [START_REF] Giovannucci | OnACID: Online analysis of calcium imaging data in real time[END_REF], the shape of the component was evaluated using the space correlation test described above. Here, we introduce a CNN classi er approach that tests candidate components by examining their spatial footprint as obtained by the average of the residual bu er across time. This online classi er (di erent from the batch classi er for quality assessment described above), is trained to be strict, minimizing the number of false positive components that enter the online processing pipeline.

Component registration across multiple sessions

CAIMAN provides a method to register components from the same FOV across di erent sessions.

The method uses a simple intersection over union metric to calculate the distance between di erent cells in di erent sessions and solving a linear assignment problem to perform the registration in a fully automated way (REGISTERPAIR, Alg. 5). To register the components between more than 2 sessions (REGISTERMULTI, Alg. 6), we order the sessions chronologically and register the components of the current session against the union of component of all the past sessions aligned to the current FOV. This allows for the tracking of components across multiple sessions without the need of pairwise registration between each pair of sessions. More details as well as discussion of other methods [START_REF] Sheintuch | Tracking the Same Neurons across Multiple Days in Ca2+ Imaging Data[END_REF] are given in Methods and Materials (Component registration).

Benchmarking against ground truth

To quantitatively evaluate CAIMAN we benchmarked its results against ground truth data.

Creating ground truth data through manual annotation

We collected manual annotations from multiple independent labelers who were instructed to nd round or donut shaped2 active neurons on 9 two-photon in vivo mouse brain datasets. The datasets were collected at various labs and from various brain areas (hippocampus, visual cortex, parietal cortex) using several GCaMP variants. A summary of the features of all the annotated datasets is given in Table 2. Details about the annotation procedure are given in Methods and Materials.

To address human variability in manual annotation each dataset was labeled by 3 or 4 independent labelers, and the nal ground truth dataset was created by having the di erent labelers reaching a consensus over their disagreements (Fig. 3a). The result of this process was de ned as ground truth for the evaluation of CAIMAN as well as each individual labeler against the consensus (Fig. 3b) 3 . More details are given in Methods and Materials (Collection of manual annotations and .

ground truth). We believe that the current database, which will be made publicly available, presents an improvement over the existing neuro nder database (http://neuro nder.codeneuro.org/) in several aspects:

Consistency: The datasets are annotated using exactly the same procedure (see Methods and Materials), and in all datasets the goal is to detect only active cells. In contrast, the annotation of the various neuro nder datasets is performed either manually or automatically by segmenting an image of a static (structural) indicator. Even though structural indicators could be used for ground truth extraction, the segmentation of such images is not a straightforward problem in the case of dense expression, and the stochastic expression of indicators can lead to mismatches between functional and structural indicators.

Uncertainty quanti cation:

By employing more than one human labeler we discovered a surprising level of disagreement between di erent annotators (see Table 1, Fig. 3b for details), which renders individual annotations somewhat unreliable for benchmarking purposes, and non-reproducible. The combination of the various annotations leads to more reliable ground truth and also quanti es the limits of human performance.

Comparing CAIMAN against ground truth

To compare CAIMAN against the consensus ground truth, the manual annotations were used as binary masks to construct the ground truth spatial and temporal components, using the SEEDEDINI-TIALIZATION procedure (Alg. 2) of CAIMAN BATCH. The set of spatial footprints obtained from CAIMAN is registered against the set of ground truth spatial footprints (derived as described above) using the REGISTERPAIR algorithm (Alg. 5) for component registration described above. Performance is then quanti ed using a precision/recall framework similar to other studies (

Results

Manual annotations show a high degree of variability

We compared the performance of each human annotator against a consensus ground truth. The performance was quanti ed with a precision/recall framework and the results of the performance of each individual labeler against the consensus ground truth for each dataset is given in Table 1.

The range of human performance in terms of F 1 score was 0.69-0.94, with average 0.83± 0.07 (mean ± STD). All annotators performed similarly on average (0.83±0.05, 0.83±0.08, 0.84±0.06, 0.85±0.08).

We also ensured that the performance of labelers was stable across time (i.e. their learning curve plateaued, data not shown). As shown in Table 1 (see also Fig 4b) the F 1 score was never 1, and in most cases it was less or equal to 0.9, demonstrating signi cant variability between annotators. 4 All future development of CAIMAN will be in Python 3, eventually rendering it incompatible with Python 2.x.

11 of 40

CAIMAN BATCH and CAIMAN ONLINE detect neurons with near-human accuracy

We rst benchmarked CAIMAN BATCH and CAIMAN ONLINE against consensus ground truth for the task of identifying neurons locations and their spatial footprints, using the same precision recall framework (Table 1). Fig. 4a shows an example dataset (K53) along with neuron-wise matches and mismatches between CAIMAN BATCH and consensus ground truth (top) and CAIMAN ONLINE vs consensus ground truth (bottom).

The results indicate a similar performance between CAIMAN BATCH and CAIMAN ONLINE; CAIMAN BATCH has F 1 scores in the range 0.68-0.79 and average performance 0.75±0.04 (mean±STD). On the other hand CAIMAN ONLINE had F 1 scores in the range 0.68-0.82 and average performance 0.76±0.04.

While the two algorithms performed similarly on average, CAIMAN BATCH tends to perform better for shorter datasets whereas online processing tends to lead to better results for longer datasets (see Table 2 for characteristics of the various datasets). CAIMAN approaches but is in most cases below the accuracy levels of human annotators (Fig. 4b). This can be attributed to a number of reasons:

First, to demonstrate the generality and ease of use of our tools, the results presented here are obtained by running CAIMAN BATCH and CAIMAN ONLINE with exactly the same parameters for each dataset (see Methods and Materials (Implementation details)): ne-tuning to each individual dataset can signi cantly increase performance. Second, CNMF detects active components regardless of their shape, and can detect non-somatic structures with signi cant transients. While non-somatic components can be ltered out to some extent using the CNN classi er, their existence degrades performance compared to the ground truth that consists only of neurons. Lastly, the ground truth is by construction a subset of the union of all individual annotations, which can bias upwards the scores of individual labelers.

Neurons with higher SNR transients are detected more accurately

While CAIMAN ONLINE had balanced performance with respect to precision and recall (mean precision 0.77±0.05, mean recall 0.76±0.07), CAIMAN BATCH showed signi cantly higher precision than recall (mean precision 0.83±0.09, mean recall 0.69±0.08). We looked into this behavior, by analyzing CAIMAN BATCH performance as a function of the SNR of the inferred and ground truth traces (Fig. 4c-d). The SNR measure of a trace corresponds to the peak-SNR averaged over the length of a typical trace (see Methods and Materials (Detecting uorescence traces with high SNR)). An example is shown in Fig. 4c where the scatter plot of SNR between matched ground truth and inferred traces is shown (false negative/positive components are shown along the x-and y-axis, respectively). To evaluate the performance we computed a precision metric as the fraction of inferred components Both distributions exhibit a sharp derivative close 1 (last bin), with the batch approach giving better results.

above a certain SNR threshold that are matched with a ground truth component (Fig. 4c, shaded blue). Similarly we computed a recall metric as the fraction of ground truth components above a SNR threshold that are detected by CAIMAN BATCH (Fig. 4c, shaded red), and an F 1 score as the harmonic mean of the two (Fig. 4d). The results indicate that the performance signi cantly grows as a function of the SNR for all datasets considered, growing on average from 0.73 when all neurons are considered to 0.92 when only neurons with traces having SNR g 9 are considered (Fig. 4d) 5 .

CAIMAN reproduces the ground truth traces with high delity

Testing the quality of the inferred traces is a more challenging task due to the complete lack of ground truth data in the context of large scale in vivo recordings. As mentioned above, we considered as ground truth the traces obtained by running the CNMF algorithm seeded with the Manuscript submitted binary masks obtained by consensus ground truth procedure. After alignment of the ground truth with the results of CAIMAN, the matched traces were compared both for CAIMAN BATCH and for CAIMAN ONLINE. Fig. 5a, shows an example of 5 of these traces for the dataset K53, showing very similar behavior of the traces in these three di erent cases.

To quantify the similarity we computed the correlation coe cients of the traces (ground truth vs CAIMAN BATCH, and ground truth vs CAIMAN ONLINE) for all the 9 datasets (Fig. 5b-c). Results indicated that for all but one dataset (Fig. 5b) CAIMAN BATCH reproduced the traces with higher delity, and in all cases the mean correlation coe cients was higher than 0.9, and the empirical histogram of correlation coe cients peaked at the maximum bin 0.99-1 (Fig. 5c). The results indicate that the batch approach extracts traces closer to the ground truth traces. This can be attributed to a number of reasons: By processing all the time points simultaneously, the batch approach can smooth the trace estimation over the entire time interval as opposed to the online approach where at each timestep only the information up to that point is considered. Moreover, CAIMAN ONLINE might not detect a neuron until it becomes strongly active. This neuron's activity before detection is unknown and has a default value of zero, resulting in a lower correlation coe cient. While this can be ameliorated to a great extent with additional passes over the data, the results indicate trade-o s between using the online and o ine versions of CAIMAN.

Online analysis of a whole brain zebra sh dataset

We tested CAIMAN ONLINE with a 380GB whole brain dataset of larval zebra sh (Danio rerio) acquired with a light-sheet microscope [START_REF] Kawashima | The serotonergic system tracks the outcomes of actions to mediate short-term motor learning[END_REF]. The imaged transgenic sh (Tg(elavl3:H2B-GCaMP6f)jf7) expressed the genetically encoded calcium indicator GCaMP6f in almost all neuronal nuclei. Data from 45 planes (FOV 820x410 m 2 , spaced at 5.5 m intervals along the dorso-ventral axis) was collected at 1Hz for 30 minutes (for details about preparation, equipment and experiment refer to Kawashima et al. (2016)). With the goal of simulating real-time analysis of the data, we run all the 45 planes in parallel on a computing cluster with 9 nodes (each node is equipped with 24

CPUs and 128-256 GB RAM). Data was not stored locally in each machine but directly accessed from a network drive.

The algorithm was initialized with CAIMAN BATCH run on 200 initial frames and looking for 500 components. The small number of frames (1885) and the large FOV size (2048 ù 1188 pixels) for this dataset motivated this choice of increased number of components during initialization. In Fig. 6 we report the results of the analysis for plane number 11 of 45. For plane 11, CAIMAN ONLINE found 1524 neurons after processing 1685 frames. Since no ground truth was available for this dataset, it was only possible to evaluate the performance of this algorithm by visual inspection. CAIMAN ONLINE identi ed all the neurons with a clear footprint in the underlying correlation image (higher SNR, Fig. 6a) and missed a small number of the fainter ones (low SNR). By visual inspection of the components the authors could nd very few false positives. Given that the parameters were not tuned and that the classi er was not trained on zebra sh neurons, we hypothesize that the algorithm is biased towards a high precision result. Spatial components displayed the expected morphological features of neurons (Fig. 6b-c). Considering all the planes (Figs 6e and 11) CAIMAN ONLINE was able to identify in a single pass of the data a total of 66108 neurons. See Supplemental Movie 3 for a summary across all planes. The analysis was performed in 21 minutes, with the rst 3 minutes allocated to the initialization and the remaining 18 to process the rest of the data in streaming mode (and in parallel for each plane). This demonstrates the ability of CAIMAN ONLINE to process large amounts of data in real-time (see also Fig. 8 for a discussion of computational performance).

Analyzing 1p microendoscopic data using CAIMAN

We tested the CNMF-E implementation of CAIMAN BATCH on in vivo microendosopic data from mouse dorsal striatum, with neurons expressing GCaMP6f. 6000 frames were acquired at 30 frames per second while the mouse was freely moving in an open eld arena (for further details refer to Zhou et al. (2018)). In Fig. 7 we report the results of the analysis using CAIMAN BATCH with patches and compare to the results of the MATLAB ® implementation of Zhou et al. (2018). Both implementations detect similar components (Fig. 7a) with an F 1 -score of 0.89. 573 neurons were found by both implementations. 106 and 31 additional components were detected by Zhou et al.

(2018) and CAIMAN BATCH respectively. The median correlation between the temporal traces of neurons detected by both implementations was 0.86. Similar results were also obtained by running CAIMAN without patches. Ten example temporal traces are plotted in Fig. 7b.

Computational performance of CAIMAN

We examined the performance of CAIMAN in terms of processing time for the various analyzed datasets presented above (Fig. 8). The processing time discussed here excludes motion correction, which is highly e cient and primarily depends on the level of the FOV discretization for non-rigid motion correction [START_REF] Pnevmatikakis | NoRMCorre: An online algorithm for piecewise rigid motion correction of calcium imaging data[END_REF]. For CAIMAN BATCH, each dataset was analyzed using three di erent computing architectures: i) a single laptop (MacBook Pro) with 8

CPUs and 16GB of RAM (blue in Fig. 8a), ii) a linux-based workstation (CentOS) with 24 CPUs and 128GB of RAM (magenta), and iii) a linux-based HPC cluster (CentOS) where 112 CPUs (4 nodes, 28

CPUs each) were allocated for the processing task (yellow). Fig. 8a shows the processing of CAIMAN BATCH as a function of dataset size on the 5 longest datasets, whose size exceeded 8GB, on log-log plot.

Results show that, as expected, employing more processing power results in faster processing.

CAIMAN BATCH on a HPC cluster processes data faster than acquisition time (Fig. 8a) even for very the cost of detecting and incorporating new neurons, and iii) the cost of periodic updates of spatial footprints. Fig. 8b-right shows that the two rst steps, which are required for each frame, can be done in real-time. In Fig. 8c the cost per frame is plotted for the analysis of the whole brain zebra sh recording. The lower imaging rate (1Hz) allows for the tracking of neural activity to be done with computational cost signi cantly lower than the 1 second between volume imaging time (Fig. 8c), even in the presence of a large number of components (typically more than 1000 per plane, Fig. 6) and the signi cantly larger FOV (2048 ù 1188 pixels). As expected the cost of updating spatial footprints can be signi cantly larger if done simultaneously for all components (Fig. 8c, bottom).

However, the average cost of updating a single spatial footprint is roughly 8ms, enabling real-time processing for each frame, when this step is evenly distributed among di erent frames/volumes, or is performed by a parallel independent process [START_REF] Giovannucci | OnACID: Online analysis of calcium imaging data in real time[END_REF].

The cost of processing 1p data in CAIMAN BATCH using the CNMF-E algorithm [START_REF] Zhou | E cient and accurate extraction of in vivo calcium signals from microendoscopic video data[END_REF] is shown (Fig. 8d) for the workstation hardware. Splitting in patches and processing in parallel can lead to computational gains at the expense of increased memory usage. This is because the CNMF-E introduces a background term that has the size of the dataset and needs to be loaded and updated in memory in two copies. This leads to processing times that are slower compared to the standard processing of 2p datasets, and higher memory requirements. However, as Fig. 8d demonstrates, memory usage can be controlled enabling scalable inference at the expense of slower processing speeds.

CAIMAN successfully tracks neurons across multiple days

Fig. 9 shows an example of tracking neurons across 6 di erent sessions corresponding to 6 di erent days of mouse cortex in vivo data using our multi-day registration algorithm REGISTERMULTI (see Methods,Alg. 6). 453,393,375,378,376, and 373 active components were found in the six sessions, respectively. Our tracking method detected a total of 686 distinct active components. Of these, 172, 108, 70, 92, 82, and 162 appeared in exactly 1, 2, 3, 4, 5, and all 6 sessions respectively. Contour plots of the 162 components that appeared in all sessions are shown in Fig. 9a, and parts of the FOV are highlighted in Fig. 9d showing that components can be tracked in the presence of non-rigid deformations of the FOV between the di erent sessions.

To test the stability of REGISTERMULTI for each subset of sessions, we repeated the same procedure running backwards in time starting from day 6 and ending at day 1, a process that also generated a total of 686 distinct active components. We identi ed the components present in at least a given subset of sessions when using the forward pass, and separately when using the backwards pass, and compared them against each other (Fig. 9b) for all possible subsets.

Results indicate a very high level of agreement between the two approaches with many of the disagreements arising near the boundaries (data not shown). Disagreements near the boundaries can arise because the forward pass aligns the union with the FOV of the last session, whereas the backwards pass with the FOV of the rst session, potentially leading to loss of information near the boundaries.

A step by step demonstration of the tracking algorithm for the rst three sessions is shown in the appendix (Fig. 10). Our approach allows for the comparison of two non-consecutive sessions through the union of components without the need of a direct pairwise registration (Fig. 10f), where it is shown that registering sessions 1 and 3 directly and through the union leads to nearly identical results. Fig. 9c compares the registrations for all pairs of sessions using the forward (red) or the backward (blue) approach, with the direct pairwise registrations. Again, the results indicate a very high level of agreement, indicating the stability and e ectiveness of the proposed approach.

Discussion

Reproducible and scalable analysis for the 99% tion). Scalability is achieved by either using a MapReduce batch approach, which employs parallel processing of spatially overlapping, memory mapped, data patches; or by integrating the online processing framework of Giovannucci et al. (2017) within our pipeline. Apart from computational gains both approaches also result in improved performance. Towards our goal of providing a single package for dealing with standard problems arising in analysis of imaging data, CAIMAN also includes an implementation of the CNMF-E algorithm of Zhou et al. (2018) for the analysis of microendoscopic data, as well as with a novel method for registering analysis results across multiple days.

Towards surpassing human neuron detection performance

To evaluate the performance of CAIMAN BATCH and CAIMAN ONLINE, we generated a corpus of multiply annotated two-photon imaging datasets. The results indicated a surprising level of disagreement between individual labelers, highlighting both the di culty of the problem, and the non-reproducibility of the laborious task of human annotation. CAIMAN reached near-human performance with respect to this ground truth, by using the same parameters for all the datasets without dataset dependent parameter tweaking. Such tweaking could for example include setting the SNR threshold based on the noise level of the recording, the complexity of the neuropil signal based on the level of background activity, or specialized treatment around the boundaries of the FOV to compensate for eventual imaging artifacts.

Apart from being used as a benchmarking tool, the set of manual annotations can also be used as labeled data for supervised learning algorithms. CAIMAN uses two CNN based classi ers trained on (a subset of) this data, one for post processing component classi cation in CAIMAN BATCH, and the other for detecting new neurons in residual images in the CAIMAN ONLINE. The deployment of these classi ers resulted in signi cant gains in terms of performance, and we expect further advances in the future. The annotations will be made freely available to the community upon publication of the paper for benchmarking and training purposes.

CAIMAN BATCH vs CAIMAN ONLINE

Future directions

While CAIMAN uses a highly scalable processing pipeline for two-photon datasets, processing of one-photon microendoscopic imaging data is less scalable due to the more complex background model that needs to be retained in memory during processing. Adapting CAIMAN ONLINE to the onephoton data processing algorithm of Zhou et al. ( 2018) is a promising way for scaling up e cient processing in this case. The continuing development and quality improvement of neural activity indicators has enabled direct imaging of neural processes (axons/dendrites), imaging of synaptic activity [START_REF] Xie | Resolution of high-frequency mesoscale intracortical maps using the genetically encoded glutamate sensor iGluSnFR[END_REF], or direct imaging of voltage activity in vivo conditions [START_REF] Piatkevich | A robotic multidimensional directed evolution approach applied to uorescent voltage reporters[END_REF]. While the approach presented here is tuned for somatic imaging through the use of various assumptions (space localized activity, CNN classi ers trained on images of somatic activity), the technology of CAIMAN is largely transferable to these domains as well. These extensions will be pursued in future work.

Methods and Materials

Memory mapping

In order to e ciently access data in parallel, CAIMAN BATCH relies on memory mapping. With memory mapped (mmap) arrays, arithmetic operations can be performed on data residing on the hard drive without explicitly loading it to RAM, and slices of data can be indexed and accessed without loading the full le in memory, enabling out-of-core processing [START_REF] Toledo | A survey of out-of-core algorithms in numerical linear algebra[END_REF]. The order in which data in a memory mapped le is stored on the hard drive can dramatically a ect the read-write performance of out-of-core operations on spinning disks, and to a lesser degree on solid state drives. On modern computers tensors are stored in linear format, no matter the number of the array dimensions. Therefore, one has to decide which elements of an array are contiguous in memory: in row-major order, consecutive elements of a row ( rst-dimension) are next to each other, whereas in column-major order consecutive elements of a column (last dimension) are contiguous.

Such decisions signi cantly a ect the speed at which data is read or written: in column-major order reading a full column is fast because memory is read in a single sequential block, whereas reading a row is ine cient since only one element can be read at a time and all the data needs to be accessed.

Therefore, the original dataset must be saved in the right order to avoid performance problems.

In the context of calcium imaging datasets, CAIMAN BATCH represents the datasets in a matrix form Y , where each row corresponds to a di erent imaged pixel, and each column to a di erent frame. As a result, a column-major order mmap le enables the fast access of individual frames at a given time, whereas a row-major order les enables the fast access of an individual pixel at all times.

To facilitate processing in patches CAIMAN BATCH stores the data in row-major order. In practice, this is opposite to the order with which the data appears, one frame at a time. In order to reduce memory usage and speed up computation CAIMAN BATCH employs a MapReduce approach, where either multiple les or multiple chunks of a big le composing the original datasets are processed and saved in mmap format in parallel. This operation includes two phases, rst the chunks/ les are saved in multiple row-major mmap format, and then chunks are simultaneously combined into a single large row-major mmap le. In order to reduce preprocessing steps, if the le(s) need to be corrected for motion artifacts, chunks of the registered data can be stored on-the-y during motion correction.

Mathematical model of the CNMF framework

The CNMF framework (Fig. 1d) for calcium imaging data representation can be expressed in mathematical terms as [START_REF] Pnevmatikakis | Simultaneous Denoising, Deconvolution, and Demixing of Calcium Imaging Data[END_REF])

Y = AC + B + E. (1)
Here, Y À R dùT denotes the observed data written in matrix form, where d is the total number of observed pixels/voxels, and T is the total number of observed timesteps (frames). A À R dùN denotes the matrix of the N spatial footprints,

A = [a 1 , a 2 , … , a N ], with a i À R dù1 being the spatial footprint of component i. C À R NùT denotes the matrix of temporal components, C = [c 1 , c 2 , … , c n ] Ò , with c i À R T ù1
B = W (Y * AC), ( 2
)
where W À R dùd is an appropriate weight matrix, where the (i, j) entry models the in uence of the neuropil signal of pixel j to the neuropil signal at pixel i.

Combining results from di erent patches

To combine results from the di erent patches we rst need to account for the overlap at the bound-

min a m g0,c m g0 ÒA old C old * a m c Ò m Ò. (3)
Prior to merging, the value of each component at each pixel is normalized by the number of patches that overlap in this pixel, to avoid counting the activity of each pixel multiple times.

We follow a similar procedure for the background/neuropil signals from the di erent patches.

For the case of two-photon data, the spatial background/neuropil components for each patch can be updated by keeping their spatial extent intact to retain a local neuropil structure, or they can be merged when they are su ciently correlated in time as described above to promote a more global structure. For the case of one-photon data, CNMF-E estimates the background using a local auto-regressive process (see Eq. ( 2)) [START_REF] Zhou | E cient and accurate extraction of in vivo calcium signals from microendoscopic video data[END_REF], a setup that cannot be immediately propagated when combining the di erent patches. To combine backgrounds from the di erent patches, we rst approximate the backgrounds B i from all the patches i with a low rank matrix using non-negative matrix factorization of rank g b to obtain global spatial, and temporal background components.

[b i , f i ] = NNMF(B i , g b ). ( 4
)
The resulting components are embedded into a large matrix B À R dùT that retains a low rank structure. After the components and backgrounds from all the patches have been combined, they are further re ned by running CNMF iteration of updating spatial footprints, temporal traces, and neuropil activity. CAIMAN BATCH implements these steps in a highly parallel fashion (as also described in Pnevmatikakis et al. (2016)): Temporal traces whose corresponding spatial traces do not overlap can be updated in parallel. Similarly, the rows of the matrix of spatial footprints A can also be updated in parallel (2b). The process is summarized in algorithmic format in Alg. 1.

Initialization strategies

Source extraction using matrix factorization requires solving a bi-convex problem where initialization plays a critical role. The CNMF/CNMF-E algorithms use initialization methods that exploit the locality of the spatial footprints to e ciently identify the locations of candidate components [START_REF] Pnevmatikakis | Simultaneous Denoising, Deconvolution, and Demixing of Calcium Imaging Data[END_REF][START_REF] Zhou | E cient and accurate extraction of in vivo calcium signals from microendoscopic video data[END_REF]. CAIMAN incorporates these methods, extending them by using the Manuscript submitted temporal locality of the calcium transient events. The available initialization methods for CAIMAN BATCH include:

GREEDYROI: This approach, introduced in Pnevmatikakis et al. (2016), rst spatially smooths the data with a Gaussian kernel of size comparable to the average neuron radius, and then initializes candidate components around locations where maximum variance (of the smoothed data) is explained. This initialization strategy is fast but requires speci cation of the number of components by the user.

ROLLINGGREEDYROI:

The approach, introduced in this paper, operates like GREEDYROI by spatially smoothing the data and looking for points of maximum variance. Instead of working across all the data, ROLLINGGREEDYROI looks for points of maximum variance on a rolling window of a xed duration, e.g., 3 seconds, and initializes components by performing a rank one NMF on a local spatial neighborhood. By focusing into smaller rolling windows, ROLLINGGREEDYROI can better isolate single transient events, and as a result detect better neurons with sparse activity. ROLLINGGREEDYROI is the default choice for processing of 2-photon data.

GREEDYCORR: This approach, introduced in Zhou et al. (2018), initializes candidate components around locations that correspond to the local maxima of an image formed by the pointwise product between the correlation image and the peak signal-to-noise ratio image. By setting a threshold for acceptance, this approach does not require the prior speci cation of number of components. This comes at the expense of a higher computational cost. GREEDYCORR is the default choice for processing of 1-photon data.

SPARSENMF: Sparse NMF approaches, when ran in small patches, can be e ective for quickly uncovering spatial structure in the imaging data, especially for neural processes (axons/dendrites) whose shape cannot be easily parametrized and/or localized.

Algorithm seeding with binary masks

Often locations of components are known either from manual annotation or from labeled data obtained in a di erent way, such as data from a static structural channel recorded concurrently with the functional indicator. CAIMAN can be seeded with binary (or real valued) masks for the spatial footprints. Apart from A, these masks can be used to initialize all the other relevant matrices C and B as well. This is performed by i) rst estimating the temporal background components f using only data from parts of the FOV not covered by any masks and, ii) then estimating the spatial background components b, and then estimating A, C (with A restricted to be non-zero only at the locations of the binary masks), using a simple NMF approach. Details are given in Alg. 2.

Details of quality assessment tests

Here we present the unsupervised and supervised quality assessment tests in more detail (Fig. 2).

Matching spatial footprints to the raw data constructs intervals around these peaks, (e.g., 50 ms in the past and 300ms in the future, to cover the main part of a possible calcium transient around that point), and then averages Y i across time over the union of these intervals to obtain a spatial image < Y i > (Fig. 2c). The Pearson's correlation over space between < Y i > and a i (both restricted on a small neighborhood around the centroid of a i ) is then computed, and component i is rejected if the correlation coe cient is below a threshold value ✓ sp , (e.g., ✓ sp < 0.5). Note that a similar test is used in the online approach of Giovannucci et al.

(2017) to accept for possible new components.

Detecting uorescence traces with high SNR For a candidate component to correspond to an active neuron its trace must exhibit dynamics reminiscent of the calcium indicator's transient. A criterion for this can be obtained by requiring the average SNR of trace c i over the course a transient to be above a certain threshold ✓ SNR , e.g., ✓ SNR = 2, (Fig. 2d). The average SNR is as a measure of how unlikely it is for the transients of c i (after some appropriate z-scoring) to have been a result of a white noise process.

To compute the SNR of a trace, let R = Y * AC * B be the residual spatiotemporal signal. We can obtain the residual signal for each component i, r i , by projecting R into the spatial footprint a i :

r i = 1 Òa i Ò 2 R Ò a i ( 5
)
Then the trace c i + r i corresponds to the non-denoised trace of component i. To calculate its SNR we rst compute a type of z-score:

z i = c i + r i * BASELINE(c i + r i ) NOISE(c i + r i ) . ( 6
)
The BASELINE( ) function determines the baseline of the trace, which can be varying in the case of long datasets exhibiting baseline trends, e.g., due to bleaching. The function NOISE( ) estimates the noise level of the trace. Since calcium transients around the baseline can only be positive, we estimate the noise level by restricting our attention only to the points t n where c i + r i is below the baseline value, i.e., t n = {t :

c i (t) + r i (t) f BASELINE(c i + r i )},
and compute the noise level as the scale parameter of a half-normal distribution (Fig. 2b):

NOISE(c i + r i ) = std([c i + r i ](t n ))_ u 1 * 2 ⇡ . (7)
We then determine how likely is that the positive excursions of z i can be attributed just to noise. We compute the probabilities p i (t) = (*z i (t)), where ( ) denotes the cumulative distribution function of a standard normal distribution, and compute the most unlikely excursion over a window of N s timesteps that corresponds to the length of a typical transient, e.g., N s = ‰0.4s ù F Â, where 0.4s could correspond to the typical length of a GCaMP6f transient, and F is the imaging rate.

p i min = min t H N s *1 « j=0 p i (t + j) I1_N s . ( 8
)
The (averaged peak) SNR of component i can then be de ned as

SNR i = *1 (1 * p i min ) = * *1 (p i min ), (9)
where *1 is the quantile function for the standard normal distribution (logit function) and a component is accepted if SNR i g ✓ SNR . Note that for numerical stability we compute p i min in the logarithmic domain and check the condition p i min f (*✓ SNR ).

We can also use a similar test for the signi cance of the time traces in the spike domain after performing deconvolution. In this case, traces can be considered as spiking if the maximum height due to a spike transient is signi cantly larger than a threshold. If we assume that the shape of each calcium transient has been normalized to have maximum amplitude 1, then this corresponds to

testing Òs i Ò ÿ g ✓ SNR i
, where s i represents the deconvolved activity trace for component i, and ✓ SNR is again an appropriate SNR threshold, e.g., ✓ SNR = 2, and i is the noise level for trace i.

Classi cation through convolutional neural networks (CNNs)

The tests described above are unsupervised but require ne-tuning of two threshold parameters (✓ sp , ✓ SNR ) that might be dataset dependent and might be sensitive to strong non-stationarities. As a third test we trained a 4-layer CNN to classify the spatial footprints into true or false components, where a true component here corresponds to a spatial footprint that resembles the soma of a neuron (See Fig. 2e and section Classi cation through convolutional networks for details). A simple threshold ✓ CNN can be used to tune the classi er (e.g., ✓ CNN = 0.5).

Collection of manual annotations and ground truth

We collected manual annotations from four independent labelers who were instructed to nd round or donut shaped neurons of similar size using the ImageJ Cell Magic Wand tool Walker (2014).

We focused on manually annotating only cells that were active within each dataset and for that reason the labelers were provided with two summary statistics: i) A movie obtained by removing a running 20th percentile (as a crude background approximation) and downsampling in time by a factor of 10, and ii) the max-correlation image. The correlation image (CI) at every pixel is equal to the average temporal correlation coe cient between that pixel and its neighbors Smith and Häusser (2010) (8 neighbors were used for our analysis). The max-correlation image is obtained by computing the CI for each batch of 33 seconds (1000 frames for a 30Hz acquisition rate), and then taking the maximum over all these images. Neurons that are inactive during the course of the dataset will be suppressed both from the baseline removed video (since their activity will always be around their baseline), and from the max-correlation image since the variation around this baseline will mostly be due to noise leading to practically uncorrelated neighboring pixels. 9 di erent mouse in vivo datasets were used from various brain areas and labs. A description is given in Table 2. To create the consensus ground truth, the labelers were asked to jointly resolve the inconsistencies with each others annotations.

The annotation procedure provides a binary mask per selected component. On the other hand, the output of CAIMAN for each component is a non-negatively valued vector over the FOV (a real-valued mask). The two sets of masks di er not only in their variable type but also in their general shape: Manual annotation through the Cell Magic Wand tool tends to produce circular shapes, whereas the output of CAIMAN will try to accurately estimate the shape of each active component. To construct ground truth that can be directly used for comparison, the binary masks from the manual annotations were used to seed the CNMF algorithm (Alg. 2). This produced a set of ground truth real valued components with spatial footprints restricted to the areas provided by the annotations, and a corresponding set of temporal components that can be used to evaluate the performance of CAIMAN (Fig. 4). Registration was performed using the REGISTERPAIR algorithm (Alg. 5) and match was counted as a true positive when the (modi ed) Jaccard distance (Eq. 11) was below 0.7. Details of the registration procedure are given below (see Component registration).

Classi cation through convolutional neural networks (CNNs)

CAIMAN uses two CNN classi ers; one for post processing component screening in CAIMAN BATCH, and a di erent one for screening candidate components in CAIMAN ONLINE. In both cases a 4 layer CNN was used, with architecture as described in Fig. 2e.

Component registration

Fluorescence microscopy methods enable imaging the same part of the brain across di erent sessions that can span multiple days or weeks. While the microscope can visit the same location in the brain with reasonably high precision, the FOV might might not precisely match due to misalignments or deformations in the brain medium. CAIMAN provides routines for FOV alignment and component registration across multiple sessions/days. Let a 1 1 , a 1 2 , … , a 1

N 1 and a 2 1 , a 2 2 , … , a 2 N 2
the sets of spatial components from sessions 1 and 2 respectively, where N 1 and N 2 denote the total number of components from each session. We rst compute the FOV displacement by aligning some summary images from the two sessions (e.g., mean or correlation image), using some nonrigid registration method, e.g., NoRMCorre [START_REF] Pnevmatikakis | NoRMCorre: An online algorithm for piecewise rigid motion correction of calcium imaging data[END_REF]. We apply the estimated displacement eld to the components of A 1 to align them with the FOV of session 2. To perform the registration, we construct a pairwise distance matrix

D À R N 1 ùN 2 with D(i, j) = d(a 1 i , a 2 j ),
where d( , ) denotes a distance metric between two components. The chosen distance corresponds to the Jaccard distance between the binarized versions of the components. A real valued component a is converted into its binary version m(x) by setting to 1 only the values of a that are above the maximum value of a times a threshold ✓ b , e.g., ✓ b = 0.2:

m(a)((x)) = T 1, a((x)) g ✓ b ÒaÒ ÿ 0, otherwise . ( 10
)
To compute the distance between two binary masks m 1 , m 2 , we use the Jaccard index (intersection over union) which is de ned as

J (m 1 , m 2 ) = m 1 " m 2  m 1 ‰ m 2  , ( 11
)
and use it to de ne the distance metric as

d(a 1 i , a 2 j ) = h n l n j 1 * J (m(a 1 i ), m(a 2 j )), 1 * J (m(a 1 i ), m(a 2 j )) f ✓ d 0, (m(a 1 i ) " m(a 2 j )) OR (m(a 2 j ) " m(a 1 i )) ÿ, otherwise. , ( 12
)
where ✓ d is a distance threshold, e.g., 0.5 above which two components are considered nonmatching and their distance is set to in nity to prevent false assignments.

After the distance matrix D has been completed, an optimal matching between the components of the two sessions is computed using the Hungarian algorithm to solve the linear assignment problem. As in nite distances are allowed, it is possible to have components from both sessions that are not matched with any other component. This process of registering components across two sessions (REGISTERPAIR) is summarized in Alg. 5.

To register components across multiple sessions, we rst order the sessions chronologically and register session 1 against session 2. From this registration we construct the union of the distinct components between the two sessions by keeping the matched components from session 2 as well as the non-matched components from both sessions aligned to the FOV of session 2. We then register this union of components to the components of session 3 and repeat the procedure until all sessions are have been registered. This process of multi session registration (REGISTERMULTI) is summarized in Alg. 6. At the end of the process the algorithm produces a list of matches between the components of each session and the union of all active distinct components, allowing for e cient tracking of components across multiple days (Fig. 9), and the comparison of non-consecutive sessions through the union without the need of direct pairwise registration (Fig. 10)). An alternative approach to the problem of multiple session registration (CELLREG) was presented recently by Sheintuch et al. (2017) where the authors register neurons across multiple days by rst constructing a similar union set of all the components which is then re ned using a clustering procedure. REGISTERMULTI di ers from the CELLREG method of Sheintuch et al. (2017) in a few key ways, that highlight its simplicity and robustness:

• REGISTERMULTI uses a very simple intersection over union metric to estimate the distance between two neighboring neurons after the FOV alignment. Cells that have a distance above a given threshold are considered di erent by default and are not tested for matching. This parameter is intuitive to set a priori for each dataset. In contrast CELLREG uses a probabilistic framework based on the joint probability distribution between the distance of two cells and the correlation of their shapes that makes speci c parametric assumptions about the distributions of centroid distances between the same and di erent cells, as well as their shape correlations. This model needs to be re-evaluated for every di erent set of sessions to be registered and potentially requires a lot of data to learn the appropriate distance metric.

• REGISTERMULTI uses the Hungarian algorithm to register two di erent set of components, a practice that solves the linear assignment problem optimally under the assumed distance function. In contrast CELLREG uses a greedy method for initializing the assignment of cells to the union superset relying on the following clustering step to re ne these estimates, and thus adding extra computational burden to the registration procedure.

Implementation details for CAIMAN BATCH

Each dataset was processed using the same set of parameters, excepting the expected size of neurons (estimated by inspecting the correlation image), the size of patches and expected number of neurons per patch (estimated by inspecting the correlation image). For the dataset N.01.01, where optical modulation was induced, the threshold for merging neurons was slightly higher (the stimulation caused clustered synchronous activity). For shorter datasets, rigid motion correction was su cient; for longer datasets K53,J115 we applied non-rigid motion correction. The parameters for the automatic selection of components were optimized using only the rst three datasets and xed for all the remaining les. For all datasets the background neuropil activity was modeled as a rank 2 matrix, and calcium dynamics were modeled as a rst order autoregressive process. The remaining parameters were optimized so that all the datasets could be run on a machine with less than 128GB RAM.

Implementation details for CAIMAN ONLINE

Datasets were processed for two epochs with the exception of the longer datasets K53,J115,J123

where only one pass of the data was performed to limit computational cost. For each dataset the online CNN classi er was used to detect new neurons, and ve candidate components were considered for each frame. The online CNN classi er had the same threshold 0.5 for all datasets, with the exception of the longest datasets J115,J123 where the threshold was set to 0.75. Setting the threshold to 0.5 for these datasets led to slightly poorer performance. Large datasets were spatially decimated by a factor of 2 to enhance processing speed, a step that did not lead to changes in detection performance. For all datasets the background neuropil activity was modeled as a rank 2 matrix, and calcium dynamics were modeled as a rst order autoregressive process. For each dataset, CAIMAN ONLINE was initialized on the rst 200 frames, using the BAREINITIALIZATION on the entire FOV with only 2 neurons, so in practice all the neurons were detected during the online mode. To highlight the truly online processing mode, no post-processing of the results was used, a step that can further enhance the performance of the algorithm. Similarly to batch processing, the expected size of neurons was chosen separately for each dataset after inspecting the correlation image.

For the analysis of the whole brain zebra sh dataset, CAIMAN ONLINE was run for 1 epoch with the same parameters as above, with only di erences appearing in the number of neurons during initialization (600 vs 2), and the value of the threshold for the online CNN classi er (0.75 vs 0.5).

The former decision was motivated by the goal of retrieving with a single pass neurons from a preparation with a denser level of activity over a larger FOV in this short dataset (1885 frames).

To this end, the number of candidate neurons at each timestep was set to 10 (per plane). The threshold choice was motivated by the fact that the classi er was trained on mouse data only, and thus a higher threshold choice would help diminish potential false positive components. Rigid motion correction was applied online to each plane.

Performance quanti cation as a function of SNR

To quantify performance as a function of SNR we approximate the ground truth traces by running CAIMAN BATCH on the datasets seeded with the "consensus" binary masks obtained from the manual annotators. After that the average peak-SNR of a trace c with corresponding residual signal r (5) is obtained as

SNR(z) = * *1 (p min ), (13)
where *1 ( ) denotes the probit function (quantile function for the standard Gaussian distribution), z is the z-scored version of c + r (6) and p min is given by (8). After we compute the SNR for both ground truth and inferred traces the performance algorithm can be quanti ed in multiple ways as a function of a SNR thresholds ✓ SNR :

Precision: Precision at level ✓ SNR , can be computed as the fraction of detected components with SNR cm > ✓ SNR that are matched with ground truth components. It quanti es the certainty that a component detected with a given SNR or above corresponds to a true component.

PREC(✓

Additional features of CAIMAN

CAIMAN contains a number of additional features that are not presented in the results section for reasons of brevity. These include:

Volumetric data processing

Apart from planar 2D data, CAIMAN BATCH is also applicable to 3D volumetric data arising either from dense raster scanning methods, or from direct volume imaging methods such as light eld microscopy [START_REF] Prevedel | Simultaneous whole-animal 3D-imaging of neuronal activity using light eld microscopy[END_REF][START_REF] Grosenick | Identi cation Of Cellular-Activity Dynamics Across Large Tissue Volumes In The Mammalian Brain[END_REF].

Segmentation of structural indicator data

Duplicate Detection

The ground truth obtained through the consensus process was screened for possible duplicate selections. To detect for duplicate components we de ne the degree of spatial overlap matrix O as

O ij = h n l n j 0, i = j m(a i
) " m(a j ) m(a j )

, i ë j , ( 14

Extraction of F _F

The uorescence trace f i of component i can be written as

f i = Òa i Ò 2 (c i + r i ). ( 15
)
The uorescence due to the component's transients overlaps with a background uorescence due to baseline uorescence of the component and neuropil activity, that can be expressed as or a for longer traces, a running percentile function, e.g., 10th percentile over a window of a hundred seconds6 . To determine the optimal percentile level an empirical histogram of the trace (or parts of it in case of long traces) is computed using a di usion kernel density estimator [START_REF] Botev | Kernel density estimation via di usion[END_REF], and the mode of this density is used to de ne the baseline and its corresponding percentile level.

f 0,i = BASELINE(f i + B Ò a i ), (16)
The F _F activity of component i can then be written as c

F _F i = f i * BASELINE(f i ) f 0,i (17)
The approach we propose here is conceptually similar to practical approaches where the F _F is computed by averaging over the spatial extent of an ROI [START_REF] Jia | In vivo two-photon imaging of sensory-evoked dendritic calcium signals in cortical neurons[END_REF] with some di erences: i) instead of averaging with a binary mask we use the a weighed average with the shape of each component, ii) signal due to overlapping components is removed from the calculation of the background uorescence, and iii) the traces have been extracted through the CNMF process prior to the F _F extraction. Note that the same approach can also be performed to the trace Òa i Ò 2 c i that does not include the residual traces for each component. In practice it can be bene cial to extract F _F traces prior to deconvolution, since the F _F transformation can alleviate the e ects of drifting baselines, e.g., due to bleaching. For the non-deconvolved traces f i some temporal smoothing can also be applied to obtain more smooth F _F traces.

Algorithm 2 SEEDEDINITIALIZATION

Require: Input data matrix Y , matrix of binary masks M, number of background components n b . 1: p = find(A1 == 0)

. Find the pixels not covered by any component.

;

  Ahrens et al. (2013); Flusberg et al. (2008); Cai et al. (2016); Prevedel et al. (2014); Grosenick et al. (2017); Bouchard et al. (2015)).

Figure 1 .

 1 Figure 1. Processing pipeline of CAIMAN for calcium imaging data. (a) The typical pre-processing steps include (i) correction for motion artifacts, (ii) extraction of the spatial footprints and uorescence traces of the imaged components, and (iii) deconvolution of the neural activity from the uorescence traces. (b) Time average of 2000 frames from a two-photon microscopy dataset (left) and magni ed illustration of three overlapping neurons (right), as detected by the CNMF algorithm. (c) Denoised temporal components of the three neurons in (b) as extracted by CNMF and matched by color (in relative uorescence change, F _F ). (d) Intuitive depiction of CNMF. The algorithm represents the movie as the sum of rank-one spatio-temporal components capturing either neurons and processes, plus additional non-sparse low-rank terms for the background uorescence and neuropil activity. (e) Flow-chart of the CAIMAN BATCH processing pipeline. From left to right: Motion correction and generation of a memory e cient data format. Initial estimate of somatic locations in parallel over FOV patches using CNMF. Re nement and merging of extracted components via seeded CNMF. Removal of low quality components. Final domain dependent processing stages. (f) Flow-chart of the CAIMAN ONLINE algorithm.After a brief mini-batch initialization phase, each frame is processed in a streaming fashion as it becomes available. From left to right: Correction for motion artifacts. Estimate of activity from existing neurons, identi cation and incorporation of new neurons. Periodically, the spatial footprints of inferred neurons are updated (dashed lines).

Figure 2 .

 2 Figure 2. Parallelized processing and component quality assessment for CAIMAN BATCH. (a) Illustration of the parallelization approach used by CAIMAN BATCH for source extraction. The data movie is partitioned into overlapping sub-tensors, each of which is processed in an embarrassingly parallel fashion using CNMF, either on local cores or across several machines in a HPC. The results are then combined. (b) Re nement after combining the results can also be parallelized both in space and in time. Temporal traces of spatially non-overlapping components can be updated in parallel (top) and the contribution of the spatial footprints for each pixel can be computed in parallel (bottom). Parallelization in combination with memory mapping enable large scale processing with moderate computing infrastructure. (c) Quality assessment in space: The spatial footprint of each real component is correlated with the data averaged over time, after removal of all other activity. (d) Quality assessment in time: A high SNR is typically maintained over the course of a calcium transient. (e) CNN based assessment. Top: A 4-layer CNN based classi er is used to classify the spatial footprint of each component into neurons or not. Bottom: Positive and negative examples for the CNN classi er, during training (left) and evaluation (right) phase. The CNN classi er can accurately classify shapes and generalizes across datasets from di erent brain areas. 7 of 40

Figure 3 .

 3 Figure 3. Ground truth generation. (a) Top: Individual manual annotations on the dataset N.04.00.t (only part of the FOV is shown) for labelers L1 (left), L2 (middle), L3(right). Bottom: Disagreements between L1 and L2 (left), and ground truth labels (right) after the consensus between all labelers has been reached. In this example, consensus considerably reduced the number of initially selected neurons. (b) Matches (top) and mismatches (bottom) between each individual labeler and consensus ground truth. Red contours on the mismatches panels denote false negative contours, i.e., components in the consensus not selected by the corresponding labeler, whereas yellow contours indicate false positive contours. Performance of each labeler is given in terms of precision/recall and F 1 score and indicates an unexpected level of variability between individual labelers.

Fig. 3 (

 3 Fig. 3 (bottom) shows an example of matches and mismatches between individual labelers and consensus ground truth for dataset K53, where the level of agreement was relatively high. The high degree of variability in human responses indicates the challenging nature of the source extraction problem and raises reproducibility concerns in studies relying heavily on manual ROI selection.

Figure 4 .

 4 Figure 4. Evaluation of CAIMAN performance against manually annotated data. (a) Comparison of CAIMAN BATCH (top) and CAIMAN ONLINE (bottom) when benchmarked against consensus ground truth for dataset K53. For a portion of the FOV, correlation image overlaid with matches (left panels, true positives red for consensus ground truth, yellow for CAIMAN) and mismatches (right panels, red for false negatives, yellow for false positives). (b) Performance of CAIMAN BATCH, CAIMAN ONLINE and all labelers (L1, L2, L3, L4) for all 9 datasets in terms of F 1 score. CAIMAN BATCH and CAIMAN ONLINE reach near-human accuracy for neuron detection. Complete results with precision and recall for each dataset are given in Table 1. (c-d) Performance of CAIMAN BATCH increases with peak SNR. (c) Example of scatter plot between SNRs of matched traces between CAIMAN BATCH and ground truth for dataset K53. False negative/positive pairs are plotted in green along the x-and y-axes respectively, perturbed as a point cloud to illustrate the density. Most false positive/negative predictions occur at low SNR values. Shaded areas represent thresholds above which components are considered for matching (blue for CAIMAN BATCH selected components and red for GT selected components) (d) F 1 score and upper/lower bounds for all datasets as a function of various peak SNR thresholds. Performance increases signi cantly for neurons with high peak SNR traces (see text for de nition of metrics and the bounds).

Figure 5 .

 5 Figure 5. Evaluation of CAIMAN extracted traces against traces derived from ground truth. (a) Examples of shapes (left) and traces (right) are shown for ve matched components extracted from dataset K53 for consensus ground truth (GT, black), CAIMAN BATCH (yellow) and CAIMAN ONLINE (red) algorithms. The dashed gray portion of the traces is also shown magni ed (bottom-right). Spatial footprints and traces for ground truth are obtained by seeding CAIMAN with the consensus binary masks. The traces extracted from both versions of CAIMAN match closely the ground truth traces. (b) Slope graph for the average correlation coe cient for matches between ground truth and CAIMAN BATCH, and between ground truth and CAIMAN ONLINE. Batch processing produces traces that match more closely the traces extracted from the ground truth data. (c) Empirical cumulative distribution functions of correlation coe cients aggregated over all the tested datasets.Both distributions exhibit a sharp derivative close 1 (last bin), with the batch approach giving better results.

Figure 6 .

 6 Figure 6. Online analysis of a 30 min long whole brain recording of the zebra sh brain. (a) Correlation image overlaid with the spatial components found by the algorithm (portion of plane 11 out of 45 planes in total). (b) Left: Spatial footprints found in the dashed region in (a), contours represent neurons displayed in (c). Right: Correlation image for the same region. (c) Spatial (left) and Temporal (right) components associated to the ten example neurons marked in panel (a). (d) Temporal traces for all the neurons found in the FOV in (a), the initialization on the rst 200 frames contained 500 neurons (present since time 0). (e) Number of neurons found per plane (See also Supplementary Fig. 11 for a summary of the results from all planes).

  Fig.8aalso shows the speed performance of CAIMAN ONLINE (red markers). Because of the low memory requirements of the streaming algorithm, this performance only mildly depends on the computing infrastructure allowing for near real-time processing speeds on a standard laptop (Fig.8a). As discussed in Giovannucci et al.(2017) processing time of CAIMAN ONLINE depends primarily on i) the computational cost of tracking the temporal activity of discovered neurons, ii)

Figure 7 .

 7 Figure 7. Analyzing microendoscopic 1p data with the CNMF-E algorithm using CAIMAN BATCH. (a) Contour plots of all neurons detected by the CNMF-E (white) implementation of Zhou et al. (2018) and CAIMAN BATCH (red) using patches. Colors match the example traces shown in (b), which illustrate the temporal components of 10 example neurons detected by both implementations. CAIMAN BATCH reproduces with reasonable delity the results of Zhou et al. (2018).

Figure 8 .

 8 Figure 8. Time performance of CAIMAN BATCH and CAIMAN ONLINE for the analyzed datasets. (a) Log-log plot of total processing time as a function of data size for CAIMAN BATCH for the 5 largest two-photon datasets using three di erent processing infrastructures: i) a laptop with 8 CPUs (blue), ii) a desktop workstation with 24 CPUs (magenta), and iii) a HPC where 112 CPUs are allocated (yellow). The results indicate a near linear scaling of the processing time with the size of dataset, with additional dependence on the number of found neurons (size of each point). Even very large datasets (> 100GB) can be processed e ciently with a single laptop, whereas access to a HPC enables processing with speed faster than the acquisition time (considered 30Hz for a 512ù512 FOV here). The results of CAIMAN ONLINE using the laptop are also plotted in red indicating near real-time processing speed. (b) Break down of processing time for CAIMAN BATCH (left) and CAIMAN ONLINE (right) (excluding motion correction). (Left) Processing with CNMF in patches and re nement takes most of the time for CAIMAN BATCH. Right: Tracking neural activity and new neuron detection can be done in real-time for CAIMAN ONLINE. (c) (Left) Cost of neural activity online tracking for the whole brain zebra sh dataset (maximum time over all planes per frame). Tracking can be done in real-time. (Right) The most expensive part during online processing occurs while updating the spatial footprints, a step that can be distributed or parallelized. Each color corresponds to the update cost for the various di erent planes. (d) Cost analysis of CNMF-E implementation for processing a 6000 frames long 1p dataset. Processing in patches in parallel induces a time/memory tradeo and can lead to speed gains (patch size in legend).

Figure 9 .

 9 Figure 9. Components registered across six di erent sessions (days). (a) Contour plots of neurons that were detected to be active in all six imaging sessions overlaid on the correlation image of the sixth imaging session. Each color corresponds to a di erent session. (b) Stability of multiday registration method. Comparisons of forward and backward registrations in terms of F 1 scores for all possible subsets of sessions. The comparisons agree to a very high level indicating the stability of the proposed approach. (c) Comparison (in terms of F 1 score) of pair-wise alignments using readouts from the union vs direct alignment. The comparison is performed for both the forward and the backwards alignment. For all pairs of sessions the alignment using the proposed method gives very similar results compared to direct pairwise alignment. (d) Magni ed version of the tracked neurons corresponding to the squares marked in panel (a). Neurons in di erent parts of the FOV exhibit di erent shift patterns over of the course of multiple days, but can nevertheless be tracked accurately by the proposed multiday registration method.

  being the temporal trace of component i. B is the background/neuropil activity matrix. For two-photon data it is modeled as a low rank matrix B = bf, where b À R dùn b , f À R n b ùT correspond to the matrices of spatial and temporal background components, and n b is the number of background components. For the case of micro-endoscopic data the integration volume is much larger and the low rank model is inadequate. A solution comes from the CNMF-E algorithm of Zhou et al. (2018) where the background is modeled as

  i denote the spatial footprint and temporal trace of component i, and the let A \i , C \i denote the matrices A, C when the component i has been removed. Similarly, let Y i = Y * A \i C \i * B denote the entire dataset when the background and the contribution of all components except i have been removed. If component i is real then Y i and a i c Ò i will look similar during the time intervals when the component i is active. As a rst test CAIMAN nds the rst N p local peaks of c i (e.g., N p = 5),

  CAIMAN BATCH classi er for post processing classi cationThe purpose of the batch classi er is to classify the components detected by CAIMAN BATCH into neuron somas or other shapes, by examining their spatial footprints. Only three annotated datasets (.03.00.t, NF.04.00.t, NF.02.00) were used to train the batch classi er. The set of estimated footprints from running CAIMAN BATCH initialized with the consensus ground truth was matched to the set of ground truth footprints. Footprints matched to ground truth components were considered positive examples, whereas the remaining components were labeled as negatives. The two sets were enriched using data augmentation (rotations, re ections, contrast manipulation etc.) through the Keras library (keras.io) and the CNN was trained on 60% of the data, leaving 20% for validation and 20% for testing. The CNN classi er reached an accuracy of 97% on test data; that also generalized to the rest of the datasets (Fig.2e) without any parameter change.Online classi er for new component detectionThe purpose of the CAIMAN ONLINE classi er is to detect new components based on their spatial footprints by looking at the mean across time of the residual bu er. To construct the ground truth data for the online classi er, CAIMAN BATCH was run on the rst ve annotated datasets seeded with the masks obtained through the manual annotations. Subsequently the activity of random subsets of found components and the background was removed from contiguous frames of the raw datasets to construct residual bu ers, which were averaged across time. From the resulting images patches were extracted corresponding to positive examples (patches around a neuron that was active during the bu er) and negative examples (patches around other positions within the FOV). A neuron was considered active if its trace attained an average peak-SNR value of 4 or higher during the bu er interval. Similarly to the batch classi er, the two sets were augmented and split into training, validation and testing sets. The resulting classi er reached a 98% accuracy on the testing set, and also generalized well when applied to di erent datasets.Di erences between the two classi ersAlthough both classi ers examine the spatial footprints of candidate components, their required performance characteristics are di erent which led us to train them separately. The batch classi er examines each component as a post-processing step to determine whether its shape corresponds to a neural cell body. As such, false positive and false negative examples are treated equally and possible mis-classi cations do not directly a ect the traces of the other components. By contrast, the online classi er operates as part of the online processing pipeline. In this case, a new component that is not detected in a residual bu er is likely to be detected later should it become more active. On the other hand, a component that is falsely detected and incorporated in the online processing pipeline will continue to a ect the future bu er residuals and the detection of future components. As such the online algorithm is more sensitive to false positives than false negatives. To ensure a small number of false positive examples under testing conditions, only components with average peak-SNR value at least 4 were considered as positive examples during training of the online classi er.

N

  be the ground truth traces and c cm 1 , c cm 2 , … , c cm N be their corresponding CAIMAN inferred traces. Here we assume that false positive and false negative components are matched with trivial components that have 0 SNR. Let also SNR gt i = SNR(c gt i ) and SNR cm i = SNR(c cm i ), respectively.

  Structural indicators expressed in the nucleus and functional indicators expressed in the cytoplasm can facilitate source extraction and help identify silent or speci c subpopulations of neurons (e.g., inhibitory). CAIMAN provides a simple adaptive thresholding ltering method for segmenting summary images of the structural channel (e.g., mean image). The obtained results can be used for seeding source extraction from the functional channel in CAIMAN BATCH or CAIMAN ONLINE as already discussed.

  ) that de nes the fraction of component i that overlap with component j, where m( ) is the thresholding function de ned in(10). Any entry of O that is above a threshold ✓ o (e.g., ✓ o = 0.7 used here) indicates a pair of duplicate components. To decide which of the two components should be removed, we use predictions of the CAIMAN BATCH CNN classi er, removing the component with the lowest score.

  where BASELINE : R T ≠ R T is a baseline extraction function, and B is the estimated background signal. Examples of the baseline extraction function are a percentile function (e.g., 10th percentile),

A

  2: [Ì, f ] } NNMF(Y [p, :], n b ) . Run NMF on these pixels just to get temporal backgrounds f 3: b } arg min bg0 ÒY * bf Ò . Obtain spatial background b. 4: C } max (M Ò M) *1 M Ò (Y * bf ), 0 . Initialize temporal traces. 5: A } arg min Ag0,A(ÌM)==0 ÒY * bf * ACÒ. . Initialize spatial footprints constrained within the masks. 6: repeat . Optionally keep updating A, C, b, f using HALS. } arg min Ag0,A(ÌM)==0 ÒY * bf * ACÒ 10: until Convergence 11: return A, C, b, f Algorithm 3 CAIMAN ONLINE (See Giovannucci et al. (2017) for explanation of routines) Require: Data matrix Y , initial estimates A, b, C, f, S, current number of components K, current timestep t ® , rest of parameters. 1: W = Y [:, 1 : t ® ]C Ò _t ® 2: M = CC Ò _t ® . Initialize su cient statistics (Giovannucci et al., 2017) 3: G = DETERMINEGROUPS([A, b], K) . Giovannucci et al. (2017), Alg. S1-S2 4: R buf = [Y * [A, b][C; f]][:, t ® * l b + 1 :

  for 23: return A, b, C, f, S

Table 1 .

 1 Results of each labeler, CAIMAN BATCH and CAIMAN ONLINE algorithms against consensus ground truth. Results are given in the form F 1 score (precision, recall), and empty entries correspond to datasets not manually annotated by the speci c labeler. In italics the datasets used to train the CNN classi ers.

	Name	L1	L2	L3	L4	CAIMAN BATCH	CAIMAN ONLINE
	N.03.00.t	X	0.90	0.85	0.78	0.78	0.76
			(0.88,0.92)	(0.78,0.93)	(0.73,0.83)	(0.77,0.79)	(0.77,0.75)
	N.04.00.t	X	0.69	0.75	0.87	0.67	0.68
			(0.54,0.97)	(0.61,0.97)	(0.78,0.98)	(0.62,0.72)	(0.65,0.71)
	N.02.00	0.89	0.87	0.84	0.82	0.79	0.77
		(0.86,0.93)	(0.88,0.85)	(0.92,0.77)	(1.00,0.70)	(0.8,0.77)	(0.79,0.76)
	N.00.00	X	0.92	0.83	0.87	0.72	0.72
			(0.93,0.91)	(0.86,0.80)	(0.96,0.80)	(0.83,0.64)	(0.83,0.64)
	N.01.01	0.80	0.89	0.78	0.75	0.77	0.73
		(0.95,0.69)	(0.96,0.83)	(0.73,0.84)	(0.80,0.70)	(0.88,0.69)	(0.78,0.68)
	YST	0.78	0.90	0.82	0.79	0.76	0.78
		(0.76,0.81)	(0.85,0.97)	(0.75,0.92)	(0.96,0.67)	(0.9,0.66)	(0.76,0.81)
	K53	0.89	0.92	0.93	0.83	0.77	0.82
		(0.96,0.83)	(0.92,0.92)	(0.95,0.91)	(1.00,0.72)	(0.83,0.72)	(0.80,0.83)
	J115	X	0.93	0.94	0.83	0.77	0.81
			(0.94,0.91)	(0.95,0.93)	(1.00,0.71)	(0.9,0.68)	(0.75,0.88)
	J123	0.85	0.83	0.90	0.91	0.68	0.80
		(0.96,0.76)	(0.73,0.96)	(0.91,0.90)	(0.92,0.89)	(0.94,0.51)	(0.82,0.79)

  aries. Neurons lying close to the boundary between neighboring patches can appear multiple times and must be merged. With this goal, we optimized the merging approach used in Pnevmatikakis

et al.

(2016)

: Groups of components with spatially overlapping footprints whose temporal traces are correlated above a threshold are replaced with a single component, that tries to explain as much of the variance already explained by the "local" components (as opposed to the variance of the data as performed in

[START_REF] Pnevmatikakis | Simultaneous Denoising, Deconvolution, and Demixing of Calcium Imaging Data[END_REF]

). If A old , C old are the matrices of components to be merged, then the merged component a m , c m are given by the solution of the rank-1 NMF problem:

Table 2 .

 2 Properties of manually annotated datasets. For each dataset the duration, imaging rate and calcium indicator are given, as well as the number of active neurons selected after consensus of the manual annotations.

	Name	Area brain	Lab	Rate (Hz) Size (TùXùY)	Indicator # labelers # neurons GT
	NF.03.00.t Hippocampus Losonczy	7	2250x498x467 GCaMP6f	3	178
	NF.04.00.t	Cortex	Harvey	7	3000x512x512 GCaMP6s	3	257
	NF.02.00	Cortex	Svoboda	30	8000x512x512 GCaMP6s	4	394
	NF.00.00	Cortex	Svoboda	7	2936x512x512 GCaMP6s	3	425
	NF.01.01 Visual Cortex Hausser	7	1825x512x512 GCaMP6s	4	333
	YST	Visual Cortex	Yuste	10	3000x200x256 GCaMP3	4	405
	K53	Parietal Cortex Tank	30	116043x512x512 GCaMP6f	4	920
	J115	Hippocampus	Tank	30	90000x463x472 GCaMP5	3	891
	J123	Hippocampus	Tank	30	41000x458x477 GCaMP5	4	183

  Recall at level ✓ SNR , can be computed as the fraction of ground truth components with SNR gt > ✓ SNR that are detected by the algorithm. It quanti es the certainty that a ground truth component with a given SNR or above is detected. The cautious reader will observe that the precision and recall quantities described above are not computed in the same set of components. This can be remedied by recomputing the quantities in two di erent ways:AND framework: Here we consider a match only if both traces have SNR above the given threshold: Here we consider a match if either trace has SNR above the given threshold and its match has SNR above 0. As demonstrated in Fig.4d, these bounds are tight.

						SNR ) =	{i : (SNR cm {i : (SNR cm i > ✓ SNR ) & (SNR gt i > ✓ SNR )}	i > 0)}
	Recall: RECALL(✓	SNR ) =	{i : (SNR gt {i : (SNR gt i > ✓ SNR ) & (SNR cm i > ✓ SNR )}	i > 0)}
	F 1 score: An overall F 1 score at level ✓	SNR , can be obtained by computing the harmonic mean
	between precision and recall
				F 1 (✓	SNR ) = 2	PREC(✓ PREC(✓ SNR ) + RECALL(✓ SNR ) ù RECALL(✓ SNR ) SNR )
	PREC	AND (✓ SNR ) =	{i : (SNR cm {i : (SNR cm i > ✓ SNR ) & (SNR gt i > ✓ SNR )} i > ✓ SNR )}
	RECALL	AND (✓	SNR ) =	{i : (SNR gt {i : (SNR gt i > ✓ SNR ) & (SNR cm i > ✓ SNR )} i > ✓ SNR )}
	OR framework: RECALL OR (✓	SNR ) =	{i : (max(SNR gt	i , SNR cm {i : (SNR cm i ) > ✓ SNR ) & (min(SNR gt i > 0)}	i , SNR cm	i ) > 0)}
	RECALL OR (✓	SNR ) =	{i : (max(SNR gt	i , SNR cm {i : (SNR gt i ) > ✓ SNR ) & (min(SNR gt i > 0)}	i , SNR cm	i ) > 0)}
	It is easy to show that
		PREC	AND (✓	SNR ) f PREC(✓ SNR ) f PREC OR (✓	SNR )
	RECALL AND (✓	SNR ) f RECALL(✓ SNR ) f RECALL OR (✓ SNR )
				F	1AND (✓ SNR ) f F 1 (✓ SNR ) f F 1OR (✓	SNR ),
	with equality holding for ✓ SNR = 0.

and Giovannucci, 2017)

  t ® ]. Initialize residual bu er 5: t = t ® 6: for i = 1, … , N epochs do FINDNEWCOMPONENTS(R buf , N comp ) .

	7:	while there is more data do
	8:	t } t + 1
	9:	y
	10:	[c
		Alg. 4
	13:

t } MOTIONCORRECT(y t , bf t*1 ) . (Pnevmatikakis t ; f t ] } UPDATETRACES([A, b], [c t*1 ; f t*1 ], y t , G) . Giovannucci et al. (2017), Alg. S3 11: C, S } OASIS(C, , s min , ) . Friedrich et al. (2017b) 12:

A new , C new } [A, b], [C, f ], K, G, R buf , W , M } INTEGRATENEWCOMPONENTS( 14: [A, b], [C, f], K, G, A new , C new , R buf , y t , W , M) . Giovannucci

et al. (2017), Alg. S4 15: R buf } [R buf [:, 2 : l b ], y t * Ac t * bf t ] . Update residual bu er 16: W , M } UPDATESUFFSTATISTICS(W , M, y t , [c t ; f t ])

Since proteins expressing the calcium indicator are con ned outside the cell nuclei, neurons will appear as ring shapes, with a dark disk in the center.

It is possible that this process generated slightly biased results in favor of each individual annotators since the ground truth was always a subset of the union of the individual annotations.

These precision and recall metrics are computed on di erent sets of neurons, and therefore strictly speaking one cannot combine them to form an F 1 score. However, they can be bound from above by being evaluated on the set of matched and non-matched components where at least one trace is above the threshold (union of blue and pink zones in Fig.4c) or below by considering only matched and non-matched components where both ground truth and inferred traces have SNR above the threshold (intersection of blue and pink zones in Fig.4c). In practice these bounds were very tight for all but one dataset (Fig.4d). More details can be found in Methods and Materials (Performance quanti cation as a function of SNR).

Computing the exact running percentile function can be computationally intensive. To reduce the complexity we compute the running percentile with a stride of W , where W is equal or smaller to the length of the window, and then linearly interpolate the values.

Acknowledgments

We thank B. Cohen, L. Myers, N. Roumelioti, and S. Villani for providing us with manual annotations.

We thank V. Staneva and B. Deverett for contributing to the early stages of CAIMAN, and L. Paninski for numerous useful discussions. We thank N. Carriero, I. Fisk, and D. Simon from the Flatiron Institute (Simons Foundation) for useful discussions and suggestions to optimize High Performance Computing code. We thank T. Kawashima and M. Ahrens for sharing the whole brain zebra sh dataset. Last but not least, we thank the active community of CAIMAN users for their great help in terms of code/method contributions, bug reporting, code testing and suggestions that have led to the growth of CAIMAN into a widely used open source package. A partial list of contributors (in the form of GitHub usernames) can be found in https://github.com/ atironinstitute/CaImAn/ graphs/contributors (Python) and https://github.com/ atironinstitute/CaImAn-MATLAB/graphs/ contributors (MATLAB ® ). The authors acknowledge support from following funding sources: A.G., E.A.P., J.F., P.G. (Simons Foundation). J.G., S.A.K., D.W.T. (NIH NRSA F32NS077840-01,5U01NS090541, 1U19NS104648 and Simons Foundation SCGB), P.Z. (NIH NIBIB R01EB022913, NSF NeuroNex DBI-1707398, Gatsby Foundation), J.T. (NIH R01-MH101198), F.N. (MURI, Simons Collaboration on the Global Brain and Pew Foundation).

Appendix 0 Figure 10. Tracking neurons across days, step-by-step description of multi session registration (Fig. 9). (a) Correlation image overlaid to contour plots of the neurons identi ed by CAIMAN BATCH in day 1 (top, 453 neurons), 2 (middle, 393 neurons) and 3 (bottom, 375 neurons). (b) Result of the pairwise registration between session 1 and 2. The union of distinct active components consists of the components that were active in i) both sessions (yellow -where only the components of session 2 are displayed), ii) only in session 2 (green), and iii) only in session 1, aligned to the FOV of session 2 (red). (c) At the next step the union of sessions 1 and 2 is registered with the results of session 3 to produce the union of all distinct components aligned to the FOV of session 3. (d) Comparison of non-consecutive sessions without pairwise registration. Keeping track of which session each component was active in, enables e cient and stable comparisons.

Algorithmic Details

In the following we present in pseudocode form several of routines introduced and used by CAIMAN.

Note that the pseudocode descriptions do not aim to present a complete picture and may refer to other work for some of the steps.

Algorithm 1 PROCESSINPATCHES

Require: Input data matrix Y , patch size , overlap size, initialization method, rest of parameters.

. Break data into memory mapped patches. 2: for i = 1, … , N p do . Process each patch 3:

. Run CNMF on each patch 4: end for 5: [A, C] = MERGECOMPONENTS[{A (i) ,

. Merge background components 7: M } (A > 0).

. Find masks of spatial footprints. 8: repeat . Optionally keep updating A, C, b, f using HALS [START_REF] Cichocki | Hierarchical ALS algorithms for nonnegative matrix and 3D tensor factorization[END_REF].

9:

. Compute motion eld between the templates.

. Turn components into binary masks.

. Compute distance matrix.

. Register A u to session i.