## Abstract

We investigated large-scale systems organization of the whole human brain using functional magnetic resonance imaging (fMRI) data acquired from healthy volunteers in a no-task or ‘resting’ state. Images were parcellated using a prior anatomical template, yielding regional mean time series for each of 90 regions (major cortical gyri and subcortical nuclei) in each subject. Significant pairwise functional connections, defined by the group mean inter-regional partial correlation matrix, were mostly either local and intrahemispheric or symmetrically interhemispheric. Low-frequency components in the time series subtended stronger inter-regional correlations than high-frequency components. Intrahemispheric connectivity was generally related to anatomical distance by an inverse square law; many symmetrical interhemispheric connections were stronger than predicted by the anatomical distance between bilaterally homologous regions. Strong interhemispheric connectivity was notably absent in data acquired from a single patient, minimally conscious following a brainstem lesion. Multivariate analysis by hierarchical clustering and multidimensional scaling consistently defined six major systems in healthy volunteers — corresponding approximately to four neocortical lobes, medial temporal lobe and subcortical nuclei — that could be further decomposed into anatomically and functionally plausible subsystems, e.g. dorsal and ventral divisions of occipital cortex. An undirected graph derived by thresholding the healthy group mean partial correlation matrix demonstrated local clustering or cliquishness of connectivity and short mean path length compatible with prior data on small world characteristics of non-human cortical anatomy. Functional MRI demonstrates a neurophysiological architecture of the normal human brain that is anatomically sensible, strongly symmetrical, disrupted by acute brain injury, subtended predominantly by low frequencies and consistent with a small world network topology.

## Introduction

Functional magnetic resonance imaging (fMRI) of the human brain reveals a complex, structured neurophysiological network, even when the person lying in the scanner has not been asked to do anything more cognitively demanding than simply to ‘rest’. For example, an fMRI time series extracted from somatosensorimotor cortex was specifically correlated with time series extracted from other components of the motor system, although the subject was not performing a motor task (Biswal *et al.*, 1995). Similar examples of resting state correlation have been reported by selective sampling of time series from one or a few brain regions, including primary visual cortex and amygdala (Lowe *et al*., 1998), hippocampus (Stein *et al*., 2000), perisylvian language centres (Hampson *et al*., 2002) and components of a putative ‘default mode network’ (Greicius *et al*., 2003). Strong, symmetric interhemispheric correlations between bilaterally homologous brain regions have been demonstrated for fusiform gyrus, basal ganglia and thalamus (Lowe *et al*., 1998; Cordes *et al*., 2002). Another consistent feature is that resting state correlations are often subtended by low-frequency components of the MR signal, typically <0.1 Hz (Biswal *et al*., 1995; Lowe *et al*., 1998; Cordes *et al*., 2000). Comparable low-frequency fluctuations in resting brain hemodynamics have been measured using other techniques such as laser Doppler flowmetry and optical imaging (Mayhew *et al*., 1996; Hudetz *et al*., 1998).

The neurobiological basis of these observations has been disputed. There are several potentially confounding sources of (co)variance in fMRI time series, including head movement and (possibly aliased) cardiac or respiratory pulsation. However, time series extracted from voxels representing blood vessels or ventricular cerebrospinal fluid (CSF) did not show the pattern of low-frequency correlation that was typical of time series extracted from diverse brain regions (Cordes *et al*., 2001; Rombouts *et al*., 2003). Moreover, using high sampling rates [repetition time (*T*_{R}) = 400 ms] to isolate cardiac and respiratory sources without aliasing, strong inter-regional correlation persisted in the low-frequency range (<0.05 Hz) and this was not attributable to physiological noise (Cordes *et al*., 2001). Further fast-sampling studies were reported to rule out a major contribution of head movement or instrumental instability (Cordes *et al*., 2002).

In addition to these experiments seeking to refute cardiorespiratory and other biologically trivial causes, there have been a number of studies more directly affirming neural mechanisms for resting state correlations. It has been shown that such correlations are stronger in BOLD than in flow-weighted echo-planar imaging data (Biswal *et al*., 1997). Several groups have shown that correlation between brain regions at rest is modulated by imposition of cognitive tasks (Lowe *et al*., 2000). For example, performance of a language task enhanced correlation between inferior frontal and superior temporal regions that were positively correlated at rest (Hampson *et al*., 2002); whereas resting state correlation between posterior cingulate cortex and other components of a default mode network was attenuated by performance of a demanding cognitive task but unaffected by simple visual stimulation (Greicius *et al*., 2003). There is also evidence that resting state correlation can be affected by nervous system disorders such as multiple sclerosis and callosal agenesis which degrade myelination or integrity of axonal tracts between correlated regions (Quigley *et al*., 2001).

Resting state correlations have often been discussed in terms of *functional connectivity*, broadly defined as the statistical association or dependency between anatomically distinct neurophysiological time series (Aertsen *et al*., 1989; Friston, 1993; Horwitz, 2003). Recent resting state fMRI studies have progressed methodologically by measuring connectivity between a larger number of brain regions *m* and by using multivariate statistical methods, such as hierarchical cluster analysis (Cordes *et al*., 2002) to summarise graphically the potentially large number of pairwise correlations or dissimilarities [(*m*^{2} − *m*)/2] that can result from a more comprehensive approach. Another methodological development has been adoption of the partial correlation coefficient to estimate specific associations between regions, factoring out the contributions to pairwise correlations that might arise due to global or third-party effects (Hampson *et al*., 2002). An analogous method in the frequency domain is the partial coherence, which provides a measure of frequency-specific association between two regions (Sun *et al*., 2004). Partial correlations (or coherences) can be used to build undirected graphs, in which connections (edges) between nodes (vertices) depict their conditional dependence (Whittaker, 1990; Stam, 2004).

Here we describe the first whole brain analysis of functional connectivity based on human fMRI data acquired at rest. A prior anatomical template image was used to extract 90 regional mean time series from each of 12 individual images acquired from healthy volunteers and the group mean inter-regional partial correlation matrix was estimated. We also acquired comparable data from a single patient with acute brain injury. We applied complementary methods of multivariate and graph theoretical analysis to characterize and visualize functional architecture of the healthy human brain, paying particular attention to the evidently important dependencies of inter-regional functional connectivity on the anatomical relationships between regions.

## Materials and Methods

### Functional MRI Data Acquisition: Sample, Scanning and Pre-processing

Twelve healthy volunteers (7 male, 5 female) with a mean age of 30 years (range 23–48 years) were recruited by advertisement locally. A single male patient in a state of minimal consciousness following acute brainstem ischaemia was also studied. Computed tomography data on the patient showed high signal intensities in the pons, midbrain and both cerebral peduncles, extending to the left subthalamic region. There was an absence of signal in the basilar artery with probable atheromatous occlusion of the left posterior cerebral artery and right middle cerebral artery.

The study was approved by the Addenbrooke's NHS Trust Local Research Ethics Committee and all healthy participants gave informed consent in writing. For the study of the patient, informed assent was obtained in writing from his next of kin. Each volunteer was scanned while lying quietly in the scanner at rest and with eyes closed.

In each scanning session, 259 gradient-echo echo-planar imaging (EPI) volumes depicting BOLD contrast were acquired using a Bruker Medspec S300 scanner operating at 3.0 T (Bruker Medical, Ettlingen, Germany) with the following acquisition parameters: number of slices = 21 (interleaved), slice thickness = 4 mm, interslice gap = 1 mm, matrix size = 64 × 64, flip angle = 90°, *T*_{R} = 1.1 s, *T*_{E} = 27.5 ms, in-plane resolution 3.125 mm. The first three volumes were discarded prior to analysis to allow for *T*_{1} saturation effects, leaving 256 volumes available for analysis.

Each image was corrected for geometrical displacements due to estimated head movement, and co-registered with the Montreal Neurological Institute (MNI) EPI template image, using SPM2 software (http://www.fil.ion.ucl.ac.uk/spm). The data were not spatially smoothed.

### Anatomical Parcellation

Regional parcellation of fMRI datasets after registration with the MNI template image was done using the anatomically labelled template image previously validated and reported by Tzourio-Mazoyer *et al*. (2002). This parcellation divides each cerebral hemisphere into 45 anatomical regions of interest (ROIs), which are listed in Table 1 together with the abbreviations used to refer to them in this study. A regional mean time series was estimated simply by averaging the fMRI time series over all voxels in each of 90 regions over the whole brain; time series were not normalized by mean subtraction before averaging in each region.

**Table 1**

Region | Abbreviation | Region | Abbreviation |
---|---|---|---|

Precentral gyrus | PreCG | Lingual gyrus | LING |

Superior frontal gyrus, dorsolateral | SFGdor | Superior occipital gyrus | SOG |

Superior frontal gyrus, orbital part | ORBsup | Middle occipital gyrus | MOG |

Middle frontal gyrus | MFG | Inferior occipital gyrus | IOG |

Middle frontal gyrus orbital part | ORBmid | Fusiform gyrus | FFG |

Inferior frontal gyrus, opercular part | IFGoperc | Postcentral gyrus | PoCG |

Inferior frontal gyrus, triangular part | IFGtriang | Superior parietal gyrus | SPG |

Inferior frontal gyrus, orbital part | ORBinf | Inferior parietal, but supramarginal and angular gyri | IPL |

Rolandic operculum | ROL | Supramarginal gyrus | SMG |

Supplementary motor area | SMA | Angular gyrus | ANG |

Olfactory cortex | OLF | Precuneus | PCUN |

Superior frontal gyrus, medial | SFGmed | Paracentral lobule | PCL |

Superior frontal gyrus, medial orbital | ORBsupmed | Caudate nucleus | CAU |

Gyrus rectus | REC | Lenticular nucleus, putamen | PUT |

Insula | INS | Lenticular nucleus, pallidum | PAL |

Anterior cingulate and paracingulate gyri | ACG | Thalamus | THA |

Median cingulate and paracingulate gyri | DCG | Heschl gyrus | HES |

Posterior cingulate gyrus | PCG | Superior temporal gyrus | STG |

Hippocampus | HIP | Temporal pole: superior temporal gyrus | TPOsup |

Parahippocampal gyrus | PHG | Middle temporal gyrus | MTG |

Amygdala | AMYG | Temporal pole: middle temporal gyrus | TPOmid |

Calcarine fissure and surrounding cortex | CAL | Inferior temporal gyrus | ITG |

Cuneus | CUN |

Region | Abbreviation | Region | Abbreviation |
---|---|---|---|

Precentral gyrus | PreCG | Lingual gyrus | LING |

Superior frontal gyrus, dorsolateral | SFGdor | Superior occipital gyrus | SOG |

Superior frontal gyrus, orbital part | ORBsup | Middle occipital gyrus | MOG |

Middle frontal gyrus | MFG | Inferior occipital gyrus | IOG |

Middle frontal gyrus orbital part | ORBmid | Fusiform gyrus | FFG |

Inferior frontal gyrus, opercular part | IFGoperc | Postcentral gyrus | PoCG |

Inferior frontal gyrus, triangular part | IFGtriang | Superior parietal gyrus | SPG |

Inferior frontal gyrus, orbital part | ORBinf | Inferior parietal, but supramarginal and angular gyri | IPL |

Rolandic operculum | ROL | Supramarginal gyrus | SMG |

Supplementary motor area | SMA | Angular gyrus | ANG |

Olfactory cortex | OLF | Precuneus | PCUN |

Superior frontal gyrus, medial | SFGmed | Paracentral lobule | PCL |

Superior frontal gyrus, medial orbital | ORBsupmed | Caudate nucleus | CAU |

Gyrus rectus | REC | Lenticular nucleus, putamen | PUT |

Insula | INS | Lenticular nucleus, pallidum | PAL |

Anterior cingulate and paracingulate gyri | ACG | Thalamus | THA |

Median cingulate and paracingulate gyri | DCG | Heschl gyrus | HES |

Posterior cingulate gyrus | PCG | Superior temporal gyrus | STG |

Hippocampus | HIP | Temporal pole: superior temporal gyrus | TPOsup |

Parahippocampal gyrus | PHG | Middle temporal gyrus | MTG |

Amygdala | AMYG | Temporal pole: middle temporal gyrus | TPOmid |

Calcarine fissure and surrounding cortex | CAL | Inferior temporal gyrus | ITG |

Cuneus | CUN |

The abbreviations listed are those used in this paper, which differ slightly from the original abbreviations by Tzourio-Mazoyer *et al*. (2002).

### Estimation of the Inter-regional Partial Correlations

Given a set of *m* random variables, the partial correlation matrix is a symmetric matrix in which each off-diagonal element is the correlation coefficient between a pair of variables after partialling out (conditioning under normality) the contributions to the pairwise correlation of all other variables included in the dataset. In this case, therefore, the partial correlation between any two brain regions partials out the effects of the 88 other brain regions defined by the template image.

To estimate the partial correlation matrix for each healthy volunteer and the patient, we used the method described in Whittaker (1990). This method is based on the fairly general assumption of multivariate normality of the observations. In the case of fMRI, however, we are dealing with a series of *t* time points and we therefore must assume jointly Gaussian stationary multivariate stochastic processes. Specifically, we will be dealing with the zero lag (instantaneous) cross-covariance and partial cross-correlation matrices, which, for simplicity, we will call covariance and partial correlation matrices.

The first step is to estimate the sample covariance matrix *S* from the (*m* = 90, *t* = 256) data matrix *Y* of observations for the individual. Each component of *S* contains the sample covariance value between two regions (say *j* and *k*)

*S*is estimated, it must be inverted and, in general, this requires time series longer than the number of regions, i.e.

*t*>

*m*. Finally, the off-diagonal elements of the inverted matrix

*S*

^{−1}must be rescaled to obtain the sample partial correlation matrix

*R*

*j,k*}th element of the inverted matrix

*S*

^{−1}, not the inverted value of the {

*j,k*}th element of matrix

*S*.

### Testing for Non-zero Partial Correlations

The partial correlation matrices individually obtained from each one of the *n* healthy volunteers (*i* = 1,…, 12) were averaged to estimate the healthy group mean inter-regional functional connectivity matrix

To test the null hypothesis that the (mean) partial correlation was zero between any pair *j,k* of regions, we conducted multiple one-sample *t*-tests, using the mean

*n*= 12). Although Fisher's

*r*-to-

*Z*transform was applied to improve normality, it had only marginal effect due to the moderate values of the observed partial correlations.

The *t*-test was performed for all 4005 possible pairs of regions and, consequently, a correction for multiple comparisons was strictly necessary. The false discovery rate (FDR) approach was applied to find a threshold that would restrict the expected proportion of type I errors to *q* < 0.05. We used the procedure described in Benjamini and Yekutieli (2001), which takes into account the lack of independence between tests. The first step involved sorting the 4005 individual *P*-values derived from the *t*-tests in ascending order

*t*tests with a

*P*-value equal or smaller than

*P*were considered significant, i.e. indicative of non-zero partial correlations.

_{k}### Frequency Dependence and Scaling Properties of Functional Connectivity

We considered the relationship between functional connectivity and frequency or scale of the time series components contributing to (partial) correlation between regions. A simple way to do this is by dividing each ‘raw’ time series into its high- and low-frequency components. To isolate the high-frequency components, we computed the discrete Fourier transform (Press *et al*., 1992) and set to zero the Fourier coefficients corresponding to frequencies <0.007 Hz; to isolate the low-frequency components we set to zero the equal number of coefficients corresponding to frequencies >0.007 Hz. In both cases the data were reconstituted in the time domain using the inverse Fourier transform before further analysis; see Figure 1 for an example.

**Figure 1.**

**Figure 1.**

This analysis can be regarded as an intermediate step towards a full analysis of functional connectivity in the frequency domain, in which an analogous quantity to the partial correlation, the partial coherence, *pcoh*, can be estimated for each of the Fourier frequencies (Brillinger, 1981; Sun *et al*., 2004):

*j,k*}th element of the inverted spectral density matrix at frequency λ (the equivalent of the covariance matrix in the frequency domain). The use of partial coherences instead of partial correlations naturally extends the concept of functional connectivity to frequency-specific coherence as a measure of association between distinct brain regions.

Another mathematical approach to the same general issue is based on the discrete wavelet transform (DWT). The DWT decomposes the variance in a time series over a hierarchy of scales (approximately frequencies) and scale-specific functional connectivity between time series can then be measured simply by estimating the covariance or correlation between the wavelet coefficients of the two time series at each scale of the decomposition. This provides a straightforward way of assessing the scaling properties of functional connectivity, i.e. the extent to which the correlation between any two time series is subtended by signal components at different scales. For greater detail on wavelet decomposition of covariance between bivariate time series, see Serroukh and Walden (2000a,b); for background on scaling and wavelet analysis of fMRI data, see Bullmore *et al*. (2004).

### Anatomical Distance and Functional Connectivity

We considered the dependency of functional connectivity on anatomical distance. Explicitly, we investigated the relationship between partial correlation *r* estimated for a pair of resting fMRI time series and the Euclidean distance *D* (in mm) between regional centroids in Talairach space:

*x*,

_{j}*y*,

_{j}*z*) and (

_{j}*x*,

_{k}*y*,

_{k}*z*) are the Talairach coordinates of the centroids for regions

_{k}*j*and

*k*.

### Multivariate Analyses of the Healthy Group Mean Partial Correlation Matrix

The mean partial correlation matrix for the group of healthy volunteers was transformed to a matrix of dissimilarities by calculating the absolute values of the partial correlations and subtracting them from 1,

This shift from partial correlations to dissimilarities allowed us to apply two complementary multivariate techniques to exploratory visualization of the positive semidefinite dissimilarity matrix: hierarchical cluster analysis with an average linkage agglomerative algorithm (Krzanowski, 1988; see also Cordes *et al*., 2002, and Stanberry *et al*., 2003, for examples of prior applications of hierarchical cluster analysis in fMRI) and classical or metric multidimensional scaling (MDS) (Gower, 1966; Everitt, 1995; see also Friston *et al*., 1996, and Welchew *et al*., 2002, for prior applications of MDS to PET and fMRI data analysis).

### Graph Theoretical Analysis

After averaging time series acquired from bilaterally homologous regions in each healthy volunteer, the unihemispheric partial correlation matrix (45 × 45) was probabilistically thresholded so that each statistically significant connection could be represented as an undirected edge between regional vertices in a graphic rendering of the whole brain functional network. For arbitrary size of the probability threshold *P* applied to the partial correlation matrix, local density or cliquishness of connections between vertices (regions) was estimated by the clustering coefficient C_{P} (Watts and Strogatz, 1998; Strogatz, 2001)

*m*regions

*j*: 1,…, 45) of a ratio. Its numerator is the number of edges connecting neighbours of region

*j*(#

*E*) where the neighbourhood of a region is defined as the set of all vertices that are directly connected to that region by an edge. The denominator of the ratio is the maximum possible number of connections between neighbours of vertex

_{j}*j*(where #

*V*is the number of neighbours of vertex

_{j}*j*). For a regular or lattice network,

*C*∼ 0.75; for random networks, C ∼

*v*/

*m (v*being the average of #

*V*and

_{j}*m*being the total number of vertices). The minimum path length,

*L*, is the average of the shortest path length over each possible pair of vertices:

_{P}*lpaths(i*,

*j*)} is the set of lengths of all possible paths connecting vertices

*i*and

*j*through one or more edges, and the path length is defined as the number of edges included in the path.

For a regular network, mean *L _{P} ∼ m/2v*; for a random network,

*L*∼ ln

_{P}*m*/ln

*v*. Thus regular networks have high cliquishness but relatively long path lengths, whereas random networks have low cliquishness but short path lengths. So-called ‘small world’ networks have intermediate properties, combining the short path lengths typical of random networks with a higher degree of local clustering than is typical of random networks.

## Results

### Inter-regional Partial Correlations

The elements of the healthy group mean partial correlation matrix were individually tested for significance and 76 pairs of regions had non-zero partial correlations, as listed in Table 2 (*P* < 0.05, FDR). Most significant connections (44; 58%) were intrahemispheric and often local, involving regions in the same lobe and/or closely adjacent to each other anatomically. Interhemispherically symmetric connections, involving bilaterally homologous regions in the two hemispheres, were also relatively numerous (29; 38%). Interhemispherically asymmetric connections (involving non-homologous regions in different hemispheres) were infrequent (3; 4%).

**Table 2**

The strength of correlation between bilaterally homologous time series sampled from a single healthy individual is evident by inspection of Figure 1. Fourier decomposition of these data into high- and low-frequency components exemplifies the prior generalization that resting state correlations are stronger at low frequencies. Wavelet analysis shows an approximately linear increase in the log of the variance of the wavelet coefficients at each scale as a function of increasing scale (decreasing frequency). This implies there is 1/*f*-like power law scaling in both resting time series, i.e. the spectrum *G* is proportional to a power law function of frequency *f:* G ∼ *f*^{α}, with α = −1.52 for left superior temporal gyrus (STG) and −1.44 for right STG. There was also an approximately linear increase in the log of the covariance between wavelet coefficients estimated from right and left STG data as a function of increasing scale (decreasing frequency), implying that power law scaling also applies to inter-regional functional connectivity.

### Relationships Between Functional Connectivity and Anatomical Distance

As shown in the upper panel of Figure 2, shorter anatomical distances (*D* < 2 cm) between centroids generally predicted stronger functional connectivity between regions in the healthy volunteers. Most possible connections at anatomical distances >4 cm had functional connectivity randomly variable around zero. The form of this non-linear relationship can be described approximately by an inverse square law: *r* ∼1/*D*^{2}. A similar relationship was observed in the patient's data, although in the context of much greater variability (Fig. 2, middle panel)

**Figure 2.**

**Figure 2.**

However, it was also notable that, for the healthy volunteers, there were a number of functional connections between regions that were much stronger than would be predicted as a function of the anatomical distance between them. These connections are defined in Figure 2 as greater than the 95% quantile of partial correlations estimated in each 1 cm bin of the distribution of Euclidean distances; these unusually strong, mostly long-distance connections are listed in Table 2. Most of them (28/36) were interhemispheric connections between bilaterally homologous regions; the remainder were intrahemispheric connections between left frontal regions (4/36), between left inferior parietal lobule and left middle or inferior frontal gyrus (2/36), and between left superior parietal gyrus and left inferior parietal lobule or inferior temporal gyrus (2/36). This pattern was entirely lost in the data acquired from the minimally conscious patient — who showed no evidence of unusually strong, long-range interhemispheric or left intrahemispheric connectivity.

When the analysis of whole brain data from the healthy volunteers is repeated after dividing each time series into high- and low-frequency components, it is clear that partial correlations consistently tend to be larger when estimated on the basis of the low-frequency components of the time series. This effect is most clearly seen for the symmetric interhemispheric connections (Fig. 2, lower panel).

### Hierarchical Cluster Analysis

The dendrogram resulting from hierarchical cluster analysis of the healthy group mean partial correlation matrix is shown in Figure 3*a*. The basic hierarchy of brain functional organization broadly respected anatomical relations between regions and could be heuristically designated lobar/sublobar/symmetrical. In other words, symmetrical links between bilaterally homologous cortical and subcortical regions were consistently represented at the lowest level of the hierarchy.

**Figure 3.**

**Figure 3.**

At the highest level of the hierarchy, six large clusters were identified, corresponding to four neocortical/paralimbic systems, one paralimbic/limbic system and one subcortical/limbic system. One of the largest neocortical/paralimbic systems comprised 22 regions of parietal, posterior cingulate and (pre)motor cortices that are specialized for motor and spatial attentional functions (Mesulam, 2000). Another large neocortical/paralimbic system comprised 22 regions of orbitofrontal, prefrontal and anterior cingulate cortices that have been implicated in diverse strategic or executive functions (Duncan and Owen, 2000). A third system comprised 14 regions of the occipital lobe specialized for visual processing. The fourth neocortical/paralimbic system comprised 12 regions of lateral temporal cortex and insula specialized for auditory-verbal functions. The paralimbic/limbic cluster comprised 10 regions of medial temporal lobe and temporal pole that have been implicated in mnemonic and affective processing. The predominantly subcortical cluster comprised 10 regions including striatum, pallidum, thalamus and olfactory cortex.

At intermediate levels of the hierarchy, symmetrical pairs of homologous regions were clustered together in configurations corresponding approximately to sub-lobar or gyral domains; the intermediate hierarchical structure of each of the six main systems is represented in anatomical space in Figure 3*b*.

### Metric Multidimensional Scaling

The functional distances between regions in the healthy group mean data were approximated by the graphical distances between them in the two-dimensional space of Figure 4: regions that are functionally similar are plotted in close proximity. This analysis confirms many of the organizational features already highlighted: symmetrical regions are often paired in the same neighbourhood of the space and the overall configuration broadly respects anatomical relations between regions.

**Figure 4.**

**Figure 4.**

Superimposed on the MDS plot are lines corresponding to the significant pairwise inter-regional partial correlations listed in Table 2. This again highlights the predominance of local intrahemispheric and symmetric interhemispheric connections.

### Small World Properties

We thresholded the unihemispheric partial correlation matrix of the healthy volunteers so that any partial correlation with *P* < 0.05 was represented by an edge between the corresponding regional vertices (and all other possible edges were set to zero). For this network, the clustering coefficient *C*_{P} was 0.25 and the mean minimum path length *L*_{P} was 2.82. Corresponding parameters for a random graph with the same number of nodes were:

## Discussion

We expected whole human brain functional architecture to be strongly conditioned by the anatomical relationships between regions and our results clearly endorse this expectation. Correlational and classical multivariate methods of analysis, applied to these neurophysiological time series data without any explicit incorporation of prior anatomical knowledge, consistently demonstrated an anatomically sensible pattern of functional organization. This is exemplified variously by the approximately lobar and sub-lobar hierarchy of the dendrogram (Fig. 3); the dorsal and ventral divergence of visual cortical regions illustrated by MDS (Fig. 4) and predicted by prior anatomical and functional neuroimaging studies (McIntosh *et al*., 1994; Ungerleider and Haxby, 1994); and the inverse square law relating anatomical distance to intrahemispheric functional connectivity (Fig. 2). Anatomy did not always predict precisely the functional relationships between regions: anterior and posterior parts of cingulate cortex had weaker functional connectivity than their anatomical contiguity might predict, and several interhemispheric and left intrahemispheric pairs of regions had stronger long-range connectivity than the distance between them would predict; but in general functional connectivity obeyed anatomical constraints.

Many of these results are compatible with prior studies of cat or monkey cortical networks defined by multivariate analyses of anatomical connectivity matrices. For example, Young (1992, 1993) used non-metric MDS to demonstrate dorsal and ventral ‘streams’ of inter-regional connectivity, a preponderance of local neighbourhood connections, and a generally hierarchical organization of primate visual cortex. Scannell *et al*. (1995, 1999) applied non-metric MDS, optimal set analysis and non-parametric cluster analysis to whole brain connectivity matrices derived from anatomical studies of the cat, and reported four major hierarchically organized systems: visual, auditory, somatosensory/motor and frontal/limbic. Hilgetag *et al*. (2000) used similar methods to demonstrate a broadly ventral–dorsal dichotomy of visual cortex, and a cluster comprising somatosensory and motor cortices, by analysis of both cat and macaque anatomical connectivity matrices.

Our data additionally highlight an aspect of functional connectivity that has not been so clearly predicted by these elegant analyses of anatomical connectivity, namely the importance of bilaterally symmetric interhemispheric connections. One reason for this discrepancy is simply that the anatomical connectivity matrices on which multivariate analyses have previously been based are unihemispheric, summarizing inter-regional connections within a single (right or left) hemisphere. Prior studies of resting state fMRI data have reported strong correlations between diverse bilaterally homologous regions, and demonstrated that interhemispheric functional connectivity depends on the integrity of the corpus callosum (Lowe *et al*., 1998; Quigley *et al*., 2003). Pioneering work by Horwitz *et al*. (1984), using a partial correlational analysis of whole brain positron emission tomography data on regional glucose metabolism, had previously demonstrated strong correlations between all bilaterally homologous regions. However, we believe that our results provide the first confirmation by fMRI of the ubiquity of symmetric inter-regional connectivity as an organizing principle of normal brain functional architecture; and they highlight other long-range intrahemispheric connections, e.g. between left prefrontal and inferior parietal cortices, that are stronger than would be predicted by the anatomical distance between regions. In a single patient, minimally conscious as a result of brainstem ischaemia, we have also shown that relatively long-range, often symmetrical connectivity was specifically attenuated, whereas short-range connectivity was less affected and tended to fall off as an inverse square function of anatomical distance in the same way as seen in the healthy volunteers. These limited clinical data confirm previous reports indicating that resting state connectivity may be a potentially useful marker of brain disease or damage (Lowe *et al*., 1998; Quigley *et al*., 2003); they also illustrate the feasibility of this approach to measuring brain function in an individual patient who would be unable to co-operate in a cognitively demanding experiment. It will be interesting in future work to explore the diagnostic and prognostic value of abnormal long-range functional connectivity in patients with acute brain injury.

The existence of many strong connections between closely neighbouring brain regions, taken together with evidence for fewer, longer-range connections, raised the possibility that human brain functional architecture might have ‘small world’ properties. Small world networks have been studied intensively since the seminal work by Watts and Strogatz (1998) demonstrated that random mutation of a few connections in a regular or lattice network substantially reduced the mean minimum path length between any pair of nodes while retaining strong local connectivity or ‘cliquishness’. Small world properties have been described for diverse social, biochemical, computational and ecological networks; they are often associated with scale-free topology, simple growth rules, economic wiring between nodes and non-linear dynamics (see Strogatz, 2001, for a recent review). The distinctive combination of high cliquishness and short minimum path length has previously been shown for non-human cortical networks derived from anatomical tract-tracing studies (Hilgetag *et al*., 2000) and studies of seizure propagation induced by direct cortical application of strychnine (Stephan *et al*., 2000); synthetic graphs selected by the criterion of high complexity of their dynamical behavior (Sporns *et al.*, 2000, 2002); and graphs derived from measures of synchronization between multiple (126) magnetoencephalography (MEG) channels of ‘resting’ human brain dynamics (Stam, 2004). Our results demonstrate for the first time that such properties may be characteristic also of large-scale human brain functional networks derived from whole brain fMRI.

Our measure of ‘small worldness’ is essentially the difference in the ratios

*et al*. (2000) reported comparable results: for macaque visual cortex

*C*

_{P}in the brain functional network are somewhat small (0.25) compared both to the limiting case of the perfectly regular network, which has

*C*

_{P}∼ 0.75, and the anatomical networks studied by Hilgetag

*et al*. (2000), which have

*C*

_{P}∼ 0.6. Moreover, we found empirically that both

*C*

_{P}and

*L*

_{P}converged quite rapidly with their expected random values as the probability threshold for significance of a partial correlation was relaxed (

*P*> 0.05), or as the number of participants in the sample (or the number of time points in each fMRI time series) was reduced. In short, we think it is likely that large-scale human brain functional networks measured by fMRI have small world characteristics, but a more conclusive and comprehensive investigation of this issue would benefit from longer time series measured in more people. On this basis it would be interesting also to explore the fMRI analog of Stam's (2004) observation in MEG data that small world properties of human brain networks are more salient in connectivity matrices derived from high- and low-frequency band-pass filtered time series than in matrices derived from intermediate frequency bands. Our wavelet decomposition of the covariance between two regional time series (Fig. 2) provides some preliminary evidence that fMRI connectivity scales log-linearly with decreasing frequency of the time series components subtending partial correlations between regions, but this observation needs to be generalized to multiresolution analysis of whole brain graphs based on longer fMRI time series.

Another important methodological consideration is the issue of anatomical parcellation. One key general advantage of parcellating fMRI datasets prior to multivariate analysis is that it is necessary to have the number of anatomical loci smaller than the number of time points, i.e. *m* < *t*, for inversion of the inter-regional covariance matrix to estimate the partial correlation matrix. This would be impossible to achieve, without exceptionally long scan times, if each voxel was treated as a separate locus. Regional averaging of time series also effects a degree of denoising and ensures that anatomical labelling is defined according to some standard, previously validated template incorporating expert neuroanatomical knowledge. However, there are many different schemes for anatomical parcellation and the template we have adopted (Tzourio-Mazoyer *et al*., 2002) is not the only one available for this purpose. Future studies of whole brain functional connectivity might usefully explore the impact of different parcellation schemes on neurophysiological architecture. It will be useful also to explore the impact of other pre-processing steps that may be implicit in anatomical parcellation, such as the effects on spatial covariance between neighbouring time series caused by interpolation and realignment of the observed time series volume to match a parcellated template image in standard space.

Can our results be explained in terms of nuisance sources of inter-regional covariance such as head movement or cardiorespiratory cycle-related pulsation? It may seem unlikely *a priori* that a global factor such as head movement could generate a set of regionally specific connectivities that conform to local anatomical relationships. And our use of the partial correlation, rather than Pearson's correlation, as a measure of functional connectivity will have attenuated the contribution of global sources of covariance to the functional connectivity between any pair of regions. Moreover, we have used a standard technique to correct head movement by co-registering each image volume with the first volume in each experimental series. We have also explored the impact of more draconian corrections for head movement, e.g. regressing the co-registered fMRI time series on the time series of estimated translations and rotations of the image center of gravity, and found that this did not make a major difference to the functional configurations demonstrated by MDS or hierarchical cluster analysis (data not shown). As noted earlier, the most successful approach to eliminate the possible effects of cardiac and respiratory cycles on functional connectivity has been to acquire images at high frequencies (*T*_{R} < 1 s) and low-pass filter the resulting time series to isolate low-frequency components (subtending functional connectivity) from the higher-frequency peaks at cardiac and respiratory cycle frequencies (Cordes *et al*., 2001, 2002; Rombouts *et al*., 2003). We did not seek to emulate this experimental approach because, at least using the MR system available to us locally, such rapid sampling implies loss of anatomical coverage and we wished to consider the functional architecture of the whole brain rather than just a few slices of it.

However, we suggest that our results linking functional connectivity in no-task fMRI data to human anatomical constraints, and the affirmative evaluation of these results in relation to prior studies of non-human anatomical connectivity, strongly implies that there are neurobiological mechanisms underpinning ‘resting’ state connectivity in human fMRI. We have replicated the prior observation that these connections are predominantly subtended by low-frequency components of the time series, and this was especially true of long-distance interhemispheric and left intrahemispheric connections that are presumably mediated by the corpus callosum and other major white matter tracts. In keeping with this interpretation, it has been shown that callosal agenesis is asssociated with loss of functional connectivity between contralaterally homologous regions of sensorimotor and auditory cortex (Quigley *et al*., 2003). There is also supportive anatomical data from non-human primate studies, e.g. demonstrating symmetry of striatal projections from left and right prefrontal and premotor regions and between left and right premotor regions (McGuire *et al*., 1991a,b). Alternative anatomical substrates for co-ordinating activity between bilaterally homologous regions are the ascending arousal systems (Robbins and Everitt, 1995). Changing levels of noradrenergic, serotonergic, cholinergic and dopaminergic input from ascending arousal systems could result in varying degrees of local neural activation and coupled cerebral blood flow (CBF) change, or affect regional CBF directly through neurotransmitter release in the microvasculature. The distribution of post-synaptic neurotransmitter receptors that serve such systems is not uniform (Zilles *et al*., 2002), and it is possible that the summed effect of varying levels of individual arousal systems may have varying effects on different parts of the brain. However, it is highly likely that these effects would be symmetrical, thus explaining the interhemispheric connectivity seen. The absence of symmetrical connectivity demonstrated in data acquired from a patient with brainstem ischaemia is arguably compatible with a role for ascending transmitter systems in maintaining coherent low-frequency oscillations in bilaterally homologous cortical regions.

## Conclusions

Functional MRI demonstrates a neurophysiological architecture of the human brain that is anatomically sensible, strongly symmetrical, subtended predominantly by low-frequency time series components and compatible with an underlying small world topology. Further work is needed to define the generative mechanisms of this architecture, to explore its scaling behavior, and to assess its utility as a marker of state changes in brain function due, for example, to drugs or diseases.

This neuroinformatics research was supported by a Human Brain Project grant from the National Institute of Biomedical Imaging & Bioengineering and the National Institute of Mental Health. The Wolfson Brain Imaging Centre is supported by a Co-operative Group grant from the Medical Research Council (UK).

## References

## Author notes

Brain Mapping Unit1 and Wolfson Brain Imaging Centre2, University of Cambridge, Addenbrooke's Hospital, Cambridge CB2 2QQ, UK