Instantaneous and Causal Connectivity in Resting State Brain Networks Derived from Functional MRI Data

Image Acquisition

Twenty six healthy volunteers (13 male, 13 female with age 15.1 ± 1.1 years) participated in this study. The study was approved by the Institutional Review Board at Emory University, and informed consent was obtained prior to scanning. The subjects were instructed to gaze at a fixation cross on the screen and not to think of anything in particular. Resting state data were acquired on a 3T Siemens Trio scanner (Siemens Medical Solutions, Erlangen, Germany). A T₂*-weighted echo planar imaging sequence was used with the following parameters: matrix of 64×64, 20 axial slices without gap and with thickness of 4 mm, flip angle of 90°, FOV of 192 cm, 210 time points, TR/TE of 2000 ms/30 ms, and a total scan time of 7 minutes for each subject.

Pre-processing and Identification of Networks

“Analysis of Functional NeuroImages” (AFNI) software package (http://afni.nimh.nih.gov/afni/) was used for pre-processing the data and extracting the time series from the ROIs. Pre-processing steps including slice-timing correction, rigid body registration, 0.009 Hz < f < 0.08 Hz band-pass filtering and 5 mm full width half max Gaussian smoothing were carried out on the data. In order to extract the time series corresponding to the four networks – DMN, HCMN, DAN, FPCN – co-ordinates of ROIs corresponding to these networks available in published literature were used. First, the ROIs corresponding to DMN – posterior cingulate cortex (PCC), left/right posterior inferior parietal lobule (L/R pIPL), orbitofrontal cortex/ventral anterior cingulate cortex (OFC/vACC), dorsomedial prefrontal cortex Brodmann area 8 (dMPFC BA8), dorsomedial prefrontal cortex Brodmann area 9 (dMPFC BA9), left dorsolateral prefrontal cortex (L DLPFC), left parahippocampal gyrus (L PHG) and left inferolateral temporal cortex (L ITC) – were defined according to Greicius' paper (Grecius, Krasnow, Reiss, & Menon, 2003). Since R DLPFC, R PHG and R ITC have been reported as part of the DMN in many subsequent studies, we included those ROIs within the DMN. Subsequently, the ROIs corresponding to HCMN, DAN and FPCN were demarcated based on co-ordinates reported by Vincent et al (Vincent, Kahn, Snyder, Raichle, & Buckner, 2008). Accordingly, DAN consisted of left /right middle temporal area (L/R MT), left/right frontal eye fields (L/R FEF) and left/right superior parietal lobule (L/R SPL), HCMN was composed of left/right hippocampal formation (L/R HF), ventromedial prefrontal cortex (vmPFC), PCC, and bilateral pIPL, and FPCN consisted of bilateral anterior prefrontal cortex (aPFC), dorsal anterior cingulate cortex (dACC), left/right dorsolateral prefrontal cortex (L/R DLPFC), left/right anterior insula (L/R aINS) and left/right anterior inferior parietal lobule (L/R aIPL). Here, we distinguished the regions common to HCMN and DMN by using the corresponding network names as their suffix (for example, PCC DMN and PCC HCMN). Table 1 lists the Talairach coordinates of the ROIs of all 4 networks (Talairach & Tournoux, 1988); it can be seen that the coordinates of regions with the same name but appearing in multiple networks (PCC and pIPL in DMN and HCMN and L DLPFC in DMN and FPCN) are different. Consequently, we considered those ROIs separately as part of distinct networks without merging them. We hypothesized that if their functional roles are similar because of anatomical proximity, then the functional similarity will be reflected in the results of our analysis.

ROIs were 2 cubic cm spheres centered at coordinates provided by the above mentioned studies and were masked by a template consisting of only those voxels which were inside the brain. Consequently, the sizes of some of the ROIs differed slightly. The coordinates given by Vincent et al (Vincent, Kahn, Snyder, Raichle, & Buckner, 2008) were transformed from MNI (Montreal Neurological Institute) space (Friston, Holmes, & Ashburner, 1999) (where they were originally identified) to Talairach space (Talairach & Tournoux, 1988) using a sub-function of the “GingerALE” software package (http://brainmap.org/ale/) called “icbm2tal,” an algorithm which has the ability to transform coordinates between SPM (Friston, Holmes, & Ashburner, 1999) (Statistical Parametric Mapping software package) MNI space and AFNI Talairach space (Turkeltaub, Eden, Jones, & Zeffiro, 2002). Four ROIs had an overlap. They were L DLPFC of DMN and FPCN (overlap of 403 voxels out of 1313 and 1623 voxels, respectively), L pIPL of DMN and HCMN (overlap of 516 voxels out of 1631 and 1935 voxels, respectively), PCC of DMN and HCMN (overlap of 14 voxels out of 1935 voxels for both ROIs) and R pIPL of DMN and HCMN (overlap of 384 voxels out of 1805 and 1935 voxels, respectively). Average time courses for the ROIs were generated for each ROI and subject. For the four ROIs which had an overlap, the average time courses were computed from the non-overlapped regions. Subsequently, for every ROI, the time series corresponding to each subject was concatenated to obtain one time series per ROI (Deshpande, Hu, Stilla, & Sathian, 2008; Deshpande, LaConte, James, Peltier, & Hu, 2009; Stilla, Deshpande, Laconte, Hu, & Sathian, 2007; Stilla, Hanna, Mariola, Deshpande, Hu, & Sathian, 2008; Deshpande, Hu, Lacey, Stilla, & Sathian, 2010).

Connectivity Analysis

Suppose x_n, n=1 .. k corresponds to the k selected ROI time series and X(t) =(x₁(t),x₂(t)… x_k(t))^T , then the traditional multivariate autoregressive (MVAR) model with model parameters A(n) of order p is given by

where E(t) is the vector corresponding to the residual error. As shown previously (Kaminski, Ding, Truccolo, & Bressler, 2001), direct causal relationship between the k selected ROI time series can be inferred from the model parameters as follows.

where a_ij are the elements of the matrix A and D_ij corresponds to the direct causal influence exerted from ROI j to ROI i. In experimental fMRI data, we can expect some interaction between the instantaneous and time-lagged relationships, in part due to the blurring of the neuronal activity with the sluggish hemodynamic response, resulting in “leakage” of the zero-lag effects into time-lagged estimates. Consequently, instantaneous correlation between time series can contribute to the causality obtained from the VAR model in Eq.(1). In order to alleviate this potential confound, we introduced the zero-lag term into Eq. (1) and modeled the effects of instantaneous correlation as shown in Eq. (3) below, which is termed as modified MVAR (or mMVAR) (Deshpande, Sathian, & Hu, 2010)

where the diagonal elements of A′(0) are zero such that only the instantaneous cross-correlation, and not the auto-correlation, between the time series are modeled. Here it is to be noted that A′(1) … A′ (p) ≠ A(1) … A(p). The causal relationship obtained from A′(1) … A′ (p) is purged of correlation leakage effects and is termed as CPGC while the off-diagonal elements of A′(0) correspond to zero-lag correlation.

The time series corresponding to all 30 ROIs of the four networks were input into a fifth order mVAR model (Deshpande, Sathian, & Hu, 2010) to obtain the instantaneous and causal connectivities between them. The order of the mVAR model was determined using Bayesian Information Criterion (BIC) (Schwartz, 1978). Surrogate data were obtained from the ROI time series by randomizing their phase, but retaining the magnitude spectrum, and input into the mVAR model. The above procedure was repeated 10,000 times to derive empirical null distributions for each path. The statistical significance of the corresponding paths was obtained by comparing the CPGC value obtained from the original data with the empirical null distributions (Deshpande, LaConte, James, Peltier, & Hu, 2009; Theiler, Eubank, Longtin, Galdrikian, & Farmer, 1992). If any pair of ROIs significantly influenced each other (p<0.01) in both directions, the two ROIs were considered bi-directionally connected. On the other hand, if only one of the ROIs in a pair had a significant causal interaction on the other but not vice versa, then the former was deemed to have a unidirectional causal influence on the latter. Finally, a graphical analysis of the thresholded network was carried out to ascertain strongly connected components within the network using Gabow's algorithm (Gabow, 2000) such that within each strongly connected component, each of the ROIs are mutually accessible to each other, taking into account the directionality of the paths.

Representation of Networks

Given the fact that even a simple haptic sensory task can involve a large number of connections in a complicated network (Deshpande, Hu, Stilla, & Sathian, 2008), the number of connections underlying four networks in the resting state is likely to be very large. Therefore rendition of all possible types of paths between all the ROIs is untrivial. In order to overcome this problem, we used a force-directed layout to represent the instantaneous correlation between the ROIs and separately overlaid the bidirectional and unidirectional causal paths, respectively, over the force-directed layout.

Force-directed Layout using Spring Embedding Algorithm

A spring embedding algorithm, based on the principle of using alternative forces of attraction and repulsion to determine the locations of data sources, was used to position the nodes in our networks on a 2-dimensional plane according to the similarity between them (Ebbels, Buxton, & Jones, 2006). With this approach, first introduced by Fruchterman and Reingold (Fruchterman & Reingold, 1991), we considered the ROI time series as the data sources and their instantaneous correlation as the similarity information driving the attractive and repulsive forces. Accordingly, the 30 nodes representing the ROIs were initially placed randomly on a 2D plane. At each iteration, the attractive force between the nodes was proportional to their zero-lag correlation and the repulsive force was determined by whether they were significantly connected or not. Therefore, the attractive force was a continuum represented by the significant instantaneous correlations between ROIs and the repulsive force was binary, represented by whether the ROIs were significantly connected or not. This situation can be visualized as having significantly correlated ROIs attached through springs whose spring constant was proportionally to the zero-lag correlation between them and no springs between ROIs which were not significantly correlated. Consequently, the movement of the nodes on the 2D plane proceeded as per Newtonian dynamics where the nodes represented the mass and the connections between them represented the spring, until all the nodes reached static equilibrium with the lowest potential energy. Every 100 iterations, the Euler solver was used to reorient the layout such that edge crossings were eliminated (Ebbels, Buxton, & Jones, 2006).

Functional Connectivity

The force-directed layout of the instantaneous correlation network is shown in Fig.1. Different networks have been represented in different colors and the significant (p<0.01) paths between the ROIs as green lines. The ROIs (or nodes) themselves are shown as yellow blobs. On the force-directed layout, the proximity of two ROIs indicates that they have a high instantaneous correlation between them, not in isolation, but in relation to their correlation with other nodes. Also, two nodes which are linked by a path have a strong direct connection between them. Consequently, the greater the number of nodes one needs to traverse to reach another node, the weaker is the instantaneous correlation between them. Therefore, the force-directed layout is an intuitive visual representation of the similarity information between the ROI time series. Such representations have been previously used to visualize brain networks (Buckner, Andrews-Hanna, & Schacter, 2008) as well as gene expression data (Schroeder, Gilbert, van Helden, & Noy, 2001) and social networks (Heer & Boyd, 2005).

Several observations can be made from Fig.1. First, ROIs which are anatomically close, and hence presumed to be functionally similar, but identified as being part of different networks, are tightly coupled. For example, L DLPFC of DMN and FPCN, and bilateral IPL and PCC of DMN and HCMN are not only proximal on the layout, but also have direct connections between them. As far as FC is concerned, the different methods employed to identify these regions and the consequent differences in their peak coordinates may not be of much consequence. In fact, DMN and HCMN may have similar functional roles or may even represent the same network as suspected before (Buckner, Andrews-Hanna, & Schacter, 2008) except for the fact that PHG and HF, which are a part of HCMN, are not identified as part of the DMN in some studies employing a PCC-based seed correlation map (Grecius, Krasnow, Reiss, & Menon, 2003). Nevertheless, it can be seen that PGH-HF complex is strongly linked to PCC HCMN in Fig.1.

Second, different types of groupings are evident in Fig.1. The ROIs of DAN are grouped together and away from the other networks. This is in line with the expectation that DAN, being active during externally directed cognition, is anti-correlated with DMN and HCMN which are active during internally directed cognition (Fox, Snyder, Vincent, Corbetta, Van Essen, & Raichle, 2005; Fransson, 2005). On the other hand, DMN and HCMN form the middle strand of the layout though not as tightly coupled as DAN. Lastly, FPCN is the least coupled network though parts of it are tightly coupled.

Third, PCC-IPL-PFC-ACC seems to be at the center of the layout and hence correlationally equidistant from most ROIs. It is interesting to note that these four regions have been consistently identified as part of DMN irrespective of the methodology or modality used (see Fig.3 in (Buckner, Andrews-Hanna, & Schacter, 2008) where these regions were identified in the convergence of all approaches) and PCC, IPL and PFC were identified as the major hubs of DMN in a previous FC-based study (Buckner, Andrews-Hanna, & Schacter, 2008). The prominent role for DMN in Fig.1 may be driven by the fact that this is a resting state study in which we expect DMN to be most active.

Effective Connectivity

Figs. 2 and 3 show significant bi-directional and unidirectional (p<0.01) causal paths, respectively, overlaid on the force-directed layout in Fig.1. The overlay helps to visualize EC within the perspectives and constraints imposed by FC. The most significant observation from the overlay is that causal connections exist between ROIs which are distant on the layout; i.e., ROIs which have a strong causal connection are not necessarily strongly connected by instantaneous correlation. This observation justifies the use of CPGC which accounted for the influence of zero-lag correlation on the derived causality and vice versa. Furthermore, EC reveals communication between ROIs whose activity may not be synchronized and are anatomically distant.

Graph analysis revealed that the causal network shown in Fig. 2 (i.e. the bi-directional causal paths represented by blue lines) was the only strongly connected component within which each ROI was causally accessible by every other ROI. However, the directional path between any given pair of ROIs had to traverse through either PCC or L pIPL. This reaffirms what is evident through visual inspection of Fig. 2, i.e. PCC and L pIPL acted as main hubs of the entire network. Given the centrality of these ROIs on the force-directed layout, it is to be expected that, functionally, they are ideally situated to act as hubs for information exchange with other ROIs. Their importance within the DMN has been reported before in the context of high resting state metabolism exhibited by them (Gusnard & Raichle, 2001) and their pivotal role in self-consciousness and self-referential processing (Bucknor & Carroll, 2007). In fact, Fransson et al came to exactly the same conclusion with their partial correlation FC analysis of DMN (Fransson & Marrelec, The precuneus/posterior cingulate cortex plays a pivotal role in the default mode network: Evidence from a partial correlation network analysis, 2008), i.e. PCC and pIPL are strongly connected to each other and to other parts of the DMN. Our results demonstrate that PCC and pIPL acted as transit hubs not only within the DMN but also for other resting state networks such as HCMN, DAN and FPCN. Due to the centrality of these ROIs on the force-directed layout, a seed-based correlation analysis using a PCC seed will only reveal DMN and HCMN ROIs which are closer to it in Fig.1 but not distant ROIs such as the ones in DAN and the frontal ROIs of FPCN. This may have led to the identification of different ROIs in different networks. This division may be valid to a certain extent because it has been corroborated by task-based studies (Fransson & Marrelec, The precuneus/posterior cingulate cortex plays a pivotal role in the default mode network: Evidence from a partial correlation network analysis, 2008; Sridharan, Levitin, & Menon, 2008). However, our results suggest that the functional segregation of different networks is supported by a soft boundary with extensive interaction between different networks as evidenced by the fact that PCC and L pIPL exhibited bidirectional causal interactions with most of the ROIs of DAN, frontal regions of FPCN and medial temporal regions of HCMN, even though they are not instantaneously correlated.

Fig. 3 shows the overlay of statistically significant (p<0.01) unidirectional influences on the force-directed layout of Fig. 1 where the width of the lines is proportional to the strength of the corresponding path. The total number of inputs/outputs and their corresponding magnitudes for each ROI are shown in Table 2. Most of the unidirectional inputs are directed at right HF, notably from SPL and MT in the DAN, IPL areas of DMN, HCMN and FPCN and medial and lateral prefrontal areas of DMN and FPCN. R aPFC is driven by bilateral insula, PCC and MT. The data shown in Table 2 also reiterates the fact that HF is driven the most, followed by pIPL, R aPFC and PCC. Since both PCC and L pIPL act as transit hubs with bi-directional paths, they have both inputs and outputs as listed in Table 2

Zero-lag Correlation and Unidirectional & Bi-directional Causality

Before we delve into the interpretation of our results, a brief discussion on the intuitive meanings of zero-lag correlation and unidirectional & bi-directional causality is apt. Zero-lag, unidirectional and bi-directional causal relationships between any given time series A and B can be intuitively explained in terms of their respective phases. The zero-lag term represents the in-phase relationship between A and B, i.e. the tendency of A and B to rise or fall together. In the VAR framework, this is the predictability of B at instant t from A at the same instant t and vice versa. A unidirectional causal relationship A→B represents a phase difference, i.e. the tendency of A to rise or fall before B does. In the VAR framework, this is the predictability of the future of B from the past of A but not vice versa. A bi-directional causal relationship A↔B is a more complex phenomenon and may arise as a result of the following scenarios. First, there may be exclusive intervals of A→B and B→A unidirectional relationships which may be inferred as a bi-directional path when the entire time series is input into the model. Second, the relationship between A and B may be a truly feedback process, i.e. A at time t may be able to predict B at time t+τ, but B at time t+τ may have a phase similar to B at time t−Δ and B at time t−Δ may help predict A at time t. The condition that B at time t+τ may have a phase similar to B at time t−Δ ensures that the future does not predict the present or past and hence does not violate the principle of causality. However, this condition also means that there is some quasi-periodicity in the time series. Distinguishing the above two scenarios requires Granger causality to be calculated as a function of time, in which case the first scenario will be easily discernable and the absence of which may imply the second scenario. While it is very much possible to compute Granger causality as a function of time (Sato, Junior, Takahashi, Felix, Brammer, & Morettin, 2006), it is beyond the scope of this work and will be investigated in future. In the third scenario, temporal aggregation of causality, i.e. if the causal events at the neuronal level occur at a much finer temporal scale which is then smoothed by the hemodynamic response, may cause a unidirectional influence to appear as being bi-directional (Wei, 1982). This can be alleviated by the correlation-purged version of Granger causality used here (Deshpande, Sathian, & Hu, 2010). However, as discussed later, evidence points to the fact that the neuronal correlates of resting state fMRI fluctuations are likely evolving at a temporal scale much coarser than the sampling period and hence we believe that the third scenario, which can be a potential confounding factor in task-based fMRI, is not relevant in resting state fMRI.

Resting State Processes

Let us first revisit the concept of resting state and the cognitive processes underlying it so that the above results can be interpreted in that context. Resting state has been often characterized as a “task-independent” state, which is not entirely accurate. In resting state, there is not a controlled task performed by the subject, however, while “mind wandering” during resting state, many different uncontrolled tasks are carried out by the brain. In that sense, resting state is really not a “task-free” state, but rather is a state of uncontrolled, but limited number of, tasks. Typically, in neuroimaging experiments designed to study a single brain process, the tasks are designed such that only the particular brain process of interest underlying the task is examined while controlling for the other processes which are not of interest. However, in resting state, multiple processes can coexist in an uncontrolled manner. For example, the perceptual experience of being inside the scanner and undergoing the experiment may be continuously encoded into the subjects' memory while simultaneously retrieving past memories which may help to put the current experience of being inside the scanner into context. In addition, cognitive integration of multiple processes is required to give the subject a unified sense of experience. Therefore, in the following subsections, our findings are discussed in the context of resting state processes.

Memory Encoding

The inputs into the hippocampal formation from the DAN ROIs – L MT, R MT and L SPL – may reflect memory encoding. Previous evidence suggests that encoding success correlates with the activity in the task-positive network (Kim, Daselaar, & Cabeza, 2010), which we call the DAN. Given that encoding success involves attention to external stimuli, this is in agreement with the role of the DAN in external-oriented processing. However, this reminds us that some amount of external oriented processing is inevitable in a conscious human being and although “resting state” is mainly supposed to involve internal oriented processing, it is not exclusive.

Memory Retrieval

Memory retrieval essentially consists of the following components – internal generation of mnemonic representations, retrieval of memories corresponding to those representations and integration of retrieved information into a form that is accessible to decision-making processes. The inputs to HF from medial frontal, posterior midline and inferior parietal regions of the task negative networks (DMN and HCMN) may convey the internal representations generated by them to the hippocampus. This is in agreement with the previously documented role of the task negative network in retrieval success (Kim, Daselaar, & Cabeza, 2010). The inputs to HF from anterior and posterior IPL may also be a result of attention to those internal mnemonic representations of interest which will be eventually retrieved. This is supported by the Attention-to-Memory hypothesis (Wagner, Shannon, Kahn, & Bucknor, 2005; Cabeza, Ciaramelli, Olson, & Moscovitch, 2008). The outputs from HF, which carry the retrieved information, were primarily directed at inferior parietal and posterior midline regions. Both the output-buffer hypothesis and mnemonic-accumulator hypothesis posit such a scenario (Wagner, Shannon, Kahn, & Bucknor, 2005; Cabeza, Ciaramelli, Olson, & Moscovitch, 2008) and are supported by the prominent anatomical connectivity between HF and parietal cortex (Kahn, Andrews-Hanna, Vincent, Snyder, & Buckner, 2008). The former hypothesis postulates that parietal regions contribute to holding retrieved information in a form that is accessible to decision-making processes while the latter hypothesis postulates that parietal regions integrate a memory-strength signal to summarize information coming from other brain regions. In either case, L pIPL and PCC are ideally situated to perform such a function in view of a large number of bi-directional causal relationships they possess with other regions (Fig.2).

Cognitive Integration

Right aPFC, the seed of FPCN, which receives most unidirectional inputs after HF, is at the apex of the control hierarchy (Vincent, Kahn, Snyder, Raichle, & Buckner, 2008) and has been implicated in integrating the outcomes of multiple cognitive operations (Ramnani & Owen, 2004) and acting as a cognitive buffer (Koechlin, Basso, Pietrini, Panzer, & Grafman, 1999). While the density of cell bodies is markedly lower in aPFC, its spine density and number of dendritic spines per cell have been reported to be the highest in the cortex (Jabobs, et al., 2001). This has led to the speculation that the primary function of aPFC is to integrate and assimilate inputs from other areas (Ramnani & Owen, 2004). Christoff and Gabrieli have suggested that aPFC may be specialized in the processing of internal states (Christoff & Gabrieli, 2000) given the fact that it was activated only during internalized problem solving in tasks such as “Tower of London Test” (Shallice, 1982) and not when the tasks demanded externalized problem solving or reasoning (Baker, et al., 1996). These previous findings regarding aPFC make it a probable candidate for a resting state role of integrating the information from the external environment with stored internal representations (Vincent, Kahn, Snyder, Raichle, & Buckner, 2008). Below, we discuss the possible roles of different inputs to aPFC with respect to types of information they bring into aPFC for cognitive integration.

We interpret the input to right aPFC from PCC as that which brings the internal representations from memory (note that PCC was bi-directionally connected to HF and the right part of aPFC has been specifically implicated in episodic memory processing (Rugg, Fletcher, Allan, Frith, Frackowiak, & Dolan, 1998; Tulving, 1983)) while the input from bilateral MT in the DAN brings the information about the external environment for integration at aPFC. Our findings regarding aPFC also support the hypothesis proposed by Vincent et al (Vincent, Kahn, Snyder, Raichle, & Buckner, 2008) using resting state fMRI and subsequently confirmed in task-based studies (Spreng, Stevens, Chamberlain, Gilmore, & Schacter, 2010), that the network seeded in aPFC (FPCN) acts as a control network. This role of integrating the external environment with stored internal representations has been sometimes ascribed to dMPFC (Buckner, Andrews-Hanna, & Schacter, 2008). However, it can be seen that dMPFC and right aPFC are close to each other in Fig.1 (and hence, were functionally connected) and aPFC has been hypothesized to process information in ways similar to other prefrontal regions such as dMPFC but in a more abstract form (Ramnani & Owen, 2004).

The strongest driver in terms of unidirectional paths in the entire network is bilateral aINS. Most of the output from aINS is directed towards aPFC. It is to be noted that both regions are part of the same network (FPCN), though, unlike previous studies based on seed-based correlation (Vincent, Kahn, Snyder, Raichle, & Buckner, 2008), instantaneous correlation between them derived from our mVAR model was not significant. This may be attributed to the fact that seed-based correlation does not take the effects of causality into account. Previous studies have elucidated the role of aINS in interoceptive awareness of homeostatic states (Craig, 2002; Critchley, Wiens, Rotshtein, Ohman, & Dolan, 2004) and reported its sensitivity to the degree of cognitive, homeostatic and emotional salience (Damasio, 2000; Naqvi, Rudrauf, Damasio, & Bechara, 2007). This information is critical to how one “feels,” more objectively termed “body ownership” (Tsakiris, Hesse, Boy, Haggard, & Fink, 2007) or “experience of agency” (Farrer & Frith, 2002), which is highlighted during resting state in the absence of an attention demanding task. The strong influence from aINS to aPFC suggests that “body ownership” information is relayed onto aPFC from aINS, to be integrated with the information received from other parts of the brain into aPFC as discussed before, in order to create the experience of resting state and the cognitive control of the underlying processes. In a recent study, Sridharan et al (Sridharan, Levitin, & Menon, 2008) showed that the fronto-insular cortex and to a certain extent the ACC drove prefrontal and cingulate regions of the DMN and parietal regions of the DAN, in an effort to achieve cognitive control both during tasks and resting state. The present findings can be reconciled with their results by the ensuing considerations. First, since they did not consider the aPFC region in their network, their analysis might have led to the conclusion that aINS to be driving other regions which ultimately drive aPFC. In fact, excluding the aPFC from our network led to many outputs from aINS being directed at other regions. Also, excluding both aPFC and aINS from our network led to ACC being the primary driver in the network, which is in agreement with Sridharan et al's (Sridharan, Levitin, & Menon, 2008) results that after aINS, ACC was the second prominent driver. Second, they employed a bivariate implementation of traditional Granger causality while we have used multivariate CPGC in our analysis; the former may have led to a somewhat incomplete picture.

Caveats

It is highly likely that networks other than the four considered in this study may have a significant bearing on our understanding of the resting state. This caveat can be alleviated only by a whole brain voxel-by-voxel multivariate analysis, which is computationally formidable when considering both FC and EC due to current methodological constraints, primarily in EC analysis. Therefore, we adopted the strategy of using networks identified with whole brain bivariate FC analysis as a template to carry out simultaneous multivariate FC and EC analysis.

Since we have interpreted our results in terms of underlying neuronal processes, it is prudent to discuss the link between neuronal processes and resting state fMRI since the latter is not a direct measure of neuronal activity. A strong consensus is yet to emerge regarding the neuronal correlates of slow fluctuations in resting state fMRI. However, there is compelling evidence from simultaneous electrophysiology-fMRI monkey experiments which suggests that slow fMRI fluctuations at rest are driven by corresponding slow neuronal fluctuations, i.e. variations in the envelope of neuronal activity (Leopold, Murayama, & Logothetis, 2003; Nir, et al., 2008; Shmuel & Leopold, 2008). Similar evidence has emerged recently in human literature as well, in the context of slow cortical potentials (He, Snyder, Zempel, Smyth, & Raichle, 2008). Therefore, even though we know that neural dynamics corresponding to sensory, motor and cognitive processes evolves rapidly, i.e. less than a second, there is evidence to show that slower neuronal fluctuations encode certain aspects of these processes (He & Raichle, 2009). Taken together, the body of evidence suggests that slow neuronal fluctuations and the corresponding fMRI fluctuations, which evolve over seconds, could be a useful indicator of the underlying neuronal processes we wish to study. Notwithstanding the above discussion, this study only illustrates preliminary evidence that a causal structure exists for resting state fMRI fluctuations which may subserve underlying neuronal functions and further studies are required to confirm and develop this idea.

Recently, questions have been raised regarding the applicability of Granger causality to fMRI in view of the spatial variability of the hemodynamic response (David, et al., 2008). However, in the study by David and colleagues (David, et al., 2008), the hemodynamic response functions (HRFs) of one of the regions showed an extreme discrepancy by peaking 5 s after, and returning to baseline more than 30 s after, the HRFs of other regions, owing to a pathological condition. This led to the wrong inference of directionality, arising due an underlying neuronal delay of 10 ms, from Granger causality analysis of fMRI data without hemodynamic deconvolution. After this study, Friston elucidated the demerits of Granger causality analysis of fMRI data (Friston, 2009) while other researchers have argued about its merits (Roebroeck, Formisano, & Goebel, in press; Roebroeck, Formisano, & Goebel, in press; Bressler & Seth, in press). Our own view is that the work by David and colleagues is not a good reference point for this discussion because: (1) their HRF variability was not in the normal physiological range (Handwerker, Ollinger, & D'Esposito, 2004), (2) their rat model of epilepsy did not represent a healthy brain, and (3) their neuronal delay of 10 ms was in the lower end of the range of expected neuronal delays. Therefore, we carried out extensive simulations by extending the simulations framework adopted by Roebroeck and colleagues (Roebroeck, Formisano, & Goebel, 2005) to reflect a more realistic scenario where in the neuronal delay and HRF variability adequately sampled their normal physiological ranges. Our results reported elsewhere (Deshpande, Sathian, & Hu, Effect of hemodynamic variability on Granger causality analysis of fMRI, 2010) suggested that the accuracy of Granger causality was well above chance and more than 80% with low noise levels and faster sampling. Further, in the worst case scenario of the hemodynamic delay opposing the neuronal delay, the simulations indicated that causality in the range of hundreds of milliseconds could be statistically inferred. Given the fact that even fast neuronal processes such as visual perception elicits an evoked potential in the primary visual area with a latency of 90-100 ms (Russo, Martínez, Sereno, Pitzalis, & Hillyard, 2001), it is quite likely that more complex cognitive processes with latencies which are hundreds of milliseconds may be detectable by Granger causality. Recent experimental evidence also supports this notion (Katwal, Gatenby, Gore, & Rogers, 2009; Rogers, Katwal, Morgan, Asplund, & Gore, 2010, in press). As discussed before, resting state fMRI most likely represents slow neuronal oscillations which evolve over seconds and hence the mismatch between neuronal and fMRI temporal resolutions is even less of an issue while investigating slow resting state BOLD fluctuations. Therefore, the application of Granger causality to resting state fMRI data is valid.

Summary

In this study, we have used functional and effective connectivity of brain networks to investigate the neural correlates of brain processes during resting state. The major findings of our study can be summarized as follows. First, based on FC, DMN forms the core resting state network around which the other networks amalgamate. Second, information exchange and transfer based on EC reveal extensive communication between brain regions whose activity may not be synchronous. Third, PCC and L pIPL in DMN act as major transit hubs for bi-directional causal exchange of information with most of the ROIs of other networks, facilitating the interaction between instantaneously dissociated ROIs, such as dMPFC and HF, for constraining the mental simulations in the former with past experiences in the latter. Fourth, HF and aPFC receive major inputs from other ROIs, likely reflecting memory encoding/retrieval and cognitive integration, respectively. Fifth, the magnitude of unidirectional output is the highest from aINS, which is mostly directed towards aPFC and may reflect cognitive, homeostatic and emotional salience information being relayed to it for integration. Finally, this study demonstrates the utility of our approach for investigating the neural correlates of brain processes in resting state.