During human brain development through infancy and childhood, microstructural and macrostructural changes take place to reshape the brain's structural networks and better adapt them to sophisticated functional and cognitive requirements. However, structural topological configuration of the human brain during this specific development period is not well understood. In this study, diffusion magnetic resonance image (dMRI) of 25 neonates, 13 toddlers, and 25 preadolescents were acquired to characterize network dynamics at these 3 landmark cross-sectional ages during early childhood. dMRI tractography was used to construct human brain structural networks, and the underlying topological properties were quantified by graph-theory approaches. Modular organization and small-world attributes are evident at birth with several important topological metrics increasing monotonically during development. Most significant increases of regional nodes occur in the posterior cingulate cortex, which plays a pivotal role in the functional default mode network. Positive correlations exist between nodal efficiencies and fractional anisotropy of the white matter traced from these nodes, while correlation slopes vary among the brain regions. These results reveal substantial topological reorganization of human brain structural networks through infancy and childhood, which is likely to be the outcome of both heterogeneous strengthening of the major white matter tracts and pruning of other axonal fibers.
Postnatal human brain development, especially from birth to the onset of adolescence, is characterized with both macroscopic and microscopic dynamic structural changes. Brain development during this period is associated with behavioral, cognitive, and emotional maturation. Structural alterations in the brain have important implications in the developmental process of brain functions. These structural changes include enhanced integrity of certain white matter (WM) axons and pruning of other axons (e.g., Innocenti 1981; Cowan et al. 1984; LaMantia and Rakic 1990, 1994; Woo et al. 1997; Innocenti and Price 2005), which may jointly contribute to significant reconfiguration of structural networks of developing brains. Specifically, axon pruning indicates the process of selective elimination of certain axons with or without the death of the parent cell during the period from birth to preadolescence.
The water molecules all over the human brain WM tend to diffuse more freely along the WM fiber bundle, instead of perpendicular to it. This diffusion property can be measured noninvasively with diffusion magnetic resonance image (dMRI) to infer the WM fiber orientations. This process is called dMRI-based tractography (e.g., Mori et al. 1999; Behrens et al. 2007), which allows us to noninvasively access the structural connectivity consisting of heterogeneous WM fiber bundles connecting different brain regions. Structural connectivity underlies the formation and evolution of human brain networks during development. Fractional anisotropy (FA; Pierpaoli and Basser 1996; Beaulieu 2002) derived from diffusion tensor imaging (DTI) (Basser et al. 1994), a type of dMRI, is capable of delineating the microstructure of WM fibers in the living human brain noninvasively. The microstructural and morphological changes of WM axons connecting different cortical areas of the developing brain have been extensively studied (e.g., Snook et al. 2005; Huang et al. 2006; Eluvathingal et al. 2007; Lebel et al. 2008; Tamnes et al. 2010; Lebel and Beaulieu 2011; Peters et al. 2012). Since the application of graph theory to studying human brain networks (Bullmore and Sporns 2009; He and Evans 2010; van den Heuvel and Sporns 2011), it has greatly advanced our understanding of the structural and functional circuits, the importance of which has been underscored by the establishment of the NIH Human Connectome Project (Sporns et al. 2005; NIH 2009). The development of functional brain networks (e.g., default mode network) or functional connectivity has been more extensively investigated with functional MRI (fMRI), especially resting-state fMRI (e.g., Fair et al. 2007; Supekar et al. 2009; Dosenbach et al. 2010; Power et al 2010; Hwang et al. 2012). Only recently have there been relatively limited studies on the development of human brain structural networks or structural connectivity constructed with dMRI tractography (Hagmann et al. 2010; Yap et al. 2011). However, so far there has been no comprehensive study investigating the structural networks covering 3 landmark cross-sectional ages during early childhood (birth, toddler, and onset of adolescence) or the relationship between structural networks and WM integrity. From birth to 2 years, the human brain undergoes several dramatic changes including rapid brain volume increases reaching 80–90% of adult volume by age 2 (Pfefferbaum et al. 1994), rapid elaboration of new synapses (Huttenlocher and Dabholkar 1997), overall very rapid gray matter volume increases (Huppi et al. 1998; Matsuzawa et al. 2001; Gilmore et al. 2007), rapid myelination of WM (Sampaio and Truwit 2001), and concurrent rapid development of a wide range of cognitive and motor functions (Kagan and Herschkowitz 2005). Preadolescence marks the beginning of a time period with dramatic changes in body and behavior. Hence, these 3 cross-sectional ages are critical for understanding brain development through infancy and childhood. With the coverage of these 3 landmark cross-sectional ages, it makes possible to determine whether the network configuration undergoes monotonic changes or whether the quantification of the network properties manifest peak values in the middle cross-sectional age.
In this study, we explored the development of human brain structural networks with high-quality dMRI data of 25 neonates, 13 toddlers, and 25 children around onset of adolescence, by applying graph theory for quantification of topological organization. High-quality dMRI data of 18 adults were also acquired to evaluate the maturity level of structural networks of preadolescents. The brain was parcellated into 80 regions or nodes with automated anatomical labeling (AAL; Tzourio-Mazoyer et al. 2002), and structural connectivity between these nodes was quantified with probabilistic tracking (Behrens et al. 2007). Topological changes during development were characterized with a comprehensive topological analysis with graph theory including computation of network metrics, such as network strength, global and local efficiency, modularity, and topological robustness. Furthermore, we investigated the relationship between network properties and microstructural integrity of WM (i.e., FA) for the developing brains. We aimed to test the hypotheses that the structural brain networks get monotonically stronger, more efficient, and more robust through infancy to childhood, and this network enhancement is associated with the increased microstructural integrity of the connected WM pathways. In addition, we also explored the structural basis for the development of the functional networks. Specifically, by surveying the nodal efficiency changes during this period, we could determine whether the important nodes involved in the functional default mode network, such as the posterior cingulate gyrus (PCG), are among those that undergo the largest increase.
Materials and Methods
dMRI and Structural MRI Acquisition for Neonates, Toddlers, and Preadolescents
Each individual volume of diffusion-weighted image (DWI) of all subjects was inspected through visual inspection by one of the authors (H.H.). After exclusion of corrupted data due to subject motion of 3 neonates and 1 toddler, high-quality dMRI and structural MRI data from 25 neonates (age range of 37–43 gestational weeks with mean and standard deviation of 39.5 ± 2.3 gestational weeks at scan; 13 M/12 F; 19 White and 6 African American), 13 toddlers (age range of 1.79–3.12 years old with mean and standard deviation of 2.3 ± 0.5 years at scan; 8 M/5 F; 13 Asian), and 25 preadolescents (age range of 10.7–13.5 years old with mean and standard deviation of 11.8 ± 1.8 years at scan; 13 M/12 F; 10 White, 3 African American, and 12 Asian) around puberty were used. To evaluate the maturity level of the structural brain networks in the preadolescent period, 18 healthy adult volunteers (age range of 25–44 years old with mean and standard deviation of 28.5 ± 5.1 years at scan; 13 M/5 F; 18 Asian) were also recruited. Due to the difficulty of reaching a good sample size for data from toddlers and preadolescents, collection of normal pediatric MR data from toddlers around 2 years- old and preadolescents around puberty was made possible by the collaboration of 2 institutions, Children's Medical Center at Dallas and Children's Hospital in Beijing. Philips 3-T Achieva MR scanners in the 2 institutions were used in this study. DTI and structural MRI data of all 25 neonates and 13 of 25 preadolescents around puberty were acquired from Children's Medical Center at Dallas. MR data of all 13 toddlers, 12 of 25 preadolescents, and 18 adults were acquired from Beijing Children's Hospital. The imaging protocols were identical except that different field of view (FOV) sizes were used for data acquisition of all participated subjects. All included neonates were part of the cohort for studying normal perinatal and prenatal development and were selected after rigorous screening procedures conducted by the neonatologist (L.C.). Exclusion criteria include mother's excessive drug or alcohol abuse during pregnancy; grade III–IV intraventricular hemorrhage; periventricular leukomalacia; hypoxic–ischemic encephalopathy; lung disease or brochopulmonary dysplasia; necrotizing enterocolitis that requires intestinal resection or complex feeding/nutritional disorders; defects or anomalies of the forebrain, brainstem, or cerebellum; brain tissue dys- or hypoplasias; abnormal meninges; alterations in the pial or ventricular surface; or WM lesions. Neonates were well fed before scanning. All included toddlers and preadolescents were healthy subjects with normal neurological and psychological records. For toddlers, scans were scheduled in the toddler's normal nap time. Besides earplugs and earphones, extra foam padding was applied to reduce the sound of the scanner while the toddler was asleep. dMRI data were acquired using a single-shot echo planar imaging with a SENSE parallel imaging scheme (SENSitivity Encoding, reduction factor = 2.5). Eight-channel SENSE head coil and consoles installed with the R2.6.3 software were used for both sites. The dMRI imaging parameters for all the 3 groups were: time echo (TE) = 78 ms, time repetition (TR) = 6850 ms, in-plane FOV = 168 × 168 to 224 × 224 mm2, in-plane imaging resolution = 2 × 2 mm2, slice thickness = 2 mm, slice number = 60–65 depending on the height of the subject brains, 30 independent diffusion-weighted directions (Jones et al. 1999) uniformly distributed in space, b-value = 1000 s/mm2, repetition = 2. The axial dMRI image dimension was 256 × 256 after reconstruction. For dMRI, the total acquisition time was 11 min. With 30 DWI volumes and 2 repetitions, we accepted those scanned dMRI datasets with <5 DWI volumes affected by motion more commonly seen in scanning of neonates and toddlers. The affected volumes were replaced by the good volumes of another DTI repetition during postprocessing. T1-weighted magnetization-prepared rapid gradient-echo (MPRAGE) image was also acquired. The imaging parameters for MPRAGE were: TR = 8.3 ms, TE = 3.8 ms, flip angle = 12°, voxel size = 1 × 1 × 1 mm3, FOV = 256 × 256 × 160 mm3, scan time = 4 min. The MPRAGE images provide superior gray and WM contrast and were used for segmentation and parcellation of the cerebral cortex of the toddler and preadolescent group. T1-weighted and dMRI images were acquired in the same session. All subjects gave informed written consent approved by the Institutional Review Board of the participated 2 institutions.
Network Construction for the Developing Brains
The workflow of structural network construction including cortical parcellation (Fig. 1a), dMRI tractography (Fig. 1b), and generation of connectivity matrix (Fig. 1c) of a typical neonate, toddler, and preadolescent brain from acquired dMRI and MPRAGE/b0 data can be found in Figure 1. After node and edge definitions described in details below, an 80 × 80 symmetric-weighted cortical network or graph representing the structural connectivity of the whole-brain was constructed for each subject, as shown in Figure 1c. Owing to the probabilistic nature of tractography in our study, the vast majority of regional pairs were assigned a nonzero probability. Unlike deterministic tracking, it is possible to include spurious connections with probabilistic tracking methods. To address this, we removed those obviously spurious connections that have extremely small probabilities by applying a threshold (Gong et al. 2009). The subsequent network analyses were repeatedly performed over the threshold range of 0.005 < wij < 0.05 with an interval of 0.0025 (19 thresholds). The comprehensive network analysis with a series of applied thresholds could help minimize the bias of the results obtained with a single threshold.
Cortical Parcellation for Network Node Definition
The procedure of defining the nodes has been previously described (Gong et al. 2009). The AAL labeling (Tzourio-Mazoyer et al. 2002) was obtained by transforming the AAL template to native space using a transformation matrix that was obtained by registering T1-weighted images of toddler and preadolescent groups and b0 images (with contrasts close to that of T2-weighted image) of the neonate group to ICBM152 template using SPM8 (http://www.fil.ion.ucl.ac.uk/spm). Two types of ICBM152 template (http://www.bic.mni.mcgill.ca/ServicesAtlases/ICBM152Lin), T1- and T2-weighted images of the ICBM152 template, were used for the preadolescent or toddler group and the neonate group, respectively. For the toddler and preadolescent groups, intrasubject affine registration was then conducted to align the AAL labeling with DTI images, resulting in 80 parcellated regions or nodes (40 for each hemisphere) for each brain. An example of transferred AAL labeling for a random subject in each group is shown in Figure 1a. Of note, discrete labeling values were preserved by the use of a nearest-neighbor interpolation method.
dMRI Tractography for Network Edge Definition
As shown in Figure 1b, with each of 80 nodes as seed, probabilistic tracking (Behrens et al. 2007) was performed by using FDT (FMRIB's Diffusion Toolbox) of FSL package (FSL, version 4.1; http://www.fmrib.ox.ac.uk/fsl). The connectivity probability from the seed voxel i to another voxel j was defined by the number of fibers passing through voxel j divided by the total number of fibers sampled from voxel i (Behrens et al. 2007). This idea could be extended from the voxel level to the regional level. For a seed region, 5000 × n fibers were sampled (5000 fibers for each voxel), where n is the number of voxels in this region. The number of fibers passing through a given region divided by 5000 × n is calculated as the connectivity probability from the seed region to this given region. In this study, each brain region was selected as the seed region and its connectivity probability to each of the other 79 regions was calculated. Notably, the probability from i to j is not necessarily equivalent to the one from j to i because of the tractography dependence on the seeding location. Thus, we defined the undirectional connectivity probability Pij between region i and j by averaging these 2 probabilities. Details on calculating Pij of each node can be found in the literature (Gong et al. 2009). Network topology could be affected by the factors such as brain size (e.g., Yan et al. 2011) and connection efficacy (Hagmann et al. 2010). Due to very different brain sizes and WM integrity of network pathways of the studied developing brains, the effects of brain size and connection efficacy were scaled by using wij = Pij×BrainSize/ADC as the weight between node i and j. ADC is the apparent diffusion coefficient of the traced WM and 1/ADC served as connection efficacy (Hagmann et al. 2010). The detailed discussion about the definition of edge weight can be found in Discussion section below.
A brain structural network (graph) G is composed by N nodes and K edges. To characterize topological organization of structural networks, the following graph measures under different thresholds were calculated: network strength, global and local efficiency, normalized shortest path length, normalized clustering coefficient, and small-worldness (Rubinov and Sporns 2010). We also investigated the modularity and robustness of the structural networks (Rubinov and Sporns 2010). For regional nodal characteristics, we considered the nodal efficiency of each node. Furthermore, we calculated the integrated area under the curve (AUC) for each network measure, which provides a summarized scalar for topological characterization of brain networks independent of single threshold selection. All network analysis was performed using the in-house GRETNA software.
Network Strength, Efficiency, and Small-Worldness
For a weighted brain network or graph, we first computed the network strength [eq. (1) in the Appendix] as the average of the strengths across all the nodes, in which the strength of a node is the sum of the edge weights (wij) linking to it. Then, we computed the weighted clustering coefficient [eq. (2) in the Appendix], which quantifies the extent of local interconnectivity or cliquishness in a network, and the weighted shortest path length [eq. (3) in the Appendix], which quantifies the ability for information propagation in a network in parallel. A network is said to be small-world if it has similar shortest path lengths but higher clustering coefficients than degree-matched random networks (Wattz and Strogatz 1998). In other words, a small-world network has the normalized clustering coefficient, , and the normalized shortest path length, . These 2 measurements can also be summarized into a simple quantitative measurement, small-worldness, [eq. (4) in the Appendix] (Humphries et al. 2005). We also computed the global efficiency as the inverse of shortest path length, the local efficiency as a measure of the efficiency of local information transfer between neighboring nodes, and the cost efficiency as the relative network efficiency normalized by its cost [eqs (5–7) in the Appendix, respectively]. For each node, we examined its nodal efficiency, which measures the average shortest path length between the given node and all the other nodes in the network [eq. (8) in the Appendix]. Besides, the normalized nodal efficiency was the nodal efficiency divided by global efficiency of the network. See Appendix for the detailed definitions and mathematical expressions of the network measures used in this work.
The modularity Q(p) for a given partition p of the brain network (Newman and Girvan 2004) is defined as:
Network robustness, characterized by the degree of tolerance against random failures and targeted attacks, is usually associated with the stability of a complex network. In the present study, we investigated the robustness (tolerance) of the WM structural networks by the removals of nodes (Albert et al. 2000; Achard et al. 2006). To address the nodal failure tolerance, we first randomly removed one node from the networks and then measured the changes in the size of the largest connected component (LCC). After this, we continued selecting and removing additional nodes from the networks at random and recomputed this measure. To evaluate the target attack tolerance, we repeated the above processes but removed the nodes in the descending order of their efficiency. Robustness is then usually visualized by a plot of the measure, LCC, versus the number of nodes removed, n (Achard et al. 2006). The robustness parameter (R) is defined as the area under this LCC (AUC of LCC) versus n curve. More robust networks retain a larger connected component even when several nodes have been knocked out, as represented by a larger AUC. The simulation procedures were performed for the individual network at each threshold T, and a comparison of the robustness parameters among 3 groups was further investigated.
Global and Regional Relationship Between Network Metrics and WM FA
Nodal FA and Whole-Brain WM FA
The individual nodal WM in this study was defined as the WM traced from a certain node. Using the traced WM as a binary mask which is mapped to the FA map (Pierpaoli and Basser 1996) derived from diffusion tensor, the averaged FA of the traced WM from a given node, namely the nodal FA, was calculated to characterize the microstructure of nodal WM. Voxels from spurious tractography which had an intensity of <0.05% of the highest value in the output connectivity distribution file of FSL were excluded from the binary mask. To ensure only WM FA was counted, an FA threshold was also applied. The FA threshold of 0.15 was adopted to include enough WM of neonates (Huang et al. 2006). For consistency of FA measurements among all the 3 groups, this lower FA threshold was applied to all the 3 age groups. FA threshold of 0.15 was relatively lower than the usually adopted 0.2 for the preadolescent group. However, our observation indicated that most of the non-WM voxels due to lower FA threshold in the preadolescent group were removed with the binary tractography mask described above. With the 2 filtering procedures described above, relatively clean and consistent nodal WM mask could be obtained for all age groups. The whole-brain WM FA was defined as the average of all nodal FA.
Global Relationship Between WM FA and Network Metrics
Network strength, global, and local efficiency of all subjects, representing the network properties of the entire brain, were correlated with the whole-brain WM FA of these subjects with Pearson's correlation.
Regional Relationship Between WM FA and Network Metrics
The node with the most significant nodal efficiency change and another node with the least significant nodal efficiency change were selected to characterize the regional relationship between nodal FA and nodal network metrics. Their absolute nodal efficiencies of all subjects were correlated with the corresponding nodal FA values with Pearson's correlations. Integrated AUC values of absolute nodal efficiency were used.
Quantitative Evaluation of Effects of Scanner Differences
To evaluate the effects of scanner differences, we have separated the preadolescent groups into 2 subgroups: one subgroup (Nsub1 = 13, subgroup 1) with data acquired from the Dallas site and another subgroup (Nsub1 = 12, subgroup 2) with data acquired from the Beijing site. The following 3 comparisons were made: (1) the between-group differences of integrated (AUC value) strength, global, local and cost efficiency, normalized Lp, normalized Cp, and small-worldness among the neonate group, toddler group, and subgroup 1 of preadolescents; (2) the between-group comparisons of above-mentioned parameters among the neonate group, toddler group, and subgroup 2 of preadolescents; (3) the between-group comparisons of above-mentioned parameters between subgroups 1 and 2 of preadolescents.
For group effects in global network measures and regional efficiency, comparisons were performed among 3 groups (neonates, toddlers, and preadolescents) using 1-way analysis of variance (ANOVA) with post hoc 2-sample t-tests when needed (P < 0.05 after correcting for multiple comparisons). Parametric ANOVA was used. For global network measures under different thresholds, P < 0.05 with Bonferroni correction (P < 0.05/19) was considered significant. Such Bonferroni correction was used for generating Figure 3a, Figure 4a, and left panels of Figure 6 below. Note that no correction was needed for statistical comparisons of integrated measures when generating Figure 3b and Supplementary Figure 1 below. Comparisons of nodal efficiency of the 3 groups were conducted with both normalized and absolute nodal efficiency. At the regional level, P < 0.05 with false discovery rate (FDR) correction for multiple comparisons across regions was considered significant when generating Figure 5 and Figure 8a,b below.
Testing Differences of Regional Correlations Between Absolute Nodal Efficiency and WM FA
Dynamics of Connectivity Matrix, Network Strength, Efficiency, and Small-World Properties During Brain Development
The averaged connectivity matrices of neonate, toddler, and preadolescent groups are shown in Figure 2a. While the connectivity pattern looks quite similar, the general stronger connectivity from the neonate to preadolescent group can be clearly observed from Figure 2a. The 3-dimensional view of the averaged network configuration of each age group is shown in Figure 2b. Stronger connectivity with increasing age, reflected by thicker blue connectivity path lines between the nodes, is visualized. As an overview of the edge weight distribution, the histogram of the edge weight for each age group is shown in Figure 2c. It is clear that larger edge weights, for example, those greater than 0.25, are more heavily distributed with age increase.
The network strength, efficiency, and small-world properties of the 3 age groups at different threshold and integrated values of these metrics are shown in Figure 3a,b, respectively. Quantitative values of integrated values of these global network properties (corresponding to Fig. 3b) are shown in Supplementary Table 1. From Figure 3a, the differences of all these network metrics at all thresholds from 0.005 to 0.05 are statistically significant with a P-value of <0.05 (Bonferroni-corrected), except normalized Lp. With large standard deviation of the normalized Lp, there is no statistical differences of this network metric among the 3 age groups at all thresholds (P > 0.1) (Fig. 3a). The network strength, global, local, and cost efficiency increase uniformly with age at all thresholds. Normalized Cp and small-worldness decrease uniformly with age at all thresholds. All the 3 groups show a small-world organization at all thresholds with a normalized Cp greater than 1 and normalized Lp close to 1 (Fig. 3a).
Correspondingly, the integrated values of connectivity strength, global, local, and cost efficiency in Figure 3b increase significantly (P < 0.001) and monotonically during brain development. The integrated values of normalized Cp and small-worldness decrease significantly and monotonically (P < 0.001) with age increase (Fig. 3b). The age-dependent and significantly decreased normalized Cp and small-worldness indicate that network segregation is reduced during development.
Network Modules of the Developing Human Brain
Figure 4a shows the number of modules, number of connectors, and modularity of the 3 developmental groups at different edge thresholds. The number of connectors and modularity are statistically different (P < 0.05, Bonferroni-corrected) among the 3 groups at most thresholds. The number of modules of the preadolescent group and that of the toddler group are similar (P > 0.1), both clearly higher than that of the neonate group at lower thresholds from 0.005 to 0.015. The number of connectors is higher with age increase at most thresholds, indicating more communications among the modules with brain development. The modularity of the neonate group with the edge weight threshold from 0.015 to around 0.05 is highest, while that of the preadolescent group lowest and that of the toddler group in the middle. With the edge weight thresholds from around 0.005 to 0.015, modularity of the toddler group is the largest, followed by that of neonate and preadolescent groups. It shows that the numbers of intramodule links of 3 age groups vary with different edge weight thresholds and intramodule links reduce at larger thresholds with development.
Three-dimensional representations of the modular configurations of the networks of neonate, toddler, and preadolescent groups are shown in Figure 4b. It can be appreciated that modular configurations of all the 3 groups are different despite that the general modular patterns are similar. Specifically, the medial fronto-parieto-occipital (purple) module is identified as a separate module in the toddler brain. However, it is combined into medial parieto-occipital (cyan) module and medial fronto-parietal (blue) module in the neonate and preadolescent groups, respectively. In addition, the left occipito-temporal module (yellow) is a separate module in the toddler and preadolescent groups, while it is integrated into the larger green module spanning frontal, parietal, temporal, and occipital cortices in the neonate group. The right occipito-temporal module (brown) is separate only in the preadolescent group.
Hub Distribution During Development
The hub distributions of the neonate, toddler, and preadolescent groups are shown in Figure 5a–c, respectively. Bilateral PCG, bilateral precuneus (PCUN), and right cuneus (CUN) are the common hubs for all the 3 groups. Bilateral dorsal cingulate gyrus (DCG) are hubs for both the neonate and toddler groups. Left anterior cingulate gyrus (ACG) and left superior occipital gyrus (SOG) become hubs only at the toddler and preadolescent group. The post hoc pair-wise comparisons of the normalized nodal efficiency of the 3 groups show that normalized nodal efficiencies change monotonically at the nodes, where significant (P < 0.01, FDR-corrected) differences were detected with the 1-way ANOVA test. Figure 5d demonstrates that significantly increased normalized nodal efficiencies take place at the bilateral cingulate, bilateral medial parietal, bilateral occipital, and left medial frontal cortical regions (including bilateral PCG, bilateral PCUN, bilateral CUN, bilateral calcarine fissure and surrounding cortex, bilateral SOG, left superior frontal gyrus, medial, right middle occipital gyrus, and right inferior occipital gyrus). Significantly decreased normalized nodal efficiencies are clustered at the left lateral temporal and parietal cortex during development. In Figure 5d, the nodal sizes indicate F values of ANOVA, and red and blue colors indicate the direction of monotonically increased and monotonically decreased normalized nodal efficiency, respectively. It is clear from Figure 5d that largest increase of normalized nodal efficiency takes place at bilateral PCG, which play a pivotal role in functional network development. To demonstrate the magnitude of the normalized efficiency changes, the values of the normalized efficiency of left Heschl's gyrus (HES) and right PCG, which has most significantly decreased (P < 0.01, FDR-corrected) and increased (P < 0.01, FDR-corrected) normalized nodal efficiency during development, are also shown in the left and right side of Figure 5d. It should be noted that although the normalized nodal efficiency can increase or decrease during development, the absolute nodal efficiency increases monotonically as demonstrated in Figure 8a below.
Topological Robustness of the Developmental Networks
Figure 6 shows the topological robustness properties of WM networks indicated by the values of the area under the ROC curve (AUC) of the LCC as a function of the removed node number by targeted attacks (Fig. 6a) or random failures (Fig. 6b). The preadolescent, toddler, and neonate groups are represented by red, blue, and green lines, respectively. Under both target attack and random failure, the AUC of LCC is highest in the preadolescent group, following in turn by the toddler and neonate groups with different edge weight thresholds (Fig. 6a,b). These differences become more apparent with higher edge weight thresholds (P < 0.05, Bonferroni-corrected). The examples at threshold 0.05 are highlighted by blue boxes in Figure 6a,b and shown in details in the right panels. With different numbers of removed nodes, the values of LCC of the preadolescent group are highest, followed in turn by the toddler and neonate groups in most node-removal situations (right panels of Fig. 6). The entire Figure 6 demonstrates that the brain is getting more robust to both target attack and random failure during both development periods from neonate to toddler and from toddler to preadolescent. In addition, the brain networks of the toddler group are approximately as robust as those of the preadolescent group in response to target attacks in the smaller range of edge weight threshold around 0–0.01 (P > 0.1) (Fig. 6a). The neonate group generally displays remarkably reduced stability against both targeted attack and random failure when compared with the other 2 groups with edge threshold from 0.025 to 0.05 and the number of removed nodes from 20 to 60, indicating special vulnerability of the neonate brain networks.
Global and Regional Relationship Between Network Metrics and WM FA
Figure 7 shows the statistically significant correlation (P < 0.001 for all global correlations) between the global network metrics (including network strength, global, and local efficiency) and whole-brain FA. It indicates a strong and positive global relationship between network strength and WM microstructural integrity and between network efficiencies and WM microstructural integrity.
Figure 8a shows the differences in absolute nodal efficiency of the 3 groups with the nodal size encoded by ANOVA F-values. The uniform red colors of all nodes in Figure 8a indicate the absolute nodal efficiency of every node increases monotonically (P < 0.0001, FDR-corrected) during development. Among all nodes, nodal efficiencies at the bilateral PCG show largest magnitude of increase. All nodal FA in Figure 8b also increases monotonically (P < 0.0001, FDR-corrected) with brain development. The statistically significant (P < 0.001) and positive correlations of the corresponding nodal efficiency and nodal FA are clear in Figure 8c, suggesting regional association of the increased nodal WM integrity and nodal efficiency. However, such correlations are heterogeneous among different nodes and whole brain. For example, with pair-wise tests, the statistically significant correlation slope differences (P < 0.001) were found between right PCG or left supramarginal gyrus (SMG) and the whole brain. The slope differences are also significant between right PCG and left SMG (P < 0.001). Figure 8c shows that the same magnitude of WM FA increases leads to heterogeneous nodal efficiency increases among different nodes and the whole brain.
Evaluation of the Maturity Level of Brain Structural Networks of Preadolescents with That of Adults as the Reference
Figure 9 shows the maturity level of brain structural networks of preadolescents with those of adults as a reference. Figure 9a indicates that the connectivity is getting even stronger from preadolescents to adults with a higher distribution of edge weights from 0.25 to 1 for adults. There is no significant difference of small-worldness, normalized Lp or normalized Cp between the 2 age groups, while global, local, and cost efficiency continues to increase from preadolescents to adults (Fig. 9b). Figure 9c shows that the hub distributions of the 2 groups are similar. For example, PCG-L/R and PCUN-L/R are the hubs in both groups. But in the adult group, the hubs at the anterior brain, for example ACG-L/R, get more prominent, whereas those in the right occipital lobe, for example right superior occipital gyrus and right cuneus, get diminished (Fig. 9c).
Evaluation of Scanner Effects
The preadolescent group was divided into 2 subgroups based on the data acquisition site. Comparisons of integrated values of the network metrics from each of these 2 subgroups with those from neonate and toddler groups were conducted to evaluate the scanner effects. The change patterns of integrated values of these network metrics are almost identical (Supplementary Fig. 1a, Supplementary Fig. 1b, and Fig. 3b), regardless of the preadolescent data from 1 of the 2 subgroups or from the entire preadolescent group. Furthermore, the comparisons of integrated values of the same network metrics between 2 preadolescent subgroups show no statistical significance (P > 0.05) for all network metrics, except integrated local efficiency (Supplementary Fig. 1c).
The human brain development from birth to adolescence is characterized by dramatic changes in connectivity and network configuration. With the coverage of 3 landmark cross-sectional ages from birth to the onset of adolescence, we found that the brain structural configuration was reshaped to be more efficient, stronger, more modularly organized, and more robust against attacks. The efficiency and strength increase continues into adulthood. Modular and hub reorganizations happen throughout the early childhood. Bilateral PCG, which play a pivotal role in functional network development, were found to undergo largest nodal efficiency increase in the structural network development. The heterogeneous nodal efficiency changes and relationship between the WM microstructural integrity and network nodal efficiency suggested that emergence of the maturing brain networks is underlined by both microstructural enhancement of some WM tracts and synaptic pruning of other fibers. To our knowledge, this study is the first comprehensive network investigation covering 3 landmark cross-sectional ages from neonates to adolescents with the weighted connectivity matrices and correlating network properties with WM microstructural integrity.
The Brain Structural Connections Grow Monotonically Stronger and More Efficient
The network quantification shows a general stronger and more efficient connectivity with development. The monotonic increase of network strength, global, and local efficiency indicates that the brain networks are consistently getting more efficient. On the other hand, monotonic decrease of clustering coefficients and small-worldness suggests that segregation decreases consistently with brain maturation from birth to onset of adolescent (Fig. 3). Figure 8a also reveals that the nodal efficiency of all nodes increases monotonically and significantly. This monotonic increase of network property during development is accompanied with that of microstructural property (Fig. 8b) also observed in the previous studies (e.g., Snook et al. 2005; Eluvathingal et al. 2007; Lebel et al. 2008; Tamnes et al. 2010). The dMRI-based binary network study (Yap et al. 2011) on infant brains has revealed that the small-worldness exists in the very early life of human brains. This small-worldness property was also found in our study of the neonate brain right around birth at 40 gestational weeks, even 2 weeks earlier than the youngest age reported in that literature (Yap et al. 2011). The dMRI-based study (Hagmann et al. 2010) on the developing human brain from 2 to 18 years old found increased efficiency, node strength, and reduced modularity of the structural networks during development through weighted network analysis. The dynamics of these network metrics are consistent with our results except that these monotonic network changes are extended further to include 0 to 2 years old in our study. Hence, our results suggest that the integration increase and segregation decrease of the human brain networks start as early as birth and continue till the onset of adolescence. These network processes are likely to be monotonic and contribute to the maturation of the human brain configuration. It should be noted that not all developmental procedures from birth to the onset of adolescence are associated with monotonic processes. For example, the previous functional study measuring local cerebral metabolic rates for glucose with positron emission tomography indicated the measured values rose rapidly to adult values by 2 years, reached a plateau to 3 and 4 years, maintained until 9 years, and began to decline then (Chugani et al. 1987). With incorporation of the 3 important developmental landmark cross-sectional ages in our study, the monotonicity of these network metric changes from birth to the onset of adolescence suggests that the brain networks are reconfigured to continuously get stronger and more efficient during this period.
Strengthened Networks with Development Involves not only Microstructural Enhancement of WM, but also Other Factors Possibly Including Pruning
Studies on development of functional connections (e.g., Fair et al. 2007,, 2008; Supekar et al. 2009; Dosenbach et al. 2010; Power et al 2010; Hwang et al. 2012) have suggested that the formation of the human brain networks is associated with some important structural developmental events taking place in parallel. Specifically, these events include myelination of certain WM fibers (e.g., Yakovlev and Lecours 1967; Benes et al. 1994), pruning of other WM fibers (e.g., Innocenti 1981; Cowan et al. 1984; LaMantia and Rakic 1990; LaMantia and Rakic 1994; Woo et al. 1997; Innocenti and Price 2005), and integration through synchronization (Varela et al. 2001). The structural connectivity and functional connectivity are linked to each other (e.g., Greicius et al. 2009; Honey et al. 2009; van den Heuvel et al. 2009). Hence, the factors reshaping the functional connectivity during brain maturation may also reshape structural connectivity. However, unlike developmental networks established with functional connectivity, the structural networks established with dMRI tractography in this study include only connectivity with axonal connections. The 2 regions connecting each other through a relay at a third region can be detected to have connection in the networks built with correlation analysis of resting-state fMRI, but not in the networks built with tractography of dMRI. Therefore, the integration through synchronization cannot be apparently manifested in the dMRI-based developmental networks in our study.
Overproduction, elimination, and reshaping of corpus callosum (Innocenti 1981; LaMantia and Rakic 1990), anterior commissure (LaMantia and Rakic 1994), and prefrontal associational circuitry (Woo et al. 1997) have been found in previous studies with monkey or cat models. Synaptic overproduction in infancy, persistence of high levels of synaptic density to adolescence and decrease after adolescence, has also been observed in the human brain (Chugani et al. 1987; Huttenlocher and Dabholkar 1997; Petanjek et al. 2011). These time-dependent changes vary across visual, auditory, and prefrontal cortices (Huttenlocher and Dabholkar 1997). The following specific findings may be related to axon pruning. First, from Figure 8c, although positive and significant correlations exist between nodal efficiencies and FA, the correlation slopes are significantly different. Of note, these FA values represent the microstructural integrity of the WM axons connected with the nodes. The largest increase of nodal efficiency takes place at PCG-R. The major WM tract connecting PCG is the cingulum bundle in cingulate gyrus. To our knowledge, there is no direct evidence so far showing cingulum bundle in cingulate gyrus undergoes most rapid increase of microstructural integrity among all major WM tracts from birth to preadolescence. It suggests that other factors, possibly including pruning of small fibers or axonal branches, in addition to the microstructural enhancement of major fiber bundles, contribute to the variability of nodal efficiency change rates among different nodes. Previously, we found that major WM tracts including the association tracts and commissural tracts connecting different cortical regions are well formed in the neonate brains, and the general morphological pattern is quite similar to that of the 5- to 6-year-old children (Huang et al. 2006). The microstructural enhancement of major WM tracts is heterogeneous during development in childhood (e.g., Lebel et al. 2008; Lebel and Beaulieu 2011). Not only increased but also decreased normalized efficiency observed in Figure 5d is probably associated with axon pruning and heterogeneous microstructural enhancement of major WM tracts, but not with formation of new major WM tracts. With Figure 5d and Figure 8c, it is likely that the strengthening of the growing structural networks is the outcome of both the heterogeneous microstructural enhancement (including myelination) of the major WM tracts and pruning of the small fibers. The analysis of integrated effects of microstructural enhancement of long-range fibers and pruning of short range fibers in the future may provide further insights on the underlying mechanism of growing human brain networks.
Formation of the Modular Organization and Reconfiguration of the Brain Hubs From Birth to Adolescence
During development, there is a general trend of increases for the numbers of modules and connectors, indicating emergence of the modular organization and more communications among the modules with brain development (Fig. 4a). Moreover, at edge weight threshold from around 0.015 to 0.05, the modularity decreases during development, indicating loss of intramodule links. This reduction of the intramodule links is possibly through the pruning process during development. The 3 plots in Figure 4a jointly suggest the strengthening of intermodule connections and the weakening of intramodule connections. We also identified modular organization of the brains at each age group (Fig. 4b). The modular organizations are different among the 3 age groups (Fig. 4b). There is no clear monotonic tread regarding with the modular organization during development. These findings suggest that the activity levels of 2 brain development processes discussed previously, major WM tract strengthening and small fiber pruning, are varied during different developing periods. It has been found that more pronounced synaptic elimination occurs in later childhood (e.g., Huttenlocher and Dabholkar 1997).
The differences of normalized nodal efficiency revealed that the largest increases of normalized nodal efficiency take place at the bilateral PCG (Fig. 5d). It is consistent to the microstructural integrity increase of the cingulum bundle in the cingulate gyrus (the major WM tract connecting PCG) during childhood observed in previous studies (e.g., Lebel et al. 2008; Lebel and Beaulieu 2011). We also found that the bilateral PCG are structural network hubs in all the 3 landmark cross-sectional ages (Fig. 5a–c) and remain to be the hubs in adulthood (Fig. 9c). PCG plays a pivotal role in default mode network in resting-state functional connectivity (e.g., Greicius et al. 2003). The relationship between structural and functional connectivity has been confirmed previously (e.g., Greicius et al. 2009; van den Heuvel et al. 2009). Our finding on PCG, therefore, provides evidence for the structural basis of the formation of the default mode network (Fair et al. 2008; Gao et al. 2009) during brain maturation. Default mode network supports mental activities such as episodic memory (Frith and Firth 1999), mentalizing (Frith and Frith 2003), and self-projection (Buckner and Carroll 2007). At early school age (7- to 9 years old), the default regions are only sparsely functionally connected; these regions integrate into a cohesive and interconnected network over development (Fair et al. 2008). Emergence of the default mode network has also been shown in neonate brains (Doria et al. 2010). As a key component of default mode network, we speculate that the prominent increase of PCG efficiency from birth to preadolescence (Fig. 5d and Fig. 8c) in the structural networks contributes to the formation of a cohesive functional default mode network (Fair et al. 2008). It is not clear why other key components of the default mode network do not manifest remarkable efficiency increases. However, it has been reported that PCG is a highly connected and metabolically active brain region (Leech and Sharp 2013). Of note, highly structurally connected PCG has been consistently found in adult brains in previous dMRI studies (Hagmann et al. 2008; Gong et al. 2009; van den Heuvel and Sporns 2011), and our study also demonstrated that it continues to be the structural network hub in adulthood (Fig. 9c).
Growing Robustness of the Brain Structural Networks with Development
Another intriguing finding is the growing robustness of the brain structural networks through topological tests after target attack or random failure. These tests of robustness against attack in the framework of developing brain revealed interesting topological advantages of the networks of the more mature brains. Figure 6 shows that, from neonates to adolescents, the younger the age, the more vulnerable the brain networks are against target attack or random failure. Increased robustness with early brain development may be related to an increased number of connectors (Fig. 4a), which play an important role in keeping the robustness and stability of the brain networks (He et al. 2009). This relatively poor topological robustness in younger brains has significant implication of autism, which develops during infancy and early childhood. Specifically, it has been speculated that excess neuron numbers may be one possible cause of early autistic brain overgrowth and produce defects in neural patterning and wiring, with exuberant local and short-distance cortical interactions impeding the function of large-scale, long-distance interactions between brain regions (e.g., Courchesne et al. 2007). Such disruptions of large-scale connections located in frontal and temporal regions and associated with corpus callosum have been observed in several studies investigating autism with dMRI (e.g., Barnea-Goraly et al. 2004; Weinstein et al. 2011; Lewis et al. 2013). These connections underlie socio-emotional and communication functions and alterations of these connections could relate to the early clinical manifestations of autism. Connector nodes identified by graph theory interconnect modules (Guimera et al. 2004; Sporns et al. 2007; Rubinov and Sporns 2010) and are likely to be associated with long-range connections. It is reasonable to extrapolate that disruptions of long-range connectivity in the autistic brain may impede increase of connectors in a normal developmental process (see Fig. 4a) and hence yield a less robust network. The other implication of this network robustness test is plasticity. Human brains become less plastic during development. It seems that network robustness and plasticity are inversely correlated to each other. The early developing brain has a poor topological robustness but more plasticity. A developing brain gives up some of its plasticity in favor of efficiency and stability (Vogel 2012), which has been jointly illustrated by our results of network efficiency (Fig. 3) and robustness (Fig. 6).
Technical Issues, Limitations, and Further Considerations
There are several technical issues related to edge weight in the weighted network analysis used in this study. As argued in a previous structural network study (Gong et al. 2009), it is impossible to choose a single threshold with the probabilistic tracking used in this study. Hence, network property measurements (Fig. 3a), modular quantifications (Fig. 4a), and network robustness measurements (Fig. 6) were conducted at different edge weight thresholds. Higher thresholds usually indicate less spurious tracing results included. However, there are also risks of removing “real” connectivity by increasing thresholds. Most of the between-group comparisons (Fig. 3a and Fig. 6) are consistent with different thresholds, while modular quantifications are more sensitive to the changes of thresholds (Fig. 4a). From birth to adolescence, the brain undergoes dramatic microstructural and macrostructural changes, both of which need to be taken as the edge weight factors. To account for the microstructural changes, 1/ADC used by Hagmann et al. (2010) was adopted as one of the edge weight scalars. For a brain with larger volume, the connectivity probability Pij tends to be smaller after WM fibers are traced for longer distance. To make up for the smaller Pij in a larger brain, brain size is another edge weight scalar. These factors were used for the fair weighted network comparisons among the 3 age groups. Currently, there has been no consensus of the weight definition for the edges derived from dMRI tractography. The variations in the brain size and microstructures of WM during development add complexity. To evaluate the effects of the factor of 1/ADC on network property changes during maturation process, we also conducted network analysis without this factor. As shown in Supplementary Figure 2, the significant differences are still preserved and order of the differences is the same as that shown in Figure 3a.
Owing to difficulty of obtaining good sample number, the datasets were acquired from 2 sites, but with the same imaging protocol and same Philips 3T scanners. Rigorous quality control of both scanners is conducted routinely. To evaluate the scanner effects, we have divided the preadolescents into 2 subgroups depending on the data acquisition site. Compared with Figure 3b, it is clear from Supplementary Figure 1a,b that the scanner difference does not impose noticeable effects on the network measurements in that the same change patterns with statistical significance can be found for all measured network metrics. Furthermore, there are no statistical differences of these network metrics between the 2 preadolescent groups (Supplementary Fig. 1c) except integrated local efficiency, which is likely to be caused by the individual differences of 2 subgroups of samples rather than the scanner differences. The FA measurements could be affected by the scanner difference too. In another study by our group, a healthy young subject (“in vivo human phantom”) was scanned in 2 Philips 3T scanners (same type of scanner as those used in this study) at 2 sites with the same dMRI sequence. Quantitative FA measurement differences caused by scanner difference were tested to be within the range of variability of scanning the same subject twice with one scanner (Saxena et al. 2012). It should also be noted that even the MR data obtained by scanning the same subject with the same scanner twice could have slight differences due to constant subtle fluctuations of the scanner. With the same type of scanners used in this study and rigorous quality control of both scanners, the effects of scanner differences may lead to relatively subtle changes of the presented results, but cannot be big enough to affect the statistical significance and overall trend of the network configurations during development found in this study. Furthermore, with the subgroup 1 of the preadolescent group from children with a relatively diverse races and the subgroup 2 exclusively from Asian children, the similar network metrics shown in Supplementary Figure 1c and dramatic differences shown in Figure 3b jointly suggest that the effects of genetic background, socio-cultural, or educational environment on brain structural networks are relatively trivial compared with those caused by normal brain development.
It should be noted that this study is not a longitudinal one as the subjects at each cross-sectional age were not the same. It could also be one limitation that the sample size (n = 13) of the toddler's group is less than those (n = 25) of the other 2 age groups due to greater difficulty of subject recruitment and data acquisition of toddlers. However, shown in Figure 3 and Figure 4, the standard deviations of the network metrics of toddlers are comparable with those of the other 2 age groups, suggesting the uneven sample numbers across the age groups have relatively limited effects on group comparisons of the network metrics.
Besides dMRI, the brain structural networks can also be investigated by measuring inter-regional covariance of gray matter morphological features (e.g., thickness or volume) derived from T1-weighted data (He et al. 2007; Bassett et al. 2008; for review, see Alexander-Bloch, Giedd, et al. 2013). Several recent studies have suggested that such structural covariance is likely to reflect developmental coordination or synchronized maturation in anatomically related regions (Raznahan et al. 2011; Alexander-Bloch, Raznahan, et al. 2013) and exhibits shared global and local topological organization with dMRI WM networks (Gong et al. 2012). Notably, a recent study (Fan et al. 2011) reported that the structural covariance networks, even in neonates, already display high efficient small-world properties and significant modular organization, and that these properties further increase during the first 2 years of life. The developmental patterns of structural covariance are highly compatible with our findings (e.g., small-worldness in Fig. 3 and modular organization in Fig. 4) based on dMRI data. Further studies combining large sample dMRI and structural MRI data would be interesting to explore developmental changes in the brain's structural covariance and WM connectivity networks.
By quantitatively characterizing the brain network properties of neonates, toddlers, and preadolescents, we found a significant and monotonic increase of network efficiency and strength through infancy and childhood with the brain's structural networks more modularly organized and more robust against attacks. We observed that the brain maturation process in this period is accompanied with heterogeneous regional nodal efficiency increases, and such increases are significantly correlated with the microstructural enhancement of the associated WM. Highest efficiency increase takes place at PCG, the pivotal component of functional default mode network and the hub of structural networks at all 3 cross-sectional ages during development and at the age of adulthood. The characterized network reconfiguration during this period could be contributed by heterogeneous enhancement of some axons and pruning of others. These findings may provide insight into pathologies of development such as autism and potentially allow us to define structural network correlates for human cognitive and behavior development.
This study is sponsored by NIH (Grant nos. EB009545, MH092535, and MH092535-S1), 973 Program (Grant no. 2013CB837300), the Natural Science Foundation of China (Grant nos. 31221003, 81000633, 81030028 and 31000499), the National Science Fund for Distinguished Young Scholars of China (Grant no. 81225012), and the 111 Project (Grant no. B07008).
Conflict of Interest: None declared.
Small-world network parameters (clustering coefficient, Cp, and shortest path length, Lp) were originally proposed by Watts and strogatz (1998). In this study, we investigated the small-world properties of the weighted brain networks. The clustering coefficient of a node i, C(i), which was defined as the likelihood whether the neighborhoods were connected with each other or not, is expressed as follows:
The path length between any pair of nodes (e.g., nodes i and j) is defined as the sum of the edge lengths along this path. For weighted networks, the length of each edge was assigned by computing the reciprocal of the edge weight, 1/wij. The shortest path length, Lij, is defined as the length of the path for nodes i and j with the shortest length. The shortest path length of a network is computed as follows:
To examine the small-world properties, the clustering coefficient, Cp, and shortest path length, Lp, of the brain networks were compared with those of random networks. In this study, we generated 100 matched random networks, which had the same number of nodes, edges, and degree distribution as the real networks (Maslov and Sneppen 2002). Of note, we retained the weight of each edge during the randomization procedure such that the weight distribution of the network was preserved. Furthermore, we computed the normalized shortest path length (lambda), , and the normalized clustering coefficient (gamma), , where and are the mean shortest path length and the mean clustering coefficient of 100 matched random networks. Of note, the 2 parameters correct the differences in the edge number and degree distribution of the networks across individuals. A real network would be considered small-world if and (Watts and Strogatz 1998). In other words, a small-world network has not only the higher local interconnectivity, but also the approximately equivalent shortest path length compared with the random networks. These 2 measurements can be summarized into a simple quantitative metric, small-worldness,
The global efficiency of G measures the global efficiency of the parallel information transfer in the network (Latora and Marchiori 2001), which can be computed as:
The local efficiency of G reveals how much the network is fault tolerant, showing how efficient the communication is among the first neighbors of the node i when it is removed. The local efficiency of a graph is defined as:
The network cost was to measure the expense for building up the connecting elements of a graph. Typically, the cost of a connection is proportional to its distance, which is the inverse of the edge weight in this study, and therefore, the overall cost of a graph is derived by taking the sum of distance, that is, Cost = sum (1/wij). The cost efficiency measures the relative network efficiency normalized by its cost.
Regional Nodal Characteristics
To determine the nodal (regional) characteristics of the WM networks, we computed the regional efficiency, Enodal(i) (Achard and Bullmore 2007):
Conflict of Interest: None declared.