Pattern-Based Granger Causality Mapping in fMRI

Abstract

Since its development, the multivoxel pattern analysis (MVPA) method has been widely used to study high-level cognitive function in the brain. The results of the MVPA indicate that the spatial pattern of functional MRI data contains useful information. In addition to the spatial pattern analysis of the brain functions, effective connectivity can also be analyzed between the spatial pattern-based information. In this article, we propose a multivoxel pattern-based causality mapping method to explore influences between the spatial pattern-based information in the brain. The method applies the Granger causality to interested regions of the brain in terms of spatiotemporal pattern-based data, which are known to play an important role in dealing with high-level functions of the brain. The method can compose a causality map throughout the entire brain for any specified region of interest. Both simulations and experiments were performed to show the performance of the proposed method, and the existence and analyzability of the connectivity between pattern-based information in the brain were verified.

Introduction

Functional magnetic resonance imaging (fMRI) has contributed to objective insights into brain functions by detecting the blood oxygenation level-dependent (BOLD) signals related to cerebral hemodynamic responses. Since the invention of the fMRI, the relationship between the brain regions and functions was the main interest of most research, and many different approaches of brain function interpretation have been proposed to understand high-level brain functions. For instance, various functional connectivity analyses have been proposed and widely used to show the interaction among brain regions in many states (Biswal et al., 1995; Boccalettia et al., 2006; Greicius et al., 2003; Hampson et al., 2002; Lowe et al., 2000; Martino et al., 2008). In particular, the effective connectivity analysis that can reveal directional information of influence between the brain regions has been highlighted as a key technology for understanding the dynamic and parallel functions of the brain (Abler et al., 2006; Deshpande et al., 2008; Friston, 2009; Gao et al., 2008; Goebel et al., 2003; Marreiros et al., 2008; Roebroeck et al., 2005). On the other hand, multivoxel pattern analysis (MVPA) has been proposed to explain virtually unlimited functions of the brain by importing spatial pattern-based brain information (Haxby et al., 2001; Kamitani and Tong, 2005; Norman et al., 2006; Peelen and Downing, 2007). The MVPA proved to be robust, which accelerated brain function studies since its development.

Inspired by the effective connectivity analysis and the MVPA, we can infer that the spatial pattern-based information in the brain is meaningful and can also have a causal relationship. However, there has been no approach to detect the causality of pattern-based information, and the necessity of the method is gathering strength (Anderson and Oates, 2010; Kriegeskorte and Kreiman, 2012; Raizada and Kriegeskorte, 2010). The causal relationship between the pattern-based brain information is expected to show dynamic features of the complex functions; hence, it can improve the understanding of brain functions. Therefore, we present a new method to analyze the causality of the pattern-based brain information.

Granger causality mapping (GCM) is one of the most widely used methods to analyze effective connectivity in the brain. The GCM imports the concept of Granger causality (Granger, 1969, 1980) to detect the influence and its direction by exploiting temporal precedence information. In the context of the Granger causality, the fMRI time series of brain regions are described using the vector autoregressive (VAR) model, which predicts current signals by linear combination of past signals. If the past signals of time series improve the prediction accuracy of the current values, it means that there is informational influence.

In conventional applications of the Granger causality into fMRI data analysis, each interested region or voxel is represented by a single fMRI time series. Recently, there have been several researches of Granger causality while taking into account the spatial information of the brain under the assumption of sparseness (Garg et al., 2011; Tang et al., 2012; Valdés-Sosa et al., 2006). In this article, to reveal the causal relationship between pattern-based information, we propose a new approach to apply the Granger causality to the interested brain regions without regressions or the assumption of sparseness. To do this, the multiple time series data are inferred as the spatiotemporal pattern-based data of the region. Through a process that is similar to the GCM, we can verify whether one pattern-based data has useful information in predicting the other by using the temporal precedence information. In other words, we can analyze the existence and direction of Granger causality between spatial pattern-based data that has been proved to be meaningful by studies using the MVPA. Similar to the GCM, the proposed method can compute directional influence between any selected two regions. The method does not require a predefined region of interest (ROI), and it can compose the Granger causality map through the entire brain for any specified ROI.

We performed simulations to confirm that the neuronal pattern Granger causality can be detected from the simulated fMRI data. We also conducted fMRI experiments using a visuomotor task to show that the pattern-based causality between visual and motor cortices can actually be detected by the proposed method.

Theory

GCM is an application of Granger causality to the brain signal, and the proposed method adopts a similar approach. Applying the VAR model, we can calculate the Granger causality that implies the directional influence using the temporal precedence (Geweke, 1982). The VAR can be applied to fit a discrete timeseries vector \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \bf x} [ n ] = ( x_1 [ n ] , \ldots , x_M [ n ] ) ^{\rm T}$$ \end{document} (Kay, 1988), with which the proposed method analyzes the dynamics of the fMRI signal. The VAR model is a prediction model in which a current vector is predicted by its past information, and the VAR model with the order p can be written as follows: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}{ \bf x} [ n ] = \mathop \sum_{i = 1}^p{ \bf A} [ i ] { \bf x} [ n - i ] + { \bf u} [ n ] , \tag{1}\end{align*} \end{document}

where u[n] is multivariate white noise. In Equation (1), x[n] represents the fMRI time series of multiple voxels and is estimated as a linear combination of its last p signals using autoregression coefficient matrices A[i].

For two arbitrary regions of X and Y in the brain, each region can be fit to the VAR independently. \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{split}{ \bf x} [ n ] = \mathop \sum_{i = 1}^p{ \bf A}_x [ i ] { \bf x} [ n - i ] + { \bf u} [ n ] \qquad & { \rm covar} ( { \bf u} [ n ] ) = \Sigma_1 \\ { \bf y} [ n ] = \mathop \sum_{i = 1}^p{ \bf A}_y [ i ] { \bf y} [ n - i ] + { \bf v} [ n ] \qquad &{ \rm covar} ( { \bf v} [ n ] ) = { \bf T}_1,\end{split} \tag{2}\end{align*} \end{document}

where x[n] and y[n] represent time-variant spatial fMRI information in the brain regions of X and Y, respectively, and u[n] and v[n] are the multivariate noises. Σ ₁ and T ₁ are covariance matrices of the respective noises quantifying the correctness of the prediction on current signals from past signals during the process. Similarly, it is possible that fMRI information from two regions is treated as one time-variant vector q[n] and fit to the model as follows: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{split}&{ \bf q} [ n ] = \left[ \begin{matrix}{ \bf x} [ n ] \\ { \bf y} [ n ] \end{matrix} \right]\\ &{ \bf q} [ n ] = \mathop \sum_{i = 1}^p{ \bf A}_q [ i ] { \bf q} [ n - i ] + { \bf w} [ n ] \quad { \rm covar} ( { \bf w} [ n ] ) = \left[ \begin{matrix} \bf \Sigma_2 \quad { \bf C} \\ { \bf C}^{\rm T} \quad { \bf T}_2\end{matrix} \right]\end{split} \tag{3}\end{align*} \end{document}

In this case, the covariance matrix of noise w[n] can be divided into four matrices, including Σ ₂ and T ₂·Σ ₂ and T ₂ are square matrices that have the same column sizes as x[n] and y[n], respectively, and can be considered indices of prediction accuracy in the two regions.

In the VAR, each element of current vector is estimated by the linear combination of its past signals and the past signals of the other elements. When the BOLD signals of two regions are considered together as a one time-variant vector, the prediction of each region can be supported by the information of the other region. In contrast, when the two regions are considered separately, each region can only use its own information. Therefore, the directional influence between the two regions can be detected, that is, if the prediction of current signals of x[n] is improved with the support of past signals of y[n], it implies that there is a directional influence from the region of Y to the region of X. Using the indices of the correctness of the prediction, the directional influences between the two regions are quantified as follows: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{split}& F_{x \rightarrow y} = \ln ( \parallel { \bf T}_1 \parallel / \parallel { \bf T}_2 \parallel ) \\ & F_{y \rightarrow x} = \ln ( \parallel { \bf \Sigma}_1 \parallel / \parallel { \bf \Sigma}_2 \parallel ) ,\end{split} \tag{4}\end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$F_{x \rightarrow y}$$ \end{document} represents the influence from region X to Y, and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$F_{y \rightarrow x}$$ \end{document} represents the influence from region Y to X. We focus our interest on these two measures among the functional aspects of Granger causality to show the directional flow of information.

Simulation Method

fMRI data can be considered a spatially low-pass filtered and temporally subsampled neuronal signal. We designed a simulation to validate that the proposed method can detect the spatial pattern-based causality of the brain regions with reference to the previous studies which validated the appliance of GCM to fMRI (Roebroeck et al., 2005). We generated the simulated fMRI data of two brain regions, between which there are causal interactions, and analyzed the data with the proposed method to find the existing causality. The simulated fMRI was generated with consideration of the hemodynamics of the brain.

The signals x[n] and y[n] of the two virtual brain regions of X and Y are generated from the VAR process as follows: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}{ \bf q} [ n ] = { \bf Aq} [ n - 1 ] + { \bf Bq} [ n - d ] + { \bf u} [ n ]\end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}{ \rm for} \quad { \bf q} [ n ] = \left[ \begin{matrix}{ \bf x} [ n ] \\ { \bf y} [ n ] \end{matrix} \right] , \tag{5}\end{align*} \end{document}

where x[n] and y[n] are the time-variant vectors with the size of a and b, respectively, and each element of noise u[n] follows a normal distribution. Matrix A implies the self-induced influence of each region, matrix B indicates the directional influence from x to y, and d is the variable for the delay of the influence. The matrices A and B are given in this simulation as follows:

where the size of the matrices A and B are (a+b)×(a+b). The bottom-left b×a part of the matrix B induces the influence from x to y, and composed with elements named ɛ, which are random variables that follow a uniform distribution between −I and I, when I is the value that indicates the strength of the influence.

Using this model, we can generate a simulated signal that has an influence from X to Y while adjusting the various parameters such as a and b (the numbers of voxels in the regions X and Y), I (the strength of the influence), and d (the delay of the influence). In this model, there is no influence from Y to X. The autoregressive coefficients on the diagonal of the matrix A were set to make the influence exist in the interested low-frequency ranges, as the hemodynamic process acts as a low-pass filter and the signal is subsampled in the time domain. The elements named ɛ in the matrix B were set to be random to exert an influence with complex interactions among the voxels in both regions.

The virtual neuronal signal of each voxel was generated with an equivalent length of 11 min with a time interval of 10 msec. The first 1 min signal was discarded for allowing steady states. Then, each resulting signal was filtered by convolution with the hemodynamic function (HRF) to simulate the BOLD signal, where the canonical HRF supplied by SPM8 (www.fil.ion.ucl.ac.uk/spm/software/spm8) was applied. After the convolution, each signal was normalized to zero mean and unit variance, and an additional 20% of white Gaussian noise was added to simulate the psychological noise in the BOLD signal. Finally, each resulted signal was subsampled with an equivalent interval of TR for simulating the data acquisition. The signal was also renormalized and another white Gaussian noise was added as the acquisition noise in the MRI system. The resulting signals that represent the fMRI data acquired from two regions X and Y were processed by the proposed method to measure \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$F_{x \rightarrow y}$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$F_{y \rightarrow x}$$ \end{document} . The order of the VAR model was set to minimize the Schwartz criterion (SC), which is suggested to reduce error variance while preventing overfitting (Luetkepohl, 1991).

We designed three sets of the simulations to evaluate the performance of the proposed method in different situations, as there are several parameters for the simulation, such as the number of voxels composing pattern-based causality (a and b), the number of voxels in each region used in analysis (V), the strength of the influence (I), the delay of the influence (D), the sampling interval of the scanner (TR), and the acquisition noise (E). In the first set of the simulation, we performed 1000 trials with random elements of B and using parameters (a=b=10, V=5, I=0.3, D=60 msec, TR=1 sec and E=20%) to examine the distributions of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$F_{x \rightarrow y}$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$F_{y \rightarrow x}$$ \end{document} . The additional 1000 trials without influence (B=0) were also performed to determine the significance threshold of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$F_{x \rightarrow y}$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$F_{y \rightarrow x}$$ \end{document} . Second, we performed different cases using combinations of the parameters (a=b=10, V={1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, I={0.1, 0.2, 0.3, 0.4, 0.5}, D={20, 40, 60, 80, 100 msec}, TR={1, 2, 3 sec}, and E={10, 20, 30, 40, 50%}), where 100 trials were performed for each of the 3750 possible combinations (i.e., 10(V)×5(I)×5(D)×3(TR)×5(E)=3750). Furthermore, the same number of additional trials, excluding influence (B=0), were performed to consider the false-positive ratio of the proposed method. In the third set of the simulation, 1000 times of trials were performed with different numbers of voxels to verify the performance of the proposed method according to the relationship between the size of regions exerting the influence and the size of each ROI used in analysis (a=b={1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, V={1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, I=0.3, D=60 msec, TR=1 sec and E=20%). When the analyzing ROI was larger than the influential region, simulated time-course data without influence were added for the analysis.

Simulation Result

The distributions of computed \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$F_{x \rightarrow y}$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$F_{y \rightarrow x}$$ \end{document} from the first set of simulation with the representative parameters are shown in Figure 1. The significance threshold can be determined using classical type I error control without any directional influence (α=0.001). The first observation is that \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$F_{x \rightarrow y}$$ \end{document} is larger than \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$F_{y \rightarrow x}$$ \end{document} on average, showing the intended true influence from X to Y. On the other hand, the average of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$F_{y \rightarrow x}$$ \end{document} is significantly larger than zero, even though there was no influence from Y to X. This is a reflection of the instantaneous influence between X and Y, which is induced because the Granger causality is applied to the subsampled time series. When the sampling interval is longer than the actual interaction that takes place, the instantaneous influence increases according to the sampling interval. Therefore, the interpretation of the influence based on individual \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$F_{x \rightarrow y}$$ \end{document} or \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$F_{y \rightarrow x}$$ \end{document} is not appropriate. The possible solution is an appliance of the differential Granger causality (Roebroeck et al., 2005), which analyzes the inference based on the difference of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$F_{x \rightarrow y}$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$F_{y \rightarrow x}$$ \end{document} . If the influence difference \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$F_{x \rightarrow y} - F_{y \rightarrow x}$$ \end{document} is positive, it can be inferred that there is a directional influence from X to Y. We observed the distribution of the difference to confirm that the differential Granger causality is appropriate for the proposed method, and the result is shown in Figure 2. The significance threshold is determined using the classical type I error control (α=0.001), and the difference is significantly positive in all trials of this simulation. The result reflects good sensitivity and specificity of the differential Granger causality. Thus, we adopted this approach for the proposed method.

FIG. 1.

Distribution of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$F_{x \rightarrow y}$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$F_{y \rightarrow x}$$ \end{document} in simulation with representative parameters: five voxels in region of interest (ROI), the influence strength of 0.3, the influence delay of 60 msec, sampling interval of 1 sec for the magnetic resonance imaging, and the acquisition noise of 20%. The distribution is caused by a random factor in the vector autoregressive model and noise in the process. Reference distribution was calculated in case of no influence between regions. Dotted line indicates the threshold controlling a fraction of type I error to be 0.001.

FIG. 2.

Distribution of difference between two measures, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$F_{x \rightarrow y} - F_{y \rightarrow x}$$ \end{document} . Reference indicates the distribution of the difference in case of no influence. Dotted line indicates the threshold controlling a fraction of type I error to be 0.001.

To examine the performance of the proposed method engaging in various parameters, we performed a z-test with hypothesis \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$F_{x \rightarrow y} - F_{y \rightarrow x} = 0$$ \end{document} and the results are shown in Figure 3. The performance improved as the strength of influence and delay of influence increased, and the TR and the acquisition noise decreased. These observations agree with the previous simulation of conventional GCM (Roebroeck et al., 2005). The performance also improved as the size of the ROI increases, and becomes significantly reliable (p<0.001) for a region larger than 6 voxels. To examine the quantified performance corresponding to the number of voxels in the ROI used in analysis, the receiver operating characteristic (ROC) curve was calculated, which shows the relationship between the true-positive ratio and false-positive ratio. The ROC curve for each case is shown in Figure 4. The performance improved according to the number of voxels in the ROI, and was almost saturated after eight voxels in the ROI.

FIG. 3.

Results of the simulation showing detecting performance of the proposed method according to the number of voxels in ROI and other parameters: (a) strength of the influence, (b) delay of the influence, (c) sampling interval of the scanner, and (d) acquisition noise. Measured z-statistics corresponding to parameters are averaged.

FIG. 4.

Receiver operating characteristic curve of the proposed method corresponding to the number of voxels in the ROI used in analysis. The results show better performance as the size of ROI becomes larger. The highest line showing the best performance indicates the case of 10, which is equal to the size of regions that contain causal information, and the lowest line indicates the case of one showing a similar performance with the reference case using the average time-course data of each region.

Finally, Figure 5 shows the performance according to the relationship between the number of voxels exerting the influence and the number of voxels in ROI used in analyzing. The overall performance is high when the influence is distributed in many voxels in the ROI. When the influential pattern size is fixed, the performance is improved as the analyzing ROI size becomes larger and saturated at the point that the analyzing ROI size becomes the influential pattern size. When the analyzing ROI size becomes larger than the influential pattern size, the performance is slightly reduced and the performance is practically maintained.

FIG. 5.

Simulation results showing the performance of the proposed method according to the relationship between the number of voxels exerting influence and analyzing ROI.

The results of the simulation suggest that the proposed method can be used to detect the pattern-based Granger causality over multiple voxels after low-pass filtering and subsampling. It was found that the performance is robust to the strength or delay of the influence, MRI acquisition noise, and TR in the case of a large number of voxels in the pattern.

Experimental Method

We applied the proposed method to actual fMRI experiments to verify the existence of the pattern-based causality in the brain and evaluated the usefulness of the proposed method for understanding brain functions. The experiment design in the previous research for GCM (Roebroeck et al., 2005) was also adopted to verify the proposed method, because the design was known to create a causal relationship between the brain regions and had desired properties. The design imported a rapid event-related paradigm. Subjects performed a visuomotor task of clicking buttons that corresponded to the displayed visual stimuli. The visual stimuli included three categories: face, house, and another object. The face and the house stimuli required a response of clicking the left hand and right hand buttons, whereas another object stimulus required no response. Since there were two responses of the left and right hand click for the two visual stimuli of face and house display, there were two possible stimulus-response mapping (S-R mapping). The S-R mapping was occasionally switched in the experiment. The face visual stimulus was known to activate the fusiform face area, and the clicking buttons caused activation in specific parts of the motor cortex. Using the rapid event-related paradigm, it could be expected that the successful performance of the task created a flow of information passing through these regions. Moreover, it was presumed that contextual information was exploited in the task, and the influences between areas involved in the control processes were persistently sustained. Therefore, the analysis would tolerate relatively low sample rates, and the whole time-course data could be analyzed to detect the causal relationship created by the task.

The experiment was composed of a slow-switching (SS) condition and fast-switching (FS) condition, and each subject participated in the experiment of both conditions. In the SS condition, the S-R mapping was switched periodically and less frequently, while it was switched randomly and frequently in the FS condition. The FS condition had a more engaging context and required a higher degree of executive control than the SS condition. Therefore, the interaction between areas was assumed to be different. If the measured influence was different according to the conditions, we could confirm that the results of the proposed method was not caused by various hemodynamic properties over the brain, as hemodynamics persist over the conditions.

The experiment was performed on 20 subjects who were right handed and had normal vision. All subjects were trained before the experiment and signed a written consent form. The images were acquired using a 3 T MRI system (Siemens Verio), a T2*-weighted echo planar imaging sequence with an echo time (TE) of 28 msec and a repetition time (TR) of 1000 msec. The image matrix size was 64×64, and the field of view was 224×224 mm². Eighteen oblique axial slices were scanned for each volume with a thickness of 5 mm and a splice gap of 1 mm. 530 volumes were scanned for both SS and FS conditions. Anatomy images were obtained using a T1-weighted MPRAGE sequence with a TE of 28 msec, a 3D image size of 256×256×192, and a voxel size of 1×1×1 mm³. The display of stimuli and the gathering of responses were guided by the E-prime 2.0 software (Psychology Software Tools). The S-R mapping was changed 24 times at every two to six trials in the FS condition, and it was changed eight times at every 15 trials in the SS condition. The stimulus onset asynchrony was 2–6 sec for both conditions. The mapping cue was indicated by changing the color of the fixation cross to magenta for mapping 1 and cyan for mapping 2, and was displayed for 500 msec at the beginning of each trial. Trial stimuli consisted of five pictures for each face, house, and object, and were displayed for 120 msec following the mapping cue. Feedback was displayed for 500 msec at the end of each trial by changing the color of the fixation cross to green for a correct response and red for an incorrect response. Each experiment continued for 45 trials, which consisted of 15 trials for each of the face, house and object stimuli. The face, house, and object stimuli were shown in pseudorandom order.

Preprocessing was performed by SPM8, and the visualization of the experimental results was guided by Brain Voyager QX (Brain Innovation). Each volume was registered to the Colin brain template supplied by Brain Voyager, and corrected for slice timing. The FFA and the visual cortex were detected by SPM8 using a general linear model (GLM). The iterative approach was applied to FFA, and the visual cortex was applied to locate specific regions that contained a distinct influence. The regions having a causal relationship with the FFA are detected by the proposed analysis. The visual cortex that significantly affects the FFA is considered a new ROI of the visual cortex. Next, the most influenced region from the new ROI of the visual cortex is detected and treated as a new ROI of the FFA again. These procedures are repeated until the two regions are not changed. For every subject, the converging 3×3×3 cubic region of FFA was analyzed, and the region was selected as ROI. The spatial smoothing was applied only for the initial GLM analysis, whereas any spatial smoothing filter was not applied for the proposed pattern-based analysis, as it may reduce the spatial pattern information (Norman et al., 2006). Every spatiotemporal data in the 3×3×3 cubic window in the brain were analyzed to map the differential pattern-based Granger causality over the entire brain. For each region, we collected 20 samples of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$F_{x \rightarrow y}$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$F_{y \rightarrow x}$$ \end{document} from all subjects to perform a statistical analysis using a t-test with the null hypothesis \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$F_{x \rightarrow y} - F_{y \rightarrow x} = 0$$ \end{document} .

Experiment Result

The experiment results of the proposed method are shown in Figure 6. The FFA was selected as ROI and is shown in the red circle. The visual cortex indicated by the green color shows a negative influence, confirming that the visual region influences the ROI. The motor cortex indicated by the blue color shows a positive influence. Thus, the motor region is influenced by the ROI. Overall, it can be inferred that there is a relay of information flowing through the visual cortex, FFA, and the motor cortex. From this result, a specific part in the frontal robe also shows a strong positive influence. Possibly the processing system involved in the given task is not direct from the visual cortex to the motor cortex, but requires additional functions, and it can be inferred that the frontal lobe is also related to the performing of the task. The frontal robe is known to be engaged in the decision making (Andersen and Cui, 2009; Barraclough et al., 2004; Bechara et al., 2000; Heekeren et al., 2004; Lee et al., 2007), and we presume the function is required to accomplish the task successfully.

FIG. 6.

Experimental results of the proposed method applied to the functional magnetic resonance imaging experiment showing the regions that interact with ROI. The circled face-selective region was selected as ROI. A significant negative differential Granger causality (p<0.001) was detected in the green regions of the visual cortex, which exert influence on the ROI. The blue regions were analyzed to have a positive differential Granger causality (p<0.01), which received influence from the ROI.

We investigated the SS and FS conditions separately to confirm that the resulting interaction is caused by the influences between regions, not by the variability of the hemodynamic response over the brain. The results from SS and FS conditions are shown in Figure 7. Both results show the connection between the visual and motor cortex, while the shapes of the detected regions have different patterns. In particular, in the frontal lobe and the motor region, a stronger positive influence was detected in the FS condition. A statistical test was performed to show that there was a significant difference between the results of the SS and FS conditions. The null hypothesis was that there was no difference between the results of the two conditions. We adopted the method controlling the false discovery rate (FDR) (Genovese et al., 2002). The FDR-based hypothesis testing controls q, which was the expected proportion of false positives within all tests. The result showed that there was a significant difference between the results of the SS and FS conditions (q=0.05), whereas the difference was not affected by the hemodynamic responses. Therefore, we can conclude that the proposed method works as intended.

FIG. 7.

The experimental results that were analyzed separately with (a) slow-switching and (b) fast-switching conditions. The circled face-selective regions were selected as ROI. Both results show the interaction between the visual and motor cortex, while the patterns of detected regions are slightly different.

It was also needed to confirm that the analyzed influences are based on pattern-based information and not on overall changes of intensity in the regions. To exclude this confound, we analyzed the data again in the same way except for spatially averaging over voxels within each region. As a result, we confirmed that there is no significant region showing an interaction with the ROI (p<0.01), as the pattern-based influence information was removed by spatial averaging.

Discussion

We proposed the multivoxel pattern-based GCM to explore the influences between high-level spatial pattern-based information in the brain. The proposed method can be applied to fMRI data using the same procedure of GCM (Roebroeck et al., 2005), except for taking the ROI as a spatial region of the brain. The simulation results showed the ability of the proposed method in detecting the influence after the hemodynamics and MRI acquisition process, and the experiment results showed the practical application of the proposed method in analyzing the existing influence.

There have been considerable debates whether the Granger causality can reveal the underlying neuronal influence, while variation in the hemodynamic latency and subsampling exists. Some simulation studies (Smith et al., 2011) insisted that the Granger causality was determined by the confounding hemodynamic latencies. Nonetheless, many studies (Schippers et al., 2011; Wen et al., 2012) have shown different results. They concluded that the influence detected by the Granger causality was likely to be true even at the neuronal level, and Granger causality was a reliable method for analyzing influence in fMRI data. In addition, Deshpande and Hu concluded that simulations should not be treated as the final results, because they were bounded by assumptions and simplifications (Deshpande and Hu, 2012). Since we assumed that there were meaningful pattern-based interactions based on the model, we adopted the VAR model to the proposed method and validated the proposed method in detecting the pattern-based Granger causality in fMRI signal under the circumstance of hemodynamics and subsampling.

The order of the VAR model in the proposed method was determined by the SC. It is known that the optimal order of the GCM is almost exclusively one when TR is 1000 msec (Roebroeck et al., 2005), and the tendency is maintained when detecting the pattern-based information.

The shape of the window in analyzing the experimental data was set to 3×3×3 cubic. The pattern-based information in the brain is presumed to be dependent on functions, and most studies using MVPA have adopted a larger window. As shown in Figure 5, there is a considerable possibility that the accuracy of the proposed method can be improved by using a larger window. Therefore, it is encouraged to apply a large window in the proposed method. However, a required temporal length of the time-course data to fit the VAR model is exponentially increased as the size of ROI increases when we apply the conventional fitting algorithm. Thus, it can generate problems related to the concentration or fatigue of subjects. There is a tradeoff between the ROI size and the required experiment time. A possible solution to fit the VAR model with the shorter time-course is the ridge-regression approach (Ahmad et al., 2012; Douglas et al., 2001), which introduces an additional parameter in order to solve an ill-posed problem. Since the ridge regression is a regularization method and does not guarantee ideal output, further research about the tradeoff between increasing the size of the window and importing the ridge-regression can be helpful in dealing with this issue.

Since any spatial smoothing was not applied to the experimental fMRI data in order to preserve the spatial pattern-based information, the noise model of the fMRI data can be more ambiguous than the conventional fMRI data analysis. To confirm that the proposed method has robust performance against the noise model, we performed a simulation that the simulated fMRI data were generated using the noise models of the uniform, skewed-normal, and bimodal-normal distribution, respectively. All noises were applied to all three stages of exerting error in the simulation, such as generation of the neuronal signal, BOLD signal, and the fMRI data. All noises had the same mean and standard deviation. Though the noise models were different, the performances were much similar in all cases of the noise models. From the result, we concluded that the proposed method had robust performance against the noise model.

Recently, it has been reported that the GCM is robust to changes in hemodynamic response properties unless the changes are accompanied by severe downsampling (Seth et al., 2013); therefore, it can be a helpful further research that investigates the properties when the proposed method is applied.

In conclusion, we proposed a method that could explore the influence between multivoxel pattern-based information, which was known to deal with high-level function in the brain, by adopting the concept of Granger causality. The existence and the analyzability of the influence were verified through the simulations and the real fMRI experiments. The proposed method can be applied to map the pattern-based influence related to ROI over the entire brain. The proposed method is expected to contribute in improving the understanding of the dynamic features of high-level information in the brain.

Footnotes

Acknowledgment

This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MEST) (No. 2010-0018833 and No. 2009-0079457).

Author Disclosure Statement

No competing financial interests exist.

References

Abler

, Roebroeck

, Goebel

, Höse

, Schönfeldt-Lecuona

, Hole

, Walter

. 2006. Investigating directed influences between activated brain areas in a motor-response task using fMRI. Magn Reson Imaging, 24:181–185.

Ahmad

, Maqbool

, Kim

, Park

, Kim

. 2012. An efficient method for effective connectivity of brain regions. Concepts Magn Reson A, 40A:14–24.

Andersen

, Cui

. 2009. Intention, action planning, and decision making in parietal-frontal circuits. Neuron, 63:568–583.

Anderson

, Oates

. A Critique of Multi-Voxel Pattern Analysis. In Proceedings of the 32nd Annual Meeting of the Cognitive Science Society, Portland, OR, 2010, pp. 1511–1516.

Barraclough

, Conroy

, Lee

. 2004. Prefrontal cortex and decision making in a mixed-strategy game. Nat Neurosci, 7:404–410.

Bechara

, Damasio

. 2000. Emotion, decision making and the orbitofrontal cortex. Cereb Cortex, 10:295–307.

Biswal

, Yetkin

, Haughton

, Hyde

. 1995. Functional connectivity in the motor cortex of resting human brain using echo planar MRI. Magn Reson Med, 34:537–541.

Boccaletti

, Latorab

, Morenod

, Chavezf

, Hwang

. 2006. Complex networks: structure and dynamics. Phys Rep, 424:175–308.

Deshpande

, Hu

. 2012. Investigating effective brain connectivity from FMRI data: past findings and current issues with reference to granger causality analysis. Brain Connect, 2:235–245.

10.

Deshpande

, LaConte

, James

, Peltier

, Hu

. 2008. Multivariate Granger causality analysis of fMRI data. Hum Brain Mapping, 30:1361–1373.

11.

Douglas

, Elizabeth

, Geoffrey

. 2001. Introduction to Linear Regression Analysis, 3rd edition. New York: Elsevier Academic Press.

12.

Friston

. 2009. Causal modelling and brain connectivity in functional magnetic resonance imaging. PLoS Biol, 7:220–225.

13.

Gao

, Chen

, Gong

. 2008. Evaluation of the effective connectivity of the dominant primary motor cortex during bimanual movement using Granger causality. Neurosci Lett, 443:1–6.

14.

Garg

, Cecchi

, Rao

. 2011. Full-brain auto-regressive modeling (FARM) using fMRI. Neuroimage, 58:416–441.

15.

Genovese

, Lazar

, Nichols

. 2002. Thresholding of statistical maps in functional neuroimaging using the false discovery rate. Neuroimage, 15:870–878.

16.

Geweke

. 1982. Measurement of linear dependence and feedback between multiple time series. J Am Stat Assoc, 77:304–313.

17.

Goebel

, Roebroeck

, Kim

, Formisano

. 2003. Investigating directed cortical interactions in time-resolved fMRI data using vector autoregressive modeling and Granger causality mapping. Magn Reson Imaging, 21:1251–1261.

18.

Granger

CWJ

. 1969. Investigating causal relations by econometric models and cross-spectral methods. Econometrica, 37:424–438.

19.

Granger

CWJ

. 1980. Testing for causality: a personal viewpoint. J Econ Dyn Control, 2:329–352.

20.

Greicius

, Krasnow

, Reiss

, Menon

. 2003. Functional connectivity in the resting brain: A network analysis of the default mode hypothesis. Proc Natl Acad Sci U S A, 100:253–258.

21.

Hampson

, Peterson

, Skudlarski

, Gatenby

, Gore

. 2002. Detection of functional connectivity using temporal correlations in MR images. Hum Brain Mapping, 15:247–262.

22.

Haxby

, Gobbini

, Furey

, Ishai

, Schouten

, Pietrini

. 2001. Distributed and overlapping representations of faces and objects in ventral temporal cortex. Science, 293:2425–2430.

23.

Heekeren

, Marrett

, Bandettini

, Ungerleider

. 2004. A general mechanism for perceptual decision-making in the human brain. Nature, 431:859–862.

24.

Kamitani

, Tong

. 2005. Decoding the visual and subjective contents of the human brain. Nat Neurosci, 8:679–685.

25.

Kay

. 1988. Modern Spectral Estimation: Theory and Application. Englewood Cliffs, NJ: Prentice Hall.

26.

Kriegeskorte

, Kreiman

. 2012. Visual Population Codes: Towards a Common Multivariate Framework for Cell Recording and Functional Imaging. Cambridge, MA: MIT press.

27.

Lee

, Rushworth

MFS

, Walton

, Watanabe

, Sakagami

. 2007. Functional specialization of the primate frontal cortex during decision making. J Neurosci, 27:8170–8173.

28.

Lowe

, Dzemidzic

, Lurito

, Mathews

, Phillips

M. D

. 2000. Correlations in low-frequency BOLD fluctuations reflect cortico-cortical connections. Neuroimage, 12:582–587.

29.

Luetkepohl

. 1991. Introduction to Multiple Time Series Analysis. Heidelberg: Springer-Verlag.

30.

Marreiros

, Kiebel

, Friston

. 2008. Dynamic causal modelling for fMRI: a two-state model. Neuroimage, 39:269–278.

31.

Martino

, Scheres

, Margulies

, Kelly

, Uddin

, Shehzad

, Biswal

, Walters

, Castellanos

, Milham

. 2008. Functional connectivity of human striatum: a resting state FMRI study. Cereb Cortex, 18:2735–2747.

32.

Norman

, Polyn

, Detre

, Haxby

. 2006. Beyond mind-reading: Multi—voxel pattern analysis of fMRI data. Trends Cogn Sci, 10:424–430.

33.

Peelen

, Downing

. 2007. Using multi-voxel pattern analysis of fMRI data to interpret overlapping functional activations. Trends Cogn Sci, 11:4–5.

34.

Raizada

RDS

, Kriegeskorte

. 2010. Pattern-information fMRI: New questions which it opens up, and challenges which face it. Int J Imaging Syst Technol, 20:31–41.

35.

Roebroeck

, Formisano

, Goebel

. 2005. Mapping directed influence over the brain using Granger causality and fMRI. Neuroimage, 25:230–242.

36.

Schippers

, Renken

, Keysers

. 2011. The effect of intra- and inter-subject variability of hemodynamic responses on group level Granger causality analyses. Neuroimage, 57:22–36.

37.

Seth

, Chorley

, Barnett

. 2013. Granger causality analysis of fMRI BOLD signals is invariant to hemodynamic convolution but not downsampling. Neuroimage, 65:540–555.

38.

Smith

, Miller

, Salimi-Khorshidi

, Webster

, Beckmann

, Nichols

, Ramsey

, Woolrich

. 2011. Network modelling methods for fMRI. Neuroimage, 54:875–891.

39.

Tang

, Bressler

, Sylvester

, Shulman

, Corbetta

. 2012. Measuring granger causality between cortical regions from voxelwise fMRI BOLD signals with LASSO. PLoS Comput Biol, 8:e1002513.

40.

Valdés-Sosa

, Bornot-Sánchez

, Vega-Hernández

, Melie-García

, Lage-Castellanos

, Canales-Rodríguez

. 2006. Granger causality on spatial manifolds: applications to neuroimaging. In: Schelter

, Winterhalder

, Timmer

(eds.) Handbook of Time Series Analysis. Berlin: Wiley-VCH; pp. 461–485.

41.

Wen

, Rangarajan

, Ding

. 2012. Linking BOLD-level Granger Causality with Neural-Level Granger Causality. In Processings of the Annual Meeting of the Organization for Human Brain Mapping, Beijing, China, p. 579.