Introduction Understanding the development of human brain organization is critical for gaining insight into brain organization and functions in adulthood as well as for investigating disorders such as autism spectrum disorders (ASD) and attention-deficit/hyperactivity disorder (ADHD), where normal developmental processes are disrupted. Neuroimaging studies of development have primarily focused on structural changes from childhood, to adolescence, and into adulthood. These studies have reported age-related changes in (1) overall brain volumes ,, (2) volumes of individual brain areas ,, (3) regional cortical thickness ,, as well as (4) regional and global grey-matter and white-matter densities –. Collectively these studies have suggested that the human brain undergoes vast developmental changes in grey and white matter structure between childhood and adulthood. These changes are thought to reflect synaptic pruning and myelination observed at the neuronal level ,. More recently, diffusion tensor imaging (DTI) studies investigating the development of white-matter pathways have shown increase in anisotropy –, decrease in overall diffusion , and maturation in major white-matter fiber tracts –, with age. In spite of growing evidence from these studies for patterned brain development, the functional organization of the human brain in childhood is not well understood and it is also not clear how the above structural changes translate to differences in functional brain organization between children and adults. Task-free (resting-state) functional connectivity MRI is a useful technique for investigating the functional organization of the human brain. This method detects interregional correlations in spontaneous blood oxygen level-dependent (BOLD) signal fluctuations ,, and has been used to investigate brain networks involved in motor , sensory , attention , salience and cognitive control , and memory , systems. However, only a small number of studies have examined developmental changes in functional brain organization. A few recent studies have examined developmental changes in functional connectivity of brain regions involved in attention and cognitive control  and the default mode network (DMN) , as well as in functional connectivity of anatomical structures such as the anterior cingulate cortex . To our knowledge, the developmental changes in the functional organization of large-scale networks at the whole-brain level have not yet been investigated. Here we use a graph theoretical approach to examine developmental changes in the large-scale functional organization of the human brain. Graph metrics such as the clustering coefficient and the characteristic path length , have been shown to be useful measures of organization of large-scale networks. Briefly, graphs are data structures that have nodes and edges between the nodes . The clustering coefficient is a measure of local network connectivity. A network with a high average clustering coefficient is characterized by densely connected local clusters. The characteristic path length is a measure of how well a network is connected. A network with a low characteristic path length is characterized by short distances between any two nodes. Many biological systems have small-world network properties, characterized by a high clustering coefficient and a low characteristic path length ,. These graph-theoretic metrics have also proven useful in modeling the large-scale functional and structural organization of the human brain –. In a graphical representation of a brain network, a node corresponds to a brain region while an edge corresponds to the functional connectivity between two brain regions. Functional connectivity networks of the human brain derived from electroencephalograms (EEGs), magnetoencephalograms (MEGs), and task-free functional magnetic resonance imaging (fMRI) data have been shown to exhibit small-world characteristics ,,. These studies suggest that small-world metrics are suited to quantify the global topological properties of large-scale organization of the human brain. Recently, in addition to small-world metrics, Bassett and colleagues used graph theoretic metrics such as hierarchy to characterize local topological properties of large-scale organization of the human brain. Using structural brain imaging data and modeling of interregional covariance in cortical thickness, they reported that hierarchical organization in anatomical human brain networks is characterized by the presence of frontal hubs . A recent study of aging by Meunier and colleagues investigated the modular organization of large-scale functional brain networks using Newman's graph-based modularity metric. They reported that while both young and older adults showed modularity of network organization, the topological roles of the specific brain regions as well as the intermodular connectivity was significantly different between the two groups . The use of small-world metrics along with more advanced graph theoretic metrics to characterize local organization of complex networks provides a new approach for investigating large-scale functional organization of the human brain at multiple levels of granularity. We investigated developmental changes in the functional organization of large-scale brain networks at multiple levels by (1) creating whole-brain functional connectivity networks from task-free fMRI data, (2) characterizing the organization of these networks using metrics of global and local brain organization (including small-worldness and hierarchy, as defined in the Methods section), and (3) comparing these metrics of global and local brain organization between healthy children (ages 7–9 y) and young-adults (ages 19–22 y). In older adults (age ±40 y) it is now well established that large-scale brain networks have a small-world architecture that reflects a robust and efficient, nonrandom, functional organization ,,. Whether children and younger adults have a similar functional brain organization is currently not known. This question is important from a developmental perspective because the brain undergoes vast changes in structural connectivity during adolescence . We hypothesized that the global functional organization of brain networks would be characterized by nonrandom, efficient, small-world characteristics in both subject groups, but that young-adults would show higher small-worldness compared to children, on the basis of previous neurobiological studies in humans and animals suggesting that developmental changes improve efficiency of information processing in the brain –. We further predicted that local organization patterns would be significantly different in children, reflecting a process of continuing structural maturation during the period between childhood and young adolescence. To further characterize developmental changes in the global and local functional organization of brain networks, we used the parcellation scheme of Mesulam  to examine functional organization in five key subdivisions: primary sensory, subcortical, limbic, paralimbic, and association areas. Additionally, developmental changes in the connectivity between these subdivisions (hereafter referred to as interregional connectivity in the manuscript) were examined. Lastly, to characterize the underlying developmental processes that produce these changes in the global and local functional organization of large-scale brain networks, we examined changes in functional connectivity as a function of DTI-based wiring distance between distinct brain regions. The formation of brain networks during development is thought to arise from a dual process of integration and segregation ,–. Accordingly, we investigated whether there is in vivo developmental evidence for the emergence of functional segregation and integration in large-scale brain networks. Results Participants Demographic and cognitive profile data for the child and young-adult groups are shown in Table 1. The two groups were well-matched and did not differ in IQ (p = 0.93) or gender (p = 0.75). 10.1371/journal.pbio.1000157.t001 Table 1 Participant characteristics. Measure Children (n = 23) Young-adults (n = 22) Age 7.95a (range: 7–9) 20.40a (range: 19–22) Gender 10 males, 13 females 11 males, 11 females IQ 112 (range: 88–137) 112 (range: 97–137) Years of education 2.52a (range: 2–3) 14.5a (range: 13–16) Age and years of education, but not IQ nor gender, are significantly different in young-adults compared with children. a Denotes significant differences between groups. Analyses of Small-World Metrics at Different Frequency Scales We first examined graph-theoretic metrics obtained for the functional brain networks constructed by thresholding (threshold values ranged from 0.01 to 0.99, with an increment of 0.01) the wavelet correlation matrix at three different frequency scales. Scale 1 encompassed 0.13–0.25 Hz, scale 2 encompassed 0.06–0.12 Hz, and scale 3 encompassed 0.01–0.05 Hz. As shown in Figure 1, for both the children and young-adult groups, the mean degree—the average number of edges incident on a node belonging to the network—was highest at scale 3 for a wide range of correlation thresholds (0.01 1) showed a monotonic increase in small-worldness as the threshold increased and the degree decreased. σ values for higher correlation thresholds are difficult to interpret because at higher threshold values, graphs of functional brain networks have fewer edges (smaller degree) and tend to split into isolated subgraphs. Graph metrics such as clustering coefficient, characteristic path length, and small-world property do not meaningfully characterize network structures that are not composed of a single, large group of interconnected nodes . 10.1371/journal.pbio.1000157.g001 Figure 1 Developmental changes in whole-brain functional connectivity network metrics. Graph metrics: degree, path length (λ), clustering coefficient (γ), small-worldness (σ), for children (Δ) and young-adults (○) at three frequency intervals. (A) For both groups, the mean degree, a measure of network connectivity, is highest at scale 3 for a wide range of correlation thresholds (0.01 0.6, mainly due to the large variance observed at higher threshold values. 10.1371/journal.pbio.1000157.g003 Figure 3 Developmental changes in network metrics for five major functional divisions of the human brain. Graph metrics—degree, path length (λ), efficiency, clustering coefficient (γ), within each of the five divisions: association, limbic, paralimbic, primary, and subcortical—are shown for children (blue) and young-adults (red), as a function of the correlation threshold. In the subcortical division, for threshold values from 0.1 to 0.6, degree and efficiency values were significantly higher and λ values significantly lower in children, compared to young-adults (p 1, σyoung-adults>1) suggesting the presence of subnetworks of densely connected nodes, mostly connected by a short path. Similar findings were observed when clustering coefficient and global efficiency were used as alternative measures of small-worldness. Small-world characterization is well-suited for analyzing functional brain networks at the systems level because these networks are complex and optimally connected to minimize information processing costs ,. Functional connectivity networks of the human brain constructed from EEG as well as MEG data have also been shown to have small-world architecture ,. Salvador et al.  examined connectivity in task-free functional MRI data with the same 90 ROI parcellation scheme used in our study and they reported small-world architecture in this network. This finding was replicated by Achard et al., who also reported that small-world properties were salient in the low frequency interval 0.03–0.06 Hz  in adults (ages 25–35 y), and by Supekar et al. in older adults (ages 37–77 y) . These findings, primarily derived from functional data obtained from middle-age to older adults, suggest that the functional organization of the brain has a small-world architecture, a characteristic that may assist in robust and dynamic information processing. Our finding that large-scale brain networks in children showed small-world properties that were very similar to young-adults, together with the above observations, suggests that key aspects of functional brain organization are conserved throughout the developmental process—from early childhood to young adulthood and into older adulthood. Critically, despite the fact that the brain undergoes vast structural reorganization at the neuronal level in the form of myelination and synaptic pruning throughout development, key global properties of functional organization appear to be conserved. Functional Brain Connectivity Patterns in Children and Young-Adults Are Significantly Different Notwithstanding similarities in global, whole-brain, small-world network properties, functional connectivity patterns in children were significantly different from those in young-adults. SVM-based pattern classification analysis showed that connectivity patterns in children could be distinguished from those in young-adults with an accuracy of over 90%. Accuracy was highest (91%) for connectivity patterns in the low frequency interval (scale 3; 0.01–0.05 Hz). Previous studies have reported that resting-state functional connectivity is most robust at frequencies below 0.1 Hz ,, and that these low frequency fMRI fluctuations are related to interregional coupling of local field potentials in the gamma band ,. Overall, these findings suggest that observed developmental changes in the functional connectivity measured by fMRI resting state signals are likely to reflect underlying differences in coupling of neuronal signals. We discuss below the nature of developmental changes in the context of hierarchical and regional organization of brain connectivity. Development of Functional Hierarchical Organization Our data provide new evidence that large-scale brain networks in children and young-adults differ in their hierarchical organization. Children showed significantly lower (p 232) was assumed to be a symmetric reflection of itself. At each of the three scales, wavelet correlations between signals in the 90 anatomical regions were determined by computing the correlation coefficient between the transformed signals at that scale. For each subject, a 90-node, scale-specific, undirected graph of the functional connectivity network was constructed by thresholding the wavelet correlation matrix computed at that scale. If the wavelet correlation value between two anatomical regions represented by nodes i and j in the network exceeded a threshold, then an edge was drawn between node i and node j. There is currently no formal consensus regarding threshold selection, so we computed networks for threshold values from 0.01 to 0.99 with an increment of 0.01. Once a whole-brain functional connectivity network was constructed from the correlation matrix, we characterized this network using graph theoretic metrics of global and local brain organization including small-worldness and hierarchy. Small-World Analysis of the Whole-Brain Functional Connectivity Small-world properties of a network are described by the clustering coefficient and the characteristic path length of the network. The clustering coefficient and characteristic path length of functional brain networks generated from the task-free fMRI data obtained from 23 children and 22 young-adults were computed. The clustering coefficient of every node was computed as the ratio of the number of connections between its neighbors divided by the maximum possible connections between its neighbors. The clustering coefficient (C) of the network was calculated as the mean of the clustering coefficients of all the nodes in the network. The mean minimum path length of a node was computed as the average of minimum distances from that node to all the remaining nodes in the network. The characteristic path length (L) of the network was the average of the mean minimum path lengths of all the nodes in the network. The clustering coefficient and path length of nodes completely disconnected with the network were set as 0 and Inf respectively, and these nodes were excluded while computing C and L. To evaluate the network for small-world properties, we compared the clustering coefficient and the characteristic path length of the network with corresponding values (C ran, L ran) obtained and averaged across 1,000 random networks with the same number of nodes and degree distribution . The degree of every node (a measure of its connectivity) was calculated by counting the number of edges incident on that node. The mean degree of the network was the average of the degree of all the nodes in the network. Small-world networks are characterized by high normalized clustering coefficient γ (C/C ran)>1 and low normalized characteristic path length λ (L/L ran)≈1 compared to random networks . A cumulative metric σ—the ratio of normalized clustering coefficient (γ) to the characteristic path length (λ), a measure of small-worldness—is thus greater than 1 for small world networks. Analysis of Global Efficiency of Whole-Brain Functional Connectivity Small-world networks are characterized by high clustering coefficient and low characteristic path length. These small-world metrics, particularly the path length, are not meaningful when the graph contains disconnected nodes. To address this issue, we ensured that only small-world metrics computed on connected graphs were considered in our analysis. Specifically, the algorithm used to choose the correlation threshold (R) guaranteed that disconnected graphs were excluded from the analysis. Also, in the node-wise clustering coefficient comparison analysis, we only considered thresholds from 0.1 to 0.6. We chose these thresholds because beyond 0.6 the network gets divided into disconnected subset of nodes. To determine if our characteristic path length findings were robust and reliable, we computed the efficiency of functional brain networks. It has been previously reported that efficiency as a graph metric (1) is not susceptible to disconnected nodes, (2) is applicable to unweighted as well as weighted graphs, and (3) is a more meaningful measure of parallel information processing than path length . Efficiency of a graph (E global−net)  is the inverse of the harmonic mean of the minimum path length between each pair of nodes, L ij, and was computed as, (1) To evaluate the network for its global efficiency of parallel information processing, we compared the global efficiency of the network (E global−net) with corresponding values (E global−ran) obtained and averaged across 1,000 random networks with the same number of nodes and degree distribution. A network with small-world properties is characterized by a global efficiency value that is lower than the random network: E global (E global−net/E global−ran)<1. Analysis of Hierarchical Organization of Whole-Brain Functional Connectivity We evaluated the hierarchical nature of the large-scale whole-brain functional connectivity network by the β parameter . β measures the extent of the power-law relationship between the clustering coefficient (C) and the degree (k): C≈k −β. The clustering coefficient (C i) and the degree (k i) of every node was computed; β of the network was calculated by fitting a linear regression line to the plot of log(C) versus log(k). Analysis of Regional Differences in Network Organization and Connectivity The human brain can be divided into five major divisions—association, limbic, paralimbic, primary, and subcortical—each of them having a distinct function . We assessed the network organization of these cortical divisions and how it differs in development by examining the regional profile of metrics (degree, λ, E global, and γ) at the divisional level. The 90 anatomical regions of our network were grouped into these five cortical divisions. The association division consists of 44 regions, the limbic division consists of 12 regions, the paralimbic division consists of 24 regions, the primary division consists of eight regions, and the subcortical division consists of eight regions (see Table S1 for region-wise division assignment). The graph metrics (degree, λ, E global, and γ) of 90 regions were aggregated into five divisions and the aggregated metrics in the two subject groups were compared using growth curve modeling, with an intercept, linear, and quadratic terms. In the aggregation step, the graph metric value at a correlation threshold of a division for a subject group was computed by averaging the corresponding metric values across regions belonging to that division. The aggregated metric values for threshold values from 0.1 to 0.6 were compared. We chose these thresholds because beyond 0.6 the network divides into disconnected subsets of nodes and small-world metrics are no longer meaningful . This analysis was performed using the Mplus software (http://www.statmodel.com). Growth curve models describe change (growth) with respect to a control variable. They are well-suited for analyzing group-level differences in biomedical data, particularly in cases where capturing and analyzing individual growth trajectories is important. Furthermore, for group comparisons, growth curve models alleviate the problem of multiple comparisons as fitted-curve coefficients are compared in contrast to traditional approaches where multiple individual points along the curve are compared. In our study, the growth trajectories of graph metric values of a subject carry important information about the variance within the group and needs to be incorporated in the model. The coefficients of growth curve models capture the baseline performance, instantaneous growth rate, and the acceleration of the variable of interest. We next examined degree, λ, E global, and γ values for each of the 90 anatomical ROIs, for the two groups, as a function of the correlation threshold. The metric values for threshold values from 0.1 to 0.6 in the two subject groups were compared using growth curve modeling, as described above. Analysis of Interregional Functional Connectivity Changes with Development To further characterize regional differences in network organization, we examined the regional connectivity at divisional level: association, limbic, paralimbic, primary, and subcortical. Differences in mean correlation coefficients for 4,005 pairs were aggregated into 15 pairs and the resulting differences were then normalized. (see also ). First, interregional pairs that showed statistically significant (p<0.01, FDR corrected) increased or decreased functional connectivity in young-adults group compared to child group were identified as (+1) or (−1), respectively. Second, the number of decreased (−1) or increased connectivities (+1) for each of the 15 pairs was counted. For example, to identify differential connectivity between the association division and the subcortical division, the number of decreased or increased connectivities between all pairs of subregions belonging to the association division and subcortical division was counted. Finally, since each brain region has a different number of subregions, the aggregated differential connectivity count was normalized by the number of possible connections between pairs of subregions belonging to the two divisions under investigation. Next we examined regional correlation values (connectivity) in the two groups. We compared regional correlation values aggregated across the 4,005 pairs of anatomical regions, between young-adults and children. No significant between-group differences in the aggregated correlation values were observed. On the basis of this observation, subsequently, individual regional correlation values were z-transformed followed by centering of the distribution around zero mean. These normalized correlation values were compared between the two subject groups. t-Test with a false discovery rate of 0.005 was used to test for significant differences. Analysis of Developmental Changes in Functional Connectivity We next examined the relationship between differences in regional correlation values (connectivity) in the two groups and the interregional wiring distance as determined using DTI. The wiring distance between two regions was computed by measuring the average length of the fiber tracks, in the MNI space, connecting those regions (see Text S1 for details). Supporting Information Figure S1 Functional connectivity in children and young-adults. Group averaged functional connectivity matrices for children and young-adults. Value of the (i,j)th element of the connectivity matrix corresponds to group averaged scale 3 wavelet correlation between the resting-state timeseries of brain region i and region j. Low correlation values are shown in darker color whereas high correlation values are shown in lighter color. Qualitatively, children, compared to young-adults, showed higher connectivity between the subcortical (caudate, globus pallidus, putamen, thalamus) and the cortical regions, and lower connectivity between the paralimbic (cingulate gyrus, orbitofrontal cortex, insula, parahippocampus gyrus, rectus gyrus, temporal pole) and the cortical regions. (2.70 MB TIF) Click here for additional data file. Table S1 Graph metrics for each anatomical region. (0.17 MB DOC) Click here for additional data file. Text S1 Experimental procedures. (0.07 MB DOC) Click here for additional data file.