1 Introduction
In the past few years, calcium imaging has been increasingly adopted in neuroscience research since it allows simultaneous measurement of the activity of a large population of neurons at the singleneuron resolution over weeks using optical imaging of living animals. For example, we recently conducted a longitudinal investigation of the neural ensemble dynamics of contextual discrimination by recording mice’s calcium imaging in the hippocampus for about 60 days [Johnston et al., 2020] using appropriate genetically encoded calcium indicators [Tian et al., 2009] and miniature fluorescence miscroscopes [Ghosh et al., 2011]. The technical advancements bring great flexibility to neuroscience research; they also create significant challenges to every aspect of data analysis  from storing the large amount of video recordings to downstream statistical modeling and inference [Pnevmatikakis, 2019]. In longitudinal recordings of freely behaving animals, after motion correction and registration of detected neurons across multiple sessions [Giovannucci et al., 2019]
, fluorescence traces of individual neurons can be extracted, for example using independent component analysis
[Mukamel et al., 2009] or nonnegative matrix factorization methods [Maruyama et al., 2014, Pnevmatikakis et al., 2016, Zhou et al., 2018]. Strategies to reduce false positive rates of neuron detection have also been proposed, such as our recent work on imposing spatial constraints on footprints of neurons [Johnston et al., 2020]. The extracted fluorescence traces of individual neurons allow for the statistical analysis of complex neural interactions at the single cell level. Fluorescence traces are often used for visualization, clustering neurons based on similar neural activity profiles, comparing neural activity levels between various experimental conditions, and studying neural encoding and decoding.
However, calcium fluorescence data provide a noisy proxy, rather than direct observations of neural activity. For many important questions, such as those involving the analysis of the precise timing of neural activity in response to stimuli, it is essential to estimate the underlying spike train from a noisy fluorescence trace. Several approaches have been developed including linear deconvolution [Yaksi and Friedrich, 2006] or nonnegative convolutions [Vogelstein et al., 2010]. Fully Bayesian methods have also been developed to obtain posterior distributions of spike trains, thus allowing uncertainty quantification of spikes’ estimates [Pnevmatikakis et al., 2013]. [Theis et al., 2016]
proposed a supervised learning approach based on a probabilistic relationship between fluorescence and spikes. This algorithm is trained on data where spike times are known. Research has shown that nonnegative deconvolution outperforms supervised algorithms and is more robust to the shape of calcium fluorescence responses
[Pachitariu et al., 2018]. One limitation of these methods is that these spike detection methods analyze only one trial at time; thus, shared information across trials is largely ignored. In a longitudinal study, neural activities are measured in multiple trials or sessions. In these settings, aggregating information across trials could increase the accuracy of spike detection.There has been limited work on how to utilize the information from multiple trials to improve accuracy in spike detection. Among the few multitrial methods we identified, Picardo et al. [2016] assumed that repeated trials share the same burst time but have trialspecific magnitude, baseline fluorescence, and noise level. Conceptually, integrating multiple trials should be beneficial if the trialtotrial variation is mainly due to randomness. In reality, however, as pointed out in Deneux et al. [2016], the gain by naively combining trials may be limited  it is likely to bring improvement for some trials but might perform worse in others when the neural activity in some trials does not follow the marginal pattern across all trials. During a learning process or when adjusting to a new environment, neurons constantly reorganize and show plasticity, which leads to varying neural dynamics across trials.
In this paper, we develop a robust multitrial spike inference method to address the challenges from longitudinal calcium imaging data. Our approach both integrates the commonality and accounts for evolving neural dynamics across trials. Efficiency is achieved by aggregating information from temporally adjacent trials whereas robustness is guaranteed by incorporating a timevarying firing rate function into a dynamic penalization framework, rather than forcing shared parameters across trials.
2 Methods
In this section, we first propose a timevarying penalized framework to analyze multitrial calcium trace data. We then present details on how to implement it by alternating a firing rate estimation step and a spike detection step.
2.1 MultiTrial timevarying penalized autoregressive model (MTVPAR)
Let be the fluorescence recorded at time point of trial , where , . In practice, multiple preprocessing steps, including normalization, [Vogelstein et al., 2010], are typically implemented to obtain . The calcium fluorescence trace is often modeled using the following firstorder autoregressive model
(1) 
where and denote the underlying calcium concentration and noise at time in the th trial, respectively; is the decay parameter for a calcium transient, which depends on multiple factors such as sampling frame, cell types, and the kinetics of genetically encoded indicators. In this model, denotes the change in calcium concentration between time and time with indicates a spike occurs at time and 0 means no spike. The main goal of spike detection to locate the time points with a positive .
Here, we propose to regularize the inference on spike detection by introducing a timevarying penalty function on the number of spikes, as the spikes tend to be sparse but not evenly distributed over time:
(2) 
We further assume that is a decreasing function of the instantaneous firing rate , which is assumed to be smooth over time within trials and change slowly over trials . Because the parameters in the AR(1) model are trialspecific but the firing rate functions across trials are assumed to change slowly to capture the evolving neural dynamics, efficiency and robustness are well balanced. In addition, our modeling approach allows the simultaneous estimation of spike trains and firing rate functions.
In the rest of this subsection, we present an iterative approach that (1) estimates the firing rate function, given the spike trains derived in the previous iterations, by using a local nonparametric method; (2) detects spike trains, given the current estimate of the firing rate function, by using a timevarying penalization.
2.2 Firing rate estimation and timevarying penalty
Estimating the firing rate function is a crucial task in the analysis of spike train data [Cunningham et al., 2009]. One commonly used descriptive approach involves building the peristimulus time histogram (PSTH) where the spike counts are averaged from multiple trials within each time bin [Gerstein and Kiang, 1960]. Kernel methods are often applied to achieve smoothness [Cunningham et al., 2008]. Bayesian methods have also been considered, e.g. by proposing Gaussian processes [Cunningham et al., 2008] and Bayesian Adaptive Regression Splines [DiMatteo et al., 2001, Olson et al., 2000, Kass et al., 2003, 2005, Behseta and Kass, 2005].
The PSTH approach and the other smoothing methods assume that the underlying firing rate function does not change over trials. The estimated firing rates are then compared between different groups to capture the association between firing rate and animal behavior [Jog et al., 1999, Wise and Murray, 1999, Wirth et al., 2003]. However, recent evidence suggests to move beyond the independent and identically distributed trial assumption and regard both neural and behavioral behavioral dynamics as smooth and continuous [Huk et al., 2018, Ombao et al., 2018]. Thus, it is essential to model the betweentrial dynamics. To account for the betweentrial dynamics, some authors have proposed a statespace framework [Czanner et al., 2008, Paninski et al., 2010]. In the statespace framework, the spike train is characterized by a point process model [Brown et al., 2003, Brown, 2005, Kass et al., 2005, Daley and VereJones, 2007] for the underlying fire rate. In this paper, to make it feasible to jointly estimate the spike trains and the firing rate function from the observed multitrial fluorescence data, we consider instead a computationally less demanding twodimensional Gaussianboxcar kernel smoothing function , which is formulated as follows
(3) 
where denotes the withintrial kernel bandwidth in the Gaussian kernel, is the indicator function, and is a bandwidth for the betweentrial sliding windows in the boxcar kernel. Thus, our estimate of is given by
where is the estimated number of spikes per second in a small time bin centered at time in trial .
For PSTH based on spike train data, common choices of the Gaussian bandwidth are around 50150 ms [Cunningham et al., 2009]. When choosing the optimal bandwidth to smooth a spike train estimated from calcium fluorescnce trace data, one should also consider the fluorescence transient kinetics. Ali and Kwan [2019] reviewed the kinetics of calcium transients with several fluorescent indicators. For example, they reported that GCaMP6f has a rise time of 42 ms and decay time of 142 ms while GCaMP6s has a rise time of 179 ms and decay time of 550 ms; a timebin of 200 ms is suggested in one of their analyses. In our analysis, is chosen between 200ms and 400ms. For the boxcar bandwidth to allow borrowing information across trials, we use a preselected window size , which is mainly determined by to the number of available trials. As discussed in the Discussion section, data driven methods, such as the one used in [Fiecas and Ombao, 2016, Ombao et al., 2018], will be considered in future work.
In spike detection, sparsity has been enforced via penalization or introducing appropriate prior distributions. Existing approaches usually adopt a tuning parameter uniformly for the whole time series of a fluorescence trace. Within a trial, the firing rate right after each stimulus are expected to be higher than baseline. Thus, using a constant penalty may not be optimal. Our proposed nonconstant penalization is inspired by prior work in the literature. For example, timevarying penalization was used for analyzing multivariate time series data [Fan et al., 2013, Yu et al., 2017]. Zbonakova et al. [2016] studied the dynamics of the penalty term in a Lasso framework for the analysis of interdependences in the stock markets. Monti et al. [2017] used varying regularization for different edges in a Gaussian graphical model to study brain connectivity in fMRI data.
Here, we expand on those approaches and consider a decreasing function of the firing rate function for the timevarying penalty, motivated by the fact that spikes are expected to be less frequent in regions with low firing rates than in regions with high firing rates. Ideally, the penalty should be small in the locations with higher firing rates. Hence, we use a negative exponential function [Tang et al., 2010] to achieve adaptive regularization. Specifically, the timevarying function is chosen as
where is the estimated firing rate (See Algorithm 2). Here the value controls how much the penalty function should depend on the firing rate function. In particular, reverts to a constant penalty. To avoid extremely large or small penalties, we scale the estimated firing rate function of each trial by dividing by its maximum value, so that it ranges between 0 and 1. This implies that the default value in our analysis leads to mildly timevarying penalties  within a trial, the penalty at the highest firing rate is about of the penalty at a firing rate of . In addition, to facilitate comparison with the case of a constant penalty, the following formula is used in each trial to scale the penalty function to have mean ,
(4) 
2.3 Timevarying penalized AR(1) model
In our proposed multitrial timevarying penalized autoregressive model (MTVPAR), a calcium fluorescence trace at the th trial , is modeled by a firstorder autoregressive model, as given in Equation (1). As previously stated, denotes the underlying calcium concentration and a positive value of implies that a spike occurs at time .
Because spikes tend to be sparse, ideally, one should penalize the number of spikes by introducing an penalty. However, penalization is computationally intractable; therefore, existing methods often impose an penalization [Vogelstein et al., 2010, Friedrich and Paninski, 2016, Friedrich et al., 2017]. Recently, Jewell and Witten [2018] found that, for spike detection, the use of an penalty brings substantial improvements over an penalty. They also showed that relaxing the nonnegative constraint of has a negligible effect on the results but the corresponding optimization problem is equivalent to a change point detection problem whose solution can be obtained efficiently using a dynamic programming algorithm. Use a similar strategy, we prove that the following timevarying penalization problem can also be solved by a dynamic programming algorithm:
(5) 
Specifically, we find that the timevarying optimization problem is equivalent to the following change point problem whose solution can be efficiently identified using a dynamic programming algorithm. For the ease of presentation, we drop the trial index and focus on the fluorescence trace of a cell at a given trial. The equivalent change point problem can be shown as (See Appendix A.1 for proof):
(6) 
where are change points, i.e., the points satisfying and
(7) 
which has the following closedform solution:
(8) 
Note that the parameter , which measures the speed at which the calcium concentration decays, is not estimated. The value of is usually close to 1 [Vogelstein et al., 2010, Yaksi and Friedrich, 2006], as a somatic calcium transient caused by an action potential is often characterized by an almost instantaneous rise but a slow decay. For computational feasibility, rather than estimating iteratively, similar to [Pnevmatikakis et al., 2016, Friedrich et al., 2017], we estimate it using the autocorrelation at lag 1.
2.4 Algorithms
Thus far, we have presented our solutions to two problems separately. The first problem is to use multitrial spike data to estimate the firing rate function for and using a Gaussianboxcar smoothing method. The second problem is to detect spikes from a single calcium fluorescence trace using a timevarying penalized approach under the assumption that the timevarying penalty function
is already known. In practice, neither firing rate functions nor spike locations are known. We therefore propose an ExpectationMaximizationtype algorithm to jointly estimate firing rate and detect spike locations. As shown in Algorithm 2, we alternate between the spike detection and firing rate estimation steps until the spike indicator function does not change anymore.
3 Simulation
3.1 Metrics for quantifying accuracy
We conducted a set of simulation studies to evaluate the performance of our MTVPAR method. Among the many competing methods, we choose the approach in Jewell and Witten [2018] because its performance has been shown better than competing methods using both simulated and benchmark data. When comparing estimated spike trains with the ground truth, we use the VictorPurpura (VP) distance [Victor and Purpura, 1996, 1997], which has been commonly used for comparing spike trains. It is defined as the minimum total cost required to transform one spike train into another using the following three basic operations:

Insert a spike into a spike train. (Cost = 1)

Delete a spike. (Cost = 1)

Shift a spike by an interval . (Cost = ) A large makes the distance more sensitive to fine timing differences. We use the default value .
Because MTVPAR estimates the firing rate function together with the spikes, we also evaluate its performance on firing rate estimation. The approach in Jewell and Witten [2018] does not estimate firing rates; hence, for a fair comparison, we apply the same Gaussianboxcar kernel smoothing in (3) to the spikes estimated using Jewell and Witten [2018] in order to estimate the firing rates. The accuracy of firing rate estimation is calibrated by the norm [Adams et al., 2009] of the difference between an estimated firing rate function and the true function:
where is the weight at and its default value is a constant.
3.2 Simulation of spike trains from inhomogeneous Poisson processes
In this simulation setting, the trials are treated as repeated trials. In other words, the trials share the same underlying firing rate function and each trial is an independent realization of an inhomogeneous Poisson process. Several forms of firing rate functions have been considered in previous work. For example, Kass et al. [2003], Behseta and Kass [2005] assumed bellshaped firing rate functions. Pachitariu et al. [2018] used a piecewise constant stimulus rate. Firing rate functions with multiple peaks from exponential stimuli functions have also been considered [ReynaudBouret et al., 2014]. In our simulation, we chose the following bimodal firing rate function
(10) 
where and the stimuli peaks are at and . Thus, reaches its maximum value 0.2 (equivalent to 10 spikes per second) at time 300 and 700.
To generate spike trains for multiple trials under an inhomogeneous Poisson process with the firing rate function in 10, we followed the idea of Adams et al. [2009]. In each trial, the spikes were first randomly drawn from a homogeneous Poisson process. A thinning process was then applied to create a realization from the desired inhomogeneous Poisson process (see Appendix B for details). After obtaining the simulated spike trains, we generated calcium fluorescence traces using the autoregressive model (1). The following parameters were used in the simulations: , , and trials in total. According to Vogelstein et al. [2010], the parameters above correspond to a sampling rate of 50 Hz and a length of 20 seconds per trial. For each simulation setting, 100 data sets were generated.
We implemented our MTVPAR using Algorithm 2. When estimating , we chose the Gaussian kernel bandwidth equal to 200 ms. Because the trials are assumed to have the same underlying firing rate function, B=50 is used, which is equivalent to averaging the spike counts across trials to estimate the firing rate . As shown in the top panel of Figure 1, the firing rate functions estimated from MTVPAR are closer to the true underlying function. Thus, adopting the proposed timevarying penalty led to an improvement in both spike detection (measured by mean VP distance, the lower panel of Figure 1) and firing rate estimation (measured by the norm, the upper panel of Figure 1) across all the values considered. For the optimal based on VP, the reductions brought by the timevarying penalty in VP and norm are and , respectively.
We also simulated data with a larger noise level (), longer series () and higher decay rate (). In all scenarios, we observed improved accuracy from using the proposed MTVPAR (results omitted).
3.3 Simulation under a dynamic firing rate function
We also simulated data under dynamic firing rate functions across trials. Figure 2 (the upper left panel) shows a dynamic firing rate function with two peaks within each trial; across trials, the peak values first increase then decrease. The twodimensional firing rate function at time point and trial is as follows:
(11) 
where and .
Figure 2 summarizes the simulation results with , , and . The sampling rate is 50 Hz and the length is 20 seconds per trial. The maximum of the firing rate is 0.2 (equivalent to 10 spikes per second). The results are based on 100 simulations. In firing rate estimation, Gaussian kernel smoothing is applied with a kernel bandwidth of 200 ms within trials and a window length of 10 across trials. Similar to the simulation in Section 3.2, MTVPAR performs much better than using a constant penalty in estimating spikes and firing rate. Specifically, at the optimal (based on minimal VP distance), the reduction brought by MTVPAR in VP distance for spike detection is and the reduction of error of firing rate estimation is .
4 Application to Calcium Imaging Data
We now pursue the scientific investigations on the neural activities in two longitudinal mice studies in which both calcium imaging data and behavioral data were recorded. In this first data set, calcium imaging data were collected within a few days after the participating mice had been trained for a discrimination task; thus, it is reasonable to treat trials with same behavior outcome as repeated trials. As a comparison, in the second study, calcium recording started in the first learning trial and lasted for a few weeks; therefore, neural dynamics are expected as a result of the learning process and neural plasticity. For this reason, incorporating neuronal dynamics in firing rate estimation is likely to improve spike detection and firing rate estimation.
4.1 Mouse task data I
The activity of neurons (labelled with GCaMP6s) from the mouse’s anterior lateral motor cortex (ALM) was recorded with twophoton calcium imaging during a headfixed whiskerbased discrimination task [Li et al., 2015]
. In the experiment, mice were supposed to discriminate the pole locations using their whiskers and report the perceived pole position by licking. Each trial is composed of three epochs: sample epoch (mice presented with a vertical pole), delay epoch (the pole was removed), response epoch (mice cued to give a response). If a mouse licked the correct lick port, it was rewarded with liquid.
The data we present here is from one mouse with 73 trials from multiple days. Each trial contained more than 100 data points at 15 Hz. The fluorescence data were obtained after typical preprocessing procedures such as correction for neuropil contamination and transformation where is the baseline averaged fluorescence within a s period right before the start of each trial. There were four possible behavioral outcomes in the experiment: correct/incorrect lick left/right. Because most of the trials were either “correct lick left” or “correct lick right”, we combined the two incorrect groups as a single group. Figure 3 shows the calcium fluorescence traces of a pyramidal tract neuron in 73 trials, including 31 trials of “correct lick left”, 21 trials of “correct lick right”, and 21 trials of “incorrect lick” .
One interesting question is whether the neuron responded differently for different outcomes. Therefore, we conducted a stratified analysis for each outcome. In the firing rate estimation, we use a Gaussian kernel smoothing with bandwidth 200 ms within trials. Since the mouse has been well trained before calcium imaging recording, firing rates across trials under the same outcome group are relatively stable, which was confirmed by available 2D visualization (data not shown). Thus, it is sensible to combine all the trials within an outcome type when estimating firing rate.
As indicated in Figure 3, in “correct lick left” trials, the neuron fired right after the cue time (when the mouse was cued to make decision); however, there was almost no neural activity under the other two outcome groups. The estimated firing rate function also confirmed this difference. It is known that the ALM brain region of mice is involved in planning directed licking [Guo et al., 2014]. The estimated spikes and firing rate functions provide convincing evidences that this neuron is likely to show neural selectivity and play a critical role in the cognitive process of making the correct decision of licking the left pole.
4.2 Mouse task data II
The second data set is from our study of longterm contextual discrimination experiment, in which mice were trained to recognize two contexts via fear conditioning (foot shock) [Johnston et al., 2020]. The research goal was to understand the behaviorassociated hippocampal neural ensemble dynamics at singlecell resolution. Viral injections were administered into mice brain to introduce a genetically encoded calcium indicator (GCaMP6f). Fluorescence signals from hippocampal CA1 excitatory neurons were then optically recorded using one photon headmounted miniscopes from freely moving mice (Figure 4, upper left). There were four stages in the experiment: habituation (mice freely exploring environment), learning (learning to freeze in a stimulus context with foot shock), extinction (no foot shock), and relearning (stimulus reinstated).
In the mouse we analyzed, the activity of 141 neurons were recorded for several weeks. We chose the 21 foot shock sessions with the first 11 sessions in the learning stage and the remaining 10 sessions in the relearning stage. In each shock session, a foot shock was administered and the fluorescence trace was recorded at 15 Hz for 2 minutes. Figure 4 (upper right panel) shows an example neuron. The calcium fluorescence traces from different trials were temporally aligned by the start of shock time. For the firing rate estimation, we apply the Gaussianboxcar kernel smoothing with a bandwidth of 400 ms for within trials and a window length of 5 for across trials. The estimated spikes (Figure 4, upper right) and firing rate functions (Figure 4, lower left) suggest that the neuron is more synchronous to the stimulus in the relearning stage than in the learning stage, which may reflect the evolving learningrelated neuronal dynamics.
It is worth noting that treating the trials as independent realizations of the same underlying process will lead to undesired results. As shown in Figure 4 (lower right panel), when assuming the trials as samples from the same underlying distribution, the estimated peak firing time is misleading. Although the stratified estimates from the two stages showed that neuronal firing took place sooner and was more frequent in the relearning stage than in the learning stage in response to the foot shock stimulus, the results based on stratified analysis were not able to completely characterize the intrinsic evolution of neural firing during the cognitive learning process.
5 Discussion
In this paper, we proposed MTVPAR, a timevarying penalized method to simultaneously conduct spike detection and firing rate estimation from longitudinal calcium fluorescence trace data. In MTVPAR, each iteration consists of two steps: the spike detection step using timevarying regularization based on the current estimate of the firing rate functions and the firing rate estimation step based on the currently detected spikes. Our approach is able to account for the intrinsic neural dynamics at both within and betweentrial levels. Here we used a timevarying penalization in the spike detection step and a Gaussianboxcar smoothing in the firing rate estimation step. In some situations, other types of regularization and smoothing methods might be preferred or more appropriate; our strategy of integrating information from multiple trials can be adopted in a similar manner.
Our work was motivated by recovering the underlying spike data from noisy calcium imaging and the main focus had been the accuracy of spike detection. Borrowing information across neighboring trials improves spike estimation; as a consequence, the firing rate estimation is also improved. In our local nonparametric estimation for firing rate functions, for computational ease, we chose the window size based on experience and the total number of available trials. The window size and other tuning parameters are best calibrated using benchmark data. In the absence of benchmark data, data driven methods might be helpful. Cunningham et al. [2007]
proposed to use a Gaussian process prior and an iterative gradient optimization algorithm to optimally choose hyperparameters.
Fiecas and Ombao [2016] proposed to use the window size that minimizes the generalized crossvalidated deviance when investigating the evolution of dynamic interactions (coherence) between brain regions. Adopting similar strategies is likely to gain further insights and efficiency when analyzing calcium imaging data.Like other existing methods, we aim to estimate the time of a neural spike or a burst of spike events. Bounded by the temporal resolution of calcium imaging data, when there are multiple spikes within the same time frame, we are not able to zoom in to estimate the times of individual spikes. In several works, the change in fluorescence intensities due to spike events at time , i.e., , was initially motivated to model spike counts [Vogelstein et al., 2010, Friedrich and Paninski, 2016, Friedrich et al., 2017]; following this line of thought, it is reasonable to assume that is positively correlated with the spike counts for rescaled fluorescence traces. In other work [Pnevmatikakis et al., 2013, Jewell and Witten, 2018], the sign of was the focus. The unknown and nonlinear relationship between spike counts and fluorescence transient measured based upon genetically encoded indicators [Vogelstein et al., 2010, Lütcke et al., 2013, Rose et al., 2014] complicates the accurate interpretation and modeling of fluorescence transient data. In the current version of MTVPAR, the magnitude of was not used in the firing rate estimation. It is worth conducting future research to examine how to best utilize in analyzing calcium fluorescence data.
Appendix A Proofs for Section 2.3
a.1 Equivalence between Problem (5) and Problem (6)
a.2 Equivalence between Problem (5) and Problem (9)
Based on the result in Appendix A.1, to show that Problems (5) and (9) are equivalent, we only need to show that Problems (6) and (9) are equivalent.
Appendix B Algorithm for Simulating Spike Trains from an Inhomogeneous Poisson Process
References

Tractable nonparametric bayesian inference in poisson processes with gaussian process intensities
. InProceedings of the 26th Annual International Conference on Machine Learning
, pp. 9–16. Cited by: §3.1, §3.2, Algorithm 3.  Interpreting in vivo calcium signals from neuronal cell bodies, axons, and dendrites: a review. Neurophotonics 7 (1), pp. 011402. Cited by: §2.2.
 Algorithms for the optimal identification of segment neighborhoods. Bulletin of mathematical biology 51 (1), pp. 39–54. Cited by: §2.3.
 Testing equality of two functions using bars. Statistics in medicine 24 (22), pp. 3523–3534. Cited by: §2.2, §3.2.
 Likelihood methods for neural spike train data analysis. Computational neuroscience: A comprehensive approach, pp. 253–286. Cited by: §2.2.
 Theory of point processes for neural systems. In Les Houches, Vol. 80, pp. 691–727. Cited by: §2.2.
 Inferring neural firing rates from spike trains using gaussian processes. In Advances in neural information processing systems, pp. 329–336. Cited by: §2.2.
 Methods for estimating neural firing rates, and their application to brain–machine interfaces. Neural Networks 22 (9), pp. 1235–1246. Cited by: §2.2, §2.2.
 Inferring neural firing rates from spike trains using gaussian processes. Advances in neural information processing systems 20, pp. 329–336. Cited by: §5.
 Analysis of betweentrial and withintrial neural spiking dynamics. Journal of neurophysiology 99 (5), pp. 2672–2693. Cited by: §2.2.
 An introduction to the theory of point processes: volume ii: general theory and structure. Springer Science & Business Media. Cited by: §2.2.
 Accurate spike estimation from noisy calcium signals for ultrafast threedimensional imaging of large neuronal populations in vivo. Nature communications 7 (1), pp. 1–17. Cited by: §1.
 Bayesian curvefitting with freeknot splines. Biometrika 88 (4), pp. 1055–1071. Cited by: §2.2.

Composite quantile regression for the singleindex model
. Preprint. Cited by: §2.2. 
Modeling the evolution of dynamic brain processes during an associative learning experiment
. Journal of the American Statistical Association 111 (516), pp. 1440–1453. Cited by: §2.2, §5.  Fast active set methods for online spike inference from calcium imaging. In Advances In Neural Information Processing Systems, pp. 1984–1992. Cited by: §2.3, §5.
 Fast online deconvolution of calcium imaging data. PLoS computational biology 13 (3), pp. e1005423. Cited by: §2.3, §2.3, §5.
 An approach to the quantitative analysis of electrophysiological data from single neurons. Biophysical Journal 1 (1), pp. 15. Cited by: §2.2.
 Miniaturized integration of a fluorescence microscope. Nature methods 8 (10), pp. 871. Cited by: §1.

CaImAn an open source tool for scalable calcium imaging data analysis
. Elife 8, pp. e38173. Cited by: §1.  Flow of cortical activity underlying a tactile decision in mice. Neuron 81 (1), pp. 179–194. Cited by: §4.1.
 Beyond trialbased paradigms: continuous behavior, ongoing neural activity, and natural stimuli. Journal of Neuroscience 38 (35), pp. 7551–7558. Cited by: §2.2.
 An algorithm for optimal partitioning of data on an interval. IEEE Signal Processing Letters 12 (2), pp. 105–108. Cited by: §2.3.
 Exact spike train inference via l0 optimization. The annals of applied statistics 12 (4), pp. 2457. Cited by: §2.3, §2.3, §3.1, §3.1, §5.
 Building neural representations of habits. Science 286 (5445), pp. 1745–1749. Cited by: §2.2.
 Robust population single neuronal calcium signal extraction using scout allows for longitudinal analysis of behaviorassociated neural ensemble dynamics. bioRxiv. Cited by: §1, Figure 4, §4.2.
 Statistical issues in the analysis of neuronal data. Journal of neurophysiology 94 (1), pp. 8–25. Cited by: §2.2, §2.2.
 Statistical smoothing of neuronal data. NetworkComputation in Neural Systems 14 (1), pp. 5–16. Cited by: §2.2, §3.2.
 Optimal detection of changepoints with a linear computational cost. Journal of the American Statistical Association 107 (500), pp. 1590–1598. Cited by: §2.3.
 A motor cortex circuit for motor planning and movement. Nature 519 (7541), pp. 51–56. Cited by: Figure 3, §4.1.
 Inference of neuronal network spike dynamics and topology from calcium imaging data. Frontiers in neural circuits 7, pp. 201. Cited by: §5.
 Detecting cells using nonnegative matrix factorization on calcium imaging data. Neural Networks 55, pp. 11–19. Cited by: §1.
 Learning population and subjectspecific brain connectivity networks via mixed neighborhood selection. The Annals of Applied Statistics 11 (4), pp. 2142–2164. Cited by: §2.2.
 Automated analysis of cellular signals from largescale calcium imaging data. Neuron 63 (6), pp. 747–760. Cited by: §1.
 Neuronal activity in macaque supplementary eye field during planning of saccades in response to pattern and spatial cues. Journal of Neurophysiology 84 (3), pp. 1369–1384. Cited by: §2.2.
 Statistical models for brain signals with properties that evolve across trials. NeuroImage 180, pp. 609–618. Cited by: §2.2, §2.2.
 Robustness of spike deconvolution for neuronal calcium imaging. Journal of Neuroscience 38 (37), pp. 7976–7985. Cited by: §1, §3.2.
 A new look at statespace models for neural data. Journal of computational neuroscience 29 (12), pp. 107–126. Cited by: §2.2.
 Populationlevel representation of a temporal sequence underlying song production in the zebra finch. Neuron 90 (4), pp. 866–876. Cited by: §1.
 Bayesian spike inference from calcium imaging data. In 2013 Asilomar Conference on Signals, Systems and Computers, pp. 349–353. Cited by: §1, §5.
 Simultaneous denoising, deconvolution, and demixing of calcium imaging data. Neuron 89 (2), pp. 285–299. Cited by: §1, §2.3.
 Analysis pipelines for calcium imaging data. Current opinion in neurobiology 55, pp. 15–21. Cited by: §1.
 Goodnessoffit tests and nonparametric adaptive estimation for spike train analysis. The Journal of Mathematical Neuroscience 4 (1), pp. 3. Cited by: §3.2.
 Putting a finishing touch on gecis. Frontiers in molecular neuroscience 7, pp. 88. Cited by: §5.
 Reliabilitydriven, spatiallyadaptive regularization for deformable registration. In International Workshop on Biomedical Image Registration, pp. 173–185. Cited by: §2.2.
 Benchmarking spike rate inference in population calcium imaging. Neuron 90 (3), pp. 471–482. Cited by: §1.
 Imaging neural activity in worms, flies and mice with improved gcamp calcium indicators. Nature methods 6 (12), pp. 875–881. Cited by: §1.
 Nature and precision of temporal coding in visual cortex: a metricspace analysis. Journal of neurophysiology 76 (2), pp. 1310–1326. Cited by: §3.1.
 Metricspace analysis of spike trains: theory, algorithms and application. Network: computation in neural systems 8 (2), pp. 127–164. Cited by: §3.1.
 Fast nonnegative deconvolution for spike train inference from population calcium imaging. Journal of neurophysiology 104 (6), pp. 3691–3704. Cited by: §1, §2.1, §2.3, §2.3, §3.2, §5.
 Single neurons in the monkey hippocampus and learning of new associations. Science 300 (5625), pp. 1578–1581. Cited by: §2.2.
 Role of the hippocampal system in conditional motor learning: mapping antecedents to action. Hippocampus 9 (2), pp. 101–117. Cited by: §2.2.
 Reconstruction of firing rate changes across neuronal populations by temporally deconvolved ca 2+ imaging. Nature methods 3 (5), pp. 377–383. Cited by: §1, §2.3.
 FRM: a financial risk meter based on penalizing tail events occurrence. Preprint. Cited by: §2.2.
 Time varying quantile lasso. Preprint. Cited by: §2.2.
 Efficient and accurate extraction of in vivo calcium signals from microendoscopic video data. Elife 7, pp. e28728. Cited by: §1.