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

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;Josic´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 s t , 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 l th trial from the j th neuron at time frame t, denoted by y ðjÞ t;l , is a noisy surrogate of the intracellular calcium concentrations. The calcium concentrations in turn are temporally blurred surrogates of the underlying spiking activity n ðjÞ t;l , as shown in Figure 1. In modeling the spiking activity, we consider two main contributions: (1) the common known stimulus s t affects the activity of the j th neuron via an unknown kernel d j , 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 s t are captured by a trial-dependent latent process x ðjÞ t;l . 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:  where fðÁÞ 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 timelocked to the repeated stimulus, s t . 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 ðS s Þ i;j and noise ðS x Þ i;j covariance between the i th neuron and j th neuron as: ðS s Þ i;j :¼d > i cov s t ; s t ð Þd j ; ðS x Þ i;j :¼cov x ðiÞ t;l ; x ðjÞ t;l ; (1) where covðÁ; ÁÞ is the empirical covariance function defined for two vector time series u t and v t as

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 S s and S x from the fluorescence observations, without requiring intermediate spike deconvolution (Figure 1, backward arrow). The signal and noise correlation matrices, denoted by S and N, can then be obtained by standard normalization of S s and S x : ; ðNÞ i;j :¼ ðS x Þ i;j ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðS x Þ i;i :ðS x Þ j;j q ; 8i; j ¼ 1; 2; Á Á Á ; N: We note that when spiking activity is directly observed using electrophysiology recordings, the conventional signal ðS con s Þ i;j and noise ðS con x Þ i;j covariances of spiking activity between the i th and j th neuron are defined as (Lyamzin et al., 2015): which after standard normalization in Equation 2 give the conventional signal ðS con Þ i;j and noise ðN con Þ i;j 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: S con » S and N con » 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 ðjÞ t;l in Equation 3 by y ðjÞ t;l , 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.
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 xt;1 , 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.
1. 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 com- , where e d j;1 and e d j;2 are akin to quadratic receptive field components. We assumed a linear stimulus integration model in our estimation framework (i.e., e d j;1 ¼ e d j;2 ¼ 0). 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. 2. 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. 3. 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 Simulation study 2: spontaneous activity Next, we present the results of a simulation study in the absence of external stimuli (i.e. s t ¼ 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 s t and the receptive field d j 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  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 twophoton 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 innetwork power (leakage) are shown below each panel.
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).
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 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 T s ðX; YÞ for two matrices X and Y. As a measure of dissimilarity, we used T d ðX; YÞ :¼ 1 À T s ðX; YÞ (see Performance evaluation for details). The comparison of T d ð b S; b NÞ 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 T d ð b S; b NÞ values. 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 < 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~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 b z ð1Þ t;1 (calcium concentration) and b n ð1Þ t;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 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.
NÞ Regression statistics ( Figure 4E) Shuffling test ( Figure 5)  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 b n ðjÞ t;l :¼ b z ðjÞ t;l À ab z ðjÞ tÀ1;l 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 b n ðjÞ t;l 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.
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 twostage 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 (ð b NÞ 1;3 ) 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 6B shows the resulting noise and signal correlation estimates under the spontaneous ( b N spon , top) and stimulus-driven ( b N stim and b S stim , 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 b N spon and b N stim are saliently similar, and distinct from b S stim .
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 T s ðX; YÞ and the dissimilarity metric T d ðX; YÞ (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  Table 2. Similarity/dissimilarity metric statistics for the estimates in Figure 6. 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). 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 timelocked 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 twophoton 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 10 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).
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).
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 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.   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 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 d j for j ¼ 1; Á Á Á ; N and the noise covariance S x . In order to simplify the treatment, we posit the following mild assumptions: Assumption (1). We assume a scalar time-varying external stimulus (i.e. s t ¼ s t , and hence . Furthermore, we set the observation noise covariance to be S w ¼ s 2 w I, 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. x t;l ¼ x Wk;l and s t ¼ s Wk , for ðk À 1ÞW þ 1 t<kW and k ¼ 1; Á Á Á ; 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> 1 64 , 0<<1=2, and 0<h 1=2 be fixed constants, Then, under Assumptions (1 -3), the bias and variance of the maximum likelihood estimators where t j and e t j denote bounded terms that are Oðs 2 w Þ or O 1 W À Á , i;j and e i;j denote bounded terms that are Oðs 2 w Þ or O 1 W 1À2 À Á and C 1 ; C 2 ; C 3 and C 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 x t;l and the true kernel d j and forms the correlation estimates. In this case, the oracle would incur zero bias and variance of order O 1=KL ð Þ in estimating d j , and outputs an estimate of S x with bias and variance in the order of O 1=KL ð Þ. Theorem 1 indeed states that for sufficiently large W and small s 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 ! ¥ (and hence the observation duration T ! ¥), the loss of performance of the proposed estimators only depends on s 2 w , the variance of the observation noise. As a result, at a given observation noise variance s 2 w , 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 s 2 w 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 ! ¥), the oracle estimates become arbitrarily close to the true parameters S x and d j . 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 mm 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 twophoton 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 < 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 t;l > 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 y t;l by the following state-space model: where A 2 R NÂN represents the scaling of the observations, w t;l is zero-mean i.i.d. Gaussian noise with covariance S w , and 0 a<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, n ðjÞ t;l 2 f0; 1g. 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): t;l > . While we assume the stimulus s t 2 R M to be common to all neurons, we model the distinct effect of this stimulus on the j th neuron via an unknown kernel d j 2 R M , 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: Finally, we assume the latent trial dependent covariates to be a Gaussian process x t;l~N ðm x ; S x Þ, with mean m x :¼ ½ ð1Þ x ; ð2Þ x ; Á Á Á ; ðNÞ x > and covariance S x . 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 S x . The signal covariance matrix, representing the covariance of the neural activity arising from the repeated application of the stimulus s t , is given by The signal and noise correlation matrices, denoted by S and N, can then be obtained by standard normalization of S s and S x : The main problem is thus to estimate fS x ; Dg from the noisy and temporally blurred data y t;l È É T;L t¼1;l¼1 .

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 p pr ðS x Þ on the noise covariance, to compensate sparsity of data. A natural estimation method to estimate fS x ; Dg in a Bayesian framework is to maximize the observed data likelihood p À fy t;l g T;L t;l¼1 S x ; 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: Inspecting this MAP problem soon reveals that estimating S x and D 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 nonlinearities 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 fx t;l g T;L t;l¼1 and S x as latent variables and fz t;l g T;L t;l¼1 and D as unknown The calcium concentration at time t is a function of the spiking activity n t;l , and the calcium activity at the previous time point z tÀ1;l . The spiking activity is driven by two independent mechanisms: latent trial-dependent covariates x t;l , and contributions from the known external stimulus s t , which we model by D > s t (in which the receptive field D is unknown). Then, we model x t;l as a Gaussian process with constant mean m x , and unknown covariance S x . Finally, we assume the covariance S x to have an inverse Wishart prior distribution with hyper-parameters c x and x . Based on this forward model, the inverse problem amounts to recovering the signal and noise correlations by directly estimating S x and D (top layer) from the fluorescence observations y t;l È É T;L t¼1;l¼1 (bottom layer).
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 y :¼ fy t;l g T;L t;l¼1 , z :¼ fz t;l g T;L t;l¼1 , and x :¼ fx t;l g T;L t;l¼1 .

Preliminary assumptions
For the sake of simplicity, we assume that the constants a, A, S w and m x are either known or can be consistently estimated from pilot trials. Next, we take p pr ðS x Þ to be an Inverse Wishart density: which turns out to be the conjugate prior in our model. Thus, c x and 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 x t;l in our model with a standard distribution is not straightforward, due to the complicated posterior arising from codependent 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, pðxjz; D; S x Þ is conditionally independent in t; l with: 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 v t;l :¼ ½! ð1Þ t;l ; ! ð2Þ t;l ; Á Á Á ; ! ðNÞ t;l > , ! ðjÞ t;l~P Gð1; 0Þ for 1 j N, 1 t T and 1 l L, to derive the complete data log-likelihood: where v :¼ v t;l È É T;L t;l¼1 and C accounts for terms not depending on y; z; x; v, S x and D. The complete data log-likelihood is notably quadratic in z t;l , 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 and S x , assuming that the parameter estimates b z and b D of the previous iteration are available. The methods that we propose to update the parameters b z and b D subsequently, will be discussed in the next section.
The objective of variational inference is to posit a family of approximate densities Q over the latent variables, and to find the member of that family that minimizes the Kullback-Leibler (KL) divergence to the exact posterior: 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): ELBO ðqÞ ¼ E½log pðx; v; S x ; yjb z; b DÞ À E½log qðx; v; S x jb z; b DÞ: Further, we assume Q to be a mean-field variational family (Blei et al., 2017), resulting in the overall variational density of the form: 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: Upon evaluation of these expectations, we derive the optimal variational distributions as: whose parameters m xt;l :¼ ½m ð1Þ xt;l ; m ð2Þ xt;l ; Á Á Á ; m ðNÞ xt;l T , Q xt;l , c ðjÞ t;l , P x , and g x can be updated given parameter estimates b D and b z:

Low-complexity parameter updates
Note that even though z is composed of the latent processes z t;l , 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 z and D using the updated variational density q Ã ðx; v; S x Þ. First, note that the log-likelihood in Equation 5 is decoupled in l, which admits independent updates to fz t;l g T t¼1 , for l ¼ 1; Á Á Á ; L. As such, given an estimate b D, we propose to estimate fz t;l g T t¼1 as: Instead, we consider an alternative solution that admits a low-complexity recursive solution by relaxing the constraints. To this end, we relax the constraint z t;l À az tÀ1;l " 1 and replace the constraint z t;l À az tÀ1;l # 0 by penalty terms proportional to jz ðjÞ t;l À az ðjÞ tÀ1;l j. The resulting relaxed problem is thus given by: 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 z t;l , followed by Fixed Interval Smoothing (FIS) (Rauch et al., 1965) to find the minimizer. At iteration k, given a current estimate z ½kÀ1 , the problem reduces to a Gaussian state-space estimation of the form: t;l , for some small constant ">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 z t;l 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 q Ã ðx; v; S x Þ and the updated b z, the estimate of d j for j ¼ 1; 2; Á Á Á ; N can be updated in closed-form by maximizing the expected complete log-likelihood E q Ã ðx;v;SxÞ log pðy; b z; x; v; S x jDÞ ½ : t;l À ab z ðjÞ tÀ1;l À 1 2 s t À ð e V t;l Þ j;j m ðjÞ xt;l s t ' ! : The VI procedure iterates between updating the variational densities and parameters until convergence, upon which we estimate the noise and signal covariances as: 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 fn t;l g T;L t;l¼1 , is directly observed in this case, the solution to the optimization problem in Equation (7) is no longer required. Thus, the parameters S x and D 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 p th tone in the stimulus set (p ¼ 1; Á Á Á ; P, where P is the number of distinct tones) is given by a binary sequence f ðpÞ t 2 f0; 1g. 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 s tÀRþ1 > 2 R R . By including all the P tones, the overall effective stimulus at the t th time frame is given by The stimulus modulation vector d j 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 < 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~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 (S w ) and scaling matrix (A)
We assumed that the observation noise covariance S w is diagonal, and estimated the diagonal elements using the background fluorescence in the absence of spiking events, for each neuron. We set A ¼ aI, where I 2 R NÂN 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 (a)
We chose a 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 a in this range through the method outlined in Hyper-parameter tuning, and accordingly chose the optimal value of a.

Mean of the latent trial-dependent process (m x )
We estimated m x as a constant that is proportional to the average firing rate. To this end, we param- The constants a and b were chosen such that À2 ðjÞ x À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 a ¼ 0:98, b ¼ 8, A ¼ 0:1I, m x ¼ À4:51 and S w ¼ 2 Â 10 À4 I (I 2 R 8Â8 represents the identity matrix and 1 2 R 8 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 6 th 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. s t ¼ ½s t ; s tÀ1 > ) in both the generative model and inference procedure. The components of the linear and quadratic stimulus modulation vectors, that is d j , e d j;1 and e d j;2 , 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 (d > j s t for the linear model, and d > j s t þ ð e d j;1 > s t Þ 2 þ ð e d j;2 > s t Þ 2 for the non-linear model) was comparable to the average power of the noise component (x ðjÞ t;l ).

Parameter choices for simulation study 2
In the second simulation study, we set a ¼ 0:98, A ¼ 0:1I, m x ¼ À4:51 and S w ¼ 10 À4 I (I 2 R 30Â30 represents the identity matrix and 1 2 R 30 represents the vector of all ones) when generating the fluorescence traces y t;l È É T;L t;l¼1 , 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: as opposed to the Bernoulli-logistic assumption in our recognition model. Then, we estimated the noise covariance b S x using the Algorithm 1, with a slight modification. Since there are no external stimuli, we set s t ¼ 0 and D ¼ 0. Accordingly, in Algorithm 1, we initialized b D ¼ 0 and did not perform the update on b D 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 a ¼ 0:95 and A ¼ I (I 2 R 16Â16 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~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 a ¼ 0:95 and A ¼ 0:75I (I 2 R 10Â10 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~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~20 neurons for the subsequent analysis. We set a ¼ 0:95 and A ¼ I.

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 X and the estimated correlations are given by the matrix b X. To quantify the similarity between X and b X, we defined the following two metrics: Normalized Mean Squared Error (NMSE): The NMSE computes the mean squared error of b X with respect to X using the Frobenius Norm: Ratio between out-of-network power and in-network power (leakage): First, we identified the innetwork and out-of-network components from the ground truth correlation matrix X. Suppose that if the true correlation between the i th neuron and the j th neuron is non-zero, then X ð Þ i;j >d x , for some d x >0. Thus, we formed a matrix X in that masks the in-network components, by setting X in À Á Likewise, we also formed a matrix X out that masks the out-ofnetwork components, by setting X out ð Þ Then, using these two matrices we quantified the leakage effect of b X comparative to X by: where ðÁÞ 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 T s ðX; YÞ for two matrices X and Y. For two vectors a and b with non-negative entries, the Tanimoto coefficient (Lipkus, 1999) is defined as: 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 X and Y as follows. Let X þ :¼ maxfX; 0Ig and X À :¼ maxfÀX; 0Ig, with the maxfÁ; Ág operator interpreted element-wise. Note that X ¼ X þ À X À , and X þ ; X À have non-negative entries. We then defined the similarity matrix by combining those of the positive and negative parts as follows: where " 2 ½0; 1 denotes the percentage of positive entries in X and Y. As a measure of dissimilarity, we used T d ðX; YÞ :¼ 1 À T s ðX; YÞ. The values of T d ð b S; b NÞ in Table 1 and T s ð b N spon ; b N stim Þ and 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 T s ð b N spon ; b N stim Þ, for each comparison and each algorithm, we fixed the first matrix (i.e. b N spon ) and randomly shuffled the entries of the second one (i.e. b N stim ) 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 T d ð b S; b NÞ and T d ð b S stim ; b N stim Þ, 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: c x and x . Given that x appears in the form of g x :¼ TL þ x , we will consider c x and g 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 b S x and weight matrix b D. 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 hyperparameter selection as follows.
For a given choice of c x and g x , we estimate b S x and b D following the proposed method. Then, based on the generative model in Proposed forward model, and using the estimated values of b S x and b D, we sample an ensemble of simulated fluorescence traces b y ¼ fb y ðlÞ t g T;L t;l¼1 , and compute the metric d c x ; g x ð Þ: dðc x ; g x Þ :¼ D frob covðb y;b yÞ; covðy; yÞ ð Þ; where covðÁ; ÁÞ denotes the empirical covariance and D frob ðX; YÞ :¼ kX À Yk 2 F . Note that D frob ðX; YÞ is strictly convex in X. Thus, minimizing D frob X; Y ð Þ over X for a given Y has a unique solution. Accordingly, we observe that d c x ; g x ð Þ is minimized when covðb y;b yÞ is nearest to covðy; yÞ. Therefore, the corresponding estimates b S x and b D that generated b y, best match the second-order statistics of y that was generated by the true parameters S x and D.
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 c x in the prior distribution. However, an exhaustive search for optimal tuning of c x 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 c x;opt .
To address this challenge, we propose the following method. First, we fix c x;init ¼ t I, where t is a scalar and I 2 R NÂN is the identity matrix. Next, given c x;init we find the optimal choice of g x as: argmin dðc x;init ; g x Þ ; where S g is a finite set of candidate solutions for g x >N À 1. Let b S x;init denote the noise covariance estimate corresponding to hyper-parameters c x;init ; g x;init . We will next use b S x;init to find a suitable choice of c x . To this end, we first fix g x;opt :¼ TL þ e x , for some N À 1<e x ( TL. Note that by choosing e x to be much smaller than TL, the final estimates become less sensitive to the choice of g x . Then, we construct a candidate set S for c x;opt by scaling b S x;init with a finite set of scalars h 2 R þ , i.
To select c x;opt , we match it with the choice of g x;opt by solving: argmin d c x ; g x;opt À Á : Finally, we use these hyper-parameters c x;opt ; g x;opt to obtain the estimators b S x and b D 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 l ¼ 920 nm and l ¼ 940 nm, respectively. The microscope was controlled by ThorImageLS software. The size of the field of view was 370 Â 370 m. Imaging frames of 512 Â 512 pixels (pixel size 0:72 m) were acquired at 30 Hz by bidirectional scanning of an 8 kHz resonant scanner. The imaging depth was around 200 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) accord-

ing to
FcellÀanFneuropilÀbaseline baseline , where a 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 y t;l È É T;L t;l¼1 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).
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Author contributions
Anuththara Rupasinghe, Ji Liu, Conceptualization, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing -original draft, Writing -review and editing; Nikolas Francis, Conceptualization, Resources, Data curation, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing -original draft, Project administration, Writing -review and editing; Zac Bowen, Data curation, Investigation, Writingreview and editing; Patrick O Kanold, Conceptualization, Resources, Supervision, Funding acquisition, Validation, Investigation, Project administration, Writing -review and editing; Behtash Babadi, Conceptualization, Resources, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing -original Suppose that the spiking events follow the forward model: where f : R 2 ! ½0; 1 is a differentiable non-linear mapping. We assume x t;l and s t to be independent. Without loss of generality, let E s t ½ ¼ 0 and E x t;l Â Ã ¼ m x . Further, we define the notation X t » Y t to denote almost sure equivalence, that is X t À ! a:s: Z and Y t À ! a:s: Z for some random variable Z.
First, let us consider ðS con Þ i;j . Noting that E n where ðiÞ t;l and ðjÞ t;l represent the higher order terms. Then, as L ! ¥, we get: where C S is a constant and t;l is typically small if the latent process x t;l and the stimulus s t are concentrated around their means. Then, the signal correlations are obtained by normalization of the signal covariance as in Equation 9, through which the scaling factor C S cancels and we get: ðS con Þ i;j » ðSÞ i;j : Thus, as T; L ! ¥, we see that S is indeed the signal correlation matrix that is aimed to be approximated by the conventional definitions.
Next, let us consider ðN con Þ i;j . Similar to foregoing analysis of the signal covariance, as T ! ¥ we get: Then, from a Taylor series expansion, we get: from the law of large numbers. Accordingly, we see that: where C N is a constant and t;l is typically small if the latent process x t;l and the stimulus s t are concentrated around their means. Then, the noise correlations are derived by normalization of the noise covariance given in Equation 9. This cancels out the scaling factor C N , and we get: Thus, we similarly conclude that as T; L ! ¥, the conventional definition of noise correlation N con indeed aims to approximate N.
As a numerical illustration, we demonstrated in Figure 2-figure supplement 2 that the conventional definitions of the correlations indeed approximate our proposed definitions, but require much larger number of trials to be accurate. More specifically, in order to achieve comparable performance to our method using L ¼ 20 trials, the conventional correlation estimates require L ¼ 1000 trials.

Proof of Theorem 1
In what follows, we present a comprehensive proof of Theorem 1. Recall the following key assumptions: Assumption (1). We assume a scalar time-varying external stimulus (i.e. s t ¼ s t , and hence d j ¼ . Furthermore, we set the observation noise covariance to be S w ¼ s 2 w I, 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 (i.e. p pr ðS x Þ / 1), 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., x t;l ¼ x Wk ;l and s t ¼ s Wk , for ðk À 1ÞW þ 1 t<kW and k ¼ 1; Á Á Á ; 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.
Proof of Theorem 1. First, recall the proposed forward model (see Proposed forward model) under Assumption (1 -3): y t;l ¼ Az t;l þ w t;l ; z t;l ¼ a z tÀ1;l þ n t;l ; n ðjÞ t;l~B ernoulli f x ðjÞ Wk ;l ; where f Á ð Þ :¼ expðÁÞ 1þexpðÁÞ , is the logistic function. Note that we have re-defined the latent process x t;l by absorbing the stimulus activity s t d to the mean of x t;l for notational convenience, without loss of generality. Hereafter, we also assume that A ¼ I without loss of generality. For a truncation level B (to be specified later), consider the event such that n Wk ;l ¼ n Wk;l ; n Wk;l ; Á Á Á ; n ðNÞ Wk;l h i > :¼ 1 W P W w¼1 n ðkÀ1ÞWþw;l . First, we derive convenient forms of the maximum likelihood estimators via the Laplace's approximations and asymptotic expansions (Wong, 2001) through the following lemma:

Lemma 1
Conditioned on event A W , the maximum likelihood estimators of the stimulus kernel of the j th neuron and the noise covariance between the i th and j th neurons take the forms: where W k ¼ ðk À 1ÞW þ 1. Then, we simplify these integrals based on the saddle point method of asymptotic expansions (Wong, 2001). To that end, first consider the numerator of Equation 10 denoted by I ð1Þ num . First, we evaluate the integration in I ð1Þ num with respect to the variable n. To that end, note: Wk ;l À ðjÞ t;l À P t k¼1 a tÀk n ðjÞ t;l 2 and dn is shorthand notation for the product measure of the discrete random vector n. Observing that rf 1 ðb nÞ ¼ 0 for b n :¼ b n t;l ¼ y t;l À ay tÀ1;l È É T;L t;l¼1 , using the method of asymptotic expansions, I ð1Þ num can be evaluated as: Wk ;l À ðjÞ x ; and A 2 ¼ W. Then, we note that the gradient of f 2 , rf 2 ðb xÞ ¼ 0 for Further, following the same sequence of reasoning, simplifying the numerator ðI ð2Þ num Þ and denominator ðI ð2Þ den Þ of Equation 11 yields: This concludes the proof of Lemma 1. & Given that f À1 ðzÞ is unbounded for z ¼ 0 or z ¼ 1, we consider another truncation: f À1 B 0 ðzÞ :¼ minfmaxff À1 ðzÞ; ÀB 0 g; B 0 g, where B 0 ¼ 2 log ð2 exp ðBÞ þ 1Þ. This choice of B 0 guarantees that over A W , f À1 B 0 n ðjÞ Wk ;l <B 0 for all j ¼ 1; Á Á Á ; N, k ¼ 1; Á Á Á ; K and l ¼ 1; Á Á Á ; L: and thus f À1 Wk ;l À ðjÞ x (i.e., observing x t;l directly) is unbiased and ðbÞ follows through the application of Jensen's inequality and triangle inequality. To simplify this bound, the triangle inequality yields: From the convexity of " e n ; The foregoing sequence of reasoning similarly follows for E f À1 Wk ;l j<B for k ¼ 1; Á Á Á ; K, l ¼ 1; Á Á Á ; L and j ¼ 1; Á Á Á ; N, conditioned on A W ). Accordingly, we derive the upper bound on the second term in Equation 19, conditioned on event A W : Then, we upper-bound the conditional second moment of e d j À ðd Oracle Þ j using the same techniques as we used in bounding the first moment. Accordingly, we get: where in ðjÞ, we have used the Cauchy-Schwarz inequality and ðkÞ follows from where ðlÞ follows from the triangle inequality, with the Oracle noise covariance estimator (i.e., observing x t;l directly), being defined as: where m ¼ km x k ¥ and we have used B 0 ¼ 2 log 2 exp ðBÞ þ 1 ð Þ. Similarly, conditioned on the event A W the second term in Equation 37 can be bounded as: Then, by combining the bounds in Equation 38 and Equation 40 and using an instance of Cauchy-Schwarz inequality P K k¼1 s Wk j j À Á 2 K P K k¼1 s 2 Wk , we see that the bound in Equation 37 conditioned on the event A W can be expressed as: B 0 e n ðiÞ Wk;l À ðiÞ x À s Wk Wk ;l À ðjÞ x À s Wk 1 L P K k 00 ¼1 sW k 00 2 P K;L k 00 ;l 00 ¼1 s W k 00 f À1 B 0 e n ðjÞ W k 00 ;l 00 À ðjÞ x À x ðiÞ Wk;l À ðiÞ x À s Wk Wk ;l À ðjÞ x À s Wk 1 L P K k 00 ¼1 sW k 00 2 P K;L k 00 ;l 00 ¼1 s W k 00 x ðjÞ W k 00 ;l 00 À ðjÞ s W k 0 s W k 00 n f À1 B 0 e n ðiÞ W k 0 ;l 0 À ðiÞ x f À1 B 0 e n ðjÞ W k 00 ;l 00 À ðjÞ x À x ðiÞ W k 0 ;l 0 À ðiÞ x x ðjÞ W k 00 ;l 00 À ðjÞ Â 1 þ s Wk L P K k 0 ¼1 s W k 0 L P K k 0 ¼1 s W k 0 2 þ s Wk L P K k 00 ¼1 s W k 00 L P K k 00 ¼1 s W k 00 2 þ s 2 Wk L P K k 0 ¼1 s W k 0 L P K k 00 ¼1 s W k 00 L 2 P K k 0 ¼1 s W k 0 2 P K k 00 ¼1 s W k 00 2 ! 2 : where the last inequality follows from an instance of the Cauchy-Schwarz inequality, that is P K k¼1 s Wk j j À Á 2 K P K k¼1 s 2 Wk . Then, following the observation KL S Oracle~I nvWish N ðS x ; KL À 1Þ, we derive the variance of S Oracle ð Þ i;j : ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where C 4 :¼ 384s m ffiffi ffi q p .
Finally, it only remains to prove that the event A W occurs with high probability for sufficiently large W: Lemma 3. The probability of occurrence of the event Next, we bound the probabilities on the right hand side using Chernoff's inequality (Boucheron et al., 2013). First, note that: where ðpÞ follows if B>2 km x k ¥ þ max k;j js Wk d j j È É À Á (which will hold under the conditions in Equation 29) and ðqÞ has been derived by applying the Chernoff's bound on the Gaussian random variable x ðjÞ Wk ;l . From the same reasoning, we see that P x