Granger causality analysis for calcium transients in neuronal networks, challenges and improvements
Abstract
One challenge in neuroscience is to understand how information flows between neurons in vivo to trigger specific behaviors. Granger causality (GC) has been proposed as a simple and effective measure for identifying dynamical interactions. At singlecell resolution however, GC analysis is rarely used compared to directionless correlation analysis. Here, we study the applicability of GC analysis for calcium imaging data in diverse contexts. We first show that despite underlying linearity assumptions, GC analysis successfully retrieves nonlinear interactions in a synthetic network simulating intracellular calcium fluctuations of spiking neurons. We highlight the potential pitfalls of applying GC analysis on real in vivo calcium signals, and offer solutions regarding the choice of GC analysis parameters. We took advantage of calcium imaging datasets from motoneurons in embryonic zebrafish to show how the improved GC can retrieve true underlying information flow. Applied to the network of brainstem neurons of larval zebrafish, our pipeline reveals strong driver neurons in the locus of the mesencephalic locomotor region (MLR), driving target neurons matching expectations from anatomical and physiological studies. Altogether, this practical toolbox can be applied on in vivo population calcium signals to increase the selectivity of GC to infer flow of information across neurons.
Editor's evaluation
This important paper provides an indepth analysis of the advantages and potential pitfalls of the application of Granger Causality (GC) to calcium imaging data, especially regarding various types of preprocessing. The authors' approach uses compelling rigor in arguing their points, and it is very clear how one would go about replicating their work. These results should be of interest to any researcher attempting to analyze calcium imaging data.
https://doi.org/10.7554/eLife.81279.sa0Introduction
The prompt integration of sensory inputs by motor command neurons in the nervous systems can become a matter of life and death, often requiring fast and precise movements in response to the stimulus. Understanding how large neuronal networks coordinate to generate cognitive activities requires knowing how information from the stimuli flows between neurons to result in the observed behavior. Recent developments in optogenetics (Deisseroth, 2011; Fenno et al., 2011), and wholebrain imaging with singlecell resolution (Ahrens et al., 2013; Nguyen et al., 2016; Venkatachalam et al., 2016; Cong et al., 2017) allow us to observe collective firing patterns of many neurons and give the basis for reconstructing a causality network that describes how neurons influence each other. However, despite the existing data, assessing the direction of neuronal communication or causality in order to understand how information flows in the brain during integrative reflexive behaviors remains a major challenge in neuroscience (Goulding, 2009).
Apart from perturbing individual neurons and observing the activity of downstream neurons, noninvasive methods for obtaining functional connectivity have been developed (Bastos and Schoffelen, 2016). The widelyused correlation analysis (Cohen and Kohn, 2011; Haesemeyer et al., 2018; Stringer et al., 2019) identifies pairs or groups of neurons that are active at the same time, suggesting either that they drive each other, or are driven by a common input. By contrast, Granger causality (GC) identifies directed interactions (Granger, 1969; Geweke, 1982; Geweke, 1984). Assuming that the two timeseries are well described by Gaussian autoregressive processes, GC quantifies the ability to predict the future values of the time series of a given neuron using prior values of the time series of other neurons. In the field of neuroscience, GC has been used to analyze data from diverse sources (Seth et al., 2015), including electroencephalography (EEG), magnetoencephalography (MEG) (Gow et al., 2008), functional magnetic resonance imaging (fMRI) (Roebroeck et al., 2005; Deshpande et al., 2009; Wen et al., 2013), and local field potentials (LFP) (Brovelli et al., 2004). Despite criticisms such as sensitivity to data preprocessing (Florin et al., 2010), internal dynamics of the neuron (Stokes and Purdon, 2017), additive noise that can be correlated among regions of interest or not (Nalatore et al., 2007; Vinck et al., 2015), and sampling frequency (Zhou et al., 2014; Barnett and Seth, 2017), GC and its extensions (Shojaie and Fox, 2021; Barnett et al., 2018) have been useful to construct functional connectomes, to identify important neurons acting as information hubs, or to classify different brain states (Nicolaou et al., 2012).
While most GC analyses in neuroscience focus on large scale ensemblelevel neuronal signals, it is unclear how well it can be adapted to address singlecell level descriptions, such as spiking (Kim et al., 2011; Gerhard et al., 2013; Sheikhattar et al., 2018) or population calcium imaging data (De Vico Fallani et al., 2014; Severi et al., 2018; Oldfield et al., 2020). Without the law of large numbers, these singlecell time series are no longer approximately Gaussian, the dynamics of signal propagation is increasingly nonlinear, and the sampling rate is typically much slower than the rate of information propagation. All of these points cast doubt on whether GC can correctly identify information flow. Recently however, De Vico Fallani et al., 2014 demonstrated in zebrafish embryos that, despite a low sampling frequency of acquisition for calcium imaging (4 Hz), GC can detect information flow consistent with the known rostrocaudal propagation of activity among motoneurons in the spinal cord. Oldfield et al., 2020 subsequently applied GC to grouped averages of calcium signals over large brain regions without single cell resolution, and detected the information flow correctly based on anatomy from the pretectum to the tectum in the visual system leading to subsequent activity of brainstem and motor circuits during prey capture in larval zebrafish. New variations of GC accounting for sparsity in the connection have also been developed and applied on singlecell level to both spiking and calcium signals in the primary auditory cortex of mice (Sheikhattar et al., 2018; Francis et al., 2018; Francis et al., 2022). Nonetheless, at the singlecell resolution, it is unclear to which extent the application of GC analysis in population calcium imaging is generalizable, and how GC analysis should be improved to accommodate intrinsic features of calcium signals.
In this manuscript, we identify a set of issues inherent to GC analysis when applied to population calcium imaging and propose solutions to address them. We first use synthetic neural signals and then exploit as examples population recordings of in vivo calcium signals in zebrafish with either correlated and noncorrelated noise. We build an improved analysis pipeline customized for calcium signals to address critical issues intrinsic to in vivo calcium imaging, such as the slow exponential decays of calcium transients, correlated noise, and nonlinear and nonGaussian statistics. We develop our method both on a synthetic network of neurons for which we know the interaction structure, as well as on real data from motoneurons in the embryonic zebrafish spinal cord that are recruited along a chain organized from rostral to caudal on either side of the spinal cord (De Vico Fallani et al., 2014). We first improve our analysis pipeline to optimally detect the information flowing between embryonic motoneurons in De Vico Fallani et al., 2014. Second, we apply our improved GC analysis to infer information flow from large neuronal populations in the brainstem of larval zebrafish performing optomotor response (Severi et al., 2018). In the context of large populations of neurons, we show that our method is able to reveal that motorcorrelated neurons within the recently identified mesencephalic locomotor region in CarboTano et al., 2022 drive the activity of downstream neurons in the reticular formation.
Results
Granger causality
Granger causality quantifies the ability to predict the future values of the discrete time series of neuron X of length $T$, $\{{x}_{t}\}$, for $t=1,2,\mathrm{\dots},T$, using prior values of the time series of neuron Y, $\{{y}_{t}\}$, for $t=1,2,\mathrm{\dots},T1$, assuming that the two timeseries are well described by Gaussian autoregressive processes (Granger, 1969; Seth et al., 2015). The method asks whether including information from Y to predict X has significant predictive power, by computing a pvalue assuming as a null model the probability that the dynamics of X can be explained without Y. The prior values of neuron X and Y itself used in that prediction are taken for a finite number of past time points, which defines a characteristic timescale.
With only the information of X itself, the dynamics of X, x_{t}, may be modeled as an autoregressive process, denoted as the ‘reduced model’,
where ${\stackrel{~}{a}}_{q}$ are the regression coefficients, the maximum lag $L$ is a hyperparameter subject to tuning, and ${\stackrel{~}{\u03f5}}_{t}$ is a residual term, assumed to be Gaussian.
Alternatively, one can also write down a linear regression model for x_{t} given the history of both neuron X and neuron Y. This “full model” is
where a_{q} and b_{q} are the regression coefficients for the contribution of neuron X and Y’s history to neuron X’s current state and ${\u03f5}_{t}$ is a residual term, also assumed to be Gaussian.
The key idea of Granger causality is to compare the residuals between the reduced and the full models. Intuitively, if neuron Y can influence X, including the dynamics of Y will improve the prediction of x_{t}, and hence reduce the residual term. As the reduced model is nested within the full model, to test whether the distributions of the residuals, $\{{\epsilon}_{t}\}$ and $\{{\stackrel{~}{\epsilon}}_{t}\}$, are significantly different requires a oneway analysis of variance (ANOVA) test carried out by computing the Ftest statistic:
where ${T}_{\text{regr}}=TL$ is the number of time points used in the regression, and ${M}_{\text{f}}=2L+1$ and ${M}_{\text{r}}=L+1$ are respectively the number of parameters in the full and reduced models. Under the null hypothesis that including y_{t} does not improve the prediction of x_{t}, and that the noise can be described by Gaussian i.i.d., the Fstatistics will follow an Fdistribution with defined degrees of freedom, $\mathcal{F}(F;{M}_{\text{f}}{M}_{\text{r}},{T}_{\text{regr}}{M}_{\text{f}})$. Hence, testing the computed Fstatistic ${F}_{Y\to X}$ against the Fdistribution will distinguish whether Y significantly Grangercauses X.
To measure the amount of reduction in the residual, the bivariate Granger causality value (BVGC) is defined as
where the sample variance requires a correction using the number of samples in the regression (${T}_{\text{regr}}=TL$) minus the number of regression parameters (${M}_{\text{r}}$ or ${M}_{\text{f}}$). Mathematically, the GC value is related to the $F$test statistic through a monotonic transformation:
A significantonly GC value, $G{C}_{\text{Y}\to \text{X}}^{\text{sig}}=G{C}_{\text{Y}\to \text{X}}$ is called if the null test of the Fstatistics is rejected, and zero otherwise.
In multicomponent systems with more than two variables, one needs to consider the influence of other neurons or of a common stimulus, which we collectively denote by $\mathbf{\text{Z}}=\{{\text{Z}}_{1},\dots ,{\text{Z}}_{N}\}$. For example, to analyze a time series of $N+2$ neurons $\{{x}_{t},{y}_{t},{\mathbf{\text{z}}}_{t}\}$, one needs to distinguish direct information flow from Y to X, from indirect information flow from Y to Z and from Z to X, or from a common but delayed drive from Z to Y and from Z to X. Geweke, 1984 extended Granger causality to a multivariate version (MVGC) to measure the conditional dependence between variables. In this case, we explicitly include the variables describing all possible other neurons and external stimuli, $\mathbf{\text{z}}}_{t}=\{{z}_{1,t},\dots ,{z}_{N,t}\$, as regressed variables in both the reduced model and the full model. The current state of neuron X, x_{t}, now also depends on the past of Z, such that the reduced and the full models become
and
where ${\stackrel{~}{c}}_{j,q}$ and ${c}_{j,q}$ are the regression coefficients from the conditioned signal z_{j}. The Ftest statistics $F}_{\text{Y}\to \text{X}\mathbf{\text{Z}}$ and the GC values $G{C}_{\text{Y}\to \text{X}\mathbf{\text{Z}}},$ have identical expression as in Equation 3 and Equation 4, with the number of parameters adjusted accordingly in the reduced and the full models. The significance test is adjusted using multiple comparisons by applying a Bonferroni correction to the pvalue threshold.
It is important to point out that, although GC is based on linear regression of time series, its use is not limited to data generated by linear dynamics. In fact, as long as the data can be approximated by a vector autoregressive process, GC can be useful in revealing the information flow (Seth et al., 2015). For analysis of calcium signals, the intrinsic nonlinearity of calcium decay means it is not obvious a priori whether GC can be useful in revealing the information structure. As we will show in the following, the success of information retrieval depends on the details of the biological network, such as the interaction strength between pairs of neurons and the connection density, which leads to a complicated interplay among different timescales, such as the firing rate, the calcium decay time constant, and the time scale of information transfer.
GC analysis of a small synthetic neuronal net
To gain some intuition about how GC identifies links, we test it on synthetic data generated using small networks of 10 neurons, where we predefine the true connectivity (Figure 1A, B). One of the simplest types of dynamics is given by the same model that underlies the GC analysis, the vector autoregressive model (VAR), where the neural activities ${f}_{i,t}$ evolve according to
where ${\mathrm{\Gamma}}_{ij,q}$ is a memory kernel capturing the effect from neuron $j$ at time $tq$ on neuron $i$ at time $t$, and ${\xi}_{i,t}$ a Gaussian noise with zero mean and variance ${\sigma}^{2}$. For simplicity, we assume the memory kernel acts only within the preceding two time steps, ${L}_{\text{true}}=2$, and set ${\mathrm{\Gamma}}_{ij,q}={A}_{ij}{s}_{j}c$ for $q\u2a7d{L}_{\text{true}}$, and zero, otherwise. ${A}_{ij}$ is the adjacency matrix: ${A}_{ij}=1$ if neuron $j$ is directly presynaptic to neuron $i$, and ${A}_{ij}=0$ otherwise. The sign of interaction s_{j} shows whether neuron $j$ is excitatory (${s}_{j}=1$) or inhibitory (${s}_{j}=1$). The parameter $c$ is the connection strength, which we set to be positive. To ensure that the synthetic time series is stationary, we set the number of excitatory cells and inhibitory cells to be equal. The simulation time step is set as the unit time, ${\mathrm{\Delta}}_{\text{sim}}=1$.
We apply the GC analysis to synthetic data generated using the VAR model (Figure 1C) to see if we can recover the true network connectivity. For each neuron pair $(i,j)$, we compute the BVGC values following Equation 4, identifying ${x}_{t}={f}_{i,t}$ and ${y}_{t}={f}_{j,t}$, and computing the regression coefficients a_{q} and b_{q} using the standard least squares method. Analogously, we also compute the MVGC values, using all the remaining neuron signals ${{f}_{k,t}}_{k\ne i,j}$ as the conditional variables $\mathit{z}}_{t$.
As shown in Figure 1D, MVGC successfully identifies all the links. However, the bivariate GC assigns spurious links, as it cannot distinguish direct links from indirect ones. In general, MVGC has a lower false positive rate (identified nonexisting links) and false negative rate (missing true links) compared to BVGC when we vary the connection strength $c$ (Figure 1I). It is interesting to point out that, at small connection strengths, both MVGC and BVGC will fail to identify true links because the signaltonoise ratio is too small to detect them. In comparison to both GC analyses, the equaltime correlation completely misses the structure (see Figure 1D and Materials and Methods). While crosscorrelation with the lag fixed to the model maximum lag ${L}_{\text{true}}$ performs better, it still fails to reproduce the true connectivity as faithfully as MVGC (see Figure 1—figure supplement 4).
To gain intuition if GC can work for more realistic singleneuron dynamics, which are typically very different from the linear VAR model, we perform the same analysis with dynamics generated with generalized linear models (GLM). We assume the spike count of neuron $i$ at time t, ${\sigma}_{i,t}\sim \text{Poiss}({\lambda}_{i,t})$ is governed by a Poisson process with spiking rate ${\lambda}_{i,t}$.
where ${\mu}_{i}$ is the base rate for spiking, and the memory kernel ${\mathrm{\Gamma}}_{ij,q}$ is identical to the VAR model. We can further convolve the spikes with an exponential decay function to simulate calcium signals (GLMcalcium; see Materials and Methods for details).
Although the GLM and GLMcalcium dynamics are highly nonlinear and nonGaussian (Figure 1E and G), MVGC successfully identifies true links (see Figure 1FGJ). In general, for large correlation strengths, the false negative rate (missing true links) is less than 0.1. Interestingly, for the GLMcalcium model, the false positive rate (identified nonexistant links) at large connection strength is smaller than in the VAR model, as is the difference of errors between the MVGC and BVGC. A possible explanation is that for neurons with nonlinear interactions, the propagation of information in indirect links is less well described by linear models in direct links. In comparison, due to large autocorrelation caused by decay of the calcium signal, the correlation and crosscorrelation performs much worse for the nonlinear dynamics than for the VAR model (see Figure 1—figure supplement 5 and Figure 1—figure supplement 6).
This analysis of synthetic data shows that GC can reveal underlying connectivity structures from calcium signals, and that the performance difference between BVGC and MVGC is smaller in nonlinear models.
In the brain, a population of neurons can encode redundant signals. If the redundancy is in the structure of the network, for example, when two neurons share the identical input and output, the stochasticity will lead to different neuronal dynamics, and MVGC is able to identify the true connection. However, if two neurons have exactly the same activity, MVGC will underestimate the causality, while BVGC still reveal these connections (see Figure 1—figure supplement 7).
In GC analysis, it is essential to choose the maximum time lag $L$. If the maximum time lag $L$ is much less than the true interaction timescale ${L}_{\text{true}}$, then the reconstructed information flow will be incomplete. On the other hand, if $L$ is too large, we will run into the risk of overfitting. Here, we briefly discuss the effect of the maximum lag $L$ on the resulting GC network, using data generated by the synthetic network given by Figure 1A, simulated with VAR dynamics with true maximum lag ${L}_{\text{true}}=3$.
As we vary $L$, the maximum time lag used in the GC regression models, we observe that increasing the maximum time lag does not always improve our ability to identify the network connectivity (Figure 2A, B). Going beyond the true maximum time lag leads to errors, increases both the false positive rate (Figure 2) and the false negative rate (Figure 2B). Compared to MVGC, BVGC is more prone to detect spurious links, even at the correct maximum time lag.
Interestingly, we observe that the average GC value as a function of the number of maximum lags can be used to identify the true maximum time lag ${L}_{\text{true}}$. As shown by Figure 2, both the average BVGC and MVGC values increase steeply as $L$ approaches the true maximum lag, and reach a plateau for a few maximum lag values before overfitting significantly increases the average GC values. Based on this observation, we select the maximum time lag to be the knee of the average GC value as a function of the maximum lag, as this maximum time lag captures enough information flow while the number of fitted parameters is kept as small as possible.
GC analysis of motoneuron data
The data
In order to test the performance of Granger causality on population calcium imaging data, we focus our analyses on previously published calcium transients in motoneurons of embryonic zebrafish (Tg(s1020t:Gal4; UAS:GCaMP3), data published in De Vico Fallani et al., 2014), where the underlying structure of information flow is known. The calcium transients were recorded using fluorescent tags in vivo for 250 s at 4 Hz and with singlecell resolution, in ventral spinal cord motoneurons while the zebrafish was showing spontaneous coiling activity of its tail (Figure 2E–G). During these fictive coiling events, the activity of motoneurons follows a wave from the rostral to the caudal side of the spinal cord (Buss and Drapeau, 2001; Masino and Fetcho, 2005; Warp et al., 2012). Upon activation along the rostrocaudal axis, spinal motoneurons that project to the axial skeletal muscle fibers elicit a rostrocaudal wave of muscle contractions that alternates between the left and right side of the body and propels fish forward. Accordingly, the calcium signals show steady oscillations alternating between the left and right sides of the spinal cord, leading to leftright muscle activation (Figure 2G). We use this data set in order to see if we recover the expected sequential information flow, from the rostral spinal cord to the caudal spinal cord (arrows in Figure 2EF, networks in Figure 2H).
To encode the underlying information structure of the neuronal network, we label the recorded $N$ neurons with index $i$, and the side of the neurons by ${s}_{i}\in \{\text{L},\text{R}\}$. Assume that there are ${N}_{l}$ neurons in the left chain and ${N}_{r}$ neurons in the right chain, we order the indices such that the odd neurons are found on the left side, with ${s}_{i}=\text{L}$ for $i=1,3,\mathrm{\dots},2{N}_{l}1$, and the neurons with $i=2,4,\mathrm{\dots},2{N}_{r}$ are found in the right chain with ${s}_{i}=\text{R}$. The indices are then ordered by the neurons position along each of the two chains, such that if $s}_{i}={s}_{j$, neuron $i$ is more rostral than neuron $j$ if $i<j$. The fluorescence signal for the $i$th cell at time $t$ is denoted by ${f}_{i}(t)$, or equivalently by $D{F}_{i}/{F}_{i}(t)$ to represent the normalized fluorescence activity.
We will compute the GC network for multiple trials of calcium transients recorded in five fishes. Following the nomenclature of the dataset published in De Vico Fallani et al., 2014, we refer to the dataset recorded for fish $A$ and trial $B$ as $fAtB$; for example, $f3t2$ denotes trial number 2 recorded for fish number 3. From the dataset, we exclude recordings for fish number 2 because the recorded time is 2 min, much less than that of the other recordings of 4 min.
Naive GC
We apply the GC algorithm as described in Granger Causality to the calcium transients observed in embryonic motoneurons, and reproduce the results obtained in De Vico Fallani et al., 2014 for the dataset f3t2. We choose the maximum time lag at the knee of average GC value, $L=3$, corresponding to a duration of 750 ms, for both BVGC and MVGC analysis (Figure 2D). This knee corresponds to the maximum lag with the best balance between accuracy and complexity: further increasing the maximum lag yields very similar results but the GC analysis becomes much more computationally expensive (Figure 2—figure supplement 1). In order to verify that the length of the motoneuron dataset is sufficient to estimate the functional connectivity using Granger causality (see Figure 1—figure supplement 2 and Figure 1—figure supplement 3 for illustrations using synthetic data), we compared the two networks obtained by applying GC analysis to each half of the calcium imaging time series. The Pearson correlation coefficient between the GC values from the two halves is 0.87 for BVGC, and 0.75 for MVGC (see Figure 2—figure supplement 2), indicating the trial length is long enough for a consistent estimation of the GC values.
The GC network is then represented by directed graphs (Figure 2I). The nodes of the graph represent the motoneurons and the circle size is proportional to the nodal delta centrality as defined by De Vico Fallani et al., 2014: the larger the circle, the more central is the neuron in terms of its tendency to act as a transmitter hub (red color, positive value) or receiver hub (blue color, negative value) of information flow. Arrows display the functional connectivity links, that is the flow of information from one neuron to another. In the expected network assuming perfect information flow (Figure 2H), all the links are ipsilateral and represent rostrocaudal information flow.
Compared to De Vico Fallani et al., 2014, we identify the same functional connectivity matrices yielding the same connectivity networks. While the previous study focuses on recovering ipsilateral links and ignores contralateral links (depicted in dim gray in the network figures, see De Vico Fallani et al., 2014), we investigate all contralateral links and ipsilateral links (Figure 2I) as we aim to study why spurious contralateral connections occur. We represented the statistically significant direct GC links by arrows, proportional in size and color intensity to the GC value. The networks resulting from the GC analyses both for BVGC and MVGC show many other links than just rostrocaudal ones on either side.
In the MVGC analysis however, one would have expected only direct links to be significant, since conditioning on the other neuronal signals should prevent indirect links from being significant. This expectation is very different from the MVGC network in Figure 2I, where spurious contralateral links prevails.
Nonetheless, it is well known that MVGC is sensitive to additive noise, lack of data points, etc. To assess the impact of these artifacts on the resulting GC networks for motoneuron data, we define measures of information flows as following.
Defining directional biases
The sequential and biaxial information flow in the motoneuron network provides an excellent testing ground for GC. Specifically for the motoneuron datasets, to quantify the success of reconstructing the information flow in the motoneuron network, and to compare the information flow across different motoneuron datasets, we define the following measures of directional biases. First, we define the ratio of normalized total weights of ipsilateral links compared to all links
where ${G}_{ij}$ is the Granger causality values between neurons $i$ and $j$. If all connections are equal, then ${W}_{\text{IC}}=0.5$. We also define the ratio of ipsilateral rostraltocaudal links compared to all ipsilateral links as
This ratio gives a measure of the directionality of information flow. A bias in rostrocaudal links results in a departure from ${W}_{\text{RC}}\sim 0.5$. Notice that based on the anatomy, the ground truth of information flow can be smaller than ${W}_{\text{IC}}=1$, because motoneuron activity has been shown to rely on command neurons in the brainstem that activate motoneurons sequentially in the spinal cord, not on motor neurons synapsing onto each other. With the two measures of directional biases of information flow, now we proceed to develop an improved pipeline for GC analysis. Using different datasets of the embryonic motoneurons as examples, we will show how to improve GC to analyze calcium transients of neuronal populations. Our goal is to recover a more biologically reasonable information flow compared to the naive GC network.
These measures are only constructed for the motoneuron networks. For more general networks with known structures, one can use the distance between the true and the GClearned connectivity matrix to test the efficacy of the GC method. For networks with unknown structure, one can compare the GClearned information flow by using different partitions of the data.
Correcting for motion artifacts
In calcium imaging experiments, small movements of the animal body or the experimental setup can result in artificial fluorescence changes in the recording. In the calcium transients of data f3t1, we observe a big drop in fluorescence at a single time point out of 1000 (Figure 3B top panel) in the calcium traces of all neurons in the recording. As seen in the top panels Figure 3C, this is sufficient to corrupt the results of the GC analyses: neuron 10 (Figure 3A) appears to drive almost all other neurons, including contralateral neurons. We correct the motion artifact by setting the outlying time point’s fluorescence value to the average of the previous and next values (Figure 3B bottom panel). The strongly dominant drive by neuron 10 disappears (Figure 3C bottom panels) and ${W}_{\text{IC}}$ increases from 0.39 to 0.66 for BVGC and from 0.43 to 0.68 for MVGC (Table 1). This shows that GC analysis is sensitive to artifacts at single time points.
The motion artifact is a global noise on the observed neurons, a source that is wellknown to create spurious GC links (Vinck et al., 2015; Nalatore et al., 2007). After careful examination of the calcium traces, we notice that the dominant neuron before motionartifact correction, neuron 10, is characterized by its small dynamic range. Because the rest of the activity is low for neuron 10 and that a single time point seems to have a strong effect on other calcium traces, GC links driven by this neuron are high. We looked at other recordings containing a motion artifact and found that removing it did not always eliminate the spurious dominant neurons. The larger the artifact, the better we eliminate the spurious links: in Figure 3, the artifact represents a deviation of about 80% of the $DF/F$ and its impact on GC links is large. When the deviation is lower, the strong spurious links are not necessarily due to the motion artifact, therefore removing the artifact does not improve the expected network results.
Removing strange neurons
In order to improve the GC analysis results, we want to keep only clean calcium signals and remove the uncharacteristically behaving whose activities do not follow the stereotypical profiles of rise and falls. This is especially important for the MVGC analysis, in which traces of all other neurons are considered in the calculation of a GC value between two given neurons. Although BVGC values are not affected by neurons outside of the observed pair, having extra neurons might also influence the network structure through the GC values of the links.
In general, cells whose calcium transients showed time decay of tens of seconds were identified as non neuronal cells and excluded from GC analysis (see overlays of calcium transients in Figure 3—figure supplement 1). In the embryonic zebrafish traces f5t2, shown in Figure 3E, we remove neuron 13 (Figure 3D) that displays only one long calcium transient lasting for tens of seconds, and neuron 21, which does not exhibit an oscillating pattern.
GC networks are drawn for the neuronal populations before (Figure 3F top panels) and after removal of these outlier neurons (Figure 3F bottom panels). Before the removal, neuron 21 appears to drive many other neurons, likely due to its lower signaltonoise ratio (SNR) and the correlated noise between calcium traces. Once removed, ${W}_{\text{RC}}$ increases, meaning that a mostly rostrocaudal propagation is observed, as information in motoneurons flows from the rostral to the caudal spinal cord (Table 2). We also observe a small decrease in ${W}_{\text{IC}}$ upon removal of neuron 21, due to the fact that the ipsilateral links going out from neuron 21 are not present anymore, but spurious contralateral links remain.
Smoothing
One known limitation of traditional Granger causality analysis is its sensitivity to noise, especially correlated noise across variables (Vinck et al., 2015; Nalatore et al., 2007). In empirical calcium imaging data, correlated noise is common due to systemwise measurement noise and spurious correlations coming from the slow sampling rate.
To illustrate the problem of the noise, we return to the synthetic models in Granger Causality, a simulated system of 10neurons with the VAR and the GLMcalcium dynamics, where the true underlying signal is ${f}_{i}(t)$, for each neuron $i$. We add a systemwide noise, such that the noisecorrupted signal is
and $\zeta (t)$ is Gaussian white noise, with $\u27e8\zeta (t)\u27e9=0$, and $\u27e8\zeta (0)\zeta (t)\u27e9={\sigma}_{\zeta}^{2}\delta (t)$. As shown in Figure 4—figure supplement 2 and Figure 4—figure supplement 3, large correlated noise significantly decreases the ability of the Granger causality analysis to identify true connections.
We exploit the stereotypical shape of calcium signals – a fast onset and a slow decay – to parametrize the structure of the noise for empirical calcium signals. We assume that any correlation in the decay phase is the correlation of the noise. Additionally, we apply total variation differentiation regularization to denoise the fluorescence time series for each neuron (Chen et al., 2019) (see Materials and Methods and Chartrand, 2011 for details). Briefly, for each neuronal signal $f(t)$, total variation differentiation regularization assumes the noise is Gaussian and white, and that the temporal variations in the derivative are exponentially distributed and only weekly correlated in time, which is wellsuited to describe calcium data with episodes of continuous firing and continuous silence. Compared to usual total variation denoising that regulates the temporal variations in the signal, regulating the variation of the derivative is better at detecting sudden change of the signal strength, and is better at keeping the causality when applied to multivariate signals. As an example, Figure 4 shows the original and the smoothed fluorescence signal for two neurons in the same fish (dataset f3t2 from De Vico Fallani et al., 2014), as well as the residual noise. The noise in the motoneuron is correlated with correlation coefficients close to 1 for some pairs of neurons (Figure 4B). The large correlated noise is likely the cause of our problematic results in the MVGC analysis (see Figure 2E and Figure 4C). Denoising the calcium time series before the GC analysis eliminates links, especially in the MVGC approach (Figure 4C). The signal as quantified in terms of the relative weight of ipsilateral links ${W}_{\text{IC}}$ is also closer to the biological expectation (${W}_{\text{IC}}\sim 1$) than before smoothing, with MVGC giving better results (Figure 4D).
We note that total variation regularization is a global denoising method that acts on the entire time series. Different from a lowpass filter, total variation regularization is able to keep the fine structure of the fast rise in the calcium signal, while smoothing out the noise. Nonetheless, the assumption is that the noise is Gaussian, which should be verified before smoothing the data by examining the highfrequency end of the signal’s power spectrum. To address the problem of highfrequency correlated noise, one may propose to use spectral Granger causality (Geweke, 1982; Geweke, 1984), and only focus on the resulting GC values at low frequencies. However, we show that spectral GC applied to the unsmoothed calcium signal is unable to recover the information flow for any of the frequencies (Figure 4figure supplement 4). Finally, in cases where correlated noise cannot be removed, there are alternative methods to reject false links, for example using Granger causality computed on the timereversed data as the null hypothesis (Vinck et al., 2015).
Slow calcium timescale
One potential pitfall of performing Granger causality analysis on singlecell calcium signals is that the timescale of the calcium indicator and the sampling are typically much slower than that of the information transfer. The onset of calcium signals can occur in few milliseconds, comparable with the timescale of information transfer, that is the time difference between the firing of a presynaptic neuron and a postsynaptic neuron downstream. Nonetheless, based on the intrinsic buffer capability of neurons and the fluorescent properties of the calcium indicator, calcium signals typically decay on the timescale of hundreds of milliseconds to seconds (Tian et al., 2012). In motoneurons, neurons can fire in bursts of spikes, which may further hinder the resolution of information propagation. In this section, we investigate how the interplay of these slow and fast timescales impacts the accuracy of Granger causality analysis to identify information flow.
In the embryonic spinal cord, motoneurons are recruited along the rostrocaudal axis sequentially. The muscle contractions alternate between the left and right side of the body, and propel fish forward. Correspondingly, the left and right sides of the motoneurons are active as a steady oscillation from the head to the tail, leading to an expected information flow from rostral to caudal (Warp et al., 2012). The spikes are followed by a plateau that lasts on the order of seconds (SaintAmant and Drapeau, 2001). Using a combined statistical and dynamical approach, we simulate artificial “motoneuron” data (see example traces in Figure 5A), as two chains of five neurons, driven externally by two binary stimuli, ${I}_{\text{L}}(t)$ and ${I}_{\text{R}}(t)$, that randomly flip between an on ($I=1$) and an off ($I=0$) state, following the same statistics as recorded in the data. The stimulus of the left side, ${I}_{\text{L}}(t)$ is generated by randomly choosing time intervals when the stimulus state is on, ${\tau}_{\text{on}}$, and when it is off, ${\tau}_{\text{off}}$, uniformly from the sets $[\mathrm{min}({\tau}_{\text{on}}^{\text{data}}),\mathrm{max}({\tau}_{\text{on}}^{\text{data}})]$ and $[\mathrm{min}({\tau}_{\text{off}}^{\text{data}}),\mathrm{max}({\tau}_{\text{off}}^{\text{data}})]$, respectively. Here, $\{{\tau}_{\text{on}}^{\text{data}}\}$ is the set of duration of all “on”states in the experiment (De Vico Fallani et al., 2014) that we find by the positive finite difference method in the fluorescent signal smoothed by totalvariational regularization (Chartrand, 2011; Chen et al., 2019). These sets of longlasting time blocks, where the neurons are always active, simulates the plateau activity of the embryonic motoneurons. After ${I}_{\text{L}}(t)$ is sampled, the stimulus for the rightside neurons, ${I}_{\text{R}}(t)$, is generated such that the time delay between the “on”state of ${I}_{\text{L}}(t)$ and ${I}_{\text{R}}(t)$ is sampled from $[\mathrm{min}({\tau}_{\text{delay}}^{\text{data}}),\mathrm{max}({\tau}_{\text{delay}}^{\text{data}})]$, where the delay time ${\tau}_{\text{delay}}$ matches the empirical stimulus distribution. We model the sequential activation of the motoneurons by imposing the spike rate as a function of ${I}_{\text{L}}(t)$ and ${I}_{\text{R}}(t)$ (see Materials and Methods) and we use a standard Poisson process with these rates to generate the spike train for each neuron, $\{{\sigma}_{i,t}\}$. The spikes are then convolved with an exponential decay with time scale ${\tau}_{\text{ca}}$ to simulate the slow decay time of calcium signals.
The timescales in this model include the sampling time ${\tau}_{\text{sampling}}$, the calcium decay time constant ${\tau}_{\text{ca}}$, and the timescale of information propagation ${\tau}_{\text{info}}$. The sampling time is given by the experimental setup for the motoneuron data we analyze (De Vico Fallani et al., 2014), with ${\tau}_{\text{sampling}}=0.25\text{s}$. As shown in Figure 3—figure supplement 1, the empirical calcium decay time constant can be found by overlaying all epochs of decays of the calcium signals from different bursts and different neurons: fitting the decay of calcium signals to exponential functions gives ${\tau}_{\text{ca}}=2.5\text{s}$. Notice this decay timescale is limited by both the buffering capability of the sensor, here GCaMP3, and the intrinsic buffering capability of the cell (Tian et al., 2012), that is having a faster sensor does not always decreases this decay time scale. The information propagation timescale is chosen such that $\tau}_{\text{info}}<{\tau}_{\text{sampling}$. The baseline spike rate is set to ${\lambda}_{0}=32{\text{s}}^{1}$. Finally, the time resolution for the simulation is chosen to be ${\mathrm{\Delta}}_{\text{sim}}=0.0125s$, such that $\lambda {\mathrm{\Delta}}_{\text{sim}}<1$ to add stochasticity in the simulation.
Applying bivariate and multivariate GC analysis to the simulated data with the above parameters, we see that when the information propagation time scale is similar to the sampling time scale, both GC methods retrieve the correct information flow (Figure 5B). In general, multivariate GC performs better in identifying information flow compared to bivariate GC. Even when the information propagation time scale is much smaller than that of the sampling rate (${\tau}_{\text{info}}\le {\tau}_{\text{sampling}}$), the information flow can still be correctly retrieved, and the multivariate GC gives the correct input ${W}_{\text{ipsi}}$ value 1, and ${W}_{\text{RC}}$ significantly larger than the null value of 0.5.
The success of retrieval of information flow also depends on the decay constant, although the effect is not strong. As shown in Figure 5C, the performance is almost perfect for short calcium decay times, with ${W}_{\text{IC}}=1$ and ${W}_{\text{RC}}=1$, and longer decay results in slightly worse performance. For the empirical calcium decay rate in the GCaMP3expressing motoneuron data (indicated by the red line), ${W}_{\text{RC}}^{\text{MVGC}}=0.81$. Naturally, using fast calcium indicators in the experiments facilitates the analysis, but the improvement may not be significant.
Adaptive threshold for Fstatistics
We now compute the Granger causality strength (Equation 4) on the smoothed calcium time series. Naively, the Ftest is used to determine whether the Granger causality is significant (Equation 3). This test is exact if the residual, or the prediction error, is Gaussiandistributed and independent. For calcium signals, these assumptions do not hold: the noise is not necessarily Gaussian, and is not independent. While subsequent versions adapting GC to nonlinear and nongaussian dynamics with assumptions of the underlying dynamics have been shown to have wellknown asymptotic null distributions (Kim et al., 2011; Sheikhattar et al., 2018), it did not account for additional irregular stimuli. To develop a general method that can be applied to datasets with diverse stimuli, we use a datadriven method to construct the null distribution.
To test if we can still use the Ftest as a test of significance, we generate data consistent with a null hypothesis by cyclically shuffling the time series of the driving neuron (see Figure 6A). Specifically, for the bivariate GC value from neuron j to another neuron i, $G{C}_{j\to i}$, and the multivariate conditioned on the activity of the rest of neurons $\left\{{f}_{k}\right\},k\ne i,j$, we compute the Fstatistics with the time series ${f}_{j}(t)$ replaced by signal shuffled by a random time constant $\mathrm{\Delta}{t}^{rand}$:
If the significance test remains valid for the calcium time series, we expect almost no significant links and the distribution of the Fstatistics of the shuffled data should be the Fdistribution. However, as Figure 6BC shows using the example dataset f3t2 from De Vico Fallani et al., 2014, for both bivariate and multivariate GC, the distribution of the Fstatistics for the shuffled data, ${\mathcal{F}}^{\text{shuffle}}$, has a shifted support that is larger than that of the Fdistribution. This means that if we use the original Ftest based on the Fdistribution as a null model, even if there are no Gcausal links between two neurons, we will classify the link as significant.
To improve specificity, we instead define the null expectation as the distribution for the Fstatistics of the cyclically shuffled data to test for the significance of GC links. A GC link from the data is significant, if it has an Fstatistics that is above the $1p$ percentile of the null Fstatistics distribution, generated from the cyclically shuffled data. The null Fstatistics distribution can be approximated nonparameterically by Monte Carlo reshuffles. However, since the significance threshold depends on the tail of the Fstatistics distribution, to get an accurate estimation requires a large number of reshuffles and is computationally undesirable.
Instead, to accelerate the computation, we notice that the bulk of the ${\mathcal{F}}^{\text{shuffle}}$ distribution is well approximated by an Fdistribution, with outliers explained by the pseudoperiodicity of the data. We parameterize the null empirical distribution by $\mathcal{F}(F;\alpha ,\beta )$, and find the parameters by maximizing the likelihood (see Materials and methods for details). Repeating the significance test with the parametrized Fdistribution, we find a resulting causality network with fewer links (Figure 6D).
Pipeline summary and GC analysis of motoneuron data
Putting together all the observations of the previous sections, we developed an improved Granger causality analysis pipeline (Figure 7A) to overcome the limitations of naive GC for the investigation of directional information flow in neuron calcium data. These improvements consist of both preprocessing the calcium traces (see GC analysis of motoneuron data) and determining a custom threshold for the significance testing of the presence of a GC link, to avoid the many false positive links that the naive GC analysis finds (see Adaptive threshold for Fstatistics). Preprocessing the data aims to remove undesirable effects of noisy fluorescence signals and of motion artifacts on the GC results. The first step consists of selecting only neurons displaying calcium traces that match their expectation (Removing strange neurons): calcium traces with a low signaltonoise ratio or inconsistent behavior should be disregarded. This step is particularly important for MVGC analysis, in which the traces of all neurons are used to determine the presence of a link between any pair of neurons. Second, we correct the instantaneous change in fluorescence caused by motion artifacts (Correcting for motion artifacts). In the final preprocessing step, the calcium signals are smoothed to remove remaining noise (Smoothing), especially noise that is correlated across neurons and that leads to spurious links in the BVGC approach, and missing links in the MVGC analysis. In the BVGC analysis, correlation in the noise might be perceived as a link between two neurons. In the MVGC case, the correlation in the noise of the receiving neuron and in the noise of neurons on which the GC is conditioned can hide the actual influence of the signal of the driving neuron on the signal of the receiving neuron.
After these preprocessing steps, one can apply the GC algorithm to these calcium traces. In order to avoid calculating GC for many maximum time lags, the choice of this parameter can be chosen to match biological delay of information flow – here 750ms. For each pair of neurons, we timeshift one of the calcium traces in the GC analysis of a pair of neurons to remove causality between the two signals and use this new null model to determine the significance threshold, as motivated in Adaptive threshold for Fstatistics. Notice that for some datasets, the resulting MVGC network has more connections than the corresponding BVGC network. This is because the adaptive thresholds for the significance test are chosen separately for BVGC and MVGC analysis.
We take advantage of anatomical knowledge of the expected network for the motoneuron data set to compare the naive GC results to our customized GC results. The improved GC pipeline significantly increases the ${W}_{\text{IC}}$ ratio from $0.62\pm 0.15$ to $1.00\pm 0.00$ for both BVGC and MVGC ($p<0.001$). ${W}_{\text{IC}}=1$ indicates that there are no contralateral flow, corresponding to what is known from anatomy. On the other hand, the ${W}_{\text{RC}}$ ratio does not significantly increase ($0.59\pm 0.09$ to $0.59\pm 0.15$ for BVGC and $0.52\pm 0.08$ to $0.53\pm 0.10$ for MVGC, $p>0.05$) (Figure 7B). ${W}_{\text{RC}}$ close to 0.5 shows that GC analysis is unable to capture the rostrocaudal information flow pattern.
The improvement of the ipsilateral ratio can be explained by the disappearance of most contralateral links that were spuriously detected by the naive GC analysis. The rostrocaudal ratio ${W}_{\text{RC}}$ is however not significantly improved. The temporal resolution of the recording (4 Hz leading to one frame every 250ms) may not be sufficient to correctly determine the directionality of the information flow as is much slower than the actual propagation speed (a few milliseconds, see Masino and Fetcho, 2005). In a similar manner, the decay timeconstant of GCaMP3 is large compared to the refractory period between two spikes. Results could be further improved by recording at higher frequency and using novel fluorescent calcium sensors such as GCaMP6f, which has faster kinetics.
GC analysis of hindbrain data
In order to investigate the information flow between neurons in a more complex dataset, we analyzed population recordings from neurons in the hindbrain of larval zebrafish performing optomotor response (OMR; data from Severi et al., 2018, also see Figure 8A and B), in which the true information flow is unknown. On each horizontal plane along the dorsoventral axis across the whole hindbrain, neurons were recorded in vivo using twophoton laser scanning calcium imaging at 5.81 Hz while Tg(elavl3:GCaMP5G) transgenic larval zebrafish (previously known as Tg(HUC:GCaMP5G)) responded to a moving grading on their right side ($n=139$ neurons in the reference plane, Figure 8A). Upon motion of the grading, zebrafish larvae responded by swimming forward, often with a left bias (Figure 8C, right). We classified hindbrain neurons in three categories: motorcorrelated neurons in the medial V2a stripe (red in Figure 8A and B, Kinkhabwala et al., 2011; Kimura et al., 2013; Pujala and Koyama, 2019), other motorcorrelated neurons (blue in Figure 8A and B) and non motorcorrelated neurons (gray in Figure 8A and B). Because V2a neurons localized in the medial hindbrain are wellknown command neurons driving locomotion (Kimura et al., 2013; Pujala and Koyama, 2019), we first focused our analyses on medial neurons whose activity was correlated to swimming (red neurons in the rectangle, ${n}_{\text{med}}=20$, in Figure 8A and B). The corresponding calcium traces (Figure 8B) were exempt of motion artifacts and smoothing was not necessary as the scanning of the infrared laser leads to noise that is not correlated between neurons. When applying the naive bivariate Granger causality using the original threshold, we found almost all links to be significant (Figure 8D), suggesting numerous spurious links due to the high overall correlation in neuron activity. Motorcorrelated neurons exhibited a higher drive (Figure 8E), defined as the sum of all its outgoing GC links:
where $j$ spans all other neurons. Critically, the neuronal drive was not correlated to the neuron’s position along the scanning direction (correlation of $r=0.15$ across all 139 neurons and $r=0.03$ across medial swimcorrelated neurons, $p>0.05$, Figure 8—figure supplement 1A), implying that a correction for the scanning delay of the laser scanning microscope was not necessary. In contrast, the drive for a given neuron was correlated to the signaltonoise ratio (SNR) of its calcium trace, computed following details in Materials and Methods, with a Pearson’s correlation coefficient of $r=0.45$ across all neurons and $r=0.69$ across medial neurons ($p<0.001$, Figure 8E). This effect could be corrected using an adaptive threshold tailored to each neuron pair to normalize the GC values. Interestingly, we noticed that the drive was higher for motorcorrelated neurons, especially the ones located in the medial stripe and assumed to likely correspond to V2a reticulospinal neurons (Kinkhabwala et al., 2011; Kimura et al., 2013; Pujala and Koyama, 2019), further justifying focusing our attention on this subset of neurons.
To remove the spurious connections, we computed, as introduced before, the custom threshold for the significance testing of GC links. The difference from motoneuron dataset is that the activities of neurons in the hindbrain are driven by a periodic stimulus. Therefore, to generate the null statistics, we randomize the neural activities of the driver neurons by shuffling the neural signals, while keeping the stimulus fixed. More specifically, we divide the entire recording into 19 epochs of consecutive onset (10 s) and offset (5 s) of the stimulus (see Figure 8B), and generate shuffled calcium traces by randomizing the order of these 19 epochs. For each neuron pair $(i,j)$, we apply $m$ such random shuffles to the calcium trace of the driver neuron $i$, and compute the Fstatistics $\{{F}_{i\to j}^{\text{shuffled}}\}=\{{F}_{i\to j}^{{s}_{1}},\mathrm{\dots},{F}_{i\to j}^{{s}_{m}}\}$. For BVGC, $m=1000$, and for MVGC, $m=100$, we denote the empirical distribution of $\{{F}_{i\to j}^{\text{shuffled}}\}$ by ${\mathcal{F}}^{\text{shuffled}}$.
As in the previous examples, we need to parameterize ${\mathcal{F}}_{i\to j}^{\text{shuffled}}$. As before, for each neuron pair $(i,j)$, we noticed that the distribution ${\mathcal{F}}_{i\to j}^{\text{shuffled}}$ was welldescribed by a constantrescaled Fdistribution (Figure 8—figure supplement 2). Applying an adaptive threshold on ${F}_{i\to j}$ reduces to applying the original threshold on the normalized ${\stackrel{~}{F}}_{i\to j}$, defined by dividing the naive $F$statistics with the expectation value of the $F$statistics generated by the shuffled data. Mathematically,
Finally, by analogy, we use Equation 5 to compute the normalized version of the GC values, directly from the normalized Fstatistics. The final BVGC matrix is shown in Figure 8F. Normalizing the BVGC values allowed the reduction of the correlation between the drive and the SNR, from 0.69 to now 0.49 for medial motorcorrelated neurons (Figure 8G). With this new threshold for significance, most spurious links disappeared. The remaining significant BVGC links between medial motorcorrelated neurons revealed strong information flow between these key command neurons. The MVGC resulted in a much smaller number of significant links (Figure 8H), suggesting a strong global correlation between the neurons, with only the strongest links remaining.
In order to quantify the direction of global information flow, we computed the distribution of the angles of the significant GC links. We normalized the raw distribution by the distribution of all possible GC links in the network, to suppress the influence of neuron positions on the shape of the histogram (see Materials and methods). The normalized distribution (Figure 8I) revealed a biased righttoleft information flow. Because the visual grading was presented on the right side of the fish, moving forward from tail to head, one expects the fish to turn more to the left. As expected, we observe in the plane of interest a biased distribution of the tail angles towards the left (Figure 8C, right). Such righttoleft flow of information was systematic across all motorcorrelated neurons in the same plane (Figure 8, blue), as well as across the 11 planes in the same reference fish, and for all planes of the 10 fish (Figure 8—figure supplement 4, blue, purple, gray respectively).
The observed flow of information from right to left in the hindbrain could reflect the recruitment of excitatory command neurons in order to recruit spinal motor neurons on the left, and start swimming with a left bend. Upstream of command neurons lies a critical high motor center that has been identified by electrical stimulations across multiple vertebrate species and is referred to as the mesencephalic locomotor region, or MLR (Shik et al., 1966). We recently identified, functionally and anatomically, the locus of the MLR in larval zebrafish (CarboTano et al., 2022). We therefore investigated the GC links across three ventral and horizontal planes where numerous command neurons are located in the caudal hindbrain (Kimura et al., 2013; Pujala and Koyama, 2019): the reference plane analyzed in Figure 8, and the adjacent planes respectively 10 µm more dorsal and more ventral than the reference plane (Figure 9A–C, left). In the behavioral recordings acquired for all three planes, we observed a biased distribution of the tail angles towards the left (Figure 9A–C, top right). Accordingly, a righttoleft flow of information was found across motorcorrelated neurons (Figure 9A–C, bottom right). Remarkably, we also detected strong driver neurons (Figure 9A–C, left) in the locus of the right MLR centered in rhombomere 1 (CarboTano et al., 2022, Figure 9B, left). As expected from anatomy (CarboTano et al., 2022) and physiology (Brocard et al., 2010, Cabelguen et al., 2003, Ryczko and Dubuc, 2013, Ryczko et al., 2016), neurons in the right MLR locus drove the activity of neurons in the left MLR locus as well as more caudally located neurons in the ipsilateral and contralateral hindbrain.
Discussion
Granger causality is a simple and effective method for information flow retrieval in large brain areas (Shojaie and Fox, 2021; Barnett et al., 2018; Nicolaou et al., 2012). However, to apply GC analysis on calcium signals from population of neurons at singlecell level requires a careful reexamination of the method, because of the nonlinearity, nonGaussianity, and the possible predominance of correlated noise in the data depending on the imaging method used. We have developed an improved analysis pipeline and demonstrated its usefulness in retrieving information flow in zebrafish either within a dozen of motoneurons (De Vico Fallani et al., 2014) and tens of hindbrain neurons (Oldfield et al., 2020). Since the properties of calcium signals in neuronal populations are comparable across animal models, the tailored analysis pipeline we are proposing is not limited to zebrafish, and can be applied to calcium imaging data across animal species.
The first concern of applying GC analysis to calcium signals is how well GC, a measure of linear dependence among signals, can retrieve information flow in these nonlinear and nonGaussian imaging data. Using synthetic data on small networks with different types of dynamics mimicking calcium signals, from linear autoregressive models to spiking dynamics with exponential decay, we have demonstrated that GC can successfully reconstruct the underlying true dynamical interactions among neurons for a wide range of dynamical variables. In the case of motoneurons that fire repeatedly in bursts, we modeled network activity with two chains of neurons using empirically identified dynamical parameters. In this simulation, we showed that GC can retrieve the expected ipsilateral flow of information with rostraltocaudal directionality, even when the information propagation occurs at a faster scale than the experimental sampling time and the time decay of the calcium transient. Altogether our results on simulations show that GC is a useful analysis method for population recordings of singlecell level calcium signals.
Throughout the examples used, we compared results from bivariate GC and multivariate GC analysis. While BVGC overestimates the true number of causal links, we find that MVGC is prone to a winnertakeall phenomenon that may represent just one of many plausible systemlevel models that can account for the observed data. This is especially problematic when the signals are strongly correlated, for example, in the presence of redundant signals. In addition, MVGC is more prone to the issue of overfitting. We therefore suggest using BVGC on datasets with large numbers (hundreds) of neurons, and MVGC on datasets with few tens of neurons. In general, starting with a small number of neurons of special interest seems wise.
After improving the confidence on the GC method, we applied it to real calcium signals. Because real data is noisy and noise can be correlated or not depending on the imaging method used, GC analysis requires careful examination of multiple pitfalls originating from these calcium signals that we addressed in our analysis pipeline. Our pipeline illustrates the preprocessing of the signals and the postprocessing of the GC results. This effort is complementary to previous studies where new variations of GC were introduced to take into account the nonlinear and nongaussian dynamics in interacting neurons, for example by modifying the underlying VAR dynamics with GLM or VAR with sparsity constraints, and subsequently using adaptive thresholds for significance tests (Kim et al., 2011; Sheikhattar et al., 2018; Francis et al., 2018; Francis et al., 2022). Future work of interest should combine these variations of GC with the pipeline we built here.
Preprocessing of calcium imaging data. As reported previously in Vinck et al., 2015; Nalatore et al., 2007, we find that GC analysis is sensitive to correlated noise, which can often occur in calcium signals, due to, for example, motion artifacts or fluctuations of the imaging laser for nonscanning imaging. We highlight that a motion artifact occurring at a single time point can create spurious GC links and should therefore be removed before proceeding to GC analysis. Similarly, correlated noise across neurons throughout the measurement will also results in spurious links, and will be especially problematic for multivariate GC analysis. The pitfalls associated with correlated noise can be resolved by identifying statistical features of the noise during the decay of calcium signals, and subsequently removing them using total variation differentiation regularization as we used here, or other nonlinear filtering methods. Future versions of the smoothing should be applied to all neurons simultaneously, assuming that the noise is correlated, to best preserve the timeordering of the signal.
Adaptive threshold for significance test in GC analysis. The naive GC analysis uses a significance test that compares the Fstatistics with the null test, a Fdistribution. While this null model is exact for Gaussian processes, for real data with nonlinear dynamics, nonGaussian noise or data with periodic stimuli, the null distribution often deviates and needs to be learned through randomization of the data. In order to develop a more stringent test of significance for GC links, we generated null models of the Fstatistics using random cyclic shuffles of the real data. By using a pairspecific threshold for the Fstatistics, we were able to correct for the effect of the signaltonoise ratio on the GC value.
When real calcium traces were acquired with a spinning disk microscope, that is using a camera chip for detection of all pixels at once illuminated by a laser, the noise was correlated because all pixels over the field of view were recorded simultaneously. In this case, smoothing the fluorescence traces was consequently necessary. In contrast, this step is not advised with twophoton laser scanning microscopy or high speed volumetric lightsheet imaging in which respectively different pixels and planes are simultaneously scanned across the sample. Furthermore, we found that the scanning delay in laser scanning microscopy acquired at 5.81 Hz did not influence the GC results: neurons recorded at the beginning of the scanning cycle did not appear as driver neurons of neurons recorded at the end of the cycle. We found however that the drive of a given neuron computed from the naive GC was correlated with the signaltonoise ratio (SNR) of the calcium trace. Using an adaptive threshold for the significance test, the correlation between the drive and SNR decreased.
Applied to motoneuron data, the new pipeline led to an improvement in removing the spurious contralateral links, but did not significantly improve the directionality of the information flow from rostral to caudal — likely due to the small temporal resolution of the data (4 Hz) while activity propagates in few tens of milliseconds between segments.
In population recordings of brainstem neurons with unknown underlying connectivity, we first focused on neurons whose activity was correlated with the motor output and were located in the area of important command neurons, referred to as the V2a medial stripe (Kinkhabwala et al., 2011) that was previously shown to be critical for generating motor output (Kimura et al., 2013). We found that V2a neurons whose activity was correlated to the motor output were more likely to drive the activity of other neurons than neurons whose activity did not encode the motor output. In the condition where the stimulus is presented on one side of the fish, resulting in the fish often performing contralateral turns, the adapted multivariate Granger causality analysis revealed that the information flows mostly from the ipsilateral to the contralateral side in the brainstem. In order to further understand the origin of the drive underlying this effect, we extended the analysis to all neurons that had correlated activity with the motor output, including neurons outside the V2a medial stripe. Remarkably, we found a cluster of neurons that contained strong driver neurons on the side where the visual stimulus was presented and in the locus of the MLR, a motor center that recruits neurons in the V2a medial stripe to generate motor output (CarboTano et al., 2022). Our observations suggest that MLR cells on the ipsilateral side mainly drive the activity of contralateral MLR cells and motorcorrelated neurons to elicit swimming. Throughout vertebrate species, the MLR has been defined functionally based on electrical stimulations effective at eliciting movement, but had not been yet identified from neuronal recordings during locomotion due to the difficulty to access this deep region in the brainstem during active locomotion. To our knowledge, our result is therefore the first evidence in vertebrates of recruitment of cells in the MLR locus during active locomotion.
We have focused in this manuscript on issues that are specific to calcium imaging. We only briefly discuss problems that are also common to other neural data types such as fMRI. Firstly, we emphasize that with current experimental techniques, only a small portion of the entire nervous system is recorded. While newer versions of GC have been adapted to take into consideration latent variables which contribute equally to all neurons (Guo et al., 2008), or through subsequent removal of causal links (Verny et al., 2017), GC does not in general allow identification of the true underlying connection in the presence of unobserved neurons. However, improved GC is able to identify information flow and an effective functional causal network for the observed neurons. Another issue specific to multivariate GC is the problem of overfitting (Seth et al., 2015), especially when the number of neurons is large compare to the number of independent samples in the data. In this situation, one would need to consider group the activities of neurons (Meshulam et al., 2019), either according to their correlation, or by adjacency, and apply GC to this coarsegrained network. One can also apply GC to subnetworks of observed neurons, and construct a distribution of information flows. Note that other approaches to identify causal links than GC have been developed and applied in genomic data (Verny et al., 2017), gene regulatory networks (Meinshausen et al., 2016) and largescale neural networks (Ke et al., 2019). Further work is now needed to explore whether and how these approaches could be adapted with our improved pipeline to calcium recordings.
Overall, comparing to common practices of using correlation to extract a functional network from calcium neural signals, our GC pipeline provides more information by informing on the direction of information flow. Additionally, compared to existing methods to overcome the pitfalls of calcium data, such as addressing the global additive noise using an inversetime GC (Vinck et al., 2015), our pipeline preprocesses calcium signals taking advantage of their slow decay, and keeps the GC approach as a simple method. By writing down GC value as a function of the Fstatistics, normalized by null models directly constructed using shuffled data, we were able to correct the GC network by adjusting the values of individual links, and thereby to reduce the influence of signaltonoise ratio in a simple yet effective way. The improved pipeline makes it possible to apply Granger causality analysis to data from experiments that simultaneously record from high numbers of neurons in larger networks.
Materials and methods
Synthetic dynamics on small networks
Request a detailed protocolWe simulate artificial neural networks with three different dynamics, evolving at a time resolution of ${\mathrm{\Delta}}_{\text{sim}}$. For simplicity, we set the information flow to occur at the same time delay. The memory kernel, or the influence of neuron $j$ on neuron $i$ separated by time lag of $q$ (in the unit of ${\mathrm{\Delta}}_{\text{sim}}$), is denoted as ${\mathrm{\Gamma}}_{ij,q}$. The three dynamics are:
Vector autoregression models (VAR)
${f}_{i,t}=\sum _{q=1}^{{L}_{\text{true}}}\sum _{j}{\mathrm{\Gamma}}_{ij,q}{f}_{j,tq}+{\xi}_{i,t},$where ${\xi}_{i,t}$ is a white noise. The vector autoregression model is identical to the underlying model for Granger causality analysis, and allows GC to best reproduce the network structure.
Spike dynamics for generalized linear model (GLM)
We simulate spiking neurons with nonlinear interaction is through the generalized linear model (GLM), where the spiking rates of the neurons are given by
${\lambda}_{i,t}=\mathrm{exp}\left({\mu}_{i}+\sum _{q=1}^{{L}_{\text{true}}}\sum _{j}{\mathrm{\Gamma}}_{ij,q}{\sigma}_{j,tq}\right)$and the spike counts
${\sigma}_{i,t}\sim \text{Poiss}({\lambda}_{i,t}).$Here, ${\mu}_{i}$ gives a base rate for spiking, ${\mathrm{\Gamma}}_{ij,q}$ is how the spike rates depend on the neurons’ past history. The spike counts ${\sigma}_{i,t}$ are what we consider as synthetic data, and what we will use for Granger causality analysis.
Spike dynamics regressed by exponential decay (GLMCalcium)
To simulate Calcium neuronal data, which has a fast onset and a slow decay, we add an exponential regression to the spike simulation with the GLM. The resulting trace is
${f}_{i,t}=\sum _{q=0}^{\mathrm{\infty}}{\sigma}_{i,tq}{e}^{q/{\tau}_{\text{ca}}}{\mathrm{\Delta}}_{\text{sim}}.$
We simulate the above three dynamics with the interaction network, ${\mathrm{\Gamma}}_{ij,q}$. The network has equal numbers of excitatory and inhibitory neurons to ensure that the dynamics is stationary. For convenience, the strength of the connection is set to be identical, $\mathrm{\Gamma}}_{ij,q}{}_{q<{L}_{\text{true}}}=c{A}_{ij$, where $c$ is the connection strength, and ${A}_{ij}$ is the signed adjacency matrix.
Crosscorrelation
Request a detailed protocolThe simplest approach of finding links in networks is through crosscorrelation,
The (equaltime) correlation matrix given by ${C}_{ij}(\mathrm{\Delta}t=0)$ is commonly used in calcium signal analysis as a functional network. Nonetheless, it is symmetric, and does not reveal directional information.
Smooth calcium trace using total variation regularization
Request a detailed protocolWe follow the same procedure as in Chen et al., 2019 to smooth calcium traces using total variation regularization. Existing MATLAB packages and custom code are used to reconstruct the smooth trace Chartrand, 2011. Here, we outline this smoothing method briefly.
Consider a raw calcium fluorescence signal $f(t)$ from neurons that are known to generate episodic spikes and silences. These episodic activities can be characterized by the underlying derivative signal $u(t)$ which governs the onset and offset of the episodes, which we assume have an exponentially distributed temporal variations, and is only weakly correlated in time. If we further assume that the noise in $f$ is Gaussian and white, the maximum likelihood reconstruction of the signal $f(t)$ is equivalent to minimizing.
where $A$ is the antiderivative operator, the combination ${\sigma}_{n}^{2}{\tau}_{n}$ is the spectral density of noise floor that we can retrieve from the power spectrum of $f$ at high frequencies, ${\sigma}_{f}$ is the total standard deviation of the signal, and ${\tau}_{f}$ is the typical time scale of these variations. We determine the one unknown parameter ${\tau}_{f}$ by asking that, after smoothing, the cumulative power spectrum of the residue $Auf$ has the least root mean square difference from the cumulative power spectrum of the extrapolated white noise.
Compute signaltonoise ratio
Request a detailed protocolWe compute the signaltonoise ratio (SNR) for the fluorescence trace of a neuron, $f(t)$, using the smoothing procedure above. Namely, for each trace $f(t)$, we first apply the smoothing procedure. Then, the power of the noise, ${P}_{\text{noise}}$, is given by the extrapolated white noise, and equivalently by the sum of the power spectrum of the residue $Auf$. The sum of the power spectrum of $f$ gives the sum of the power of the noise and the power of the signal, ${P}_{\text{signal}}+{P}_{\text{noise}}$. Finally, the SNR is computed as the ratio of ${P}_{\text{signal}}/{P}_{\text{noise}}$.
Synthetic statistical data generation
Request a detailed protocolTo study the effect of slow calcium decays on the accuracy of GC analysis, we generate synthetic data statistically. Specifically, we consider two chains of five neurons, driven externally by block stimuli with progressing timedelay across each chain, and alternates between the left and the right side. The stimulus of the leftside is generated by choosing intervals, ${\tau}_{\text{on}}$, for consecutive onstate (and ${\tau}_{\text{off}}$ for consecutive offstates) uniformly from the set $[\mathrm{min}({\tau}_{\text{on}}^{\text{data}}),\mathrm{max}({\tau}_{\text{on}}^{\text{data}})]$, where the $\{{\tau}_{\text{on}}^{\text{data}}\}$ is the set of duration of all ‘on’states in a real zebrafish experiment, founded by positive finite difference in the fluorescent signal smoothened by totalvariational regularization. The stimulus of the rightside follows with the delay time ${\tau}_{\text{delay}}$ also matches the uniform distribution on the support of the empirical distribution.
We model the sequential activation of the motoneurons by imposing the spike rate as
where $k$ is the ordered index of neurons on each side, ${\tau}_{\text{info}}$ is the time scale of information propagation, and ${\lambda}_{0}$ is a large basal spike rate to simulate the plateau in the neural activity. We use a standard Poisson process with these rates to generate the spike train for each neuron $\{{\sigma}_{i},t\}$. Finally, the observed calcium signal is computed by convoluting the spikes with exponential decays for each neuron $i$,
and downsampling at a sampling frequency ${\tau}_{\text{sampling}}^{1}$.
Neuron identification in the hindbrain dataset
Request a detailed protocolWe analyzed the raw fluorescence movies from Severi et al., 2018 using the suite2p software (Pachitariu et al., 2017), a pipeline for processing twophoton calcium imaging data. The pipeline first performs 2D registration of the plane recording to remove motion artifacts. Regions of interest (ROIs) are detected automatically, classified into cell and noncell, and their calcium traces are extracted by the software. Following the automatic neuron classification, we performed a manual verification to remove spurious cells and add missing cells. These were selected based on the shape of the detected ROI and on the shape of the extracted fluorescence trace. The V2a neurons are identified by their positions and motorcorrelated calcium activitites. The MLR neurons are identified based on their location within the MLR locus as determined in CarboTano et al., 2022, and specifically not using any information about their shape nor the fact that their activity is motorcorrelated.
Normalized polar histogram of Granger causality links
Request a detailed protocolWe computed the distribution of the angles of the significant GC links to quantify the direction of global information flow. To suppress the influence of neuron positions on the direction distribution, we normalized the distribution by dividing it over the distribution of the angles of all possible connections given the position of the neurons.
Data availability
The current manuscript is a computational study, so no data have been generated for this manuscript. All data used in our manuscript have been previously published. Modelling code is available on the github: https://github.com/statbiophys/zebrafishGC, (copy archived at Statistical Biophysics Consortium, 2023). Data files are available at: https://doi.org/10.5281/zenodo.6774389.

ZenodoData access for figures of Chen, Ginoux, Wyart, Mora & Walczak.https://doi.org/10.5281/zenodo.6774389
References

The MVGC multivariate granger causality toolbox: a new approach to grangercausal inferenceJournal of Neuroscience Methods 223:50–68.https://doi.org/10.1016/j.jneumeth.2013.10.018

Detectability of granger causality for subsampled continuoustime neurophysiological processesJournal of Neuroscience Methods 275:93–121.https://doi.org/10.1016/j.jneumeth.2016.10.016

A tutorial review of functional connectivity analysis methods and their interpretational pitfallsFrontiers in Systems Neuroscience 9:175.https://doi.org/10.3389/fnsys.2015.00175

Synaptic drive to motoneurons during fictive swimming in the developing zebrafishJournal of Neurophysiology 86:197–210.https://doi.org/10.1152/jn.2001.86.1.197

Numerical differentiation of noisy, nonsmooth dataISRN Applied Mathematics 2011:1–11.https://doi.org/10.5402/2011/164564

Searching for collective behavior in a small brainPhysical Review. E 99:052418.https://doi.org/10.1103/PhysRevE.99.052418

Measuring and interpreting neuronal correlationsNature Neuroscience 14:811–819.https://doi.org/10.1038/nn.2842

Hierarchy of neural organization in the embryonic spinal cord: grangercausality graph analysis of in vivo calcium imaging dataIEEE Transactions on Neural Systems and Rehabilitation Engineering 23:333–341.https://doi.org/10.1109/TNSRE.2014.2341632

Multivariate granger causality analysis of fMRI dataHuman Brain Mapping 30:1361–1373.https://doi.org/10.1002/hbm.20606

The development and application of optogeneticsAnnual Review of Neuroscience 34:389–412.https://doi.org/10.1146/annurevneuro061010113817

Successful reconstruction of a physiological circuit with known connectivity from spiking activity alonePLOS Computational Biology 9:e1003138.https://doi.org/10.1371/journal.pcbi.1003138

Measurement of linear dependence and feedback between multiple time seriesJournal of the American Statistical Association 77:304–313.https://doi.org/10.1080/01621459.1982.10477803

Measures of conditional linear dependence and feedback between time seriesJournal of the American Statistical Association 79:907–915.https://doi.org/10.1080/01621459.1984.10477110

Circuits controlling vertebrate locomotion: moving in a new directionNature Reviews. Neuroscience 10:507–518.https://doi.org/10.1038/nrn2608

Partial granger causality  eliminating exogenous inputs and latent variablesJournal of Neuroscience Methods 172:79–93.https://doi.org/10.1016/j.jneumeth.2008.04.011

A granger causality measure for point process models of ensemble neural spiking activityPLOS Computational Biology 7:e1001110.https://doi.org/10.1371/journal.pcbi.1001110

Fictive swimming motor patterns in wild type and mutant larval zebrafishJournal of Neurophysiology 93:3177–3188.https://doi.org/10.1152/jn.01248.2004

Coarse graining, fixed points, and scaling in a large population of neuronsPhysical Review Letters 123:178103.https://doi.org/10.1103/PhysRevLett.123.178103

Mitigating the effects of measurement noise on granger causalityPhysical Review. E, Statistical, Nonlinear, and Soft Matter Physics 75:031123.https://doi.org/10.1103/PhysRevE.75.031123

The multifunctional mesencephalic locomotor regionCurrent Pharmaceutical Design 19:4448–4470.https://doi.org/10.2174/1381612811319240011

The mesencephalic locomotor region sends a bilateral glutamatergic drive to hindbrain reticulospinal neurons in a tetrapodThe Journal of Comparative Neurology 524:1361–1383.https://doi.org/10.1002/cne.23911

Granger causality analysis in neuroscience and neuroimagingThe Journal of Neuroscience 35:3293–3297.https://doi.org/10.1523/JNEUROSCI.439914.2015

Upravlenie khod’boĭ i begom posredstvom elektricheskoĭ stimulatsii srednego mozgaBiofizika 11:659–666.

Granger causality: a review and recent advancesAnnual Review of Statistics and Its Application 9:289–319.https://doi.org/10.1146/annurevstatistics040120010930

Imaging neuronal activity with genetically encoded calcium indicatorsCold Spring Harbor Protocols 2012:647–656.https://doi.org/10.1101/pdb.top069609

Learning causal networks with latent variables from multivariate information in genomic dataPLOS Computational Biology 13:e1005662.https://doi.org/10.1371/journal.pcbi.1005662

Analysis of sampling artifacts on the granger causality analysis for topology extraction of neuronal dynamicsFrontiers in Computational Neuroscience 8:75.https://doi.org/10.3389/fncom.2014.00075
Article and author information
Author details
Funding
European Research Council (COG 724208)
 Aleksandra M Walczak
European Research Council (COG 101002870)
 Claire Wyart
New York Stem Cell Foundation (NYSCFRNI39)
 Claire Wyart
Human Frontier Science Program (RGP0063/2018)
 Claire Wyart
Fondation pour la Recherche Médicale (FRM EQU202003010612)
 Claire Wyart
Fondation BettencourtSchueller (FBSdon0031)
 Claire Wyart
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Urs Böhm and Kristen Severi for providing feedback on the hindbrain data, Jenna Sternberg for the data on motoneurons, and Fabrizio De Vico Fallani, Mario Chavez (Paris Brain Institute) and Moritz GrosseWentrup (University of Wien) for insightful feedback. This work was supported by the European Research Council Consolidator Grants n. 724208 (AMW) and n. 101002870 (CW), the New York Stem Cell Foundation (NYSCF) Robertson Award 2016 #NYSCFRNI39 (CW), the Human Frontier Science Program (HFSP) Research Grant #RGP0063/2018 (CW), the Fondation pour la Recherche Médicale Team grant (FRM EQU202003010612) and the Fondation BettencourtSchueller #FBSdon0031 (CW).
Copyright
© 2023, Chen, Ginoux et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 2,140
 views

 300
 downloads

 7
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Neuroscience
Brain water homeostasis not only provides a physical protection, but also determines the diffusion of chemical molecules key for information processing and metabolic stability. As a major type of glia in brain parenchyma, astrocytes are the dominant cell type expressing aquaporin water channel. How astrocyte aquaporin contributes to brain water homeostasis in basal physiology remains to be understood. We report that astrocyte aquaporin 4 (AQP4) mediates a tonic water efflux in basal conditions. Acute inhibition of astrocyte AQP4 leads to intracellular water accumulation as optically resolved by fluorescencetranslated imaging in acute brain slices, and in vivo by fiber photometry in mobile mice. We then show that aquaporinmediated constant water efflux maintains astrocyte volume and osmotic equilibrium, astrocyte and neuron Ca^{2+} signaling, and extracellular space remodeling during optogenetically induced cortical spreading depression. Using diffusionweighted magnetic resonance imaging (DWMRI), we observed that in vivo inhibition of AQP4 water efflux heterogeneously disturbs brain water homeostasis in a regiondependent manner. Our data suggest that astrocyte aquaporin, though bidirectional in nature, mediates a tonic water outflow to sustain cellular and environmental equilibrium in brain parenchyma.

 Neuroscience
Neuromodulatory inputs to the hippocampus play pivotal roles in modulating synaptic plasticity, shaping neuronal activity, and influencing learning and memory. Recently, it has been shown that the main sources of catecholamines to the hippocampus, ventral tegmental area (VTA) and locus coeruleus (LC), may have overlapping release of neurotransmitters and effects on the hippocampus. Therefore, to dissect the impacts of both VTA and LC circuits on hippocampal function, a thorough examination of how these pathways might differentially operate during behavior and learning is necessary. We therefore utilized twophoton microscopy to functionally image the activity of VTA and LC axons within the CA1 region of the dorsal hippocampus in headfixed male mice navigating linear paths within virtual reality (VR) environments. We found that within familiar environments some VTA axons and the vast majority of LC axons showed a correlation with the animals’ running speed. However, as mice approached previously learned rewarded locations, a large majority of VTA axons exhibited a gradual rampingup of activity, peaking at the reward location. In contrast, LC axons displayed a premovement signal predictive of the animal’s transition from immobility to movement. Interestingly, a marked divergence emerged following a switch from the familiar to novel VR environments. Many LC axons showed large increases in activity that remained elevated for over a minute, while the previously observed VTA axon rampingtoreward dynamics disappeared during the same period. In conclusion, these findings highlight distinct roles of VTA and LC catecholaminergic inputs in the dorsal CA1 hippocampal region. These inputs encode unique information, with reward information in VTA inputs and novelty and kinematic information in LC inputs, likely contributing to differential modulation of hippocampal activity during behavior and learning.