## Abstract

Small-world properties have been demonstrated for many complex networks. Here, we applied the discrete wavelet transform to functional magnetic resonance imaging (fMRI) time series, acquired from healthy volunteers in the resting state, to estimate frequency-dependent correlation matrices characterizing functional connectivity between 90 cortical and subcortical regions. After thresholding the wavelet correlation matrices to create undirected graphs of brain functional networks, we found a small-world topology of sparse connections most salient in the low-frequency interval 0.03–0.06 Hz. Global mean path length (2.49) was approximately equivalent to a comparable random network, whereas clustering (0.53) was two times greater; similar parameters have been reported for the network of anatomical connections in the macaque cortex. The human functional network was dominated by a neocortical core of highly connected hubs and had an exponentially truncated power law degree distribution. Hubs included recently evolved regions of the heteromodal association cortex, with long-distance connections to other regions, and more cliquishly connected regions of the unimodal association and primary cortices; paralimbic and limbic regions were topologically more peripheral. The network was more resilient to targeted attack on its hubs than a comparable scale-free network, but about equally resilient to random error. We conclude that correlated, low-frequency oscillations in human fMRI data have a small-world architecture that probably reflects underlying anatomical connectivity of the cortex. Because the major hubs of this network are critical for cognition, its slow dynamics could provide a physiological substrate for segregated and distributed information processing.

## Introduction

A remarkable variety of social, economic, and biological networks demonstrate “small-world” properties (Strogatz, 2001), meaning the nodes of the network have greater local interconnectivity or cliquishness than a random network, but the minimum path length between any pair of nodes is smaller than would be expected in a regular network or lattice (Watts and Strogatz, 1998). Some small-world networks, such as the worldwide web or the Hollywood network of movie actors, include a small number of “hubs,” nodes with an unusually large number of connections (or large degree) (Barabási and Albert, 1999; Barabási, 2003). The degree distribution of the worldwide web obeys a power law, implying no meaningful “average” number of links to each site, and for this reason, it has been described as a scale-free network (Barabási and Albert, 1999). Many other small-world networks, including Hollywood, have exponential or exponentially truncated power law distributions, implying relatively reduced probabilities of huge hubs (Amaral et al., 2000; Albert and Barabási, 2002).

Small-worlds are attractive models for connectivity of nervous systems because the combination of high clustering and short path length confers a capability for both specialized or modular processing in local neighborhoods and distributed or integrated processing over the entire network (Sporns et al., 2004). It has been demonstrated that the neuronal network of *Caenorhabditis elegans* has small-world topology at a microscopic anatomical scale (Watts and Strogatz, 1998). Anatomical connectivity matrices from tract-tracing studies of cat and monkey cortices show small-world properties at a macroscopic scale (Hilgetag et al., 2000). There is evidence also for small-world properties of graphs inferred from functional connectivity matrices measured at a macroscopic (regional) scale in monkey and human neurophysiological data (Stephan et al., 2000; Stam, 2004; Salvador et al., 2005b) and at a mesoscopic (voxel) scale in human functional magnetic resonance imaging (fMRI) data (Eguíluz et al., 2005).

Here, we used wavelets to decompose the pairwise correlations between fMRI time series measured in multiple cortical and subcortical human brain regions in normal human volunteers scanned during “rest” or a no-task condition. Analysis of restingstate data has the advantage that it may focus attention on endogenous or background neurophysiological processes; however, it is disadvantaged by incomplete understanding of the generative mechanisms and cognitive significance of endogenous, correlated oscillations in fMRI data (we will return to this issue in the Discussion). Wavelet analysis of these data results in a set of inter-regional correlation matrices, each matrix describing the functional connectivity between regions that is subtended by signal components at a distinct scale or frequency interval (Table 1). We thresholded the scale-specific wavelet correlation matrices to estimate small-world properties of the resulting undirected graphs, and we explored in detail their statistical and anatomical properties. This allows us for the first time to evaluate the frequency dependency and degree distribution of a small-world functional network in the entire human brain, to identify and contextualize the hubs of this network anatomically, and to test the resilience of the network to random error and targeted attack.

## Materials and Methods

*Sample.* Five healthy human volunteers (age, 25–35 years; two males) were recruited by advertisement. Volunteers had no personal history of neurologic or psychiatric disorder, were not taking medication, and were not abusing alcohol or illicit drugs. All participants gave informed consent in writing. The study was approved by the Addenbrooke's NHS Trust Local Research Ethics Committee (Cambridge, UK).

*fMRI data acquisition.* Each participant was scanned on a single occasion, lying quietly at rest with eyes closed for 37 min, 44 s. Gradient-echo echoplanar imaging (EPI) data depicting blood oxygen level-dependent contrast were acquired using a Medspec S300 scanner (Bruker Medical, Ettlingen, Germany) operating at 3.0T in the Wolfson Brain Imaging Centre (Cambridge, UK). We acquired 2058 volumes with the following parameters: number of slices, 21 (interleaved); slice thickness, 4 mm; interslice gap, 1 mm; matrix size, 64 × 64; flip angle, 90°; repetition time (TR), 1100 ms; echo time, 27.5 ms; in-plane resolution, 3.125 mm. The first 10 volumes were discarded to allow for T1 saturation effects, leaving 2048 volumes available for analysis of resting state connectivity in each subject.

*fMRI data preprocessing.* Each data set was corrected initially for geometrical displacements because of estimated head movement and coregistered with the MNI (Montreal Neurological Institute) EPI template image, using an affine transform implemented in SPM2 software (www.fil.ion.ucl.ac.uk/spm). The data were not spatially smoothed before regional parcellation using the anatomically labeled template image validated previously by Tzourio-Mazoyer et al. (2002). This parcellation divides each cerebral hemisphere into 45 anatomical regions of interest, which are listed in Table 2 together with the abbreviated regional labels used in this study. Regional mean time series were estimated for each individual by averaging the fMRI time series over all voxels in each of 90 regions. Each regional mean time series was further corrected for the effects of head movement by regression on the time series of translations and rotations of the head estimated in the course of initial movement correction by image realignment. The residuals of these regressions constituted the set of regional mean time series used for wavelet correlation analysis.

*Wavelet correlation analysis.* Wavelet transforms effect a time-scale decomposition that partitions the total energy of a signal over a set of compactly supported basis functions, or little waves, each of which is uniquely scaled in frequency and located in time. Wavelet analysis is particularly well suited to analysis of signals that have fractal scaling or 1/*f* properties, as is typical of cortical fMRI time series in the resting state (Maxim et al., 2005) [see Bullmore et al. (2004) for a review of wavelet methods for fMRI data analysis and Percival and Walden (2000) for a general text on wavelet analysis of time series].

Here, we use the wavelet transform as a basis to estimate the scale-dependent correlations between each of 4005 possible pairs of the 90 cortical and subcortical fMRI time series extracted from each individual data set. A scale- or frequency-specific analysis of multivariate neurophysiological time series is motivated by previous observations from electroencephalographic studies that coherence between electrodes is not equal at all frequencies and anatomically distinct systems of brain regions may be most coherent at different frequencies. A comparable observation from previous fMRI studies of functional connectivity (Horwitz, 2003) is that resting-state correlations are predominantly subtended by very low-frequency (<0.1 Hz) signal components (Biswal et al., 1995; Lowe et al., 1998; Cordes et al., 2000), and long-distance connections (e.g., between the prefrontal and posterior parietal cortices) are especially dependent on low-frequency coherence (Salvador et al., 2005a).

We applied the maximal overlap discrete wavelet transform (MODWT; see Appendix) to each regional mean time series and estimated the pairwise inter-regional correlations between wavelet coefficients at each of the first six scales. This resulted in a set of six {90 × 90} inter-regional wavelet correlation matrices for each subject, which were averaged over subjects at each scale to produce six group mean wavelet correlation matrices (Whitcher et al., 2000). These matrices can be understood as representing the inter-regional functional connectivity that is subtended by time series components in the frequency bands defined approximately by scales 1–6. For example, the scale 1 wavelet correlation matrix comprises the functional connectivity subtended by the relatively high frequency interval 0.23–0.45 Hz, and the scale 6 matrix comprises the functional connectivity subtended by the relatively low-frequency interval 0.007–0.01 Hz (Table 1).

*Small-world analysis of wavelet correlation matrices.* Small-world properties have been described for graphs defined by a number of nodes or vertices, *n*, each connected by a number of undirected edges, *k*, to other nodes in the network. The key metrics are the degree *k*, the clustering coefficient *C*, and the mean minimum path length *L*, all defined for each of the *i* = 1,2,3,...,*n* nodes in the network. The clustering coefficient 0 < *C*_{i} < 1 is a ratio that defines the proportion of possible connections that actually exist between the nearest neighbors of a node (Watts and Strogatz, 1998); high values of *C*_{i} imply that most of the nearest neighbors of that node are also nearest neighbors of each other, or that the node is located in a cliquish local neighborhood. The minimum path length is the number of edges comprising the shortest path between any pair of nodes; the mean minimum path length *L*_{i} is the average of the *n* – 1 minimum path lengths between the index node and all other nodes in the network. The hubs of a network are the nodes with the smallest values of *L*_{i} or largest values of *k*_{i}. All of these quantities can be averaged over nodes to estimate the network means *k*_{net} and *C*_{net} and the characteristic path length *L*_{net}. To diagnose small-world properties, the characteristic path length and clustering coefficient were compared with the same metrics estimated in random networks configured with the same number of nodes, mean degree *k*_{ran}, and degree distribution as the network of interest, under the constraint that *k*_{ran} > log(*n*). Typically, in a small-world network, we expect the ratio γ = *C*_{net}/*C*_{ran} > 1 and the ratio λ = *L*_{net}/*L*_{ran} ∼ 1 (Watts and Strogatz, 1998; Montoya and Solé, 2002). A scalar summary of small-worldness is therefore the ratio σ = γ/λ, which is typically >1 (Humphries et al., 2005).

To apply these tools to analysis of scale-dependent brain functional networks, we first had to draw an undirected graph such that the regions were only connected by an edge in the graph if the wavelet correlation between regions *i* and *j* exceeded a threshold. We defined the threshold primarily in terms of the probability of the observed correlation P(*r _{i,j}* >

*R*) under the null hypothesis that

*r*

_{i,j}is less than an arbitrary value

*R*. Because of the multiple, nonindependent tests entailed by thresholding each of 4005 inter-regional correlations, we controlled the false discovery rate (FDR) at the 5% level. So a functional connection between regions was not regarded as significant unless P(

*r*>

_{i,j}*R*) was less than α

_{FDR}, the

*p*value for an individual test that controlled the FDR at 5% over multiple dependent tests (Benjamini and Yekutieli, 2001).

The group mean wavelet correlation matrices are shown for all six scales of the MODWT (Fig. 1), and the effects of thresholding are illustrated for the scale 4 matrix using three different values of the correlation threshold *R*. Because the value of *R* is increased, fewer connections survive thresholding, and the resulting graphs become correspondingly sparser. If we apply a range of thresholds to the correlation matrices at all six scales, we can see that the mean degree of the graphs becomes monotonically smaller as the threshold becomes higher (Fig. 2). Increasing the threshold corresponds to eliminating the weaker connections that are more likely to be noisy. In contrast, if the threshold is too high, the mean degree *k*_{net} of the resulting graph will be less than the log of the number of nodes [i.e., *k*_{net} < log(*n*) = 4.5 in these data] and, under these circumstances, small-world properties are not estimable (Watts and Strogatz, 1998). Here, we therefore report primarily the results of thresholding wavelet correlation matrices with the maximum value of *R* for which small-world properties are estimable [i.e., the value of *R* at which *k*_{net} = log(*n*)].

*Anatomical and topological visualization of brain functional networks.* To visualize the thresholded correlation matrices in anatomical space, we simply located each regional node according to the coordinates of its centroid in the *y–z* plane of Talairach space and drew a line between connected nodes (Figs. 1, 3). The three-dimensional anatomical distance between connected nodes was estimated by the Euclidean distance between regional centroids *Di,j* ∼ (*x*_{i} – *x*_{j})^{2} + (*y*_{i} – *y*_{j})^{2} +(*z*_{i} – *z*_{j})^{2}, and long-distance connections, defined as those edges connecting regions separated by a distance >7.5 cm, could then be distinguished graphically by line color (Fig. 3). The Euclidean distance here serves merely as an approximation to the true physical distance (axonal length) of connections between regions.

To clarify visualization of the regional hubs, we also mapped the thresholded correlation matrices topologically using NEATO implemented in Graphviz (Koutsofios and North, 1991) (Fig. 4). This software instantiates an algorithm proposed by Kamada and Kawai (1989): for any pair of nodes, the graph is drawn such that the graphical distance between nodes approximates the path length between them; thus, the hubs of the network, which have a short mean path length to other regions, will tend to be clustered in the center of the plot, whereas nodes with a relatively long path length to other regions will tend to be located peripherally. This method can be regarded as a variant of multidimensional scaling for binary (thresholded) matrices rather than continuous (correlation or distance) matrices.

*Resilience of brain functional networks to random error and targeted attack.* To assess the resilience of the brain network to random error and targeted attack, we adopted the approach described by Albert and Barabási (2002). Random error of the network was simulated by selecting one regional node at random and removing it (and all its connections) from the graph before recalculating the size of the largest connected cluster in the network and its mean path length *L*_{net}. We then repeated this process, incrementally eliminating additional nodes from the network at random, until the size of the largest cluster was 1. The same process was applied to simulate a targeted attack, but in this case, the first node to be eliminated was the hub with the largest degree, and nodes were subsequently eliminated in rank order of decreasing degree. The curves describing change in the largest cluster size and mean path length as a function of random error and targeted attack on the brain network were compared with equivalent curves calculated by random and targeted elimination of nodes in an Erdös-Rényi random network and a scale-free network, both of which had the same number of nodes and mean degree as the brain network.

Additionally, to explore the global impact of a single isolated “lesion,” we systematically eliminated each regional node (and all its connections) in turn and re-estimated the mean path length of the network in its absence. The difference in path length Δ*L*_{net} before and after isolated elimination of a single node provides a measure of the centrality of that region to global network topology.

## Results

### Scale dependency of functional connectivity

Brain functional connectivity was most salient in the frequency interval 0.03–0.06 Hz represented by the scale 4 wavelet correlation matrix. The mean wavelet correlation coefficient was greatest at scale 4 (Table 1), and the degree of the thresholded graphs was consistently greatest at scale 4 over a wide range of correlation thresholds (0.2 < *R* < 0.44) (Fig. 2). The maximum correlation threshold for scale 4 [i.e., the value of *R* (0.44) corresponding to mean degree = log(*n*)] was greater than for any other scale. The mean path length of the network obtained by thresholding the scale 4 correlation matrix was *L*_{net} = 2.49, with clustering coefficient *C*_{net} = 0.525. Compared with random graphs matched for number of nodes, mean degree, and degree distribution, the scale 4 brain network had an almost identical path length (λ = 1.09) but was more locally clustered (γ = 2.37), resulting in a small-world scalar σ = 2.18. Small-worldness decreased monotonically with *R*, because lower thresholds make the network topology less clearly distinguishable from matched random graphs, but small-world properties indicated by σ > 1 were evident over a range of values of *R* at all scales. In view of the greater strength of functional connectivity at scale 4, here we focus detailed attention primarily on the scale 4 brain network corresponding to 0.03–0.06 Hz. However, we note that the brain network derived from the scale 5 correlation matrix also had σ ∼ 2 at high values of *R*, implying that small-world properties were salient also in the frequency interval 0.01–0.03 Hz (Table 1).

### Anatomical and topological maps of low-frequency network

An anatomical map of this low-frequency brain network (Fig. 3) had sparse connectivity overall [comprising 405 edges or ∼10% of the total number of possible edges (4005)] but showed neighborhoods of dense local connectivity and clustering, especially in the occipital and parietal cortex, as well as a number of long-distance connections (with Euclidean distance >7.5 cm), especially between regions of the frontal cortex [middle frontal gyrus (MFG), dorsal superior frontal gyrus (SFGdor), inferior frontal gyrus (IFG), supplementary motor area (SMA), superior orbitofrontal cortex] and heteromodal regions of the lateral temporal and parietal association cortex [inferior temporal gyrus (ITG), middle temporal gyrus (MTG), superior temporal gyrus (STG), inferior parietal lobule (IPL), precuneus (PCUN), angular gyrus (ANG), supramarginal gyrus (SMG)]. The network also had marked bilateral symmetry: many regions (43 of 90) were connected to their contralateral homologs, and the right and left intrahemispheric connections were similar. Most brain regions (79 of 90) were connected to a single giant cluster: isolated regions were the left and right olfactory cortex, anterior cingulate cortex, pallidum, putamen, and amygdala (all of which had symmetrical connectivity between bilateral homologs).

A topological map (Fig. 4) clarified the regional identity of the hubs by clustering nodes with a large degree at the center of the plot so that the small geometric distances between the plotted symbols approximated the short path lengths between these highly connected cortical regions. This representation also highlights the relative paucity of connections to limbic/paralimbic regions and the concentration of long-distance connections between hub regions of the association cortex.

### Degree distribution and hubs

The degree distribution was somewhat heavy-tailed or heterogeneous (Fig. 5). Goodness-of-fit was compared using Akaike's information criterion (AIC) for three possible forms of the degree distribution P(*k*): a power law, P(*k*) ∼ *k*^{–α}; an exponential, P(*k*) ∼ *e*^{–α}* ^{k}*; and an exponentially truncated power law, P(

*k*) =

*k*

^{α–1}

*e*. Of these, the exponentially truncated power law was the best-fitting model for the degree distribution (AIC = 545 compared with 578 and 623 for the exponential and power law, respectively) with estimated exponent α = 1.80 and cutoff degree

^{k/kc}*k*= 5.

_{c}Twenty hubs of this network, defined in terms of a regionally characteristic path length averaged over both hemispheres that was less than the network mean path length, *L _{net}* = 2.49, included 15 regions of the heteromodal or unimodal (especially visual) association cortex, 3 regions of the primary sensory or motor cortex, and two regions of the paralimbic cortex (for details, see Table 2). The less well connected nodes, in contrast, included several regions of the limbic or paralimbic cortex such as the hippocampus, parahippocampal gyrus, insula, most temporal pole, and orbitofrontal regions.

Clustering was negatively correlated with connection distance over all neocortical nodes in the network (*r* =–0.58; df = 25; *p* = 0.0046) (Fig. 6). Regions of the parietal, temporal, and frontal heteromodal association cortex (IPL, PCUN, superior parietal gyrus, ANG, SMG, STG, MTG, ITG, MFG, SFGdor, triangular IFG) tended to have both low clustering and long connection distances, indicating that they were strongly connected to remote regions that were not otherwise connected to each other. Regions of the primary and unimodal association cortex [calcarine cortex (CAL), inferior occipital gyrus, middle occipital gyrus (MOG), fusiform gyrus (FFG), lingual gyrus (LING), Heschl's gyrus, SMA] tended to show a contrasting pattern of high clustering and short connection distances, indicating that they were strongly connected to a clique of richly interconnected, neighboring regions. Limbic nodes did not show an association between clustering and distance.

### Resilience of brain functional network

The brain network was approximately as resilient to random error as the random and scale-free networks. However, it was considerably more resilient to targeted attack than the scale-free network (Fig. 7). The size of the largest connected cluster in the scale-free network was reduced by 50% by targeted attack on the most highly connected 20% of hubs; the brain network did not disintegrate to the same extent until ∼40% of the most connected nodes had been attacked. The global mean path length of the brain network increased by a factor of 2 as the top 40% of nodes were eliminated by targeted attack (Fig. 7).

The percentage of change in global mean path length was estimated after individually attacking each regional node. As shown in Table 2, the most damaging attacks were generally on hub regions such as the MFG and SFGdor; PCUN; STG, MTG, and ITG; MOG, FFG, and LING; primary visual, somatosensory, and motor cortex (CAL, postcentral gyrus, precentral gyrus); and thalamus. Elimination of these nodes was associated with a few percentage points increase in global mean path length estimated in the remaining network.

## Discussion

The picture that emerges overall is that of a sparse, resilient, low-frequency, small-world human brain functional network with a heterogeneous degree distribution reflecting a highly connected neocortico-cortical “core” and a less well connected paralimbiclimbic “periphery.”

### Anatomy of human and primate brain small-world networks

At a global level of analysis, the small-world properties of the human brain functional network demonstrated here are quantitatively similar to those reported previously by anatomical studies of cat and macaque cortex using equally carefully matched random networks to assess path length and clustering (Sporns and Zwi, 2004). For example, the anatomical network derived from tract-tracing studies of the entire macaque cortex had λ = 1.16, γ = 3.07, and σ = 2.65, whereas the human brain functional network reported here had λ = 1.09, γ = 2.37, and σ = 2.18.

The anatomical identities of the regional hubs of this markedly symmetric network are arguably not too surprising: large-scale neocortico-thalamic circuits are key physiological substrates for human cognition and consciousness (Tononi et al., 1998; Mesulam, 2000), and the association cortex, by definition, receives rich and convergent inputs from multiple other cortical regions (Mesulam, 2000). More specifically, the only previous study to have measured characteristic path length and clustering for individual cortical regions [based on anatomical connectivity in cat and monkey cortex (Sporns and Zwi, 2004)] also found that many hubs were regions with multimodal or integrative functions, low path lengths, and low clustering coefficients: “In the entire connection system of macaque cortex, areas with markedly low path length and clustering include A7b, LIP, TPT, FEF, and A46, all of which mediate interactions between various sensory and motor systems” (Sporns and Zwi, 2004). The relatively marginal or isolated topological location of the limbic and paralimbic cortex is also not unprecedented anatomically. In macaque cortex, afferent and efferent connectivity was approximately equally balanced, but the anterior cingulate cortex (area 24) and amygdala showed remarkably few afferent connections (Kötter and Stephan, 2003). One might ask whether low functional connectivity in limbic regions demonstrated by fMRI was attributable to susceptibility artifact, because regions such as the hippocampus and orbitofrontal cortex are in close juxtaposition to CSF, air, or bone. However, we found high connectivity in nonlimbic regions also prone to susceptibility effects, such as the inferior temporal cortex, and low connectivity in subcortical regions, such as the pallidum, where susceptibility effects will be small, so we consider this possible interpretation to be unlikely.

In short, we propose that the human brain functional network described here is likely to reflect an underlying anatomical architecture with small-world properties and hubs of long-distance connectivity in regions of the heteromodal association cortex. This prediction could perhaps be tested by future tractography studies of human brain connectivity using diffusion tensor imaging.

### Why the human brain network is not scale-free

If the topology of the human brain anatomical network, which we assume strongly conditions the observed topology of the functional network, had developed or evolved by preferential attachment, or a “rich-get-richer” rule, whereby nodes joining the network were most likely to connect to the nodes that already had the largest number of connections (Albert and Barabási, 2002), one would not expect relatively late-developing regions such as the dorsolateral prefrontal cortex (MFG, SFGdor) to be among the hubs of the functional network. The observation that it is regions of the association cortex, rather than the more primitive limbic or paralimbic regions, that dominate the highly connected core of the entire brain network is one piece of evidence against a scale-free model for this system. Two other pieces of evidence against a scale-free model are (1) that the brain network was more resilient to targeted attack than a comparable scale-free network and (2) that its degree distribution did not follow a power law.

The resilience of the brain network was assessed in response to both random error and targeted attack. Compared with a scale-free network, it was about equally resilient to random error but considerably more resilient to targeted attack. Whereas a scale-free network tends to disintegrate rapidly when the hubs are selectively attacked (Albert et al., 2000), up to 40% of the most connected nodes in the brain network could be eliminated before precipitating a 50% reduction in size (and twofold increase in path length) of the largest connected cluster. It seems that the small-world architecture of the brain may confer distinctive benefits in terms of robustness to both random elimination of nodes and selective attack on hubs, and, of course, one can speculate that this robustness might have considerable fitness value in mitigating the loss of network functionality in the face of developmental aberration or disease.

A log–log plot of the degree distribution of the brain network clearly was not linear (Fig. 5), as it would be for a scale-free network with a power law distribution. Indeed the best-fitting distribution was an exponentially truncated power law that defined a scaling regimen, followed by an exponential decay in probability of hubs with a degree greater than a cutoff value of ∼5. Comparable exponential or truncated power law distributions have been shown previously for the neural network of *C. elegans* as well as the Hollywood actor network, the western United States electrical power supply grid, and the United States airport network (Strogatz, 2001; Amaral et al., 2000). These diverse networks are all physically constrained in one or more ways that are likely to make the emergence of very highly connected hubs less probable than a power law would predict: Hollywood actors grow old and eventually stop making movies; the wiring costs of adding a connection between power generators or neurons separated by long distances may be prohibitive; the dynamical consequence of scheduling another flight into an already congested hub airport may be to disrupt services throughout the network of the airline. It seems likely a priori (Laughlin and Sejnowski, 2003; Sporns et al., 2004; Sporns and Zwi, 2004; Kaiser and Hilgetag, 2004a,b) that similar physical constraints (aging or connection costs) might apply to the formation of whole-brain functional networks and thus account for the empirical form of this degree distribution.

### Comparable previous fMRI studies

Small-world properties have been reported previously for networks derived from time-domain analysis of inter-regional correlations in fMRI (Salvador et al., 2005b; Eguíluz et al., 2005), and there have been previous fMRI studies of frequency dependency of functional coherence between regions in the Fourier domain (Sun et al., 2004; Salvador et al., 2005a; Wink et al., 2005). However, this is the first study both to use wavelets for frequency decomposition of inter-regional correlations or functional connectivity and to describe the anatomy and small-world topology of the resulting low-frequency networks. Wavelet correlation analysis parsimoniously summarizes the frequency dependency of brain functional connectivity in terms of a small number of scale-specific correlations (compared with a larger number of frequency-specific coherences in the Fourier domain) and the statistical properties of covariance estimators based on the MODWT have been well described (see Appendix for detail). Eguíluz et al. (2005) reported a small-world and scale-free topology in networks derived from experimentally activated fMRI time series at voxel resolution. The difference between the power law degree distribution of their network, based on activated voxel data, and the truncated power law degree distribution of the network reported here, based on resting-state fMRI data at regional resolution, suggests that small-world properties of brain functional networks could be conditioned by anatomical resolution of analysis and/or experimental stimulation of the subjects.

### Neural and cognitive correlates of correlated low-frequency oscillations in fMRI

Our finding that small-world networks exist at low frequencies (0.01–0.11 Hz), and most saliently in the 0.03–0.06 Hz interval, is consistent with many previous studies of resting-state fMRI data, which have typically reported strongest resting-state inter-regional correlations at low frequencies <0.1 Hz (Biswal et al., 1995; Lowe et al., 1998; Cordes et al., 2000; Salvador et al., 2005a). There are also comparable observations of low-frequency oscillations in optical imaging data (Mayhew et al., 1996). Both age and muscarinic receptor blockade increased the spectral power of low frequencies in hippocampal fMRI time series recorded at rest, and these changes were associated with increased frontohippocampal coherence (Wink et al., 2005). Experimental studies manipulating the sampling rate for fMRI (e.g., using short TR sequences) have shown that low-frequency power or coherence is unlikely to represent (possibly aliased) signals related to the cardiac or respiratory cycles (Cordes et al., 2000; Kiviniemi et al., 2005). Moreover, long period or very low-frequency coherent oscillations have been described recently in both monkey and human electrophysiological data (Linkenkaer-Hansen et al., 2001; Leopold et al., 2003). For example, slow (< 0.1 Hz) and coherent fluctuations in band-limited power of multielectrode local field potential recordings from monkey visual cortex were clearly of neuronal origin and might be related to low-frequency fluctuations in fMRI time series (Leopold et al., 2003). Together, these data strongly suggest that low-frequency oscillations and correlations in resting-state fMRI data are neurobiologically interesting and may reflect endogenously coordinated dynamics in large-scale neuronal populations.

## Appendix

We used the MODWT (Percival and Walden, 2000), instead of the orthogonal DWT, to estimate wavelet correlations. The MODWT is a redundant transform that is translation invariant and straightforward to compute using the pyramid algorithm. MODWT-based estimators of wavelet correlation are superior to those derived from the orthogonal DWT (Whitcher et al., 2000).

### MODWT

The MODWT is defined formally as follows. Let **X** be a time series of length *N*. Let {*h*_{j,l};*l* = 0,...,*L*_{j} – 1} and {*g*_{j,l};*l* = 0,...,*L*_{j} – 1} be, respectively, a *j*th level wavelet filter and scaling filter; here, *L*_{j} = (2* ^{j}* – 1)(

*L*– 1) + 1, and

*L*denotes the width of the initial filter. The corresponding

*j*th level MODWT wavelet and scaling filters are defined, respectively, by and and have the same width

*L*

_{j}. Then the

*j*th level MODWT wavelet and scaling coefficients are

*N*dimensional vectors denoted, respectively,

**W**

_{j}and

**V**

_{j}and defined, for

*t*= 1,...,

*N*– 1, by the following: (A1) and (A2) The MODWT yields an energy decomposition: (A3)

### MODWT estimator of wavelet correlation

Let us denote by {*X*_{t}}_{t = 0,...,N – 1} and {*Y*_{t}}_{t = 0,...,N – 1}, realizations of the processes {*X*_{t}}_{t} and {*Y*_{t}}_{t} that are Gaussian processes with stationary increments. Then for *N* ≥ *L*_{j}, the scale-dependent covariance between {*X*_{t}} and {*Y*_{t}} is estimated by the following: (A4) which is an unbiased estimator of (A5) where are the scale λ_{j} = 2^{j}^{–1} MODWT coefficients for {*X*_{t}} and {*Y*_{t}}, respectively, and *N*_{j} = *N* – *L*_{j} + 1 (Whitcher et al., 2000).

An estimator of the wavelet correlation is then defined (Whitcher et al., 2000) as follows: (A6) where is the wavelet variance for the time series **X**. Assuming that *L* > 2*d*, and is a bivariate Gaussian weakly stationary process with square integrable autospectra, then the MODWT estimator of the wavelet correlation ρ_{XY}(λ_{j}) is asymptotically normally distributed with mean ρ_{XY}(λ_{j}) and a large sample variance (Gençay et al., 2001).

## Footnotes

This neuroinformatics research was supported by a Human Brain Project grant from the National Institute of Bioengineering and Biomedical Imaging and the National Institute of Mental Health and was conducted in the Medical Research Council (MRC)/Wellcome Trust Behavioral and Clinical Neurosciences Institute (Cambridge, UK). The Wolfson Brain Imaging Centre was supported by an MRC Cooperative Group grant. This work was presented in part at the Brain Connectivity 2005 Workshop at Florida Atlantic University (Boca Raton, FL), April 15–16, 2005.

Correspondence should be addressed to Prof. Ed Bullmore, Brain Mapping Unit, University of Cambridge, Department of Psychiatry, Addenbrooke's Hospital, Cambridge CB2 2QQ, UK. E-mail: etb23{at}cam.ac.uk.

DOI:10.1523/JNEUROSCI.3874-05.2006

Copyright © 2006 Society for Neuroscience 0270-6474/06/260063-10$15.00/0