Tracking calcium dynamics from individual neurons in behaving animals.

Measuring the activity of neuronal populations with calcium imaging can capture emergent functional properties of neuronal circuits with single cell resolution. However, the motion of freely behaving animals, together with the intermittent detectability of calcium sensors, can hinder automatic monitoring of neuronal activity and their subsequent functional characterization. We report the development and open-source implementation of a multi-step cellular tracking algorithm (Elastic Motion Correction and Concatenation or EMC2) that compensates for the intermittent disappearance of moving neurons by integrating local deformation information from detectable neurons. We demonstrate the accuracy and versatility of our algorithm using calcium imaging data from two-photon volumetric microscopy in visual cortex of awake mice, and from confocal microscopy in behaving Hydra, which experiences major body deformation during its contractions. We quantify the performance of our algorithm using ground truth manual tracking of neurons, along with synthetic time-lapse sequences, covering a wide range of particle motions and detectability parameters. As a demonstration of the utility of the algorithm, we monitor for several days calcium activity of the same neurons in layer 2/3 of mouse visual cortex in vivo, finding significant turnover within the active neurons across days, with only few neurons that remained active across days. Also, combining automatic tracking of single neuron activity with statistical clustering, we characterize and map neuronal ensembles in behaving Hydra, finding three major non-overlapping ensembles of neurons (CB, RP1 and RP2) whose activity correlates with contractions and elongations. Our results show that the EMC2 algorithm can be used as a robust and versatile platform for neuronal tracking in behaving animals.

This is a PLOS Computational Biology Methods paper.

Introduction

Measuring the activity of neuronal populations in freely behaving animals can help a detailed understanding of how neural circuits integrate external information, compute, learn and control animal behavior. Calcium imaging has become widespread for measuring single neuron activity as it is non-invasive and allows the simultaneous measurement of hundreds to thousands of cells, with single cell resolution [1]. Moreover, monitoring single neuron activity in freely moving animals such as rodents can be achieved with miniaturized microscopes attached to the head [2]. However, technical limitations of current microscopy techniques and of mathematical analysis hinder a more complete imaging and analysis of entire brains. Alternative strategies consist of monitoring single neuron activity in targeted brain regions of rodents using two-photon microscopy [3], or imaging the nervous system of a smaller animal, one that can fit entirely within a microscope’s field of view, such as Caenorhabditis elegans [4], Hydra [5], Zebrafish [6] or Drosophila larvae [7]. An advantage of simple model organisms is that they contain many fewer neurons than mammals and have a limited repertoire of behaviors [8] that may be entirely characterized in the near future.

Aside from the difficulties in imaging entire nervous systems with high temporal and spatial resolution [3, 6, 9], an important bottleneck in analyzing calcium imaging data is to achieve robust and automatic tracking of individual cells’ position while the animal is behaving. Single-cell tracking is challenging for three main reasons: First, there may be a large number of cells in a cluttered environment. Therefore, false positives and negatives during single cell detection and localization impede the association of detected neurons between successive time frames and call for more elaborate tracking algorithms. Second, neurons can remain undetectable over large periods of time because calcium sensors may be significantly brighter than background only when neurons are firing. Third, tracking methods in behaving animals have to handle animal motion and body deformations [10]. To tackle these issues, numerous hardware solutions have been proposed, including fixation of animal (usually the head) [3, 6, 11], high-speed motorization of microscopes to track animal movements [4, 12] and dual-color labeling of neurons with a calcium-insensitive probe that can be detected and localized even when neurons are not firing [4, 9]. However, even when (some of) these solutions are implemented, the residual motion of neurons, the limited spatial resolution of the microscope and the intermittency of the fluorescence signal (when single-channel calcium imaging is used) hinder the robust tracking of single neuron activity, particularly over long periods of time. Therefore, elaborate post-processing of acquired movies is required. Most current image processing methods consist of registering images (volumes) to reference image(s) (volume(s)) using either the local fluorescence intensity of images [13, 14], or the extracted neuron positions directly [4, 10, 15, 16]. Then, additional image processing, such as non-negative matrix factorization [17, 18], might be required for demixing and denoising cellular calcium activity.

Despite these software developments, the difficult implementation of hardware solutions, such as the dual-labeling and imaging of calcium-insensitive probes, together with the significant deformability of challenging animal models such as Hydra, prevent the robust automatic tracking of single neurons in many experiments. Tracking has then to be performed manually [5] or semi-manually [19]. This limits the analysis to a few hundred frames, introduces operator bias and, ultimately, hinders our understanding of the functional organization of nervous systems.

To robustly track particles with intermittent detectability (neurons) in a cluttered and deforming environment, we report the development of an algorithm named Elastic Motion Correction and Concatenation (EMC2). EMC2 is based on the versatile framework of single-particle-tracking (SPT), enabling the robust monitoring of single neuron activity in most, if not all, animal models. In contrast to traditional SPT, EMC2 does not set expected priors for particle motion (typically diffusion and/or linear motion). Instead, it uses information about local motion and deformation from detectable and tracked particles in the neighborhood of undetectable particles. EMC2 is therefore more versatile, and does not require motion priors or heuristics to close tracking gaps. In addition, for the local tracking of detectable particles, EMC2 uses a probabilistic method and is therefore robust to very cluttered conditions. We validate the robustness and accuracy of EMC2 with manual tracking of neurons in two calcium imaging datasets from behaving animals. Our first dataset consists of two-photon calcium imaging of few tens of neurons in the visual cortex of awake mice. We show that our algorithm accurately tracks the limited motion of single neurons in the two-photon field-of-view, enabling the fast analysis of long recordings of individual neuron activity. We then monitor single neurons of Hydra’s nervous system while the animal is behaving and deforming. Hydra imaging datasets represent perhaps the worst possible scenario for tracking purposes, since animals can have major changes in body size with non-isometric deformations. We also quantify the performance of EMC2 using simulations of fluorescence time-lapse sequences with different types of motion (confined diffusion, linear displacement and elastic deformation), and show that EMC2 outperforms state-of-the art tracking algorithms.

After integrating EMC2 in an open-source and freely available platform Icy [20] (icy.bioimageanalysis.org), we explore the utility of the algorithm in two experimental scenarios. We first monitor the activity of single neurons in two-photon calcium imaging of layer 2/3 of mouse visual cortex over multiple days and show that there is an important turnover of active neurons and that very few neurons remain active across days. We then track complete neuron activity in behaving Hydra, and find functional clustering of individual neurons into co-active ensembles [21]. Consistent with previous observations [5], we find that Hydra contains three main neuronal ensembles (CB, RP1 and RP2); and, after mapping the positions of individual neurons from each ensemble, we also confirm that these ensembles are not overlapping, i.e., they do not share neurons, and that they are correlated with contraction bursts and elongation behaviors.

These results demonstrate that EMC2 is an effective and versatile tracking algorithm for the tracking of single neuron activity in calcium imaging of living animals. Robust tracking constitutes a prerequisite for the statistical analysis of the functional organization of neural circuits, the description of emergent computational units such as neuronal ensembles, and ultimately, for the understanding and prediction of animals’ adaptive behavior.

Results

1 Statistical mapping of neuron positions versus single-particle-tracking

Most of the methods that have been developed for tracking single neuron activity in the well-documented animal model C. Elegans are based on the statistical mapping of neuron positions to a reference set of coordinates within the animal [10, 16, 22]. An important advantage of these methods is their robustness to the length of the analyzed time-lapse sequence, as the different images are registered independently from each other to the reference set of positions. In these methods, reference positions are extracted from reference frames and cell identities can be obtained either from stereotyped fluorescent color maps of all neurons (NeuroPAL [23]) or online cell atlases (e.g. WormAtlas [24] and OpenWorm [25]). Therefore, mapping methods heavily rely on the stable repertoire of neurons (position and/or type) within single worms, when mapping is used for tracking neurons along time [10, 18, 26], and even across different worms when different animals or strains are compared [16, 22].

Unfortunately, in Hydra, the number and position of neurons differ from one animal to another [27]. Moreover, even within a single animal, the accurate mapping of a reference set of neuron positions is prevented by the important and continual deformations of the animal (typically, the length of the animal is reduced by more than half during longitudinal contraction), and the 2D imaging of transparent 3D tissues that induces apparent changes in neuron positions even between stereotypical poses of the animal. Finally, while the aforementioned mapping methods used in C. Elegans can accommodate moderate changes in the total number of neurons between frames due to missing or spurious detections (counting noise), the intermittent and sparse detectability of neurons in calcium imaging, without a reference fluorophore like calcium-insensitive red fluorescent proteins (RFP) [10], definitely hinders the applicability of mapping methods.

Therefore, to track neurons in Hydra, we used the SPT framework. SPT methods sequentially link cell detections through time [8], with no need for a reference set of positions. They are more versatile than statistical mapping, and can handle the large deformations of Hydra and the intermittent detectability of cells in calcium imaging.

2 Limitations of standard SPT Algorithms

Most SPT algorithms rely on the automatic detection of particles (cells, molecules…) that are significantly brighter than the noisy background in each frame of the time-lapse sequence (see [28] and [29] for review) and, subsequently, the linking of detections between frames corresponding to the reconstruction of coherent particle trajectories (see Table 1). The prevalence of false positives (i.e. background signal) and negatives (i.e. missing detections) in the detection of particles, together with the influence of high particle density and stochastic dynamics have, over the last two decades, motivated the development of algorithms that go beyond naïve tracking methods that simply associate nearest-neighbor detections between consecutive time frames (see Table 1 and [30, 31] for a review of existing methods). Indeed, the erroneous association of one detection to a track, or the premature ending of a track due to missing detection(s), can lead to important error propagation as detections are sequentially associated to existing tracks.

Table 1

Tracking methods in bio-imaging.

Algorithm	Type	Detection	Linking	Gap closing	Pros	Cons	Freely available	Ref.
Sage et al.	Global	None	Energy minimization	Yes	Global.	Designed for one or few sparse particles	ImageJ plugin	[53]
Bonneau et al.	Global	None	Energy minimization	Yes	Global. Robust gap closing with minimal-path algorithm	Designed for few sparse particles. High computational load.	No	[54]
NeRVE	Detect & Mapping	Watershed Segmentation	Point-set registration & clustering–Animal deformation estimated with elastic transformations	Yes	Robust to dense packing of particles. Handles non-linear deformations.Time-independent (i.e. robust even in long time-lapse sequences)	High computational load. Not robust to many missing detections and long gaps in highly deforming environments.	Matlab GUI	[10]
fDLC	Detect & Mapping	Watershed Segmentation	Point-set registration to reference set of positions–learning of animal deformation	Yes	Robust to dense packing of particles. Handles non-linear deformations.Time-independent	Not robust to many missing detections and long gaps in highly deforming environments.	Python (Github repository)	[16]
CRF_ID	Detect & Mapping	Gaussian mixture model fitting	Point-set registration to reference & temporally-nearby frames–Use of graphical model (neighbors) to predict identities	Yes	Robust to dense packing of particles. Handles non-linear deformations.Time-independent	High computational load. Not robust to many missing detections and long gaps in highly deforming environments.	Matlab (Github repository)	[22]
Mosaic	Detect & Link	Gaussian Convolution & Thresholding	Global distance minimization	Yes	Fast. Accounts for spot intensity and size in distance computation.	Gap closing does not handle large motion. Not robust in very cluttered conditions.	Particle Tracker plugin (ImageJ)	[32]
TrackMate	Detect & Link	Wavelet transformation or Gaussian convolution & thresholding	Global distance minimization	Yes	Fast. Handles split & merge events.	Gap closing does not handle large, non-linear motion. Many user-defined parameters	TrackMate plugin (ImageJ)	[33]
eMHT	Detect & Link	Wavelet transformation & thresholding	Probabilistic (Multiple Hypothesis)	Yes	Robust in cluttered environment. Few user-defined parameters	Slower than global distance minimization. Cannot close large gaps (> ~5 frames) due to computational load	Spot Tracking plugin (Icy)	[30]
MAP-4D-DAE	Detect & Link	Not specified	Probabilistic (Multiple Hypothesis) + Autoencoding for particle motion modeling	Yes	Robust in cluttered environment. Few user-defined parameters. Handles non-linear deformations	Slower than global distance minimization. Cannot close large gaps (> ~5 frames) due to computational load	No	[10]

One category of elaborated tracking algorithms is based on global distance minimization (GDM) between all pairs of detections in consecutive time frames. The distance measure between detections can be simply the Euclidean distance, or can additionally take into account the similarity of the intensity and/or shape of the detected particle [32]. To handle possible missing and false detections, heuristics for track termination and initiation are defined by the user [33]. GDM methods are fast and robust, but user-defined parameters that regulate the ending, the initiation and the fusion of tracks hinder their applicability in very cluttered conditions with numerous missing and spurious detections [30]. Moreover, the limitations of current particle motion models (confined diffusion and linear directed motion [32-34]) prevent the robust estimation of particle positions when they remain undetectable over long periods of time, as for sparsely firing neurons in calcium imaging. This limits the ability to close gaps in the trajectories of tracked particles to a very few frames.

Probabilistic methods are an attractive alternative to GDM methods, even if their computational cost is higher. Probabilistic algorithms model both the stochastic motion of the particles and their detectability, and then compute the optimal tracking solution by maximizing the model likelihood of observed detections [30, 35, 36]. The gold standard of probabilistic association is multiple hypothesis tracking (MHT) in which one computes all the possible tracking solutions over the entire time-lapse sequence, before inferring the tracking solution that maximizes the observation’s likelihood. However, MHT is generally not computationally tractable. Approximate solutions that iteratively compute a nearly optimal solution over a limited number of frames (typically up to 5) have been proposed [30]. As probabilistic methods model the particles’ detectability, they are usually more robust than GDM methods in cluttered conditions. However, the small number of frames considered when approximating the MHT solution, together with limited particle motion models (again diffusion and/or linear displacements [30, 35]), reduce the capability of probabilistic algorithms to close large tracking gaps and keep track of particles’ putative position when they remain undetectable over many frames.

3 EMC2 Algorithm

To increase the capability of SPT algorithms to track single particles with intermittent detectability in cluttered and deforming environments, we have developed EMC2. This multi-step algorithm and software is particularly well-suited for tracking single neuron activity with calcium imaging in a behaving animal. The EMC2 algorithm can be decomposed into four main steps (Fig 1).

Fig 1

Multi-step EMC2 for tracking neuron activity in calcium imaging data.

a- Time-lapse imaging (N frames) of intermittent fluorescence activity of a neuron in a deforming environment (e.g. behaving animal). b- Fluorescent spots (neurons), that are significantly brighter than background, are automatically detected with a wavelet-based algorithm. c- Tracklets of detectable neurons are robustly reconstructed using probabilistic tracking algorithm (eMHT). d- Short tracklets of detectable particles are used to compute the elastic deformation of the field of view at each time frame. Associated detections in neuron tracklets are used as fiducials, and the whole deformation is interpolated using a poly-harmonic thin-plate spline function. Forward- and backward-propagated positions of tracklet particle positions are shown with a thin blue line. e- After having corrected for the deformation of the field-of-view where neurons are embedded, gaps between the end- and starting-points of tracklets are closed by minimizing the global Euclidean distance between points (dotted line). f- Finally, complete single neuron tracks over the time-lapse sequence are obtained by applying the elastic transformation of the field-of-view to concatenated tracklets.

First, bright spots (e.g. firing neurons in calcium imaging sequence) are automatically detected with a robust method based on wavelet decomposition of time-lapse sequence and statistical thresholding of wavelet coefficients (Materials and Methods). Second, detected spots are linked into single particle trajectories with a state-of-the-art probabilistic algorithm, a variant of multiple hypothesis tracking (eMHT [30]), which is particularly robust in cluttered conditions. Obtained track(let)s correspond to trajectories of detectable particles. However, in many time-lapse sequences such as calcium imaging of neuron activity, tracks would be terminated prematurely when particles switch to an undetectable state (e.g. non-firing neuron) and new tracks would be generated when particles can be detected again (e.g. firing neurons). This creates time-gaps in individual tracks that need to be closed to allow the accurate tracking of each particle’s identity over the time-lapse sequence. Thus, the two last steps of our method aim to close gaps in trajectories using information about the motion and deformation of the field of view along the time-lapse sequence. We considered that tracked particles are embedded in a deformable medium (e.g. neurons within tissue) and that local estimation of the deformation of the field-of-view should allow the inference of particles’ positions even when they are undetectable.

The third step of EMC2 is therefore the computation of the elastic deformation of the field-of-view at each time using the information contained in tracklets of detectable particles. For this, we used the positions of tracked particles between consecutive time frames as fiducial source and target points. We then computed the forward and backward elastic deformations of the whole field by interpolating the deformation at any position between fiducials with a thin-plate-spline function. The thin-plate spline is a popular poly-harmonic spline whose robustness in image alignment and point-set matching has been demonstrated [37], and which has been recently applied for automatic neuron registration in time-lapse sequences [10]. In our hands, the Neuron Registration Vector Encoding (NeRVE) method developed to map single neurons in C. elegans [10] is unfortunately not sufficiently robust for tracking neurons in calcium imaging of behaving Hydra, as shown in [15]. In addition to the inherent limitations of mapping methods when applied to calcium imaging of Hydra (see paragraph “Statistical Mapping of Neuron Positions versus Single-Particle-Tracking”), the poor tracking accuracy of the NeRVE method when applied to Hydra is also due to the fact that, contrary to C. elegans, one cannot map the Hydra neuron positions in fixed cylindrical coordinates along the principal axis of the animal. This reflects a much lower level of effective deformation in C. elegans and consequently a lower complexity of neuron registration. To increase the robustness and accuracy of single neuron tracking, we rather implement, in EMC2, a local concatenation of short tracklets after having propagated forward and backward the estimated deformation of the field of view.

Therefore, the fourth and last step of our method is the iterative estimation and correction of the elastic deformation of the field-of-view, followed by the optimal concatenation of tracklets. In this last step, we used the elastic transformation computed with thin-plate splines to propagate forward (and backward) the putative positions of undetectable particles, following the termination of their detectable tracklets (or preceding the initiation of novel tracklets). After having corrected for the elastic deformation of the field-of-view, we then linked short tracklets by minimizing the global distance between the end-points of prematurely terminated tracklets with the starting-points of newly appearing tracklets (Materials and Methods). Finally, single-particle tracks over the time-lapse sequence are obtained by applying the computed elastic transformation to the concatenated tracklets.

Contrary to gap-closing GDM approaches [33], EMC2 contains only two free parameters: the maximal distance between forward-propagated ending-points and backward-propagated starting-points of short tracklets for concatenation, and, to avoid important error propagation, a maximum time-lag between ending- and starting-point candidates. We highlight that concatenated tracklets do not necessarily span the entire time lapse sequence: each track begins with the first detection of its first concatenated tracklet and ends with its last detection of the last concatenated tracklet. Moreover, EMC2 algorithm handles complex natural motions and deformations, contrary to GDM methods that only account for confined or linear motion. EMC2 is therefore more robust and versatile. For the sake of reproducibility and dissemination of our method, we implemented the EMC2 multi-step procedure in the bio-image analysis software suite Icy [20] (http://icy.bioimageanalysis.org/). Icy is an open-source platform that is particularly well-suited for multi-step analysis thanks to graphical programming (plugin protocols) where each step of the analysis can be implemented as a block with inputs and outputs that can be linked to the other blocks (Fig 2 and Materials and Methods). Our method builds on well-established Icy preprocessing functions for spot detection and tracking.

Fig 2

Implementation of the EMC2 algorithm in Icy platform.

Time-lapse sequence of fluorescent particles is the input to a multi-step, automatic protocol in Icy. A first series of blocks, highlighted in blue, detects the position of fluorescent neurons (spots) in each frame of the time-lapse sequence. Block 1 uses the wavelet transform of each image and statistical thresholding of wavelet coefficients to determine spots that are significantly brighter than background. To separate close spots that form clusters in the wavelet-based mask of the image, the thresholded sequence is convolved with a log-gaussian kernel to enhance single spots (block 2), and local maxima algorithm is applied (block 3). A second series of blocks, highlighted in red, computes single particle tracks from computed spot positions. First, the Bayesian tracking algorithm (eMHT) computes tracklets of detectable particles (block 4). Due to fluctuating particle detectability, many Bayesian tracklets are terminated prematurely and new tracklets are created when particles can be detected again. To close detection gaps in single particle tracks, block 5 applies the EMC2 algorithm. Final output of the Icy protocol is the collection of single particle tracks over the time-lapse sequence. Tracking protocol can be found here: http://icy.bioimageanalysis.org/protocol/detection-with-cluster-un-mixing-and-tracking-of-neurons-with-emc2/ and is also directly accessible through the search bar of the Icy software (see step-by-step Supplementary Icy tutorial).

4 Validation of EMC2

a Manual tracking of calcium dataset in behaving animals

To validate the capabilities of EMC2, we first compared the results of our algorithm with manual tracking in calcium imaging sequences of neuron activity in behaving Hydra [5]. We used the first 250 frames (25 seconds at 10 Hz) of a time-lapse sequence previously acquired in a genetically-engineered animal [5] and automatically detected the active neurons (bright spots) using the multi-step detection process described in the Materials and Methods. We tracked the detected particles with the eMHT algorithm ([30], implemented in Spot tracking plugin in Icy) and obtained short Bayesian tracklets (n = 784 tracklets) for the detected neurons (step 4 of the Icy protocol in Fig 2). We then manually concatenated all the corresponding tracklets, i.e. we closed gaps, and obtained complete neuron tracks (n = 444 tracks). We observed that, before gap closing, tracklets were significantly shorter than concatenated tracks, meaning that many tracklets are indeed terminated prematurely by the undetectability of silent (non-firing) neurons. We measured the accuracy of EMC2 by comparing the computed tracks with those obtained after manual association of tracklets, which we took as an approximate ground truth (see Materials and Methods and S2 Fig). We also measured how tracks obtained with TrackMate [33] implemented in Fiji [34] matched this manual ground truth. TrackMate is a GDM method, based on optimal linear assignment between closest detections. To handle the gaps in tracking when particles are undetectable, TrackMate again uses a GDM algorithm to compute the optimal linear assignment between ending- and starting-points of previously computed tracks. As a result, TrackMate uses the same type of gap closing algorithm as EMC2 but without correcting for potential elastic deformation of the field-of-view. Finally, to evaluate how well the algorithm performs to generate the initial set of short tracklets, i.e. to compare the probabilistic eMHT used in EMC2 with the GDM algorithm used in TrackMate, we also measured the accuracy of EMC2, but without elastic motion correction before gap closing. Compared algorithms are summarized in Table 2. First, we found that EMC2 (n = 453 tracks, with 410 (90.5%) matched tracks) outperformed TrackMate (n = 474 tracks, with 259 (54.6%) matched tracks) and EMC2 without elastic correction (n = 514 tracks, with 279 (54.3%) matched tracks). The similar capabilities of TrackMate and EMC2 without elastic motion correction indicate that Bayesian eMHT and the GDM tracking method perform similarly for local association of detectable spots, but fail at closing longer tracking gaps in deformable media. This highlights the importance of elastic motion correction before the optimal concatenation of short tracks.

Table 2

Tested tracking algorithms in manual validation.

Name	Local association of detected particles	Elastic Motion Correction?	Gap closing	% match
Manual	Bayesian (eMHT)	No	Manual	100% (ground truth)
EMC 2	Bayesian (eMHT)	Yes	GDM	90.5%
TrackMate	Global distance minimization (GDM)	No	GDM	54.6%
EMC2 without Correction	Bayesian (eMHT)	No	GDM	54.3%

Using the same methodology, we compared manual and automatic (EMC2) tracking of single neurons in the less challenging case of two-photon imaging of the visual cortex (layer 2/3) of mice (we used the first day calcium recording from the first animal [38]). Here, the deformation of the field-of-view is much more limited and the motion of embedded neurons can be assimilated to confined diffusion. As expected, we obtained an EMC2 accuracy that was close to 100% (64 correct tracks over 65, i.e. 98.5% accuracy) for a time lapse sequence of 3,700 frames (5 minutes).

b Synthetic time-lapse sequences

Manual gap closing in time-lapse sequences is tedious and prone to operator bias. Moreover, the ground truth, i.e. the identity of each individual neuron along the whole time-lapse sequence, is unknown. Therefore, we designed a reproducible, synthetic approach where we simulated individual neurons’ activity and animal motion with different sets of parameters.

We modeled three different types of motion and/or deformation (Materials & Methods): Confined diffusion, where blinking neurons diffuse within a confined area (as in two-photon imaging of targeted brain areas), Linear motion, where neurons all move together in the same direction at constant velocity, and finally, using deformation fields measured in Hydra experimental data, we simulated naturalistic Hydra deformations. We further modeled the intermittent activity of neuronal ensembles with a probabilistic Poisson model. We also modeled the fluorescence dynamics of individual spikes using a parametric curve that we fitted to experimental data. Finally, using neuron positions, firing activity and fluorescence dynamics, we generated synthetic time-lapse sequences using a mixed Poisson-Gaussian noise model ([30] and Materials & Methods).

For confined diffusion (150 simulated tracks, n = 10 simulations), both EMC2 and TrackMate gave excellent results, with track matches of 93.5% ± 1.6% (144.4 ± 1.8 correct tracks over 154.6 ± 0.9 reconstructed tracks) for EMC2 and 92.9% ± 0.8% for TrackMate (140.3 ± 1.0 correct tracks over 151.0 ± 0.3 reconstructed tracks) (Fig 3A). The good performance of TrackMate was expected as this algorithm was initially designed to track confined endocytic spots at the cell membrane [33]. In addition to confined motion, TrackMate can also model linear motion of particles when computing the optimal gap closing between short tracks. Therefore, in linear motion simulations (337.6 ± 1.0 simulated tracks, n = 10 simulations), we used TrackMate with linear motion correction instead of standard confined motion correction. However, even with linear correction, the performance of Trackmate (76.3% ± 0.6% (262.6 ± 1.4 correct tracks over 358.2.0 ± 1.8 reconstructed tracks)) was significantly worse than EMC2’s performance (97.7 ± 0.5% (326.1 ± 1.6 correct tracks over 337.3 ± 1.1 reconstructed tracks)). This difference is due to the different estimation methods that are used in the two tracking algorithms to estimate the direction of tracks: in TrackMate, the estimation of track directions is local, based on the last detection within each short track, whereas the estimation of track direction in EMC2 uses global information provided by neighbouring short tracks and is therefore more robust. Finally, in the third case, we used the deformation field that we estimated over 250 frames within a time-lapse experimental sequence of Hydra (500 simulated tracks, n = 10 simulations, see section 2.a and Materials & Methods). We found that Trackmate had similar performances with (matching score 69.4 ± 1.1%, or 414.8 ± 3.7 correct tracks over 598.2 ± 4.8 reconstructed tracks) or without (72.4 ± 0.9%, or 427.2 ± 3.0 correct tracks over 590.7 ± 3.5 reconstructed tracks) linear motion correction, and that both were outperformed by EMC2 (98.6 ± 0.3%, or 493.4 ± 1.6 correct tracks over 500.4 ± 0.2 reconstructed tracks).

Fig 3

Testing EMC2 robustness with synthetic simulations.

For each simulated type of motion (confined diffusion (a), linear motion (b) and “Hydra-like” elastic deformation (c)), we simulated the stochastic firing of neuronal ensembles and corresponding fluorescence dynamics in synthetic time-lapse sequences (see Materials and Methods for details). We then compared the performances of EMC2 (blue) with TrackMate (no motion correction (red) or linear motion correction (magenta)). P-values are obtained with the Wilcoxon rank sum test over n = 10 simulations in each case. (d-e) Using “Hydra-like” synthetic deformation, we measured the accuracy of EMC2 for increasing proportion of stable (i.e. non-blinking) particles (neuronal cells) and increasing number of simulated particles. (f) After having estimated the deformation-field in three different animals (animal 1 (black), 4 (blue) and 6 (green)), we measured the accuracy of EMC2 for simulated sequences with increasing length (25, 50, 100 and 240 seconds. Imaging and simulations were performed at 10 Hz). For comparison purposes, the performance of TrackMate algorithm for 25 seconds (animal 1), extracted from (c), is shown.

We also measured the robustness of EMC2 to parameter change in synthetic motion simulations. In particular, we measured the performance of the algorithm for an increased percentage of stable cells (Fig 3D) (see Material and Methods), an increased number of neurons (Fig 3E) and an increased length of simulated sequences (Fig 3F and Table 3). First, we found that even without stable cells (α = 0), the accuracy of EMC2 remained high (77.7% ± 3.6% (398.8 ± 18.4 correct tracks over 513.3 ± 1.7 reconstructed tracks), and rapidly increased to 91.8% ± 1.2% for α = 5%, before reaching a plateau above 95% accuracy for α≥10%. Conversely, we found that EMC2 was very sensitive to the number of neurons, with poor performance for very few neurons (for simulations with only 10 neurons, the accuracy was 11.6% ± 0.6%, or 5.6 ± 0.3 correct tracks over 49.8 ± 0.8 reconstructed tracks). The accuracy rapidly increased to > 90% when more than 100 neurons were simulated. Tracking errors in simulations with few neurons are due to the inaccurate estimation of the deformation field, and iterative error propagation, when few fiducial source and target points are used. Finally, we measured how the performance of the EMC2 evolved with the length of the simulated sequence, for three different animals (animals 1, 4 and 6, see Table 4). We chose animals with a high and homogeneous density of neurons so that we could accurately estimate, and therefore simulate, their body deformation. We found that the EMC2 accuracy decreased with the sequence length, but remained high (i.e. >80%) for up to 2 to 4 minutes of simulation at 10 Hz, depending on the animal (animal 1 ~ 1200 frames = 2 min., animal 4 ~ 1500 frames = 2 min. 30s and animal 6 ~ 2400 frames = 4 min.). For longer simulations, the accuracy dropped below 80% and reached a mean of 72% at the end of the simulation (4 minutes). The number of reconstructed tracks remained close to 500 even for longer sequences (Fig 3F and Table 3). The decreased performance of EMC2 for longer time-lapse sequences is an expected drawback of SPT-based tracking algorithms. Indeed, the probability of false associations between newly detected particles and existing tracks increases with time. Yet, the EMC2 tracking algorithm, in contrast with state-of-the-art SPT algorithms, allows the robust, automatic tracking of individual neurons over few (~2–3) minutes in the highly deforming Hydra model. This sustained performance is particularly desirable as it allows the robust analysis of single neuron activity and functional coupling during different animal behaviors (section 4).

Table 3

Results of synthetic simulations (Hydra-like deformation, increasing length).

	250 frames (25 s.)		500 frames (50 s.)		1000 frames (100 s.)		2400 frames (240 s.)
	Accuracy	Nb. Tracks	Accuracy	Nb. Tracks	Accuracy	Nb. Tracks	Accuracy	Nb. Tracks
Animal 1	98.6±0.3%	500.4±0.2	87.5±0.8%	511.7±1.0	83.8±1.8%	512.0±3.1	67.4±4.1%	538.8±8.0
Animal 4	98.0±0.2%	503.5±0.4	95.9±0.4%	507.9±1.1	88.7±1.2%	521.5±2.6	66.7±2.5%	571.9±8.5
Animal 6	97.8±0.2%	502.9±0.7	96.9±0.3%	505.6±0.7	92.0±1.3%	515.7±2.3	81.7±1.9%	540.7±5.4

Table 4

Results of statistical extraction of neuronal ensembles in Hydra (n = 8 animals).

Animal	Movie length (frames)	Tracks (>150 frames)	Activity peaks	Ensembles	CB neurons	RP1 neurons	RP2 neurons
1	900	411	11	2		139 (33.8%)	37 (9.0%)
1	1636	723	37	3	202(27.9%)	121(16.7%)	49(6.8%)
1	1371	607	27	3	222(36.6%)	78(12.9%)	41(6.8%)
1	988	516	23	3	220(42.6%)	49(9.5%)	55(10.7%)
2	1000	531	17	2	34(6.4%)	29(5.5%)
2	1000	518	17	2	59(11.4%)	29(5.6%)
3	1000	191	20	2	16(8.4%)	11(5.8%)
4	1000	243	13	2	16(6.6%)	17(7.0%)
4	1000	173	16	2		24(13.9%)	8(4.6%)
5	1000	147	15	2	14(9.5%)	7(4.8%)
6	1000	379	11	2	14(3.7%)	10(2.6%)
7	1000	173	9	2	59(34.1%)	25(14.5%)
8	1000	244	10	2	44(18.0%)	18(7.4%)
Mean	1067	374	17.3		82(18.7%)	43(10.8%)	38(7.6%)
Standard error	56	53	2.2		26(4.4%)	12(2.3%)	8(1.8%)

Altogether these simulations show that EMC2 is a robust tracking algorithm for different types of particle motions, and is therefore a versatile method for single particle tracking.

5 Tracking and analyzing single neuron activity in behaving animals

a. Monitoring neuron activity with two-photon microscopy in mouse cortex

Two-photon calcium imaging is widely used to monitor single neuron activity in targeted brain regions of awake and behaving animals [3]. The head of the animal is usually fixed under the objective of the microscope, which limits the motion of neurons and facilitates their individual tracking. However, residual neuron motion due to animal movements, breathing or heartbeats requires computational post-processing of acquired time-lapse sequences to robustly monitor calcium activity of single neurons. The most popular technique for (slight) motion correction in time-lapse sequences is the elastic registration of fluorescence images with respect to one (or multiple) reference frame(s) using image intensity [13]. The computational load of image registration is important, especially for long sequences with large images, and algorithms have been developed to speed up the registration process and decrease the computational time to a few minutes for ~2,000 frame time-lapse sequences (with ~256x256 images) [14]. After image registration, the segmentation of neuronal masks for calcium fluorescence monitoring is then performed with standard intensity thresholding [38] or more elaborate techniques when neuronal masks are overlapping, such as non-negative matrix factorization [17].

To simultaneously localize and correct for neuron motion, we applied the EMC2 Icy tracking protocol (Fig 2) to two-photon time-lapse calcium images of mouse visual cortex (Fig 4 and Material and Methods). Imaging was performed for 5 minutes at 12.3 Hz (~3,700 frame sequences) at days 1, 2 and 46. The entire EMC2 protocol with neuron spot detection and tracking for each ~3,700 frame sequence ran in ~3 minutes with a 2.7 GHz Intel Core i7 processor. Consistent with [38], we observed a significant turn-over of active neurons (Fig 4A) with few neurons (median = 22% (15/68 neurons), n = 4 animals) that remained active across all days. We observed a similar number of active neurons on day 1 and 2, but a decreased number at day 46 which is probably due to several factors, such as decreased transgene expression or repeated experimental procedures [38]. We then analyzed single neuron trajectories obtained with EMC2 (Fig 4B). Even if the animal’s head was fixed under the microscope, residual motion of the field-of-view led to confined stochastic trajectories for single neurons. The positions of single neurons at each time were either computed with the intensity center of detected spots when neurons were detectable, or estimated using the computed elastic deformation of the field of view when neurons were silent and undetectable. We measured a median neuron displacement between frames of ~ 0.25 pixels, and a median maximum distance of excursion (relative to the center point of the trajectory) of ~ 4 pixels (Fig 4B).

Fig 4

Monitoring the activity of individual neurons in two-photon calcium imaging of mouse visual cortex with EMC2.

a- Two-photon calcium imaging of single neuron activity in visual cortex of awake mice is performed at Day 1, Day 2 and Day 46 during 5 minutes at 12.3 Hz. Tracking of neuron positions is performed with EMC2 and reveals an important turn-over of active neurons across days. Examples of neurons that are active at Day 1, Day 2 or Day 46 are respectively highlighted with red, green or blue arrows. Neurons active at Day 1&2 are highlighted with yellow arrows, at Day 2&46 with cyan arrows, and at Day 1,2&46 with white arrows. The median number of active neurons each day is also plotted (n = 4 animals). The number of neurons that are active from Day 1, 2 or 3 are respectively represented in red, green or blue. b- Single neuron trajectories can be modeled with confined stochastic motion. Two example trajectories are shown (green & blue trajectories) with a maximum excursion distance of ~ 4 pixels. Boxplots of single neuron displacement between two consecutive frames, and maximum excursion distance (in pixels) are plotted (n = 590 trajectories).

Altogether, these results show that EMC2 tracking protocol is a robust and fast method to post-process two-photon calcium imaging from awake mice. Trajectory analysis revealed the stochastic confined motion of single neuron positions, even in head-fixed animals. This residual motion is partly due to the animal’s movements, but also to the uncertainty of sub-pixel localization of neuron spots at each time frame [39]. Moreover, the analysis of single neuron activity across days in the same cortical region showed a significant turn-over within the pool of active neurons each day, with few neurons remaining active over all days. Statistical analysis recently showed that this latter pool of active neurons could be a stable neuronal ensemble [38].

b. Characterizing neuronal ensembles in behaving Hydra

There is increasing experimental evidence that neurons are organized into neuronal ensembles composed of a few tens of highly coupled neurons, and that these co-active ensembles are the fundamental computational units of the brain rather than single neurons themselves [21, 40]. Using manual annotation, it has been shown that Hydra’s nervous system, one of the simplest of the animal kingdom, may be dominated by three main functional networks that are distributed through the entire animal [5, 41]. To confirm (or refute) these observations, we used EMC2 and automatically tracked single neurons in n = 13 time-lapse sequences (length 1067±56 frames at 10 Hz (Materials and Methods) from 8 different animals (Table 4 and Fig 5). Movies used in this study (avi and tiff files) can be downloaded from the BioStudies platform (https://www.ebi.ac.uk/biostudies/studies/S-BSST428). Analyzed movies were significantly longer than the manually annotated one (length = 200 frames at 10 Hz [5]). Ensemble activity, corresponding to the co-firing of neurons, can be detected as significant peaks within the raster plot of single neuron activity (Fig 5 and Materials and Methods). We detected a mean number of 1.7±2.2 peaks per movie, corresponding to a mean rate of 1 activity peak every 63 frames (6.3 s) which corresponds well with the 4 peaks observed previously in the 200 frame movie [5]. To associate each peak with a putative neuronal ensemble, we adopted a similar approach as in [42] and measured the similarity between these events in terms of the identities of the participating neurons (Fig 5 and Materials and Methods). We then performed k-means clustering of peak similarity, and used the Silhouette criterion [43] to determine the most likely number of neuronal ensembles causing the detected peaks of activity. We found 2 or 3 neuronal ensembles in each movie (3 neuronal ensembles were detected in 3 out of the 13 movies, or 23%). We then categorized each detected ensemble into one of the previously defined ensembles [5, 41]: Contraction Burst (CB) neurons that fire during longitudinal contraction of the animal, Rhythmic Potential 1 (RP1) that fire during the longitudinal elongation of the animal, and Rhythmic Potential 2 (RP2) neurons that fire independently of RP1 and CB activity. We found CB neuronal ensembles in almost all movies (11/13, or 85% of movies) and all animals (8/8, 100%), RP1 ensembles in all movies and animals and RP2 ensembles in fewer movies (5/13, or 38%) and only 2/8 (25%) animals. The absence of detected CB ensembles in 2 movies corresponds to the observed absence of contraction cycles within these movies. On the other hand, we hypothesize that the absence of RP2 ensembles in 6/8 (75%) of the animals is due to the limited depth of the field-of-view in confocal microscopy (see Materials and Methods). Indeed, RP2 neurons lie in the thin ectoderm of the animal [5] that may not have been imaged in some animals. Finally, we classified and mapped each individual neuron in the detected ensembles (see Materials and Methods). As in [5], we found that the CB ensemble was the most important group of neurons with a mean number of 82±26 neurons representing 18.7%±4.4% of the total number of neurons (a mean of 374±53 neurons were tracked over >150 frames in the different movies). The RP1 ensemble, with 43±12 neurons (10.8%±2.3%), was the next largest ensemble and RP2, with 38±8 (7.6±1.8%), was the third. The relative number of neurons in the different ensembles is in agreement with previous observations [5]. However, the overall size of each ensemble is smaller than the size reported previously. This is due to the fact that automatic classification of each individual neuron in an ensemble is more stringent than the manual classification that had been previously performed. Finally, we measured the overlap between ensembles, i.e. the proportion of single neurons that belonged to more than one ensemble. For each pair of ensembles, we computed the ratio between the number of shared neurons and the total number of neurons in both ensembles. We obtained ratios of 2.1%±0.7%, 2.4%±0.4% and 1.6%±0.6% respectively for CB-RP1, CB-RP2 and RP1-RP2 ensemble overlap. These very low values confirm the near non-overlap of the main neuronal ensembles in Hydra [5]. By coupling the automatic tracking of individual neurons in multiple Hydra with robust clustering analysis of neuronal activity, we were therefore able to confirm the previous observations that the Hydra nervous system is dominated by three main non-overlapping ensembles that are involved in different animal behaviors.

Fig 5

Neuron tracking and mapping of neuronal ensembles in behaving Hydra.

a- Calcium imaging of single neuron activity in behaving Hydra. The images and analysis of the 3rd movie of animal 1 are given as representative examples. b- Single neuron tracks and fluorescence intensity are obtained with EMC2 algorithm. c- For each neuron, spikes are extracted from fluorescence traces. Peaks of activity (highlighted with red stars) correspond to significant co-activity of individual neurons (sum of individual activities (solid red line) > statistical threshold (dashed red line), p = 0.001 see Materials and Methods). Each peak putatively corresponds to the activation of one neuronal ensemble. d- Similarity between activity peaks is computed using the identities of individual neurons that are firing at each peak (see Materials and Methods). e- The optimal number of peak classes (that putatively corresponds to the number of neuronal ensembles) is computed using the Silhouette index on k-means clustering of the similarity matrix (see Material and Methods). Median fluorescence trace of each neuronal ensemble and corresponding activity peaks are shown. The classification of individual neurons in each ensemble is determined based on their firing at ensemble peaks (see Materials and Methods). f- Individual neurons of each ensemble can be dynamically mapped in the original time-lapse sequence.

Altogether, these results show that EMC2 is sufficiently robust to monitor single neuron activity in behaving animals. This constitutes a fundamental prerequisite for the analysis of neurons’ functional organization and, ultimately, for our understanding of their emergent computational properties.

Discussion

When using calcium imaging in living animals, an important challenge is the sustained tracking of neuronal positions over extended times in a moving, and potentially also deformable, environment. To tackle this issue, we have developed an algorithm, EMC2, that tracks detectable particles with a state-of-the-art probabilistic tracking algorithm and uses the information contained in reconstructed tracks about the local deformation of the field-of-view (e.g. animal) to estimate the position of undetectable particles and potentially close tracking gaps. We validated the performance and versatility of EMC2 by comparing its performance with a state-of-the-art tracking algorithm on manually tracked neurons in time-lapse calcium-imaging of behaving animals, including imaging over days of the same neurons from mouse visual cortex, the challenging deformable Hydra, and also on synthetic time-lapse sequences that modeled different types of motion/deformation of the field-of-view (confined diffusion, linear motion and Hydra-like elastic deformation). In all cases, EMC2 showed high accuracy and outperformed state-of-the-art tracking methods.

Compared to traditional tracking approaches, composed only of particle detection and linking (see Table 1), our hybrid algorithm is better equipped to handle tracking gaps and general particle motion. Recently, tracking methods based on artificial neural networks [44, 45] have been introduced to handle more general particle motion, not only diffusion and/or linear motion. However, these methods are designed to link two sets of particle detections in consecutive time frames and, even if they can handle some missing or false detections, they do not appear well suited for tracking gaps such as those encountered in calcium imaging of neuronal activity.

On-going development of high-speed three-dimensional microscopes enables the imaging of neuronal activity with high temporal and spatial resolution in an increasing number of animal models [6, 9]. In the near future, the extension of EMC2 tracking method to three-dimensional movies should not present important technical issues. Indeed, the thin-plate-spline transform used here to estimate the deformation of the field-of-view in two dimensions is actually a special case of polyharmonic splines that have been specifically designed for robust interpolation between data points in any dimensional space[46].

After having implemented EMC2 in the open-source bio-imaging platform Icy [20], we tracked single neuron activity in behaving mouse visual cortex and also Hydra. In mouse cortex, our algorithm performed essentially flawlessly, and reveals the presence of neurons that continue to respond to the same stimulus over several days. In addition, using statistical clustering, we confirmed the previous observations that neural activity of Hydra is dominated by three major, non-overlapping neuronal ensembles that are involved in the animal’s repetitive contraction and elongation. At the same time, automatic clustering of neuronal activity was not able to extract other smaller ensembles of the animal’s nervous system, such as the tentacle and sub-tentacle ensembles [5] that we could observe by eye. Indeed, the coordinated activity of small ensembles is difficult to detect from individual neuronal spiking because of sparse activity and noise. Moreover, the activation of these smaller ensembles is less frequent than the activation of the three major ensembles and happened in only few movies. Therefore, the complete mapping of neuronal ensembles of Hydra, with the characterization of even the smallest neuronal ensembles with less frequent activity, will require in the near future further development of better imaging and tracking of single neuron activity over longer periods of time. For tracking neurons over longer time (several minutes or even hours), hardware implementations such as dual-color imaging with calcium-insensitive dye or partial immobilization of the animal [47] will be required. Another possible software strategy for increasing the robustness of EMC2 tracking over long times would be the characterization and use of stereotypical poses of the animal [8] as reference frames for mapping the (almost) stable positions of neuronal sub-sets.

To conclude, our results show that EMC2 is a robust and versatile tracking algorithm that allows monitoring and quantification of single neuron activity in behaving animals. Robust tracking of neural activity is a first step towards a better understanding of the neural code, i.e. how connected neuronal ensembles integrate information, underlie adaptive behavior and, more generally, compute the animal’s behavioral or internal states.

Materials and methods

1 Elastic motion correction and concatenation (EMC2) of short tracks

The eMHT algorithm returns a set of N short tracks (tracklets) τ, for 1≤i≤N, of detectable particles; the ith tracklet τ starts at position at time and ends at position at time . To concatenate tracklets, i.e. to link the tracklets that putatively correspond to the same neuron, we then estimate the backward position of the particle (neuron) corresponding to the i tracklet before the starting-point () of the tracklet. To estimate for , we iteratively apply backward thin-plate-spline transformation , with the initial condition . Similarly, to estimate the forward position after the ending-point of the tracklet, i.e. for , we iteratively apply the forward transformation with the initial condition From the estimated positions of particles corresponding to each tracklet, we then compute the distance d between tracklets τ≠τ with where gap is a user-defined maximum time gap that EMC2 is allowed to close (typically a few hundreds of frames). We need to apply a maximum time gap for time-lapse sequences with particles that remain undetectable over long periods of time due to the error growth during backward/forward estimation of the putative position of undetectable particles.

Using the computed distances between tracklets, we then define an association cost matrix, with ϕ = d if d

This minimization problem, known as the assignment problem, is similar to the problem solved in GDM methods of tracking, where algorithms determine the optimal association between particle detections by minimizing the global distance between detections in consecutive time frames of the sequence. One of the first and most popular algorithms to solve assignment problems is the Hungarian algorithm [48]. However, due to its computational load, faster algorithms have been proposed over the years. We used here the Jonker-Volgenant algorithm [49] implemented in the TrackMate plugin in ImageJ (see Table 1).

Finally each particle (neuron) track T for 1≤i≤N′, with N′≤N the number of “long” tracks, results from the concatenation of n tracklets: with . We highlight that, by construction, each long track T does not necessarily span the whole time-lapse imaging sequence, but begins at the starting-time of the first tracklet τ and ends at the ending-time of the last tracklet .

2 Icy protocol

a Detection of spots (e.g. neurons) in time-lapse sequences

To detect automatically the positions of fluorescent spots, corresponding to detectable particles, in each frame of the time lapse sequence, we designed a multi-step algorithm (see Fig 2) where we first detected fluorescent spots that are significantly brighter than background with a fast and robust algorithm based on a wavelet transformation of the image and statistical thresholding of the wavelet coefficients (block number 1 in Icy protocol (Fig 2)) [50]. These spots correspond to individual particles or clusters of particles (e.g. neurons). To separate individual particles in the detected clusters, we then multiplied the original sequence with the binary mask obtained with wavelet thresholding and convolved the result of the multiplication with a log-Gaussian transformation (block number 2). The log-Gaussian convolution is similar to the point-spread function of microscopes and thus enhances individual particles [51]. Finally, we extracted the positions of single particles by applying a local-maxima algorithm (block number 3) to the convolved sequence.

b Tracking (EMC2)

After having detected the positions of fluorescent spots in each time frame, a second series of blocks computed the tracks of each single particle. First, block number 4 used the positions of spots and a robust Bayesian algorithm (eMHT [30]) to compute single tracks of detectable particles. Due to fluctuating detectability, many computed tracks are terminated prematurely and new tracks are created when particles are detectable again. We thus applied the EMC2 algorithm (block number 5) to close gaps and reconstruct single-particle tracks over the entire time lapse sequence.

The tracking protocol can be found here: http://icy.bioimageanalysis.org/protocol/detection-with-cluster-un-mixing-and-tracking-of-neurons-with-emc2/ and is also directly accessible through the search bar of the Icy software. A step-by-step tutorial for tracking neurons and exporting track intensity with Icy is provided as Supplementary material.

3 Validation of EMC2

a Comparison metric

To compare the tracks obtained with EMC2 and other automatic tracking algorithms with ground truth tracks, we first considered the whole set of detections x(t), 1≤i≤N(t), with N(t) the number of detections at time 1≤t≤T (T being the length of the time sequence) and assigned each detection to the closest active track at time t. Therefore, for each reference track , 1≤j≤|Θ|, with |Θ| the total number of reference tracks, and for each test track , 1≤k≤|Θ|, with |Θ| the total number of test tracks, we obtained a set of associated detections. We then considered that a test track matched a reference track if it shared at least 80% of common detections. Finally, for each reference track, we either obtained no test track that matched, exactly one test track that matched or more than one.

b Synthetic motions

To validate the robustness and accuracy of EMC2 in different scenarios, we simulated three classes of motions: confined diffusion, linear motion and elastic deformation. For confined diffusion, each simulated spot (e.g. fluorescent neuron) can diffuse with coefficient D = 1 pixel per frame and is confined to a 10 pixel disk area. For linear motion, each simulated spot moves linearly at speed v = 1 pixel per frame. When a track reaches the boundary of the field-of-view (a rectangle of 200x200 pixels), it is terminated and another track is initiated at the other side of the FOV. Finally, for elastic deformation, we used the experimental tracks in Hydra to estimate iteratively (i.e. from one frame to the following one) the local deformation for each synthetic track position.

c Firing rates of individual neurons

To model the stochastic firing rates of individual neurons (total number of neurons n), we first determined a proportion (α) of stable spots, i.e. non-blinking cells, with constant fluorescent intensity. In Hydra, stable cells typically correspond to nematocytes or other cell types that also express fluorescent proteins after the genetic editing of the animal, but that don’t fire as neurons do [5]. To simulate the correlated activity patterns observed in Hydra, we then divided the (1−α)n firing neurons, with intermittent activity and detectability, into n ensembles. All neurons in each ensemble fire simultaneously with Poisson rate λ = size λ, with λ the firing rate of individual neurons and size = (1−α)n/n the number of neurons in each group. Parameters for each simulation used for the validation of EMC2 algorithm are summarized in Table 5.

Table 5

Parameters used for synthetic simulations.

Parameters	n neuron	α stable	n group	λ individual	A	τ decay	β	μ	τ rise	σ PSF	σ n	B	SNR=AB+σn2
Name	Total number of neurons	% of non-firing spots	Number of neuron groups	Individual firing rate	Amplitude	Fluorescence decay time constant	Decay power	Median rising time	Rising time constant	St. Dev. of the PSF	St. Dev. of the Gaussian noise	Poisson background	Signal-to-noise ratio
Confined diffusion & Linear motion	150	20%	10	0.01 frame-1	100	3 frames	1	1 frame	0.5 frames	1 pixel	5	10	≈3
Elastic deformation (Hydra)	500	20%	10	0.0002 frame-1	100	15 frames	2	2 frames	0.5 frames	1 pixel	5	10	≈3

d Generation of synthetic images

To generate synthetic fluorescence time-lapse sequences, we used a mixed Poisson-Gaussian model [31]. In this model, the intensity I[x,y] at pixel location [x,y] is equal to I[x,y] = U[x,y]+N(0,σ) where U is a random Poisson variable and N(0,σ) is additive white Gaussian noise with standard deviation σ. The intensity λ[x,y] of the Poisson variable varies spatially because it depends on the presence or not of particle spots (neurons). Therefore, λ[x,y] = P[x,y]+B with B a constant background value and P[x,y] the spots’ intensity at position [x,y]. Assuming an additive model for the intensity of the spots, , where P[x,y] is the signal originating from the i spot. We approximated the point-spread-function (PSF) of the microscope with a Gaussian profile. Thus, for a particle located at position , its intensity at position [x,y] is given by , with A the amplitude of the i particle and σ the standard deviation of the 2D Gaussian profile of the PSF. We chose a constant amplitude for each particle A = A, for all 1≤i≤n. Parameters used in simulations are summarized in Table 5.

e Fluorescence kinetics

When a neuron fires at time t0, we model its fluorescence time course with the general kinetics equation where the numerator models a power-law exponential decay of the fluorescence (β = 1 models a standard single exponential decay), with a decay time constant τ, and the denominator models a sigmoidal increase of fluorescence with median τ and time constant μ. Kinetic parameters for each synthetic simulation are summarized in Table 5. These parameters were obtained by fitting n = 3075 individual spikes from 444 individual neuron tracks in an experimental time-lapse sequence (250 frames at 10 Hz) of GCAMP-labeled Hydra [5] (S1 Fig). We highlight that, for Hydra elastic simulations, we used a long decay time constant and a power index β = 2 instead of 1 for standard confined diffusion and linear motion simulations.

4 Tracking and analyzing single neuron activity in living animals

a Two-photon calcium imaging of mouse visual cortex

Movies used in this study are issued from [38], and experimental protocol for two-photon volumetric imaging of targeted brain regions in mouse visual cortex can be found in the Methods section of this manuscript.

b Hydra maintenance

were cultured using standard methods [52] in Hydra medium at 18°C in the dark. They were fed freshly hatched twice per week.

c Hydra Imaging

Transgenic Hydra expressing GCaMP6s in the interstitial cell lineage were used and prepared for imaging studies as previously described [5]. Calcium imaging was performed using a custom spinning disc confocal microscope (Solamere Yokogawa CSU-X1). Samples were illuminated with a 488 nm laser (Coherent OBIS) and emission light was detected with an ICCD camera (Stanford Photonics XR-MEGA10). Images were captured with a frame rate of 10 frames per second using either a 6X objective (Navitar HRPlanApo 6X/0.3) or a 10X objective (Olympus UMPlanFl 10x/0.30 W).

All movies used in this study can be downloaded from the BioStudies website https://www.ebi.ac.uk/biostudies/studies/S-BSST428.

d Extracting single neuron activity

We extracted the fluorescence trace of each individual neuron using the Track Processor Intensity profile within the TrackManager plugin in Icy [20] (a step-by-step tutorial for tracking neurons and exporting track intensity with Icy can be found on the EMC2 plugin web documentation: http://icy.bioimageanalysis.org/plugin/elastic-motion-correction-concatenation-emc2-of-tracks/). For each detection within the track, the extracted intensity corresponded to the mean intensity over a disk centered at the detection’s position, with a 2 pixels diameter. When detections are missing (tracking gap), the intensity was set to 0. Then, for each individual neuron, we denoised its non-zero fluorescence trace using wavelet denoising (wdenoise) in Matlab. We then automatically extracted spikes with a custom procedure where we first computed the discrete derivative of the smoothed fluorescence signal, then set to 0 all negative variations and finally, we detected significant positive variations of the signal (discrete positive derivative > quantile at 98% of all empirical positive variations) that putatively corresponded to spikes.

e Statistical characterization of neuronal ensembles in Hydra

Neuronal ensembles are groups of neurons that repeatedly fire together. Therefore, the activity of neuronal ensembles can be detected as significant co-activity peaks in the raster plot of single neuron firing. To detect significant peaks of activity, we applied the procedure described in [42], and identified as peaks times at which the sum of single neuron activity within a time step of 100 ms fell within the quantile at 0.999% obtained empirically by circularly shuffling the individual spikes in the activity raster plot.

Then, to relate detected peaks of activity to putative neuronal ensembles, we constructed a vector describing the activity of each individual neuron at the detected peaks with entries 1 if the neuron fires at that peak, and 0 otherwise. We then computed the similarity between these vectors for each of the activity peaks, using the Jaccard index:

Then, to estimate the number of neuronal ensembles underlying the detected peaks of activity, we clustered the peaks with a k-means algorithm based on their similarity for different numbers of classes (from 1 to 5 classes typically). K-means clustering was performed using the cosine distance. The optimal number of classes in the k-means clustering algorithms, and therefore the putative number of neuronal ensembles, was computed using the Silhouette index [43]. For a clustering of N peaks into k classes, the Silhouette index is given by where a is the average distance from the i peak to the other peaks in the same cluster as i, and b is the minimum average distance from the i peak to peaks in a different cluster, minimized over clusters. An advantage of the Silhouette evaluation criterion over other clustering criteria is its versatility, as it can be used with any distance (cosine distance was used here for the k-means clustering). Finally, we assigned each individual neuron to an ensemble if that neuron fired in more than 50% of the activity peaks of the identified ensemble.

Fitting the fluorescence kinetics of firing neurons in calcium imaging of Hydra.

Fluorescence time course of a firing neuron (at time t) is modeled with the kinetics equation where the numerator models a power-law exponential decay of the fluorescence with a decay time constant τ, and the denominator models a sigmoidal increase of fluorescence with median τ and time constant μ. Boxplot for fitted parameters (n = 3075 individual spikes from 444 individual neuron tracks) are shown.

(EPS)

Click here for additional data file.

Comparing tracking methods with manual ground truth on calcium imaging data in behaving Hydra.

Tracks obtained with TrackMate (red), EMC2 without motion correction (purple) and EMC2 are compared to manual ground truth obtained by manual concatenation of tracklets in calcium imaging of behaving Hydra (see Material and Methods).

(EPS)

Click here for additional data file.

21 Oct 2020

Dear Dr. Lagache,

Thank you very much for submitting your manuscript "Robust tracking of single neuron calcium dynamics in behaving Hydra" for consideration at PLOS Computational Biology.

As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Boris S. Gutkin

Associate Editor

PLOS Computational Biology

Lyle Graham

Deputy Editor

PLOS Computational Biology

***********************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: The review is uploaded as a PDF attachment.

Reviewer #2: Lagache and colleagues present a complete user-friendly software package for tracking neurons during short calcium recordings in moving and deforming hydra. They demonstrate that their method can track neurons even when those neurons are inactive, without the need for an additional fluorescent marker. The development of effective algorithms that are integrated into user-friendly and open source software packages like Icy is very valuable for the community.

There are two substantive concerns.

1) It is unclear how the method performs on recordings longer than 25 seconds.

I could not find any description of the effective recording length for the validation performed against synthetic data with empirically derived noise, so I don't know if that was greater than 25s. But it is important to show that the method can work over longer recordings because calcium recordings are often many minutes long, and the sequential nature of this method would seem to be particularly susceptible to accumulating errors over time. For example, I might reasonably expect to accumulate errors in stitching tracks together such that after a few minutes a significant fraction of neurons may now be assigned different labels at the end of the recording compared to the beginning. One way the authors could address this concern is by showing how performance degrades over recording length when validated against realistic synthetic data.

2) Authors should clarify novelty and correct description of a previous method.

A key selling point of the method is that it uses deformation defined by tracked neurons to infer the location of missing neurons during gaps. The text highlights this method extensively, as if it were novel. But method [9], Ngueyn et al, explicitly performs the same step. Like this method, [9] uses thin plate splines of tracked neurons to infer the position of missing neurons, see main text accompanying Fig 5 in that reference [9].

The authors make a special point to explicitly compare to reference [9] in detail (for example on p. 7) and they also compare to many other methods in Table 1 and throughout the manuscript, all of which is commendable. However, they incorrectly claim that [9] is " designed to link two sets of particle detections in consecutive time frames." But this is incorrect. The method in [9] is time independent, in the sense that all the recorded volumes could be shuffled in time without changing tracking performance. Similarly, Table 1 incorrectly states that ref [9] is not robust to long gaps. On the contrary, because it is time independent, ref [9] is robust to exceptionally long gaps.

The authors should clarify novelty and correct descriptions of the previous method.

Additional technical concern:

There seems to be an inconsistency with Figure 4. It shows calcium activity lasting 1,400 s or 23 mins, whereas the main text and Table 3 suggests that recordings lasted less than 2 mins (1000 frames @ 10hz). Please resolve this inconsistency.

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biology data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: None

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, PLOS recommends that you deposit laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions, please see

Submitted filename: review_PCB_Lagache_2020.pdf

Click here for additional data file.

24 Feb 2021

Submitted filename: Response to reviewers.docx

Click here for additional data file.

30 Mar 2021

Dear Dr. Lagache,

Thank you very much for submitting your manuscript "EMC2: A versatile algorithm for robust tracking of calcium dynamics from individual neurons in behaving animals" for consideration at PLOS Computational Biology.

As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

In particular, it is necessary that you address the concerns of Reviewer #2 regarding the relevant time scale for the method, which is a critical factor for it's experimental application.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Boris Gutkin

Associate Editor

PLOS Computational Biology

Lyle Graham

Deputy Editor

PLOS Computational Biology

***********************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: The authors have more than adequately adressed all my relevant comments. The addition of a second test case and the exploration of the performances of EMC2 with respect to the fraction of stable cells, cell density and time of recording both help in demonstrating the versatility of the method. I believe the paper will acquire greater visibility and applicability as it now stands. I also appreciated the clarifications the authors made on the stitching procedure both in the Methods section and in modifying Fig1. Finally, quantification of performances of the different methods in terms of reconstructed tracks is very instructive and allows for a rational comparison. I strongly support publication of this paper.

Reviewer #2: In this revision Lagache and colleagues have made important revisions and clarifications. It remains unclear, however, whether the method can track neurons in hydra over a sufficiently long timescale to be useful.

The authors specifically motivate their work by the need to "achieve robust and automatic tracking of individual cells over long time-lapse sequences." Most calcium activity experiments are for at least a few minutes, and recent published work by this group recorded calcium activity in moving hydra for 2 hours (Yamamoto and Yuste 2020, see Fig 3). It is therefore imperative that the authors show what accuracy to expect on the 5 min timescale and beyond.

With synthetic data, the authors now show that their method's accuracy drops from 97% for 250 frames (presumably 25s) to an accuracy of ~80% for 1000 frames (presumably 100s). It seems possible that accuracy will continue to drop dramatically as the timing approaches 5 mins.

Related: synthetic data is reported in terms of frames, but the authors should instead report the equivalent time in seconds, based on the timescale of the motion in the synthetic data. The authors should make this easy to find in the text and figure or caption.

If the method is unable to achieve reasonable accuracy on the multiple minutes timescale, the authors should scale back their claims, and it may also reduce the significance of the work.

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biology data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: None

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

27 Apr 2021

Submitted filename: Cover_Letter_Revision_2.docx

Click here for additional data file.

13 May 2021

Dear Dr. Lagache,

Thank you very much for submitting your manuscript "EMC2: A versatile algorithm for robust tracking of calcium dynamics from individual neurons in behaving animals" for consideration at PLOS Computational Biology.

As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments. In particular the authors must address the issues brought up by reviewer 1 on the limitations of the method for extended recording times that had come up after the last revision. Indeed, in order to meet the editorial requirements of PLoS CB it is necessary to show that the method presents a significant advance on existing methods

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Boris S. Gutkin

Associate Editor

PLOS Computational Biology

Lyle Graham

Deputy Editor

PLOS Computational Biology

***********************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: I still support the publication of this paper.

Reviewer #2: With further simulation studies, the authors now show that their method's performance degrades as recording duration increases. For Hydra, the method is only suitable for very short recordings (<5 min timescale). This contradicts claims made in the abstract and throughout the text related to "long time lapse sequences" or "long-term" monitoring. Claims about closing potentially "long" gaps also lack evidence. In short there is a disconnect between many of the authors' claims and the evidence.

Evidence does not support the authors claim that that this method "greatly supersedes existing approaches." In contrast, other methods (NeRVE from Ngueyn et al. 2017; Conditional Random Fields from Chaudhary et al. elife 2021) are fundamentally time-independent and thus maintain the same performance regardless of recording duration. The authors need to revise their comparisons to other methods with a frank assessment of this limitation, which seems to arise from the sequential time-dependent nature of the algorithm. For example, it would be relevant to report in Table 1 whether each algorithm is time-dependent or not.

Crucially, the authors need to better justify why this method is meritorious despite being limited to only very short recordings. To be compelling the authors should provide specific use cases, and also include a discussion of why and how errors accumulate with time and what fundamentally about the algorithm imposes this limitation.

Minor, but important: it is imprecise to claim that accuracy "remained above ~70% for all animals" when Table 3 reports accuracy of 66.7% and 67.4% for two of the three animals.

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: None

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Olivier Cochet-Escartin

Reviewer #2: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

25 Jun 2021

Submitted filename: Cover_Letter_Revision_3.doc

Click here for additional data file.

7 Jul 2021

Dear Dr. Lagache,

Thank you very much for submitting your manuscript "EMC2: A versatile algorithm for robust tracking of calcium dynamics from individual neurons in behaving animals" for consideration at PLOS Computational Biology. As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. The reviewers appreciated the attention to an important topic. Based on the reviews, we are likely to accept this manuscript for publication, providing that you modify the manuscript according to the review recommendations. Please follow the suggestions laid out by the reviewer 2.

Please prepare and submit your revised manuscript within 30 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to all review comments, and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Thank you again for your submission to our journal. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Boris S. Gutkin

Associate Editor

PLOS Computational Biology

Lyle Graham

Deputy Editor

PLOS Computational Biology

***********************

A link appears below if there are any accompanying review attachments. If you believe any reviews to be missing, please contact ploscompbiol@plos.org immediately:

[LINK]

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: I thank the authors for the clarifications of time-dependent and time-independent methods following Reviewer 2's comments. I feel these clarifications now clearly define both the advantages and limitations of the EMC2 method in comparison with previously developped solutions. By clearly discussing in which situations, such as Hydra neuronal tracking, EMC2 provides improvement from existing methods, the authors now clearly provide the reader with both a novel tracking method and rational arguments for using it, or not, depending on her/his own experimental data.

Reviewer #2: The authors have taken important steps to bring their claims more into alignment with the evidence they present, but some problematic language remains.

• The phrase “Long time lapse sequences” still appears on p. 4 This should be removed because it implies that the method will address the “long time lapse” case, when it does not.

• The placement of new language about “…premature ending of a track due to missing detection(s), can lead to important error propagation” on p. 8 is problematic. That text appears in subsection “2 - Limitations of Standard SPT Algorithms,” but the authors claim incorrectly that EMC^2 overcomes those limitations of standard SPT algorithms in subsection “3 – EMC^2 Algorithm”: “To overcome these limitations, we have developed EMC2…” The authors need to revise this framing so that they do not incorrectly claim that they have overcome the limitations on accumulation of errors over time.

• The claim (p.16) that EMC^2 is robust “independent of the type and complexity of particle motion” is not supported by the data. Indeed if this were true it would be extraordinary. I suspect that the performance of EMC^2, like all time-dependent tracking algorithms before it, is dependent on the amount and type of motion per frame, and with enough motion per frame its performance would degrade precipitously. Please qualify this claim or provide evidence to support it.

Minor, but important:

• This statement on p. 7 incorrectly implies that the sparse detectability of neurons is uniquely problematic for Hydra and not C. elegans or other organisms, when in fact this depends only on whether the animal is expressing a reference fluorophore like RFP, and is not an inherent property of the organism. Please reword and clarify.

“while the aforementioned mapping methods used in C. Elegans can accommodate moderate changes in the total number of neurons between frames due to missing or spurious detections (counting noise), the intermittent and sparse detectability of neurons in calcium imaging definitely hinders the applicability of mapping methods to Hydra [7]”

• Subsection “4 - Validation of EMC^2”: “We used the first 250 frames” p.13 The authors should specify in the main text the duration in seconds of these 250 frames (Perhaps 25 seconds?) This information is important for the reader to interpret the method’s performance.

• In the abstract the authors write: “Our results prove…”. The statement of proof is unnecessarily, and in my opinion, inappropriately strong. Best usage would be to reserve proofs for mathematical proofs. Please revise.

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org.

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

References:

Review your reference list to ensure that it is complete and correct. If you have cited papers that have been retracted, please include the rationale for doing so in the manuscript text, or remove these references and replace them with relevant current references. Any changes to the reference list should be mentioned in the rebuttal letter that accompanies your revised manuscript.

If you need to cite a retracted article, indicate the article’s retracted status in the References list and also include a citation and full reference for the retraction notice.

16 Jul 2021

Submitted filename: Cover_Letter_Revision_4.doc

Click here for additional data file.

8 Sep 2021

Dear Dr. Lagache,

We are pleased to inform you that your manuscript 'Tracking calcium dynamics from individual neurons in behaving animals' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

Boris S. Gutkin

Associate Editor

PLOS Computational Biology

Lyle Graham

Deputy Editor

PLOS Computational Biology

***********************************************************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #2: The revisions have adressed my concerns. Thank you.

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #2: None

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #2: No

4 Oct 2021

PCOMPBIOL-D-20-01324R4

Tracking calcium dynamics from individual neurons in behaving animals

Dear Dr Lagache,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Olena Szabo

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol