Granger causality analysis for calcium transients in neuronal networks, challenges and improvements

Version of Record

Accepted for publication after peer review and revision.

Download
Cite
Share
CommentOpen annotations (there are currently 0 annotations on this page).

Version of Record published: March 15, 2023 (This version)
Accepted Manuscript published: February 7, 2023 (Go to version)
Accepted: February 6, 2023
Preprint posted: June 29, 2022 (Go to version)
Received: June 22, 2022

1. Of interest
Dynamic organization of cerebellar climbing fiber response and synchrony in multiple functional components reduces dimensions for reinforcement learning

Huu Hoang, Shinichiro Tsutsumi ... Keisuke Toyama

Research Article Sep 27, 2023
Further reading

Abstract
Editor's evaluation
Introduction
Results
Discussion
Materials and methods
Data availability
References
Decision letter
Author response
Article and author information
Metrics

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 single-cell 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 non-linear 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 in-depth analysis of the advantages and potential pitfalls of the application of Granger Causality (GC) to calcium imaging data, especially regarding various types of pre-processing. 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.sa0

Introduction

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 whole-brain imaging with single-cell 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, non-invasive methods for obtaining functional connectivity have been developed (Bastos and Schoffelen, 2016). The widely-used 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 time-series 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 ensemble-level neuronal signals, it is unclear how well it can be adapted to address single-cell 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 single-cell 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 single-cell 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 single-cell 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 non-correlated 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 non-linear and non-Gaussian 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 motor-correlated neurons within the recently identified mesencephalic locomotor region in Carbo-Tano 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, \dots, T$ , using prior values of the time series of neuron Y, ${y_{t}}$ , for $t = 1, 2, \dots, T - 1$ , assuming that the two time-series 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 p-value 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’,

x_{t} = {\tilde{a}}_{0} + \sum_{q = 1}^{L} {\tilde{a}}_{q} x_{t - q} + {\tilde{ε}}_{t},

where ${\tilde{a}}_{q}$ are the regression coefficients, the maximum lag $L$ is a hyper-parameter subject to tuning, and ${\tilde{ϵ}}_{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

x_{t} = a_{0} + \sum_{q = 1}^{L} a_{q} x_{t - q} + \sum_{q = 1}^{L} b_{q} y_{t - q} + ε_{t},

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 $ϵ_{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, ${ε_{t}}$ and ${{\tilde{ε}}_{t}}$ , are significantly different requires a one-way analysis of variance (ANOVA) test carried out by computing the F-test statistic:

F_{Y \to X} = \frac{\sum_{t = L + 1}^{T} ({\tilde{ε}}^{2} (t) - ε^{2} (t)) / (M_{f} - M_{r})}{\sum_{t = L + 1}^{T} ε^{2} (t) / (T_{regr} - M_{f})},

where $T_{regr} = T - L$ is the number of time points used in the regression, and $M_{f} = 2 L + 1$ and $M_{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 F-statistics will follow an F-distribution with defined degrees of freedom, $ℱ (F; M_{f} - M_{r}, T_{regr} - M_{f})$ . Hence, testing the computed F-statistic $F_{Y \to X}$ against the F-distribution will distinguish whether Y significantly Granger-causes X.

To measure the amount of reduction in the residual, the bivariate Granger causality value (BVGC) is defined as

\begin{array}{ll} G C_{Y \to X} & \equiv max (\ln \frac{var (\tilde{ε})}{var (ε)}, 0) \\ = max (\ln \frac{\sum_{t = L + 1}^{T} {\tilde{ε}}_{t}^{2} / (T_{regr} - M_{r})}{\sum_{t = L + 1}^{T} ε_{t}^{2} / (T_{regr} - M_{f})}, 0), \end{array}

where the sample variance requires a correction using the number of samples in the regression ( $T_{regr} = T - L$ ) minus the number of regression parameters ( $M_{r}$ or $M_{f}$ ). Mathematically, the GC value is related to the $F$ -test statistic through a monotonic transformation:

G C_{Y \to X} = \max (\ln [(\frac{M_{f} - M_{r}}{T_{regr} - M_{f}} F_{Y \to X} + 1) \frac{T_{regr} - M_{f}}{T_{regr} - M_{r}}], 0) .

A significant-only GC value, $G C_{Y \to X}^{sig} = G C_{Y \to X}$ is called if the null test of the F-statistics is rejected, and zero otherwise.

In multi-component 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 $Z = {Z_{1}, \dots, Z_{N}}$ . For example, to analyze a time series of $N + 2$ neurons ${x_{t}, y_{t}, 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, $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

x_{t} = {\tilde{a}}_{0} + \sum_{q = 1}^{L} {\tilde{a}}_{q} x_{t - q} + \sum_{j = 1}^{N} \sum_{q = 1}^{L} {\tilde{c}}_{j, q} z_{j, t - q} + {\tilde{ε}}_{t},

and

x_{t} = a_{0} + \sum_{q = 1}^{L} a_{q} x_{t - q} + \sum_{q = 1}^{L} b_{q} y_{t - q} + \sum_{j = 1}^{N} \sum_{q = 1}^{L} c_{j, q} z_{j, t - q} + ε_{t},

where ${\tilde{c}}_{j, q}$ and $c_{j, q}$ are the regression coefficients from the conditioned signal z_j. The F-test statistics $F_{Y \to X | Z}$ and the GC values $G C_{Y \to X | 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 p-value 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 pre-define 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

f_{i, t} = \sum_{q = 1}^{L_{true}} \sum_{j} Γ_{i j, q} f_{j, t - q} + ξ_{i, t},

Figure 1 with 7 supplements see all

Download asset Open asset

Granger causality (GC) analysis on synthetic data generated from N = 10 interconnected neurons of specific connectivity.

(A) The network diagram and (B) the connectivity matrix of the model. The excitatory links are represented by arrows in (A) and red dots in (B). The inhibitory links are represented by short bars in (A) and blue dots in (B). Time series are generated and analyzed using correlation analysis and GC, using either a vector autoregressive (VAR) process with connection strength $c = 0.1265$ (**C–D**), a generalized linear model (GLM) with connection strength $c = 0.9$ (**E–F**), and a GLM with calcium filtering as in the experiments (**G–H**). (**D,F,H**) Comparison performance of correlation analysis using trajectory simulated with T = 5000 timepointes, bivariate (BVGC) and multivariate (MVGC) GC in predicting the connectivity matrix. True connectivity are indicated with red dots for excitatory links and blue dots for inhibitory links. The error rate of connectivity identification varies as a function of interaction strength $c$ for both the VAR dynamics (I) and the GLM-calcium-regressed dynamics (J). The multivariate GC always does better at predicting true connectivities than the bivariate (pairwise) GC and than correlation analysis. Error bars in (I) and (J) are standard deviation across 10 random realizations of the dynamics, each lasting T = 5000 timepoints.

where $Γ_{i j, q}$ is a memory kernel capturing the effect from neuron $j$ at time $t - q$ on neuron $i$ at time $t$ , and $ξ_{i, t}$ a Gaussian noise with zero mean and variance $σ^{2}$ . For simplicity, we assume the memory kernel acts only within the preceding two time steps, $L_{true} = 2$ , and set $Γ_{i j, q} = A_{i j} s_{j} c$ for $q ⩽ L_{true}$ , and zero, otherwise. $A_{i j}$ is the adjacency matrix: $A_{i j} = 1$ if neuron $j$ is directly pre-synaptic to neuron $i$ , and $A_{i j} = 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, $Δ_{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 \neq i, j}$ as the conditional variables $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 non-existing 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 signal-to-noise ratio is too small to detect them. In comparison to both GC analyses, the equal-time correlation completely misses the structure (see Figure 1D and Materials and Methods). While cross-correlation with the lag fixed to the model maximum lag $L_{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 single-neuron 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, $σ_{i, t} \sim Poiss (λ_{i, t})$ is governed by a Poisson process with spiking rate $λ_{i, t}$ .

λ_{i, t} = \exp (μ_{i} + \sum_{q = 1}^{L_{true}} \sum_{j} Γ_{i j, q} σ_{j, t - q}),

where $μ_{i}$ is the base rate for spiking, and the memory kernel $Γ_{i j, q}$ is identical to the VAR model. We can further convolve the spikes with an exponential decay function to simulate calcium signals (GLM-calcium; see Materials and Methods for details).

Although the GLM and GLM-calcium dynamics are highly non-linear and non-Gaussian (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 GLM-calcium model, the false positive rate (identified non-existant 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 auto-correlation caused by decay of the calcium signal, the correlation and cross-correlation 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_{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_{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.

Figure 2 with 2 supplements see all

Download asset Open asset

Bivariate and multivariate GC analyses of motoneuron calcium imaging recordings in a zebrafish embryo do not match expected ipsilateral and rostral to caudal information flow.

(A) Effect of the maximum lag $L$ on the false positive rate, that is the proportion of spurious positive links over all negative links, on synthetic data (VAR) with a true lag of 3, for BVGC (*purple*) and MVGC *(green*). BVGC is more prone to detect spurious links, even at the correct lag. (B) The false negative rate, that is the proportion of missed positive links over all positive links, is very sensitive to the choice of the lag. If the chosen lag is very different compared to the true lag, many links will not be detected. (**C,D**) Effect of the lag on the average GC values for synthetic (C) and real motoneuron (D) data (N = 11 neurons, De Vico Fallani et al., 2014). Using the average GC value, we choose a lag of 3 for the motoneuron data analysis, corresponding to 750ms. Further supporting evidence for choosing $L = 3$ is given by the correlation analysis of GC values using different lags in Figure 1. Error bars in (A), (B) and (C) are standard deviation across 10 random realizations of the simulated dynamics with T = 4000 on the network (N = 10 neurons) specified by Figure1A. Error bars in (D) correspond to the standard deviation of the average GC value across the ten recordings of the motoneuron dataset. (E) Dorsal view of the embryonic zebrafish and positions of the left (L) and right (R) motoneurons targeted in calcium imaging (*colored circles*). The arrows represent the wave of activity, propagating ipsilaterally from the rostral (Ro) spinal cord to the caudal (C) spinal cord. (F) Mean image of a 30-somite stage embryonic zebrafish *Tg(s1020t:Gal4; UAS:GCaMP3*), dataset *f3t2* in De Vico Fallani et al., 2014, with motoneurons selected for analysis shown in color. (G) Fluorescence signals (proxy for calcium activity) of selected neurons. We observe an oscillating pattern of activation of the left vs right motoneurons. (H) Representation of the motorneuron flow of information as a directed graph as would be ideally defined by bivariate (*left*) and multivariate (*right*) GC analyses. (I) Resulting directed networks for bivariate (*left*) and multivariate (*right*) GC analyses of fluorescence traces of (G) at maximum lag $L = 3$ : we observe different networks than expected, both for BVGC and MVGC. This is consistent with the effect observed in (A), that BVGC detects many spurious links.

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_{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 single-cell 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 left-right 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 {L, 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} = L$ for $i = 1, 3, \dots, 2 N_{l} - 1$ , and the neurons with $i = 2, 4, \dots, 2 N_{r}$ are found in the right chain with $s_{i} = 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 $f A t B$ ; for example, $f 3 t 2$ 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 bi-axial 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

W_{IC} = \frac{{⟨ G_{i j} ⟩}_{(i, j) | s_{i} = s_{j}}}{{⟨ G_{i j} ⟩}_{(i, j) | s_{i} = s_{j}} + {⟨ G_{i j} ⟩}_{(i, j) | s_{i} \neq s_{j}}},

where $G_{i j}$ is the Granger causality values between neurons $i$ and $j$ . If all connections are equal, then $W_{IC} = 0.5$ . We also define the ratio of ipsilateral rostral-to-caudal links compared to all ipsilateral links as

W_{RC} = \frac{⟨ G_{i j} ⟩_{i < j | s_{i} = s_{j}}}{⟨ G_{i j} ⟩_{i < j | s_{i} = s_{j}} + ⟨ G_{i j} ⟩_{i > j | s_{i} = s_{j}}} .

This ratio gives a measure of the directionality of information flow. A bias in rostrocaudal links results in a departure from $W_{RC} \sim 0.5$ . Notice that based on the anatomy, the ground truth of information flow can be smaller than $W_{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 GC-learned connectivity matrix to test the efficacy of the GC method. For networks with unknown structure, one can compare the GC-learned 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_{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.

Figure 3 with 1 supplement see all

Download asset Open asset

One time-point artifact or one noisy neuron fluorescence trace is sufficient to corrupt the GC algorithm.

(A) Mean image of a 30-somite stage larval zebrafish *Tg(s1020t:Gal4; UAS:GCaMP3*), dataset *f3t1* in De Vico Fallani et al., 2014, with neurons selected for analysis shown in color. (B) Fluorescences traces of neurons of (A). *Top:* A motion artifact consisting of a drop in fluorescence at a single time point is present (star). *Bottom:* the motion artifact is corrected as being the mean of the two surrounding time points. (C) Directed graph resulting from the bivariate (*left*) and multivariate (*right*) GC analysis, before (*top*) and after (*bottom*) motion artifact correction. Before the correction, the neuron with the smallest activity appears to drive many other neurons, perhaps due to its lower SNR. This spurious dominant drive disappears once the motion artifact is corrected, demonstrating that GC is not resilient to artifacts as small as one time point out of one thousand. (D) Mean image of a 30-somite stage larval zebrafish *Tg(s1020t:Gal4; UAS:GCaMP3*), dataset *f5t2* in De Vico Fallani et al., 2014, with neurons selected for analysis shown in color. The fluorescence traces (E) of neurons displayed in red exhibit activity patterns different from the oscillatory pattern expected in motorneurons. GC analysis was run after removal of these neurons. (F) Directed graph resulting from the bivariate (*left*) and multivariate (*right*) GC analysis, before (*top*) and after (*bottom*) removing neuron 13 and neuron 21.

Table 1

The directional measures for the bivariate and multivariate GC network for the calcium transients f3t1 before and after motion-artifact correction, a pre-processing step for calcium signals.

The GC network after motion artifact correction has clear ipsilateral structures.

	BVGC		MVGC
	$W_{IC}$	$W_{RC}$	$W_{IC}$	$W_{RC}$
Before MA correction	0.39	0.53	0.43	0.52
After MA correction	0.66	0.56	0.68	0.44

The motion artifact is a global noise on the observed neurons, a source that is well-known 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 motion-artifact 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 $D F / 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 signal-to-noise ratio (SNR) and the correlated noise between calcium traces. Once removed, $W_{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_{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.

Table 2

The directional measures for the bivariate and multivariate GC network for the calcium transients f5t2 before and after removal of the atypical neurons.

	BVGC		MVGC
	$W_{IC}$	$W_{RC}$	$W_{IC}$	$W_{RC}$
With strange neurons	0.56	0.51	0.59	0.38
Without strange neurons	0.54	0.67	0.53	0.53

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 system-wise 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 10-neurons with the VAR and the GLM-calcium dynamics, where the true underlying signal is $f_{i} (t)$ , for each neuron $i$ . We add a system-wide noise, such that the noise-corrupted signal is

g_{i} (t) = f_{i} (t) + ζ (t),

and $ζ (t)$ is Gaussian white noise, with $⟨ ζ (t) ⟩ = 0$ , and $⟨ ζ (0) ζ (t) ⟩ = σ_{ζ}^{2} δ (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 well-suited 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_{IC}$ is also closer to the biological expectation ( $W_{IC} \sim 1$ ) than before smoothing, with MVGC giving better results (Figure 4D).

Figure 4 with 4 supplements see all

Download asset Open asset

Smoothing the calcium imaging signal can improve the accuracy of GC.

(A) We smooth the noisy calcium imaging signal $f = D F / F$ using total-variational regularization (see Materials and Methods), and plot the example traces of the original and the smoothened neuronal signals, and the residual noise, using the dataset *f3t2* (N = 11 neurons) from De Vico Fallani et al., 2014. (B) The Pearson correlation coefficient shows the residual noise is correlated. (C) GC networks constructed using the original noisy calcium signal compared to ones using the smoothed signal for motorneurons. (D) Weight of ipsilateral GC links, $W_{IC}$ , and the weight of rostal-to-cordal links, $W_{RC}$ , for the original calcium signals and smoothened signals. We expect $W_{IC} \sim 1$ and a null model with no bias for where the links are placed will have $W_{IC} = 0.5$ . The x-axis is sorted using the $W_{IC}$ value for bivariate GC analysis with the original calcium data. Data points are connected with lines to guide the eye. The bivariate GC results are plotted with purple circles, and the multivariate GC results with green diamonds. GC analysis results using the original noisy fluroscence signals are shown with empty markers, while the results from first smoothed data are shown with filled markers.

We note that total variation regularization is a global denoising method that acts on the entire time series. Different from a low-pass 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 high-frequency end of the signal’s power spectrum. To address the problem of high-frequency 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 4-figure 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 time-reversed data as the null hypothesis (Vinck et al., 2015).

Slow calcium timescale

One potential pitfall of performing Granger causality analysis on single-cell 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 (Saint-Amant 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_{L} (t)$ and $I_{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_{L} (t)$ is generated by randomly choosing time intervals when the stimulus state is on, $τ_{on}$ , and when it is off, $τ_{off}$ , uniformly from the sets $[\min (τ_{on}^{data}), \max (τ_{on}^{data})]$ and $[\min (τ_{off}^{data}), \max (τ_{off}^{data})]$ , respectively. Here, ${τ_{on}^{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 total-variational regularization (Chartrand, 2011; Chen et al., 2019). These sets of long-lasting time blocks, where the neurons are always active, simulates the plateau activity of the embryonic motoneurons. After $I_{L} (t)$ is sampled, the stimulus for the right-side neurons, $I_{R} (t)$ , is generated such that the time delay between the “on”-state of $I_{L} (t)$ and $I_{R} (t)$ is sampled from $[\min (τ_{delay}^{data}), \max (τ_{delay}^{data})]$ , where the delay time $τ_{delay}$ matches the empirical stimulus distribution. We model the sequential activation of the motoneurons by imposing the spike rate as a function of $I_{L} (t)$ and $I_{R} (t)$ (see Materials and Methods) and we use a standard Poisson process with these rates to generate the spike train for each neuron, ${σ_{i, t}}$ . The spikes are then convolved with an exponential decay with time scale $τ_{ca}$ to simulate the slow decay time of calcium signals.

Figure 5

Download asset Open asset

Synthetic data mimicking bursting motoneurons shows the slow timescale of calcium signal decay decreases GC performance.

(A) Example traces of the synthetic data for two chains of five neurons (N = 10 neurons for the whole system). The common stimuli, $I_{L}$ and $I_{R}$ , are sampled using the empirically observed on- and off-durations for the neuron bursts ( $τ_{on}$ and $τ_{off}$ ). The time delay of the onset of the right stimulus $I_{R}$ from the offset of the left stimulus $I_{L}$ is also sampled using empirical distribution. Each synthetic data is generated for the duration of 1000s. (B) Success of information flow retrieval, measured by the weight of the ipsilateral GC links, $W_{IC}$ , and the weight of rostral-to-caudal links, $W_{RC}$ , as a function of the information propagation time scale $τ_{info}$ , evaluated at the empirical time scales $τ_{ca} = 2.5 s$ , $τ_{sampling} = 0.25 s$ , and biologically reasonable base spike rate $λ_{0} = 32 s^{- 1}$ . Results from bivariate GC are plotted in *purple*, and results from multivariate GC are in *green*. (C) Success of information flow retrieval for $τ_{info} = 0.025 s$ and $τ_{info} = 0.25 s$ for synthetic data with different calcium decay time constants. Greater calcium decay constants lead to worse information retrivals. Errorbars represents standard deviations across 10 realizations of synthetic data. The empirical calcium time scale $τ_{ca} = 2.5 s$ is indicated by the red vertical line. Error bars in panel (B) and (C) are standard deviation across 10 random realizations of the synthetic data.

The timescales in this model include the sampling time $τ_{sampling}$ , the calcium decay time constant $τ_{ca}$ , and the timescale of information propagation $τ_{info}$ . The sampling time is given by the experimental setup for the motoneuron data we analyze (De Vico Fallani et al., 2014), with $τ_{sampling} = 0.25 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 $τ_{ca} = 2.5 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 $τ_{info} < τ_{sampling}$ . The baseline spike rate is set to $λ_{0} = 32 s^{- 1}$ . Finally, the time resolution for the simulation is chosen to be $Δ_{sim} = 0.0125 s$ , such that $λ Δ_{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 ( $τ_{info} \leq τ_{sampling}$ ), the information flow can still be correctly retrieved, and the multivariate GC gives the correct input $W_{ipsi}$ value 1, and $W_{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_{IC} = 1$ and $W_{RC} = 1$ , and longer decay results in slightly worse performance. For the empirical calcium decay rate in the GCaMP3-expressing motoneuron data (indicated by the red line), $W_{RC}^{MVGC} = 0.81$ . Naturally, using fast calcium indicators in the experiments facilitates the analysis, but the improvement may not be significant.

Adaptive threshold for F-statistics

We now compute the Granger causality strength (Equation 4) on the smoothed calcium time series. Naively, the F-test is used to determine whether the Granger causality is significant (Equation 3). This test is exact if the residual, or the prediction error, is Gaussian-distributed 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 non-gaussian dynamics with assumptions of the underlying dynamics have been shown to have well-known 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 data-driven method to construct the null distribution.

To test if we can still use the F-test 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 ${f_{k}}, k \neq i, j$ , we compute the F-statistics with the time series $f_{j} (t)$ replaced by signal shuffled by a random time constant $Δ t^{r a n d}$ :

f_{j}^{CS} (t) = {\begin{cases} f_{j} (t - Δ t^{rand}), & 0 < t < T - Δ t^{rand} \\ f_{j} (t + T - Δ t^{rand}), & T - Δ t^{rand} < t < T \end{cases} .

Figure 6

Download asset Open asset

Significance tests for Granger causality on real calcium data require new thresholds, generated using null models.

(A) Schematics for generating null models by randomly shuffling the data. The *blue* dots and curves correspond to the original data, the effect neuron $f_{X}$ and the cause neuron $f_{Y}$ ; the *green* ones are generated by random cyclic shuffling of the original $f_{Y} (t)$ , with gray rectangles indicating matching time points before and after the cyclic shuffle. (B) The probability density of the F-statistics of the bivariate Granger causality for an example experiment (dataset *f3t2* from De Vico Fallani et al., 2014 with N = 10 neurons). Blue circles are computed using the smoothed Calcium signals, and the green crosses are computed using the null model generated with cyclic shuffles. The black curve is the F-distribution used in the naive discrimination of GC statistics, and the gray curve is the best fitted F-distribution, fitting an effective number of samples. The gray and the black dashed vertical lines indicate the significance threshold for Granger causality, using the original threshold and the fitted threshold, respectively. (C) Same as (B), for multivariate Granger causality analysis. (D) Comparison between the significant GC network for the experiment *f3t2* using the naive vs. fitted F-statistics threshold.

If the significance test remains valid for the calcium time series, we expect almost no significant links and the distribution of the F-statistics of the shuffled data should be the F-distribution. 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 F-statistics for the shuffled data, $ℱ^{shuffle}$ , has a shifted support that is larger than that of the F-distribution. This means that if we use the original F-test based on the F-distribution as a null model, even if there are no G-causal 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 F-statistics 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 F-statistics that is above the $1 - p$ percentile of the null F-statistics distribution, generated from the cyclically shuffled data. The null F-statistics distribution can be approximated non-parameterically by Monte Carlo reshuffles. However, since the significance threshold depends on the tail of the F-statistics 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 $ℱ^{shuffle}$ distribution is well approximated by an F-distribution, with outliers explained by the pseudo-periodicity of the data. We parameterize the null empirical distribution by $ℱ (F; α, β)$ , and find the parameters by maximizing the likelihood (see Materials and methods for details). Repeating the significance test with the parametrized F-distribution, 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 pre-processing 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 F-statistics). Pre-processing 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 signal-to-noise 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 pre-processing 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.

Figure 7

Download asset Open asset

Comparison of GC analysis results with and without processing the calcium transients and adapting the threshold.

(A) Description of the whole pipeline. Before running the GC analysis on the raw calcium transients ( $D F / F$ ), the lag parameter must be chosen. Then, before running the GC analysis again, motion artifacts are corrected, the fluorescence is smoothed and a new threshold is calculated. (B) Comparison of $W_{IC}$ and $W_{RC}$ before and after applying the pipeline to the GC analyses. Points in the upper right triangle gray-shaded area represent the ratios that have increased in the final GC results. Each fish (n=9) is represented by a color, BVGC by circles and MVGC by triangles. Means are shown in black. We find that using our pipeline, $W_{IC}$ strongly increases and we can clear out the spurious contra-lateral links present in the original GC. $W_{RC}$ is not significantly improved: the recording frequency and GCaMP decay are likely too low and slow compared to the speed of rostro-caudal propagation of the information flow. (C) Network results before (*top row*) and after (*bottom row*) applying our pipeline for computing bivariate Granger causality (BVGC), for all fish. For consistency in the network representation and better ability to compare, we removed the uncharacteristically behaving neurons of the GC analysis before applying the pipeline. Note that more links are found on the side that was better in focus in the recording. (D) Same as (C) but for the multivariate Granger causality (MVGC).

After these pre-processing 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 time-shift 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 F-statistics. 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_{IC}$ ratio from $0.62 \pm 0.15$ to $1.00 \pm 0.00$ for both BVGC and MVGC ( $p < 0.001$ ). $W_{IC} = 1$ indicates that there are no contralateral flow, corresponding to what is known from anatomy. On the other hand, the $W_{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_{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_{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 time-constant 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 two-photon 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: motor-correlated 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 motor-correlated neurons (blue in Figure 8A and B) and non motor-correlated neurons (gray in Figure 8A and B). Because V2a neurons localized in the medial hindbrain are well-known 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_{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. Motor-correlated neurons exhibited a higher drive (Figure 8E), defined as the sum of all its outgoing GC links:

{d r i v e}_{i} = \sum_{j \neq i} G C_{i \to j}

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 swim-correlated 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 signal-to-noise 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 motor-correlated 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.

Figure 8 with 5 supplements see all

Download asset Open asset

Our Granger causality pipeline revealed revealed right-to-left information flow among motor-correlated hindbrain neurons in larval zebrafish performing optomotor (OMR) responses.

(A) Two-photon laser scanning calcium imaging *in vivo* was acquired at 5.81Hz in the hindbrain of *Tg(elavl3:GCaMP5G*) transgenic larvae performing OMR (data from Severi et al., 2018). Forty-one (*red and blue*) out of 139 neurons from a reference plane of a larval zebrafish hindbrain were motor-correlated and selected for subsequent analysis. Neurons located in the medial area (*red*, N = 20 neurons) are in the location of the V2a reticulospinal stripe (Kinkhabwala et al., 2011), distinct from other motor-correlated neurons (*blue*, N = 21). The activity of all other neurons (*gray*, N = 98 neurons) were not correlated to the motor output. (B) Sampled calcium traces of the three categories of neurons mentioned above were ordered by decreasing signal-to-noise ratio (SNR) from top to bottom. The moving grading (*top trace*) was presented to the right of the fish from tail to head to induce swimming, depicted as the tail angle (*bottom trace*, with a zoom over 4s showing two subsequent bouts, positive tail angles correspond to left bends). (C) The distribution of the tail angles for turns (tail angle $α$ with $45^{\circ} < | α | < 90^{\circ}$ ) was biased to the left. (D) GC matrix (*left*) and map of information flow (*right*, the arrows indicate the directionality of the GC links) for the naive BVGC on medial motor-correlated neurons. Of the 380 possible pairs, 351 were found to have a significant drive, suggesting that the naive BVGC algorithm is too permissive to false positive links, and justifying an adapted pipeline to remove spurious connections. (E) The neuronal drive, calculated as the sum of the strength of all outgoing GC links for each neuron, was correlated to the SNR of the calcium traces, especially high for medial neurons ( $r = 0.69$ , compared to $r = 0.45$ across all 139 neurons). (F) Customization of the threshold and normalization of the GC values reduced the number of significant GC links and revealed numerous commissural links that were both rostrocaudal and caudorostral between medial motor-correlated neurons. (G) The new neuronal drive was less correlated with the SNR (from 0.69 to 0.49 for medial neurons). (H) Due to the high correlation in the selected neurons and the relatively high number of neurons (20), the new MVGC pipeline highlighted only the strongest GC links. (I) The density distribution of the significant GC link angles from the new MVGC analyses revealed a biased right-to-left information flow. The directionality of the information flow is consistent with the presence of the stimulus on the right side of the fish and the fish swimming forward with a bias towards the left (C, *right*). These results were conserved across planes of the same fish as well as across fish (Figure 4).

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 F-statistics ${F_{i \to j}^{shuffled}} = {F_{i \to j}^{s_{1}}, \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}^{shuffled}}$ by $ℱ^{shuffled}$ .

As in the previous examples, we need to parameterize $ℱ_{i \to j}^{shuffled}$ . As before, for each neuron pair $(i, j)$ , we noticed that the distribution $ℱ_{i \to j}^{shuffled}$ was well-described by a constant-rescaled F-distribution (Figure 8—figure supplement 2). Applying an adaptive threshold on $F_{i \to j}$ reduces to applying the original threshold on the normalized ${\tilde{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,

{\tilde{F}}_{i \to j} = F_{i \to j} / ⟨ F_{i \to j}^{shuffled} ⟩ .

Finally, by analogy, we use Equation 5 to compute the normalized version of the GC values, directly from the normalized F-statistics. 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 motor-correlated neurons (Figure 8G). With this new threshold for significance, most spurious links disappeared. The remaining significant BVGC links between medial motor-correlated 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 right-to-left 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 right-to-left flow of information was systematic across all motor-correlated 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 (Carbo-Tano 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 right-to-left flow of information was found across motor-correlated 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 (Carbo-Tano et al., 2022, Figure 9B, left). As expected from anatomy (Carbo-Tano 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.

Figure 9

Download asset Open asset

New multivariate Granger causality analysis of motor-correlated hindbrain neurons revealed neurons in the locus of Mesencephalic Locomotor Region (MLR) on the side of the visual grading driving activity of neurons in the contralateral MLR and on both sides in the medulla.

The distributions of the tail angles for turns (tail angle $α$ with $45^{\circ} < | α | < 90^{\circ}$ , *gray, top left* of **A-C**) and the distributions of the new MVGC link angles (*blue, bottom left* of **A-C**) indicate a bias when the visual grading stimulus was presented on the right towards the left across 3 hindbrain planes spaced by 10µm (A: plane 5, dorsal to B: plane 7, analyzed in Figure 8; C: plane 9, ventral to plane 7). Across these three planes, a group of strong driver neurons (*middle and right panels* of **A-C**) was found located on the right, in the locus of the MLR. This is to our knowledge the first evidence of recruitment of cells in the MLR locus during locomotion. As expected from anatomy (Carbo-Tano et al., 2022), right MLR neurons drive contralateral neurons in the left MLR and in the medulla where the reticular formation is located, both on the right and left sides.

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 single-cell level requires a careful re-examination of the method, because of the non-linearity, non-Gaussianity, 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 non-linear and non-Gaussian 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 rostral-to-caudal 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 single-cell 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 winner-take-all phenomenon that may represent just one of many plausible system-level 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 pre-processing of the signals and the post-processing of the GC results. This effort is complementary to previous studies where new variations of GC were introduced to take into account the non-linear and non-gaussian 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 non-scanning 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 time-ordering of the signal.
Adaptive threshold for significance test in GC analysis. The naive GC analysis uses a significance test that compares the F-statistics with the null test, a F-distribution. While this null model is exact for Gaussian processes, for real data with nonlinear dynamics, non-Gaussian 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 F-statistics using random cyclic shuffles of the real data. By using a pair-specific threshold for the F-statistics, we were able to correct for the effect of the signal-to-noise 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 two-photon laser scanning microscopy or high speed volumetric light-sheet 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 signal-to-noise 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 (Carbo-Tano et al., 2022). Our observations suggest that MLR cells on the ipsilateral side mainly drive the activity of contralateral MLR cells and motor-correlated 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 coarse-grained 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 large-scale 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 inverse-time 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 F-statistics, 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 signal-to-noise 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 protocol

We simulate artificial neural networks with three different dynamics, evolving at a time resolution of $Δ_{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 $Δ_{sim}$ ), is denoted as $Γ_{i j, q}$ . The three dynamics are:

Vector autoregression models (VAR)
$f_{i, t} = \sum_{q = 1}^{L_{true}} \sum_{j} Γ_{i j, q} f_{j, t - q} + ξ_{i, t},$
where $ξ_{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 non-linear interaction is through the generalized linear model (GLM), where the spiking rates of the neurons are given by
$λ_{i, t} = \exp (μ_{i} + \sum_{q = 1}^{L_{true}} \sum_{j} Γ_{i j, q} σ_{j, t - q})$
and the spike counts
$σ_{i, t} \sim Poiss (λ_{i, t}) .$
Here, $μ_{i}$ gives a base rate for spiking, $Γ_{i j, q}$ is how the spike rates depend on the neurons’ past history. The spike counts $σ_{i, t}$ are what we consider as synthetic data, and what we will use for Granger causality analysis.
Spike dynamics regressed by exponential decay (GLM-Calcium)

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}^{\infty} σ_{i, t - q} e^{- q / τ_{ca}} Δ_{sim} .$

We simulate the above three dynamics with the interaction network, $Γ_{i j, 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, $Γ_{i j, q} |_{q < L_{true}} = c A_{i j}$ , where $c$ is the connection strength, and $A_{i j}$ is the signed adjacency matrix.

Cross-correlation

Request a detailed protocol

The simplest approach of finding links in networks is through cross-correlation,

C_{i j} (Δ t) = \frac{\sum_{t} (f_{i} (t) - \bar{f_{i}}) (f_{j} (t + Δ t) - \bar{f_{j}})}{σ_{f_{i}} σ_{f_{j}}} .

The (equal-time) correlation matrix given by $C_{i j} (Δ 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 protocol

We 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.

F (u) = \frac{τ_{f}}{σ_{f}} \int_{0}^{T} d t | \dot{u} | + \frac{1}{2 σ_{n}^{2} τ_{n}} \int_{0}^{T} d t | A u - f |^{2},

where $A$ is the antiderivative operator, the combination $σ_{n}^{2} τ_{n}$ is the spectral density of noise floor that we can retrieve from the power spectrum of $f$ at high frequencies, $σ_{f}$ is the total standard deviation of the signal, and $τ_{f}$ is the typical time scale of these variations. We determine the one unknown parameter $τ_{f}$ by asking that, after smoothing, the cumulative power spectrum of the residue $A u - f$ has the least root mean square difference from the cumulative power spectrum of the extrapolated white noise.

Compute signal-to-noise ratio

Request a detailed protocol

We compute the signal-to-noise 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_{noise}$ , is given by the extrapolated white noise, and equivalently by the sum of the power spectrum of the residue $A u - f$ . The sum of the power spectrum of $f$ gives the sum of the power of the noise and the power of the signal, $P_{signal} + P_{noise}$ . Finally, the SNR is computed as the ratio of $P_{signal} / P_{noise}$ .

Synthetic statistical data generation

Request a detailed protocol

To 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 time-delay across each chain, and alternates between the left and the right side. The stimulus of the left-side is generated by choosing intervals, $τ_{on}$ , for consecutive on-state (and $τ_{off}$ for consecutive off-states) uniformly from the set $[\min (τ_{on}^{data}), \max (τ_{on}^{data})]$ , where the ${τ_{on}^{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 total-variational regularization. The stimulus of the right-side follows with the delay time $τ_{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

λ_{L_{k}} (t) = I_{L} (t - k τ_{info}) λ_{0}, λ_{R_{k}} (t) = I_{R} (t - k τ_{info}) λ_{0},

where $k$ is the ordered index of neurons on each side, $τ_{info}$ is the time scale of information propagation, and $λ_{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 ${σ_{i}, t}$ . Finally, the observed calcium signal is computed by convoluting the spikes with exponential decays for each neuron $i$ ,

f_{i, t} = \sum_{q = 0}^{\infty} σ_{i, t - q} e^{- q / τ_{ca}} Δ_{sim}

and downsampling at a sampling frequency $τ_{sampling}^{- 1}$ .

Neuron identification in the hindbrain dataset

Request a detailed protocol

We analyzed the raw fluorescence movies from Severi et al., 2018 using the suite2p software (Pachitariu et al., 2017), a pipeline for processing two-photon 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 non-cell, 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 motor-correlated calcium activitites. The MLR neurons are identified based on their location within the MLR locus as determined in Carbo-Tano et al., 2022, and specifically not using any information about their shape nor the fact that their activity is motor-correlated.

Normalized polar histogram of Granger causality links

Request a detailed protocol

We 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.

The following previously published data sets were used

(2022) Zenodo
Data access for figures of Chen, Ginoux, Wyart, Mora & Walczak.
https://doi.org/10.5281/zenodo.6774389

References

1. Ahrens MB
2. Orger MB
3. Robson DN
4. Li JM
5. Keller PJ
(2013) Whole-Brain functional imaging at cellular resolution using light-sheet microscopy
Nature Methods 10:413–420.
https://doi.org/10.1038/nmeth.2434
- PubMed
- Google Scholar
1. Barnett L
2. Seth AK
(2014) The MVGC multivariate granger causality toolbox: a new approach to granger-causal inference
Journal of Neuroscience Methods 223:50–68.
https://doi.org/10.1016/j.jneumeth.2013.10.018
- PubMed
- Google Scholar
1. Barnett L
2. Seth AK
(2017) Detectability of granger causality for subsampled continuous-time neurophysiological processes
Journal of Neuroscience Methods 275:93–121.
https://doi.org/10.1016/j.jneumeth.2016.10.016
- PubMed
- Google Scholar
(2018) Solved problems for granger causality in neuroscience: a response to Stokes and purdon
NeuroImage 178:744–748.
https://doi.org/10.1016/j.neuroimage.2018.05.067
- PubMed
- Google Scholar
1. Bastos AM
2. Schoffelen JM
(2016) A tutorial review of functional connectivity analysis methods and their interpretational pitfalls
Frontiers in Systems Neuroscience 9:175.
https://doi.org/10.3389/fnsys.2015.00175
- PubMed
- Google Scholar
1. Brocard F
2. Ryczko D
3. Fénelon K
4. Hatem R
5. Gonzales D
6. Auclair F
7. Dubuc R
(2010) The transformation of a unilateral locomotor command into a symmetrical bilateral activation in the brainstem
The Journal of Neuroscience 30:523–533.
https://doi.org/10.1523/JNEUROSCI.3433-09.2010
- PubMed
- Google Scholar
1. Brovelli A
2. Ding M
3. Ledberg A
4. Chen Y
5. Nakamura R
6. Bressler SL
(2004) Beta oscillations in a large-scale sensorimotor cortical network: directional influences revealed by granger causality
PNAS 101:9849–9854.
https://doi.org/10.1073/pnas.0308538101
- PubMed
- Google Scholar
1. Buss RR
2. Drapeau P
(2001) Synaptic drive to motoneurons during fictive swimming in the developing zebrafish
Journal of Neurophysiology 86:197–210.
https://doi.org/10.1152/jn.2001.86.1.197
- PubMed
- Google Scholar
(2003) Bimodal locomotion elicited by electrical stimulation of the midbrain in the salamander Notophthalmus viridescens
The Journal of Neuroscience 23:2434–2439.
https://doi.org/10.1523/JNEUROSCI.23-06-02434.2003
- PubMed
- Google Scholar
Preprint
1. Carbo-Tano M
2. Lapoix M
3. Jia X
4. Auclair F
5. Dubuc R
6. Wyart C
(2022) Functional Coupling of the Mesencephalic Locomotor Region and V2a Reticulospinal Neurons Driving Forward Locomotion
bioRxiv.
https://doi.org/10.1101/2022.04.01.486703
- Google Scholar
1. Chartrand R
(2011) Numerical differentiation of noisy, nonsmooth data
ISRN Applied Mathematics 2011:1–11.
https://doi.org/10.5402/2011/164564
- Google Scholar
1. Chen X
2. Randi F
3. Leifer AM
4. Bialek W
(2019) Searching for collective behavior in a small brain
Physical Review. E 99:052418.
https://doi.org/10.1103/PhysRevE.99.052418
- PubMed
- Google Scholar
1. Cohen MR
2. Kohn A
(2011) Measuring and interpreting neuronal correlations
Nature Neuroscience 14:811–819.
https://doi.org/10.1038/nn.2842
- PubMed
- Google Scholar
1. Cong L
2. Wang Z
3. Chai Y
4. Hang W
5. Shang C
6. Yang W
7. Bai L
8. Du J
9. Wang K
10. Wen Q
(2017) Rapid whole brain imaging of neural activity in freely behaving larval zebrafish (Danio rerio)
eLife 6:e28158.
https://doi.org/10.7554/eLife.28158
- PubMed
- Google Scholar
(2014) Hierarchy of neural organization in the embryonic spinal cord: granger-causality graph analysis of in vivo calcium imaging data
IEEE Transactions on Neural Systems and Rehabilitation Engineering 23:333–341.
https://doi.org/10.1109/TNSRE.2014.2341632
- PubMed
- Google Scholar
1. Deisseroth K
(2011) Optogenetics
Nature Methods 8:26–29.
https://doi.org/10.1038/nmeth.f.324
- PubMed
- Google Scholar
1. Deshpande G
2. LaConte S
3. James GA
4. Peltier S
5. Hu X
(2009) Multivariate granger causality analysis of fMRI data
Human Brain Mapping 30:1361–1373.
https://doi.org/10.1002/hbm.20606
- PubMed
- Google Scholar
(2011) The development and application of optogenetics
Annual Review of Neuroscience 34:389–412.
https://doi.org/10.1146/annurev-neuro-061010-113817
- PubMed
- Google Scholar
1. Florin E
2. Gross J
3. Pfeifer J
4. Fink GR
5. Timmermann L
(2010) The effect of filtering on granger causality based multivariate causality measures
NeuroImage 50:577–588.
https://doi.org/10.1016/j.neuroimage.2009.12.050
- PubMed
- Google Scholar
(2018) Small networks encode decision-making in primary auditory cortex
Neuron 97:885–897.
https://doi.org/10.1016/j.neuron.2018.01.019
- PubMed
- Google Scholar
(2022) Sequential transmission of task-relevant information in cortical neuronal networks
Cell Reports 39:110878.
https://doi.org/10.1016/j.celrep.2022.110878
- PubMed
- Google Scholar
1. Gerhard F
2. Kispersky T
3. Gutierrez GJ
4. Marder E
5. Kramer M
6. Eden U
(2013) Successful reconstruction of a physiological circuit with known connectivity from spiking activity alone
PLOS Computational Biology 9:e1003138.
https://doi.org/10.1371/journal.pcbi.1003138
- PubMed
- Google Scholar
1. Geweke J
(1982) Measurement of linear dependence and feedback between multiple time series
Journal of the American Statistical Association 77:304–313.
https://doi.org/10.1080/01621459.1982.10477803
- Google Scholar
1. Geweke JF
(1984) Measures of conditional linear dependence and feedback between time series
Journal of the American Statistical Association 79:907–915.
https://doi.org/10.1080/01621459.1984.10477110
- Google Scholar
1. Goulding M
(2009) Circuits controlling vertebrate locomotion: moving in a new direction
Nature Reviews. Neuroscience 10:507–518.
https://doi.org/10.1038/nrn2608
- PubMed
- Google Scholar
1. Gow DW
2. Segawa JA
3. Ahlfors SP
4. Lin FH
(2008) Lexical influences on speech perception: a granger causality analysis of MEG and EEG source estimates
NeuroImage 43:614–623.
https://doi.org/10.1016/j.neuroimage.2008.07.027
- PubMed
- Google Scholar
1. Granger CWJ
(1969) Investigating causal relations by econometric models and cross-spectral methods
Econometrica 37:424.
https://doi.org/10.2307/1912791
- Google Scholar
1. Guo S
2. Seth AK
3. Kendrick KM
4. Zhou C
5. Feng J
(2008) Partial granger causality -- eliminating exogenous inputs and latent variables
Journal of Neuroscience Methods 172:79–93.
https://doi.org/10.1016/j.jneumeth.2008.04.011
- PubMed
- Google Scholar
1. Haesemeyer M
2. Robson DN
3. Li JM
4. Schier AF
5. Engert F
(2018) A brain-wide circuit model of heat-evoked swimming behavior in larval zebrafish
Neuron 98:817–831.
https://doi.org/10.1016/j.neuron.2018.04.013
- PubMed
- Google Scholar
Preprint
1. Ke NR
2. Bilaniuk O
3. Goyal A
4. Bauer S
5. Larochelle H
6. Schölkopf B
7. Mozer MC
8. Pal C
9. Bengio Y
(2019) Learning Neural Causal Models from Unknown Interventions
arXiv.
https://doi.org/10.48550/arXiv.1910.01075
- Google Scholar
1. Kim S
2. Putrino D
3. Ghosh S
4. Brown EN
(2011) A granger causality measure for point process models of ensemble neural spiking activity
PLOS Computational Biology 7:e1001110.
https://doi.org/10.1371/journal.pcbi.1001110
- PubMed
- Google Scholar
1. Kimura Y
2. Satou C
3. Fujioka S
4. Shoji W
5. Umeda K
6. Ishizuka T
7. Yawo H
8. Higashijima S
(2013) Hindbrain v2a neurons in the excitation of spinal locomotor circuits during zebrafish swimming
Current Biology 23:843–849.
https://doi.org/10.1016/j.cub.2013.03.066
- PubMed
- Google Scholar
1. Kinkhabwala A
2. Riley M
3. Koyama M
4. Monen J
5. Satou C
6. Kimura Y
7. Higashijima SI
8. Fetcho J
(2011) A structural and functional ground plan for neurons in the hindbrain of zebrafish
PNAS 108:1164–1169.
https://doi.org/10.1073/pnas.1012185108
- PubMed
- Google Scholar
1. Masino MA
2. Fetcho JR
(2005) Fictive swimming motor patterns in wild type and mutant larval zebrafish
Journal of Neurophysiology 93:3177–3188.
https://doi.org/10.1152/jn.01248.2004
- PubMed
- Google Scholar
(2016) Methods for causal inference from gene perturbation experiments and validation
PNAS 113:7361–7368.
https://doi.org/10.1073/pnas.1510493113
- PubMed
- Google Scholar
1. Meshulam L
2. Gauthier JL
3. Brody CD
4. Tank DW
5. Bialek W
(2019) Coarse graining, fixed points, and scaling in a large population of neurons
Physical Review Letters 123:178103.
https://doi.org/10.1103/PhysRevLett.123.178103
- PubMed
- Google Scholar
(2007) Mitigating the effects of measurement noise on granger causality
Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics 75:031123.
https://doi.org/10.1103/PhysRevE.75.031123
- PubMed
- Google Scholar
1. Nguyen JP
2. Shipley FB
3. Linder AN
4. Plummer GS
5. Liu M
6. Setru SU
7. Shaevitz JW
8. Leifer AM
(2016) Whole-brain calcium imaging with cellular resolution in freely behaving Caenorhabditis elegans
PNAS 113:E1074–E1081.
https://doi.org/10.1073/pnas.1507110112
- PubMed
- Google Scholar
(2012) Eeg-Based automatic classification of “ awake ” versus “ anesthetized ” state in general anesthesia using granger causality
PLOS ONE 7:e33869.
https://doi.org/10.1371/journal.pone.0033869
- PubMed
- Google Scholar
1. Oldfield CS
2. Grossrubatscher I
3. Chávez M
4. Hoagland A
5. Huth AR
6. Carroll EC
7. Prendergast A
8. Qu T
9. Gallant JL
10. Wyart C
11. Isacoff EY
(2020) Experience, circuit dynamics, and forebrain recruitment in larval zebrafish prey capture
eLife 9:e56619.
https://doi.org/10.7554/eLife.56619
- PubMed
- Google Scholar
Preprint
(2017) Suite2p: Beyond 10,000 Neurons with Standard Two-Photon Microscopy
bioRxiv.
https://doi.org/10.1101/061507
- Google Scholar
1. Pujala A
2. Koyama M
(2019) Chronology-based architecture of descending circuits that underlie the development of locomotor repertoire after birth
eLife 8:e42135.
https://doi.org/10.7554/eLife.42135
- PubMed
- Google Scholar
(2005) Mapping directed influence over the brain using granger causality and fMRI
NeuroImage 25:230–242.
https://doi.org/10.1016/j.neuroimage.2004.11.017
- PubMed
- Google Scholar
1. Ryczko D
2. Dubuc R
(2013) The multifunctional mesencephalic locomotor region
Current Pharmaceutical Design 19:4448–4470.
https://doi.org/10.2174/1381612811319240011
- PubMed
- Google Scholar
(2016) The mesencephalic locomotor region sends a bilateral glutamatergic drive to hindbrain reticulospinal neurons in a tetrapod
The Journal of Comparative Neurology 524:1361–1383.
https://doi.org/10.1002/cne.23911
- PubMed
- Google Scholar
1. Saint-Amant L
2. Drapeau P
(2001) Synchronization of an embryonic network of identified spinal interneurons solely by electrical coupling
Neuron 31:1035–1046.
https://doi.org/10.1016/s0896-6273(01)00416-0
- PubMed
- Google Scholar
(2015) Granger causality analysis in neuroscience and neuroimaging
The Journal of Neuroscience 35:3293–3297.
https://doi.org/10.1523/JNEUROSCI.4399-14.2015
- PubMed
- Google Scholar
(2018) Investigation of hindbrain activity during active locomotion reveals inhibitory neurons involved in sensorimotor processing
Scientific Reports 8:13615.
https://doi.org/10.1038/s41598-018-31968-4
- PubMed
- Google Scholar
1. Sheikhattar A
2. Miran S
3. Liu J
4. Fritz JB
5. Shamma SA
6. Kanold PO
7. Babadi B
(2018) Extracting neuronal functional network dynamics via adaptive granger causality analysis
PNAS 115:E3869–E3878.
https://doi.org/10.1073/pnas.1718154115
- PubMed
- Google Scholar
(1966)
Upravlenie khod’boĭ i begom posredstvom elektricheskoĭ stimulatsii srednego mozga

Biofizika 11:659–666.
- Google Scholar
1. Shojaie A
2. Fox EB
(2021) Granger causality: a review and recent advances
Annual Review of Statistics and Its Application 9:289–319.
https://doi.org/10.1146/annurev-statistics-040120-010930
- Google Scholar
Software
1. Statistical Biophysics Consortium
(2023) ZebrafishGC, version swh:1:rev:5d480ffadde77f224ffdb0c70d51a0d9a5b646bc
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:71327fe1d9d2117718f2e1498a6be399d519a3aa;origin=https://github.com/statbiophys/zebrafishGC;visit=swh:1:snp:4e84e444b7684c721710e1ce28e9379c38e7f2b5;anchor=swh:1:rev:5d480ffadde77f224ffdb0c70d51a0d9a5b646bc
1. Stokes PA
2. Purdon PL
(2017) A study of problems encountered in granger causality analysis from A neuroscience perspective
PNAS 114:E7063–E7072.
https://doi.org/10.1073/pnas.1704663114
- PubMed
- Google Scholar
(2019) Spontaneous behaviors drive multidimensional, brainwide activity
Science 364:255.
https://doi.org/10.1126/science.aav7893
- PubMed
- Google Scholar
(2012) Imaging neuronal activity with genetically encoded calcium indicators
Cold Spring Harbor Protocols 2012:647–656.
https://doi.org/10.1101/pdb.top069609
- PubMed
- Google Scholar
1. Venkatachalam V
2. Ji N
3. Wang X
4. Clark C
5. Mitchell JK
6. Klein M
7. Tabone CJ
8. Florman J
9. Ji H
10. Greenwood J
11. Chisholm AD
12. Srinivasan J
13. Alkema M
14. Zhen M
15. Samuel ADT
(2016) Pan-neuronal imaging in roaming Caenorhabditis elegans
PNAS 113:E1082–E1088.
https://doi.org/10.1073/pnas.1507109113
- PubMed
- Google Scholar
1. Verny L
2. Sella N
3. Affeldt S
4. Singh PP
5. Isambert H
(2017) Learning causal networks with latent variables from multivariate information in genomic data
PLOS Computational Biology 13:e1005662.
https://doi.org/10.1371/journal.pcbi.1005662
- PubMed
- Google Scholar
(2015) How to detect the granger-causal flow direction in the presence of additive noise?
NeuroImage 108:301–318.
https://doi.org/10.1016/j.neuroimage.2014.12.017
- PubMed
- Google Scholar
1. Warp E
2. Agarwal G
3. Wyart C
4. Friedmann D
5. Oldfield CS
6. Conner A
7. Del Bene F
8. Arrenberg AB
9. Baier H
10. Isacoff EY
(2012) Emergence of patterned activity in the developing zebrafish spinal cord
Current Biology 22:93–102.
https://doi.org/10.1016/j.cub.2011.12.002
- PubMed
- Google Scholar
(2013) Is granger causality a viable technique for analyzing fMRI data?
PLOS ONE 8:e67428.
https://doi.org/10.1371/journal.pone.0067428
- PubMed
- Google Scholar
1. Zhou D
2. Zhang Y
3. Xiao Y
4. Cai D
(2014) Analysis of sampling artifacts on the granger causality analysis for topology extraction of neuronal dynamics
Frontiers in Computational Neuroscience 8:75.
https://doi.org/10.3389/fncom.2014.00075
- PubMed
- Google Scholar

Decision letter

Gordon J Berman

Reviewing Editor; Emory University, United States
Timothy E Behrens

Senior Editor; University of Oxford, United Kingdom

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Granger causality analysis for calcium transients in neuronal networks: challenges and improvements" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Ronald Calabrese as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

1) It was not clear to the reviewers what lessons are specific to the system they are studying and which ones are to be taken as more general lessons. Certainly, dealing with slow calcium dynamics, motion artifacts, and smoothing, are general problems in calcium imaging, but the reviewers were puzzled a bit about how to decide which neurons are "strange" without a lot of system-specific knowledge. This seems to be a rather important effect, and having a bit more guidance as to this point in the discussion would be useful.

2) It was not entirely sure what results one should take home from the hindbrain analysis. It is clear that there is a more-or-less global signal modulating all neural activity, but this is a common occurrence in population recordings (often, one subtracts this off via PCA or another means before proceeding). Is the general lack of causal links (via the MVGC at least) a generic phenomenon in recurrent networks, or is there something more system-specific here? Accordingly, it might be interesting to run a recurrent neural network simulation with similar properties to the hindbrain (and perhaps with correlated driving) to see what GC/MVGC would predict. Is there any hope of these methods finding information flow in recurrent networks, or should we restrict the method to networks where we expect the primary mode of information transmission to be feedforward? Either performing simulations or limiting the stated results to feedforward networks would be important.

3) The reviewers pointed out several studies that adapt the classical GC for both electrophysiology data and calcium imaging data (e.g., [1] A. Sheikhattar et al., "Extracting Neuronal Functional Network Dynamics via Adaptive Granger Causality Analysis", PNAS, Vol. 115, No. 17, E3869-E3878, 2018. [2] N. A. Francis et al., "Small Networks Encode Decision-Making in Primary Auditory Cortex", Neuron, Vol. 97, No. 4, 2018. [3] N. A. Francis et al., "Sequential Transmission of Task-Relevant Information in Cortical Neuronal Networks", Cell Reports, Vol. 39, No. 9, 110878, 2022.). In reference [1], a variation of GC based on GLM log-likelihoods is proposed that addresses the issues of non-linearity, non-stationarity, and non-Gaussianity of electrophysiology data. In [2] and [3], a variation of GC using sparse multi-variate models is introduced with application to calcium imaging data. In particular, all three references use the sparse estimation of the MVAR parameters in order to mitigate overfitting and also use corrections for multiple comparisons that also reduce the number of spurious links. We suggest discussing these relevant references in the introduction (paragraphs 2 and 3) and discussion.

4) A major issue of GC applied to calcium imaging data is that the trials are typically limited in duration, which results in overfitting of the MVAR parameters when using least squares (See references [2] and [3] above, for example). The authors mention on page 4 that they use least squares to estimate the parameters. However, for the networks of ~10 neurons considered in this work, stationary trials of a long enough duration are required to estimate the parameters correctly. We suggest that the authors discuss this point and explicitly mention the trial durations and test whether the trial durations suffice for stable estimation of the MVAR parameters (this can be done by repeating some of the results on the synthetic data and using different trial lengths and then assessing the consistency of the detected GC links).

5) The definition of the "knee" of the average GC values as a function of the lag L needs to be a bit more formalized. In Figure 2H using the synthetic data, the "knee" effect is more clear, but in the real data shown in Figure 2I, the knee is not obvious, given that the confidence intervals are quite wide. Is there a way to quantify the "knee" by comparing the average GC values as well as their confidence bounds along the lag axis?

6) Another key element of existing GC methods applied to large-scale networks is dealing with the issue of multiple comparisons: for instance, in Figures 2, 3, 4, 6, 7, and 8, it seems like all arrows corresponding to all possible links are shown, where the colormap indicates the GC value. However, when performing multiple statistical tests, many of these links can be removed by a correction such as the Benjamini-Hochberg procedure. It seems that the authors did not consider any correction of multiple comparisons; we suggest doing so and adding this to your pipeline.

7) In Figure 5C, the values of W_IC for the MV cases seem to be more than 1, whereas by definition they should be less than or equal to 1. Please clarify.

8) Is there evidence that the lateralized and rostrocaudal connectivity of the motoneurons occurs at the time-scale of ~750 ms? Given that this time scale is long enough for multiple synapses, it could be the case that some contralateral and non-rostrocaudal connections could be "real", as they reflect multi-hop synaptic connections. Please clarify.

9) Redundant signals: throughout the brain, it's expected that a population of neurons can encode the same information. It's unclear how GC (both the original and the modified versions) can handle this redundancy. Given how pervasive redundant signals are in the brain, this should be addressed in both simulation and experimental data. For example, in one of the manuscript's simulated networks, replace one neuron with 10 copies of it, each with identical inputs and outputs but with the weights scaled by 1/10. Such a network is functionally equivalent to the original but may pose some challenges for the various versions of GC. This issue may also account for the MVGC results in the hindbrain dataset. It might be more appropriate to apply GC to groups of neurons (as indeed the authors cited), instead of applying it at the single-cell level with redundant signals.

10) Both BVGC and MVGC appear to be extremely sensitive to any outlier signals. The most worrying aspect is that the authors developed their corrections and pipelines with the benefit of knowing the structure of the underlying system, whereas in the case where GC would be most useful, the user would be unable to rely on prior knowledge of the underlying structure. For instance, the motion artifact in Figure 3a-c was a helpful example of a vulnerability of naive GC, but one could easily imagine scenarios involving an unmeasured disturbance (e.g. the table is bumped) causing a similar artifact, but if the experimenter is unaware of such unmeasured disturbances then they will not be included in Z, and hence can result in the detection of widespread spurious links. There is a circularity here that's concerning. It seems that one already needs to have the answer (e.g. circuit connectivity) in order to clean up the data sufficiently for BVGC or MVGC to work effectively. Perhaps the authors would be interested in incorporating ideas from the systems identification literature, which can include the estimation of unmeasured disturbances, perhaps in conjunction with L1 regularization on the GC links. This is certainly out of scope for the present work, but it would be worth acknowledging the difficulties of unmeasured disturbances and deferring a general solution to future work. Similar considerations apply to a common unmeasured neuronal input (e.g. from a brain region not included in the field of view of the imaging).

11) Interpretation – would it be correct to state that BVGC identifies plausible causal links, while MVGC identifies a plausible system-level model? We think these interpretations, carefully stated, might provide a helpful way of thinking about the two GC approaches. Taking the results of the paper together, neither BVGC nor MVGC is definitive – BVGC may overestimate the true number of causal links but MVGC is prone to a winner-take-all phenomenon that may represent just one of many plausible system-level models that can account for the observed data. This should be more clearly stated in the manuscript.

Reviewer #3 (Recommendations for the authors):

There are numerous typos in the manuscript and some odd choices in the citations in the introduction. Given that there are still some major conceptual issues to address, I'll defer these to the future.

https://doi.org/10.7554/eLife.81279.sa1

Author response

Essential revisions:

1) It was not clear to the reviewers what lessons are specific to the system they are studying and which ones are to be taken as more general lessons. Certainly, dealing with slow calcium dynamics, motion artifacts, and smoothing, are general problems in calcium imaging, but the reviewers were puzzled a bit about how to decide which neurons are "strange" without a lot of system-specific knowledge. This seems to be a rather important effect, and having a bit more guidance as to this point in the discussion would be useful.

We thank the reviewer for pointing out the lack of quantitative definition of the “strange” neurons. In this revision, we introduce a quantitative criteria to identify atypical calcium traces to exclude from the analysis, by recognizing that calcium transients of neurons decay typically with the same profile, i.e. similar decaying time constant when fitted to an exponential function.

If a calcium trace has a decay constant much longer (> tens of s) than an average trace (~ 0.1-5 s), we infer that it’s not the trace of a neuron but likely a glial cell.

This method can be applied generally to systems with neurons expected to have the same characteristics in terms of excitability and calcium buffer. This is the method we already used to identify the time constant of calcium decay (τ = 2.5s) in Figure 3—figure supplement 1 (Figure S7 from first submission) and Results: GC analysis of motoneuron data: slow calcium timescale. In this figure, we now include the overlay of calcium transients from typical neurons as well as from atypical traces (the “strange” neurons). We have also changed the text in Results: GC analysis of motoneuron data: Removing strange neurons as follows:

“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)”.

We also included in the discussion that this method can be generally applied to any other systems relying on calcium imaging to infer neuronal activity.

2) It was not entirely sure what results one should take home from the hindbrain analysis. It is clear that there is a more-or-less global signal modulating all neural activity, but this is a common occurrence in population recordings (often, one subtracts this off via PCA or another means before proceeding). Is the general lack of causal links (via the MVGC at least) a generic phenomenon in recurrent networks, or is there something more system-specific here? Accordingly, it might be interesting to run a recurrent neural network simulation with similar properties to the hindbrain (and perhaps with correlated driving) to see what GC/MVGC would predict. Is there any hope of these methods finding information flow in recurrent networks, or should we restrict the method to networks where we expect the primary mode of information transmission to be feedforward? Either performing simulations or limiting the stated results to feedforward networks would be important.

The reviewers are right. We have now replaced the network in Figure 1 with a more realistic network with feedback. We have also added a figure in the Supplementary Information (Figure 1—figure supplement 1), that tests the performance of GC on VAR dynamics run for randomly connected neurons (N = 10). Each underlying true interaction structure is sampled based on, P_connect, the probability that one neuron has a synaptic connection with the other, and also the connection strength c. As shown in Figure 1— figure supplement 1, MVGC performs well in identifying the true and only the true links when the interaction strength is large.

3) The reviewers pointed out several studies that adapt the classical GC for both electrophysiology data and calcium imaging data (e.g., [1] A. Sheikhattar et al., "Extracting Neuronal Functional Network Dynamics via Adaptive Granger Causality Analysis", PNAS, Vol. 115, No. 17, E3869-E3878, 2018. [2] N. A. Francis et al., "Small Networks Encode Decision-Making in Primary Auditory Cortex", Neuron, Vol. 97, No. 4, 2018. [3] N. A. Francis et al., "Sequential Transmission of Task-Relevant Information in Cortical Neuronal Networks", Cell Reports, Vol. 39, No. 9, 110878, 2022.). In reference [1], a variation of GC based on GLM log-likelihoods is proposed that addresses the issues of non-linearity, non-stationarity, and non-Gaussianity of electrophysiology data. In [2] and [3], a variation of GC using sparse multi-variate models is introduced with application to calcium imaging data. In particular, all three references use the sparse estimation of the MVAR parameters in order to mitigate overfitting and also use corrections for multiple comparisons that also reduce the number of spurious links. We suggest discussing these relevant references in the introduction (paragraphs 2 and 3) and discussion.

We thank the reviewers for pointing out these references. While our GC pipeline focuses on preprocessing the calcium signal by denoising, and post-processing the resulting GC matrix by using an adaptive threshold, these previous studies mentioned by the reviewers emphasize more on modifying the underlying dynamics for GC in order to address the non-linear and non-gaussian dynamics in interacting neurons. In future work, it is of our interest to combine these two approaches. To bring this understanding to the readers, we have incorporated the following change of text. We have included in the introduction (paragraph 3) the following:

“New variations of GC accounting for sparsity in the connection have also been developed and applied on single-cell level to both spiking and calcium signals in the primary auditory cortex of mice (Sheikhattar et al. 2018; Francis et al., 2018, 2022).”

We also modified the Discussion section to include the above reference and added in the text, after we review the two important steps of our GC pipeline, the paragraph connecting our approach better to the body of literature:

“Our pipeline focuses on the pre-processing of the signals and the post-processing of the GC results. This effort is complementary to previous studies where new variations of GC were introduced to take into account the non-linear and non-gaussian 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, 2022). Future work of interest includes combining these variations of GC with our pipeline.”

4) A major issue of GC applied to calcium imaging data is that the trials are typically limited in duration, which results in overfitting of the MVAR parameters when using least squares (See references [2] and [3] above, for example). The authors mention on page 4 that they use least squares to estimate the parameters. However, for the networks of ~10 neurons considered in this work, stationary trials of a long enough duration are required to estimate the parameters correctly. We suggest that the authors discuss this point and explicitly mention the trial durations and test whether the trial durations suffice for stable estimation of the MVAR parameters (this can be done by repeating some of the results on the synthetic data and using different trial lengths and then assessing the consistency of the detected GC links).

We thank the reviewers for pointing out the importance of the length of the time series in the GC calculation and the information missing to describe our dataset lengths. We address this comment with three additional pieces of supporting evidence.

First, we perform the simulation suggested by the reviewer. With the same underlying structure, we simulate two trials with length T, measure the GC values using each trial, and compute the Pearson’s correlation coefficient between the two sets of GC values as a measure of consistency of detecting GC links. The results are included in the Figure 1—figure supplement 3, showing that the stronger the coupling is, the smaller the sample size is required to reconstruct a consistent map of information flow.

Second, to compare with the ground truth, we use synthetic data to check the error rate of link detection as a function of trial lengths, the result of which is included in Figure 1—figure supplement 2. Specifically, we simulate GLM-Calcium dynamics on randomly connected networks, with a fixed probability of connection and a range of connection strength, for different trial durations, and we compute the accuracy of link detection. For MVGC, the accuracy depends on the connection strength. If the connection strength is large, then the true links are mostly detected. If the connection strength is small, then more true links are missed.

These numerical results prompt us to test how consistent the inferred GC matrix is for the zebrafish data, if we use only half of the data. If half of the existing time points are already sufficient to give consistent GC structures, then using all data (and more data) suffices for stable estimation. This is the case for the motoneuron data, which we show in the scatter plot of GC values inferred using the first halves vs. the second halves of each experimental trial. The GC values from the two halves have a Pearson’s correlation coefficient of 0.75 for MVGC, and 0.87 for BVGC. The results for motor neurons in the embryonic spinal cord are shown in an additional supplementary figure, Figure 2—figure Supplement 2.

For hindbrain data, the consistency is not as good. The Pearson’s correlation coefficient is 0.41 for MVGC and 0.76 for BVGC. We note that the similarity is higher for the motoneuron dataset due to the cyclic activity pattern leading to two similar halves.

We added in Results: GC analysis of motoneuron data: Naive GC the above results, specifically addressing the point of trial duration:

“… 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's 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.”

We also added a supplementary figure to discuss the overfitting issue (Figure 8—figure Supplement 3). Specifically, we compare the resulting GC links for the hindbrain data, using VAR dynamics, and VAR plus regularization. We show that in MVGC, not observing many links is not purely a problem of overfitting. And that if including regularization, one still needs to use the adaptive threshold to increase the connection.

5) The definition of the "knee" of the average GC values as a function of the lag L needs to be a bit more formalized. In Figure 2H using the synthetic data, the "knee" effect is more clear, but in the real data shown in Figure 2I, the knee is not obvious, given that the confidence intervals are quite wide. Is there a way to quantify the "knee" by comparing the average GC values as well as their confidence bounds along the lag axis?

We thank the reviewers for asking for clarification. We can compare the average GC values as suggested, or the individual pairwise GC values when adding an extra lag. Using 3 lags or more yields very similar results, so we took the decision to keep the simplest model with fewest lags. To clarify, we added in Results: GC analysis of motoneuron data: Naive GC:

“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).”

The additional Figure 2—figure supplement 1 shows the “Comparison of Granger causality analysis results at different maximum lags. Correlation between the GC values obtained for each pair of neurons using max lag L and max lag L + 1. All fish (n=9) of the motoneuron data set are used and represented by color (same as Figure 7). The similarity in Granger causality values, measured by the Pearson correlation coefficient r, increases as L increases and then stabilizes after maximum lag L = 3 (r = 0.97 both for BVGC and MVGC between L = 3 and L = 4, as well as between L = 4 and L = 5), justifying our decision to use L = 3. Using a larger maximum lag does not bring more information but brings more complexity, and increases the computation time (especially for MVGC).”

6) Another key element of existing GC methods applied to large-scale networks is dealing with the issue of multiple comparisons: for instance, in Figures 2, 3, 4, 6, 7, and 8, it seems like all arrows corresponding to all possible links are shown, where the colormap indicates the GC value. However, when performing multiple statistical tests, many of these links can be removed by a correction such as the Benjamini-Hochberg procedure. It seems that the authors did not consider any correction of multiple comparisons; we suggest doing so and adding this to your pipeline.

In our original analysis, we had performed multiple comparisons by applying a Bonferroni correction to the p-value threshold. We now clarify this in the main text by adding in Results: Granger Causality, after Equation 7,

“…with the number of parameters adjusted accordingly in the reduced and the full model. The significant test is adjusted using multiple comparisons by applying a Bonferroni correction to the p-value threshold.”

7) In Figure 5C, the values of W_IC for the MV cases seem to be more than 1, whereas by definition they should be less than or equal to 1. Please clarify.

We thank the reviewers for carefully reading our submission. In our first submission, W_IC can be more than 1 when the averaged GC value for contralateral links is negative (Equation 10). This negative GC value can happen when the residuals of fitting from the full and from the reduced model are similar, and after dividing by the degrees of freedom, (T_{regr –} M_r) and (T_{regr –} M_f), respectively, leads to $v a r (\tilde{ε}) < v a r (ε)$ . We did not impose that GC value has to be nonnegative in our definition.

To improve clarity, we have introduced in the revisions the non-negativeness in the definition of GC in Equation 4 and 5 such that GC is defined as the maximum between 0 and the log of ratio of residuals. This way, the W_IC and W_RC will always be less than or equal to 1. We have changed Figure 5C accordingly. The problem of negative GC value does not occur in the results of GC applied to real motoneuron and hindbrain data.

8) Is there evidence that the lateralized and rostrocaudal connectivity of the motoneurons occurs at the time-scale of ~750 ms? Given that this time scale is long enough for multiple synapses, it could be the case that some contralateral and non-rostrocaudal connections could be "real", as they reflect multi-hop synaptic connections. Please clarify.

In the spinal cord, motor neurons form synapses onto muscle fibers and are not synaptically connected among themselves. At the embryonic stage, motor neurons are part of an electronically-coupled network relying on GAP junctions. During twitches, motor neurons are recruited by excitatory neurons located in the hindbrain and enable their sequential activation from head to tail. We therefore use this dataset to infer flow of information, not synaptic connectivity. This important point is now added to the result section.

Physiologically, the two chains of motor neurons are not directly connected to each other. However, all motor neurons are controlled on each side by common driver neurons, in order to be sequentially activated (St Amant and Drapeau, Neuron 2001). Overall, each chain of neurons on one side can be recruited at a frequency up to 2hz, so with a typical 250ms time difference (Warp et al., Current Biology 2012; Wan et al., Cell 2019 https://pubmed.ncbi.nlm.nih.gov/31564455/). Depending on the position of a neuron in the focal plane of the spinning disk microscope, we have good reasons to think that the kinetics of the calcium transient could appear fast (in the plane) versus slow (out of the focal plane). This artifact may lead to variations in the apparent pattern of recruitment of motor neurons in the rostral caudal axis. This could explain why in the data, we sometimes observe differences in the recruitment pattern leading to a caudal motor neuron being activated a few hundreds milliseconds before a rostral one.

This physiological structure is explained in Results: GC analysis for motoneuron data: The data, when the motoneuron dataset is introduced. To further clarify that the ground truth is not strictly W_RC = 1, we added in Results: GC analysis for motoneuron data: Defining directional biases, after the definition of

W_RC:

“Notice that based on the anatomy, the ground truth of information flow can be smaller than W_RC = 1, because motoneuron activity has been shown to rely on command neurons in the brainstem that activates motoneurons sequentially in the spinal cord, not that the motor neurons synapse onto each other.”

9) Redundant signals: throughout the brain, it's expected that a population of neurons can encode the same information. It's unclear how GC (both the original and the modified versions) can handle this redundancy. Given how pervasive redundant signals are in the brain, this should be addressed in both simulation and experimental data. For example, in one of the manuscript's simulated networks, replace one neuron with 10 copies of it, each with identical inputs and outputs but with the weights scaled by 1/10. Such a network is functionally equivalent to the original but may pose some challenges for the various versions of GC. This issue may also account for the MVGC results in the hindbrain dataset. It might be more appropriate to apply GC to groups of neurons (as indeed the authors cited), instead of applying it at the single-cell level with redundant signals.

We agree with the reviewer that the problem of redundant signals is subtle. Redundant signals are not equivalent to repeated signals. A repeated signal is an exact copy of another signal, which experimentally could come from a mis-partition of a single neuron. In this case, GC is unable to detect the correct links, especially the connection between the repeated neuron and the others. In contrast, redundant encoding is more in the structure of the network. Because of stochasticity, the dynamics of two neurons – with identical input and output weights to other neurons – can be different, as we test it on the network in Figure 1 with VAR and with GLM-Calcium dynamics. Such a difference in the dynamics allows MVGC to detect the true structure despite redundant encoding. We have added a paragraph in Results: GC analysis of a small synthetic neuronal net, and an additional supplementary figure (Figure 1 —figure supplement 7) to clarify this issue. The added text is:

“Redundant signals. 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 can still reveal these connections (Figure 1— figure Supplement 7).”

The additional Figures 1 —figure supplement 7 shows:

Presence of redundant signals harms the performance of MVGC. GC analysis performed on synthetic data generated using VAR dynamics on the network structure in the main figure, for T = 5000 data points. GC reveals all links for networks with redundant structures, where an “artificial” neuron 11 is created with identical input and output strength as neuron 1 (redundant structure, middle column). However, MVGC fails to identify causal links when the signal from neuron 1 is copied to create another “artificial” neuron 11 (redundant signal, right most column). BVGC is able to identify the underlying true connectivity.

For the hindbrain data in the last figure, we previously tried grouping the neurons. However, we find that the most correlated neurons are spatially located far from each other, as one would expect given that pairs of pre- and postsynaptic neurons may be connected through a very long axon. This makes grouping difficult: grouping neurons by correlation will disable us finding the spatial information flow, while grouping neurons by spatial proximity washes away meaningful signals.

10) Both BVGC and MVGC appear to be extremely sensitive to any outlier signals. The most worrying aspect is that the authors developed their corrections and pipelines with the benefit of knowing the structure of the underlying system, whereas in the case where GC would be most useful, the user would be unable to rely on prior knowledge of the underlying structure. For instance, the motion artifact in Figure 3a-c was a helpful example of a vulnerability of naive GC, but one could easily imagine scenarios involving an unmeasured disturbance (e.g. the table is bumped) causing a similar artifact, but if the experimenter is unaware of such unmeasured disturbances then they will not be included in Z, and hence can result in the detection of widespread spurious links. There is a circularity here that's concerning. It seems that one already needs to have the answer (e.g. circuit connectivity) in order to clean up the data sufficiently for BVGC or MVGC to work effectively. Perhaps the authors would be interested in incorporating ideas from the systems identification literature, which can include the estimation of unmeasured disturbances, perhaps in conjunction with L1 regularization on the GC links. This is certainly out of scope for the present work, but it would be worth acknowledging the difficulties of unmeasured disturbances and deferring a general solution to future work. Similar considerations apply to a common unmeasured neuronal input (e.g. from a brain region not included in the field of view of the imaging).

We corrected for unmeasured inputs by preprocessing the data using characteristic features for calcium transients, and not assuming any underlying network structures. As we expect, all calcium transients follow a fast rise and a slow decay that can be fitted by exponentials, disturbances such as the table being bumped can be identified and corrected. We also emphasize that knowing the underlying structure of the motoneuron network only prompts us to clean up the correlated noise, and that this method developed is applicable to general datasets without a priori known structures, such as the hindbrain dataset.

We thank the reviewers for pointing out literature in system identification theory. To provide more detailed context, we added in the Discussion section sentences on recent development of GC and other causality measures that considers hidden input, and discussed their advantages and disadvantages:

“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 [49], or through subsequent removal of causal links [50], GC does not in general allow identification of the true underlying connection in the presence of unobserved neurons. However, it is able to identify information flow and an effective functional causal network for the observed neurons.”

11) Interpretation – would it be correct to state that BVGC identifies plausible causal links, while MVGC identifies a plausible system-level model? We think these interpretations, carefully stated, might provide a helpful way of thinking about the two GC approaches. Taking the results of the paper together, neither BVGC nor MVGC is definitive – BVGC may overestimate the true number of causal links but MVGC is prone to a winner-take-all phenomenon that may represent just one of many plausible system-level models that can account for the observed data. This should be more clearly stated in the manuscript.

We agree with the referees, and have added the following sentences to Discussion:

“We also compared results from bivariate GC and multivariate GC analysis. While BVGC overestimate the true number of causal links, MVGC is prone to a winner-take-all phenomenon that may represent just one of many plausible system-level 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 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 interest seems wise.”

https://doi.org/10.7554/eLife.81279.sa2

Article and author information

Author details

Xiaowen Chen

Laboratoire de physique de l'École normale supérieure, CNRS, PSL University, Paris, France

Contribution
Conceptualization, Data curation, Software, Formal analysis, Investigation, Visualization, Methodology, Writing – original draft, Writing – review and editing

Contributed equally with
Faustine Ginoux

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-4029-1805
Faustine Ginoux

Spinal Sensory Signaling team, Sorbonne Université, Paris Brain Institute (Institut du Cerveau, ICM), Paris, France

Contribution
Conceptualization, Data curation, Software, Investigation, Visualization, Methodology, Writing – original draft, Writing – review and editing

Contributed equally with
Xiaowen Chen

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-5031-7346
Martin Carbo-Tano

Spinal Sensory Signaling team, Sorbonne Université, Paris Brain Institute (Institut du Cerveau, ICM), Paris, France

Contribution
Validation, Investigation

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-1936-7174
Thierry Mora

Laboratoire de physique de l'École normale supérieure, CNRS, PSL University, Paris, France

Contribution
Conceptualization, Supervision, Investigation, Methodology, Writing – original draft, Project administration, Writing – review and editing

Contributed equally with
Aleksandra M Walczak and Claire Wyart

For correspondence
thierry.mora@phys.ens.fr

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-5456-9361
Aleksandra M Walczak

Laboratoire de physique de l'École normale supérieure, CNRS, PSL University, Paris, France

Contribution
Conceptualization, Supervision, Funding acquisition, Investigation, Methodology, Writing – original draft, Project administration, Writing – review and editing

Contributed equally with
Thierry Mora and Claire Wyart

For correspondence
aleksandra.walczak@phys.ens.fr

Competing interests
Senior editor, eLife

"This ORCID iD identifies the author of this article:" 0000-0002-2686-5702
Claire Wyart

Spinal Sensory Signaling team, Sorbonne Université, Paris Brain Institute (Institut du Cerveau, ICM), Paris, France

Contribution
Conceptualization, Resources, Data curation, Supervision, Funding acquisition, Investigation, Visualization, Methodology, Writing – original draft, Writing – review and editing

Contributed equally with
Thierry Mora and Aleksandra M Walczak

For correspondence
claire.wyart@icm-institute.org

Competing interests
Reviewing editor, eLife

"This ORCID iD identifies the author of this article:" 0000-0002-1668-4975

Funding

European Research Council (COG 724208)

Aleksandra M Walczak

European Research Council (COG 101002870)

Claire Wyart

New York Stem Cell Foundation (NYSCF-R-NI39)

Claire Wyart

Human Frontier Science Program (RGP0063/2018)

Claire Wyart

Fondation pour la Recherche Médicale (FRM- EQU202003010612)

Claire Wyart

Fondation Bettencourt-Schueller (FBS-don-0031)

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 Grosse-Wentrup (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 #NYSCF-R-NI39 (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 Bettencourt-Schueller #FBS-don-0031 (CW).

Senior Editor

Timothy E Behrens, University of Oxford, United Kingdom

Reviewing Editor

Gordon J Berman, Emory University, United States

Version history

Received: June 22, 2022
Preprint posted: June 29, 2022 (view preprint)
Accepted: February 6, 2023
Accepted Manuscript published: February 7, 2023 (version 1)
Version of Record published: March 15, 2023 (version 2)

Copyright

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

1,235

Page views
178

Downloads
1

Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Mendeley

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Xiaowen Chen
Faustine Ginoux
Martin Carbo-Tano
Thierry Mora
Aleksandra M Walczak
Claire Wyart

(2023)

Granger causality analysis for calcium transients in neuronal networks, challenges and improvements

eLife 12:e81279.

https://doi.org/10.7554/eLife.81279

Categories and tags

Research organism

Zebrafish

Share this article

Cite this article

Granger causality (GC) analysis on synthetic data generated from N = 10 interconnected neurons of specific connectivity.

Bivariate and multivariate GC analyses of motoneuron calcium imaging recordings in a zebrafish embryo do not match expected ipsilateral and rostral to caudal information flow.

One time-point artifact or one noisy neuron fluorescence trace is sufficient to corrupt the GC algorithm.

The directional measures for the bivariate and multivariate GC network for the calcium transients f3t1 before and after motion-artifact correction, a pre-processing step for calcium signals.

The directional measures for the bivariate and multivariate GC network for the calcium transients f5t2 before and after removal of the atypical neurons.

Smoothing the calcium imaging signal can improve the accuracy of GC.

Synthetic data mimicking bursting motoneurons shows the slow timescale of calcium signal decay decreases GC performance.

Significance tests for Granger causality on real calcium data require new thresholds, generated using null models.

Comparison of GC analysis results with and without processing the calcium transients and adapting the threshold.

Our Granger causality pipeline revealed revealed right-to-left information flow among motor-correlated hindbrain neurons in larval zebrafish performing optomotor (OMR) responses.

New multivariate Granger causality analysis of motor-correlated hindbrain neurons revealed neurons in the locus of Mesencephalic Locomotor Region (MLR) on the side of the visual grading driving activity of neurons in the contralateral MLR and on both sides in the medulla.

Author details

Xiaowen Chen

Contribution

Contributed equally with

Competing interests

Faustine Ginoux

Contribution

Contributed equally with

Competing interests

Martin Carbo-Tano

Contribution

Competing interests

Thierry Mora

Contribution

Contributed equally with

For correspondence

Competing interests

Aleksandra M Walczak

Contribution

Contributed equally with

For correspondence

Competing interests

Claire Wyart

Contribution

Contributed equally with

For correspondence

Competing interests

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)

Categories and tags

Research organism

Further reading