Direct extraction of signal and noise correlations from two-photon calcium imaging of ensemble neuronal activity

Abstract

Neuronal activity correlations are key to understanding how populations of neurons collectively encode information. While two-photon calcium imaging has created a unique opportunity to record the activity of large populations of neurons, existing methods for inferring correlations from these data face several challenges. First, the observations of spiking activity produced by two-photon imaging are temporally blurred and noisy. Secondly, even if the spiking data were perfectly recovered via deconvolution, inferring network-level features from binary spiking data is a challenging task due to the non-linear relation of neuronal spiking to endogenous and exogenous inputs. In this work, we propose a methodology to explicitly model and directly estimate signal and noise correlations from two-photon fluorescence observations, without requiring intermediate spike deconvolution. We provide theoretical guarantees on the performance of the proposed estimator and demonstrate its utility through applications to simulated and experimentally recorded data from the mouse auditory cortex.

Introduction

Neuronal activity correlations are essential in understanding how populations of neurons encode information. Correlations provide insights into the functional architecture and computations carried out by neuronal networks (Abbott and Dayan, 1999; Averbeck et al., 2006; Cohen and Kohn, 2011; Hansen et al., 2012; Kohn et al., 2016; Kohn and Smith, 2005; Lyamzin et al., 2015; Montijn et al., 2014; Smith and Sommer, 2013; Sompolinsky et al., 2001; Yatsenko et al., 2015). Neuronal activity correlations are often categorized in two groups: signal correlations and noise correlations (Cohen and Kohn, 2011; Cohen and Maunsell, 2009; Gawne and Richmond, 1993; Josić et al., 2009; Lyamzin et al., 2015; Vinci et al., 2016). Given two neurons, signal correlation quantifies the similarity of neural responses that are time-locked to a repeated stimulus across trials, whereas noise correlation quantifies the stimulus-independent trial-to-trial variability shared by neural responses that are believed to arise from common latent inputs.

Two-photon calcium imaging has become increasingly popular in recent years to record in vivo neural activity simultaneously from hundreds of neurons (Ahrens et al., 2013; Romano et al., 2017; Stosiek et al., 2003; Svoboda and Yasuda, 2006). This technology takes advantage of intracellular calcium flux mostly arising from spiking activity and captures calcium signaling in neurons in living animals using fluorescence microscopy. The observed fluorescence traces of calcium concentrations, however, are indirectly related to neuronal spiking activity. Extracting spiking activity from fluorescence traces is a challenging signal deconvolution problem and has been the focus of active research (Deneux et al., 2016; Friedrich et al., 2017; Grewe et al., 2010; Jewell et al., 2020; Jewell and Witten, 2018; Kazemipour et al., 2018; Pachitariu et al., 2018; Pnevmatikakis et al., 2016; Stringer and Pachitariu, 2019; Theis et al., 2016; Vogelstein et al., 2009; Vogelstein et al., 2010).

The most commonly used approach to infer signal and noise correlations from two-photon data is to directly apply the classical definitions of correlations for firing rates (Lyamzin et al., 2015), to fluorescence traces (Fallani et al., 2015; Francis et al., 2018; Rothschild et al., 2010; Winkowski and Kanold, 2013). However, it is well known that fluorescence observations are noisy and blurred surrogates of spiking activity, because of dependence on observation noise, calcium dynamics and the temporal properties of calcium indicators. Due to temporal blurring, the resulting signal and noise correlation estimates are highly biased. An alternative approach is to carry out the inference in a two-stage fashion: first, infer spikes using a deconvolution technique, and then compute firing rates and evaluate the correlations (Kerlin et al., 2019; Najafi et al., 2020; Ramesh et al., 2018; Soudry et al., 2015; Yatsenko et al., 2015). These two-stage estimates are highly sensitive to the accuracy of spike deconvolution, and require high temporal resolution and signal-to-noise ratios (Lütcke et al., 2013; Pachitariu et al., 2018). Furthermore, these deconvolution techniques are biased toward obtaining accurate first-order statistics (i.e. spike timings) via spatiotemporal priors, which may be detrimental to recovering second-order statistics (i.e. correlations). Finally, both approaches also undermine the non-linear dynamics of spiking activity as governed by stimuli, past activity and other latent processes (Truccolo et al., 2005). There are a few existing studies that aim at improving estimation of neuronal correlations, but they either do not consider signal correlations (Rupasinghe and Babadi, 2020; Yatsenko et al., 2015), or aim at estimating surrogates of correlations from spikes such as the connectivity/coupling matrix (Aitchison et al., 2017; Mishchenko et al., 2011; Soudry et al., 2015; Keeley et al., 2020).

Here, we propose a methodology to directly estimate both signal and noise correlations from two-photon imaging observations, without requiring an intermediate step of spike deconvolution. We pose the problem under the commonly used experimental paradigm in which neuronal activity is recorded during trials of a repeated stimulus. We avoid the need to perform spike deconvolution by integrating techniques from point processes and state-space modeling that explicitly relate the signal and noise correlations to the observed fluorescence traces in a multi-tier model. Thus, we cast signal and noise correlations within a parameter estimation setting. To solve the resulting estimation problem in an efficient fashion, we develop a solution method based on variational inference (Jordan et al., 1999; Blei et al., 2017), by combining techniques from Pólya-Gamma augmentation (Polson et al., 2013) and compressible state-space estimation (Rauch et al., 1965; Kazemipour et al., 2018; Ba et al., 2014). We also provide theoretical guarantees on the bias and variance performance of the resulting estimator.

We demonstrate the utility of our proposed estimation framework through application to simulated and real data from the mouse auditory cortex during presentations of tones and acoustic noise. In application to repeated trials under spontaneous and stimulus-driven conditions within the same experiment, our method reliably provides noise correlation structures that are invariant across the two conditions. In addition, our joint analysis of signal and noise correlations corroborates existing hypotheses regarding the distinction between their structures (Keeley et al., 2020; Rumyantsev et al., 2020; Bartolo et al., 2020). Moreover, while application of our proposed method to spatial analysis of signal and noise correlations in the mouse auditory cortex is consistent with existing work (Winkowski and Kanold, 2013), it reveals novel and distinct spatial trends in the correlation structure of layers 2/3 and 4. In summary, our method improves on existing work by: (1) explicitly modeling the fluorescence observation process and the non-linearities involved in spiking activity, as governed by both the stimulus and latent processes, through a multi-tier Bayesian forward model, (2) joint estimation of signal and noise correlations directly from two-photon fluorescence observations through an efficient iterative procedure, without requiring intermediate spike deconvolution, (3) providing theoretical guarantees on the performance of the proposed estimator, and (4) gaining access to closed-form posterior approximations, with low-complexity and iterative update rules and minimal dependence on training data. Our proposed method can thus be used as a robust and scalable alternative to existing approaches for extracting signal and noise correlations from two-photon imaging data.

Results

In this section, we first demonstrate the utility of our proposed estimation framework through simulation studies as well as applications on experimentally recorded data from the mouse auditory cortex. Then, we present theoretical performance bounds on the proposed estimator. Before presenting the results, we will give an overview of the proposed signal and noise correlation inference framework, and outline our contributions and their relationship to existing work. For the ease of reproducibility, we have archived a MATLAB implementation of our proposed method in GitHub (Rupasinghe, 2020) and have deposited the data used in this work in the Digital Repository at the University of Maryland (Rupasinghe et al., 2021).

Signal and noise correlations

We consider a canonical experimental setting in which the same external stimulus, denoted by [Formula/expression omitted. See PDF], is repeatedly presented across L independent trials and the spiking activity of a population of $N$ neurons are indirectly measured using two-photon calcium fluorescence imaging. Figure 1 (forward arrow) shows the generative model that is used to quantify this procedure. The fluorescence observation in the [Formula/expression omitted. See PDF] trial from the [Formula/expression omitted. See PDF] neuron at time frame $t$, denoted by $y_{t,l}^{(j)}$, is a noisy surrogate of the intracellular calcium concentrations. The calcium concentrations in turn are temporally blurred surrogates of the underlying spiking activity $n_{t,l}^{(j)}$, as shown in Figure 1.

Figure 1.

The proposed generative model and inverse problem.

Observed (green) and latent (orange) variables pertinent to the [Formula/expression omitted. See PDF] neuron are indicated, according to the proposed model for estimating the signal (blue) and noise (red) correlations from two-photon calcium fluorescence observations. Calcium fluorescence traces $\left(y_{t,l}^{(j)}\right)$ of $L$ trials are observed, in which the repeated external stimulus [Formula/expression omitted. See PDF] is known. The underlying spiking activity $\left(n_{t,l}^{(j)}\right)$ , trial-to-trial variability and other intrinsic/extrinsic neural covariates that are not time-locked with the external stimulus $\left(x_{t,l}^{(j)}\right)$ , and the stimulus kernel [Formula/expression omitted. See PDF] are latent. Our main contribution is to solve the inverse problem: recovering the underlying latent signal [Formula/expression omitted. See PDF] and noise [Formula/expression omitted. See PDF] correlations directly from the fluorescence observations, without requiring intermediate spike deconvolution.

In modeling the spiking activity, we consider two main contributions: (1) the common known stimulus [Formula/expression omitted. See PDF] affects the activity of the [Formula/expression omitted. See PDF] neuron via an unknown kernel [Formula/expression omitted. See PDF], akin to the receptive field; (2) the trial-to-trial variability and other intrinsic/extrinsic neural covariates that are not time-locked to the stimulus [Formula/expression omitted. See PDF] are captured by a trial-dependent latent process $x_{t,l}^{(j)}$. Then, we use a Generalized Linear Model to link these underlying neural covariates to spiking activity (Truccolo et al., 2005). More specifically, we model spiking activity as a Bernoulli process:

[Formula omitted. See PDF]

where $\varphi (\cdot )$ is a mapping function, which could in general be non-linear.

The signal correlations aim to measure the correlations in the temporal response that are time-locked to the repeated stimulus, [Formula/expression omitted. See PDF]. On the other hand, noise correlations in our setting quantify connectivity arising from covariates that are unrelated to the stimulus, including the trial-to-trial variability (Keeley et al., 2020). Based on the foregoing model, we propose to formulate the signal [Formula/expression omitted. See PDF] and noise [Formula/expression omitted. See PDF] covariance between the [Formula/expression omitted. See PDF] neuron and [Formula/expression omitted. See PDF] neuron as:

[Formula omitted. See PDF]

where ${\displaystyle \mathrm{cov}(\cdot ,\cdot )}$ is the empirical covariance function defined for two vector time series ${\displaystyle {\mathbf{u}}_t}$ and ${\displaystyle {\mathbf{v}}_t}$ as ${\displaystyle \mathrm{cov}\left({\mathbf{u}}_t,{\mathbf{v}}_t\right):=\frac{1}{T}\sum _{t=1}^T\left({\mathbf{u}}_t-\frac{1}{T}\sum _{t^{\mathrm{\prime }}=1}^T{\mathbf{u}}_{t^{\mathrm{\prime }}}\right){\left({\mathbf{v}}_t-\frac{1}{T}\sum _{t^{\mathrm{\prime }}=1}^T{\mathbf{v}}_{t^{\mathrm{\prime }}}\right)}^{\mathrm{\top }}}$, for a total observation duration of $T$ time frames.

Our main contribution is to provide an efficient solution for the so-called inverse problem: direct estimation of [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF] from the fluorescence observations, without requiring intermediate spike deconvolution (Figure 1, backward arrow). The signal and noise correlation matrices, denoted by [Character omitted. See PDF] and [Character omitted. See PDF], can then be obtained by standard normalization of [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF]:

[Formula omitted. See PDF]

We note that when spiking activity is directly observed using electrophysiology recordings, the conventional signal [Formula/expression omitted. See PDF] and noise [Formula/expression omitted. See PDF] covariances of spiking activity between the [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF] neuron are defined as (Lyamzin et al., 2015):

(3)

${\displaystyle ({\mathbf{\Sigma }}_s^{\mathsf{c}\mathsf{o}\mathsf{n}}{)}_{i,j}:=\mathrm{cov}\left(\frac{1}{L}\sum _{l=1}^L{n}_{t,l}^{(i)},\frac{1}{L}\sum _{l=1}^L{n}_{t,l}^{(j)}\right),({\mathbf{\Sigma }}_x^{\mathsf{c}\mathsf{o}\mathsf{n}}{)}_{i,j}:=\frac{1}{L}\sum _{l=1}^L\mathrm{cov}\left(n_{t,l}^{(i)}-\frac{1}{L}\sum _{l^{\mathrm{\prime }}=1}^L{n}_{t,l^{\mathrm{\prime }}}^{(i)},n_{t,l}^{(j)}-\frac{1}{L}\sum _{l^{\mathrm{\prime }}=1}^L{n}_{t,l^{\mathrm{\prime }}}^{(j)}\right),}$which after standard normalization in Equation 2 give the conventional signal [Formula/expression omitted. See PDF] and noise [Formula/expression omitted. See PDF] correlations. While at first glance, our definitions of signal and noise covariances in Equation 1 seem to be a far departure from the conventional ones in Equation 3, we show that the conventional notions of correlation indeed approximate the same quantities as in our definitions:${\displaystyle {\mathbf{S}}^{\mathsf{c}\mathsf{o}\mathsf{n}}\approx \mathbf{S}\text{and}{\mathbf{N}}^{\mathsf{c}\mathsf{o}\mathsf{n}}\approx \mathbf{N},}$under asymptotic conditions (i.e. $T$ and $L$ sufficiently large). We prove this assertion of asymptotic equivalence in Appendix 1, which highlights another facet of our contributions: our proposed estimators are designed to robustly operate in the regime of finite (and typically small) $T$ and $L$, aiming for the very same quantities that the conventional estimators could only recover accurately under ideal asymptotic conditions.

Existing methods used for performance comparison

In order to compare the performance of our proposed method with existing work, we consider three widely available methods for extracting neuronal correlations. In simulation studies, we additionally benchmark these estimates with respect to the known ground truth. The existing methods considered are the following:

Pearson correlations from the two-photon data

In this method, fluorescence observations are assumed to be the direct measurements of spiking activity, and thus empirical Pearson correlations of the two-photon data are used to compute the signal and noise correlations (Rothschild et al., 2010; Winkowski and Kanold, 2013; Francis et al., 2018; Bowen et al., 2020). Explicitly, these estimates are obtained by simply replacing $n_{t,l}^{(j)}$ in Equation 3 by $y_{t,l}^{(j)}$, without performing spike deconvolution.

Two-stage Pearson estimation

Unlike the previous method, in this case spikes are first inferred using a deconvolution technique. Then, following temporal smoothing via a narrow Gaussian kernel the Pearson correlations are computed using the conventional definitions of Equation 3. For spike deconvolution, we primarily used the FCSS algorithm (Kazemipour et al., 2018). In order to also demonstrate the sensitivity of these estimates to the deconvolution technique that is used, we provide a comparison with the f-oopsi deconvolution algorithm (Pnevmatikakis et al., 2016) in Figure 2—figure supplement 1.

Two-stage GPFA estimation

Similar to the previous method, spikes are first inferred using a deconvolution technique. Then, a latent variable model called Gaussian Process Factor Analysis (GPFA) (Yu et al., 2009) is applied to the inferred spikes in order to estimate the latent covariates and receptive fields. Based on those estimates, the signal and residual noise correlations are derived through a formulation similar to Equation 1 and Equation 2 (Ecker et al., 2014).

Simulation study 1: neuronal ensemble driven by external stimulus

We simulated calcium fluorescence observations according to the proposed generative model given in Proposed forward model, from an ensemble of $N=8$ neurons for a duration of $T=5000$ time frames. We considered $L=20$ repeated trials driven by the same external stimulus, which we modeled by an autoregressive process (see Guidelines for model parameter settings for details). Figure 2 shows the corresponding estimation results.

Figure 2.

Results of simulation study 1.

(A) Estimated noise and signal correlation matrices from different methods. Rows from left to right: ground truth, proposed method, Pearson correlations from two-photon recordings, two-stage Pearson estimates and two-stage GPFA estimates. The normalized mean squared error (NMSE) of each estimate with respect to the ground truth and the leakage effect quantified by the ratio between out-of-network and in-network power (leakage) are indicated below each panel. (B) Simulated external stimulus (orange), latent trial-dependent process (red), fluorescence observations (black), estimated calcium concentrations (purple), putative spikes (green), and estimated mean of the latent state (blue) by the proposed method, for the first trial of neuron 1.

Figure 2—figure supplement 1.

Sensitivity of two-stage estimates to the choice of the underlying spike deconvolution technique.

(A) Noise (first row) and signal (second row) correlations corresponding to the ground truth (first column), estimated by the two-stage Pearson method using the FCSS (Kazemipour et al., 2018) (second column) and constrained f-oopsi (Pnevmatikakis et al., 2016) (third column) spike deconvolution techniques, for the simulation study in Figure 2. The NMSE and leakage ratios of the estimates are indicated below each panel. While the correlation estimates based on these two methods are comparable, there exist notable differences between them, as a result of the slight discrepancies in the deconvolved spikes. This demonstrates that the two-stage estimates are sensitive to minor differences in the estimated spikes obtained by different deconvolution techniques. In addition, both two-stage Pearson estimates fail to capture the ground truth correlations (as is also evident from the high NMSE and leakage values). (B) Simulated observations (black, re-scaled for ease of visual comparison) and ground truth spikes (blue), as well as the estimated calcium concentrations (purple) and putative spikes (green) for the 1st trial of neuron one in the simulation study of Figure 2, using the FCSS (Kazemipour et al., 2018) (second row) and constrained f-oopsi (Pnevmatikakis et al., 2016) (third row) spike deconvolution methods.

Figure 2—figure supplement 2.

Performance of two-stage estimates based on ground truth spikes.

Performance of two stage estimates based on ground truth spikes. Noise (first row) and signal (second row) correlations corresponding to the ground truth (first column) are repeated from Figure 2. The second and third columns show the results of two-stage GPFA and two-stage Pearson methods using $L=20$ trials, respectively. The fourth column shows the results of the two-stage Pearson method using $L=1000$ trials. All estimates were obtained using the ground truth spikes, as opposed to extracting the spikes via a deconvolution technique. Thus, these results isolate the effect of the non-linearities involved in spike generation on the estimation performance. The NMSE and leakage ratios of the estimates are indicated below each panel. Even though the ground truth spikes are used, the NMSE and leakage ratios indicated in the second and third columns are remarkably high. This further shows that the usage of conventional definitions and GPFA estimates is not optimal for the recovery of signal and noise correlations. In accordance with our theoretical analysis in Direct Extraction of Signal and Noise Correlations from Two-Photon Calcium Imaging of Ensemble Neuronal Activity, the performance of the two-stage Pearson method significantly improves as the number of trials is increased to $L=1000$, a number that is unrealistic in the context of typical two-photon imaging experiments. However, our proposed method shown in Figure 2 achieves comparable performance with number of trials as low as $L=20$. In summary, these results suggest that the two-stage methods produce highly biased estimates under limited number of trials, even if the ground truth spikes were ideally deconvolved from the two-photon data.

Figure 2—figure supplement 3.

Performance comparison under stimulus integration model mismatch.

Estimated noise and signal correlation matrices from different methods based on data generated with non-linear stimulus integration. Spikes were generated by replacing the linear receptive field model [Formula/expression omitted. See PDF] with a non-linear one given by [Formula/expression omitted. See PDF], but a linear stimulus model was used for estimation (i.e., [Formula/expression omitted. See PDF]). Rows from left to right: ground truth, proposed method, Pearson correlations from two-photon recordings, two-stage Pearson estimates and two-stage GPFA estimates. The normalized mean squared error (NMSE) of each estimate with respect to the ground truth and the leakage effect quantified by the ratio between out-of-network and in-network power (leakage) are indicated below each panel. While the NMSE in our proposed signal correlation estimates under this setting is greater than that in Figure 2 with no model mismatch, our proposed estimates still outperform existing methods. In addition, model mismatch in the stimulus integration component does not affect the accuracy of noise correlations estimated by our method.

Figure 2—figure supplement 4.

Performance under calcium decay model mismatch.

(A) Proposed noise and signal correlation estimates for data simulated at lower SNR than the setting of Figure 2 and model mismatch introduced by using a second-order autoregressive model for the calcium decay. The ground truth correlations are the same as those in Figure 2. The NMSE and leakage ratio are given at the bottom. (B) Putative spikes (green) and estimated calcium concentrations (purple). The model mismatch and lower SNR result in slight performance degradation compared to Figure 2 (in terms of NMSE and leakage), and our method is capable of recovering the underlying correlations faithfully.

Figure 2—figure supplement 5.

Performance comparison under varying SNR levels and firing rates.

Performance comparison with respect to varying SNR levels and average firing rates. (A) NMSE (top) and leakage ratios (bottom) for the noise (left) and signal (right) correlation estimates vs. SNR (in dB), for the proposed method, Pearson correlations from two-photon data and two-stage Pearson method. The SNR setting corresponding to Figure 2 is indicated by a dashed vertical line. The mean and standard deviation (std) of the normalized performance gain of the proposed method in comparison to the two existing methods are indicates as insets in each panel. (B) Same organization as panel A, but with respect to varying firing rates (in Hz). (C) Sample simulated white observation noise (red), two-photon observations (black, re-scaled for ease of visual comparison) and ground truth spikes (blue), as well as the estimated calcium concentrations (purple) and putative spikes (green) for the 1st trial of neuron 1. While the performance of all methods degrade at low SNR levels or firing rates (SNR lt₁₀ dB, firing rate lt_0.5 Hz), our proposed method outperforms the existing methods for almost all SNR and firing rate settings considered.

Figure 2—figure supplement 6.

Performance comparison under observation noise model mismatch.

Performance comparison with respect to varying SNR levels and average firing rates, with additional observation noise model mismatch. (A) NMSE (top) and leakage ratios (bottom) for the noise (left) and signal (right) correlation estimates vs. SNR (in dB), for the proposed method, Pearson correlations from two-photon data and two-stage Pearson method. The observation noise is generated by a white noise signal with an additive drift component from a low-frequency auto-regressive process. The mean and standard deviation (std) of the normalized performance gain of the proposed method in comparison to the two existing methods are indicates as insets in each panel. (B) Same organization as panel A, but with respect to varying firing rates (in Hz). (C) Sample simulated observation noise (red), two-photon observations (black, re-scaled for ease of visual comparison) and ground truth spikes (blue), as well as the estimated calcium concentrations (purple) and putative spikes (green) for the 1st trial of neuron 1. Panels (D), (E), and (F) are respectively in the same organization as panels (A), (B), and (C), but the observation noise is generated by a pink noise process. Our proposed method outperforms the existing methods for a wide range of SNR and firing rate values and under both observation noise model mismatch conditions.

The first column of Figure 2A shows the ground truth noise (top) and signal (bottom) correlations (diagonal elements are all equal to one and omitted for visual convenience). The second column shows estimates of the noise and signal correlations using our proposed method, which closely match the ground truth. The third, fourth and fifth columns, respectively, show the results of the Pearson correlations from the two-photon data, two-stage Pearson, and two-stage GPFA estimation methods. Through a qualitative visual inspection, it is evident that these methods incur high false alarms and mis-detections of the ground truth correlations.

To quantify these comparisons, the normalized mean square error (NMSE) of different estimates with respect to the ground truth are shown below each of the subplots (Figure 2A). Our proposed method achieves the lowest NMSE compared to the others. Furthermore, we observed a significant mixing between signal and noise correlations in these other estimates. To quantify this leakage effect, we first classified each of the correlation entries as in-network or out-of-network, based on being non-zero or zero in the ground truth, respectively (see Performance evaluation). We then computed the ratio between the power of out-of-network components and the power of in-network components as a measure of leakage. The leakage ratios are also reported in Figure 2A. The leakage of our proposed estimates is the lowest of all four techniques, in estimating both the signal and noise correlations. In order to further probe the performance of our proposed method, the simulated external stimulus s_t, latent trial-dependent process $x_{t,1}^{(1)}$ , simulated observations $y_{t,1}^{(1)}$ , estimated calcium concentration ${\widehat{z}}_{t,1}^{(1)}$ , the putative spikes ${\widehat{n}}_{t,1}^{(1)}:={\widehat{z}}_{t,1}^{(1)}-\alpha {\widehat{z}}_{t-1,1}^{(1)}$ , and the estimated mean of the latent state [Formula/expression omitted. See PDF] , for the first trial of the first neuron are shown in Figure 2B. These results demonstrate the ability of the proposed estimation framework in accurately identifying the latent processes, which in turn leads to an accurate estimation of the signal and noise correlations as shown in Figure 2B.

The main sources of the observed performance gap between our proposed method and the existing ones are the bias incurred by treating the fluorescence traces as spikes, low spiking rates, non-linearity of spike generation with respect to intrinsic and external covariates, and sensitivity to spike deconvolution. For the latter, we demonstrated the sensitivity of the two-stage Pearson estimates to the choice of the deconvolution technique in Figure 2—figure supplement 1. Furthermore, in order to isolate the effect of said non-linearities on the estimation performance, we applied the two-stage methods to ground truth spikes in Figure 2—figure supplement 2. Our analysis showed that both two-stage estimates incur significant estimation errors even if the spikes were recovered perfectly, mainly due to the limited number of trials ($L=20$ here). In accordance with our theoretical analysis of the asymptotic behavior of the conventional signal and noise correlation estimates given in Appendix 1, we also showed in Figure 2—figure supplement 2 that the performance of the two-stage Pearson estimates based on ground truth spikes, but using $L=1000$ trials, dramatically improves. Our proposed method, however, was capable of producing reliable estimates with the number of trials as low as $L=20$, which is typical in two-photon imaging experiments.

Analysis of robustness with respect to modeling assumptions

While the preceding results are quite favorable to our proposed method, the underlying generative models were the same as those used to estimate signal and noise correlations, which is in contrast to conventional real data validation with known ground truth. Access to ground truth correlations in two-photon imaging experimental settings, however, is quite challenging. In order to further probe the robustness of our proposed method in the absence of ground truth data, we utilized surrogate data that parallel the setting of Figure 2, but deviate from our modeling assumptions.

Robustness to stimulus integration model mismatch. First, we considered surrogate data generated with a non-linear stimulus integration model by replacing the linear receptive field component [Formula/expression omitted. See PDF] with [Formula/expression omitted. See PDF], where [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF] are akin to quadratic receptive field components. We assumed a linear stimulus integration model in our estimation framework (i.e., [Formula/expression omitted. See PDF]). Figure 2—figure supplement 3 shows the resulting correlation estimates. While the performance of our proposed signal correlation estimates degrade under this setting as compared to Figure 2, our proposed estimates still outperform existing methods. In addition, the model mismatch in the stimulus integration component does not affect the accuracy of noise correlation estimation in our method.

Robustness to calcium decay model mismatch. Next, we tested our proposed estimation framework on data simulated with a different calcium decay model. Specifically, we simulated data with second-order autoregressive calcium dynamics, and at a lower signal-to-noise ratio (SNR) compared to the setting of Figure 2, and used our inference framework which assumes first order calcium dynamics for estimation. Figure 2—figure supplement 4 shows the corresponding noise and signal correlations estimated by the proposed method under these conditions. Even though the performance slightly degrades (in terms of NMSE and leakage ratio), our method is able to recover the underlying correlations faithfully under this setting.

Robustness to SNR level and firing rate. Next, we compared the performance of Pearson and Two-Stage Pearson methods with our proposed method under varying SNR levels and average firing rates, as shown in Figure 2—figure supplement 5. While the performance of all methods degrades at low SNR levels or firing rates (SNR < 10 dB, firing rate < 0.5 Hz), our proposed method outperforms the existing methods for a wide range of SNR and firing rate values. To quantify this comparison, we have also indicated the mean and standard deviation of the relative performance gain of our proposed estimates across SNR levels and firing rates as insets in Figure 2—figure supplement 5.

Robustness to observation noise model mismatch. Finally, we repeated the foregoing comparisons under varying SNR levels and firing rates, only now we included an additional observation noise model mismatch. Similar to the treatment in Deneux et al., 2016, we considered two temporally correlated observation noise models: white noise with a low frequency drift (Figure 2—figure supplement 6, top panels) and pink noise (Figure 2—figure supplement 6, bottom panels). In accordance with the results in Figure 2—figure supplement 5, our proposed method outperforms the existing ones for a wide range of SNR and firing rate values and under both observation noise model mismatch conditions. From Figure 2—figure supplement 6C and F, it can be observed that the ground truth spikes are favorably recovered as a byproduct of our method, even though the estimated calcium concentrations are contaminated by the temporally correlated fluctuations in observation noise. This in turn results in accurate signal and noise correlation estimates.

Simulation study 2: spontaneous activity

Next, we present the results of a simulation study in the absence of external stimuli (i.e. ${\displaystyle {\mathbf{s}}_t=\mathbf{0}}$), pertaining to the spontaneous activity condition. It is noteworthy that the proposed method can readily be applied to estimate noise correlations during spontaneous activity, by simply setting the external stimulus [Formula/expression omitted. See PDF] and the receptive field [Formula/expression omitted. See PDF] to zero in the update rules (see Proposed forward model for details). We simulated the ensemble spiking activity based on a Poisson process (Smith and Brown, 2003) using a discrete time-rescaling procedure (Brown et al., 2002; Smith and Brown, 2003), so that the data are generated using a different model than that used in our inference framework (i.e. Bernoulli process with a logistic link as outlined in Proposed forward model). As such, we eliminated potential performance biases in favor of our proposed method by introducing the aforementioned model mismatch. We simulated $L=20$ independent trials of spontaneous activity of $N=30$ neurons, observed for a time duration of $T=5000$ time frames. The number of neurons in this study is notably larger than that used in the previous one, to examine the scalability of our proposed approach with respect to the ensemble size.

Figure 3 shows the comparison of the noise correlation matrices estimated by our proposed method, Pearson correlations from two-photon recordings, two-stage Pearson, and two-stage GPFA estimates, with respect to the ground truth. The Pearson and the two-stage estimates are highly variable and result in excessive false detections. Our proposed estimate, however, closely follows the ground truth, which is also reflected by the comparatively lower NMSE and leakage ratios, in spite of the mismatch between the models used for data generation and inference. In addition, our proposed method exhibits favorable scaling with respect to the ensemble size, thanks to the underlying low-complexity variational updates (see Low-complexity parameter updates for details).

Figure 3.

Results of simulation study 2.

Estimated noise correlation matrices using different methods based from spontaneous activity data. Rows from left to right: ground truth, proposed method, Pearson correlations from two-photon recordings, two-stage Pearson and two-stage GPFA estimates. The normalized mean squared error (NMSE) of each estimate with respect to the ground truth and the ratio between out-of-network power and in-network power (leakage) are shown below each panel.

Real data study 1: mouse auditory cortex under random tone presentation

We next applied our proposed method to experimentally recorded two-photon observations from the mouse primary auditory cortex (A1). The dataset consisted of recordings from 371 excitatory neurons in layer 2/3 A1, from which we selected $N=16$ responsive neurons (i.e. neurons that exhibited at least one spiking event in at least half of the trials considered; see Guidelines for model parameter settings). A random sequence of four tones was presented to the mouse, with the same sequence being repeated for $L=10$ trials. Each trial consisted of $T=3600$ time frames, and each tone was 2 s long followed by a 4 s silent period (see Experimental procedures for details). We considered an integration window of $R=25$ frames for stimulus encoding (see Guidelines for model parameter settings for details). The comparison of the noise and signal correlation estimates obtained by our proposed method, Pearson correlations from two-photon recordings, two-stage Pearson and two-stage GPFA methods is shown in Figure 4A. The spatial map of the 16 neurons considered in the analysis in the field of view is shown in Figure 4B. Figure 4C shows the stimulus tone sequence s_t , two-photon observations $y_{t,1}^{(1)}$ , estimated calcium concentration ${\widehat{z}}_{t,1}^{(1)}$ , putative spikes ${\widehat{n}}_{t,1}^{(1)}:={\widehat{z}}_{t,1}^{(1)}-\alpha {\widehat{z}}_{t-1,1}^{(1)}$ and the estimated mean of the latent state [Formula/expression omitted. See PDF] , for the first trial of the first neuron.

Figure 4.

Application to experimentally-recorded data from the mouse A1.

(A) Estimated noise (top) and signal (bottom) correlation matrices using different methods. Rows from left to right: proposed method, Pearson correlations from two-photon data, two-stage Pearson and two-stage GPFA estimates. (B) Location of the selected neurons with the highest activity in the field of view. (C) Presented tone sequence (orange), observations (black), estimated calcium concentrations (purple), putative spikes (green) and estimated mean latent state (blue) in the first trial of the first neuron. (D) Null distributions of chance occurrence of dissimilarities between signal and noise correlation estimates using different methods. The observed test statistic in each case is indicated by a dashed vertical line. (E) Scatter plots of signal vs. noise correlations for individual cell pairs (blue dots) corresponding to each method. Data were normalized for comparison by computing z-scores. For each case, the linear regression model fit is shown in red, and the slope and p-value of the t-test are indicated as insets.

Figure 4—figure supplement 1.

Probing the effect of stimulus integration window length on the performance of the proposed estimates.

Proposed noise and signal correlation estimates under different settings of the stimulus integration window length ($R$). (A) Proposed noise correlation (top) and signal correlation (bottom) estimates under different settings of $R$, from left to right: $R=1$, $R=10$, $R=25$ and $R=50$. (B) Null distributions of dissimilarities between proposed signal and noise correlation estimates corresponding to different choice of $R$. The observed test statistic in each case is indicated by a dashed vertical line, and the p-values are indicated above each panel. These results show that small values of $R=1$ and $R=10$ are not adequate to capture stimulus effect. However, both signal and noise correlation estimates exhibit consistency for $R=25$ and $R=50$.

Figure 4—figure supplement 2.

Inspecting the inferred latent processes under high fluorescence activity due to rapid increase in firing rate.

The fluorescence observations (black), inferred calcium concentrations (purple) and putative spikes (green) by our proposed method, for a sample data segment with high fluorescence activity due to successive closely spaced spikes. The rise onset of the fluorescence activity is marked by the vertical dashed line and spiking magnitude level of 1 is indicated by the horizontal dashed line. The proposed method favorably recovers the underlying calcium concentrations by predicting putative spikes in successive windows following the rapid rise of the fluorescence and with magnitudes possibly larger than 1.

We estimated the Best Frequency (BF) of each neuron as the tone that resulted in the highest level of fluorescence activity. The results in Figure 4A are organized such that the neurons with the same BF are neighboring, with the BF increasing along the diagonal. Thus, expectedly (Bowen et al., 2020) our proposed method as well as the Pearson and two-stage Pearson estimates show high signal correlations along the diagonal. However, the two-stage GPFA estimates do not reveal such a structure. By visual inspection, as also observed in the simulation studies, the Pearson correlations from two-photon recordings, two-stage Pearson and two-stage GPFA estimates have significant leakage between the signal and noise correlations, whereas our proposed signal and noise correlation estimates in Figure 4A suggest distinct spatial structures.

To quantify this visual comparison, we used a statistic based on the Tanimoto similarity metric (Lipkus, 1999), denoted by [Formula/expression omitted. See PDF] for two matrices [Character omitted. See PDF] and [Character omitted. See PDF]. As a measure of dissimilarity, we used [Formula/expression omitted. See PDF] (see Performance evaluation for details). The comparison of [Formula/expression omitted. See PDF] for the four estimates is presented in the second column of Table 1. To assess statistical significance, for each comparison we obtained null distributions corresponding to chance occurrence of dissimilarities using a shuffling procedure as shown in Figure 4D, and then computed one-tailed $p$-values from those distributions (see Performance evaluation for details). Table 1 and Figure 4D includes these p-values, which show that the proposed estimates (boldface numbers in Table 1, second column) indeed have the highest dissimilarity between signal and noise correlations. The higher leakage effect in the other three estimates is also reflected in their smaller [Formula/expression omitted. See PDF] values.

Table 1.

Dissimilarity metric statistics for the estimates in Figure 4A (also illustrated in Figure 4D), linear regression statistics of the comparison between signal and noise correlations in Figure 4E, and the average NMSE across 50 trials used in the shuffling procedure illustrated in Figure 5A.

	Dissimilarity [Formula/expression omitted. See PDF]	Regression statistics (Figure 4E)		Shuffling test (Figure 5)
Estimate	Figure 4D	Slope (p-value)	${\mathrm{R}}^2$ Value	NMSE in [Formula/expression omitted. See PDF]	NMSE in [Formula/expression omitted. See PDF]
Proposed	$\mathbf{0.8725}(p<{10}^{-4})$	${\displaystyle 0.02(p=0.84)}$	$4\times {10}^{-4}$	$\mathbf{1.07}\pm \mathbf{0.16}$	$\mathbf{1.32}\pm \mathbf{0.19}$
Pearson	$0.6675(p=0.71)$	$0.33(p=2\times {10}^{-4})$	0.11	0	0
Two-stage Pearson	$0.7325(p=0.09)$	$0.15(p=0.10)$	0.02	$1.84\pm 0.34$	$0.55\pm 0.12$
Two-stage GPFA	$0.7625(p<{10}^{-4})$	$0.02(p=0.86)$	$3\times {10}^{-4}$	$2.32\pm 0.52$	$2.26\pm 0.51$

To further investigate this effect, we have depicted the scatter plots of signal vs. noise correlations estimated by each method in Figure 4E. To examine the possibility of the leakage effect on a pairwise basis, we performed linear regression in each case. The slope of the model fit, the p-value for the corresponding t-test, and the $R^2$ values are reported in the third and fourth columns of Table 1 (the slope and p-values are also shown as insets in Figure 4E). Consistent with the results of Winkowski and Kanold, 2013, the Pearson estimates suggest a significant correlation between the signal and noise correlation pairs (as indicated by the higher slope in Figure 4E). However, none of the other estimates (including the proposed estimates) in Figure 4E register a significant trend between signal and noise correlations. This further corroborates our assessment of the high leakage between signal and noise correlations in Pearson estimates, since such a leakage effect could result in overestimation of the trend between the signal and noise correlation pairs. The signal and noise correlations estimated by our proposed method show no pairwise trend, suggesting distinct patterns of stimulus-dependent and stimulus-independent functional connectivity (Kohn et al., 2016; Montijn et al., 2014; Rothschild et al., 2010; Keeley et al., 2020).

A key advantage of our proposed method over the Pearson and two-stage approaches is the explicit modeling of stimulus integration. The relevant parameter in this regard is the length of the stimulus integration window $R$. While in our simulation studies the value of $R$ was known, it needs to be set by the user in real data applications. To this end, domain knowledge or data-driven methods such as cross-validation and model order selection can be utilized (see Guidelines for model parameter settings for details). Noting that the number of parameters to be estimated linearly scales with $R$, it must be chosen large enough to capture the stimulus effects, yet small enough to result in favorable computational complexity. Here, given that the typical tone response duration of mouse A1 neurons is ${\displaystyle <1}$ s (Linden et al., 2003; DeWeese et al., 2003; Petrus et al., 2014), with a sampling frequency of $f_s=30$ Hz, we surmised that a choice of $R\sim 30$ suffices to capture the stimulus effects. We further examined the effect of varying $R$ on the proposed correlation estimates in Figure 4—figure supplement 1. As shown, small values of $R$ (e.g. $R=1$ or 10) may not be adequate to fully capture stimulus integration effects. By considering values of $R$ in the range $25-50$, we observed that the correlation estimates remain stable. We thus chose $R=25$ for our analysis.

Careful inspection of the second panel in Figure 4C shows that the fluorescence activity often saturates to ∼4 times its baseline value. This effect is due to successive closely spaced spikes, which implies the occurrence of more than one spike per frame and thus violates our Bernoulli modeling assumption. To inspect the performance of our method more carefully under this scenario, we show in Figure 4—figure supplement 2 a zoomed-in view of the estimated latent processes ${\widehat{z}}_{t,1}^{(1)}$ (calcium concentration) and ${\widehat{n}}_{t,1}^{(1)}$ (putative spikes) for a sample data segment with high fluorescence activity. The estimated latent processes reveal two mechanisms leveraged by our inference method to mitigate the aforementioned model mismatch: first, our proposed method predicts spiking events in adjacent time frames to compensate for rapid increase in firing rate and thus infers calcium concentration levels that match the observed fluorescence; secondly, even though our generative model assumes that there is only one spiking event in a given time frame, this restriction is mitigated in our inference framework by relaxing the constraint ${\widehat{n}}_{t,l}^{(j)}:={\widehat{z}}_{t,l}^{(j)}-\alpha {\widehat{z}}_{t-1,l}^{(j)}\le 1$, as explained in Low-complexity parameter updates. While this relaxation was performed for the sake of tractability of the inverse solution, it in fact leads to improved estimation results under episodes of rapid increase in firing rate, by allowing the putative spike magnitudes ${\widehat{n}}_{t,l}^{(j)}$ to be greater than 1. The latter is evident in the magnitude of the inferred spikes in Figure 4—figure supplement 2, following the rise of fluorescence activity.

Given that the ground truth correlations are not available for a direct comparison, we instead performed a test of specificity that reveals another key limitation of existing methods. Fluorescence observations exhibit structured dynamics due to the exponential intracellular calcium concentration decay (as shown in Figure 4C, for example), which are in turn related to the underlying spikes that are driven non-linearly by intrinsic/extrinsic stimuli as well as the properties of the indicator used. As such, an accurate inference method is expected to be specific to this temporal structure. To test this, we randomly shuffled the $T$ time frames consistently in the same order in all trials, in order to fully break the temporal structure governing calcium decay dynamics, and then estimated correlations from these shuffled data using the different methods. The resulting estimates of noise correlations are shown in Figure 5A for one instance of such shuffled data. The average NMSE for a total of 50 shuffled samples with respect to the original un-shuffled estimates (in Figure 4A) are tabulated in the fifth and sixth columns of Table 1, and are also indicated below each panel in Figure 5A.

Figure 5.

Assessing the specificity of different estimation results shown in Figure 4.

Rows from left to right: proposed method, Pearson correlations from two-photon data, two-stage Pearson and two-stage GPFA estimates. (A) The estimated noise correlations using different methods after random temporal shuffling of the observations. The mean and standard deviation of the NMSE across 50 trials are indicated below each panel. (B) Histograms of the noise correlation estimates between the first and third neurons over the 50 temporal shuffling trials. The estimate based on the original (un-shuffled) data in each case is indicated by a dashed vertical line.

A visual inspection of Figure 5A shows that the Pearson correlations from two-photon recordings expectedly remain unchanged. Since this method treats each time frame to be independent, temporal shuffling does not impact the correlations in anyway. On the other extreme, both of the two-stage estimates seem to detect highly variable and large correlation values, despite operating on data that lacks any relevant temporal structure. Our proposed method, however, remarkably produces negligible correlation estimates. Although both the two-stage and proposed estimates show variability with respect to the shuffled data (Table 1, fifth column), the standard deviation of the NMSE values of our proposed method are considerably smaller than those of the two-stage methods (Table 1, fifth column). For further inspection, the histograms of a single element ([Formula/expression omitted. See PDF]) of the estimated correlation matrices across the 50 shuffling trials are shown in Figure 5B. The original un-shuffled estimates are marked by the dashed vertical lines in each case. The proposed estimate in Figure 5B is highly concentrated around zero, even though the un-shuffled estimate is non-zero. However, the two-stage estimates produce correlations that are widely variable across the shuffling trials. This analysis demonstrates that our proposed method is highly specific to the temporal structure of fluorescence observations, whereas the Pearson correlations from two-photon recordings, two-stage Pearson and two-stage GPFA methods fail to be specific.

Real data study 2: spontaneous vs. stimulus-driven activity in the mouse A1

To further validate the utility of our proposed methodology, we applied it to another experimentally-recorded dataset from the mouse A1 layer 2/3. This experiment pertained to trials of presenting a sequence of short white noise stimuli, randomly interleaved with silent trials of the same duration. Figure 6A shows a sample trial sequence. The two-photon recordings thus contained episodes of stimulus-driven and spontaneous activity (see Experimental procedures for details). Under this experimental setup, it is expected that the noise correlations are invariant across the spontaneous and stimulus-driven conditions. In accordance with the foregoing results of real data study 1, we also expect the signal and noise correlation patterns to be distinct. Each trial considered in the analysis consisted of $T=765$ frames (see Experimental procedures for details). We selected $N=10$ responsive neurons (according to the criterion described in Guidelines for model parameter settings), each with $L=10$ trials. Similar to real data study 1, we chose a stimulus integration window of length $R=25$ frames.

Figure 6.

Comparison of spontaneous and stimulus-driven activity in the mouse A1.

(A) A sample trial sequence in the experiment. Stimulus-driven (stim) trials were recorded with randomly interleaved spontaneous (spon) trials of the same duration. (B) Estimated noise and signal correlation matrices under spontaneous (top) and stimulus-driven (bottom) conditions. Rows from left to right: proposed method, Pearson correlations from two-photon data, two-stage Pearson and two-stage GPFA estimates. (C) Location of the selected neurons with highest activity in the field of view. (D) Stimulus onsets (orange), observations (black), estimated calcium concentrations (purple) and putative spikes (green) for the first trial from two pairs of neurons with high signal correlation (top) and high noise correlation (bottom), as identified by the proposed estimates.

Figure 6—figure supplement 1.

Histograms of the similarity/dissimilarity metrics under the shuffling procedure.

Null distributions of (A) the similarities between [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF] (top: [Formula/expression omitted. See PDF]) and (B) the dissimilarities between [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF] (bottom: [Formula/expression omitted. See PDF]), obtained by the shuffling procedure applied to the results of real data study 2 in Figure 6. The observed test statistic in each case is indicated by a dashed vertical line. Rows from left to right: proposed method, Pearson correlations from two-photon data, two-stage Pearson correlations and two-stage GPFA estimates. These results show that the only statistically significant outcomes (with $p\le 0.05$) are the similarities and dissimilarities obtained by our proposed method.

Figure 6B shows the resulting noise and signal correlation estimates under the spontaneous ([Formula/expression omitted. See PDF], top) and stimulus-driven ([Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF], bottom) conditions. Figure 6C shows the spatial map of the 10 neurons considered in the analysis in the field of view. A visual inspection of the first column of Figure 6B indeed suggests that [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF] are saliently similar, and distinct from [Formula/expression omitted. See PDF]. The Pearson correlations obtained from two-photon data (second column) and the two-stage Pearson and GPFA estimates (third and fourth columns, respectively), however, evidently lack this structure. As in the previous study, we quantified this visual comparison using the similarity metric [Formula/expression omitted. See PDF] and the dissimilarity metric [Formula/expression omitted. See PDF] (see Performance evaluation for details). These statistics are reported in Table 2 along with the p-values (null distributions are shown in Figure 6—figure supplement 1), which show that the only significant outcomes (boldface numbers) are those of our proposed method. While it is expected from the experiment design for the noise correlations under the two settings to be similar, the only method that detects this expected outcome with statistical significance is our proposed method. Moreover, the statistically significant dissimilarity between the signal and noise correlations of our proposed estimates corroborate the hypothesis that signal and noise are encoded by distinct functional networks (Kohn et al., 2016; Montijn et al., 2014; Rothschild et al., 2010; Keeley et al., 2020).

Table 2.

Similarity/dissimilarity metric statistics for the estimates in Figure 6.

Estimation method	[Formula/expression omitted. See PDF]	[Formula/expression omitted. See PDF]
Proposed	0.5716 ($p=0.003$)	0.7946 ($p=0.004$)
Pearson	0.3031 ($p=0.61$)	0.5032 ($p=0.92$)
Two-stage Pearson	0.2790 ($p=0.05$)	0.7862 ($p=0.39$)
Two-stage GPFA	0.2008 ($p=0.50$)	0.7792 ($p=0.22$)

Furthermore, Figure 6D shows the time course of the stimulus, observations, estimated calcium concentrations and putative spikes for the first trial from two pairs of neurons with high signal correlation ($j=2,8$, top) and high noise correlation ($j=3,5$, bottom). As expected, the putative spiking activity of the neurons with high signal correlation (top) are closely time-locked to the stimulus onsets. The activity of the two neurons with high noise correlation (bottom), however, is not time-locked to the stimulus onsets, even though the two neurons exhibit highly correlated activity. The correlations estimated via the proposed method thus encode substantial information about the inter-dependencies of the spiking activity of the neuronal ensemble.

Real data study 3: spatial analysis of signal and noise correlations in the mouse A1

Lastly, we applied our proposed method to examine the spatial distribution of signal and noise correlations in the mouse A1 layers 2/3 and 4 (data from Bowen et al., 2020). The dataset included fluorescence activity recorded during multiple experiments of presenting sinusoidal amplitude-modulated tones, with each stimulus being repeated across several trials (see Experimental procedures and Bowen et al., 2020 for experimental details). In each experiment, we selected on average around 20 responsive neurons for subsequent analysis (according to the criterion described in Guidelines for model parameter settings). For brevity, we compare the estimates of signal and noise correlations using our proposed method only with those obtained by Pearson correlations from the two-photon data. The latter method was also used in previous analyses of data from this experimental paradigm (Winkowski and Kanold, 2013).

In parallel to the results reported in Winkowski and Kanold, 2013, Figure 7A and Figure 7B illustrate the correlation between the signal and noise correlations in layers 2/3 and 4, respectively. Consistent with the results of Winkowski and Kanold, 2013, the signal and noise correlations exhibit positive correlation in both layers, regardless of the method used. However, the correlation coefficients (i.e. slopes in the insets) identified by our proposed method are notably smaller than those obtained from Pearson correlations, in both layer 2/3 (Figure 7A) and layer 4 (Figure 7B). Comparing this result with our simulation studies suggests that the stronger linear trend between the signal and noise correlations observed using the Pearson correlation estimates is likely due to the mixing between the estimates of signal and noise correlations. As such, our method suggests that the signal and noise correlations may not be as highly correlated with one another as indicated in previous studies of layer 2/3 and 4 in mouse A1 (Winkowski and Kanold, 2013).

Next, to evaluate the spatial distribution of signal and noise correlations, we plotted the correlation values for pairs of neurons as a function of their distance for layer 2/3 (Figure 7C) and layer 4 (Figure 7D). The distances were discretized using bins of length ${\displaystyle 10\mu m}$. The scatter of the correlations along with their median at each bin are shown in all panels. Then, to examine the spatial trend of the correlations, we performed linear regression in each case. The slope of the model fit, the p-value for the corresponding t-test, and the $R^2$ values are reported in Table 3 (the slope and p-values are also shown as insets in Figure 7C and D).

Table 3.

Linear regression statistics for the analysis of correlations vs. cell-pair distance.

	Statistics of layer 2/3 correlations		Statistics of layer 4 correlations
Correlations	Slope (p-value)	${\mathrm{R}}^2$ Value	Slope (p-value)	${\mathrm{R}}^2$ Value
Proposed Signal Corr.	[Formula/expression omitted. See PDF] ($p=0.002$)	0.012	[Formula/expression omitted. See PDF] ($p=3\times {10}^{-6}$)	0.023
Pearson Signal Corr.	$-5\times {10}^{-5}$ ($p=0.02$)	0.007	$-3\times {10}^{-5}$ ($p=0.02$)	0.005
Proposed Noise Corr.	[Formula/expression omitted. See PDF] ($p=0.005$)	0.010	[Formula/expression omitted. See PDF] ($p=0.06$)	0.004
Pearson Noise Corr.	$-4\times {10}^{-5}$ ($p=0.1$)	0.003	$-5\times {10}^{-5}$ ($p=0.02$)	0.005

From Table 3 and Figure 7C and D (upper panels), it is evident that the signal correlations show a significant negative trend with respect to distance, using both methods and in both layers. However, the slope of these negative trends identified by our method (boldface numbers in Table 3) is notably steeper than those identified by Pearson correlations. On the other hand, the trends of the noise correlations with distance (bottom panels) are different between our proposed method and Pearson correlations: our proposed method shows a significant negative trend in layer 2/3, but not in layer 4, whereas the Pearson correlations of the two-photon data suggest a significant negative trend in layer 4, but not in layer 2/3. In addition, the slopes of these negative trends identified by our method (boldface numbers in Table 3) are steeper than or equal to those identified by Pearson correlations.

Our proposed estimates also indicate that noise correlations are sparser and less widespread in layer 4 (Figure 7D) than in layer 2/3 (Figure 7C). To further investigate this observation, we depicted the two-dimensional spatial spread of signal and noise correlations in both layers and for both methods in Figure 7E and F, by centering each neuron at the origin and overlaying the individual spatial spreads. The horizontal and vertical axes in each panel represent the relative dorsoventral and rostrocaudal distances, respectively, and the heat-maps represent the magnitude of correlations. Comparing the proposed noise correlation spread in Figure 7E with the corresponding spread in Figure 7F, we observe that the noise correlations in layer 2/3 are indeed more widespread and abundant than in layer 4, as can be expected by more extensive intralaminar connections in layer 2/3 vs. 4 (Watkins et al., 2014; Meng et al., 2017a; Meng et al., 2017b; Kratz and Manis, 2015).

Figure 7.

Comparison of signal and noise correlations across layers 2/3 and 4.

(A) Scatter-plot of noise vs. signal correlations (blue) for individual cell-pairs in layer 2/3, based on the proposed (left) and Pearson estimates (right). Data were normalized for comparison by computing z-scores. The linear model fits are shown in red, and the slope and p-value of the t-tests are indicated as insets. Panel (B) corresponds to layer 4 in the same organization as panel A. (C) Signal (top) and noise (bottom) correlations vs. cell-pair distance in layer 2/3, based on the proposed (left) and Pearson estimates (right). Distances were binned to ${\displaystyle 10\mu m}$ intervals. The median of the distributions (black) and the linear model fit (red) are shown in each panel. The slope of the linear model fit, and the p-value of the t-test are also indicated as insets. Dashed horizontal lines indicate the zero-slope line for ease of visual comparison. Panel D corresponds to layer 4 in the same organization as panel C. (E) Spatial spread of signal (top) and noise (bottom) correlations in layer 2/3, based on the proposed (left) and Pearson estimates (right). The horizontal and vertical axes in each panel respectively represent the relative dorsoventral and rostrocaudal distances between each cell-pair, and the heat-map indicates the magnitude of correlations. Marginal distributions of the signal (blue) and noise (red) correlations along the dorsoventral and rostrocaudal axes for the proposed method (darker colors) and Pearson method (lighter colors) are shown at the top and right sides of the sub-panels. Panel F corresponds to layer 4 in the same organization as panel E.

Figure 7—figure supplement 1.

Comparing the marginal distributions of signal and noise correlations along the dorsoventral and rostrocaudal axes.

Comparison of marginal distributions of signal and noise correlations. (A) Cumulative marginal probability distributions of signal (blue) and noise (red) correlations along the rostrocaudal (top) and dorsoventral (bottom) directions, as estimated by the proposed method (left) and Pearson correlations from two-photon data (right), in layer 2/3 neurons. The Kolmogorov–Smirnov (KS) test statistic along with the corresponding p-values are indicated as insets in each panel. Panel B shows the results for layer 4 in the same organization as panel A. These results show that along both directions and in both layers, the signal correlation distributions are significantly different from the corresponding noise correlation distributions, consistently for both methods. However, the KS statistics (i.e. effect sizes) for the proposed estimate are remarkably larger than those obtained from the Pearson estimates.

Figure 7—figure supplement 2.

Marginal angular distributions of signal and noise correlations.

Polar plots of the angular marginal distributions of correlations. (A) Polar histograms indicating the distribution of signal (top) and noise (bottom) correlations as a function of relative angle (in the dorsoventral-rostrocaudal coordinate system) between pairs of neurons in layer 2/3, as estimated by the proposed method (left) and Pearson correlations from two-photon data (right). The KS test statistic comparing each polar distribution with a uniform distribution (shown in magenta), along with the corresponding p-values are indicated below each polar plot. The mode of each probability distribution is also indicated in blue fonts. Panel B shows the results for layer four in the same organization as panel A. All distributions are significantly non-uniform, and particularly indicate a rostrocaudal directionality in layer 4 (as indicated by the mode angles in panel B).

The spatial spreads of signal and noise correlations based on the Pearson estimates are remarkably similar in both layers (Figure 7E and F, right panels), whereas they are saliently different for our proposed estimates (Figure 7E and F, left panels). This further corroborates our hypothesis on the possibility of high mixing between the signal and noise correlation estimates obtained by the Pearson correlation of two-photon data. To further examine the differences between the signal and noise correlations, the marginal distributions along the dorsoventral and rostrocaudal axes are shown in Figure 7E and F, selectively overlaid for ease of visual comparison. To quantify the differences between the spatial distributions of signal and noise correlations estimated by each method, we performed Kolmogorov-Smirnov (KS) tests on each pair of marginal distributions, which are summarized in Figure 7—figure supplement 1. Although the marginal distributions of signal and noise correlations are significantly different in all cases from both methods, the effect sizes of their difference (KS statistics) are higher for our proposed estimates compared to those of the Pearson estimates.

Finally, the spatial spreads of correlations for either method and in each layer suggest non-uniform angular distributions with possibly directional bias. To test this effect, we computed the angular marginal distributions and performed KS tests for non-uniformity, which are reported in Figure 7—figure supplement 2. These tests indicate that all distributions are significantly non-uniform. In addition, the angular distributions of both signal and noise correlations in layer 4 exhibit salient modes in the rostrocaudal direction, whereas they are less directionally selective in layer 2/3 (Figure 7—figure supplement 2).

In summary, the spatial trends identified by our proposed method are consistent with empirical observations of spatially heterogeneous pure-tone frequency tuning by individual neurons in auditory cortex (Winkowski and Kanold, 2013). The improved correspondence of our proposed method compared to results obtained using Pearson correlations could be the result of the demixing of signal and noise correlations in our method. As a result of the demixing, our proposed method also suggests that noise correlations have a negative trend with distance in layer 2/3, but are much sparser and spatially flat in layer 4. In addition, the spatial spread patterns of signal and noise correlations are more structured and remarkably more distinct for our proposed method than those obtained by the Pearson estimates.

Theoretical analysis of the bias and variance of the proposed estimators

Finally, we present a theoretical analysis of the bias and variance of the proposed estimator. Note that our proposed estimation method has been developed as a scalable alternative to the intractable maximum likelihood (ML) estimation of the signal and noise covariances (see Overview of the proposed estimation method). In order to benchmark our estimates, we thus need to evaluate the quality of said ML estimates. To this end, we derived bounds on the bias and variance of the ML estimators of the kernel [Formula/expression omitted. See PDF] for $j=1,\mathrm{\cdots },N$ and the noise covariance [Formula/expression omitted. See PDF]. In order to simplify the treatment, we posit the following mild assumptions:

Assumption (1). We assume a scalar time-varying external stimulus (i.e. [Formula/expression omitted. See PDF], and hence [Formula/expression omitted. See PDF] ). Furthermore, we set the observation noise covariance to be [Formula/expression omitted. See PDF], for notational convenience.

Assumption (2). We derive the performance bounds in the regime where $T$ and $L$ are large, and thus do not impose any prior distribution on the correlations, which are otherwise needed to mitigate overfitting (see Preliminary assumptions).

Assumption (3). We assume the latent trial-dependent process and stimulus to be slowly varying signals, and thus adopt a piece-wise constant model in which these processes are constant within consecutive windows of length $W$ (i.e. [Formula/expression omitted. See PDF] and $s_t=s_{W_k}$, for $(k-1)W+1\le t<kW$ and $k=1,\mathrm{\cdots },K$ with $W_k=(k-1)W+1$ and $KW=T$) for our theoretical analysis, as is usually done in spike count calculations for conventional noise correlation estimates.

Our main theoretical result is as follows:

Theorem 1 (Performance Bounds) Let $q\mathrm{>}\frac{\mathrm{1}}{\mathrm{64}}$, $\mathrm{0}\mathrm{<}ϵ\mathrm{<}\mathrm{1}\mathrm{/}\mathrm{2}$, and $\mathrm{0}\mathrm{<}\eta \mathrm{\le }\mathrm{1}\mathrm{/}\mathrm{2}$ be fixed constants, ${\sigma }_m^{\mathrm{2}}\mathrm{:=}{\max }_i{\mathrm{(}{\mathrm{\Sigma }}_x\mathrm{)}}_{i\mathrm{,}i}$ and ${\sigma }_s^{\mathrm{2}}\mathrm{:=}\frac{\mathrm{1}}{K}\mathit{}{\mathrm{\sum }}_{k\mathrm{=}\mathrm{1}}^K{s}_{W_k}^{\mathrm{2}}$. Then, under Assumptions (1 - 3), the bias and variance of the maximum likelihood estimators $\widehat{\mathrm{d}}$ and ${\widehat{\mathrm{\Sigma }}}_x$, conditioned on an event ${\mathrm{A}}_W$ with $\mathrm{P}\mathit{}\mathrm{\left(}{\mathrm{A}}_W\mathrm{\right)}\mathrm{\ge }\mathrm{1}\mathrm{-}\eta$ satisfy:

[Formula omitted. See PDF]

for all $i\mathrm{,}j\mathrm{=}\mathrm{1}\mathrm{,}\mathrm{2}\mathrm{,}\mathrm{\cdots }\mathrm{,}N$, if${\displaystyle \log W\ge max{\displaystyle \left\{\frac{\log (8KLN/\eta )}{q},\frac{32{\sigma }_m^{2}q}{{ϵ}^2},\frac{2\log (64q)}{1-2ϵ},\frac{max\{6.25,4{\left(\Vert {\mathit{\mu }}_x{\Vert }_{\mathrm{\infty }}+\underset{k,j}{max}\left\{|s_{W_k}{d}_j|\right\}\right)}^2\}}{8q{\sigma }_m^2},\log 2\right\}},}$where ${\tau }_j$ and ${\stackrel{\mathrm{~}}{\tau }}_j$ denote bounded terms that are $\mathrm{O}\mathit{}\mathrm{(}{\sigma }_w^{\mathrm{2}}\mathrm{)}$ or $\mathrm{O}\mathit{}\mathrm{\left(}\frac{\mathrm{1}}{W}\mathrm{\right)}$, ${\xi }_{i\mathrm{,}j}$ and ${\stackrel{\mathrm{~}}{\xi }}_{i\mathrm{,}j}$ denote bounded terms that are $\mathrm{O}\mathit{}\mathrm{(}{\sigma }_w^{\mathrm{2}}\mathrm{)}$ or $\mathrm{O}\mathit{}\mathrm{\left(}\frac{\mathrm{1}}{W^{\mathrm{1}\mathrm{-}\mathrm{2}\mathit{}ϵ}}\mathrm{\right)}$ and $C_{\mathrm{1}}\mathrm{,}{C}_{\mathrm{2}}\mathrm{,}{C}_{\mathrm{3}}$ and $C_{\mathrm{4}}$ are bounded constants given in Appendix 2.

Proof. The proof of Theorem 1 is provided in Appendix 2.

■

In order to discuss the implications of this theoretical result, several remarks are in order:

Remark 1: Achieving near Oracle performance

A common benchmark in estimation theory is the performance of the idealistic oracle estimator, in which an oracle directly observes the true latent process [Formula/expression omitted. See PDF] and the true kernel d_j and forms the correlation estimates. In this case, the oracle would incur zero bias and variance of order [Formula/expression omitted. See PDF] in estimating d_j, and outputs an estimate of [Formula/expression omitted. See PDF] with bias and variance in the order of [Formula/expression omitted. See PDF]. Theorem 1 indeed states that for sufficiently large $W$ and small ${\sigma }_w$, the bias and variance of the ML estimators are arbitrarily close to those of the oracle estimator. Recall that our variational inference framework is in fact a solution technique for the regularized ML problem. Hence, the bounds in Theorem 1 provide a benchmark for the expected performance of the proposed estimators, by quantifying the excess bias and variance over the performance of the oracle estimator.

Remark 2: Effect of the observation noise and observation duration

As the assumed window of stationarity $W\to \mathrm{\infty }$ (and hence the observation duration $T\to \mathrm{\infty }$), the loss of performance of the proposed estimators only depends on ${\sigma }_w^2$, the variance of the observation noise. As a result, at a given observation noise variance ${\sigma }_w^2$, these bounds provide a sufficient upper bound on the time duration of the observations required for attaining a desired level of estimation accuracy. It is noteworthy that ${\sigma }_w^2$ is typically small in practice, as it pertains to the effective observation noise and is significantly diminished by pixel averaging of the fluorescence traces following cell segmentation.

Remark 3: Effect of the number of trials

Finally, note that the bounds in Theorem 1 have terms that also drop as the number of trials $L$ grows. These terms in fact pertain to the performance of the oracle estimator. As the number of trials grows ($L\to \mathrm{\infty }$), the oracle estimates become arbitrarily close to the true parameters [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF]. Thus, our theoretical performance bounds also provide a sufficient upper bound on the number of trials $L$ required for the oracle estimator to attain a desired level of estimation accuracy.

Discussion

We developed a novel approach for the joint estimation of signal and noise correlations of neuronal activities directly from two-photon calcium imaging observations and tested our method with experimental data. Existing widely used methods either take the fluorescence traces as surrogates of spiking activity, or first recover the unobserved spikes using deconvolution techniques, both followed by computing Pearson correlations or connectivity matrices. As such, they typically result in estimates that are highly biased and are heavily dependent on the choice of the spike deconvolution technique. We addressed these issues by explicitly relating the signal and noise covariances to the observed two-photon data via a multi-tier Bayesian model that accounts for the observation process and non-linearities involved in spiking activity. We developed an efficient estimation framework by integrating techniques from variational inference and state-space estimation. We also established performance bounds on the bias and variance of the proposed estimators, which revealed favorable scaling with respect to the observation noise and trial length.

We demonstrated the utility of our proposed estimation framework on both simulated and experimentally recorded data from the mouse auditory cortex. In our simulation studies, we evaluated the robustness of our proposed method with respect to several model mismatch conditions induced by the stimulus integration model, calcium decay, SNR level, firing rate, and temporally correlated observation noise. In all cases, we observed that our proposed estimates outperform the existing methods in recovering the signal and noise correlations.

There are two main sources for the observed performance gap between our proposed method and existing approaches. The first source is the favorable soft decisions on the timing of spikes achieved by our method as a byproduct of the iterative variational inference procedure. An accurate probabilistic decoding of spikes results in better estimates of the signal and noise correlations, and conversely having more accurate estimates of the signal and noise covariances improves the probabilistic characterization of spiking events. This is in contrast with both the Pearson correlations computed from two-photon data and two-stage methods: in computing the Pearson correlations from two-photon data, spike timing is heavily blurred by the calcium decay; in the two-stage methods, erroneous hard decisions on the timing of spikes result in biases that propagate to and contaminate the downstream signal and noise correlation estimation and thus results in significant errors.

The second source of performance improvement is the explicit modeling of the non-linear mapping from stimulus and latent covariates to spiking through a canonical point process model, which is in turn tied to a two-photon observation model in a multi-tier Bayesian fashion. Our theoretical analysis in Theorem 1 corroborates that this virtue of our proposed methodology results in robust performance under limited number of trials. As we have shown in Appendix 1, as the number of trials $L$ and trial duration $T$ tend to infinity, conventional notions of signal and noise correlation indeed recover the ground truth signal and noise correlations, as the biases induced by non-linearities average out across trial repetitions. However, as exemplified in Figure 2—figure supplement 2, in order to achieve comparable performance to our method using few trials (e.g. $L=20$), the conventional correlation estimates require considerably more trials (e.g. $L=1000$).

Application to two-photon data recorded from the mouse primary auditory cortex showed that unlike the aforementioned existing methods, our estimates provide noise correlation structures that are expectedly invariant across spontaneous and stimulus-driven conditions within an experiment, while producing signal correlation structures that are largely distinct from those given by noise correlation. These results provide evidence for the involvement of distinct functional neuronal network structures in encoding the stimulus-dependent and stimulus-independent information.

Our analysis of the relationship between the signal and noise correlations in layers 2/3 and 4 in mouse A1 indicated a smaller correlation between signal and noise correlations than previously reported (Winkowski and Kanold, 2013). Thus, our proposed method suggests that the signal and noise correlations reflect distinct circuit mechanisms of sound processing in layers 2/3 vs 4. The spatial distribution of signal correlations obtained by our method was consistent with previous work showing significant negative trends with distance (Winkowski and Kanold, 2013). However, in addition, our proposed method revealed a significant negative trend of noise correlations with distance in layer 2/3, but not in layer 4, in contrast to the outcome of Pearson correlation analysis. The lack of a negative trend in layer 4 could be attributed to the sparse nature of the noise correlation spread in layer 4, as revealed by our analysis of two-dimensional spatial spreads. The latter analysis indeed revealed that the noise correlations in layer 2/3 are more widespread than those in layer 4, consistent with existing work based on whole-cell patch recordings (Meng et al., 2017a; Meng et al., 2017b).

The two-dimensional spatial spreads of signal and noise correlations obtained by our method are more distinct than those obtained by Pearson correlations. The spatial spreads also allude to directionality of the functional connectivity patterns, with a notable rostrocaudal preference in layer 4. This result seems surprising in light of existing evidence for quasi-rostrocaudal organization of the tonotopic axis in mouse A1 (Romero et al., 2020). However, given the heterogeneity of tuning in both layers 2/3 and 4 with a best frequency interqartile range of ∼1–1.5 octaves over the imaging field (Bowen et al., 2020) and using supra-threshold tones, we expect that the tones will drive not only neurons with the corresponding best frequency, but also neurons tuned to neighboring frequencies. Moreover, there is high connectivity between layer 4 cells within a few 100 μm across the tonotopic axis (Kratz and Manis, 2015; Meng et al., 2017a), potentially amplifying and broadening the effect of supra-threshold tones.

Our proposed method can scale up favorably to larger populations of neurons, thanks to the underlying low-complexity variational updates in the inference procedure. Due to its minimal dependence on training data, our estimation framework is also applicable to single-session analysis of two-photon data with limited number of trials and duration. Another useful byproduct of the proposed framework is gaining access to approximate posterior densities in closed-form, which allows further statistical analyses such as construction of confidence intervals. Our proposed methodology can thus be used as a robust and scalable alternative to existing approaches for extracting neuronal correlations from two-photon calcium imaging data.

A potential limitation of our proposed generative model is the assumption that there is at most one spiking event per time frame for each neuron, in light of the fact that typical two-photon imaging frame durations are in the range of 30–100 ms. Average spike rates of excitatory neurons in mouse A1 layers 2/3 and 4 are of the order of ${\displaystyle <10}$ Hz (Petrus et al., 2014; Forli et al., 2018) and thus our model is reasonable for the current study, although it might not be optimal during bursting activity. It is noteworthy that we relax this assumption in the inference framework by allowing the magnitude of putative spikes to be greater than one, thus alleviating the model mismatch during episodes of rapid increase in firing rate. This assumption can also be made more precise by adopting a Poisson model, but that would render closed-form variational density updates intractable.

Furthermore, in the regime of extremely low spiking rate and high observation noise, the proposed method may fail to capture the underlying correlations faithfully and its performance degrades to those of existing methods based on Pearson correlations, as we have shown through our simulation studies. Nevertheless, our method addresses key limitations of conventional signal and noise correlation estimators that persist even in high spiking rate and high SNR conditions.

Our proposed estimation framework can be used as groundwork for incorporating other notions of correlation such as the connected correlation function (Martin et al., 2020), and to account for non-Gaussian and higher order structures arising from spatiotemporal interactions (Kadirvelu et al., 2017; Yu et al., 2011). Other possible extensions of this work include leveraging variational inference beyond the mean-field regime (Wang and Blei, 2013), extension to time-varying correlations that underlie rapid task-dependent dynamics, and extension to non-linear models such as those parameterized by neural networks (Aitchison et al., 2017). In the spirit of easing reproducibility, a MATLAB implementation of our proposed method as well as the data used in this work are made publicly available (Rupasinghe, 2020; Rupasinghe et al., 2021).

Materials and methods

Proposed forward model

Suppose we observe fluorescence traces of $N$ neurons, for a total duration of $T$ discrete-time frames, corresponding to $L$ independent trials of repeated stimulus. Let ${\displaystyle {\mathbf{y}}_{t,l}:=[y_{t,l}^{(1)},y_{t,l}^{(2)},\cdots ,y_{t,l}^{(N)}{]}^{\mathrm{\top }}}$, ${\displaystyle {\mathbf{z}}_{t,l}:=[z_{t,l}^{(1)},z_{t,l}^{(2)},\cdots ,z_{t,l}^{(N)}{]}^{\mathrm{\top }}}$, and ${\displaystyle {\mathbf{n}}_{t,l}:=[n_{t,l}^{(1)},n_{t,l}^{(2)},\cdots ,n_{t,l}^{(N)}{]}^{\mathrm{\top }}}$ be the vectors of noisy observations, intracellular calcium concentrations, and ensemble spiking activities, respectively, at trial $l$ and frame $t$. We capture the dynamics of [Formula/expression omitted. See PDF] by the following state-space model:

[Formula omitted. See PDF]

where [Formula/expression omitted. See PDF] represents the scaling of the observations, [Formula/expression omitted. See PDF] is zero-mean i.i.d. Gaussian noise with covariance [Formula/expression omitted. See PDF], and ${\displaystyle 0\le \alpha <1}$ is the state transition parameter capturing the calcium dynamics through a first order model. Note that this state-space is non-Gaussian due to the binary nature of the spiking activity, that is, ${\displaystyle n_{t,l}^{(j)}\in \{0,1\}}$. We model the spiking data as a point process or Generalized Linear Model with Bernoulli statistics (Eden et al., 2004; Paninski, 2004; Smith and Brown, 2003; Truccolo et al., 2005):${\displaystyle n_{t,l}^{(j)}\sim \mathrm{Bernoulli}\left({\lambda }_{t,l}^{(j)}\right),{\lambda }_{t,l}^{(j)}=\varphi \left(x_{t,l}^{(j)},{{\mathbf{d}}_j}^{\mathrm{\top }}{\mathbf{s}}_t\right),}$where ${\lambda }_{t,l}^{(j)}$ is the conditional intensity function (Truccolo et al., 2005), which we model as a non-linear function of the known external stimulus [Formula/expression omitted. See PDF] and the other latent intrinsic and extrinsic trial-dependent covariates, ${\displaystyle {\mathbf{x}}_{t,l}:={\displaystyle [x_{t,l}^{(1)},x_{t,l}^{(2)},\cdots ,x_{t,l}^{(N)}{]}^{\mathrm{\top }}}}$. While we assume the stimulus [Formula/expression omitted. See PDF] to be common to all neurons, we model the distinct effect of this stimulus on the [Formula/expression omitted. See PDF] neuron via an unknown kernel [Formula/expression omitted. See PDF], akin to the receptive field.

The non-linear mapping of our choice is the logistic link, which is also the canonical link for a Bernoulli process in the point process and Generalized Linear Model frameworks (Truccolo et al., 2005). Thus, we assume:

[Formula omitted. See PDF]

Finally, we assume the latent trial dependent covariates to be a Gaussian process [Formula/expression omitted. See PDF], with mean ${\displaystyle {\mathit{\mu }}_x:={\displaystyle [{\mu }_x^{(1)},{\mu }_x^{(2)},\cdots ,{\mu }_x^{(N)}{]}^{\mathrm{\top }}}}$ and covariance [Formula/expression omitted. See PDF].

The probabilistic graphical model in Figure 8 summarizes the main components of the aforementioned forward model. According to this forward model, the underlying noise covariance matrix that captures trial-to-trial variability can be identified as [Formula/expression omitted. See PDF]. The signal covariance matrix, representing the covariance of the neural activity arising from the repeated application of the stimulus [Formula/expression omitted. See PDF], is given by [Formula/expression omitted. See PDF], where ${\displaystyle \mathbf{D}:=[{\displaystyle {\mathbf{d}}_1,{\mathbf{d}}_2,\cdots ,{\mathbf{d}}_N}]\in {\mathbb{R}}^{M\times N}}$. The signal and noise correlation matrices, denoted by [Character omitted. See PDF] and [Character omitted. See PDF], can then be obtained by standard normalization of [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF]:

[Formula omitted. See PDF]

The main problem is thus to estimate [Formula/expression omitted. See PDF] from the noisy and temporally blurred data [Formula/expression omitted. See PDF] .

Figure 8.

Probabilistic graphical model of the proposed forward model.

The fluorescence observations at the [Formula/expression omitted. See PDF] time frame and [Formula/expression omitted. See PDF] trial: [Formula/expression omitted. See PDF], are noisy surrogates of the intracellular calcium concentrations: [Formula/expression omitted. See PDF]. The calcium concentration at time $t$ is a function of the spiking activity [Formula/expression omitted. See PDF], and the calcium activity at the previous time point [Formula/expression omitted. See PDF]. The spiking activity is driven by two independent mechanisms: latent trial-dependent covariates [Formula/expression omitted. See PDF], and contributions from the known external stimulus [Formula/expression omitted. See PDF], which we model by [Formula/expression omitted. See PDF] (in which the receptive field [Character omitted. See PDF] is unknown). Then, we model [Formula/expression omitted. See PDF] as a Gaussian process with constant mean [Formula/expression omitted. See PDF], and unknown covariance [Formula/expression omitted. See PDF]. Finally, we assume the covariance [Formula/expression omitted. See PDF] to have an inverse Wishart prior distribution with hyper-parameters ${\displaystyle {\mathit{\psi }}_x}$ and ${\rho }_x$. Based on this forward model, the inverse problem amounts to recovering the signal and noise correlations by directly estimating [Formula/expression omitted. See PDF] and [Character omitted. See PDF] (top layer) from the fluorescence observations [Formula/expression omitted. See PDF] (bottom layer).

Overview of the proposed estimation method

First, given a limited number of trials $L$ from an ensemble with typically low spiking rates, we need to incorporate suitable prior assumptions to avoid overfitting. Thus, we impose a prior [Formula/expression omitted. See PDF] on the noise covariance, to compensate sparsity of data. A natural estimation method to estimate [Formula/expression omitted. See PDF] in a Bayesian framework is to maximize the observed data likelihood ${\displaystyle p({\displaystyle \{{\mathbf{y}}_{t,l}{\}}_{t,l=1}^{T,L}|{\mathbf{\Sigma }}_x,\mathbf{D}})}$, that is maximum likelihood (ML). Thus, we consider the joint likelihood of the observed data and latent processes to perform Maximum a Posteriori (MAP) estimation:

(4)

${\displaystyle \begin{array}{ll}{\displaystyle p(\mathbf{y},\mathbf{z},\mathbf{x},{\mathbf{\Sigma }}_x|\mathbf{D})}& ={\displaystyle p_{\mathsf{p}\mathsf{r}}({\mathbf{\Sigma }}_x)\prod _{t,l=1}^{T,L}\frac{1}{\sqrt{(2\pi {)}^N|{\mathbf{\Sigma }}_w|}}\exp \left(-\frac{1}{2}({\mathbf{y}}_{t,l}-\mathbf{A}{\mathbf{z}}_{t,l}{)}^{\mathrm{\top }}{\mathbf{\Sigma }}_w^{-1}({\mathbf{y}}_{t,l}-\mathbf{A}{\mathbf{z}}_{t,l})\right)}\\ & \times {\displaystyle \prod _{t,l,j=1}^{T,L,N}\frac{{\left(\exp \left(x_{t,l}^{(j)}+{{\mathbf{d}}_j}^{\mathrm{\top }}{\mathbf{s}}_t\right)\right)}^{z_{t,l}^{(j)}-\alpha z_{t-1,l}^{(j)}}}{1+\exp \left(x_{t,l}^{(j)}+{{\mathbf{d}}_j}^{\mathrm{\top }}{\mathbf{s}}_t\right)}\prod _{t,l=1}^{T,L}\frac{1}{\sqrt{(2\pi {)}^N|{\mathbf{\Sigma }}_x|}}\exp \left(-\frac{1}{2}({\mathbf{x}}_{t,l}-{\mathit{\mu }}_x{)}^{\mathrm{\top }}{\mathbf{\Sigma }}_x^{-1}({\mathbf{x}}_{t,l}-{\mathit{\mu }}_x)\right).}\end{array}}$

Inspecting this MAP problem soon reveals that estimating [Formula/expression omitted. See PDF] and [Character omitted. See PDF] is a challenging task: (1) standard approaches such as Expectation-Maximization (EM) (Shumway and Stoffer, 1982) are intractable due to the complexity of the model, arising from the hierarchy of latent processes and the non-linearities involved in their mappings and (2) the temporal coupling of the likelihood in the calcium concentrations makes any potential direct solver scale poorly with $T$.

Thus, we propose an alternative solution based on Variational Inference (VI) (Beal, 2003; Blei et al., 2017; Jordan et al., 1999). VI is a method widely used in Bayesian statistics to approximate unwieldy posterior densities using optimization techniques, as a low-complexity alternative strategy to Markov Chain Monte Carlo sampling (Hastings, 1970) or empirical Bayes techniques such as EM. To this end, we treat [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF] as latent variables and [Formula/expression omitted. See PDF] and [Character omitted. See PDF] as unknown parameters to be estimated. We introduce a framework to update the latent variables and parameters sequentially, with straightforward update rules. We will describe the main ingredients of the proposed framework in the following subsections. Hereafter, we use the shorthand notations [Formula/expression omitted. See PDF], [Formula/expression omitted. See PDF], and [Formula/expression omitted. See PDF].

Preliminary assumptions

For the sake of simplicity, we assume that the constants α, [Character omitted. See PDF], [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF] are either known or can be consistently estimated from pilot trials. Next, we take [Formula/expression omitted. See PDF] to be an Inverse Wishart density:

[Formula omitted. See PDF]

which turns out to be the conjugate prior in our model. Thus, [Formula/expression omitted. See PDF] and ${\rho }_x$ will be the hyper-parameters of our model. Procedures for hyper-parameter tuning and choosing the key model parameters are given in subsections Hyper-parameter tuning and Guidelines for model parameter settings, respectively.

Decoupling via Pólya-Gamma augmentation

Direct application of VI to problems containing both discrete and continuous random variables results in intractable densities. Specifically, finding a variational distribution for [Formula/expression omitted. See PDF] in our model with a standard distribution is not straightforward, due to the complicated posterior arising from co-dependent Bernoulli and Gaussian random variables. In order to overcome this difficulty, we employ Pólya-Gamma (PG) latent variables (Pillow and Scott, 2012; Polson et al., 2013; Linderman et al., 2016). We observe from Equation 4 that the posterior density, [Formula/expression omitted. See PDF] is conditionally independent in $t,l$ with: 

[Formula omitted. See PDF]

Thus, upon careful inspection, we see that this density has the desired form for the PG augmentation scheme (Polson et al., 2013). Accordingly, we introduce a set of auxiliary PG-distributed i.i.d. latent random variables ${\displaystyle {\mathit{\omega }}_{t,l}:=[{\omega }_{t,l}^{(1)},{\omega }_{t,l}^{(2)},\cdots ,{\omega }_{t,l}^{(N)}{]}^{\mathrm{\top }}}$, ${\displaystyle {\omega }_{t,l}^{(j)}\sim \mathrm{PG}(1,0)}$ for $1\le j\le N$, $1\le t\le T$ and $1\le l\le L$, to derive the complete data log-likelihood:

(5)

${\displaystyle \begin{array}{ll}& {\displaystyle \log p(\mathbf{y},\mathbf{z},\mathbf{x},\mathit{\omega },{\mathbf{\Sigma }}_x|\mathbf{D})}\\ & ={\displaystyle -\frac{TL}{2}\log |{\mathbf{\Sigma }}_x|+\log p_{\mathsf{p}\mathsf{r}}\left({\mathbf{\Sigma }}_x\right)+\sum _{t,l=1}^{T,L}\{-\frac{1}{2}{\left({\mathbf{y}}_{t,l}-\mathbf{A}{\mathbf{z}}_{t,l}\right)}^{\mathrm{\top }}{\mathbf{\Sigma }}_w^{-1}\left({\mathbf{y}}_{t,l}-\mathbf{A}{\mathbf{z}}_{t,l}\right)-\frac{1}{2}({\mathbf{x}}_{t,l}-{\mathit{\mu }}_x{)}^{\mathrm{\top }}{\mathbf{\Sigma }}_x^{-1}({\mathbf{x}}_{t,l}-{\mathit{\mu }}_x)}\\ & {\displaystyle +\sum _{j=1}^N\{(z_{t,l}^{(j)}-\alpha z_{t-1,l}^{(j)}-\frac{1}{2})\left(x_{t,l}^{(j)}+{{\mathbf{d}}_j}^{\mathrm{\top }}{\mathbf{s}}_t\right)-\frac{1}{2}{\omega }_{t,l}^{(j)}{\left(x_{t,l}^{(j)}+{{\mathbf{d}}_j}^{\mathrm{\top }}{\mathbf{s}}_t\right)}^2+\log p_{{}_{\mathsf{P}\mathsf{G}\mathsf{(}\mathsf{1}\mathsf{,}\mathsf{0}\mathsf{)}}}\left({\omega }_{t,l}^{(j)}\right)\}\}+C,}\end{array}}$where [Formula/expression omitted. See PDF] and $C$ accounts for terms not depending on [Formula/expression omitted. See PDF], [Formula/expression omitted. See PDF] and [Character omitted. See PDF]. The complete data log-likelihood is notably quadratic in [Formula/expression omitted. See PDF], which as we show later admits efficient estimation procedures with favorable scaling in $T$.

Deriving the optimal variational densities

In this section, we will outline the procedure of applying VI to the latent variables [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF], assuming that the parameter estimates [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF] of the previous iteration are available. The methods that we propose to update the parameters [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF] subsequently, will be discussed in the next section.

The objective of variational inference is to posit a family of approximate densities [Character omitted. See PDF] over the latent variables, and to find the member of that family that minimizes the Kullback-Leibler (KL) divergence to the exact posterior:

[Formula omitted. See PDF]

However, evaluating the KL divergence is intractable, and it has been shown (Blei et al., 2017) that an equivalent result to this minimization can be obtained by maximizing the alternative objective function, called the evidence lower bound (ELBO):${\displaystyle \mathrm{ELBO}(q)=\mathbb{E}[\log p(\mathbf{x},\mathit{\omega },{\mathbf{\Sigma }}_x,\mathbf{y}|\hat{\mathbf{z}},\hat{\mathbf{D}})]-\mathbb{E}[\log q(\mathbf{x},\mathit{\omega },{\mathbf{\Sigma }}_x|\hat{\mathbf{z}},\hat{\mathbf{D}})].}$

Further, we assume [Character omitted. See PDF] to be a mean-field variational family (Blei et al., 2017), resulting in the overall variational density of the form:

[Formula omitted. See PDF]

Under the mean field assumptions, the maximization of the ELBO can be derived using the optimization algorithm ‘Coordinate Ascent Variational Inference’ (CAVI) (Bishop, 2006; Blei et al., 2017). Accordingly, we see that the optimal variational densities in Equation 6 take the forms:${\displaystyle \begin{array}{ll}{\displaystyle \log q^{\ast }\left({\mathbf{x}}_{t,l}\right)}& \propto {\mathbb{E}}_{q^{\ast }\left({\mathbf{\Sigma }}_x\right)q^{\ast }\left({\mathit{\omega }}_{t,l}\right)}\left[\log p\left({\mathbf{x}}_{t,l}|{\mathit{\omega }}_{t,l},{\mathbf{\Sigma }}_x,\mathbf{y},\hat{\mathbf{z}},\hat{\mathbf{D}}\right)\right],\\ \log q^{\ast }\left({\omega }_{t,l}^{(j)}\right)& \propto {\mathbb{E}}_{q^{\ast }\left({\mathbf{x}}_{t,l}\right)}\left[\log p\left({\omega }_{t,l}^{(j)}|{\mathbf{x}}_{t,l},{\mathbf{\Sigma }}_x,\mathbf{y},\hat{\mathbf{z}},\hat{\mathbf{D}}\right)\right],\\ \log q^{\ast }\left({\mathbf{\Sigma }}_x\right)& \propto {\mathbb{E}}_{q^{\ast }\left(\mathbf{x}\right)}\left[\log p\left({\mathbf{\Sigma }}_x|\mathbf{x},\mathbf{y},\hat{\mathbf{z}},\hat{\mathbf{D}}\right)\right].\end{array}}$

Upon evaluation of these expectations, we derive the optimal variational distributions as:

[Formula omitted. See PDF]

whose parameters ${\displaystyle {\mathbf{m}}_{{\mathbf{x}}_{t,l}}:={{\displaystyle [m_{{\mathbf{x}}_{t,l}}^{(1)},m_{{\mathbf{x}}_{t,l}}^{(2)},\cdots ,m_{{\mathbf{x}}_{t,l}}^{(N)}]}}^T}$, [Formula/expression omitted. See PDF], $c_{t,l}^{(j)}$, [Formula/expression omitted. See PDF], and ${\gamma }_x$ can be updated given parameter estimates [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF]:${\displaystyle \begin{array}{l}{\displaystyle {\displaystyle {\mathbf{Q}}_{{\mathbf{x}}_{t,l}}=({\tilde{\mathbf{\Omega }}}_{t,l}+{\gamma }_x{\mathbf{P}}_x^{-1}{)}^{-1}},{\displaystyle {\mathbf{m}}_{{\mathbf{x}}_{t,l}}={\mathbf{Q}}_{{\mathbf{x}}_{t,l}}({\hat{\mathbf{z}}}_{t,l}-\alpha {\hat{\mathbf{z}}}_{t-1,l}-\frac{1}{2}\mathbf{1}-{\tilde{\mathbf{\Omega }}}_{t,l}{\hat{\mathbf{D}}}^{\mathrm{\top }}{\mathbf{s}}_t+{\gamma }_x{\mathbf{P}}_x^{-1}{\mathit{\mu }}_x)},}\\ {\displaystyle {\mathbf{P}}_x:={\mathit{\psi }}_x+\sum _{t,l=1}^{T,L}\left\{{\mathbf{Q}}_{{\mathbf{x}}_{t,l}}+{\mathbf{m}}_{{\mathbf{x}}_{t,l}}{\mathbf{m}}_{{\mathbf{x}}_{t,l}}^{\mathrm{\top }}-{\mathit{\mu }}_x{\mathbf{m}}_{{\mathbf{x}}_{t,l}}^{\mathrm{\top }}-{\mathbf{m}}_{{\mathbf{x}}_{t,l}}{\mathit{\mu }}_x^{\mathrm{\top }}+{\mathit{\mu }}_x{\mathit{\mu }}_x^{\mathrm{\top }}\right\}},c_{t,l}^{(j)}={\displaystyle \sqrt{{\left({\mathbf{Q}}_{{\mathbf{x}}_{t,l}}\right)}_{j,j}+{\left(m_{{\mathbf{x}}_{t,l}}^{(j)}+{{\hat{\mathbf{d}}}_j}^{\mathrm{\top }}{\mathbf{s}}_t\right)}^2}},\end{array}}$and ${\gamma }_x:={\rho }_x+TL$, with ${\displaystyle {\tilde{\mathbf{\Omega }}}_{t,l}\in {\mathbb{R}}^{N\times N}}$ denoting a diagonal matrix with entries ${\displaystyle ({\tilde{\mathbf{\Omega }}}_{t,l}{)}_{j,j}:={\displaystyle \frac{1}{2c_{t,l}^{(j)}}\tanh (\frac{c_{t,l}^{(j)}}{2})}}$ and [Formula/expression omitted. See PDF] denoting the vector of all ones.

Low-complexity parameter updates

Note that even though [Character omitted. See PDF] is composed of the latent processes [Formula/expression omitted. See PDF], we do not use VI for its inference, and instead consider it as an unknown parameter. This choice is due to the temporal dependencies arising from the underlying state-space model in Equation 4, which hinders a proper assignment of variational densities under the mean field assumption. We thus seek to estimate both [Character omitted. See PDF] and [Character omitted. See PDF] using the updated variational density [Formula/expression omitted. See PDF].

First, note that the log-likelihood in Equation 5 is decoupled in $l$, which admits independent updates to [Formula/expression omitted. See PDF], for $l=1,\mathrm{\cdots },L$. As such, given an estimate [Formula/expression omitted. See PDF], we propose to estimate [Formula/expression omitted. See PDF] as:

${\displaystyle \begin{array}{ll}{\displaystyle \{{\hat{\mathbf{z}}}_{t,l}{\}}_{t=1}^T}& =\underset{\{{\mathbf{z}}_{t,l}{\}}_{t=1}^T}{\text{argmax}}{\mathbb{E}}_{q^{\ast }(\mathbf{x},\mathit{\omega },{\mathbf{\Sigma }}_x)}\left[\log p\left(\mathbf{y},\mathbf{z},\mathbf{x},\mathit{\omega },{\mathbf{\Sigma }}_x|\hat{\mathbf{D}}\right)\right]\\ & =\underset{\{{\mathbf{z}}_{t,l}{\}}_{t=1}^T}{\text{argmin}}\sum _{t=1}^T\left\{\frac{1}{2}{\left({\mathbf{y}}_{t,l}-\mathbf{A}{\mathbf{z}}_{t,l}\right)}^{\mathrm{\top }}{\mathbf{\Sigma }}_w^{-1}\left({\mathbf{y}}_{t,l}-\mathbf{A}{\mathbf{z}}_{t,l}\right)-\sum _{j=1}^N\left(m_{{\mathbf{x}}_{t,l}}^{(j)}+{{\hat{\mathbf{d}}}_j}^{\mathrm{\top }}{\mathbf{s}}_t\right)\left(z_{t,l}^{(j)}-\alpha z_{t-1,l}^{(j)}\right)\right\},\end{array}}$ under the constraints ${\displaystyle 0\le {\displaystyle z_{t,l}^{(j)}-\alpha z_{t-1,l}^{(j)}}\le 1}$, for $t=1,\mathrm{\cdots },T$ and $j=1,\mathrm{\cdots },N$. These constraints are a direct consequence of $n_{t,l}^{(j)}=z_{t,l}^{(j)}-\alpha z_{t-1,l}^{(j)}$ being a Bernoulli random variable with [Formula/expression omitted. See PDF]. While this problem is a quadratic program and can be solved using standard techniques, it is not readily decoupled in $t$, and thus standard solvers would not scale favorably in $T$.

Instead, we consider an alternative solution that admits a low-complexity recursive solution by relaxing the constraints. To this end, we relax the constraint [Formula/expression omitted. See PDF] and replace the constraint [Formula/expression omitted. See PDF] by penalty terms proportional to $\left|z_{t,l}^{(j)}-\alpha z_{t-1,l}^{(j)}\right|$. The resulting relaxed problem is thus given by:

[Formula omitted. See PDF]

where [Formula/expression omitted. See PDF] with $\beta \ge 1$ being a hyper-parameter. Given that the typical spiking rates are quite low in practice, [Formula/expression omitted. See PDF] is expected to be a negative number. Thus, we have assumed that [Formula/expression omitted. See PDF].

The problem of Equation 7 pertains to compressible state-space estimation, for which fast recursive solvers are available (Kazemipour et al., 2018). The solver utilizes the Iteratively Re-weighted Least Squares (IRLS) (Ba et al., 2014) framework to transform the absolute value in the second term of the cost function into a quadratic form in [Formula/expression omitted. See PDF], followed by Fixed Interval Smoothing (FIS) (Rauch et al., 1965) to find the minimizer. At iteration $k$, given a current estimate [Formula/expression omitted. See PDF], the problem reduces to a Gaussian state-space estimation of the form:

[Formula omitted. See PDF]

with [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF], where [Formula/expression omitted. See PDF] is a diagonal matrix with [Formula/expression omitted. See PDF], for some small constant $\epsilon >0$. This problem can be efficiently solved using FIS, and the iterations proceed for a total of $K$ times or until a standard convergence criterion is met (Kazemipour et al., 2018). It is noteworthy that our proposed estimator of the calcium concentration [Formula/expression omitted. See PDF] can be thought of as soft spike deconvolution, which naturally arises from our variational framework, as opposed to the hard spike deconvolution step used in two-stage estimators.

Finally, given [Formula/expression omitted. See PDF] and the updated [Formula/expression omitted. See PDF], the estimate of [Formula/expression omitted. See PDF] for $j=1,2,\mathrm{\cdots },N$ can be updated in closed-form by maximizing the expected complete log-likelihood [Formula/expression omitted. See PDF]:${\displaystyle {\displaystyle {\hat{\mathbf{d}}}_j={\displaystyle {\left(\sum _{t,l=1}^{T,L}\left(({\tilde{\mathbf{\Omega }}}_{t,l}{)}_{j,j}{\mathbf{s}}_t{{\mathbf{s}}_t}^{\mathrm{\top }}\right)\right)}^{-1}\left(\sum _{t,l=1}^{T,L}\{\left({\hat{z}}_{t,l}^{(j)}-\alpha {\hat{z}}_{t-1,l}^{(j)}-\frac{1}{2}\right){\mathbf{s}}_t-({\tilde{\mathbf{\Omega }}}_{t,l}{)}_{j,j}{m}_{{\mathbf{x}}_{t,l}}^{(j)}{\mathbf{s}}_t\}\right)}.}}$

The VI procedure iterates between updating the variational densities and parameters until convergence, upon which we estimate the noise and signal covariances as:

[Formula omitted. See PDF]

The overall combined iterative procedure is outlined in Algorithm 1. Furthermore, a MATLAB implementation of this algorithm is publicly available in Rupasinghe, 2020. It is worth noting that a special case of our proposed variational inference procedure can be used to estimate signal and noise correlations from electrophysiology recordings. Given that spiking activity, that is [Formula/expression omitted. See PDF], is directly observed in this case, the solution to the optimization problem in Equation (7) is no longer required. Thus, the parameters [Formula/expression omitted. See PDF] and [Character omitted. See PDF] can be estimated using a simplified variational procedure, which is outlined in Algorithm 2 in Appendix 3.

Guidelines for model parameter settings

There are several key model parameters that need to be set by the user prior to the application of our proposed method. Here, we provide our rationale and criteria for choosing these parameters, which could also serve as guidelines in facilitating the applicability and adoption of our method by future users. We will also provide the specific choices of these parameters used in our simulation studies and real data analyses.

Number of neurons selected for the analysis ($N$)

While our proposed method scales-up well with the population size due to low-complexity update rules involved, including neurons with negligible spiking activity in the analysis would only increase the complexity and potentially contaminate the correlation estimates. Thus, we performed an initial pre-processing step to extract $N$ neurons that exhibited at least one spiking event in at least half of the trials considered.

Stimulus integration window length ($R$)

The number of lags $R$ considered in stimulus integration is a key parameter that can be set through data-driven approaches or using prior domain knowledge. Examples of common data-driven criteria include cross-validation, Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC), which balance the estimation accuracy and model complexity (Arlot and Celisse, 2010; Ding et al., 2018).

To quantify the effect of $R$ on model complexity, we first describe the stimulus encoding model in our framework. Suppose that the onset of the [Formula/expression omitted. See PDF] tone in the stimulus set ($p=1,\mathrm{\cdots },P$, where $P$ is the number of distinct tones) is given by a binary sequence ${\displaystyle f_t^{(p)}\in \{0,1\}}$. The choice of $R$ implies that the response at time $t$ post-stimulus depends only on the $R$ most recent time lags. As such, the effective stimulus at time $t$ corresponding to tone $p$ is given by [Formula/expression omitted. See PDF]. By including all the $P$ tones, the overall effective stimulus at the $t^{th}$ time frame is given by [Formula/expression omitted. See PDF]. The stimulus modulation vector [Formula/expression omitted. See PDF] would thus be $RP$-dimensional. As a result, the number of parameters ($M=RP$) to be estimated linearly increases with $R$. By using additional domain knowledge, we chose $R$ to be large enough to capture the stimulus effects, and at the same time to be small enough to control the complexity of the algorithm.

As an example, given that the typical tone response duration of mouse primary auditory neurons is ${\displaystyle <1}$ s (Linden et al., 2003; DeWeese et al., 2003; Petrus et al., 2014), with a sampling frequency of $f_s=30$ Hz, a choice of $R\sim 30$ would suffice to capture the stimulus effects. By further examining the effect of varying $R$ on the proposed correlation estimates in Figure 4—figure supplement 1, we chose $R=25$ for our real data analyses.

Observation noise covariance ([Formula/expression omitted. See PDF]) and scaling matrix ([Character omitted. See PDF])

We assumed that the observation noise covariance [Formula/expression omitted. See PDF] is diagonal, and estimated the diagonal elements using the background fluorescence in the absence of spiking events, for each neuron. We set [Formula/expression omitted. See PDF], where [Formula/expression omitted. See PDF] represents the identity matrix, and estimated $a$ by considering the average increase in fluorescence after the occurrence of isolated spiking events. Specifically, we derived the average fluorescence activity of multiple trials triggered to the fluorescence rise onset, and set $a$ as the increment in the magnitude of this average fluorescence immediately following the rise onset.

State transition parameter (α)

We chose α in the range $[0.95,0.98]$, which match the slow dynamics of the calcium indicator in our data. We tested the robustness of our estimates under different choices of α in this range through the method outlined in Hyper-parameter tuning, and accordingly chose the optimal value of α.

Mean of the latent trial-dependent process ([Formula/expression omitted. See PDF])

We estimated [Formula/expression omitted. See PDF] as a constant that is proportional to the average firing rate. To this end, we parametrized each component of [Formula/expression omitted. See PDF] as ${\mu }_x^{(j)}=-a_{\mu }+b_{\mu }\frac{1}{LT}{\sum }_{t,l=1}^{T,L}{y}_{t,l}^{(j)}$, for $j=1,\mathrm{\cdots },N$. The constants $a_{\mu }$ and $b_{\mu }$ were chosen such that $-2\le {\mu }_x^{(j)}\le -10$, which gives the range of baseline parameters compatible with observed firing rates in our experimental data.

Parameter choices for simulation study 1

In the first simulation study, we set $\alpha =0.98$, $\beta =8$, [Formula/expression omitted. See PDF], [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF] ([Formula/expression omitted. See PDF] represents the identity matrix and [Formula/expression omitted. See PDF] represents the vector of all ones), so that the SNR of simulated data was in the same range as that of experimentally-recorded data. We used a [Formula/expression omitted. See PDF] order autoregressive process with a mean of -1 as the stimulus (s_t), and considered $R=2$ ($M=2$) lags of the stimulus (i.e. [Formula/expression omitted. See PDF]) in both the generative model and inference procedure. The components of the linear and quadratic stimulus modulation vectors, that is [Formula/expression omitted. See PDF], [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF], were chosen at random uniformly in the range $[-0.5,0.5]$. The variance of s_t was set in each case such that the average power of the overall signal component ([Formula/expression omitted. See PDF] for the linear model, and [Formula/expression omitted. See PDF] for the non-linear model) was comparable to the average power of the noise component ($x_{t,l}^{(j)}$).

Algorithm 1 Estimation of ${\displaystyle {\mathbf{\Sigma }}_x}$ and ${\displaystyle \mathbf{D}}$ through the proposed iterative procedure

Inputs: Ensemble of fluorescence measurements [Formula/expression omitted. See PDF], constants [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF], hyper-parameters [Formula/expression omitted. See PDF], ${\rho }_x$, β and $ϵ$, tolerance at convergence δ and the external stimulus [Formula/expression omitted. See PDF] Outputs: [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF] Initialization: Initial choice of [Formula/expression omitted. See PDF], ${\displaystyle {\tilde{\mathbf{\Omega }}}_t}$, [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF], $\mathrm{\text{residual}}\mathrm{=}\mathrm{10}\mathrm{}\delta$, ${\gamma }_x\mathrm{=}{\rho }_x\mathrm{+}L\mathrm{}T$ 1: while $\text{residual}\ge \delta$ do Estimate calcium concentrations using Fixed Interval Smoothing 2:    for ${\displaystyle l=1,\cdots ,L}$ do Forward filter: 3:      for ${\displaystyle t=1,...,T}$ do 4:              ${\displaystyle {\mathbf{z}}_{(t|t-1),l}=\alpha {\mathbf{z}}_{(t-1|t-1),l}}$ 5:              ${\displaystyle {\mathbf{P}}_{(t|t-1),l}={\alpha }^2{\mathbf{P}}_{(t-1|t-1),l}+{\mathbf{\Sigma }}_{{\mathbf{v}}_{t,l}}}$ 6:              ${\displaystyle {\mathbf{B}}_{t,l}={\mathbf{P}}_{(t|t-1),l}{\mathbf{A}}^{\mathrm{\top }}(\mathbf{A}{\mathbf{P}}_{(t|t-1),l}{\mathbf{A}}^{\mathrm{\top }}+{\mathbf{\Sigma }}_w{)}^{-1}}$ 7:              ${\displaystyle {\mathbf{z}}_{(t|t),l}={\mathbf{z}}_{(t|t-1),l}+{\mathbf{B}}_{t,l}(\mathbf{y}{}_{t,l}-\mathbf{A}{\mathbf{z}}_{(t|t-1),l})}$ 8:              ${\displaystyle {\mathbf{P}}_{(t|t),l}=(\mathbf{I}-{\mathbf{B}}_{t,l}\mathbf{A}){\mathbf{P}}_{(t|t-1),l}}$ 9:      end for Backward smoother: 10:      for ${\displaystyle t=T-1,...,1}$ do 11:              ${\displaystyle {\hat{\mathbf{z}}}_{t,l}={\mathbf{z}}_{(t|t),l}+\alpha {\mathbf{P}}_{(t|t),l}{\mathbf{P}}_{(t+1|t),l}^{-1}({\hat{\mathbf{z}}}_{t+1,l}-{\mathbf{z}}_{(t+1|t),l})}$ 12:      end for 13:   end for Update variational parameters 14:      for ${\displaystyle t=1,\cdots ,T}$ and ${\displaystyle l=1,\cdots ,L}$ do 15:              [Formula/expression omitted. See PDF] 16:              ${\displaystyle {\mathbf{m}}_{{\mathbf{x}}_{t,l}}={\mathbf{Q}}_{{\mathbf{x}}_{t,l}}({\hat{\mathbf{z}}}_{t,l}-\alpha {\hat{\mathbf{z}}}_{t-1,l}-\frac{1}{2}\mathbf{1}-{\tilde{\mathbf{\Omega }}}_{t,l}{\hat{\mathbf{D}}}^{\mathrm{\top }}{\mathbf{s}}_t+{\gamma }_x{\mathbf{P}}_x^{-1}{\mathit{\mu }}_x)}$ 17: ${\displaystyle v_{t,l}^{(j)}:=\beta |m_{{\mathbf{x}}_{t,l}}^{(j)}+{\hat{\mathbf{d}}}_j^{\mathrm{T}}{\mathbf{s}}_t|}$ 18: for ${\displaystyle j=1,\cdots ,N}$ do 19:              [Formula/expression omitted. See PDF] 20:              ${\displaystyle ({\tilde{\mathbf{\Omega }}}_{t,l}{)}_{j,j}:={\displaystyle \frac{1}{2c_{t,l}^{(j)}}\tanh (\frac{c_{t,l}^{(j)}}{2})}}$ 21:      end for 22:   end for 23:   ${\displaystyle {\mathbf{P}}_x:={\mathit{\psi }}_x+\sum _{t,l=1}^{T,L}\left\{{\mathbf{Q}}_{{\mathbf{x}}_{t,l}}+{\mathbf{m}}_{{\mathbf{x}}_{t,l}}{\mathbf{m}}_{{\mathbf{x}}_{t,l}}^{\mathrm{\top }}-{\mathit{\mu }}_x{\mathbf{m}}_{{\mathbf{x}}_{t,l}}^{\mathrm{\top }}-{\mathbf{m}}_{{\mathbf{x}}_{t,l}}{\mathit{\mu }}_x^{\mathrm{\top }}+{\mathit{\mu }}_x{\mathit{\mu }}_x^{\mathrm{\top }}\right\}}$ Update IRLS covariance approximation 24:      for ${\displaystyle l=1,\cdots ,L,t=1,\cdots ,T}$ and ${\displaystyle j=1,\cdots ,N}$ do 25:              ${\displaystyle ({\mathbf{\Sigma }}_{{\mathbf{v}}_{t,l}}{)}_{j,j}:=\sqrt{({\hat{z}}_{t,l}^{(j)}-\alpha {\hat{z}}_{t-1,l}^{(j)}{)}^2+{\epsilon }^2}/v_{t,l}^{(j)}}$ 26:      end for Update outputs and the convergence criterion 27:      for ${\displaystyle j=1,\cdots ,N}$ do 28:              ${\displaystyle {\hat{\mathbf{d}}}_j={\displaystyle {\left(\sum _{t,l=1}^{T,L}\left(({\tilde{\mathbf{\Omega }}}_{t,l}{)}_{j,j}{\mathbf{s}}_t{{\mathbf{s}}_t}^{\mathrm{\top }}\right)\right)}^{-1}\left(\sum _{t,l=1}^{T,L}\left\{\left({\hat{z}}_{t,l}^{(j)}-\alpha {\hat{z}}_{t-1,l}^{(j)}-\frac{1}{2}\right){\mathbf{s}}_t-({\tilde{\mathbf{\Omega }}}_{t,l}{)}_{j,j}{m}_{{\mathbf{x}}_{t,l}}^{(j)}{\mathbf{s}}_t\right\}\right)}}$ 29:      end for 30:      [Formula/expression omitted. See PDF], ${\displaystyle \hat{\mathbf{D}}=\left[{\hat{\mathbf{d}}}_1,{\hat{\mathbf{d}}}_2,\cdots ,{\hat{\mathbf{d}}}_N\right]}$ 31:      [Formula/expression omitted. See PDF], ${\displaystyle {\hat{\mathbf{\Sigma }}}_x=\frac{{\mathbf{P}}_x}{{\gamma }_x+N+1}}$ 32:      [Formula/expression omitted. See PDF] 33:   end while 34:   Return [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF]

Parameter choices for simulation study 2

In the second simulation study, we set $\alpha =0.98$, [Formula/expression omitted. See PDF], [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF] ([Formula/expression omitted. See PDF] represents the identity matrix and [Formula/expression omitted. See PDF] represents the vector of all ones) when generating the fluorescence traces [Formula/expression omitted. See PDF], so that the SNR of the simulated data was in the same range as of real calcium imaging observations. Furthermore, we simulated the spike trains based on a Poisson process (Smith and Brown, 2003) using the discrete time re-scaling procedure (Brown et al., 2002; Smith and Brown, 2003). Following the assumptions in Brown et al., 2002, we used an exponential link to simulate the observations:${\displaystyle n_{t,l}^{(j)}\sim \mathrm{Poisson}\left({\lambda }_{t,l}^{(j)}\right),{\lambda }_{t,l}^{(j)}=\exp \left(x_{t,l}^{(j)}\right),}$as opposed to the Bernoulli-logistic assumption in our recognition model. Then, we estimated the noise covariance [Formula/expression omitted. See PDF] using the Algorithm 1, with a slight modification. Since there are no external stimuli, we set [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF]. Accordingly, in Algorithm 1, we initialized [Formula/expression omitted. See PDF] and did not perform the update on [Formula/expression omitted. See PDF] in the subsequent iterations.

Parameter choices for real data study 1

The dataset consisted of recordings from 371 excitatory neurons, from which we selected $N=16$ responsive neurons for the analysis. Each trial consisted of $T=3600$ time frames (the sampling frequency was 30 Hz, and each trial had a duration of 120 s), with the presentation of a random sequence of four tones ($P=4$). The spiking events were very sparse and infrequent, and hence this dataset fits our model with at most one spiking event in a time frame.

We considered $R=25$ ($M=100$) time lags in this analysis and further examined the effect of varying $R$ in Figure 4—figure supplement 1. We set $\alpha =0.95$ and [Formula/expression omitted. See PDF] ([Formula/expression omitted. See PDF] represents the identity matrix).

Parameter choices for real data study 2

Each trial consisted of $T=765$ frames (25.5 s) at a sampling frequency of 30 Hz. The A1 neurons studied here had low response rates (in both time and space), with only $\sim 10$ neurons exhibiting spiking activity in at least half of the trials. Thus, we selected $N=10$ neurons and $L=10$ trials for the analysis, and chose $R=25$ lags of the stimulus ($M=25$) in the model for the stimulus-driven condition. We set $\alpha =0.95$ and [Formula/expression omitted. See PDF] ([Formula/expression omitted. See PDF] represents the identity matrix).

Parameter choices for real data study 3

Each experiment consisted of $L=5$ trials of $P=9$ different tone frequencies repeated at four different amplitude levels, resulting in each concatenated trial being $\sim 180$ s long (see Bowen et al., 2020 for more details). We set the number of stimulus time lags considered to be $R=25$ ($M=225$). For each layer, we analyzed fluorescence observations from six experiments. In each experiment, we selected the most responsive $N\sim 20$ neurons for the subsequent analysis. We set $\alpha =0.95$ and [Formula/expression omitted. See PDF].

Performance evaluation

Simulation studies

Since the ground truth is known in simulations, we directly compared the performance of each signal and noise correlation estimate with the ground truth signal and noise correlations, respectively. Suppose the ground truth correlations are given by the matrix [Character omitted. See PDF] and the estimated correlations are given by the matrix [Formula/expression omitted. See PDF]. To quantify the similarity between [Character omitted. See PDF] and [Formula/expression omitted. See PDF], we defined the following two metrics:

Normalized Mean Squared Error (NMSE): The NMSE computes the mean squared error of [Formula/expression omitted. See PDF] with respect to [Character omitted. See PDF] using the Frobenius Norm:${\displaystyle \mathrm{N}\mathrm{M}\mathrm{S}\mathrm{E}:=\frac{\Vert \mathbf{X}-\hat{\mathbf{X}}{\Vert }_F^2}{\Vert \mathbf{X}{\Vert }_F^2}.}$

Ratio between out-of-network power and in-network power (leakage): First, we identified the in-network and out-of-network components from the ground truth correlation matrix [Character omitted. See PDF]. Suppose that if the true correlation between the [Formula/expression omitted. See PDF] neuron and the [Formula/expression omitted. See PDF] neuron is non-zero, then [Formula/expression omitted. See PDF], for some ${\delta }_x>0$. Thus, we formed a matrix [Formula/expression omitted. See PDF] that masks the in-network components, by setting [Formula/expression omitted. See PDF] if [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF] if [Formula/expression omitted. See PDF]. Likewise, we also formed a matrix [Formula/expression omitted. See PDF] that masks the out-of-network components, by setting [Formula/expression omitted. See PDF] if [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF] if [Formula/expression omitted. See PDF]. Then, using these two matrices we quantified the leakage effect of [Formula/expression omitted. See PDF] comparative to [Character omitted. See PDF] by:

[Formula omitted. See PDF]

where $(\cdot )$ denotes element-wise multiplication.

Real data studies

To quantify the similarity and dissimilarity between signal and noise correlation estimates, we used a statistic based on the Tanimoto similarity metric (Lipkus, 1999), denoted by [Formula/expression omitted. See PDF] for two matrices [Character omitted. See PDF] and [Character omitted. See PDF]. For two vectors [Character omitted. See PDF] and [Character omitted. See PDF] with non-negative entries, the Tanimoto coefficient (Lipkus, 1999) is defined as:${\displaystyle T(\mathbf{a},\mathbf{b}):=\frac{{\mathbf{a}}^{\mathrm{T}}\mathbf{b}}{{\mathbf{a}}^{\mathrm{T}}\mathbf{a}+{\mathbf{b}}^{\mathrm{T}}\mathbf{b}-{\mathbf{a}}^{\mathrm{T}}\mathbf{b}}.}$

The Tanimoto similarly metric between two matrices can be defined in a similar manner, by vectorizing the matrices. Thus, we formulated a similarity metric between two correlation matrices [Character omitted. See PDF] and [Character omitted. See PDF] as follows. Let [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF], with the $\max \{\cdot ,\cdot \}$ operator interpreted element-wise. Note that [Formula/expression omitted. See PDF], and [Formula/expression omitted. See PDF] have non-negative entries. We then defined the similarity matrix by combining those of the positive and negative parts as follows:

[Formula omitted. See PDF]

where $\epsilon \in [0,1]$ denotes the percentage of positive entries in [Character omitted. See PDF] and [Character omitted. See PDF]. As a measure of dissimilarity, we used [Formula/expression omitted. See PDF]. The values of [Formula/expression omitted. See PDF] in Table 1 and [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF] reported in Table 2 were obtained based on the foregoing definitions.

To further assess the statistical significance of these results, we performed following randomized tests. To test the significance of [Formula/expression omitted. See PDF], for each comparison and each algorithm, we fixed the first matrix (i.e. [Formula/expression omitted. See PDF]) and randomly shuffled the entries of the second one (i.e. [Formula/expression omitted. See PDF]) while respecting symmetry. We repeated this procedure for 10000 trials, to derive the null distributions that represented the probabilities of chance occurrence of similarities between two random groups of neurons.

To test the significance of [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF], for each comparison and each algorithm, again we fixed the first matrix (i.e. signal correlations). Then, we formed the elements of the second matrix (akin to noise correlations) as follows. For each element of the second matrix, we assigned either the same element as the signal correlations (in order to model the leakage effect) or a random noise (with same variance as the elements in the noise correlation matrix) with equal probability. As before, we repeated this procedure for 10,000 trials, to derive the null distributions that represent the probabilities of chance occurrence of dissimilarities between two matrices that have some leakage between them.

Hyper-parameter tuning

The hyper-parameters that directly affect the proposed estimation are the inverse Wishart prior hyper-parameters: ${\displaystyle {\mathit{\psi }}_x}$ and ${\rho }_x$. Given that ${\rho }_x$ appears in the form of ${\gamma }_x:=TL+{\rho }_x$, we will consider ${\displaystyle {\mathit{\psi }}_x}$ and ${\gamma }_x$ as the main hyper-parameters for simplicity. Here, we propose a criterion for choosing these two hyper-parameters in a data-driven fashion, which will then be used to construct the estimates of the noise covariance matrix [Formula/expression omitted. See PDF] and weight matrix [Formula/expression omitted. See PDF]. Due to the hierarchy of hidden layers in our model, an empirical Bayes approach for hyper-parameter selection using a likelihood-based performance metric is not straightforward. Hence, we propose an alternative empirical method for hyper-parameter selection as follows.

For a given choice of [Formula/expression omitted. See PDF] and ${\gamma }_x$, we estimate [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF] following the proposed method. Then, based on the generative model in Proposed forward model, and using the estimated values of [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF], we sample an ensemble of simulated fluorescence traces [Formula/expression omitted. See PDF], and compute the metric [Formula/expression omitted. See PDF]:${\displaystyle d({\mathit{\psi }}_x,{\gamma }_x):=D_{\mathsf{f}\mathsf{r}\mathsf{o}\mathsf{b}}\left(\mathrm{cov}(\hat{\mathbf{y}},\hat{\mathbf{y}}),\mathrm{cov}(\mathbf{y},\mathbf{y})\right),}$where ${\displaystyle \mathrm{cov}(\cdot ,\cdot )}$ denotes the empirical covariance and [Formula/expression omitted. See PDF]. Note that [Formula/expression omitted. See PDF] is strictly convex in [Character omitted. See PDF]. Thus, minimizing [Formula/expression omitted. See PDF] over [Character omitted. See PDF] for a given [Character omitted. See PDF] has a unique solution. Accordingly, we observe that [Formula/expression omitted. See PDF] is minimized when [Formula/expression omitted. See PDF] is nearest to [Formula/expression omitted. See PDF]. Therefore, the corresponding estimates [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF] that generated [Formula/expression omitted. See PDF], best match the second-order statistics of [Character omitted. See PDF] that was generated by the true parameters [Formula/expression omitted. See PDF] and [Character omitted. See PDF].

The typically low spiking rate of sensory neurons observed in practice may render the estimation problem ill-posed. It is thus important to have an accurate choice of the scale matrix [Formula/expression omitted. See PDF] in the prior distribution. However, an exhaustive search for optimal tuning of [Formula/expression omitted. See PDF] is not computationally feasible, given that it has $N(N+1)/2$ free variables. Thus, the main challenge here is finding the optimal choice of the scale matrix [Formula/expression omitted. See PDF].

To address this challenge, we propose the following method. First, we fix [Formula/expression omitted. See PDF], where τ is a scalar and [Formula/expression omitted. See PDF] is the identity matrix. Next, given [Formula/expression omitted. See PDF] we find the optimal choice of ${\gamma }_x$ as:

[Formula omitted. See PDF]

where [Formula/expression omitted. See PDF] is a finite set of candidate solutions for ${\gamma }_x>N-1$. Let [Formula/expression omitted. See PDF] denote the noise covariance estimate corresponding to hyper-parameters [Formula/expression omitted. See PDF]. We will next use [Formula/expression omitted. See PDF] to find a suitable choice of [Formula/expression omitted. See PDF]. To this end, we first fix [Formula/expression omitted. See PDF], for some $N-1<{\tilde{\rho }}_x\ll TL$. Note that by choosing ${\tilde{\rho }}_x$ to be much smaller than $TL$, the final estimates become less sensitive to the choice of ${\gamma }_x$. Then, we construct a candidate set [Formula/expression omitted. See PDF] for [Formula/expression omitted. See PDF] by scaling ${\displaystyle {\hat{\mathbf{\Sigma }}}_{x,\mathsf{i}\mathsf{n}\mathsf{i}\mathsf{t}}}$ with a finite set of scalars ${\displaystyle \eta \in {\mathbb{R}}^{+}}$, i.e. [Formula/expression omitted. See PDF]. To select [Formula/expression omitted. See PDF], we match it with the choice of [Formula/expression omitted. See PDF] by solving:

[Formula omitted. See PDF]

Finally, we use these hyper-parameters [Formula/expression omitted. See PDF] to obtain the estimators [Formula/expression omitted. See PDF] and [Formula/expression omitted. See PDF] as the output of the algorithm.

Experimental procedures

All procedures were approved by the University of Maryland Institutional Animal Care and Use Committee. Imaging experiments were performed on a P60 (for real data study 1) and P83 (for real data study 2) female F1 offspring of the CBA/CaJ strain (The Jackson Laboratory; stock #000654) crossed with transgenic C57BL/6J-Tg(Thy1-GCaMP6s)GP4.3Dkim/J mice (The Jackson Laboratory; stock #024275) (CBAxThy1), and F1 (CBAxC57). The third real data study was performed on data from P66-P93 and P166-P178 mice (see Bowen et al., 2020 for more details). We used the F1 generation of the crossed mice because they have good hearing into adulthood (Frisina et al., 2011).

We performed cranial window implantation and two-photon imaging as previously described in Francis et al., 2018; Liu et al., 2019; Bowen et al., 2019. Briefly, we implanted a cranial window of 3 mm in diameter over the left auditory cortex. We used a scanning microscope (Bergamo II series, B248, Thorlabs) coupled to Insight X3 laser (Spectra-physics) (study 1) or pulsed femtosecond Ti:Sapphire two-photon laser with dispersion compensation (Vision S, Coherent) (studies 2 and 3) to image GCaMP6s fluorescence from individual neurons in awake head-fixed mice with an excitation wavelengths of $\lambda =920$ nm and $\lambda =940$ nm, respectively. The microscope was controlled by ThorImageLS software. The size of the field of view was ${\displaystyle 370\times 370\mu m}$. Imaging frames of $512\times 512$ pixels (pixel size ${\displaystyle 0.72\mu m}$) were acquired at 30 Hz by bidirectional scanning of an 8 kHz resonant scanner. The imaging depth was around ${\displaystyle 200\mu m}$ below pia.

Data pre-processing

A circular ROI was manually drawn over each cell body to extract raw fluorescence traces from individual cells. Neuropil contamination subtraction and baseline correction were performed on the raw fluorescence traces of each cell (Francis et al., 2018; Liu et al., 2019; Bowen et al., 2020) according to [Formula/expression omitted. See PDF], where ${\alpha }_n$ was set to 0.7 in real data study 1 (Francis et al., 2018), 0.8 in real data study 2 (Liu et al., 2019) and 0.9 in real data study 3 (Bowen et al., 2020). The two-photon observations [Formula/expression omitted. See PDF] used in our analyses are the output of this pre-processing step.

Stimuli for real data study 1

During imaging experiments, we presented four tones (4, 8, 16, and 32 kHz) at 70 dB SPL. The tones were 2 s in duration with an inter-trial silence of 4 s. For the sequence of tones, we first generated a randomized sequence that consisted of five repeats for each tone (20 tones in total) and then the same sequence was repeated for 10 trials.

Stimuli for real data study 2

During imaging experiments, we presented a 75 dB SPL 100 ms broadband noise (4–48 kHz) as the auditory stimulus. Each trial was 5.1 s long (1 s pre-stimulus silence + 0.1 s stimulus + 3 s post-stimulus silence), and the inter-trial duration was 3 s. Spontaneous neuronal activity was collected from activity during randomly interleaved no-stimuli trials of the same duration, and these trials had complete silence throughout the trial duration (5.1 s long).

Then, we extracted 50 such trials from each type, and formed 10 ($L=10$) trials each of 25.5 s duration ($T=765$ frames) for the subsequent analysis, by concatenating five 5.1 s trials. This final step was performed to increase the effective trial duration.

Stimuli for real data study 3

During imaging experiments, sounds were played at four sound levels (20, 40, 60, and 80 dB SPL). Auditory stimuli consisted of sinusoidal amplitude-modulated (SAM) tones (20 Hz modulation, cosine phase), ranging from 3 to 48 kHz. The frequency resolution was two tones/octave (0.5 octave spacing) and each of these tonal stimuli was 1 s long, repeated five times with a 4−6 s inter-stimulus interval (see Bowen et al., 2020 for details).