Dimensionality reduction of calcium-imaged neuronal population activity

Model overview

Our central goal is to develop dimensionality reduction methods appropriate for calcium imaging recordings to capture the shared, time-varying activity among the population of recorded neurons. To do so, we systematically compare three approaches.

For the first approach, we applied a standard dimensionality reduction method directly to the recorded fluorescence traces from calcium imaging (Fig. 2a, top). Here we chose to use a latent Linear Dynamical System (LDS), which is among the most basic methods for extracting neural trajectories. Conceptually, an LDS seeks to explain the temporal structure in the data using latent variables that vary smoothly over time.

Each time a neuron spikes, intracellular free calcium increases, then decays slowly over time. The calcium indicator kinetics influence the measured decay time, resulting in fluorescence traces whose intensity decays over hundreds of milliseconds to seconds, depending on the particular calcium indicator used [13]. This decay transient induces temporal correlations in the fluorescence measurements which are input to the LDS, which might attempt to capture these correlations in its latent variable estimates. This motivates our second approach, which first deconvolves each fluorescence trace separately, and then applies a LDS to the resulting estimated spiking activity. The deconvolution serves to remove a substantial portion of the calcium decay transient, producing activity traces similar to spike trains (or time varying firing rates). We term this two-stage method deconv-LDS (Fig. 2a, middle). Here we chose to deconvolve the fluorescence traces with OASIS [34, 32], which has been widely used in calcium imaging studies [23, 20, 35, 24].

In deconv-LDS, each neuron is deconvolved independently, and the stages of deconvolution and dimensionality reduction are performed sequentially. We asked if performing these stages jointly would lead to more accurate recovery of the latent variables (Fig. 2a, bottom). More specifically, we hypothesized that a sequential method (e.g., deconv-LDS) may inadvertently discard some of the shared activity amongst neurons due to the independent deconvolution of each neuron (Fig. 2b, left). Since the latent variables are intended to capture the shared spiking activity among neurons (and not the calcium decay, which is independent across neurons), it might be possible to better separate the calcium decay from the latent dynamics by considering all the neurons together, and performing the two stages of deconvolution and dimensionality reduction jointly. This allows the dimensionality reduction component to influence the deconvolution estimates, and vice versa, to more accurately estimate the ground truth latent variable than deconv-LDS (Fig. 2b, right). Thus, for our third method we developed a unified approach, CILDS, in which dimensionality reduction and deconvolution are performed jointly (Fig. 2a, bottom).

Deconvolution with dimensionality reduction

We first asked whether it is beneficial to first deconvolve fluorescence to estimate spiking activity in order to recover the underlying latent variables[27]. Indeed, multiple studies have applied dimensionality reduction to deconvolved spiking estimates [19, 20, 21, 22, 35, 23, 24, 36]. However, deconvolution is subject to particular statistical modeling assumptions (as is any statistical method) and usually does not recover the underlying spikes exactly. Therefore it is unclear how, or even if, deconvolution should be used with dimensionality reduction. We addressed this question by comparing the three approaches described above (Fig. 2a).

With calcium imaging recordings, the “ground truth” latent variables are unknown. Hence to directly compare each method’s ability to extract latent variables, we designed a simulation framework in which we created known ground truth latent variables with smoothly-varying time courses. These latent variables were used to generate spike trains which, in turn, were used to generate fluorescence traces (see Methods). We then applied each of the three approaches to these simulated fluorescence traces to assess how accurately they recovered the ground truth latent variables. Examples of two combinations of experimental variables are illustrated in Fig. 3a. In the simulations, we used 10 latent variables, of which two are shown.

We found that in both settings, CILDS outperformed the other two methods, returning more accurate estimates of the ground truth latent variables (Fig. 3b–c, points above the diagonal). By operating on the entire population of neurons together, CILDS was better able to separate the calcium transients from the shared activity among neurons (i.e., the latent variables) compared to deconv-LDS, which deconvolves the activity of each neuron individually, and LDS, which makes no attempt at this separation. Taken together, when extracting single-trial neural trajectories, one should use deconvolution jointly with dimensionality reduction, as in CILDS.

Impact of experimental variables on latent variable recovery

We next asked how different experimental variables affect the accuracy of the recovered latent variables. The simulation framework enables us to systematically vary experimentally relevant variables. These variables comprise four axes along which we can explore different combinations of values that might mimic a particular experimental paradigm. We varied the timescales of the latent time courses from 50ms to 5000ms, the number of neurons from 20 to around 100, the calcium decay timescales to match GCaMP6f, GCaMP6m, and GCaMP6s (fast, medium, slow) [13], and the variance of the added imaging noise that was independent of the calcium and spiking activity (Supplementary Table 1). We assessed how the accuracy of the dimensionality reduction methods changed as we systematically changed these variables.

As we increased the timescales of the latent variables, all dimensionality reduction methods improved their accuracy in estimating the ground truth latent variables (Fig. 3d). This occurs because with slower latent fluctuations, the latent variables become less independent across time. As a result, all methods can leverage future and past time points to better estimate the latent variables at the current time point. When the latent timescales are slow (order of seconds), the calcium indicator decay (which also has a timescale on the order of seconds) tends not to blur the underlying neural activity. Therefore a method that is not able to disambiguate between the time courses of the latent variables and the calcium decay (e.g., LDS) can still accurately recover the latent variables (see latent timescale of 5000 ms in Fig. 3d). However, at faster latent timescales (tens to hundreds of milliseconds) which reflect the timescales of many sensory, cognitive, and motor functions [37, 35], it is critical to use a method that accounts for the calcium decay, as both CILDS and deconv-LDS do (see latent timescale of 50 ms in Fig. 3d).

Next, when we increased the number of “recorded” neurons, all three methods improved in their ability to reconstruct the ground truth latent variables (Fig. 3e). This makes sense because each neuron provides a different, noisy view of the underlying latent variables. With more neurons, all methods are better able to “triangulate” the values of the latent variables.

We also varied the time constant of the GCaMP calcium indicator decay to match GCaMP6f, 7 6m, and 6s (from fast to slow) (Fig. 3f). All methods performed worse as the decay time constant increased. This occurred because the slower the calcium indicator is, the less the resulting fluorescence signal resembles the original spike train, which increases the difficulty in disambiguating between the latent time courses and the calcium decay (Fig. 3f).

Finally, we varied the amount of noise added to the fluorescence, which reflects imaging noise independent of calcium and spiking activity (Fig. 3g). As the variance of the noise increased, all three methods performed worse as expected, although the extent of the performance degradation differed across methods. As the fluorescence noise variance increased, CILDS continued to outperform the other two methods (Fig. 3g). This indicates that leveraging the population of neurons for simultaneous deconvolution and dimensionality reduction, as done by CILDS, provides statistical power to mitigate a loss in accuracy due to increased fluorescence noise (Fig. 2b).

Overall, CILDS performed as well or better than the other two methods in every simulated setting we tested (Fig. 3d–g, orange higher than purple and cyan, also see Supplementary Fig. 1 for additional combinations of simulation parameter settings and Supplementary Fig. 2 for how the results vary with mean firing rates). We additionally found that deconv-LDS usually outperformed LDS in accuracy of recovered latent variables, consistent with Wei et al., which applied PCA to trial-averaged activity [27] (Fig. 3d–g, purple higher than cyan). Using deconvolution is particularly important in regimes where the time scales of neural activity (i.e., the latent timescales) are faster than that of the calcium decay, which is the case for many commonly-studied brain functions.

Need for latent dynamical model in dimensionality reduction

All three dimensionality reduction methods considered so far explicitly attempt to extract latent variables that evolve smoothly over time via a latent dynamical model. We asked if including a latent dynamical model was necessary for the accurate recovery of latent variables.

To address this, we compared three methods: 1) our method that incorporates deconvolution and latent dynamics (CILDS), 2) a method we created that incorporates deconvolution but no latent dynamics (Calcium Imaging Factor Analysis, CIFA, see Methods), and 3) an off-the-shelf method that does not incorporate deconvolution or latent dynamics (Factor Analysis, FA). CIFA is identical to CILDS, except that CIFA does not have a latent dynamical model, as with conventional FA. Across a range of latent timescales, we measured the ability of each method to accurately recover the ground truth latent variables (Fig. 4a).

We found that when the latent timescales are fast, the calcium transient blurs the neural activity and deconvolution is necessary for more accurate latent variable recovery (Fig. 4a, left, consistent with Fig. 3d). This is why the methods that involve deconvolution (CILDS, CIFA) outperform the method that does not involve deconvolution (FA). CILDS outperforms CIFA in this regime because CILDS additionally smooths the latent variables over time via the latent dynamical model in a way that is determined by the data (see Methods).

When the latent timescales are slow, it becomes important that the latent variables are smoothed temporally (Fig. 4a, right), either through the use of a latent dynamical model (CILDS) or by retaining the smoothing provided by the calcium indicator dynamics when they are not deconvolved away (FA). These methods outperform the method for which the latent variables do not have temporal smoothing (CIFA). Furthermore, the fact that CIFA does not have a latent dynamical model influences its estimate of the calcium decay constant. CIFA can attribute temporal smoothness in the fluorescence to only one possible source – calcium decay. For this reason, as the simulated latent timescale increases (left to right in Fig. 4b, Supplementary Fig. 3), CIFA erroneously attributes the slower varying fluorescence to a slower calcium decay, rather than to a longer latent timescale. By contrast, CILDS can attribute smoothness in the fluorescence to two possible sources – latent variables that evolve smoothly over time and calcium decay. As a result, as the simulated latent timescales increase (left to right in Fig. 4b, Supplementary Fig. 3), CILDS correctly attributes the slower varying fluorescence to the latent variables, and not to a longer calcium decay time constant.

In summary, we found that rather than a latent dynamical model being necessary per se, temporal smoothing of the latent variables is necessary for their accurate recovery, to a degree that is appropriate for the underlying latent or neural process. The trends that we observed for CILDS, CIFA, and FA are true across a variety of related methods (Supplementary Fig. 4). Overall, CILDS performs well in all regimes, regardless of whether the latent timescales are fast or slow. For this reason, we consider CILDS to be an excellent tool for exploratory data analysis for which we do not know a priori whether the underlying neural processes are fast or slow.

Application to calcium imaging recordings

To assess whether the advantages of CILDS also hold in real data, we applied each of the dimensionality reduction methods described above (CILDS, deconv-LDS, LDS, and CIFA) to calcium imaging recordings in two experimental contexts: larval zebrafish and mice. To emphasize the generality of our findings, these two experimental settings involve not only different animal species, but also different brain areas, behavioral tasks, and properties of the recorded fluorescence (see below). In these experiments, the ground truth latent variables are unknown. To quantify the accuracy of each method, we adopted a leave-neuron-out procedure used in previous studies [38, 39]. We estimate the latent variables using all-but-one neuron, and then assess how well these latent variables predict the recorded fluorescence of the held-out neuron (see Methods). This procedure assesses the relative ability of each method to identify a set of latent variables that captures the variability that is shared amongst the neurons. We used between 3 to 17 latent variables, depending on the dimensionality reduction method, for the zebrafish and 50 latent variables for the mice (see Methods).

The first experimental context involves larval zebrafish engaged in a “fictive swimming” moto-sensory gain adaptation task (Fig. 5a) [40]. Calcium imaging was performed on neurons expressing GCaMP6f in dorsal raphe nucleus (DRN) using light-sheet microscopy in three fish (19 to 22 neurons; mean: 20) (Fig. 5a). We applied each dimensionality reduction method to these recordings and assessed their performance using the leave-neuron-out prediction procedure (Fig. 5b). We found that CILDS more accurately predicted the fluorescence of the held-out neurons than the other methods (Fig. 5c–e), as quantified by the correlation between the predicted and recorded fluorescence. Similar results were obtained when examining each fish individually (Supplementary Fig. 5), as well as for a larval zebrafish experiment in which hundreds of neurons across multiple brain areas were analyzed together (300 neurons) (Supplementary Fig. 6).

The second experimental context involves awake, head-fixed mice passively viewing static visual gratings (Fig. 6a) [41]. Two-photon calcium imaging was performed on neurons expressing GCaMP6f in the primary visual cortex (V1) of three mice (133 to 319 neurons; mean: 234.7). Comparing the raw recorded fluorescence of the two experimental contexts (Fig. 5a versus Fig. 6a), the neurons in the DRN of the larval zebrafish tend to exhibit slower fluctuations in fluorescence that are more correlated across neurons (Supplementary Fig. 7). In contrast, the neurons in mouse V1 exhibit faster changes in fluorescence, which are less correlated across neurons. We applied the same leave-neuron-out fluorescence prediction analysis from Fig. 5 to these mouse recordings (Fig. 6b). Despite the stark differences between the fish DRN and mouse V1 fluorescence traces, we again found that CILDS more accurately predicted the fluorescence of held-out neurons than LDS (Fig. 6c), deconv-LDS (Fig. 6d), and CIFA (Fig. 6e). These results were also true for each mouse individually (Supplementary Fig. 8).

Thus far, we have shown that CILDS extracts latent variables that more accurately predict the fluorescence of held-out neurons than the other methods. Another way to assess how meaningful the latent variables extracted by each method are is to measure how strongly they reflect external variables [42]. We thus performed a decoding analysis, whereby we classified the orientation and spatial frequency of the presented grating using the latent variables extracted by each of three methods (Fig. 6f). We found that CILDS and CIFA performed similarly, and both outperformed deconv-LDS (by 1.16 times) and LDS (by 1.43 times) (Fig. 6g). This demonstrates that the joint methods are better at capturing the shared activity changes among neurons that are relevant to the visual stimulus than the other methods. To understand why CILDS and CIFA performed similarly, recall from the simulations that the performance between CILDS and CIFA becomes more similar as latent timescales get faster (Fig. 4a). The latent timescales of these mouse recordings are at least as fast as the rate of change in the visual stimulus (every 250ms), as seen by the abrupt transients in the raw fluorescence traces (Fig. 6a). Furthermore, the higher decoding accuracy for deconv-LDS compared to LDS indicates that accounting for the calcium transients is important when the variable of interest changes on the timescale of tens to hundreds of milliseconds.

Taken together, the results based on calcium imaging recordings from two different recording regimes are consistent with what we identified in simulation. Namely, deconvolution should be used jointly with a dimensionality reduction method that provides temporal smoothing (as in CILDS), particularly if the neural process of interest changes on a timescale of tens to hundreds of milliseconds.

Dimensionality reduction methods

Here we describe mathematically the dimensionality reduction methods used in this work: LDS, deconv-LDS, CILDS, and CIFA. For the purposes of this work, we assume that the spatial footprint of each neuron has already been identified from the raw calcium imaging data (a procedure known as image segmentation[62, 55]), resulting in a fluorescence time course for each neuron. The dimensionality reduction methods presented here are applied to these fluorescence time courses.

Linear Dynamical System (LDS)

We first considered a standard dimensionality reduction method for summarizing the time course of spiking activity, a latent Linear Dynamical System (LDS), here applied to calcium imaging recordings. Let ${\mathit{y}}_t\in {ℝ}^{q\times 1}$ be a high-dimensional vector of fluorescence values recorded at time point t, where q is the number of neurons imaged simultaneously. The goal is to extract a corresponding low-dimensional latent variable ${\mathit{z}}_t\in {ℝ}^{p\times 1}$ at each time point, where p is the number of latent dimensions $(p<q)$. The observation model is

where $A\in {ℝ}^{q\times q}$ is the loading matrix that specifies how each neuron’s activity is related to the latent variables, $\mathit{b}\in {ℝ}^{q\times 1}$ is an offset vector that accounts for constant background fluorescence, $R\in {ℝ}^{q\times q}$ is the observation noise covariance, and $t=1,\dots ,T$. We constrained R to be diagonal, thereby capturing activity variability and imaging noise independent to each neuron. The time-evolution of the latent variables is described as a linear dynamical system

where $D\in {ℝ}^{p\times p}$ is the dynamics matrix that determines the timescale of the latent variables, $P\in {ℝ}^{p\times p}$ is the dynamics noise covariance, ${\mathit{h}}_1\in {ℝ}^{p\times 1}$ and $G_1\in {ℝ}^{p\times p}$ describe the mean and covariance of the latent variable at the first time point, and $t=2,\dots ,T$. We constrained D, P, and G₁ here to be diagonal as a form of regularization, although a general LDS with these parameters unconstrained could also be used.

Equations (1), (2), and (3) together define the latent LDS. We fit the model parameters (A, b, R, D, P, h₁, G₁) using the expectation-maximization (EM) algorithm. To initialize the model parameters, we first performed factor analysis (FA) on y_t to obtain A, b, and R. We initialized D to be 0.999I (a stable system), which we found to work well in practice in the simulations. We ran the EM algorithm until convergence (defined as a log data likelihood increase of < 10⁻⁶ or 1500 iterations, whichever came first). This maximum number of iterations was chosen heuristically, by noting empirically that the latent variables do not change significantly beyond this point.

Deconvolution - Linear Dynamical System (deconv-LDS)

Since fluorescence traces are an indirect measure of spiking activity, we also considered a two-stage approach, whereby we first deconvolve the fluorescence traces one neuron at a time to estimate the underlying spiking activity. We then applied a standard dimensionality reduction method, in this case LDS, to those deconvolved estimates. We refer to this two-stage method as deconv-LDS.

For the deconvolution stage of deconv-LDS, we used the Online Active Set method to Infer Spikes (OASIS), developed by Friedrich et al. [34], using their L₁ regularization. We also tested their L₀ regularization and found the L₁ regularization to work better for our datasets. As per Friedrich et al., this first order autoregressive model for OASIS is described for each neuron as

where $y_t$ is the recorded fluorescence at at time $t,c_t$ represents the calcium concentration at time $t,{ϵ}_t$ captures imaging noise independent of the calcium and spiking activity, and s_t is the spiking activity. The parameter a relates the calcium concentration to fluorescence, b accounts for the baseline fluorescence, ${\sigma }^2$ captures the variance of the imaging noise, and γ specifies how quickly the calcium trace decays, which depends on the calcium indicator. Additionally, there is a hyperparameter in the OASIS model, minimum spike size, that sets the minimum value of s_t that would be identified. In this model, all variables are scalars and a is constrained to be non-negative. We initialized OASIS with γ values that are typical for the calcium indicators used [13] (see Supplementary Table 1), and allowed OASIS to optimize $a,b,\gamma ,{\sigma }^2$ and the minimum spike size.

Applying deconvolution to the recorded fluorescence traces returned estimates of the time course of spiking activity for each neuron, $s_t$. We then used the estimated spiking activity of all the neurons as the observations ${\mathit{y}}_t\in {ℝ}^{q\times 1}$ for time points $t=1,\dots ,T$ in the LDS model defined in equations (1)–(3).

Calcium Imaging Linear Dynamical System (CILDS)

CILDS unifies the approaches described above by allowing estimates of shared activity among neurons (i.e., the latent variables) to influence the estimates of deconvolved spiking activity, and vice versa. In other words CILDS performs deconvolution for all neurons and dimensionality reduction jointly, in a unified framework. This is in contrast to deconv-LDS, which deconvolves the activity of each neuron independently. With low-dimensional latent variables that are jointly estimated with the model of calcium decay, CILDS is better able to peer through the calcium decay to more clearly identify the shared activity among neurons, as compared to deconv-LDS and LDS applied directly on fluorescence.

Let ${\mathit{y}}_t\in {ℝ}^{q\times 1}$ be the high-dimensional vector of fluorescence traces recorded at each time point t, where q is the number of neurons imaged simultaneously. The goal is to extract a corresponding low-dimensional latent variable ${\mathit{z}}_t\in {ℝ}^{p\times 1}$ at each time point, where p is the number of latent dimensions $(p<q)$. The observation model follows the multivariate form of equation (4) which was used for deconvolution

where $B\in {ℝ}^{q\times q}$ maps the calcium concentration to the recorded fluorescence, $R\in {ℝ}^{q\times q}$ is the fluorescence noise covariance, and $t=1,\dots ,T$. We constrained B and R to be diagonal to allow each dimension of c_t to represent the calcium concentration of one neuron. B accounts for all experimental variables influencing the scale of the signal from each neuron, such as the amplification of the imaging system [32], and R accounts for fluorescence fluctuations independent of calcium concentration. Here we omit the additive offset found in equation (4) without loss of generality due to the offset included in equation (7).

Similar to equation (5), the calcium decay for each neuron is described using a first order autoregressive model

where $\text{Γ}\in {ℝ}^{q\times q}$ captures the calcium decay, $A\in {ℝ}^{q\times p}$ is the loading matrix that describes how latent variable maps to calcium concentrations, $\mathit{b}\in {ℝ}^{q\times 1}$ is a constant vector, $Q\in {ℝ}^{q\times q}$ captures the spiking variability independent to each neuron, ${\mathit{\mu }}_1\in {ℝ}^{q\times 1}$ and $V_1\in {ℝ}^{q\times q}$ describe the mean and variance of the calcium concentration at the first time point, respectively, and $t=2,\dots ,T$. Our key innovation is to replace the spiking activity s_t from equation (5) with $A{\mathit{z}}_t+\mathit{b}+{\mathit{w}}_t$. By analogy to factor analysis, $A{\mathit{z}}_t$ describes the shared activity among neurons, and ${\mathit{w}}_t$ describes the activity independent to each neuron. We constrain $\text{Γ},Q$, and $V_1$ to be diagonal to prevent intermixing among neurons outside of the loading matrix. The form of Γ allows each neuron to have a different calcium decay as determined by the fitting procedure (see below), which can depend on the extent of calcium buffering within a cell and the calcium indicator used [32, 13, 34].

Similar to LDS, the low-dimensional latent variables z_t evolve over time according to a linear dynamical system

where $D\in {ℝ}^{p\times p}$ is the dynamics matrix, $P\in {ℝ}^{p\times p}$ is the dynamics noise covariance, ${\mathit{h}}_2\in {ℝ}^{p\times 1}$ and $G_2\in {ℝ}^{p\times p}$ are the mean and covariance of the latent variable at the first time point respectively, and $t=3,\dots ,T$. We constrained D, P, and $G_2$ here to be diagonal as a form of regularization, although a model with these parameters unconstrained could also be used. Note that according to equation (7), ${\mathit{z}}_2$ is the first latent variable in the time series (i.e., there is no ${\mathit{z}}_1$).

Equations (6) – (10) define CILDS. The joint estimation of the parameters B, R, Γ, A, b, Q, μ₁, V₁, D, P, h₂, and G₂ allows CILDS to leverage the entire recorded neural population to perform deconvolution and estimate latent variables in a unified fashion. CILDS can be viewed as a special case of the standard LDS, where the parameters are constrained in a specific way. We fit CILDS using the EM algorithm, initialized using deconv-LDS run for 100 EM iterations. The EM algorithm was run until convergence, defined as a log data likelihood increase of < 10⁻⁶ or 1500 iterations, whichever came first. The maximum number of iterations was chosen by noting empirically that the latent variables do not change substantially beyond this point. The EM equations for CILDS are provided in Supplemental Information.

Calcium Imaging Factor Analysis (CIFA)

To evaluate the importance of incorporating a latent dynamical system in dimensionality reduction methods for calcium imaging, we developed a method (CIFA) identical to CILDS, with the only difference being that CIFA does not enforce latent dynamics. Like CILDS, CIFA also uses equations (6), (7), and (8). The only difference between CIFA and CILDS is that we replaced the latent dynamical system equations (9) and (10) with

for $t=2,\dots ,T$. In other words, CIFA defines latent variables that are independent across time, whereas CILDS defines latent variables that evolve smoothly over time. Like CILDS, we fit CIFA using the EM algorithm, initialized using parameters from deconv-FA run for 100 EM iterations. The EM algorithm was run until convergence, defined as a log data likelihood increase of < 10⁻⁶ or 1500 iterations, whichever came first. The EM equations for CIFA are the same as for CILDS (see Supplemental Information), without the equations estimating D, P, h₂, G₂. Note that there is no loss of generality by setting the prior distribution of ${\mathit{z}}_t$ to $𝒩(0,I)$, compared to a general Gaussian distribution. Additionally, although CIFA has FA in its name, there is one key difference from FA. Whereas FA can be fit on data with no concept of time, CIFA requires a time series for deconvolution (equations (7) and (8)).

Simulation framework

We created a framework to simulate fluorescence recordings from calcium imaging for two reasons. First, in calcium imaging recordings, the ground truth latent variables are unknown. A simulation of fluorescence traces from known latent variables allows us to directly evaluate our dimensionality reduction methods by comparing the estimated latent variables with the ground truth latent variables. Second, we wanted a simulation framework in which we could systematically vary various experimentally relevant parameters to see their effects on the estimated latent variables. Specifically, we evaluated our ability to recover the ground truth latent variables as a function of the timescale of the latent variables, the number of neurons, the calcium decay rate, and the size of the fluorescence noise. The simulation procedure consists of first generating fluorescence traces from known latent variables while varying the experimentally relevant parameters listed above. Then, we applied each dimensionality reduction method to estimate latent variables from the simulated fluorescence traces. The estimated latent variables were then compared to the ground truth latent variables. The steps of this procedure are detailed below.

Generating fluorescence traces

To simulate fluorescence traces, we first drew p Gaussian processes (GP), where each GP has T time points at 1 ms time resolution. We denote the ith GP as ${\mathit{z}}_{i,:}\in {ℝ}^{1\times T}$, where $i=1,\dots ,p$. The GP allows us to specify the covariance $K_i\in {ℝ}^{T\times T}$ for the ith GP across the T time points as

and $t_1,t_2=1,\dots ,T$. Here we chose the commonly used squared exponential covariance function. The squared exponential covariance is defined by its signal variance ${\sigma }_{f,i}^2\in {ℝ}_{+}$, characteristic timescale ${\tau }_i\in {ℝ}_{+}$, and noise variance ${\sigma }_{n,i}^2\in {ℝ}_{+}$. The Kroneker delta ${\delta }_{t_1,t_2}$ equals 1 if $t_1=t_2$ and 0 otherwise. We set ${\sigma }_{f,i}^2=1-{\sigma }_{n,i}^2$ so that the latent variable ${\mathit{z}}_t\in {ℝ}^{p\times 1}$ at every time point has mean 0 and a variance I. The vector ${\mathit{z}}_t$ comprises the tth time point from each of the p Gaussian processes. The noise variance ${\sigma }_{n,i}^2$ must be nonzero to ensure that $K_i$ is invertible, hence for practical purposes we set ${\sigma }_{n,i}^2={10}^{-9}$. We note that an LDS can also be used to introduce latent dynamics, but we chose to use a GP due to the ease by which we can generate stationary time series with specified time scales. Furthermore, using a GP introduced model mismatch for all of the dimensionality reduction methods, which enables a more meaningful comparison across methods.

Next, we projected the low-dimensional latent variables ${\mathit{z}}_t$ into the high-dimensional neural space to obtain neural firing rates for each of q neurons at each time point $t=1,\dots ,T$. We then imposed a rectifying nonlinearity applied element-by-element using $\log (1+\exp (W{\mathit{z}}_t+\mathit{\mu }))$, where $W\in {ℝ}^{q\times p}$ is the loading matrix and $\mathit{\mu }\in {ℝ}^{q\times 1}$ is a constant offset, to ensure that firing rates are non-negative. We generated binary spikes ${\mathit{s}}_t\in {ℝ}^{q\times 1}$ at each time point using an inhomogeneous Poisson process with time-varying rates defined by the output of this rectifying nonlinearity. This allows us to assess how the dimensionality reduction methods, which assume Gaussian noise for mathematical tractability, would perform on realistic spiking variability that is not Gaussian. Finally, we obtained the calcium concentration ${\mathit{c}}_t\in {ℝ}^{q\times 1}$ using a first order autoregressive model and fluorescence ${\mathit{y}}_t\in {ℝ}^{q\times 1}$, as in Friedrich et al. 2017 [34]

where $\text{Γ}\in {ℝ}^{q\times q}$ determines how quickly the calcium decays, $B\in {ℝ}^{q\times q}$ relates the calcium concentration to the fluorescence, $\mathit{b}\in {ℝ}^{q\times 1}$ is the baseline fluorescence, $R\in {ℝ}^{q\times q}$ is the imaging noise covariance, and $t=1,\dots ,T$, at the same (1ms) resolution as the GP. We set Γ, B, and R to be diagonal. We specified Γ such that the decay constants approximately matched the decay constants of GCaMP6f, GCaMP6m, and GCaMP6s, found in Chen et al. [13] (See Supplementary Table 1). B represents experimental variables influencing the scale of the calcium signal of each neuron, such as the amplification of the imaging system. B is set as the identity matrix and b is set to be 0 in our simulations. We varied the signal-to-noise ratio by varying R to simulate low, medium, and high noise regimes (Supplementary Table 1).

For each simulation run, we simulated p = 10 latent variables, where each latent dimension had the same latent timescale for ease of interpretation. Across different settings of experimentally relevant parameters, we explored a range of timescales $\tau \in \{50,100,200,1000,2000,5000\}\text{ms}$. To define loading matrices W that were realistic for spiking activity, we used electrophysiological recordings with 94 neurons [6]. Across runs, we tested $q\in \{20,50,94\}$ neurons. For the q = 20 and q = 50 cases, we randomly subsampled the electrophysiological recordings. To avoid a start-up transient at the start of every trial from the calcium model, we generated two long fluorescence traces for each neuron, each 6,000,000 time points long (1 ms resolution). One fluorescence trace was used for training, and the other for testing. We divided each trace into 100 trials. Each trial was 60 s long, comprising 60,000 time points. Finally, having generated these fluorescence traces at 1 kHz, we down-sampled the fluorescence rate to a more typical imaging rate of 40 Hz (2,400 time points per trial). The simulation parameters are summarized in Supplementary Table 1. See Practical Considerations below for how we generated GPs with 6, 000, 000 time points within the computer’s memory constraints.

Evaluation of dimensionality reduction methods

We applied the four dimensionality reduction methods to the simulated fluorescence traces to evaluate how accurately each method reconstructed the ground truth latent variables. Rather than simulating fluorescence traces directly from the generative models of the dimensionality reduction methods, we simulated fluorescence traces in a way that introduced model mismatch, which is present when applying any statistical method to real-world data. We first generated the latent variables using Gaussian processes, rather than using an LDS. Then, we simulated spike trains based on the latent variables using an inhomogeneous Poisson process, rather than adding Gaussian noise to the firing rates. We believe the presence of model mismatch strengthens our results with simulated data.

We used two-fold cross validation in estimate the accuracy of latent variable recovery, so that the model parameters were fit using half of the data (100 training trials, each with 2,400 time points) and those parameters were then used to estimate the latent variables in the other half of the data (100 held-out trials). For all methods, the latent variables are only unique up to an arbitrary linear transformation. In order to compare the estimated and ground truth latent variables, we aligned them using the following procedure. First we split the estimated latent variables for each test fold into two further inner halves. We concatenated the trials over time, such that the estimated latent variables are defined as $\tilde{Z}\in {ℝ}^{p\times \tilde{T}}$, and the ground truth latent variables are defined as $Z\in {ℝ}^{p\times \tilde{T}}$, where p = 10 is the dimensionality of the latent variables and $\tilde{T}=120,000$ for one inner half. We then applied linear regression to relate $\tilde{Z}$ and the Z corresponding to that inner half. This yielded a transformation matrix W, where $W=(Z{\tilde{Z}}^{\prime }){(\tilde{Z}{\tilde{Z}}^{\prime })}^{-1}$. Finally we used W to align the $\tilde{Z}$ corresponding to the other inner half to obtain the transformed estimated latent variables $\widehat{Z}=W\tilde{Z}$, where $\widehat{Z}\in {ℝ}^{p\times \tilde{T}}$. We performed the same procedure to align each inner half of each cross-validation fold. In Figs. 3 and 4, we provide each dimensionality reduction method with the number of ground truth latent variables (p = 10) to be extracted. We also performed an analysis where we varied the number of extracted latent variables for a fixed number of ground truth latent variables (Supplementary Fig. 10).

We computed the accuracy of the ith transformed estimated latent variables

where $z_t^{(i)}\in ℝ$ is the $(i,t)$ entry of Z, and ${\widehat{z}}_t^{(i)}\in ℝ$ is the $(i,t)$ entry of $\widehat{Z}$. The mean of the ith ground truth latent variable across time is defined as ${\overline{z}}^{(i)}\in ℝ$, where $i=1,\dots ,p$. A larger R² means a better match between the estimated latent variables and the ground truth latent variables, based on the proportion of total variance (of the ground truth latent variables) explained. R² has an upper limit of 1, and any value below 0 indicates that the estimate is poorer than the using the ground truth mean. We repeated this process for each cross-validation fold and averaged the results across all p latent dimensions and folds. It is important to note that since the results are cross-validated, a dimensionality reduction method with more parameters will not necessarily outperform a method with fewer parameters.

Larval Zebrafish

Neurons were imaged from the dorsal raphe nucleus (DRN) of larval zebrafish while they engaged in a “fictive swimming” motosensory gain adaptation task [40]. Calcium imaging was performed using light-sheet microscopy at 30 Hz on a single plane of narrow area around the DRN. This was performed with Tg(elavl3:GCaMP6f)^jf1 fish expressing GCaMP6f in the cytosol [63]. We analyzed baseline corrected fluorescence recordings from three fish, the same data used in Kawashima et al. 2016[40]. In the task, the fish underwent an initialization period of 20 seconds to increase locomotor drive, a training period in which the fish attenuated their locomotor drive, and a delay period of 10 seconds which stopped the fish from swimming. Finally, there was a test period of 5 seconds to probe the extent to which the attenuated locomotor drive persisted throughout the delay period. There were three different training period lengths of 7, 15, or 30 seconds. Here we combined the data from the different training periods. We included only the 20-second initialization period and the first 7 seconds of the training period of each trial. Thus, the analyzed portion of each trial is nominally identical. For each fish, this yielded 15 trials, each 27 seconds long. We analysed 22, 19, and 19 neurons imaged from the DRN of each fish, respectively. All experiments presented in this study were conducted according to the animal research guidelines from NIH and were approved by the Institutional Animal Care and Use Committee and Institutional Biosafety Committee of Janelia Research Campus.

Mouse

Two-photon calcium imaging was performed in the binocular zone of V1 in awake head-fixed mice resting atop a floating spherical treadmill [41]. GCaMP6f was expressed in excitatory neurons and imaged at 15.5 Hz. We analyzed recordings from three mice (2 male, 1 female, age 42-80 days). These were homozygous Emx1cre mice (Jackson Laboratories, stock number 005628) crossed with homozygous Ai93/heterozygous Camk2a-tTA mice (Jackson Laboratories, stock number 024108). Experimental mice were heterozygous for all three alleles. Mice were positioned to passively view static sinusoidal gratings, without reward. There were 180 gratings presented, comprising 12 different orientations equally spaced with range {0 – 165}° and 15 different spatial frequencies equally spaced with range {0.02 – 0.30} cycles/°. Each “trial” was a 196.7 seconds long recording (3049 time points), comprising 4 presentations of each of the 180 possible gratings in random order. Each presentation lasted 250 ms without an intervening grey screen. The onset time of the first stimulus relative to the beginning of the recording was varied. This means that there is a short period of time recorded before the first stimulus is shown, and a short period of time recorded after the last stimulus is shown. The experiment comprised 15 trials for each mouse. We analysed 133, 252, and 319 neurons from V1 of mouse (labeled mouse 1-3, respectively). These mice correspond to mouse2317, mouse2320, and mouse2209 in the experiments. All experimental procedures were compliant with the guidelines established by the Institutional Animal Care and Use Committee of Carnegie Mellon University and the National Institute of Health, and all experimental protocols were approved by the Institutional Animal Care and Use Committee of Carnegie Mellon University.

Leave-neuron-out fluorescence prediction

We sought to compare the four dimensionality reduction methods using experimental data. In experimental data, ground truth latent variables are unknown, and so we could not use the same evaluation procedure of comparing estimated and ground truth latent variables as in the simulations. Furthermore, the cross-validated data likelihoods are not comparable across all methods. Hence to compare the four methods, we performed a leave-neuron-out fluorescence prediction test to determine which method best summarizes the neuronal activity with low-dimensional latent variables [38, 39]. The intuition is that a method that provides a better summary of the population activity using the latent variables would be better able to predict the activity of held-out neuronal fluorescence traces.

For the leave-neuron-out fluorescence prediction, we performed 5-fold cross-validation. We first split the trials into five equal-sized folds. For a given dimensionality reduction method, we fit the model parameters using four of the folds. With the remaining validation fold, we estimated the latent variables using all but one neuron, and then predicted the activity of that held-out neuron using those estimated latent variables. We did this for every neuron in the validation fold. In the same manner, we performed leave-neuron-out fluorescence predictions for the remaining folds and thus obtained predictions of the fluorescence activity for all the neurons at all time points.

After predicting the leave-neuron-out fluorescence for each neuron using each dimensionality reduction method, we computed the Pearson’s correlation coefficient between the predicted fluorescence and the ground truth fluorescence for every neuron and every trial. For each neuron, we then use the correlation values obtained for all the trials to apply a paired two-tailed t-test (across 15 values per neuron for the larval zebrafish DRN recordings, and 10 values per neuron for the mouse V1 recordings). We indicated the neurons for which the difference in leave-neuron-out fluorescence prediction correlation values between CILDS and each of the other methods was statistically significant (p < 0.05; Figs. 5, 6, Supplementary Figs. 5, 6, 8). For the fish (Fig. 5), we found that CILDS outperforms LDS (59% of red points above the diagonal), deconv-LDS (100% of red points above the diagonal), and CIFA (93% of red points above the diagonal). For the mice (Fig 6), we found that CILDS outperforms LDS (68% of red points above the diagonal), deconv-LDS (84% of red points above the diagonal), and CIFA (99% of red points above the diagonal).

Performing the leave-neuron-out prediction procedure requires selecting the latent dimensionality p for each method. The classic way of finding the optimal latent dimensionality is through nested cross-validation. Nested cross-validation uses, in this case, the 5-fold cross-validation described earlier as the outer folds. Within each outer fold, the training portions are then used to perform an inner 4-fold cross-validation to determine the latent dimensionality that maximizes the cross-validated data likelihood. We used this latent dimensionality to estimate the model parameters using the same four training folds. These model parameters were then used for the leave-neuron-out prediction procedure of the remaining (validation) outer fold. Note that although the cross-validated data likelihood is not comparable across methods, it is comparable across different latent dimensionalities for the same method. To illustrate the model selection process, we show how the cross-validated data likelihood varies with dimensionality for both fish and mice, which is representative of how the cross-validation procedure is carried out for each outer fold (Supplementary Fig. 11).

For the larval zebrafish DRN recordings, we selected the latent dimensionality for each dimensionality reduction method (LDS: 17, deconv-LDS: 17, CILDS: 7 to 17, CIFA: 3 to 5) using nested cross-validation, as described above. For the mouse V1 recordings, nested cross-validation would have required us to fit the model 25 times, for which the running time would have been prohibitively long. For this reason, we selected dimensionality for the mouse V1 leave-neuron-out analysis by assessing the cross-validated decoding performance. We found that the performance of all dimensionality reduction methods increased with latent dimensionality across the range of dimensionalities tested (5-50). Thus we used a latent dimensionality of 50 for all methods for cross-validated leave-neuron-out prediction in mouse V1.

Decoding analysis

Another way to assess how meaningful the extracted latent variables are is by decoding external variables from the latent variables [42]. For the larval zebrafish recordings, there is no moment by moment behavior that can be decoded from the neural activity. In the mouse recordings, we can decode the orientation and spatial frequency of the grating stimuli. To begin, we applied the dimensionality reduction methods to the mouse recordings. We then applied a linear Gaussian Naïve Bayes classifier [64] to the extracted latent variables

where $\mathit{z}\in {ℝ}^{p\times 1}$ is the latent variable averaged across the time points in a 250 ms window corresponding to a given stimulus $C_k$, k = 1, …, 180. The parameters of the classifier are ${\mathit{\mu }}_k$, the mean of the latent variables corresponding to class k, and $\sum$, the covariance of the latent variables across trials. We constrain $\sum$ to be diagonal and the same across all classes. The parameters ${\mathit{\mu }}_k$ and $\sum$ are fit by maximizing the likelihood of the training data using 5-fold cross-validation.

We then used the parameters found from the training data to classify the held-out data according to

where $\widehat{y}$ is the predicted class label (1,…,180). We then computed the accuracy of the predicted class labels against the true class labels (chance level is $\frac{1}{180}$). A higher accuracy indicates that the latent variables are better at capturing the shared modulations among neurons that are relevant to the visual stimulus. For latent dimensionalities {5, 10, 20, 30, 40, 50}, we fit each dimensionality reduction method using all data. Then, we performed 5-fold cross-validation in the decoding stage.

To decode visual stimuli, there are two considerations. First, visual information takes time to arrive in the visual cortex, hence there is a need to shift the window of neural activity relative to the stimulus presentation. Second, there is an additional time delay introduced by the calcium indicators. Unlike deconv-LDS and CILDS, LDS does not attempt to remove the calcium decay. Hence, we expected a longer latency for LDS than for deconv-LDS and CILDS. To determine the appropriate window, we considered a range of time lags and evaluated their cross-validated classification accuracy. We found that the best cross-validated accuracy was obtained for a 4 time point (260 ms) shift for LDS, and a 3 time point (190 ms) shift for deconv-LDS and CILDS. Thus we used these time lags to report the classification accuracy of our decoding analysis (Fig. 6g).