Eighteen right-handed healthy volunteers (15 male, 3 female) were recruited from the GlaxoSmithKline (GSK) Clinical Unit Cambridge, a clinical research facility in Addenbrooke's Hospital, Cambridge, UK. All volunteers (mean age 32.7 years ± 6.9 SD) had a satisfactory medical examination prior to study enrolment and were screened for any other current Axis I psychiatric disorder using the Structured Clinical Interview for the DSM-IV-TR Axis I Disorders (SCID). Participants were also screened for normal radiological appearance of structural MRI scans by a consultant neuroradiologist, and female participants were screened for pregnancy. Urine samples were used to confirm abstinence from illicit drugs and breath was analysed to ensure that no participant was under the influence of acute alcohol intoxication. All volunteers provided written informed consent and received monetary compensation for participation. The study was reviewed and approved by the Cambridge Local Research Ethics Committee (REC06/Q0108/130; PI: TW Robbins).

Whole-brain echoplanar imaging (EPI) data depicting BOLD contrast were acquired at the Wolfson Brain Imaging Centre, University of Cambridge, UK, using a Siemens Magnetom Tim Trio whole body scanner operating at 3 T with a birdcage head transmit/receive coil. Gradient-echo, EPI data were acquired for the whole brain with the following parameters: repetition time = 2000 ms; echo time = 30 ms, flip angle = 78°, slice thickness = 3 mm plus 0.75 mm interslice gap, 32 slices parallel to the inter-commissural (AC-PC) line, image matrix size = 64 × 64, within-plane voxel dimensions = 3.0 mm × 3.0 mm.

Participants were asked to lie quietly in the scanner with eyes closed during the acquisition of 300 images. The first four EPI images were discarded to account for T1 equilibration effects, resulting in a series of 296 images, of which the first 256 images were used to estimate wavelet correlations.

The images were corrected for motion and registered to the standard stereotactic space of the Montreal Neurological Institute EPI template image using an affine transform (Suckling et al., 2006). Time series were then extracted using a whole brain, high resolution, regional parcellation of the images, implemented in the following manner; see Figure 1A. First, a binarized version of a commonly used template image (Tzourio-Mazoyer et al., 2002) was used as a broad grey matter mask. Second, each 8 mm³ voxel in this mask was downsampled by a factor of 4 such that each equally sized region in the parcellation comprised 4 × 4 × 4 voxels of the original image. This initial parcellation included some regions of the image which were not largely representative of grey matter: these were excluded from further analysis by applying the criteria that each region must be at least 50% overlapping with the grey matter mask and must contain at least 80% voxels having BOLD signal (defined operationally as mean signal intensity >50). To be included in the definitive parcellation scheme (which comprised 1808 regional nodes), a region had to satisfy these two inclusion criteria for every individual dataset in the sample.

The mean time series of each region was extracted and wavelet-filtered using Brainwaver R package¹ (Achard et al., 2006; Achard and Bullmore, 2007). The wavelet correlation coefficient was estimated for each of four wavelet scales between each pair of nodes, resulting in a {1808 × 1808} association matrix, or frequency-dependent functional connectivity matrix, for each wavelet scale in the overall frequency range 0.25–0.015 Hz. In what follows, we will focus on results at wavelet scale 3, subtending a frequency interval of 0.06–0.03 Hz.

This choice of frequency interval was guided by the fact that prior work on resting-state fMRI functional connectivity has found that the greatest power in connectivity occurs in frequency bands lower than 0.1 Hz (Cordes et al., 2001). However, analysing very low frequency scales in limited time series such as those acquired with fMRI can reduce precision in estimating inter-regional wavelet correlations (Achard et al., 2006). So scale 3 was chosen for the focus of this study as representing a reasonable compromise between retaining sufficient estimation precision while measuring low frequency network properties.

Each association matrix was thresholded to create an adjacency matrix A, the a_i,jth element of which is either 1, if the absolute value of the wavelet correlation between nodes i and j, w_i,j, exceeds a threshold value τ; or 0, if it does not. We have chosen here to take a high threshold, leading to very sparse networks comprising 8000 edges, i.e. with a connection density of 0.5% of all possible edges in a network of this size. Modularity of neuroimaging networks is typically greater (Meunier et al., 2009), and computational costs are lower, when the networks are more sparsely thresholded. Up to 10% of nodes were disconnected from the rest of the network at this threshold.

In recent years, many methods have been proposed to discover the modular organization of complex networks. A key step was taken when Girvan and Newman popularized graph-partitioning problems by introducing the concept of modularity. Modularity is by far the most widespread quantity for measuring the quality of a partition 𝒫 of a network. In its original definition, an unweighted and undirected network that has been partitioned into communities has modularity (Newman and Girvan, 2004): (1)Q=12m∑C∈P ∑i,j∈C [Aij−kikj2m]  where A is the adjacency matrix of the network; m is the total number of edges; and ki=∑jAij is the degree of node i. The indices i and j run over the N nodes of the graph. The index C runs over the modules of the partition 𝒫. Modularity counts the number of edges between all pairs of nodes belonging to the same community or module, and compares it to the expected number of such edges for an equivalent random graph. Modularity therefore evaluates how well a given partition concentrates the edges within the modules.

A popular method for discovering the modules of a network consists in optimizing modularity, namely in finding the partition having the largest value of 𝒬. However, it is typically impossible computationally to sample modularity exhaustively by enumerating all the possible partitions of a network into communities. Several heuristic algorithms have therefore been proposed to provide good approximations, and so to allow for the analysis of large networks in reasonable times. The computational expediency of the algorithm has become a crucial factor due to the increasing size of the networks to be analysed.

More recently, methods to study hierarchical modularity, also called nested modularity, have been introduced (Sales-Pardo et al., 2007; Arenas et al., 2008; Rosvall and Bergstrom, 2008). In this case, each module obtained at the partition of the highest level can further be decomposed into sub-modules, which in turn can be decomposed into sub-submodules, and so on. Here, we will use a multi-level method which was introduced very recently in order to optimize modularity (Blondel et al., 2008); see Figure 1B. The primary advantages of this method are that it unfolds a complete hierarchical community structure for the network and outperforms previous methods with respect to computation time. This so-called “Louvain method” takes advantage of the hierarchical organization of complex networks in order to facilitate the optimization. The algorithm starts by assigning a different module to each node of the network. The initial partition of the network is therefore made of N communities. It then consists of two phases that are repeated iteratively. The first phase consists in a greedy optimization where nodes are selected sequentially in an order that has been randomly assigned. When a node is selected, it may leave its community and adopt a community which is in its direct neighbourhood, but only if this change of community leads to an increase of modularity (GO on Figure 1B). The second phase builds a new network whose meta-nodes are the communities found in the first phase (MN on Figure 1B). Let us denote by “pass” a combination of these two phases. These passes are repeated iteratively until a maximum of modularity is attained and an optimal partition of the network into communities is found. By construction, the meta-nodes, or intermediate communities, are made of more nodes at subsequent passes. The optimization is therefore done in a multi-scale way: among adjacent nodes at the first pass, among adjacent meta-nodes at the second pass, etc. The output of the algorithm is a set of partitions, one for each pass. The optimal partition is the one found at the last pass. It has been shown on several examples that modularity estimated by this method is very close to the optimal value obtained from slower methods (Blondel et al., 2008). Intermediate partitions can also be shown to be meaningful and to correspond to communities at intermediate resolutions (see Section “Discussion”). In the following, we will call “lowest non-trivial level” the partition found after the first pass.

Once a maximally modular partition of the network has been identified, it is possible to assign topological roles to each node based on its density of intra- and inter-modular connections (Guimerà and Amaral, 2005a,b; Guimerà et al., 2005; Sales-Pardo et al., 2007).

Intra-modular connectivity is measured by the normalized within-module degree, (2)zi=κni−κ¯nσκn where κni is the number of edges connecting the ith node to other nodes in the nth module, κ¯n is the average of κni over all nodes in the module n, and σκn is the standard deviation of the intra-modular degrees in the nth module. Thus z_i will be large for a node that has a large number of intra-modular connections.

Inter-modular connectivity is measured by the participation coefficient, (3)Pi=1−∑n=1N(κniki)2 where κni is the number of edges linking the ith node to other nodes in the nth module, and k_i is the total degree of the ith node. Thus P_i will be close to 1 if the ith node is extensively linked to all other modules in the community and 0 if it is linked exclusively to other nodes in its own module.

The two-dimensional space defined by these parameters, the {P, z} plane, can be partitioned to assign categorical roles to the nodes of the network. Contrarily to our previous study (Meunier et al., 2009), where we used a simplified definition of node roles, the higher number of nodes examined in the current study allowed us to adopt the original definitions of node roles as described for large metabolic (Guimerà et al., 2005) and transportation networks (Guimerà and Amaral, 2005b):

These different categories allowed us to classify the nodes according to their topological functions in the network. For example, a provincial hub is a hub with greater intra- vs inter-modular connectivity, thus having a pivotal role in the function realized by its module, whereas a connector hub will play a central role in transferring information from its module to the rest of the network.

The results of modular decomposition were visualized in anatomical space using Caret software for cortical surface mapping², and in topological space using Pajek software³.

To compare the different modularity partitions obtained at different hierarchical levels in the same subject, or at the same hierarchical level in different subjects, we used the normalized mutual information, as defined in Kuncheva and Hadjitodorov (2004). For two given partitions A and B, with a number of communities denoted C_A and C_B: (4)I(A,B)=−2∑i=1CA∑j=1CBNijlog(NijNNi.N.j) ∑i=1CANi.log(Ni.N) +∑j=1CBN.jlog(N.jN)  where N_ij is the number of nodes in common between modules i and j, the sum over row i of matrix N_ij is denoted N_i, and the sum over column j is denoted N_.j. If the two partitions are identical then I(A,B) takes its maximum value of 1. If the two partitions are totally independent, I(A,B) = 0.

The initial application of this quantity was to evaluate different modularity partition algorithms (Danon et al., 2005). The similarity index was used to compute how closely the partitions obtained from different algorithms matched the “target” partition of a given test network, i.e. a network whose modular structure was known a priori. Here the application was different, since we wanted to compare partitions obtained for different subjects in a group. Since the equation is symmetric in A and B, it is however possible to use the index without a target partition.

The networks constructed for each individual had the same number of nodes N, so the partitions of each subject have the same number of nodes. However, due to the high threshold applied to construct the adjacency matrix, the number of disconnected nodes in the networks can be different for each subject. One solution is to consider each disconnected node as a single module. In this case, each node (disconnected or not) of the network will be in the set of modules of each subject. However, it introduces artificially high values in the similarity values, especially if the networks of two subjects have similar sets of disconnected edges. So we have chosen to remove the disconnected nodes from the partitions and study only the partitions obtained on the giant component of each network, but keeping the value of N in the equation as the total number of nodes. This leads to a value of similarity slightly lower than if the disconnected nodes were included in the partitions, but is more representative of the relevant set of connected modules.

A first consideration concerns the choice of the Louvain method (LM) in order to uncover nested modules in the brain networks. LM was first proposed in order to uncover optimal partitions of a graph by maximising modularity. This is a greedy method which is known to be very fast and very precise (Blondel et al., 2008), albeit less precise than much slower methods such as simulated annealing (SA). It is interesting to note, however, that this lack of precision may be an advantage, in practice, as it may avoid some of the pitfalls of modularity analysis such as its resolution limit (Fortunato and Barthélemy, 2007). For instance, it has been recently shown that LM performs much better than SA when applied to benchmark networks with unbalanced modules comprising different numbers of nodes (Lancichinetti and Fortunato, 2009). We therefore believe that there is good evidence that the top level partitions uncovered by LM are valid. The validity of the intermediate hierarchical levels identified by the algorithm is, however, more arguable, as it has not been studied in detail yet. In order to show the validity of these intermediate levels, we need to verify that the method uncovers all the significant partitions present in the network and only those.

To do so, we have tested LM on a benchmark network with known hierarchical structure (Sales-Pardo et al., 2007); Figure 5A). This benchmark network is made of 640 nodes with three levels of organization: small modules comprising 10 nodes, medium-size modules comprising 40 nodes and large modules comprising 160 nodes. The cohesiveness of the hierarchy between levels is tuned by a single parameter σ, i.e. the larger the value of σ, the more difficult it is to find the sub-modules. When applied on this benchmark network, the algorithm finds with an excellent precision the first two levels (16 modules and 64 modules), but does not uncover the partition into 4 modules. This result is to be expected because this partition into four modules is sub-optimal in terms of modularity and can therefore not be uncovered by an aggregative method. This shows that the method can at best uncover the partition optimising modularity and finer partitions. In order to uncover coarser partitions, one needs to decrease the resolution of the method, which can be done by following Reichardt and Bornholdt (2006), or Sales-Pardo et al. (2007), for instance.

On the same benchmark network, the algorithm typically finds two levels (one corresponding to 64 modules and one corresponding to 16 modules) but it may occasionally find three levels (one level corresponding to 64 modules and two levels similar to the partition into 16 modules). When σ = 1.0, for instance, over 100 realizations of the graph, the algorithm finds two levels on 86 realizations, and three levels on 14 realizations. This result is encouraging as it suggests that the algorithm only produces significant partitions. However, it is possible to find situations where it is not the case, e.g. random graphs. It is therefore still necessary to verify the significance of intermediate partitions, as we will discuss below.

A second consideration concerns the comparison of the partition of the original network with randomized data, as the algorithm also gives a hierarchical decomposition for randomly generated networks. To show that the representative brain network under study (subject ID 2) displays a non-random hierarchical modular structure, we have randomized the original data and processed the hierarchical structure of randomized networks, with two kind of randomization. First, by computing 100 randomizations of the time points in the original time-series (in green on Figure 5B) and, second, by randomising the original adjacency matrix 100 times (in blue on Figure 5B). Note that the two kinds of randomization lead to networks with different sizes: in the randomized time-series networks, almost all the nodes are connected, thus leading to networks with 1808 nodes and 8000 edges. Whereas starting from the original adjacency matrix leads to networks of 844 nodes and 8000 edges. The modularity obtained for the lowest and highest partitions of the original network are displayed in Figure 5B. The modularity values are clearly reduced in the randomized networks, relative to the original data, indicating that our results on real brain networks are not trivially reproduced in random networks.

In order to show that the intermediate levels considered in this paper are significant, we have followed the argument that significant partitions should be robust, in the sense that they should only be weakly altered by a modification of the optimization algorithm. As argued by Ronhovde and Nussinov (2009), comparing the optimal partitions found by the algorithm for different orders of the nodes is a way to test their robustness and therefore their validity. We have therefore optimized the modularity of the representative brain network 100 times by choosing the nodes in a different order, and focused on the first non-trivial partition found by the algorithm. The mutual information between pairs of partitions obtained for each different order is then computed. The average mutual information among those pairs is very high (0.89) compared to what is obtained for a comparable random network (0.44), thereby suggesting that partitions obtained at the lowest non-trivial levels are relevant for the network under study.

A third consideration concerns the number m of edges that we have chosen in order to map the correlation matrices onto unweighted graphs. This is a known problem when dealing with fMRI data and building brain networks. If m is too small, i.e. keeping the top most significant links, the network will be so sparsely connected that it will be made of several disconnected clusters. If m is too large, in contrast, the network will be very densely connected, but mainly made of unsignificant links. In these two extremes, the network structure is a bad representation of the correlation matrix. This is still an open problem that requires the right trade-off between these two competing factors. In order to show the robustness of our results, we propose to look at the resilience of the hierarchical modular organization under the tuning of the value of m. Meaningful values of m are identified by intervals over which the structure of the network is preserved. We have applied this scheme to the optimal partitions of the most representative subject (Subject 2), over a wide range of threshold (2000–14000 edges, with a step of 200 edges). Our results show that partitions are very similar (in terms of mutual information) over the range (6000–11000) for both highest level (left on Figure 5C) and non-trivial lowest level (right on Figure 5C), indicating our results are robust to the specific choice of threshold.