realSEUDO for real-time calcium imaging analysis

Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 14:25:07

Abstract: Closed-loop neuroscience experimentation, where recorded neural activity is used to modify the experiment on-the-fly, is critical for deducing causal connections and optimizing experimental time. Thus while new optical methods permit on-line recording (via Multi-photon calcium imaging) and stimulation (via holographic stimulation) of large neural populations, a critical barrier in creating closed-loop experiments that can target and modulate single neurons is the real-time inference of neural activity from streaming recordings. In particular, while multi-photon calcium imaging (CI) is crucial in monitoring neural populations, extracting a single neuron's activity from the fluorescence videos often requires batch processing of the video data. Without batch processing, dimmer neurons and events are harder to identify and unrecognized neurons can create false positives when computing the activity of known neurons. We solve these issues by adapting a recently proposed robust time-trace estimator---Sparse Emulation of Unused Dictionary Objects (SEUDO) algorithm---as a basis for a new on-line processing algorithm that simultaneously identifies neurons in the fluorescence video and infers their time traces in a way that is robust to as-yet unidentified neurons. To achieve real-time SEUDO (realSEUDO), we introduce a combination of new algorithmic improvements, a fast C-based implementation, and a new cell finding loop to enable realSEUDO to identify new cells on-the-fly with no "warm-up" period. We demonstrate comparable performance to offline algorithms (e.g., CNMF), and improved performance over the current on-line approach (OnACID) at speeds of 120 Hz on average. This speed is faster than the typical 30 Hz framerate, leaving critical computation time for the computation of feedback in a closed-loop setting.

Main paper: 
 introduction: 
 Closed loop experiments enable neuroscientists to modify or adapt presented stimuli or introduce perturbations, such as optogenetic stimulation, in real-time based on new observations of the neural activity. Such experiments are critical for both optimizing experimental time, e.g., by optimally selecting stimuli to fit neural response models (Charles et al., 2018) , or by deducing causality instead of only functional connectivity. Despite this critical need, closed loop experiments are incredibly difficult as they require real-time processing of neural data, which can be computationally intensive to process. In particular, population-level recordings using modern technologies often require significant computation to even extract individual neuronal activity traces, e.g., in high-density recordings or in fluorescence microscopy.One particularly challenging recording technology is fluorescence microscopy, in particular multiphoton calcium imaging (CI). CI has progressed significantly since its inception with optical advances enabling larger fields of view, and therefore higher data throughput. While neuroscientists now have access to hundreds-to-thousands of neurons at a time (up to even millions), the neuronal time traces are only indirectly observed as fluorescing objects in the full image. To extract each neuron's activity, a number of methods have been developed, including matrix factorization approaches, deep learning approaches, and others (we refer to a recent review for a more complete coverage of all the available methods and their nuances (Benisty et al., 2022) ).Almost all current calcium image processing methods rely on batch processing: i.e., using a full video all at once to simultaneously identify the neurons in the data and their time traces. For example, a common approach is to identify cells in a mean image (the image containing the average fluorescence per pixel over all time) and then to extract the time-trace from the video given the neuron's location. Real-time processing does not afford such luxury. Instead, frames must be processed as they are collected. Furthermore minimal data can be stored and used, as large image batches reduce algorithmic speed. Finally, the incomplete knowledge of the full set of cells in the video can cause unintended cross-talk as new cells firing may overlap with known cells, causing a well-documented effect of false transients (Gauthier et al., 2022) .We thus present an algorithm capable of demixing CI data frame-by-frame in real-time. Our design goals are to operate at > 30 Hz with minimal temporary data storage. Our primary contributions are: 1) A fast implementation of the SEUDO algorithm, 2) A new feedback loop to identify cells in real time using the SEUDO output, and 3) patch-based parallelization that enables high-throughput calcium trace estimation across larger fields of view. 
background: 
 Traditionally, CI analysis has been performed on full imaging videos. The goal of these algorithms is to extract from a pixel-by-time data matrix Y 2 R M ⇥T , where M is the number of pixels in each frame and T is the number of frames, a set of neural profiles X 2 R M ⇥N (one for each of N neurons) and a corresponding set of time traces 2 R N ⇥T . The former of these has, as each column of X, a single component profile depicting which pixels constitute that fluorescing object, and how strong that pixel is fluorescing. The latter has as each row the corresponding time traces that represent how bright that object was at each frame. These time-traces are particularly important for relating neural activity to each other via population dynamics, or to stimuli and behavior.In typical approaches, full videos are required to either 1) identify summary images (e.g., mean or max images (Pachitariu et al., 2013; Diego & Hamprecht, 2014; Petersen et al., 2018; Klibisz et al., 2017) ) to identify cells in, 2) to create a dataset within which points are clustered into cells (Kaifosh et al., 2014; Spaen et al., 2019; Mishne et al., 2018; Reynolds et al., 2017; Apthorpe et al., 2016; Soltanian-Zadeh et al., 2019; Bao et al., 2021; Soltanian-Zadeh et al., 2019; Kirschbaum et al., 2020) , or 3) to perform simultaneous cell identification and demixing (Pnevmatikakis et al., 2016; Pachitariu et al., 2016; Charles et al., 2022; Haeffele & Vidal, 2019; Inan et al., 2017; Maruyama et al., 2014; Mishne & Charles, 2019; Song et al., 2017; Giovannucci et al., 2019 ) (e.g., via matrix factorization or dictionary learning). For example, in the latter example, the data decomposition is solved via a regularized optimization, e.g.,EQUATIONwhere k•k 2 F is the Frobenius norm (sum of squares of all matrix elements), and R X (X) and R ( ) are regularization terms for the profiles and time-traces, respectively. While many regularization combinations exist, common terms include sparse neural firing, sparse overlaps, non-negativity, and spatial locality. Regardless, all methods require all frames to identify the fluorescing components, with the exception of OnACID (Giovannucci et al., 2017) and FIOLA Cai et al. (2023) .OnACID and FIOLA operate in an on-line manner, utilizing the buffer of last l b residuals r t = y t Xc t Bf t where X and c represent the spatial and temporal profiles of already recognized cells and B and f represent the spatial and temporal profiles of the known background signal. Both methods use a local Constrained Non-negative Matrix Factorization (CNMF) (Pnevmatikakis et al., 2016) in the spatial and temporal vicinity of that point. CNMF is an off-line algorithm that repeatedly performs alternating optimizations on [X, B] and on [c, f ] using the full dataset, until it converges to a designated precision. Both methods require initialization periods, and FIOLA further requires GPU and CPU optimizaiton, raising the computational infrastructure costs. We seek a solution that does not need any initialization data and can be run on simpler CPU machines for easier incporporation into user's workflows.One primary challenge in fully on-line settings is the incomplete knowledge of all fluorescing components at the experiment onset. Even in off-line methods, incomplete identification of cells can create scientifically impactful cross talk-termed false transients-in inferred activity (Gauthier et al., 2022; Inan et al., 2017) . Another challenge is identifying new components from few frames: ideally from individual frames to reduce memory usage. Recent work has provided an algorithm that solves both challenges: The Sparse Emulation of Unused Dictionary Objects (SEUDO) algorithm (Gauthier et al., 2022) . SEUDO is a robust time-trace estimator for neuronal time traces. Given a single fluorescence video frame y t , and a set of known profiles X, SEUDO models contamination as W c where W is a basis of small Gaussian bumps that linearly construct the interfering components, weighted by the sparse coefficients c (i.e., most c values are zero). SEUDO then solves the optimizationEQUATIONwhere and are model parameters and the internal min selects from the two internal expression that which has the minimal value. Since SEUDO operates per-frame, y t , P hi, c are all vectors. While SEUDO has demonstrated the ability to remove false transients (Gauthier et al., 2022), SEUDO's application has been limited to off-line post-processing due to: (1) slow computational speed, and (2) the need for pre-defined profiles X.The computational bottleneck in SEUDO is a weighted LASSO (Tibshirani, 1996) optimization, which can be implemented with the Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) that implements a momentum gradient descent (Beck & Teboulle, 2009) . FISTA optimizesEQUATIONwheref (x) is a smooth convex function with Lipschitz constant L > 0, such that krf (x) rf (y)k  Lkx yk8x, y 2 R n , e.g. in SEUDO f (x) = ky t X t W ck 2 2, and g(x) is a continuous convex function that is typically non-smooth, e.g., g(x) = kxk 1 . Each descent step of FISTA consists of an ISTA descent step and a momentum step:EQUATIONwhere the parameter ⌘ k gradually reduces with⌘ k = t k 1 t k+1 , t k+1 = 1+ p 1+4t 2 k 2 , t 1 = 1. 
real-time seudo: 
 Here we develop Real-time SEUDO (realSEUDO) that resolves the primary limitations of SEUDO and extends the algorithm significantly from a time-trace estimator to a real-time cell identification method. Specifically we improve the computationally intensive momentum descent algorithm used to solve Equation equation 2 by reducing the number of steps of momentum descent, implementing parallelism, optimizing internal computations (e.g., of smoothness parameters), and reducing the complexity of the original fitting problem without a substantial loss of quality by manipulating it inputs. Moreover we add a new algorithm that automatically recognizes the neurons that have not previously been seen and adds them to the SEUDO dictionary. Finally, we present a patch-based parallelism that avoids the slowdown of running LASSO in higher dimensions.realSEUDO 1 operates at high level as follows (Fig. 1 ): For simplicity we assume that we begin with no known profiles. When significant activity occurs, triggering a SEUDO cell, that cell is stored in a temporary array. Profiles are moved from the temporary profiles to the static set of profiles when they disappear from the current frames. The static profile set is then used to identify the activity of those components in future frames, with unexplained components becoming SEUDO cells and cycling back into the temporary profiles, followed by an update of the static profile set.Our first step in creating realSEUDO is to vastly increase the speed of SEUDO to operate at > 30 fps.The original implementation of SEUDO relied on the TFOCS library (Becker et al., 2011) in MAT-LAB. We began by improving the run-time by switching from MATLAB's interpreted programming language running TFOCS to a fast implementation of LASSO (Tibshirani, 1996) via FISTA (Beck & Teboulle, 2009) in the C++ compiled language. To further improve performance, we optimized the C++ code with the use of templates to eliminate the function call overhead in tight loops, and also employ parallelism, based on the POSIX threads with TPOPP library wrapper (Babkin, 2010).To prevent a bottleneck from the passing of data through Matlab's OOP API at the MATLAB/C++ interface, we switched to the non-OOP version of the MATLAB-to-C API.While beneficial, simply coding SEUDO in C++ was not sufficient to achieve the desired processsing rate. We further improved runtimes with optimizations to the computation of the cost function and derivatives. The partial LASSO component of SEUDO that performs optimization on a single frame y using the FISTA algorithm can be written as arg min f ( ) + g( ), where f = ky k 2 2 is the least-squares portion and g = k k 1 is the `1 penalty. In the FISTA algorithm the spaces of time traces and Gaussian kernels are unified, with representing a concatenation of X and W , representing a concatenation of and c. With the number of pixels in the frame M , number of neurons N , and number of Gaussian kernels K,y 2 R M , 2 R M ⇥(N +K) , 2 R N +K , 2 R N +K ,with the first N elements of corresponding to being equal to 0.In working with this optimization program, a number of internal computations become bottlenecks: specifically computing the gradients r f and r (f + g), the Lipshitz smoothness estimation, and the momentum/stopping criteria.We reduce the burden of the gradient computations by both reducing the number of times the gradient must be used, and by improving the internal gradient computation. For the former, we note that naive implementations compute both a step in the direction of r f and then in the direction of r (f + g). Moreover, we noted that these two steps in slightly different directions cause the gradient to dither around the optimum. We thus instead only take a step in the direction of r (f + g) (similar to (Beck & Teboulle, 2009) ). For the latter, computing r (f + g) requires matrix vector multiplications with and its transpose. Since is sparse, to save memory and computation, we avoided storing the matrix explicitly, instead generating it on the fly in an efficient convolution. Specifically, the gradient r f requires computing T , which we reorganize to compute in two passes: The first pass:EQUATIONthe first computes intermediary variables, and the second uses these values to compute the gradient dimensions.This split into two passes factors out repeated computations, and reduces the time complexity from O(n 3 ) to O(n 2 ). Moreover, the computation of each pass becomes highly parallelizable, with partitioning of the first pass by j and of the second pass by m, and efficiently skipping the iterations over the zero elements in the sparse matrix .To improve the efficiency of estimating the Lipschitz constant L, note that L = max(krf (x 1 ) rf (x 2 )k/kx 1 x 2 k). We found that for our cost function, we could approximate L with independent computations in each dimension. This estimation reduced the number of steps by as much as 30% over typical computations of L before each step based on the local gradient.The momentum descent central to solving the partial LASSO tends to spend many steps on stopping the momentum, especially with the large values of the Lipshitz constant L (i.e. the non-momentum steps are small). One such case is "circling the drain" around the minimum, with the momentum causing the overshoot in one dimension while another dimension is stopping. Another case is when a dimension is moved past the boundary (e.g., x 0 for SEUDO). The momentum continues to push the solution to negative values, producing suboptimal solutions and increasing the number of steps necessary. FISTA includes a parameter ⌘ that progressively limits the top speed of descent to reduce such problems. We improved these cases by resetting the momentum to zero on a dimension when it either attempts to cross into the negative values or when its gradient changes sign. The dimensional momentum stopping stops abruptly at the right time, obviating the need for slowing and thus we can simplify FISTA by fixing ⌘ = 1.We have also experimented with using the dimensional momentum stopping on a different kind of optimization problems, on momentum optimization of the neural network training, where it provided a substantial improvement. While FISTA performed worse than simple gradient descent, modified FISTA produced the same error rate and squared mean error in about 10 times fewer training passes than simple descent, or about 10 times lower error rate and 1.5 times lower mean square error in the same number of passes. The FISTA algorithm incorporating all these extensions is presented in Algorithm 2 in Supplement.As a third step, we modified the optimization program (the SEUDO model) itself to achieve the final speedups. In the original SEUDO, Gaussian components in W are spaced apart by one pixel. We instead found that this tight spacing is redundant. The kernels with radius r cover (⇡ ⇤r 2 ) pixels, and thus each pixel is covered by (⇡⇤r 2 ) kernels with varying weights. This redundancy results in FISTA continuing to adjust the kernel coefficients c after the neural activations x converge. Reducing the number of kernels thus reduces both the number of gradient descent steps and the per-step cost, accelerating the computation more than quadratically. For a kernel with diameter of 30 pixels this improves the performance by a factor of over 100 without substantial degradation of false transients removal or the recognition of the interfering components' shape W c.We tested our speedups against the original SEUDO on 45000 frames across 50 cells from Gauthier et al. (2022) . SEUDO ran at 5.8-6.9 s/cell on a Macintosh M1. The optimized C++ implementation without the MATLAB C++ API reduced the runtime to 0.9-1.1 s/cell (a 6-7x improvement). Sparse SEUDO provided further acceleration to a run time of 0.2 s, with 0.1 s for the computation and 0.1 s for the overhead of converting data between Matlab and native code; a total of a 29-34.5x speed-up.With our improvements, the core SEUDO algorithm was able to run in real time, and we next moved on to the cell recognition feedback loop that forms the realSEUDO algorithm. As realSEUDO must rely only on past video frames to identify new profiles, we designed realSEUDO to run on a frame-by-frame basis. At a high level, our automatic cell recognition first runs SEUDO on the current incoming (denoised) frame given the known profiles X stab identified thus far. The loop then identifies contiguous bright areas in the residual frames, i.e., the SEUDO cells resulting from SEUDO, and places them in a 'temporary profile' array X temp . These temporary profiles are then updated (via merging with new potential profiles) given new, incoming, frames until they are stable and moved to the stable, known profile list X stab that is updated less frequently by addition, merging and splitting of temporary profiles.Procedually, we first preprocess each incoming frame with the goal of reducing noise to improve profile recognition. Typical calcium imaging analysis often uses running averages in space and time for noise reduction. We thus implement both a spatial Gaussian filter, as well as a running average of several sequential frames. We keep the window length as a tunable parameter that can be set to one to perform frame-by-frame processing and minimize temporal blurring of neural activity.After each frame is denoised, we estimate the activation level for each of the stable profiles X stab using the fast version of the SEUDO estimator. We denote the activation level for the k th profile at time t as kt . SEUDO also returns a robust residual that contains all the fluorescence that was not captured by X stab . We then run the residual through SEUDO a second time using the temporary profiles X temp to test if any temporary profile represents the fluorescence well and should be moved to X stab from X temp . The residual after this second application of SEUDO represents completely unknown profiles and are analyzed separately to determine if a new member of X temp should be created.Starting from the bottom up, the detection of new temporary profiles is based on finding the areas of the image that stand above the noise level. The noise level is evaluated by noting that most of each video frame has no activity, indicating that the median pixel value will be very close to the median value of the pixels in an all-dark frame containing the same noise. The half-amplitude of the noise ha can be estimated as: ha = median(y t ) min(y t ).(7) In some cases different areas of the image may contain a different amount of background lighting (e.g., changes in neuropil), which can skew the noise estimate. To overcome this challenge, we split the larger image into sections, with each section computing a local median which is smoothly interpolated between section of the image. This can be thought as either a krigging procedure or a cheap approximation to a local median evaluated independently for each pixel in the image.To merge new profiles into the exiting profiles when adding new profiles to X temp , we compute an overlap score. The overlap computation is not a plain spatial overlap but includes a heuristics that identifies when one profile is mostly contained in another. The condition for merging two profiles in X temp is based on the comparison of numbers of common and unique pixels between profiles, where P 1 and P 2 are numbers of pixels in each of two profiles, U 1 and U 2 are the numbers of unique pixels in each profile, C is the number of common pixels, B 1 and B 2 are the perimeters of bounding boxes for each profile, and k temp is a constant with an empirically chosen value of 0.75:U 1  B 1 ⇤ 0.5 or U 2  B 2 ⇤ 0.5 or C k temp ⇤ min(P 1 , P 2 ) (8)After a profile is moved from X temp to X stab , a different score is computed pair-wise between the new profile and each existing profile, to decide whether they should be left separate, or merged, or one of profiles split. If a merge or split is performed, the original profiles are removed from X stab and the results entered recursively into X stab as new profiles. The score for two profiles A and B in X stab is computed based on the brightness and measures of least-squares fit of the cells into each other 1) as whole cells (i.e., ↵ AB = hA, Bi/hA, Ai) and 2) using only the overlapping region (i.e., AB = hA ol , Bi/hA ol , A ol i where A ol is the profile A restricted to the region overlapping with B).AB is used as a measure of difference in brightness, against which the fit of the whole cells ↵ AB is compared as a measure of proximity in shape. Specifically we compute two ratios ⇢ AB and ⇢ BA are computed asEQUATIONA higher value of ⇢ AB (always  1) means that cell A fits better inside cell B. The value of 1 means that it fits entirely inside cell B. The same principle applies symmetrically to ⇢ BA .The scores express the following logic: If two temporal profiles are a close match in cross-section, they likely represent the same cell and should be merged. If they overlap only partially, they likely represent separate cells. If one cross-section is inside the other, look at relative brightness: if the smaller cross-section is also weaker, it's likely a weaker partial activation of the same cell and should be merged, if the smaller cross-section has a close or higher brightness, the larger cross-section likely represents an intertwining of two cells that has to be split.realSEUDO, although highly efficient for smaller patches, is still based on the LASSO algorithm that reduces in efficiency with much larger frames. Thus we adopt a patching scheme that breaks each frame into small patches that can be parallelized to maintain the high framerate by utilizing the multithreading in many modern processors. Patching, however, requires matching profiles across patches.Traditionally profiles discovered in data split into patches is to add overlap margins to the patches and to use profiles overlap in this region to determine matchings in neighboring patches. Additional margins, however, introduces redundant computation and decreases computational efficiency.In realSEUDO we note that the logic behind the scoring we use to merge profiles in X stab and X temp within each patch can also be used to score the match of profiles across patches. Specifically, we extend the matching to include the profile temporal activity as an additional dimension to find matching cells in neighboring patches via consecutive gluing of the profiles. The highest score is assigned to the bidirectional match of both spatial and temporal dimensions, a lower score to symmetrical match of spatial dimensions and asymmetrical match of temporal dimension, a yet lower score to an asymmetrical match of both kinds of dimensions. We have observed successful matches even with zero margin, using the neighboring strips of pixels around the perimeters of the patches as the spatial dimension for matching.We first applied all three algorithms to a simulated video created with Neural Anatomy and Optical Microscopy simulation (NAOMi) (Song et al., 2021) . Specifically we simulated the neural activity over 20000 frames at 30Hz with fame size of 500x500 pixels. There are approximately 450 cells visible in this dataset (i.e. fluorescing cells intersecting the plane of imaging). We benchmarked the patch-based parallel processing of realSEUDO 80x80 pixel patches.For comparison we ran the off-line CNMF (a staple batch-based calcium imaging demixing algorithm) and OnACID, the computationally similar on-line method. realSEUDO found 201 true cells, identified as strongly correlated with ground-truth cells, while OnACID found 152 and CNMF (the offline method) found 308 (Fig. 2A-B ). Furthermore, both realSEUDO and OnACID found many fewer false positives than CNMF, presumably because they cannot be fooled by small fluctuations integrated over the full recording (Fig. 2C-D ). Note that to remove the confound of post-processing we followed prior work Song et al. (2021) in using the CNMF raw fluorescence traces instead of the model-based denoised traces. On average, realSEUDO processed 67.8 frames per second endto-end, while OnACID ran at 9.2 fps.Applications to in-vivo mouse CA1 recordings We applied realSEUDO to an in vivo calcium imaging recordings from mouse hippocampal area CA1 previously described in Gauthier et al. 2022 (Gauthier et al., 2022) . They consisted of 36 videos, each sized 90x90 pixels with 41750 frames sampled at 30 Hz. The outputs had previously been verified manually by Gauthier et al. 2022 (Gauthier et al., 2022) with human labeling of CNMF outputs. We applied realSEUDO to all videos and compared the outputs with the current online cell demixing algorithm OnACID (Giovannucci et al., 2017) , as well as a popular offline algorithm, CNMF (Giovannucci et al., 2019) , as an additional baseline.We benchmarked realSEUDO against OnACID (an online analysis tool) and CNMF (an offline analysis tool) on real in-vivo calcium imaging movies (Fig. 3 ). Algorithmic performance was measured on an x86-64 computer with 48 CPU cores (Intel Xeon 6248R), 2 hyperthreads per core, 78 GB of memory, and without the use of a GPU. The initialization times were not included. On average, realSEUDO processed 162 frames per second compared to 26 processed by OnACID and 13 by CNMF: an improvement of 6.5x and 12.5x respectively (Fig. 3C ). Quality-wise, we found that OnACID exhibited difficulties with adapting to larger ranges of pixel brightness, sometimes missing bright cells. Scaling pixel values improved OnACID results, but only mildly (one additional cell). Numerically OnACID and CNMF appear similar but they identified different components. SEUDO results were most similar to CNMF, and with additional cells identified, and less false positives (Fig. 3B ). Finally, the per-transient manual classification provided by (Gauthier et al., 2022) enabled us to assess if realSEUDO inherited the false transient removal properties of SEUDO. For a reasonable value of = 0.15, realSEUDO had a true positive rate of 75% and a false positive rate of 24%. While these numbers are a bit lower than the numbers reported in (Gauthier et al., 2022) , in that study the authors average N=3 frames to reduce noise, while maintained single-frame analysis. This can be evident by the fact that missed transients were very small: realSEUDO kept 98% of real fluorescence and only 15% of false fluorescence.Additional in-vivo tests: As final test we applied all three algorithms (realSEUDO, OnACID, CNMF) to a 2000-frame mesoscope video example collected by the Yuste lab at Columbia Uni- versity and provided with the OnACID github package as a demo. For this example we similarly saw improved cell detection and runtime improvement in terms of fps (Fig. 3B ). 
discussion: 
 We present here an online method for cell detection and fluorescence time-trace estimation from streaming CI data: realSEUDO. realSEUDO is based on the SEUDO robust time-trace estimator that reduces bias due to unknown cells while also providing approximate shapes of the unknown fluorescing objects. To build realSEUDO we 1) improved SEUDO's slow runtime via significant modifications at the code, algorithmic, and model levels 2) built a new feedback loop that allowed SEUDO (which has no cell finding component currently) to identify cells in real-time. Overall, realSEUDO can achieve frame processing rates of 80-200 fps, depending on cell density. While our goal was to exceed the typical 30 Hz data collection rate common to many experiments, the high processing efficiency leaves additional time to compute feedback in future closed loop systems. Moreover, realSEUDO can scale with faster recording rates as calcium indicators become faster, e.g., GCaMP8 (Zhang et al., 2023) . realSEUDO's implementation exhibits a higher degree of parallelism than OnACID, however both are likely constrained by the employed tools. Specifically, the measurements in Fig. 3C show very little fps fluctuation with respect to cell count for OnACID, to the point where the error bars are invisible, which might be due to a bottleneck in a single thread. realSEUDO has a lower latency to the first events for a new cell, OnACID requires a history of 100 frames to recognize a cell, while realSEUDO would produce the first events starting with the first frame that passes the low brightness threshold. Alternatively, OnACID has a higher native scalability with respect to the frame area, conditioned on similar cell counts. This likely result from OnACID's algorithm restriction of processing to areas in the immediate vicinity of known cells. realSEUDO can still achieve a high processing efficiency with reasonable hardware with our parallelization. Future research may further improve realSEUDO runtime by blanking out entire patches until activity is detected via simpler activity detectors.One strength of realSEUDO is that it can be initialized with either an empty profile set or a precomputed set, e.g., via prior analysis or via anatomical (z-stack) imaging. The ability to start from nothing will be useful in mesoscope settings when the field-of-view may be changed on-the-fly, and not requiring an initialization step will reduce start-up overhead. Furthermore, we note that many of our speed adaptions deviate from the traditional gradient descent approach. In particular the approximation of the Lipshitz constant and the changes in the momentum and stopping criteria. We have observed that in other domains, in particular for training deep neural networks, these deviations appear to also provide significant speedups. We thus will consider future research directions in analyzing and quantifying the extent of these improvements across broader applications.Limitations: While our results achieve the design criteria we initially set out, there are some potential barriers. For one, as with many real-time systems, the compute environment is very important to configure correctly. We have found the importance of explicitly setting the Linux CPU manager to enable the performance mode. The default automatic adaptable mode does not react properly to CPU loads of less than 100% of the whole capacity, and significantly skews the benchmarking by running the CPUs at low frequency.While our core algorithm is written completely in C++, and thus open source, we have found MAT-LAB convenient and efficient as a wrapper for prototyping wrappers for our core functions. Further work will add Python wrappers to allow for seamless integration into both MATLAB and Python pipelines, enabling realSEUDO to be more widely used. Finally, we focused here only on cell detection, assuming access to the on-line motion correction algorithm from OnAcid. Future work should more holistically incorporate the full pipeline in C++. 
 
 Appendix:  


Decision: This paper extends a previously developed algorithm called SEUDO for extracting time traces from two-photon imaging data robustly also in challenging real-time close-loop setting. The algorithm can simultaneously extract ROIs and time traces. The algorithm performs competitively with state-of-the-art offline algorithms and better than state-of-the-art online algorithms. The results and concepts are well explained, broadly accessible and potentially very impactful for closed-loop neuroscience experiments.

Reviews: review: paper_summary: paper_summary: The paper presents a novel online method for cell detection and fluorescence estimation from streaming calcium imaging (CI) data, called realSEUDO. Building upon the SEUDO algorithm, realSEUDO introduces significant advancements in real-time cell identification and fluorescence estimation. The authors effectively address the challenge of high-speed processing by leveraging parallel computing and profile matching techniques across patches and temporal dimensions.
strength_weakness: strength_weakness: While the paper presents a novel and efficient approach with realSEUDO, several limitations are apparent. The reliance on specific hardware configurations, such as Linux CPU performance modes, could limit accessibility and generalizability. Additionally, the implementation primarily in C++ with MATLAB wrappers, while effective, could benefit from broader integration into Python pipelines to enhance usability. The paper also assumes access to external motion correction algorithms and does not address the full pipeline holistically, which may impact its applicability in scenarios without these additional tools.
questions: questions: N/A
limitations: limitations: The paper presents notable advancements with realSEUDO but has some limitations. Firstly, the performance of the system is highly dependent on specific hardware configurations and the proper setting of Linux CPU manager modes, which may not be universally applicable. Additionally, the reliance on MATLAB for prototyping, combined with limited Python integration, may restrict its accessibility and ease of use for a broader audience. Moreover, the paper focuses on cell detection without fully integrating or addressing the motion correction and complete processing pipeline, which could limit its applicability in real-world scenarios where such components are not readily available.
paper_summary: paper_summary: The authors present an online method for neuron profile and trace estimation in multi-photon calcium imaging. Their approach enhances the SEUDO algorithm by running it on a frame-by-frame basis and periodically updating profiles using a patching mechanism across space and time. The implementation leverages different degrees of stability, namely temporary and stable. The proposed method is benchmarked against both an online method and an offline method using simulated and real data.
strength_weakness: strength_weakness: - The work largely relies on previously described algorithms like FISTA and SEUDO. Since a significant portion of the method describes these algorithms in detail, it is unclear what the key contributions of this work are.
- The evaluation was conducted on synthetic data (Naomi) and a real dataset (36 videos). Including additional sources like 2p movies from the Allen Brain Observatory would have enhanced the credibility of the results.
- The given metrics should be improved. While the images are good, providing tables with quantitative metrics can enhance the overall presentation and support the claims made. In addition to correlation, metrics like IoU could be evaluated on neuron profiles.
- The experimental setting lacks an ablation study on the main components of the work, such as automatic cell recognition (e.g. with and w/o X_temp) and patching.
questions: questions: 1. What are the key methodological differences between SEUDO and realSEUDO?
2. How did you choose the benchmark datasets?
3. Why correlation was selected as the main (and only) metric?
4. How do the authors justify the lack of an ablation study on the main components of the work? Can the authors describe the quantitative improvement of their method compared to existing methods they rely on?
limitations: limitations: The authors have adequately addressed the limitations.
paper_summary: paper_summary: In this study, the authors created a real-time version of the “Sparse Emulation of Unused Dictionary Objects” (SEUDO) algorithm termed realSEUDO, that can process multi-photon calcium imaging (CI) data in real time. The goal of this approach is to overcome the limitations of batch processing methods commonly employed in CI analysis by processing brain activity data in real-time, thus enabling closed-loop neuroscience investigations. The authors identified several places that the original SEUDO algorithm warranted improvements, implemented the changes and compared their method to other current online methods, with both simulated and real datasets.
strength_weakness: strength_weakness: While the performance with a small number of neurons is impressive there are two drawbacks. 
First, the testing datasets only include a very limited number of neurons and the scalability with much larger datasets are not presented. For example, current methods such as with a mesoscope  (Sofroniew et al. 2016), allow for much larger recording samples in the of thousands of neurons. From the figure 3c we can see that the algorithm drops in performance quite fast and when much larger samples are included, it could potentially perform worse than existing methods, something that remains to be tested. 
Lastly, calcium sensors already limit our ability to extract with fine temporal precision the neural activity thus the usefulness of a method that far exceeds that temporal resolution remains limited. 
Another potential issue is to validate the model with more typical cortical recordings that have higher level of background fluorescence, and how it deals with overlapping sources of fluorescent signal.
questions: questions: By facilitating the real-time processing of CI data, the paper shows the potential to improve the accuracy and efficiency of closed-loop studies. The claims of broad applicability would be strengthened by more extensive testing under various experimental situations and with diverse kinds of brain data. Additionally, it would be great if they could characterize and show the performance with different patch overlaps. Finally, while the paper discusses the integration with motion correction algorithms from OnACID, integration of all preprocessing steps such as motion correction within the realSEUDO pipeline would further widen its potential usability.
limitations: limitations: The authors have a limitation section that addresses potential issues with their algorithm.
paper_summary: paper_summary: The authors extend an available algorithm called SEUDO for extracting time traces from two-photon imaging data robustly to the real-time close-loop setting. The algorithm can simultaneously extract ROIs and time traces. The algorithm performs on par with state-of-the-art offline algorithms and better than state-of-the-art online algorithms.
strength_weakness: strength_weakness: * The main improvement come from engineering the implementation of SEUDO such as to run in close to real time (is that a weakness or strength?).
* Lines 228: The merging procedure seems quite ad-hoc and heuristic. The logic is only laid out at the end (lines 247), maybe it would help to discuss this earlier.
* Only two algorithms are chosen for comparison. 
* Fig. 2B is very hard to read and understand. Does this mean, OnACID produces no false alarms and no ambiguous units?
questions: questions: * A lot of engineering work has gone into the fast implementation of the online algorithm using C++ and various bells and whistles. Have the authors considered a GPU based implementation e.g. in jax and how would that compare in performance?
* Where does equation 7 come from? Min operations are quite sensitive to outliers. The details won’t matter, but in the spike sorting literature, a threshold of n*0.6745(median (abs(x))) is often used for detecting spikes, which is indeed a multiple of a robust estimate of the noise sd without the spikes.
* Why was the tool not compared to the other online algorithm FIOLA?
* Why was the plain CNMF algorithm used and not one of the more advanced implementations/methods or an additional algorithm?
* Doesn’t Fig. 2C speak very clearly for the results of the OnACID algorithm?
* Lines 348: Why do the authors say that a c++ implementation makes the code open source?
limitations: limitations: Adequately adressed.

Comments: title: Thank you
Thank you
comment: Thank you for pointing us to the availability of the raw Allen Visual Coding - 2P dataset on AWS. We will certainly be making use of this resource as we continue to assess and develop new methods for 2P data processing. We are also thankful that our clarifications were helpful in understanding our choices in metrics and will make sure these points are clear in the final manuscript.
I would like to thank the authors for the response to my comments and the clarifications that will appear in the text.
I thank the authors for their clarification. Regarding the evaluation, I acknowledge the authors' efforts in choosing the right datasets. However, I would find the results more convincing if another source of real data were used. In this regard, the Allen Visual Coding - 2P dataset, in addition to extracted traces and ROIs, also provides raw movies on AWS. 

Aside from this, I also acknowledge that all technical aspects were well addressed, and I find the authors' motivation for not using segmentation metrics convincing. I will update my original assessment to reflect my new perspective.
Dear reviewer,
Thank you and I am glad that we were able to clarify! Your point about the 5th percentile is well taken and we will definitely test that out.
Thanks for the clarification. I maintain my score, and would vote clearly in favor of "accept".

Regarding the use of min: Consider using the 5th percentile or similar.