Database and deep learning toolbox for noise-optimized, generalized spike inference from calcium imaging

<jats:title>ABSTRACT</jats:title><jats:p>Calcium imaging is a key method to record patterns of neuronal activity across populations of identified neurons. Inference of temporal patterns of action potentials (‘spikes’) from calcium signals is, however, challenging and often limited by the scarcity of ground truth data containing simultaneous measurements of action potentials and calcium signals. To overcome this problem, we compiled a large and diverse ground truth database from publicly available and newly performed recordings. This database covers various types of calcium indicators, cell types, and signal-to-noise ratios and comprises a total of &gt;35 hours from 298 neurons. We then developed a novel algorithm for spike inference (CASCADE) that is based on supervised deep networks, takes advantage of the ground truth database, infers absolute spike rates, and outperforms existing model-based algorithms. To optimize performance for unseen imaging data, CASCADE retrains itself by resampling ground truth data to match the respective sampling rate and noise level. As a consequence, no parameters need to be adjusted by the user. To facilitate routine application of CASCADE we developed systematic performance assessments for unseen data, we openly release all resources, and we provide a user-friendly cloud-based implementation.</jats:p>

section_name: Supplementary Figures

section_name: Supplementary Table

section_content: To ensure reliable inference of spike rates it is advantageous to train the supervised deep network with a training dataset that matches the noise level of the test neuron from the population imaging data. In practice, the noise level of each neuron from the population imaging data is determined and an existing model trained with approximately the same noise levels is loaded for spike inference. Noise-suppression is only effective if the noise statistics of the ground truth used for training the network resemble the noise statistics of the calcium data used for testing. To generate a ground truth that matches this requirement, we tested two different approaches to increase the noise of ground truth recordings to match population recordings (Fig. S2 ). 

First, we used the raw imaging data to extract not only the mean fluorescence trace, but the fluorescence trace of each pixel of the region of interest (ROI) that defines the neuron. To achieve a given noise level ν, a random subset of pixels was drawn from the ROI pixels until the average fluorescence trace of this sub-ROI reached the desired noise level (Fig. S2a ). This method generates realistic noise characteristics through spatial subsampling but is computationally costly and cannot be applied to ground truth datasets when raw fluorescence movies are not available. 

Alternatively, we extracted only the mean fluorescence trace of the ROI of a ground truth neuron and added artificial noise until the overall high-frequency noise matched the test calcium dataset (Fig. S2a ). This procedure can in theory be repeated to produce an infinite number of examples (replicas) from a simple ground truth recording. However, since the mean ΔF/F of a ground truth recording already is associated with a certain noise level, these noise patterns would be correlated across replicas. To avoid this undesired effect, which could lead to overfitting of correlated noise during training, the number of replicas was restricted to a number n that was computed with the noise level of the mean ΔF/F of a neuronal ROI, ν, and the target noise level, ν, by ν ν ⁄ and thresholded at a maximum of n 500. We tested both simple Gaussian noise as well as Poisson noise, where the variance of the noise is proportional to the signal amplitude, as is typical for photon shot noise. 

To test whether artificial noise enables the network to learn and suppress natural noise patterns we generated ground truth with either natural (spatial sub-sampling) or artificial (Gaussian or Poisson) noise for a subset of the available ground truth datasets (Fig. S2a-c ). Using models trained on artificial rather than natural noise slightly but significantly decreased the correlation of predictions with the ground truth when applied to test datasets based on natural noise (decrease for Gaussian noise: Δ 1 = 0. 010 ± 0. 004, pseudo-median ± 95% C. I., p < 1e-4, paired Wilcoxon test; decrease for Poisson noise: Δ 1 = 0. 006 ± 0. 004, p < 0. 005). This decrement indicates how much better a model trained with naturally sampled noise would likely perform. Conversely, testing models trained with artificially generated ground truth on artificially instead of naturally sampled ground truth increased the correlation of predictions slightly (Gaussian noise: Δ 2 = 0. 006 ± 0. 008, p = 0. 13; Poisson noise: Δ 2 = 0. 005 ± 0. 009, p = 0. 28). This increment indicates how much the correlations with the ground truth will be overestimated when training and testing is done with artificially generated ground truth only. Differences were generally more pronounced for larger noise levels, but overall remained very small compared to the absolute values (not shown). We therefore conclude that artificial noise allows deep networks to effectively learn noise patterns that can be applied to natural noise, with only minor performance loss compared to spatially subsampled recordings, and we further conclude that using Poissondistributed noise yields slightly improved performance compared to artificial Gaussian noise.
section_name: NOISE-MATCHING OF RESAMPLED GROUND TRUTH DATA

section_name: DEPENDENCE OF PERFORMANCE ON HYPER-PARAMETERS, OVERFITTING AND NETWORK

section_content: We tested how spike inference performance depends on the choice of hyper-parameters and network architecture. Networks were trained on a specific ground truth dataset using all neurons except one, which was held out for testing (leave-one-out strategy). The algorithm turned out highly robust with respect to changes of the optimizer for gradient descent, the batch size during learning, the number of convolutional features per layer, the number of neurons in the dense layer, and the extent of the temporal window of the receptive field (Fig. S4a-e ). All of these observations were confirmed over a surprisingly large range, indicating that the network performance is highly robust with respect to any hyper-parameter choices. 

In addition, we investigated potential overfitting of the training dataset. The performance of the network was high already after one training epoch (i. e., as soon as every sample had been seen once by the network), then reached a maximum after 10-30 training epochs, and slightly decreased thereafter (Fig. S4f ). This learning behavior suggests only moderate overfitting. At the same time the training loss decreased monotonically (Fig. S4g ). We believe that the high abundance of noise and sparseness of events acts as a natural regularizer that prevents overfitting. More importantly, while the learning curve was smooth on average (Fig. S4f ), individual network instances sometimes reached unfavorable states. As expected from known properties of deep networks 77, this effect could be easily eliminated by ensemble averaging over 5 networks (Fig. S4h ). 

Finally, we also tested the performance when employing a network architecture different from the standard convolutional architecture ('default'). We tested a large variety of standard deep learning architectures, including recurrent LSTM networks and non-convolutional deep networks with only the input, output and loss function remaining unchanged (see Methods for detailed descriptions and explanations). Most of these networks performed well, and the performance of several networks was statistically indistinguishable: the default convolutional network, a convolutional network with reduced filter size, a locally connected network, and a bi-directional LSTM network (Fig. S5a ). Much smaller networks (single convolutional layer, Fig. S8a ) tended to underfit the ground truth. On the other hand, networks with larger numbers of parameters (the deeper convolutional networks and the locally connected network) overfitted the data when training continued (dashed lines in Fig. S5b ), consistent with previous observations 16. 

The locally connected network and the bi-directional LSTM network performed equally well compared to the default convolutional network despite very different architectures. However, some architectures that were not adapted to spike inference showed lower performance, for example the naïve LSTM network, which by its recurrent design prevents the network from looking precisely at the time point of interest (see Methods for details). Another example is a network identical to the default convolutional network but with purely linear activation functions (Fig. S5a, b ), which prevents the algorithm from non-linearly adjusting decision boundaries.
section_name: ARCHITECTURE

section_content: Many existing spike inference methods do not aim at the inference of spike rates, but rather of discrete spikes 11, 12, 15. Previous publications reported that the average ΔF/F value triggered by a single spike is larger than zero, and that the calcium transients corresponding to single spikes can be detected in selected neurons 34. However, the precise identification of individual spikes in practical situations is more challenging since a detection scheme should also work in unseen data. In addition, the task may depend in unknown ways on the variable expression of calcium indicators, on shot noise, on other noise sources, on low sampling rates and on the non-linear response of calcium indicators 3, 22, 30, 31, 78. It is therefore not clear whether discrete spikes can be reliably inferred in realistic scenarios. We therefore devised two approaches to test whether single-spike precision can be achieved. First, we focused on single, isolated spikes in the ground truth and compared them with inferred spike rates for the same time window. This approach is discussed in the main text. Second, we transformed spike rates into discrete spikes and analyzed whether the discretization improved spike inference. 

With respect to the second approach, we argue that the spike rates inferred by the deep network will exhibit a tendency to quantize if the predictions are close to single-spike precision. Therefore, a procedure that takes into account the prior about discretized spiking could improve the inferred spike rates. We therefore applied an algorithmic procedure that uses prior knowledge about the spike rate (spiking probability) waveform associated with a single spike to fill up the inferred probability trace with discrete spikes using an optimization procedure based on Monte-Carlo importance sampling (Methods). 

We found that spiking probabilities that were almost correctly inferred by the deep network were optimized by suppression of noise or by rounding of close matches (blue arrows in Fig. S7a-d ). However, this procedure can also enhance small false positive or small negative errors (red arrows in Fig. S7a-d ). Over all ground truth datasets, the correlation metric slightly but consistently decreased when spike rates were discretized (Fig. S7e ), indicating that the data quality did not allow for efficient use of the prior. Although discretized spike rates tended to decrease the error (Fig. S7e ), this effect was primarily due to suppression of small noise events in the absence of spiking, and we found that this positive effect could be achieved without reduction of the correlation by thresholding the inferred spike rates (Fig. S7d, f ). Together, this suggests that the available datasets do not permit discretization of predicted spike rates without performance loss. Despite these caveats, discretization of spike rates might still be useful for two reasons. First, discrete spike events may be more intuitive visualizations of activity than smooth probabilities. Second, while spike rates are smoothed with a Gaussian kernel for each spike, the detection of a single spike that optimally explains this spike provides better temporal resolution. Here, we see a potentially useful application of discrete spike inference, which is, however, beyond the scope of this study. We include the algorithm to discretize spike rates as a script in the public repository (https://git. io/JtZe4). for an isolated action potential in the ground truth across all datasets. Across all datasets, the estimate (white) of a single ground truth action potential is clearly lower than expected from ground truth (yellow). h, Distribution of the overall area under the curve associated with a single isolated ground truth action potential. In an ideal case, values would be narrowly distributed around 1 (dashed line). Instead, the distribution is broad and biased with a central value <1. i, Spiking probability for a single dataset that should be ideally qualified for single-spike precision (DS#04; no movement artifacts due to ex vivo preparation, synthetic indicator OGB-1). Even here, the median prediction is on average systematically lower than a single action potential, indicating that the network trades off single-spike precision in order to prevent false positive detection of spikes amidst of shot noise.
section_name: DISCRETE SPIKES AND SINGLE-SPIKE PRECISION

section_content: We considered several characteristics of ground truth datasets and evaluated whether they could improve the mutual predictability across ground truth datasets. This includes characteristics that are accessible without ground truth (indicator species, i. e., synthetic dyes vs. genetically encoded indicators; animal species, i. e., mouse vs. zebrafish; median spike rate across neurons; the burstinesss; the Fano factor) and characteristics that are only accessible with available ground truth (area under the curve of the linear kernel, cf. Extended Data Fig. 1 ; decay constant of the linear kernel; the correlation between the kernels of the training and test datasets), with detailed descriptions in the Methods section. These 8 predictors were used as regressors for a GLM to fit the mutual predictability matrix (correlations) among datasets. a, Correlations predicted from the GLM vs. measured correlations (see Fig. 3a ). b, A GLM based on "all predictors" results in a correlation that recovers panel (a). In addition, this plot also shows values for fits to the bias and error matrices. 

Using only one of the predictors reduces the correlation, often very significantly. c, Measured cross-dataset predictability (reproducing Fig. 3a ). d, Cross-dataset predictability as computed with the GLM. Together, this supplementary figure shows that the GLM is not able to explain a large fraction of the variance of the original matrix, and that the main predictor is the correlation between kernel time courses, which is not accessible without available ground truth. As consequence, we chose the use of a "global EXC model" that was trained on all reliable datasets (Fig. 3a ). In principle, both re-assignment to past and future is equally undesirable. In some scenarios, however, re-assignment of the activity into the past of the spike can be particularly adverse, for example, when a precise external stimulus (e. g., an auditory tone) is used to generate a peri-stimulus activity trace. While the ground truth used for training (a) is typically smoothed symmetrically in time with a Gaussian function (b), by design activity that happened briefly after the stimulus will be deconvolved to time points both after and before the stimulus. To circumvent these a-causal events, it is possible to use a more causal filter (c), using a highly skewed inverse Gaussian distribution. d, The deep network can use the 'causal' ground truth, resulting in inferred spike rates that re-assign activity in a more causal way (lower prediction trace; see zoom-in). However, noise sources result in re-assignment of activity to time bins prior to ground truth spikes. e, The predictions become necessarily more sloppy and smeared out for higher noise levels. This shows that even when using a causal kernel to smooth the ground truth, a-causal effects can be induced by spike inference. Depending on the desired temporal resolution of the algorithm, different width of the smoothing kernels can be chosen. A Gaussian with σ = 0. 2 s will lead to a FWHM of the smoothing Gaussian of ca. 0. 5 s (b). For panel (d), a smoothing kernel with σ = 0. 2 s was used.
section_name: Figure S9 | Predicting cross-dataset predictability with a generalized linear model (GLM).

section_name: Supplementary Notes

section_content: Supplementary Note 1 | Noise-matching of resampled ground truth data. 

 Supplementary Note 2 | Dependence of performance on hyper-parameters, overfitting and network architecture. 

 Supplementary Note 3 | Discrete spikes and single-spike precision.
section_name: 