Log In Sign Up

Time-varying ℓ_0 optimization for Spike Inference from Multi-Trial Calcium Recordings

Optical imaging of genetically encoded calcium indicators is a powerful tool to record the activity of a large number of neurons simultaneously over a long period of time from freely behaving animals. However, determining the exact time at which a neuron spikes and estimating the underlying firing rate from calcium fluorescence data remains challenging, especially for calcium imaging data obtained from a longitudinal study. We propose a multi-trial time-varying ℓ_0 penalized method to jointly detect spikes and estimate firing rates by robustly integrating evolving neural dynamics across trials. Our simulation study shows that the proposed method performs well in both spike detection and firing rate estimation. We demonstrate the usefulness of our method on calcium fluorescence trace data from two studies, with the first study showing differential firing rate functions between two behaviors and the second study showing evolving firing rate function across trials due to learning.


Estimating a Separably-Markov Random Field (SMuRF) from Binary Observations

A fundamental problem in neuroscience is to characterize the dynamics of...

Scalable Spike Source Localization in Extracellular Recordings using Amortized Variational Inference

Extracellular recordings using modern, dense probes provide detailed foo...

Analyzing second order stochasticity of neural spiking under stimuli-bundle exposure

Conventional analysis of neuroscience data involves computing average ne...

Fast Nonconvex Deconvolution of Calcium Imaging Data

Calcium imaging data promises to transform the field of neuroscience by ...

Cross-population coupling of neural activity based on Gaussian process current source densities

Because LFPs arise from multiple sources in different spatial locations,...

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 single-neuron 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 multi-trial methods we identified, Picardo et al. [2016] assumed that repeated trials share the same burst time but have trial-specific magnitude, baseline fluorescence, and noise level. Conceptually, integrating multiple trials should be beneficial if the trial-to-trial 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 multi-trial 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 time-varying firing rate function into a dynamic penalization framework, rather than forcing shared parameters across trials.

2 Methods

In this section, we first propose a time-varying penalized framework to analyze multi-trial 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 Multi-Trial time-varying penalized auto-regressive model (MTV-PAR)

Let be the fluorescence recorded at time point of trial , where , . In practice, multiple pre-processing steps, including normalization, [Vogelstein et al., 2010], are typically implemented to obtain . The calcium fluorescence trace is often modeled using the following first-order auto-regressive model


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 time-varying penalty function on the number of spikes, as the spikes tend to be sparse but not evenly distributed over time:


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 trial-specific 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 time-varying penalization.

2.2 Firing rate estimation and time-varying 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 between-trial dynamics. To account for the between-trial dynamics, some authors have proposed a state-space framework [Czanner et al., 2008, Paninski et al., 2010]. In the state-space framework, the spike train is characterized by a point process model [Brown et al., 2003, Brown, 2005, Kass et al., 2005, Daley and Vere-Jones, 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 multi-trial fluorescence data, we consider instead a computationally less demanding two-dimensional Gaussian-boxcar kernel smoothing function , which is formulated as follows


where denotes the within-trial kernel bandwidth in the Gaussian kernel, is the indicator function, and is a bandwidth for the between-trial 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 50-150 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 time-bin of 200 ms is suggested in one of their analyses. In our analysis, is chosen between 200ms and 400ms. For the box-car bandwidth to allow borrowing information across trials, we use a pre-selected 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 non-constant penalization is inspired by prior work in the literature. For example, time-varying 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 time-varying 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 time-varying 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 time-varying 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 ,


2.3 Time-varying penalized AR(1) model

In our proposed multi-trial time-varying penalized auto-regressive model (MTV-PAR), a calcium fluorescence trace at the th trial , is modeled by a first-order auto-regressive 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 time-varying penalization problem can also be solved by a dynamic programming algorithm:


Specifically, we find that the time-varying 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):


where are change points, i.e., the points satisfying and


which has the following closed-form solution:


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 auto-correlation at lag 1.

Finally, as shown in Appendix A.2, the optimization problem (5) can be solved by computing recursively:


The resulting algorithm has a time complexity . This optimization problem can be solved in time (Algorithm 1) with a dynamic programming algorithm [Auger and Lawrence, 1989, Jackson et al., 2005, Jewell and Witten, 2018, Killick et al., 2012] .

1:Input: Time-varying penalty function , single-trial fluorescence with time points.
3:for  do
4:     Calculate
5:     Find
6:     Let
7:     if   then
8:         Add to
9:     end if
10:end for
Algorithm 1 Dynamic programming algorithm to detect spikes with a time-varying penalization function

2.4 Algorithms

Thus far, we have presented our solutions to two problems separately. The first problem is to use multi-trial spike data to estimate the firing rate function for and using a Gaussian-boxcar smoothing method. The second problem is to detect spikes from a single calcium fluorescence trace using a time-varying penalized approach under the assumption that the time-varying penalty function

is already known. In practice, neither firing rate functions nor spike locations are known. We therefore propose an Expectation-Maximization-type 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.

1:Input: Multi-trial fluorescence data with trials and time points. Penalty term .
2:Initialize: Penalty function .
3:Set the binary spike indicator .
4:while  the indicator matrix changes do
5:     for  do
6:         Apply Algorithm 1 to detect spikes and let denote the times at which spikes were detected.
7:         Set
8:     end for
9:     Estimate the firing rate , for and using Gaussian-Sliding-Window kernel smoothing.
10:     Calculate the weight function
11:     Update the time-varying penalty function
12:end while
Algorithm 2 Simultaneous spike detection and firing rate estimation

3 Simulation

3.1 Metrics for quantifying accuracy

We conducted a set of simulation studies to evaluate the performance of our MTV-PAR 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 Victor-Purpura (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 MTV-PAR 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 Gaussian-boxcar 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 bell-shaped 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 [Reynaud-Bouret et al., 2014]. In our simulation, we chose the following bi-modal firing rate function


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 auto-regressive 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 MTV-PAR 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 MTV-PAR are closer to the true underlying function. Thus, adopting the proposed time-varying 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 time-varying 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 MTV-PAR (results omitted).

Figure 1: Top: firing rate estimates using MTV-PAR (time-varying) and constant penalty. Bottom: VP distance of spike trains and norm of firing rates between the truth and estimators.

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 two-dimensional firing rate function at time point and trial is as follows:


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, MTV-PAR 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 MTV-PAR in VP distance for spike detection is and the reduction of error of firing rate estimation is .

Figure 2: Top Left: True firing rates; Top Middle: Estimated firing rates using MTV-PAR; Top Right: Estimated firing rates using constant penalty. Bottom: VP distance of spike trains and L2 norm of firing rates between the truth and estimators.

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 two-photon calcium imaging during a head-fixed whisker-based 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 pre-processing 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.

Figure 3: An example cell from the water lick experiment [Li et al., 2015]. Top Left: The calcium fluorescence trace under correctly lick left outcome (black dots: estimated spikes). Top Right: The calcium fluorescence trace under correctly lick right outcome. Bottom Left: The calcium fluorescence trace under incorrect lick outcomes. Bottom Right: Estimated firing rate functions of under different conditions. Vertical dashed line: the start of response epoch.

4.2 Mouse task data II

The second data set is from our study of long-term 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 behavior-associated hippocampal neural ensemble dynamics at single-cell 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 head-mounted 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 Gaussian-boxcar 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 learning-related 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.

Figure 4: An example neuron from the foot shock study. Top Left: Calcium imaging of mouse hippocampus using a miniaturized scope and the four stages of the experimental design: habituation, learning, extinction, and relearning [Johnston et al., 2020]. Top Right: Calcium fluorescence traces and estimated spikes of a sample neuron. The red and blue traces were obtained during the learning and relearning stages, respectively. The long dashed vertical line denotes the time when a foot shock was applied and the short black vertical lines are the time locations at which spikes were detected. Bottom Left: the estimated firing rate function using MTV-PAR. Bottom Right: the estimated firing rate function using stratified analysis. Vertical dashed line: the start time of foot shock.

5 Discussion

In this paper, we proposed MTV-PAR, a time-varying penalized method to simultaneously conduct spike detection and firing rate estimation from longitudinal calcium fluorescence trace data. In MTV-PAR, each iteration consists of two steps: the spike detection step using time-varying 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 between-trial levels. Here we used a time-varying penalization in the spike detection step and a Gaussian-boxcar 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 cross-validated 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 re-scaled 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 MTV-PAR, 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

1:Input Length of data , upper bound of rate intensity , firing rate function
2:Draw the number of spikes:
3:Uniformly distribute the spikes:
4:for  do
5:     if  then
6:         exclude from
7:     end if
8:end for
9:Return the set of spikes
Algorithm 3 Simulate Spike Trains from an Inhomogeneous Poisson Process ( Modified from Algorithm 1 of Adams et al. [2009])


  • R. P. Adams, I. Murray, and D. J. MacKay (2009)

    Tractable nonparametric bayesian inference in poisson processes with gaussian process intensities


    Proceedings of the 26th Annual International Conference on Machine Learning

    pp. 9–16. Cited by: §3.1, §3.2, Algorithm 3.
  • F. Ali and A. C. Kwan (2019) Interpreting in vivo calcium signals from neuronal cell bodies, axons, and dendrites: a review. Neurophotonics 7 (1), pp. 011402. Cited by: §2.2.
  • I. E. Auger and C. E. Lawrence (1989) Algorithms for the optimal identification of segment neighborhoods. Bulletin of mathematical biology 51 (1), pp. 39–54. Cited by: §2.3.
  • S. Behseta and R. E. Kass (2005) Testing equality of two functions using bars. Statistics in medicine 24 (22), pp. 3523–3534. Cited by: §2.2, §3.2.
  • E. N. Brown, R. Barbieri, U. T. Eden, and L. M. Frank (2003) Likelihood methods for neural spike train data analysis. Computational neuroscience: A comprehensive approach, pp. 253–286. Cited by: §2.2.
  • E. N. Brown (2005) Theory of point processes for neural systems. In Les Houches, Vol. 80, pp. 691–727. Cited by: §2.2.
  • J. P. Cunningham, M. Y. Byron, K. V. Shenoy, and M. Sahani (2008) Inferring neural firing rates from spike trains using gaussian processes. In Advances in neural information processing systems, pp. 329–336. Cited by: §2.2.
  • J. P. Cunningham, V. Gilja, S. I. Ryu, and K. V. Shenoy (2009) 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.
  • J. P. Cunningham, B. M. Yu, K. V. Shenoy, and M. Sahani (2007) Inferring neural firing rates from spike trains using gaussian processes. Advances in neural information processing systems 20, pp. 329–336. Cited by: §5.
  • G. Czanner, U. T. Eden, S. Wirth, M. Yanike, W. A. Suzuki, and E. N. Brown (2008) Analysis of between-trial and within-trial neural spiking dynamics. Journal of neurophysiology 99 (5), pp. 2672–2693. Cited by: §2.2.
  • D. J. Daley and D. Vere-Jones (2007) An introduction to the theory of point processes: volume ii: general theory and structure. Springer Science & Business Media. Cited by: §2.2.
  • T. Deneux, A. Kaszas, G. Szalay, G. Katona, T. Lakner, A. Grinvald, B. Rózsa, and I. Vanzetta (2016) Accurate spike estimation from noisy calcium signals for ultrafast three-dimensional imaging of large neuronal populations in vivo. Nature communications 7 (1), pp. 1–17. Cited by: §1.
  • I. DiMatteo, C. R. Genovese, and R. E. Kass (2001) Bayesian curve-fitting with free-knot splines. Biometrika 88 (4), pp. 1055–1071. Cited by: §2.2.
  • Y. Fan, W. K. Härdle, W. Wang, and L. Zhu (2013)

    Composite quantile regression for the single-index model

    Preprint. Cited by: §2.2.
  • M. Fiecas and H. Ombao (2016)

    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.
  • J. Friedrich and L. Paninski (2016) 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.
  • J. Friedrich, P. Zhou, and L. Paninski (2017) Fast online deconvolution of calcium imaging data. PLoS computational biology 13 (3), pp. e1005423. Cited by: §2.3, §2.3, §5.
  • G. L. Gerstein and N. Y. Kiang (1960) An approach to the quantitative analysis of electrophysiological data from single neurons. Biophysical Journal 1 (1), pp. 15. Cited by: §2.2.
  • K. K. Ghosh, L. D. Burns, E. D. Cocker, A. Nimmerjahn, Y. Ziv, A. El Gamal, and M. J. Schnitzer (2011) Miniaturized integration of a fluorescence microscope. Nature methods 8 (10), pp. 871. Cited by: §1.
  • A. Giovannucci, J. Friedrich, P. Gunn, J. Kalfon, B. L. Brown, S. A. Koay, J. Taxidis, F. Najafi, J. L. Gauthier, P. Zhou, et al. (2019)

    CaImAn an open source tool for scalable calcium imaging data analysis

    Elife 8, pp. e38173. Cited by: §1.
  • Z. V. Guo, N. Li, D. Huber, E. Ophir, D. Gutnisky, J. T. Ting, G. Feng, and K. Svoboda (2014) Flow of cortical activity underlying a tactile decision in mice. Neuron 81 (1), pp. 179–194. Cited by: §4.1.
  • A. Huk, K. Bonnen, and B. J. He (2018) Beyond trial-based paradigms: continuous behavior, ongoing neural activity, and natural stimuli. Journal of Neuroscience 38 (35), pp. 7551–7558. Cited by: §2.2.
  • B. Jackson, J. D. Scargle, D. Barnes, S. Arabhi, A. Alt, P. Gioumousis, E. Gwin, P. Sangtrakulcharoen, L. Tan, and T. T. Tsai (2005) An algorithm for optimal partitioning of data on an interval. IEEE Signal Processing Letters 12 (2), pp. 105–108. Cited by: §2.3.
  • S. Jewell and D. Witten (2018) 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.
  • M. S. Jog, Y. Kubota, C. I. Connolly, V. Hillegaart, and A. M. Graybiel (1999) Building neural representations of habits. Science 286 (5445), pp. 1745–1749. Cited by: §2.2.
  • K. G. Johnston, S. F. Grieco, Z. Yu, S. Jin, T. Shen, R. Crary, J. Guzowski, T. Holmes, Q. C. Nie, and X. Xu (2020) Robust population single neuronal calcium signal extraction using scout allows for longitudinal analysis of behavior-associated neural ensemble dynamics. bioRxiv. Cited by: §1, Figure 4, §4.2.
  • R. E. Kass, V. Ventura, and E. N. Brown (2005) Statistical issues in the analysis of neuronal data. Journal of neurophysiology 94 (1), pp. 8–25. Cited by: §2.2, §2.2.
  • R. E. Kass, V. Ventura, and C. Cai (2003) Statistical smoothing of neuronal data. Network-Computation in Neural Systems 14 (1), pp. 5–16. Cited by: §2.2, §3.2.
  • R. Killick, P. Fearnhead, and I. A. Eckley (2012) Optimal detection of changepoints with a linear computational cost. Journal of the American Statistical Association 107 (500), pp. 1590–1598. Cited by: §2.3.
  • N. Li, T. Chen, Z. V. Guo, C. R. Gerfen, and K. Svoboda (2015) A motor cortex circuit for motor planning and movement. Nature 519 (7541), pp. 51–56. Cited by: Figure 3, §4.1.
  • H. Lütcke, F. Gerhard, F. Zenke, W. Gerstner, and F. Helmchen (2013) Inference of neuronal network spike dynamics and topology from calcium imaging data. Frontiers in neural circuits 7, pp. 201. Cited by: §5.
  • R. Maruyama, K. Maeda, H. Moroda, I. Kato, M. Inoue, H. Miyakawa, and T. Aonishi (2014) Detecting cells using non-negative matrix factorization on calcium imaging data. Neural Networks 55, pp. 11–19. Cited by: §1.
  • R. P. Monti, C. Anagnostopoulos, G. Montana, et al. (2017) Learning population and subject-specific brain connectivity networks via mixed neighborhood selection. The Annals of Applied Statistics 11 (4), pp. 2142–2164. Cited by: §2.2.
  • E. A. Mukamel, A. Nimmerjahn, and M. J. Schnitzer (2009) Automated analysis of cellular signals from large-scale calcium imaging data. Neuron 63 (6), pp. 747–760. Cited by: §1.
  • C. R. Olson, S. N. Gettner, V. Ventura, R. Carta, and R. E. Kass (2000) 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.
  • H. Ombao, M. Fiecas, C. Ting, and Y. F. Low (2018) Statistical models for brain signals with properties that evolve across trials. NeuroImage 180, pp. 609–618. Cited by: §2.2, §2.2.
  • M. Pachitariu, C. Stringer, and K. D. Harris (2018) Robustness of spike deconvolution for neuronal calcium imaging. Journal of Neuroscience 38 (37), pp. 7976–7985. Cited by: §1, §3.2.
  • L. Paninski, Y. Ahmadian, D. G. Ferreira, S. Koyama, K. R. Rad, M. Vidne, J. Vogelstein, and W. Wu (2010) A new look at state-space models for neural data. Journal of computational neuroscience 29 (1-2), pp. 107–126. Cited by: §2.2.
  • M. A. Picardo, J. Merel, K. A. Katlowitz, D. Vallentin, D. E. Okobi, S. E. Benezra, R. C. Clary, E. A. Pnevmatikakis, L. Paninski, and M. A. Long (2016) Population-level representation of a temporal sequence underlying song production in the zebra finch. Neuron 90 (4), pp. 866–876. Cited by: §1.
  • E. A. Pnevmatikakis, J. Merel, A. Pakman, and L. Paninski (2013) Bayesian spike inference from calcium imaging data. In 2013 Asilomar Conference on Signals, Systems and Computers, pp. 349–353. Cited by: §1, §5.
  • E. A. Pnevmatikakis, D. Soudry, Y. Gao, T. A. Machado, J. Merel, D. Pfau, T. Reardon, Y. Mu, C. Lacefield, W. Yang, et al. (2016) Simultaneous denoising, deconvolution, and demixing of calcium imaging data. Neuron 89 (2), pp. 285–299. Cited by: §1, §2.3.
  • E. A. Pnevmatikakis (2019) Analysis pipelines for calcium imaging data. Current opinion in neurobiology 55, pp. 15–21. Cited by: §1.
  • P. Reynaud-Bouret, V. Rivoirard, F. Grammont, and C. Tuleau-Malot (2014) Goodness-of-fit tests and nonparametric adaptive estimation for spike train analysis. The Journal of Mathematical Neuroscience 4 (1), pp. 3. Cited by: §3.2.
  • T. Rose, P. M. Goltstein, R. Portugues, and O. Griesbeck (2014) Putting a finishing touch on gecis. Frontiers in molecular neuroscience 7, pp. 88. Cited by: §5.
  • L. Tang, G. Hamarneh, and R. Abugharbieh (2010) Reliability-driven, spatially-adaptive regularization for deformable registration. In International Workshop on Biomedical Image Registration, pp. 173–185. Cited by: §2.2.
  • L. Theis, P. Berens, E. Froudarakis, J. Reimer, M. R. Rosón, T. Baden, T. Euler, A. S. Tolias, and M. Bethge (2016) Benchmarking spike rate inference in population calcium imaging. Neuron 90 (3), pp. 471–482. Cited by: §1.
  • L. Tian, S. A. Hires, T. Mao, D. Huber, M. E. Chiappe, S. H. Chalasani, L. Petreanu, J. Akerboom, S. A. McKinney, E. R. Schreiter, et al. (2009) Imaging neural activity in worms, flies and mice with improved gcamp calcium indicators. Nature methods 6 (12), pp. 875–881. Cited by: §1.
  • J. D. Victor and K. P. Purpura (1996) Nature and precision of temporal coding in visual cortex: a metric-space analysis. Journal of neurophysiology 76 (2), pp. 1310–1326. Cited by: §3.1.
  • J. D. Victor and K. P. Purpura (1997) Metric-space analysis of spike trains: theory, algorithms and application. Network: computation in neural systems 8 (2), pp. 127–164. Cited by: §3.1.
  • J. T. Vogelstein, A. M. Packer, T. A. Machado, T. Sippy, B. Babadi, R. Yuste, and L. Paninski (2010) 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.
  • S. Wirth, M. Yanike, L. M. Frank, A. C. Smith, E. N. Brown, and W. A. Suzuki (2003) Single neurons in the monkey hippocampus and learning of new associations. Science 300 (5625), pp. 1578–1581. Cited by: §2.2.
  • S. P. Wise and E. A. Murray (1999) Role of the hippocampal system in conditional motor learning: mapping antecedents to action. Hippocampus 9 (2), pp. 101–117. Cited by: §2.2.
  • E. Yaksi and R. W. Friedrich (2006) 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.
  • L. Yu, W. K. Härdle, L. Borke, and T. Benschop (2017) FRM: a financial risk meter based on penalizing tail events occurrence. Preprint. Cited by: §2.2.
  • L. Zbonakova, W. K. Härdle, and W. Wang (2016) Time varying quantile lasso. Preprint. Cited by: §2.2.
  • P. Zhou, S. L. Resendez, J. Rodriguez-Romaguera, J. C. Jimenez, S. Q. Neufeld, A. Giovannucci, J. Friedrich, E. A. Pnevmatikakis, G. D. Stuber, R. Hen, et al. (2018) Efficient and accurate extraction of in vivo calcium signals from microendoscopic video data. Elife 7, pp. e28728. Cited by: §1.