1 Introduction

Fluorescence indicators of calcium activity allow us to monitor the dynamics of neuronal populations both in vivo and in vitro. In the last decade there has been a proliferation of new methods to identify single spikes from fluorescence time series using template matching (Dyer et al., 2013; Grewe et al., 2010; Kerr et al., 2005; Quan et al., 2010), linear deconvolution (Holekamp et al., 2008; Tubiana et al., 2020; Wei et al., 2020; Yaksi and Friedrich, 2006), finite rate of innovation (Oñativia and Dragotti, 2015; Oñativia et al., 2013), independent component analysis (Mukamel et al., 2009), non model-based signal processing (Sebastian et al., 2019), supervised learning (Hoang et al., 2020; Rupprecht et al., 2021; Sasaki et al., 2008; Sebastian et al., 2021; Theis et al., 2016), constrained non-negative matrix factorization (Pnevmatikakis et al., 2014, 2016; Zhou et al., 2018), active set methods (Friedrich and Paninski, 2016; Friedrich et al., 2017), convex and non-convex optimization methods (Jewell and Witten, 2018; Jewell et al., 2020; Malik et al., 2011; Ranganathan and Koester, 2010; Stern et al., 2020), and interior point method (Vogelstein et al., 2010).

Model-based approaches aim to describe the mechanistic link between spikes and fluorescence. This is typically achieved by introducing a set of time-dependent “state” variables such as calcium level and baseline modulation, and their temporal dynamics. The temporal evolution of the state variables usually depends on additional parameters that remain constant over time and define the dynamics. Depending on the model, state variables might not be fully observable, in which case we need to infer them from the data. In statistics, this class of models is referred to as state-space models, and they are used in time series analysis to describe the probabilistic dependence between data and unobserved variables. Previous works have used tractable state-space models to derive maximum-a-posteriori estimates of spike times (Deneux et al., 2016; Fletcher and Rangan, 2014; Kazemipour et al., 2017; Tsunoda et al., 2010b).

The vast majority of spike inference methods provide single estimates of spike times by minimizing a cost function defined by the underlying model and constraints. This optimization approach does not provide information about the statistical uncertainty associated with single-point estimates. To address this issue, previous studies proposed Bayesian inference methods (Greenberg et al., 2018; Huys and Paninski, 2009; Im et al., 2019; Mishchencko et al., 2011; Mishchenko and Paninski, 2011; Pnevmatikakis et al., 2013; Rahmati et al., 2016; Shibue and Komaki, 2020; Speiser et al., 2017; Theis et al., 2013; Tsunoda et al., 2010a; Vogelstein et al., 2009) which give access to the full probability distribution of the unknowns given the data.

However, the state-space models used in these approaches do not take into account the possibility of burst firing and slow changes in the fluorescence baseline, which is known to be an important issue for the analysis of in-vivo recordings. Moreover, current Bayesian methods do not treat time-independent model parameters (e.g. rate constants) and dynamic variables equally. Instead, they require additional optimization procedures to calibrate model parameters, thus neglecting how the uncertainty about model parameters propagates into spike times.

Inference on non-linear and non-Gaussian state-space models is analytically intractable, requiring the application of Monte Carlo methods to obtain unbiased approximations of the posterior distributions. Because the number of unknowns in these models is of the order of the number of time steps in the fluorescence time-series, the analysis of long time series requires efficient strategies to sample from high-dimensional spaces. A major breakthrough in analyzing state-space models has been the introduction of sequential Monte Carlo methods (Chopin and Papaspiliopoulos, 2020). These algorithms can sample efficiently from the latent space by approximating sequentially the target posterior distribution by combining importance sampling and resampling techniques. In particular, the particle Gibbs algorithm can be used to obtain unbiased estimates of the joint distribution of time-independent model parameters and dynamical variables, but it has never been applied in the context of spike inference.

In this work we employ the particle Gibbs (PG) sampler on a bursting autoregressive (BAR) model of time series calcium-dependent fluorescence. Our generative model accounts for periods of high firing rates between periods of baseline (lower) firing rates. By quantifying the performance of our method (PGBAR) on the CASCADE benchmark dataset (Rupprecht et al., 2021) we showed that our approach is competitive with existing techniques. Finally, we tested PGBAR on in-vitro recordings of cerebellar granule cells using the ultrafast GCaMP8f calcium indicator, showing that our method permits the detection of spikes reliably even for high firing rates (∼ 200Hz).

2 Results

2.1 The model

We used a generative model of time-dependent fluorescence changes that will be used for spike inference. We model the normalized fluorescence F_1:T as the sum of a calcium-dependent fluorescence level c_t (hereinafter referred to as calcium level for brevity), a time-varying baseline b_t and Gaussian noise η_t

where the fluorescence noise η_t is normally distributed with zero mean and variance σ². To describe bursts of calcium transients, we introduce firing states q_t = 0, 1, associated respectively with low and high firing rates, r₀ and r₁. We allow for stochastic transitions between these two states with rates w_0→1 and w_1→0. The probability of switching from q to q^l within a sampling period ∆ is given by the transition matrix

The number of spikes at time t, s_t, is modeled by a Poisson distribution with rate r₁ when q_t = 1 otherwise with baseline firing rate r₀ when q_t = 0. The dynamics of the calcium level in response to a spike train is modeled as a second order autoregressive process

where c₀ is the initial calcium level and A controls the calcium increase upon single action potential. Note that in Eq. (3) the calcium level at time t depends on the previous calcium levels up to c_t−2. The dynamics of c_t in response to a single spike (kernel response) is characterized by a finite rise time (time to peak response) and exponential decay (see Section Response kernel for a derivation).

We can recast our model as a first order Markov process, where the state at time t only depends on the state at time t − 1. This can be done by introducing a calcium vector and a spike count vector (see S2.2.3 in Ref. Pnevmatikakis et al. (2016))

where in particular the calcium vector at time t is constructed by combining the calcium levels at current and previous time. With this definition, the calcium vector C_t satisfies the first-order Markov dynamics

In this model the calcium dynamics is a deterministic process given the sequence of spikes s_1:T.

We reparameterize the autoregressive model in terms of phenomenological parameters: peak response (A^(max)), rise time (time to peak response τ_r) and decay time (τ_d) of unitary fluorescence response. Thanks to previous characterization of GCaMP probes (Chen et al., 2013), we can more easily design prior distributions for such phenomenological parameters rather than the autoregressive model parameters γ’s, A and γ_1,2. A^(max), τ_r and τ_d, referred to as kernel parameters in the upcoming sections, can be derived from γ_1,2 and A as (see Section Reparameterization for a derivation)

Finally, the fluorescence baseline B_t is described by a Gaussian random walk with normally distributed initial condition

where ∆ is the sampling period of the time series.

In the language of state-space models, the latent state of our model is the combination of the bursting state q_t, the spike count s_t, the calcium vector C_t and the baseline b_t, whereas the fluorescence F_t, defined in Eq. (1) is our observation. The time-independent parameters of our model are the firing rate constants r_0,1, the transition rates of the 2-state bursting process W_0→1, W_1→0, the kernel parameters of the calcium indicators (peak amplitude, rise and decay constants), the initial calcium level c₀ and the fluorescence noise σ. To simplify the notation we will denote the latent space as X = {q_t, s_t, C_t, b_t} and the combination of time-independent parameters as θ.

2.2 State-space model formulation

The joint probability of the latent state trajectory X_1:T and the fluorescence observations F_1:T conditional to the time-independent parameters θ can be expressed as

where is the transition probability of the latent state, is the probability of the observed fluorescence conditional to the latent state at time t and μ^θ(X₁) is the probability distribution of the initial latent transitions and the Poisson probability of spike counts, namely

By assuming the fluorescence noise to be normally distributed, we have

To infer latent states and time-independent parameters from fluorescence observations, we need to compute the posterior probability

where P (θ) denotes the prior probability on model parameters and P (F_1:T) is the normalization factor of the posterior distribution, also known as marginal likelihood. This distribution encodes all the information about the statistics of the latent state trajectory and the model parameters. We can use it to compute point estimates but also to quantify uncertainties. The posterior distribution for general state-space models is not analytically tractable. However, we can use Monte Carlo methods to generate random samples from the posterior distribution and use them to obtain unbiased approximations of any posterior average.

2.3 Sequential Monte Carlo

In this work, we use Monte Carlo methods to approximate the posterior distribution in Eq. (13). There are two critical issues when applying these methods to state-space models. First the high-dimensionality of the posterior distribution and, second, the joint inference of state variables and constant parameters. In the general setting of a state-space model we need to estimate state variables at each time point. Consequently, the support of the target posterior distribution is a high-dimensional space. The high dimensionality of state-space models is an issue when applying standard sampling methods such as Markov Chain Monte Carlo, as their performance rapidly decreases with increasing dimensionality. Sequential Monte Carlo (SMC) methods address this issue by providing efficient sampling strategies from the latent space. The typical approach is to construct a sequential approximation of the posterior distribution where observations are accounted for iteratively. SMCs solve the problem of estimating dynamical variables at fixed model parameters (filtering).

The second issue, related to the joint estimation of parameters and dynamic variables, was addressed by Andrieu et al. (2010) with the introduction of the particle Gibbs algorithm. This algorithm alternates the sampling of model parameters and state variables as in the Gibbs sampler, with the difference that state variables are sampled from an SMC-based transition kernel that leaves the filtering distribution P (X_1:T | F_1:T, θ) invariant. In this work, we employ a version of the particle Gibbs algorithm developed by Lindsten et al. (2014), named particle Gibbs with ancestor sampling (PGAS), with better mixing properties (see Algorithm 2 and the Methods section).

To carry out inference of time-independent and dynamical variables for our model, we employed Algorithm 1, which alternates the two steps mentioned above: (1) run the PGAS transition kernel to sample a new trajectory of the state variables, and (2) sample model parameters from their full conditional distributions when analytically available, otherwise from a Metropolis-Hastings kernel.

Algorithm 1

Gibbs sampler

2.4 Validation and performance of PGBAR on simulated time-series fluorescence data

To test the performance of our inference method, we generated latent state variables and fluorescence time series from our model and compared the spikes inferred using our sampling algorithm against the ground truth simulations. In Fig. 2A we show a fluorescence time series simulated from our model. The firing pattern displays periods of increased firing rate separated by quiet time windows. By using this trace as input to Algorithm 1, we can generate a latent state trajectory and a set of model parameters at each iteration. In Fig. 2B we show 1000 samples of spike counts obtained by fitting the normalized fluorescence in Fig. 2A. The average spike counts over the random samples at each time frame (Fig. 2C) can be interpreted as the instantaneous firing rate multiplied by the sampling period. To illustrate the accuracy of our method, we calculated the spike counts within 1s time intervals for each random sample, providing the posterior distribution of the number of spikes in each time bin. As shown in Fig. 2D, the ground-truth spike counts are well within the range of the posterior. Our method allows us to infer not only spike times but also the time windows of high and low firing states (q_t = 0, 1 in the model) of the neuron. In particular, the probability of a burst-firing state (q = 1) can be obtained by averaging the firing state across the Monte Carlo samples. As shown in Fig. 2E, this probability is close to one during the ground-truth bursting periods and zero otherwise, with some degree of uncertainty at the onset and offset of the bursting period. Fig. 2F compares the ground-truth baseline and the sample average. One of the key advantages of our sampling algorithm is the joint estimation of latent states and time-independent model parameters. Figure 2G illustrates the posterior distributions associated with the fluorescence probe’s peak amplitude, noise level, decay, and rise time. The ground truth parameters used to simulate the testing time series are always close to the peak of the corresponding posterior distributions, showing that our model is identifiable.

Figure 1.

Figure 2.

To quantify the importance of having two firing states on the accuracy of the inference, we compared the performance of our method against a variant with only one global Poisson firing rate. We simulated fluorescence traces at different signal-to-noise (SNR) levels (1.1, 2 and 10), defined as the ratio between peak response A^(max) and the fluorescence noise parameter σ, and the burst firing frequency parameter, r₁ (5Hz, 10Hz, 20Hz, 50Hz). Then, we used our algorithm and its non-bursting variant to infer spikes from the fluorescence trace. To quantify the inference performance, we calculated the correlation between ground truth and estimated spikes (downsampled at 7.5 Hz for consistency with other analyses in this work), the average absolute error and the bias (average error) per time point (see Methods). The top panel of Figure 3A shows two example traces with stimulation times generated by a 5 and 50Hz Poisson distribution of spike times. For the 5Hz firing trace, the bursting model did not improve the inference accuracy compared to the variant with a single firing rate(Fig. 3A, middle and bottom). The two analyses produced comparable correlations, errors and biases (Fig. 3C-E). In contrast, the single firing rate model induced a systematic underestimate of the number of spikes for the 50 Hz trace, the original model captures the ground-truth spikes reliably. The lower performance of the non-bursting version of the model is due to the bias induced by forcing a single firing rate across the time series. While for all conditions of noise and frequency, inference using the bursting model gives unbiased spike counts (Fig. 3E), in the case of the non-bursting model, the single Poisson firing rate leads necessarily to an underestimation of the spike count during bursting time windows and an overestimation during low activity windows. At increasing noise levels and firing frequency, the performance difference between bursting and non-bursting versions of our algorithm become more pronounced (Fig. 3D), with a clear advantage of our original bursting model in increasing correlation with ground truth and reducing error.

Figure 3.

2.5 Validation of PGBAR on the CASCADE benchmark data and comparison to previous methods

To test our method on experimental data we analyzed neuronal recordings from the CASCADE benchmark dataset (Rupprecht et al., 2021), which allowed us to quantify the performance of our algorithm on different calcium indicators. We compared our method to CASCADE (Rupprecht et al., 2021), MLSpike (Deneux et al., 2016), Peeling (Lütcke et al., 2013), CaImAn (Giovannucci et al., 2019), Suite2p (Pachitariu et al., 2017) and JewellWitten (Jewell et al., 2020) by using their previously benchmarked performance on the same datasets obtained from extensive parameter optimization (Rupprecht et al., 2021) (Fig. 4C). As a metric for inference performance, we used the Pearson’s correlation coefficient between ground-truth spikes and predicted spikes after filtering with a Gaussian kernel with 200 ms bandwidth. This metric allowed us to directly compare PGBAR with previous analyses Rupprecht et al. (2021). The correlation obtained with PGBAR averaged across cells and datasets is 0.75. As for the other methods, we observed a large variability of performance across recordings, however we did not find a particular condition where our method performed systematically better or worse across indicators (Fig. 4A) and standardized noise level (Fig. 4B), which is a noise index introduced by Rupprecht et al. (2021) to account for different sampling frequencies (ratio between the standard deviation of the normalized fluorescence and the square root of the sampling frequency). In particular, we did not observe a statistically significant difference between mean correlation coefficient across standardized noise levels in the range 0.5-3.5 (Fig. 4B), although lower noise levels seem to be associated to a larger range of correlation.

Figure 4.

The performance range of PGBAR across the CASCADE database was comparable to previous methodologies. Among model-based approaches, our method slightly underperformed MLSpike (PGBAR median is within 6% from MLSpike, Mann-Whitney test p-value 0.0043), likely due to its description of the non-linear properties of the calcium indicator. The top performance is achieved by the supervised CASCADE method, however, it is a supervised method (i.e. needs training on ground truth data) and does not provide posterior distributions of spike times and model parameters.

The previous analysis provides an overview of the general performance of PGBAR in comparison with previous methods. To illustrate the advantages of PGBAR in estimating statistical uncertainties, we now take a closer look at a representative recording in the CASCADE dataset (Fig. 5A, CASCADE dataset 9 Chen et al. (2013)), a GCaMP6f fluorescence time series from a pyramidal neuron in the mouse visual cortex. The comparison between ground-truth spikes and those inferred using PGBAR in Fig. 5A shows differences outside the 1st-3rd interquartile range in 30% of the 1s time intervals. This statistical discrepancy between the posterior distributions and the ground truth can be attributed to the limitations of the autoregressive model to reproduce the biophysical properties of GCaMP6. In spite of such discrepancy, the estimated bursting pattern shown in Fig. 5A (lower panel) captures faithfully the overall periods of increased neuronal activity. The posterior distributions of model parameters are shown in Fig. 5B. Some distributions (burst firing rate and the rise time) shift significantly from their corresponding priors. The mismatch between prior and posterior also highlights the shortcomings of the model. In this particular case, the posterior distribution of the burst firing rate appears shifted to larger values. This result could be associated to the use of Gaussian noise in the model. Indeed, non Gaussian fluorescence fluctuations might be erroneously interpreted by the sampler as true spikes, pushing the burst firing rate to larger values. This example illustrates that prior distributions can be used to partially compensate this bias by penalizing unrealistic parameter configurations.

Figure 5.

2.6 Validation and performance of PGBAR on simulations of short-interval stimuli

We tested the accuracy of PGBAR in resolving pairs of stimuli with inter-spike intervals (ISI) below 10ms by performing model simulations of different signal-to-noise levels. Figure 6A shows two simulated fluorescence traces with two spikes 10ms apart at low (1.4) and high (3.4) SNR levels. By running our algorithm on these simulated data PGBAR was accurately detecting 2 spikes with over 95% confidence across all SNR values. Next we focus on the temporal accuracy of the ISI estimation. Increasing SNR or sampling frequency has the effect of narrowing the ISI posterior distribution around the ground-truth value (Fig. 6B). When examining multiple stimulus intervals, SNRs, and sampling frequencies we show that high sampling frequeuncies enable more reliable extraction of spike intervals. Figure 6C shows that the posterior probability of the ISI falls within 3-6ms of the ground-truth ISI for low (1kHz) and high (3kHz) sampling frequencies, respectively (averaged over 30 trials). When this probability is above 0.5, the ground-truth value is equal to the mode of the posterior ISI distribution (using 3ms binsize), therefore we can use this probability as a metric to quantify the temporal accuracy of our inference. For ground-truth values above 3ms, the contours of constant probability (iso-probability) show that increasing the temporal separation between stimuli has no effect on the inference accuracy. On the contrary, increasing the sampling frequency shifts the iso-probability contours towards lower SNR, indicating that for high sampling rates, lower SNRs are required to reach the same accuracy. Figure 6D shows the posterior ISI distributions obtained from ten independent simulations with ground-truth ISI of 5ms at different SNR levels. At low SNR the posterior distributions have larger variance while at higher SNR levels they shrink around the ground-truth value. This analysis illustrates the variability expected when analyzing multiple trials of the same neuron. We have shown that PGBAR can provide accurate estimates of inter-spike intervals down to 5 ms but the trade-off between temporal accuracy, SNR and sampling frequency must be considered.

Figure 6.

2.7 PGBAR spike inference from fluorescence transients recorded using the fast calcium indicator GCaMP8f

We tested our approach on the fast calcium indicator GCaMP8f by performing high-speed (≈ 2.8kHz) 2-photon linescan calcium imaging of cerebellar granule cells in vitro. We used adeno-associated viruses (AAV) to express GCaMP8f in Crus I of the cerebellum (Fig. 7A). Fluorescence was recorded for both granule cell somata and synaptic boutons while ground truth spikes were evoked by extracellular stimulation of granule cell axons in the molecular layer (Fig. 7B).

Figure 7.

First, we obtained prior distributions of amplitudes, rise and decay times of single action potential-evoked GCaMP8f fluorescence responses (Fig. 7C) to constrain the inference procedure. Next, we recorded fluorescence responses in a granule cell soma in response to a 20 Hz Poisson stimulation protocol; both a single trial and several trials are shown. (Fig. 7E). We used PGBAR to analyze independently each trial in Figure 7E. By generating thousands of posterior samples of spike time patterns, we obtained the spike probability for all time frames and trials (Fig. 7F). The spike patterns obtained using our method are very similar across trials, showing that PGBAR can reliably detect single-trial action potential-evoked GCaMP8f fluorescence transients.

To illustrate the temporal accuracy of PGBAR we focused on a part of the stimulus train with a short 5 ms interval between two spikes (Fig. 7G). Despite the relatively low SNR (A^(max)/σ ≈ 2.4), we could reliably identify the two spikes in each trial. The ground-truth ISI is well within the posterior distribution of each trial, with posterior modes symmetrically distributed around it (Fig. 7H). To better understand the variability of the posterior modes, we simulated 10 fluorescence traces with two spikes separated by 5 ms and sampling frequency matching our recordings. By analysing these simulations with PGBAR we obtained a distribution of ISI posterior modes that is statistically compatible with the real data (Fig. 7I, Mann-Whitney test p-value = 0.97), which provides further evidence that our results are statistically unbiased and that PGBAR can be used to infer spike times for intervals down to 5 ms.

Genetically encoded calcium indicators have been used to monitor neuronal activity in somata and boutons. We examine the accuracy of of PGBAR to infer spike times in 2.8 kHz linescans. Fluorescence transients were evoked using a 20Hz Poisson stimulation protocol containing 29 stimulation events. Each trial was analyzed independently using PGBAR. The distribution of the posterior modes of the total number of spikes across experiments and trials is centered around the ground-truth value for both somas and boutons with a relative error of 15% (Fig. 8A). The time to nearest spike detection, averaged over all stimulation time points, is 0.43 ms (± 0.10 ms, SEM) and 1.39 ms (± 0.11 ms, SEM) in boutons and somas respectively, highlighting a significant delay of 0.95 ms (± 0.15 ms SEM, Mann-Whitney test p-value = 4.8 × 10⁻⁷) from the average time of detection in boutons (Fig. 8B). To examine the temporal accuracy of our approach we calculated the Pearson’s correlation between detected spike times and stimulation times over different time-scales. We binned detected spike times and the Poisson stimulation train using bin sizes ranging from 2 ms to 50 ms and calculated their correlation. If we consider as reference correlation the average across the CASCADE dataset (0.75), this analysis shows that the correlation is larger than 0.75 up to 10ms bin size, confirming that our method provides accurate inference of short intervals.

Figure 8.

3 Discussion

Fluorescence indicators provide essential tools for monitoring the activity of neuronal populations in model organisms. However, the extraction of underlying firing patterns from fluorescence time series is challenging due to low signal-to-noise ratio, incomplete knowledge of the indicator dynamics, complex firing statistics, and unknown fluorescence modulation. In spite of the proliferation of methodologies developed to address this issue, limited attention has been devoted to the estimation of the statistical uncertainties associated to spike inference. The vast majority of the spike detection algorithms are indeed based on optimization techniques, providing only point estimates of the detected spikes. Quantifying statistical uncertainties is key to compare firing patterns across neurons (Diana et al., 2019) and establishing their causal relationships. The work of Pnevmatikakis et al. (2013) and Vogelstein et al. (2009), addressed this issue by using Monte Carlo methods to approximate the full posterior distributions of spiking patterns. Building upon their work, here we improve the generative model used to infer spiking patterns from fluorescence time series and describe efficient Monte Carlo strategies to infer spike times and their statistical uncertainties.

Bursting dynamics and baseline modulation

Neural activity patterns are not always well described by a simple Poisson spiking process. The PGBAR is the first Monte Carlo method (Greenberg et al., 2018; Pnevmatikakis et al., 2013; Vogelstein et al., 2009) to perform statistical inference based on a non-homogeneous Poisson firing and baseline modulation model. PGBAR uses a two-state process to enable transitions between low and high firing rates. This feature is used to mimic the alternation between periods of low baseline firing and bursting activity transients where the firing rate is significantly increased. We have shown that not taking into account bursting activity in the model used for inference can lead to biased results, especially at low SNR levels and high firing frequency. By explicitly modeling the bursting dynamics, PGBAR produced unbiased results at all levels of noise and frequency when tested on simulated data (Figure 3).

Although the model used in Pnevmatikakis et al. (2013) did not account for fluorescence baseline modulation, they acknowledged that it would be important to account for baseline fluctuations in in vivo recordings. PGBAR uses a Gaussian random walk, as MLSpike Deneux et al. (2016)), to describe slow changes in baseline fluorescence across the recording. While this is a simple Markov model of fluorescence baseline it introduces additional noise. To avoid this effect, it is possible to employ alternative baseline models, such as the integrated random walk.

Joint estimation of time-independent parameters and dynamic variables

The estimation of time-independent model parameters is a well known issue in spike detection algorithms, requiring ad-hoc calibration procedures and manual settings. In particular, unknown firing rate and peak amplitude (in response to a single action potential) can lead to unidentifiable parameters. PGBAR employs a fully Bayesian approach where model parameters and dynamic variables are treated equally. This enables users to constrain the detection of action potential by controlling mean and variance of phenomenological parameters (e.g. rise and decay constants, firing rates, bursting frequencies) by constraining their prior distributions. We employed for the first time state-of-the-art particle Gibbs algorithms to infer spikes from noisy fluorescence. This is a key novelty compared to previous SMC-based methods (Greenberg et al., 2018; Vogelstein et al., 2009), allowing for a joint estimation of time-independent parameters and dynamic variables.

Comparison with benchmark datasets

The proliferation of spike inference methodologies led to the development of community-based initiatives (Berens et al., 2018; Theis et al., 2016) to rank the performance of available methods. We applied our approach on the CASCADE dataset (Rupprecht et al., 2021) which provides a curated database of neuronal recordings from mouse and zebrafish using different calcium indicators. The performance of PGBAR, measured by the correlation coefficient with ground-truth spikes, is within one quartile from the correlation distributions of existing unsupervised approaches. In addition, it provides information about the statistical uncertainty associated to spike detection that is not currently possible with state-of-the-art techniques.

PGBAR detection of short high-frequency bursts using an ultrafast calcium indicator

PGBAR employs a second-order autoregressive process to link spiking activity to fluorescence. This simple model accounts for the basic qualitative aspects of calcium transients, and it is well-suited for linear indicators. For this reason, we have tested the performance of PGBAR on the ultrafast GCaMP8f (Zhang et al., 2021) that exhibits improved linearity in comparison to previous calcium probes. We showed that the combination of PGBAR with the GCaMP8f enables the detection of inter-spike intervals of 5 ms with an accuracy of 2.5 ms from single trials, thus offering a statistical tool for estimating high-frequency neural activity patterns.

PGBAR limitations and future perspectives

Although full Bayesian inference is known to be computationally expensive, SMC algorithms are highly parallelizable. Posterior distributions are represented by particles that are simultaneously propagated through time. In particular, GPU parallelization of SMC methods is an active field of research in computational statistics. Future advances in this directions might dramatically boost these methods and offer tools for online processing of fluorescence time series.

Many commonly used indicators exhibit a nonlinear fluorescence response to increases in intracellular calcium (Chen et al., 2013). The autoregressive model used in this study does not account for that nonlinearity. This work provides a statistical framework generalizable to more specific biophysical models of calcium indicators to account for non-linear effects. In addition, these models depend on rate constants that are difficult to measure directly. Not being able to constrain model parameters is an issue for spike inference as different parameter combinations might describe equally well the observed fluorescence by adjustments of the inferred spike pattern. For this reason, when model parameters are not identifiable, the problem of spike inference from fluorescence time series is ill-defined. Our Bayesian framework offers a systematic approach to address this issue by exploiting current and future data on the kinetics of calcium indicators through informative priors. Thanks to our joint sampling strategy of model parameters and spike patterns, our Monte Carlo method can keep the inference within biophysically relevant regions.

4 Materials and Methods

4.1 Particle Gibbs with ancestor sampling

The PGAS step used in Algorithm 1 to sample latent state trajectories was introduced in Lindsten et al. (2014) to improve the performance of the original particle Gibbs sampler (Andrieu et al., 2010). We refer to their original works for details on convergence and mixing properties of the method. The PGAS step is a SMC algorithm that generate a new latent state trajectory starting from a reference trajectory and the model parameters. We initialize the algorithm with a set of N latent states (particles) at time t = 1,

where the first N − 1 are sampled from a proposal distribution equal to the probability of the initial state μ^θ(X₁)

where bursting and baseline firing states q₁ = 0, 1 have equal probability and the calcium vector is constrained by the initial condition of Eq. (5), Ĉ₁ ≡ (c₀ + AS₁, 0) and the model parameter c₀. The last particle is constrained by the reference trajectory X^l. To conclude the initialization stage, we assign importance weights to all particles

Next, for t > 1, we evolve the particle system through time by assigning ancestor particles to time t to get a new set of particles and assigning weights . For the first N − 1 particles the ancestors are obtained by multinomial resampling from the particle system at time t 1 with probability proportional to the importance weights . The ancestor J of the last particle is drawn from the distribution

exploiting the fact that we know its state at time t.

In order to evolve the particle system through time we need to set a proposal distribution ρ_t(X_t | X_t−1) to sample new latent states. This proposal is then taken into account in the reweighting stage. Although the choice of the proposal distribution is arbitrary, it can be shown that the conditional distribution P (X_t | X_t−1, F_t) reduces the variance of the importance weights. We can express this optimal proposal as

where Z^θ(X_t−1, F_t) is the normalization factor as a function of the latent state at time t −1 and current fluorescence F_t. We can now use the expressions of and in Eqs. (11,12) for our model to compute the optimal proposal distribution. To do so we decompose P (X_t|X_t−1, F_t) as the product

where we used the fact that C_t is deterministic and only depends on C_t−1 and the spike count at time t. The idea is to use this chain decomposition to sample first the firing state q_t and the spike count s_t, then calculate C_t from its deterministic evolution and finally sample the baseline b_t from its distribution conditional to the other variables. The first term P (q_t, s_t|X_t−1) can be obtained by integrating the product over b_t and C_t and then normalizing the result. The integration over b_t and C_t leads to

where we introduced the function as the integral

The normalization factor Z^θ(X_t−1, F_t) is obtained by taking the sum of Eq. (20) over firing state and spike count:

To draw a combination of q_t and s_t from this distribution we applied a cutoff to the number of spikes per time step S^(max) = 20 and constructed a probability matrix of size 2 S^(max) for all combinations of firing state and spike count.

To obtain the full conditional distribution of b_t we consider again the product and by keeping only terms in b_t and normalizing we obtained a Gaussian distribution with mean μ_prop and variance σ_prop given by

The final step is to reweight all particles according using the importance weight

However, due to the form of the optimal proposal in Eq. (18) the importance weights reduce to the normalization factor calculated in Eq. (22)

Algorithm 2

PGAS kernel

4.2 Prior distributions

As discussed in the text, we use a reparameterization of the autoregressive model in terms of the maximal amplitude A^(max), rise and decay times τ_r and τ_d, for which it is easier to design realistic prior distributions based on previous empirical estimates of the kinetics of calcium indicators. We have used truncated normal priors for the maximal amplitude, the initial condition of the autoregressive model c₀, rise and decay time. In order to calculate the full conditional distributions on bursting/baseline firing rates r_0,1 and the transition matrix parameters w_{q→q ′} we have used a gamma distribution, whereas for the noise level σ² an inverse gamma distribution.

4.3 Sampling rules for time-independent parameters

In Algorithm 1, after a new latent state trajectory is sampled from the PGAS kernel, we draw time-independent model parameters from the conditional distribution P (θ_i | X_1:T, F_1:T). We use a mixed approach where the parameters r_0,1, w_q→q′ and σ² are sampled from their full conditional distribution, which can be obtained analytically by using gamma priors, while kernel parameters, A^(max) and τ_r,d, are sampled using the Metropolis-Hastings acceptance rule.

The full conditional distribution of a given parameter can be obtained from the joint probability of model parameters, latent state and fluorescence trajectories

where P (θ) is the prior distribution. We will now calculate the full conditionals for firing rates r_q, transition parameters w_{q→ q′} and noise variance σ². For simplicity we will use the same symbols for shape, α, and rate, β of all prior distributions, although they differ numerically for each parameter. By combining the expressions in Eqs.(10), (11) and (12) with the gamma prior gamma(α, β) and by keeping only terms proportional to r_0,1 we obtain

therefore the full conditional is a gamma distribution with updated parameters

By applying the same method to the transition rates w_{q→ q′} we have

where the approximation holds when the transition rate between firing states is much slower than the sampling frequency (w_{q→ q′} « ∆⁻¹). Therefore the full conditional is again a gamma distribution with parameters

For the noise variance parameter we used an inverse gamma prior and by applying the same method we can compute the full conditional as

therefore the updated shape and rate of the inverse gamma are

4.4 Response kernel

The response to a single spike can be obtained by writing Eq. (3) in the form of a first order Markov process in terms of the new variables C_t ≡ [c_t, c_t−1] and S_t ≡ [s_t, 0] so that

with the initial condition C₁ = [c₀ + AS₁, 0]. We can now write the solution at time t as

If s_t = δ_t,1 and c₀ = 0 then Eq. (38) simplifies to

By introducing eigenvectors and eigenvalues of M

we can express Eq. (39) as

therefore we have

By setting the time derivative of c_t to zero we obtain the time to reach the maximal response τ_r as

whereas if we take the long time limit of Eq. (39) we obtain

4.5 Reparameterization

To reparameterize the autoregressive model in terms of kinetic parameters we need to find the inverse map τ_r,d → γ_1,2 to obtain the original autoregressive parameters given rise and decay times. By combining the expressions of τ_r and τ_d in terms of g_± we have

which shows that the ratio g₋/g₊ can be expressed as

where the inverse function f ⁻¹(x) can be determined by numerical interpolation in the range [0,1]. To obtain the original autoregressive parameteres γ_1,2 we first obtain g_± as

and then

4.6 Experimental methods

GCaMP8f virus injection in the Crus I lobule of the cerebellum

Virus injection was targeted to Crus I (6.00 mm posterior to the Bregma, 3.00 mm lateral to the midline; vertical depth of 0.50 mm from pial surface) under deep isofluorane anesthesia. 100 nl of adeno-associated virus encoding GCaMP8f (AAV-DJ-CAG-FLEx-jGCaMP8f, Janelia Research Campus) was injected with Nanoject III (Drummond Scientific) using thin glass pipette (diameter 30 μm). Injection was made when the C57BL/6J-Gabra6tm2(cre)Wwis mice were 7 weeks old. After injections the mice were returned to their home cage for 6 weeks to allow time for expression.

Slice Preparation

Acute coronal slices (200 μm) of the Crus I lobule of the cerebellum was prepared from adult C57BL/6J-Gabra6tm2(cre)Wwis mice, aged 97 and 102 days. Following transcardial perfusion with an ice-cold solution containing (in mM): 2.5 KCl, 0.5 CaCl2, 4 MgCl2, 1.25 NaH2PO4, 24 NaHCO3, 25 glucose, 230 sucrose, and 0.5 ascorbic acid, the brains were removed and placed in the same solution. The solution was bubbled with 95% O2 and 5% CO2. Slices were cut from the dissected lateral cerebellum using a vibratome (Leica VT1200S), and incubated at room temperature for 30 minutes in a solution containing (in mM): 85 NaCl, 2.5 KCl, 0.5 CaCl2, 4 MgCl2, 1.25 NaH2PO4, 24 NaHCO3, 25 glucose, 75 sucrose and 0.5 ascorbic acid. Slices were then transferred to an external recording solution containing (in mM): 125 NaCl, 2.5 KCl, 1.5 CaCl2, 1.5 MgCl2, 1.25 NaH2PO4, 24 NaHCO3, 25 glucose and 0.5 ascorbic acid, and maintained at room temperature for up to 6 hr.

Cellular Imaging

The brain slices containing GCaMP8f expressing cells were identified with a 4x objective lens (Olympus UplanFI 4x, 0.13 NA) using very brief illumination with 470 nm light to excite GcaMP8f fluorescence. GCs were identified using infrared Dodt-gradient contrast and a QlClick digital CCD camera (QImaging, Surrey, BC, Canada) mounted on an Ultima multiphoton microscopy system (Bruker Nano Surfaces Division, Middleton, WI, USA) that was mounted on an Olympus BX61W1 microscope, equipped with a water-immersion objective (Olympus 60x, 1.1 NA). Two-photon excitation was performed with a Ti-sapphire laser (Spectraphysics). To visualize GCs expressing GcaMP8f, two-photon excitation was performed at 920 nm. Infrared Dodt-gradient contrast was used to position the stimulation pipette in molecular layer targeting PF to directly activate GC bodies. Linescan imaging of GC bodies was performed by scanning through the whole cell membrane marked by freehand linescan mode (Prairie View). Total laser illumination per sweep lasted 2000 ms. Fluorescence was detected using both proximal epifluorescence and substage photomultiplier tube gallium arsenide phosphide (H7422PA-40 SEL, Hamamatsu).

5 Code and Data

In this work we employed the publicly available CASCADE dataset in Rupprecht et al. (2021). Linescan imaging data and the source code code of PGBAR are available through the GitHub repository https://github.com/giovannidiana/pgbar.

Acknowledgements

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 896051. This work was also supported by a Pasteur-Roux-Cantarini fellowship of the Institut Pasteur. We would like to thank Peter Rupprecht for his help with the analysis of the CASCADE dataset and for thorough discussions. We also thank Nicolas Chopin for suggesting the use of backward steps methods in particle Gibbs, Andrea Giovannucci, Marco Banterle and Diana Passaro for helpful discussions.