*3 Department of Radiology, Capital Medical University, Beijing, China*

*4 State Key Laboratory of Cognitive Neuroscience and Learning & IDG/McGovern Institute for Brain Research, Beijing Normal University, Beijing, China*

*5 Department of Radiology, University of Texas Southwestern Medical Center, Dallas, TX, USA*

#### *Edited by:*

*Alan Evans, McGill University, Canada*

#### *Reviewed by:*

*Pratik Mukherjee, University of California, USA Xiaobo Li, Albert Einstein College of Medicine, USA*

#### *\*Correspondence:*

*Hao Huang, Advanced Imaging Research Center, University of Texas Southwestern Medical Center, 5323 Harry Hines Blvd., Dallas, TX 75390-8542, USA e-mail: hao.huang@ utsouthwestern.edu*

The human brain development is a complicated yet well-organized process. Metrics derived from diffusion tensor imaging (DTI), including fractional anisotropy (FA), radial (RD), axial (AxD), and mean diffusivity (MD), have been used to noninvasively access the microstructural development of human brain white matter (WM). At birth, most of the major WM tracts are apparent but in a relatively disorganized pattern. Brain maturation is a process of establishing an organized pattern of these major WM tracts. However, how the linkage pattern of major WM tracts changes during development remains unclear. In this study, DTI data of 26 neonates and 28 children around puberty were acquired. 10 major WM tracts, representing four major tract groups involved in distinctive brain functions, were traced with DTI tractography for all 54 subjects. With the 10 by 10 correlation matrices constructed with Spearman's pairwise inter-tract correlations and based on tract-level measurements of FA, RD, AxD, and MD of both age groups, we assessed if the inter-tract correlations become stronger from birth to puberty. In addition, hierarchical clustering was performed based on the pairwise correlations of WM tracts to reveal the clustering pattern for each age group and pattern shift from birth to puberty. Stronger and enhanced microstructural inter-tract correlations were found during development from birth to puberty. The linkage patterns of two age groups differ due to brain development. These changes of microstructural correlations from birth to puberty suggest inhomogeneous but organized myelination processes which cause the reshuffled inter-tract correlation pattern and make homologous tracts tightly clustered. It opens a new window to study WM tract development and can be potentially used to investigate atypical brain development due to neurological or psychiatric disorders.

 The major cerebral white matter (WM) tracts connecting different brain regions are involved in different brain functions. These major cerebral WM tracts are often categorized into different tract groups based on their distinct functions. There are roughly four tract groups, namely limbic, projection, callosal, and association tract groups (e.g., Wakana et al., 2004; Huang et al., 2012a,b), for cerebral WM tracts. The WM tracts within a tract group perform similar functions. For example, limbic tracts underlie the connectivity in the limbic system and association tracts connect between cerebral cortical areas. The pair of tracts in both cerebral hemispheres belongs to the same tract group and is considered as homologous tracts. At birth, most of major WM tracts are well formed (e.g., Huang et al., 2006; Oishi et al., 2011), except the arcuate fasciculus which is a part of the superior longitudinal fasciculus (SLF) and related to language function.

The water molecules all over the human brain tend to diffuse more freely along the WM fiber bundle, instead of perpendicular to it. This diffusion property can be measured noninvasively with diffusion MRI (dMRI), a modality of MRI. The widely used diffusion tensor imaging (DTI) (Basser et al., 1994) characterizes the water diffusion properties in the brain voxels with a tensor model. Fractional anisotropy (FA) (Pierpaoli and Basser, 1996; Beaulieu, 2002) and mean diffusivity (MD), derived from DTI, have been widely used to quantify the microstructural properties of the WM voxels. Other than FA or MD, the other two DTI-derived metrics, radial diffusivity (RD) and axial diffusivity (AxD), convey unique information related to myelination and axonal integrity, respectively (Song et al., 2002). The four DTI-derived metrics, FA, RD, AxD, and MD, characterize different aspects of diffusion tensor and are highly sensitive to WM microstructural changes.

Compared to voxel-based morphometry (VBM), recent tract analyses (Yushkevich et al., 2008; Goodlett et al., 2009; O'Donnell et al., 2009; Zhang et al., 2010; Colby et al., 2012) including ours (Huang et al., 2011, 2012a,b) have become important due to great functional and clinical significance of the tracts.  And the heterogeneous WM tracts can be noninvasively segmented with the traced fibers. With these segmented WM tracts as binary masks for the maps of DTI-derived metrics, the microstructural properties of the WM tracts can be quantified.

Dramatic microstructural changes take place during normal human brain development from birth to puberty, which are two landmark time points in early brain development. Birth marks the beginning time point of postnatal brain development. Puberty marks the end of the child development and beginning of adolescence. The four DTI derived metrics, FA, RD, AxD and MD, have been incorporated in numerous studies investigating WM microstructural changes for infants, children and adolescents during development.  On the other hand, it was shown with DTI-derived metrics of a cohort of adults that specific WM tracts involved in similar functions vary in a similar pattern with each other across different individuals (Wahl et al., 2010; Li et al., 2012), while hemispheric asymmetries of DTI-derived metrics in homologous pairs of WM tracts (Bonekamp et al., 2007; Wilde et al., 2009) have been reported. However, from perspective of brain development, whether or not the significant inter-tract correlations exist at birth and around puberty is still unclear. Furthermore, it remains elusive if these inter-tract correlations will be strengthened and how the correlation patterns change from birth to puberty.

During the development from birth to puberty, the human brain is likely to change from a more randomized state to a more balanced and organized state. In this study, we hypothesized that inter-tract correlations become stronger and the correlation patterns are reshuffled from birth to puberty. Specifically, the reshuffling process will cause more homologous tracts to form tight relationship. DTI data were acquired from 26 normal neonates and 28 normal children around puberty. The following 10 major WM tracts covering limbic, association, commissural and projection tract groups were selected for tract-level measurements of DTI metrics of each subject: left and right corticospinal tract (CST\_L and CST\_R), left and right inferior fronto-occipital fasciculus (IFO\_L and IFO\_R), left and right cingulate part of cingulum tract (CGC\_L and CGC\_R), left and right hippocampal part of cingulum tract (CGH\_L and CGH\_R), forceps major (FMajor) and forceps minor (FMinor). The tract level comparisons of all four DTI-derived metrics were conducted between the two age groups. Spearman's pairwise inter-tract correlations were performed. We tested if significant correlations of homologous WM tracts exist in neonates and children around puberty. After obtaining four 10 by 10 inter-tract correlation matrices corresponding to four DTI-derived metrics, FA, RD, AxD, and MD, for each age group, we tested these correlation matrices against the identity matrix or a matrix with equal non-diagonal entries. We then assessed if the inter-tract correlations become statistically stronger from birth to puberty. In addition, hierarchical clustering was performed with the pairwise correlations based on FA, RD, AxD, and MD measurements for each age group to reveal the pattern of clustering in either age group and reveal the pattern shift from birth to puberty.

## **MATERIALS AND METHODS**

## **SUBJECTS AND DATA ACQUISITION**

Twenty six normal neonates (14 males; age: 37 to 43 gestational weeks with mean and standard deviation 40*.*1 ± 2*.*0 gestational weeks) and 28 normal children around puberty (15 males; age: 9.5–15 years with mean and standard deviation 12*.*0 ± 2*.*3 years), free of current and past neurological or psychiatric disorders, were recruited at Children's Medical Center (CMC) at Dallas and Advanced Imaging Research Center (AIRC) of the University of Texas Southwestern Medical Center (UTSW), respectively. The parents of all the subjects gave written informed consents approved by Institutional Review Board of UTSW.

Two 3T Philips Achieva MR systems at CMC and AIRC were used to acquire dMRI of neonate and child group, respectively. dMRI data were acquired using a single shot echo planar imaging (EPI) with SENSE parallel imaging scheme (SENsitivity Encoding, reduction factor = 2.3). dMRI parameters for neonates were: *FOV* = 200/200/100 mm, in-plane imaging matrix = 100 × 100, axial slice thickness = 2 mm. dMRI parameters for children around puberty were: *FOV* = 224/224/143 mm, in-plane imaging matrix = 112 × 112, axial slice thickness = 2*.*2 mm. The common parameters for dMRI acquisition of both neonate and child group were: *b*-value = 1000 s/mm2, *TE* = 97 ms, *TR* = 7.6 s, 30 independent diffusion-weighted directions (Jones et al., 1999b) and 2 repetitions to increase signal-to-noise ratio (SNR).

## **DTI PREPROCESSING**

dMRI acquired from all the subjects was processed offline using DTIStudio (mristudio.org; Jiang et al., 2006). dMRI images for each subject were corrected for motion and eddy current by registering all the diffusion weighted images to the b0 image using a 12-parameter (affine) linear image registration with automated image registration (AIR) algorithm (Woods et al., 1998). After the registration, six independent elements of the 3 × 3 diffusion tensor (Basser et al., 1994) were determined by multivariate least-square fitting of diffusion weighted images. The tensor was diagonalized to obtain three eigenvalues (λ1 − 3) and eigenvectors (ν1 − 3). FA, MD, AxD and RD, derived from DTI, were obtained for all the subjects with the following equations of eigenvalues:

$$\begin{aligned} FA &= \frac{\sqrt{(\lambda\_1 - \lambda\_2)^2 + (\lambda\_1 - \lambda\_3)^2 + (\lambda\_2 - \lambda\_3)^2}}{\sqrt{2}\sqrt{\lambda\_1^2 + \lambda\_2^2 + \lambda\_3^2}} \\ MD &= (\lambda\_1 + \lambda\_2 + \lambda\_3)/3 \\ AxD &= \lambda\_1 \\ RD &= (\lambda\_2 + \lambda\_3)/2 \end{aligned}$$

## **TRACT-LEVEL MEASUREMENTS OF DTI METRICS**

The following 10 major WM tracts were selected for tract-level measurements of DTI metrics, left and right corticospinal tract (CST\_L and CST\_R), left and right inferior fronto-occipital fasciculus (IFO\_L and IFO\_R), left and right cingulate part of cingulum tract (CGC\_L and CGC\_R), left and right hippocampal part of cingulum tract (CGH\_L and CGH\_R), forceps major (FMajor) and forceps minor (FMinor). These tracts could be reproducibly traced with DTI of all neonates and children and cover all four major tract groups, namely projection, limbic, commissural and association tract group. Other major WM tracts such as SLF could not be traced reproducibly with our cohort of neonate DTI dataset. Following the literature (Wakana et al., 2007), the tractography protocol described in details below was used to trace all these tracts. DTIstudio (mristudio.org) was used to conduct the tractography. The binary masks of the individually traced tracts were used to compute the tract-level FA, RD, AxD, and MD. The test-retest reliability was quantified by coefficient of variation (CV) and κ values of variability shown in Supplemental Table 1, after tracing the tracts below 3 times with the data from 3 subjects randomly selected from each group. All CV values are less than 2% and κ values are greater than 95% for both neonate an child group, indicating almost perfect test-retest reliability and almost perfect agreement of measurements among different tests in both neonate and child group.

## *CST-L and CST-R*

For the first ROI, the entire cerebral peduncle of the desired hemisphere was delineated at the level of the decussation of the superior cerebellar peduncle using an axial slice. "OR" operation was used to select all the CST fibers in this hemisphere that reach the primary motor cortex. The second ROI was then drawn at the most ventral axial slice that identifies the cleavage of the central sulcus. "AND" operation is performed at this axial slice to select all the CST fibers in this hemisphere. The fibers running through to the opposite hemisphere were removed using the "NOT" operation.

#### *CGC-L and CGC-R*

For the first ROI, a coronal plane was selected at the middle of the splenium of the corpus callosum (CC) using the mid-sagittal plane and the region containing the entire cingulum in the desired hemisphere is selected. All the fibers in this coronal plane passing through the cingulum were selected using the "OR" operation. The second ROI was drawn by selecting a coronal plane in the middle of the genu of the CC and all the CGC fibers were selected using the "AND" operation.

## *CGH-L and CGH-R*

For the first ROI, a coronal plane in the middle of the splenium of the CC was selected using the mid-sagittal plane and the cingulum below the CC of the desired hemisphere was delineated. All the fibers in this coronal plane passing through the cingulum were selected using the "OR" operation. The second ROI was drawn at a coronal slice anterior to the pons using the mid-sagittal plane and the fibers passing through the cingulum in this hemisphere were selected using the "AND" operation.

## *IFO-L and IFO-R*

For the first ROI, a coronal slice at the middle point between the posterior edge of the cingulum and the posterior edge of the parieto-occipital sulcus was selected and the entire occipital lobe of the desired hemisphere was delineated. All the fibers in this hemisphere were selected using the "OR" operation. The second ROI was drawn at the anterior edge of the CC using a coronal slice and all the fibers in this hemisphere were selected using the "AND" operation. The fibers running through to the thalamus were removed by using the "NOT" operation.

## *FMajor*

For the first ROI, a coronal plane including only the left occipital lobe was selected at the most posterior edge of the parietooccipital sulcus. The "OR" operation in this coronal plane delineated all the fibers of FMajor. The second ROI was drawn in the same coronal plane on the right hemisphere using the "AND" operation such that all the fibers in the right occipital lobe was selected.

#### *FMinor*

For the first ROI, a coronal plane at the middle point between the anterior tip of the frontal lobe and the anterior edge of the genu of the CC was selected using the mid-sagittal plane. The "OR" operation was used to select all the fibers in the entire left hemisphere. The second ROI was drawn in the same coronal plane on the right hemisphere and all the fibers in the right hemisphere were selected using the "AND" operation.

#### **INTER-TRACT CORRELATION ANALYSIS**

Shapiro-Wilk normality test was performed with DTI-derived metrics of all the 10 WM tracts of 26 neonate brains and 28 child brains. Distributions of DTI-derived metrics for most of the tracts in the neonate group did not show significant difference (*p >* 0*.*05) from normality. However, distributions of DTI-derived metrics of most of the tracts in the child group differed significantly (*p <* 0*.*05) from normality. Hence, following the method in the literature (Wahl et al., 2010), non-parametric Spearman's rank correlation coefficient ρ was used to measure all correlations. Subsequently, a correlation matrix was constructed for each of the 4 DTI-derived metrics using pairwise correlation values between any two tracts. Symmetric correlation matrices were obtained with a value of unity along the diagonal of the correlation matrix. The diagonal element represents perfect correlation of the DTI-derived parameter of the tract with itself and the offdiagonal element represents the correlation of the DTI derived parameter of one tract with that of another tract. There were 10∗*(*10 − 1*)/*2 = 45 nontrivial independent correlation values in each correlation matrix.

#### **STATISTICAL ANALYSIS**

Two independent tests of correlation matrices were performed to evaluate if the correlation matrices were significantly different from identity and homogeneous matrix (Rencher, 2002; Wahl et al., 2010). The null hypothesis to test for identity was the

correlation matrix is an identity matrix, *-*0 = ⎛ ⎜ ⎝ 100 *. . . ... . . .* 0 ··· 1 ⎞ ⎟ ⎠ . The

size of the correlation matrix was 10 × 10 as we were testing for pairwise correlations of 10 independent tracts. To test the correlation matrices for homogeneity, the null hypothesis was that the correlation matrix was homogeneous and that the non-diagonal

elements of the matrix were equal, i.e.,: μ0 = ⎛ ⎜ ⎝ 1 ρ ρ *. . . ... . . .* ρ ··· 1 ⎞ ⎟ ⎠ . This

null hypothesized homogeneous correlation matrix was derived by following the procedure outlined in the literature (Rencher, 2002; Wahl et al., 2010). Bonferroni correction was conducted for both test of correlation matrix against identify matrix and test of correlation matrix against matrix with equal non-diagonal elements. Once the DTI-derived correlation matrices were found to be significantly different from identity and homogeneity within each group, correlation matrices of each DTI-derived metric were compared between the two age groups. Spearman's rank correlation coefficients were converted to *z* values by using Fisher's r-to-z-transform (Fisher, 1915). Note that *r* in Fisher's r-to-ztransform is the spearman's rank correlation coefficient ρ in this study. *Z*-statistics was then performed to identify the pair of tracts that showed significant change in correlation strength from birth to puberty.

## **HIERARCHICAL CLUSTERING ANALYSIS**

Hierarchical clustering methods were used to characterize the patterns of inter-tract correlation in each matrix among the groups. We used 1-ρ, where ρ is the Spearman's rank correlation coefficient, as a measure of distance or dissimilarity between the WM tracts for the purpose of clustering. Hierarchical clustering was performed using hclust function in R version 2.15.2. Depending on the correlation coefficient, different WM tracts were successively grouped into larger groups and the results were visualized as a dendrogram. WM tracts that had stronger correlation among themselves were fused together and were linked together. To characterize the uncertainty of the linkage among WM tracts and reduce the standard error in the percentage of confidence level for each cluster, a multi-scale bootstrap with 1000 repetitions of the analysis was performed using pvclust function in R (Suzuki and Shimodaira, 2006). The multi-scale bootstrap analysis yields an approximately unbiased *p*-value of each linkage in hierarchical clustering and has been applied in various other studies (Shimodaira, 2002, 2004; Wahl et al., 2010). The threshold for determining statistical significance for the grouping of tracts was set at an unbiased *p*-value of 0.05 or 95% confidence interval.

### **RESULTS**

## **CHANGES OF WHITE MATTER MICROSTRUCTURE FROM NEONATES TO CHILDREN AROUND PUBERTY**

**Figure 1** shows the three-dimensional (3D) visualization of the traced 10 major WM tracts for a typical neonate and a typical child around puberty. These 10 major tracts, namely CST\_L, CST\_R, IFO\_L, IFO\_R, CGC\_L, CGC\_R, CGH\_L, CGH\_R, FMajor and FMinor, cover projection, limbic, association and commissural tract groups involved in distinct brain functions. The microstructural changes from neonates to children at puberty and measured with FA, MD, AxD, and RD are shown in **Figure 2**. For all 10 major WM tracts, FA values are higher in the child group than those of the neonate group while MD, AxD, and RD values of the child group are less. FA of CST\_L and CST\_R of both age groups are highest among all tracts, followed by FMajor and FMinor. RD of CST\_L and CST\_R of both age groups are

**FIGURE 1 | 3-D visualization of the traced WM tracts overlaid on mid-sagittal slice of the FA image of a typical neonate (A) and a typical child around puberty (B).** Different colors represent different tracts traced for both subjects. CGC\_L/R, CGH\_L/R, CST\_L/R, FMajor/FMinor and IFO\_L/R are painted by red, orange, green, blue, and yellow color, respectively.

lowest among all tracts, indicating better myelination of CST compared to all other tracts. MD and AxD of FMajor and FMinor are highest among all tracts for child group. The values of these DTI-derived metrics of all tracts for the two age groups are shown in Supplemental Table 2.

## **ENHANCED INTER-TRACT CORRELATION FROM NEONATES TO CHILDREN AROUND PUBERTY**

**Figure 3** shows the scatterplot of the FA (**Figure 3A**), RD (**Figure 3B**), AxD (**Figure 3C**), and MD (**Figure 3D**) values for 2 pairs of homologous tracts and 2 pairs of non-homologous tracts, CGC\_L vs. CGC\_R, CST\_L vs. CST\_R, FMajor vs. FMinor, and CST-R vs. CGC-R. The 2 homologous tracts and FMajor/FMinor represent tract pairs of the three tract groups. Significant correlations (*p <* 0*.*05) can be observed for these 4 pairs of tracts in both age groups.

Significant differences (*p <* 0*.*05, Bonferroni-corrected) were observed for tests of all inter-tract correlation matrices of both groups against identity matrix or matrix with equal non-diagonal entries. The inter-tract correlation matrices of neonate and child group and the differences of these correlation matrices for FA, RD, AxD and MD are shown in **Figures 4A–D**, respectively. General stronger inter-tract correlations can be observed in the child group compared to those in the neonate group for all DTI-derived metrics, represented by the warmer colors in correlation matrices of child group. The statistics with *z*-scores (right panel of **Figure 4**) show that 64.4% (29/45), 84.4% (38/45), 73.3% (33/45), and 73.3% (33/45) of independent entries in the correlation matrix of child group are significantly higher than the corresponding entries of neonate group for FA, RD, AxD and MD, respectively. The correlation matrix from RD shows highest percentage changes (84.4%) among correlation matrices from all DTI-derived metrics, indicating more widespread enhanced inter-tract correlations with RD measurements. Note that the denominator 45 above indicates the number of all independent entries in the correlation matrix. The absolute values of correlation coefficients for AxD and MD, represented by the warmer colors in **Figures 4C,D**, are higher in both neonate and child group than those for FA (**Figure 4A**) or RD (**Figure 4B**). Much smaller percentages of independent entries of the correlation matrices are associated with the situation where correlation coefficients are significantly higher in the neonates than the children. Specifically, these percentages are 8.9% (4/45, namely CGC\_R vs. FMinor, CGC\_R vs. FMajor, CGH\_L vs. IFO\_L and IFO\_L vs. IFO\_R), 6.7% (3/45, namely CGC\_R vs. FMinor, CGC\_R vs. FMajor and IFO\_L vs. FMajor), 4.4% (2/45, namely CGC\_R vs. FMinor and CGC\_R vs. FMajor) and 4.4% (2/45, namely CGC\_R vs. IFO\_R and CGC\_R vs. FMajor) for FA, RD, AxD and MD, respectively. The inter-tract correlation coefficient values based on FA, RD, AxD and MD measurements for both age groups are shown in Supplemental Table 3.

**from 2 homologous and 2 non-homologous tracts.** The 2 homologous pairs and FMajor/FMinor represent tract pairs of three tract groups for both age groups. Top panel shows the inter-tract scatter plots for the neonates and plots represents the data from an individual subject in that group. ρ is the Spearman's rank correlation coefficient of the tract pair while *p*-value shows the statistical confidence of the inter-tract correlation strength.

**from tract-level FA (A), RD (B), AxD (C) and MD (D) measurements of both age groups.** The left and middle panels show the inter-tract correlation matrices for neonates and children around puberty, respectively. The right panel shows the *z*-scores of the changes in correlation strength between the two age groups. In the *z*-score plots, the entries with significant (*p <* 0*.*05) changes in inter-tract correlation strengths are shown in red color while entries with non-significant (*p >* 0*.*05) change in inter-tract correlation strengths are shown in green color. Color bar encoding the correlation strengths in the left and middle columns is also shown.

### **RESHUFFLED INTER-TRACT CORRELATION PATTERNS FROM NEONATES TO CHILDREN AROUND PUBERTY**

The inter-tract correlations based on each of DTI-derived metrics are reshuffled from neonates to children at puberty. These reshuffled inter-tract correlation patterns can be appreciated from dendrograms based on FA, RD, AxD, and MD measurements in **Figures 5A–D**, respectively. In general, the pairing of homologous tracts is more prominent for child group compared to that of the neonate group based on measurements of all DTI-derived metrics. The reshuffling leading to a more organized pairing among WM tracts is most prominent with the dendrograms based on tract-level RD measurements (**Figure 5B**). For inter-tract correlation based on FA measurement (**Figure 5A**), the tract pairs of CST\_L/R becomes clear in the dendrogram of child group, while this pair is not as apparent in neonate group. The IFO\_L/R pair is prominent for dendrograms of both neonate and child group based on FA measurements (**Figure 5A**). For dendrograms based on RD measurements (**Figure 5B**), it is clear that all homologous tracts and FMajor/FMinor get paired for child group while the correlation patterns are more random for neonate group. For dendrograms based on AxD (**Figure 5C**) or MD (**Figure 5D**) measurements, stronger and more clusters of the homologous tracts can be observed. The clustering pattern obtained from MD (**Figure 5D**) is similar to that obtained from RD (**Figure 5B**). However, the homologous tracts of child group are not well paired in MD-based dendrogram (**Figure 5D**), compared to those in RD-based dendrogram (**Figure 5B**).

The ranks of the clustering are shown in bold. The confidence intervals for the clustering (the percentage values) are shown in italics.

#### **DISCUSSION**

In this study, dynamics of inter-tract correlations from birth to onset of adolescence was investigated with DTI-based tract-level microstructural measurements from 10 major WM tracts, CST\_L, CST\_R, IFO\_L, IFO\_R, CGC\_L, CGC\_R, CGH\_L, CGH\_R, FMajor, and FMinor. Higher WM tract integrity, reflected by higher FA, lower MD, AxD, and RD, were found for the corresponding tracts from neonates to children around puberty. It is clear that even at birth, nearly all major WM tracts demonstrate similar morphology as those in children around puberty. Significant correlations of homologous tracts are shown for both neonate and child groups. The comparisons of the inter-tract correlation matrices between the neonate and child group indicated that stronger inter-tract correlations are established during development. Using data-driven hierarchical clustering algorithm with no *a priori* information, we were able to reveal that the linkage patterns of the major tracts differ between the two age groups. Specifically, homologous tracts involved in similar brain functions tend to cluster together for children around puberty especially with tract-level RD measurements. Such clustering patterns of homologous tracts become more prominent from birth to puberty. These changes of inter-tract correlations between neonates and children around puberty suggest inhomogeneous but organized axonal development which causes the reshuffled inter-tract correlation pattern while keeping homologous tracts tightly correlated. To the best of our knowledge, this is the first study investigating dynamics of inter-tract correlations with DTIbased microstructural measurements during early human brain development.

#### **HETEROGENEOUS WM GROWTH**

The WM development is heterogeneous among different tract groups, but more homogeneous among homologous tracts. The heterogeneity among different tract groups includes heterogeneous tract-level measurements of DTI-derived metrics at each time point and heterogeneous changes of these DTI-based tractlevel measurements from birth to puberty. The general FA increase and general MD, AxD, and RD decrease for all WM tracts in early brain development shown in **Figure 2** are consistent with the previous findings (e.g., Barnea-Goraly et al., 2005; Snook et al., 2005; Eluvathingal et al., 2007; Dubois et al., 2008; Lebel et al., 2008; Gao et al., 2009; Giorgio et al., 2010; Schmithorst and Yuan, 2010; Tamnes et al., 2010; Westlye et al., 2010). The heterogeneity among different tract groups is most prominent with FA measurements for children around puberty (**Figure 2**). From **Figure 2A**, it is clear that FA of CST-L/R (project tract group) is highest in children around puberty, followed by FMajor/FMinor (commissural tract group), IFO-L/R (association tract group) and CGC/H\_L/R (limbic tract group). This order of tract group FA measurements is preserved back at birth, but the differences of FA among the tract groups are smaller (**Figure 2A**) for neonatal brains. The heterogeneous changes of tract-level DTI metrics in **Figure 2** may cause the reshuffling of inter-tract correlation patterns. Quite similar FA values (**Figure 2A**) as well as MD, AxD, and RD values (**Figures 2B,D**) of the homologous tracts can be found for both age groups.

## **STRENGTHENED AND RESHUFFLED INTER-TRACT CORRELATION DURING DEVELOPMENT**

At birth, significant correlations can be observed based on microstructural measurements of homologous tracts, as shown in upper panels of **Figures 3A,D**. From **Figure 4**, the general intertract correlations from birth to puberty are clearly stronger for all four DTI-derived microstructural measurements. With FA measurement as an example, each entry of correlation matrix in **Figure 4A** can be expanded to the correlation scatterplots such as those shown in panels of **Figure 3A**. Our results also indicated that significantly increased inter-tract correlations are more widespread with correlation coefficients based on RD measurements. For both age groups, correlation coefficients are highest based on AxD measurements (**Figure 4C**) and lowest based on FA measurements (FA). This pattern is especially clear for children around puberty. The overall higher inter-tract correlation coefficients based on AxD measurements suggest that axonal integrity is coherent among the WM tracts within an individual child's brain but varies among different children. To help understand this finding, we could assume a situation when each of the children around puberty has identical AxD for all 10 major WM tracts and the AxD values vary among these children. Under such situation, all inter-tract correlation coefficients based on AxD would be the perfect value 1. The relatively low inter-tract correlation coefficients based on FA measurements may be caused by heterogeneity of FA values among the WM tracts within each individual subject's brain. Larger variability of FA values among WM tracts can also be observed in **Figure 2A**.

From **Figure 5A**, the FA-based dendrograms show that several homologous tracts are clustered together even for neonates. The left and right IFO are clustered with rank 1 for dendrograms based on tract-level FA of neonatal brains. IFO is thought to play a role in integrating the information from auditory and visual cortices to the prefrontal cortex (e.g., Martino et al., 2010). Resting-state fMRI studies have shown consistent pattern of activation in auditory and visual networks in neonates (Fransson et al., 2007; Doria et al., 2010). It is noteworthy that in the restingstate fMRI studies, the connectivity was also identified through correlation of brain-oxygen-level-dependent (BOLD) signal fluctuations in the homologous brain regions. The strong correlations of BOLD signal time courses in visual and auditory networks between left and right hemisphere in neonatal brains may be related to tight cluster of IFO-L and IFO-R involved in these brain functions. Both neonatal and child dendrograms based on FA measurements (**Figure 5A**) have two more separate clusters of projection tracts (CST) and limbic tracts (CGC). All these clusters can also be found in adult brains (Wahl et al., 2010). **Figure 5B** demonstrates one of the most compelling findings with dendrograms based on RD measurements. Although two pairs of homologous WM tracts are clustered together with rank 1 (IFO\_L/R) and rank 3 (CGC\_L/R), the other two pairs of homologous tracts and FMajor/FMinor are spread all over the dendrogram for neonatal brains (**Figure 5B**). It is striking that all 4 pairs of homologous tracts and FMajor/FMinor are tightly clustered together for children around puberty (**Figure 5B**). RD is closely related to myelination of WM tracts (Song et al., 2002). We should be careful to associate the RD values with myelination due to crossing-fiber and pathological situations (Wheeler-Kingshott and Cercignani, 2009). Nevertheless, during normal brain development with no implication of pathology and with the assumption that tract-level RD measurement for major WM tracts is much less affected by crossing fiber compared to voxelwise RD measurement, tract-level RD is still considered as an important index reflecting the degree of myelination (e.g., Snook et al., 2005; Eluvathingal et al., 2007; Gao et al., 2009; Tamnes et al., 2010; Westlye et al., 2010). The results in **Figure 5B** suggest that organized myelination from birth to puberty plays an important role to reshuffle the inter-tract correlations and result in clustered homologous tracts and clustered functionally similar tracts. The dendrogram patterns based on AxD and MD measurements are different than those based on FA or RD measurements. Unlike RD-based dendrogram (**Figure 5B**), neither of the dendrograms based on AxD or MD measurements for children around puberty shows well organized clusters of all 4 pairs of homologous tracts and FMajor/FMinor (**Figures 5C,D**). Previous studies (Mukherjee et al., 2002; Snook et al., 2005; Gao et al., 2009) found that the measurements of RD decrease dramatically with relatively little changes in AxD of the major WM tracts during brain development. Our results demonstrated in **Figure 2** also indicate smaller and more homogeneous changes of AxD from birth to puberty compared to those of RD. These AxD change features may explain why relatively disorganized AxD-based dendrograms remain for children around puberty. From Equation (1), MD is the linear combination of AxD and RD. The relatively disorganized dendrograms based on MD measurement for child group could also be originated from smaller and more homogeneous changes of AxD measurements during development. Nevertheless, with the brain development, there are still trends for homologous and functionally similar tracts to cluster together in dendrograms based on AxD and MD in **Figures 5C,D**.

## **POSSIBLE MECHANISMS OF THE INTER-TRACT CORRELATION CHANGES FROM BIRTH TO ONSET OF ADOLESCENCE**

This study on dynamics of inter-tract correlations from birth to onset of adolescence provides unique insight on the wellorganized cerebral WM development. The development of human cerebral WM tracts is characterized with enhanced myelination and axonal integrity. From the results in this study, inhomogeneous RD decreases among the tracts take place during development. More (84.4% of all independent correlation coefficients) inter-tract correlations become stronger with RD measurements, compared to those with any other DTI metric measurements. In addition, the dendrograms based on RD (**Figure 5B**) demonstrate that all 4 pairs of homologous tracts and FMajor/FMinor are tightly clustered at puberty while only 2 pairs of homologous tracts are clustered at birth. The cerebral WM at birth is likely to be in a relatively random and disorganized status. Due to close relationship of RD with the myelination of WM tract, these results suggest that inhomogeneous enhancement of myelination rather than strengthening of axonal integrity plays a key role in reshaping the WM configuration during development. Both genes and experiences could also play a role to adjust WM microstructures so that the homologous WM tracts reach to a coherent status to meet the needs of certain brain functions by the time of onset of adolescence. Although it is not known by what mechanism inhomogeneous myelination is modulated during WM maturation, our results suggest that the myelination process is precisely controlled so that all 4 pairs of homologous tracts and FMajor/FMinor are clustered together around puberty (**Figure 5B**).

#### **LIMITATIONS OF THIS STUDY AND FUTURE DIRECTIONS**

There are several issues which may affect the results in this study. Five pairs of major WM tracts were chosen for this study due to the fact that only these five pairs of major WM tracts could be reproducibly traced with neonate DTI. The numbers of participated subjects, 26 for neonates and 28 for children, just exceeds 25 which is needed for correlation analysis of five pairs of WM tracts. Higher sample numbers could increase the confidence level in analysis of hierarchical clustering. Corrections for multiple comparisons were only performed on testing the correlation matrices against identity matrix and matrix with equal non-diagonal elements. No correction was performed on the hierarchical clustering results due to lack of any known methods to perform such a correction on dendrograms. The accuracy of measuring tract-level DTI-derived metrics plays a key role in inter-tract correlations. This accuracy is affected by three major factors. They are the crossing-fiber factor, SNR of the data and partial volume effects. Both the tractography and the DTI-derived measurements are biased at crossing-fiber regions (Wheeler-Kingshott and Cercignani, 2009). With single tensor model and tractography method of fiber assignment by continuous tracking (FACT) (Mori et al., 1999), it is apparent that the tracing method adopted in this study cannot resolve the crossingfiber issue, resulting in imperfect binary mask of the traced tracts for tract-level measurements of DTI metrics. Nevertheless, the tracing protocol (Wakana et al., 2007) based on FACT tractography captures the core of the major WM tracts and is still widely used in the field. We conducted two repetitions of diffusion MRI and the SNR is sufficient for data acquisition with 3T magnet. With diffusion imaging resolution 2 × 2 × 2 mm3 for neonates and 2 × 2 × 2*.*2 mm3 for children around puberty, the partial volume effects are inevitable. However, it seems the effects of imperfect WM fiber tracing offsets the partial volume effects for obtaining accurate tract-level DTI metrics, in that the tracing algorithm adopted in this study cannot trace the small branches of the fibers where the problem of partial volume effects is most prominent. Ranks of the DTI measurements, instead of measured metrics themselves, were used for correlation pattern analysis. Despite possibly different levels of biases of the DTI metric measurements caused by partial volume effects due to different head sizes between neonates and children, the similar shifts of measurements in the same age group could have minimum effects on the rank of metric measurements and therefore minimum effects on the Spearman's rank correlation patterns. Although same types of scanners were used in this study, systematic differences of the scanners may affect the DTI analysis results. To make sure that the effects of systematic differences caused by two different scanners are minimal, a healthy young subject ("*in vivo* human phantom") was scanned in both scanners used in this study and with the same DTI sequence. The quantitative comparisons were also conducted. Quantitative DTI 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). With the same type of scanners used in this study and rigorous quality control of both scanners, the effects of scanner differences on the presented results in this study are thus negligible. The group of children around puberty included both pre-puberty and post-puberty subjects. The age difference between the two groups, newborns and children around puberty, is much larger than the age difference within the group of around puberty. Therefore, we hypothesized that the intra-group WM developmental heterogeneity for the children group exists, but is not big enough to affect the intergroup results presented in this study. This has been tested and proved by **Supplemental Figure 1**. In **Supplemental Figure 1**, we separated the children around puberty into two subgroups, prepuberty (9.5–12 years) and post-puberty (12–15 years). With reduced sample number for each subgroup, we could conduct the correlation analysis for 6 tracts (3 pairs of homologous tracts). It is shown in **Supplemental Figure 1** that homologous tracts are still clustered together in both subgroups, like the cluster patterns shown in the right column of **Figure 5**. Moreover, the linkage patterns are very similar between the two subgroups, with slight dendrogram rank change.

In the future, several improvements can be made to address the issues affecting accuracy of measuring tract-level DTI metrics. The tracking methods (e.g., Tournier et al., 2004; Behrens et al., 2007) which are capable of resolving the crossing-fiber issue can be adopted. Although there is no general consensus on which parameters can replace these four DTI-derived metrics and better characterize the WM microstructure, there have been a few metrics such as general fractional anisotropy (GFA) (Tuch, 2004; Fritzsche et al., 2010; Zhan et al., 2010), generalized anisotropy (GA) (Ozarslan et al., 2005), mode of anisotropy (Ennis and Kindlmann, 2006; Douaud et al., 2011) and fractions (Hosey et al., 2008; Jbabdi et al., 2010) which are less sensitive to crossing-fiber problem. In addition, tract-based spatial statistics (TBSS) from FSL (http://www*.*fmrib*.*ox*.*ac*.*uk/fsl) (Smith et al., 2006) can be used to alleviate the partial volume effects, as shown in the study of FA correlations in adults (Li et al., 2012).

## **CONCLUSION**

In conclusion, inter-tract correlation changes during development from birth to onset of adolescence were investigated with tract-level FA, RD, AxD, and MD measurements. Stronger and enhanced microstructural inter-tract correlations were found during development. The linkage patterns of the major tracts also differ with the dendrograms of two age groups due to brain development. These changes of microstructural correlations from birth to puberty suggest inhomogeneous but organized myelination processes which cause the reshuffled inter-tract correlation pattern and make homologous tracts tightly clustered. Especially RD-based dendrograms reveal that all 4 pairs of homologous tracts and FMajor/FMinor investigated in this study are tightly clustered for children around puberty while only 2 out of these 5 pairs of tracts are clustered at birth, indicating important role of myelination to reshape the WM configuration. It opens a new window to study WM tract development and can be potentially used to investigate atypical brain development due to neurological or psychiatric disorders.

#### **ACKNOWLEDGMENTS**

This study is sponsored by NIH (Grant Nos. EB009545, MH092535, and MH092535-S1, Hao Huang), the Natural Science Foundation of China (Grant Nos. 31221003 and 81030028, Yong He), the National Science Fund for Distinguished Young Scholars of China (Grant No. 81225012, Yong He) and the 111 Project (Grant No. B07008, Qi Dong).

#### **SUPPLEMENTARY MATERIAL**

The Supplementary Material for this article can be found online at: http://www.frontiersin.org/journal/10.3389/fnhum.2013. 00721/abstract

#### **Supplementary Figure 1 | Dendrograms depicting the hierarchical clustering pattern obtained from tract-level FA (A), RD (B), AxD (C) and MD (D) measurements for two subgroups of children around puberty.**

Dendrograms on the left column are for the subgroup of children younger than 12-years-old (9.5–12 years); those on the right column are for the subgroup of children older than 12-year-old (12–15 years).

 *Magn. Reson. Imaging* 26, 236–245. doi: 10.1016/j.mri.2007.07.002

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 18 June 2013; accepted: 11 October 2013; published online: 29 October 2013. Citation: Mishra V, Cheng H, Gong G, He Y, Dong Q and Huang H (2013) Differences of inter-tract correlations between neonates and children around puberty: a study based on microstructural measurements with DTI. Front. Hum. Neurosci. 7:721. doi: 10.3389/fnhum.2013.00721*

*This article was submitted to the journal Frontiers in Human Neuroscience.*

*Copyright © 2013 Mishra, Cheng, Gong, He, Dong and Huang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Graph theoretical analysis of developmental patterns of the white matter network

#### *Zhang Chen1, Min Liu1, Donald W. Gross 2 and Christian Beaulieu1 \**

*1 Department of Biomedical Engineering, Faculty of Medicine and Dentistry, University of Alberta, Edmonton, AB, Canada 2 Division of Neurology, Department of Medicine, Faculty of Medicine and Dentistry, University of Alberta, Edmonton, AB, Canada*

#### *Edited by:*

*Yong He, Beijing Normal University, China*

#### *Reviewed by:*

*Hao Huang, University of Texas Southwetern Medical Center, USA Pew-Thian Yap, University of North Carolina, USA*

#### *\*Correspondence:*

*Christian Beaulieu, Department of Biomedical Engineering, 1098 Research Transition Facility, University of Alberta, 8308-114 Street, Edmonton, AB T6G 2V2, Canada e-mail: christian.beaulieu@ ualberta.ca*

Understanding the development of human brain organization is critical for gaining insight into how the enhancement of cognitive processes is related to the fine-tuning of the brain network. However, the developmental trajectory of the large-scale white matter (WM) network is not fully understood. Here, using graph theory, we examine developmental changes in the organization of WM networks in 180 typically-developing participants. WM networks were constructed using whole brain tractography and 78 cortical regions of interest were extracted from each participant. The subjects were first divided into 5 equal sample size (*n* = 36) groups (early childhood: 6.0–9.7 years; late childhood: 9.8–12.7 years; adolescence: 12.9–17.5 years; young adult: 17.6–21.8 years; adult: 21.9–29.6 years). Most prominent changes in the topological properties of developing brain networks occur at late childhood and adolescence. During late childhood period, the structural brain network showed significant increase in the global efficiency but decrease in modularity, suggesting a shift of topological organization toward a more randomized configuration. However, while preserving most topological features, there was a significant increase in the local efficiency at adolescence, suggesting the dynamic process of rewiring and rebalancing brain connections at different growth stages. In addition, several pivotal hubs were identified that are vital for the global coordination of information flow over the whole brain network across all age groups. Significant increases of nodal efficiency were present in several regions such as precuneus at late childhood. Finally, a stable and functionally/anatomically related modular organization was identified throughout the development of the WM network. This study used network analysis to elucidate the topological changes in brain maturation, paving the way for developing novel methods for analyzing disrupted brain connectivity in neurodevelopmental disorders.

**Keywords: graph theory, neurodevelopment, anatomical connectivity, modular networks, small world network**

## **INTRODUCTION**

Neuroimaging studies have demonstrated widespread and regionally specific structural and functional brain changes during development from infancy to adulthood. Structural magnetic resonance imaging (MRI) studies have reported age-related changes in brain volumes (Giedd et al., 1999; Good et al., 2001), areas (Thompson et al., 2000), cortical thickness (Sowell et al., 2004; Shaw et al., 2008), and regional gray matter (GM) and white matter (WM) density (Paus et al., 1999; Gogtay et al., 2004). The developmental changes in GM and WM on gross scale MRI may reflect synaptic pruning and myelination occurring at the neuronal level (Gogtay et al., 2004). Functional neuroimaging studies have demonstrated increased connectivity among distant regions and decreased connectivity among neighboring regions in brain maturation which suggests a mechanism of segregation of local regions and integration of distant regions into disparate subnetworks for the developing brain (Fair et al., 2008, 2009; Vogel et al., 2010). Diffusion tensor imaging (DTI) studies of WM have shown age-related increases in fractional anisotropy (FA) and decrease in overall diffusion with development [many studies

but some include (Snook et al., 2005; Lebel et al., 2008; Tamnes et al., 2010)], including into young adulthood (Giorgio et al., 2008; Lebel and Beaulieu, 2011).

The recent advent of modern network analysis based on graph theory (Strogatz, 2001), has enabled the investigation of the largescale topological organization of various structural and functional brain networks such as the small-world property, network efficiency and modularity (He et al., 2007; Bullmore and Sporns, 2009; He and Evans, 2010). The network metrics have also proven useful in modeling the large-scale functional and structural organization of the developing brain. Several functional brain network studies have reported age-related increases in the small-worldness (Wu et al., 2013) and a progression from local to distributed organization (Fair et al., 2009) in brain development. The analysis of the structural brain network constructed from regional cortical thickness correlations has revealed a non-linear developmental pattern in network metrics and that most topological changes happen at the late childhood stage (Khundrakpam et al., 2013).

Recently, there has been an increasing interest in the study of how graph metrics of the anatomical brain network change during development. Using DTI, Yap et al. (2011) examined WM networks of 39 healthy pediatric subjects with longitudinal data collected at average ages of 2 weeks, 1 year, and 2 years and demonstrated that the small-world architecture exists at birth with efficiency that increases in later stages of development. Two recent brain connectivity studies of WM maturation pattern using diffusion MRI tractography demonstrated linear and non-linear patterns of increasing structural efficiency with age between ages 2 and 18 years in 30 patients scanned clinically and otherwise deemed normal post-MRI (Hagmann et al., 2010) and between ages 12 and 30 years in 439 healthy subjects (Dennis et al., 2013). However, those studies were limited by different constraints such as a binarized brain network, limited sample size, or restricted age range (early adolescence to early adulthood), thus the developmental trajectory of the WM network from early childhood to adulthood remains unclear.

Therefore, the main goal of this study was to map the developmental changes of the structural brain network based on WM connectivity in 180 typically-developing subjects from 6 to 30 years of age. We hypothesized (i) non-linear age-related developmental trajectories of network metrics as most changes would be expected to happen at late childhood and adolescent stages, and (ii) altered modular organization in different age groups that reflects a process of fine-tuning in structural brain development.

## **MATERIALS AND METHODS**

#### **SUBJECTS**

This study included 180 healthy right-handed subjects aged from 6 to 30 years. Health was verified by asking participants a series of questions to ensure there was no history of neurological or psychiatric disease or brain injury. All subjects gave informed consent; child assent and parent/guardian consent was obtained for volunteers under 18 years. The subjects were divided into 5 age groups with equal numbers of subjects and demographics of all groups are shown in **Table 1**.

#### **IMAGE ACQUISITION**

All data were acquired on a 1.5-T Siemens Sonata MRI scanner. Standard DTI was acquired using a dual spin-echo, single shot echo-planar imaging sequence with the following parameters: 40 3-mm-thick slices with no inter-slice gap, *TR* = 6400 ms, *TE* = 88 ms, 6 non-collinear diffusion sensitizing gradient directions

with *b* = 1000 s/mm2, 8 averages, field-of-view 220 × 220 mm2, matrix of 96 × 128 zero-filled to 256 × 256, and scan time of 6:06 min. T1-weighted images were also acquired using MPRAGE with *TE* = 4*.*38 ms, *TR* = 1870 ms, *TI* = 1100 ms, and scan time of 4:29 min.

#### **STRUCTURAL BRAIN NETWORK CONSTRUCTION**

Image preprocessing steps including motion and eddy current corrections were performed using FSL 5.0 for all DTI images (http://fsl*.*fmrib*.*ox*.*ac*.*uk/fsl/fslwiki). The T1-weighted (MPRAGE) image of each subject was first linearly coregistered (**Figures 1A,B**) to its corresponding b0 image. Each transformed T1 image was then non-linearly registered to a pre-segmented and validated volumetric template, the automated anatomical labeling (AAL) atlas (Tzourio-Mazoyer et al., 2002) as shown in **Figures 1B,C**. This parcellation divided the cortical surface into 78 regions (39 per hemisphere). See **Table 2** for the name of the regions and their corresponding abbreviations. The resulting inverse deformation map (*T*−1) for each subject was then applied to warp the AAL template to the DTI native space of each subject using nearest neighbor interpolation method (**Figures 1E,F**) as each AAL region was defined as a brain network node. Whole brain WM tractography was performed using a brute-force streamline-tracking method (Basser et al., 2000) with a FA threshold of 0.2 and primary eigenvector turning angle of 45 degrees (**Figures 1A,D**). To limit false positive connections, two cortical regions were deemed connected if at least 10 connecting fibers with two end points were located between them; the same threshold was also applied in a recent brain network study (van den Heuvel et al., 2012). The effects of different nodeconnecting fiber number (FN) thresholds ranging from 3 to 10 were determined for our network analysis. We quantified the weight of each valid connection between two cortical regions (*i* and *j*) as the product of the connecting FN and mean FA of the connecting fiber, normalized by dividing the average volume of the two connecting regions to counteract the bias where larger cortical regions inherently project/receive more "virtual" fibers (*wij* = FN∗FA/Volume). Several diffusion brain network studies have applied this weighting function (Lo et al., 2010; Brown et al., 2011). As a result, the structural brain network of each participant was represented by a symmetric 78 × 78 matrix (**Figure 1G**).

To examine the small-worldness and modular organization of the WM networks for all different age groups, one weighted backbone network for each age group was generated to capture the underlying anatomical connectivity patterns using a previously published method by our group (Gong et al., 2009). In summary, to identify the highly consistent cortical connections, a non-parametric one-tailed sign test was applied. For each pair of cortical regions, the sign test was performed with the null hypothesis that there is no existing connection. The Bonferroni method was applied to correct for multiple comparisons at *P <* 0*.*05. The use of this conservative statistical criterion generated a symmetric weighted matrix as each edge weight was calculated as the mean of all existing edges in all subjects that captured underlying anatomical connectivity patterns in the human cerebral cortex (Gong et al., 2009).

#### **NETWORK ANALYSIS**

Several network topological properties were applied for the weighted anatomical brain network derived from each participant, including small-worldness, efficiency and modularity (Watts and Strogatz, 1998; Latora and Marchiori, 2001; Newman, 2006). The connection weights of all edges (*wij*) were normalized by the mean weight of the whole network to keep network cost at the same level for all subjects.

For a weighted network *G* with *N* nodes and *K* edges, the total strength *S* was defined as the mean of all edge weights in the network, *S(G)* = 1 *N N* -*M i* = *j* ∈*G wij* where *i* and *j* are two distinct nodes in graph *G*. The clustering coefficient (*CC*) of the weighted network *G* quantifies the likelihood whether the neighboring nodes of any network nodes are connected with each other (Onnela et al., 2005), which was defined as: *CC* = 1 *N* - *j,k* ∈*G (wij wjk wik)*1*/*3*/(k*∗ *i (ki* − 1*)/*2*)*, where *ki* is the number of

connected neighbors of node *i*. The weighted characteristic path length *L* of a network is the average minimum connectional weights that link any two nodes of the network. To avoid the issue of disconnected nodes, *L* was measured here by using a "harmonic mean" distance between any pair of nodes such as the reciprocal of the average of the reciprocals (Newman, 2003). A real network would be considered small world if it meets the following criteria: γ = *Creal p /Crand p >* 1 and λ = *Lreal p /Lrand p* ∼ 1 (Watts and Strogatz, 1998), where *Crand p* and *Lrand p* are the mean *CC* and *L* of matched random networks that preserve the same number of nodes, edges and degree distribution as the real network (Maslov and Sneppen, 2002). In this study, we generated 1000 matched random networks for each group network.

The global efficiency *Eglob* of a weighted network *G* is defined as *Eglob(G)* = 1 *N(N*−1*)* - *i* = *j*∈*G* 1 *wij* , where *wij* is the smallest connectional weight between node *i* and *j* and *N* is the number of nodes. It characterizes the efficiency of a system transporting information in parallel (Latora and Marchiori, 2003). The local efficiency *Eloc* of a weighted network *G* is defined as: *Eloc(G)* = 1 *N* - *i* ∈*G Eglob(Gi)*, where *Gi* denotes the subgraph composed of the nearest neighbors of node *i*. The local efficiency represents the fault tolerance level of the network in response to the removal of a node (Latora and Marchiori, 2003).

The regional global efficiency *Ereg* of a given node *i* is defined as: *Enodal(i)* = 1 *N*−1 - *i* = *j* ∈*G* 1 *wij* , as it measures the average smallest path weight between given node *i* and all other nodes in the network. The node *i* was considered as a hub if its regional global efficiency was at least one standard deviation (*SD*) greater than the mean nodal efficiency of the whole network.

The modularity *Q(p)*for a given partition *p* of a weighted brain structural network is defined as *Q(p)* = *N* -*M s* = 1 *ws W* − *Ws* 2*W* 2 , where *NM* is the number of modules, *W* is the total weight of the network, *ws* is the sum of the connectional weights between all nodes in module *s*, and *Ws* is the sum of the nodal weights in module *s*. The modularity index quantifies the difference between the weight of intra-module links of the actual network and that of a random network in which connections are weighted at random (Newman, 2004). The aim of this module identification process is to find a specific partition *p* which yields the largest network modularity, *Q(p)*. Here a modified greedy optimization algorithm (Clauset et al., 2004; Danon et al., 2006) is used to find the modules of the network. The advantage of this modularity optimization method is that it takes into account the heterogeneity of module size observed in real networks (Danon et al., 2006).

We also determined the participation coefficient (PC) for each cortical region in terms of their inter-modular connection density (Guimera and Amaral, 2005; Guimera and Nunes Amaral, 2005; Sales-Pardo et al., 2007). The PC, *Pi*, measures the intermodular connectivity of node *i*, and is defined as: *Pi* = 1 − *N* -*M s* = 1 *wis wi* 2 , where *NM* is the number of modules and *wis* is intermodular connectional weight between the node *i* and module

*s*. *wi* is the total weight of node *i* in the network. The PC of node *i* will be close to 0 if all weights are within its module.

The node *i* was considered as an inter-modular hub if its *PC* value was at least one *SD* greater than the mean PC of the whole network.

#### **STATISTICAL ANALYSIS**

Between-group differences analysis of all the global network metrics (*S*, *CC*, *L*, *Eglob*, *Eloc*, *Ereg , Q*) was performed between adjacent age groups using the General linear model (GLM) with age and gender included as covariates. The nodal properties (*Ereg* , *z, P*) were corrected by false discovery rate at *q* = 0*.*05 (Genovese et al., 2002; Zeisel et al., 2011).

### **RESULTS**

#### **AGE-RELATED CHANGES IN FIBER NUMBER AND NETWORK SPARSITY**

To examine the age effect on the tractography results, we mapped age-related changes in the FN and sparsity of white matter network as shown in **Figures 2A,B**. We found that age has an incremental effect on both the FN and sparsity, where both increase by a factor of ∼1.6 and ∼1.2, respectively, from age 6 to 30 years. These increases are presumably due to the known elevations of FA in WM with age. Given the fact that our network edge weighting function depends on FN and FA, it is expected that the connectivity strength of the network would also increase with age.

#### **SMALL-WORLD EFFICIENCY OF DEVELOPING WM NETWORKS**

To examine the small-worldness of the WM networks for all different age groups, using a previously published method by our group (Gong et al., 2009), one weighted backbone network for each age group was generated to capture the underlying anatomical connectivity patterns as shown in **Figure 3**. Compared with their corresponding 1000 random networks, all

five age groups showed strong small-worldness (σearly childhood = 3*.*54, σlate childhood = 3*.*19, σadolescence = 3*.*25, σyoung adult = 3*.*12, σadult = 3*.*19).

### **GLOBAL NETWORK PROPERTIES AND THEIR AGE-RELATED TRAJECTORIES**

Over all subjects in each age group, the total network weight, CC, Lp, modularity (Q), Eglob and Eloc was calculated for the WM network and the age-related trajectories are shown in **Figure 4**. The total network weight displayed significant increases in three of the four developing stages, whereas the other metrics such as L*p*, Q, Eglob, and Eloc demonstrated non-linear alteration patterns where most changes happened from young childhood to late childhood that then leveled off. Both L*p* and Q decreased significantly from young childhood to late childhood but stabilized at older ages. Global network efficiency increased significantly from young childhood to late childhood but also stabilized

later. Local network efficiency increased significantly from late childhood to adolescence and stabilized afterwards.

#### **REGIONAL EFFICIENCY OF THE DEVELOPING WM NETWORKS**

We found consistent hubs regions, measured here as the AAL areas with highest regional global efficiency, such as bilateral PCUN, SFGdor, and SFGmed, that are shared by all age groups as shown in **Figure 5**. Comparing the regional efficiency changes from group to group in these hubs, seven regions had increased nodal efficiency (*P <* 0*.*05, FDR corrected) from early childhood to late childhood and two regions from late childhood to adolescence (**Figure 6**). Most regional changes from early to late childhood are in the default-mode system, including bilateral PCUN and left DCG. Left STG and right INS were found to have increased efficiency from late childhood to adolescence.

**FIGURE 4 | Age-related changes in different network metrics for the developing WM network from early childhood to adult. (A)** Total network strength (S), **(B)** Clustering coefficient (CC), **(C)** Shortest path length (L*p*), **(D)** Modularity (Q), **(E)** Global efficiency (Eglob), and **(F)** Local efficiency (Eloc). Significant changes between any adjacent age groups are indicated by their *p* value. An increase with age is observed in S over the full age span. Eglob increases only between the two youngest age groups and Eloc only between late childhood and adolescence; in both cases, the efficiency values then stay elevated. Reductions are observed in L*p* and Q from early to late childhood that is then maintained low. There is no change in CC between any adjacent age groups. The + signs indicate outliers.

### **MODULAR ORGANIZATION AND CONNECTOR HUBS OF THE DEVELOPING WM NETWORKS**

The modular organization of the developing structural brain networks for the five different age groups is shown in **Table 3** and **Figure 7**. Six modules (1–6) were detected in all age groups indicating strong stability (Greicius et al., 2003) in the modularity of the developing brain network. Despite decreased modularity from early to late childhood, the modular structures of both groups were almost identical. Module 1 was mostly composed of bilateral orbitofrontal regions (ORBsup, ORBsupmed, ORBmid, REC) in early and late childhood that becomes more lateralized in adolescence. Right orbitofrontal regions become connected with right temporal and occipital regions that resembles the ventral visual system (Grill-Spector et al., 2008) and left orbitofrontal regions become part of lateral frontal system. Module 2 consists of mostly occipital regions (SOG, CAL, CUN) throughout the youngest age groups except at adulthood when the left occipital regions become part of left ventral visual system (Grill-Spector et al., 2008). Lateralized modules 3 (left hemisphere) and 4 (right hemisphere) consist of regions mostly across frontal, parietal and temporal lobes from each hemisphere from early to late childhood. However, module 4

indicated in **Table 2**.

SFGmed are consistent over all age groups.

**FIGURE 6 | Regions with significantly increased nodal efficiency from early childhood to late childhood and late childhood to adolescence.** Seven regions had increased nodal efficiency (*P <* 0*.*05, FDR corrected) from early childhood to late childhood and two regions from late childhood to adolescence. Most regional changes from early to late childhood are in the default-mode system, including bilateral PCUN and left DCG. Left STG and right INS have increased efficiency from late childhood to adolescence.

is pruned to a mainly frontal-parietal system from adolescence onwards and module 3 doesn't reach a similar outcome until adulthood. Modules 5 and 6 are two of the most consistent modules during development and include mostly bilateral frontal (SFG, MFG) and posterior parietal (PCUN, DCG, SMA) regions, respectively.

The distribution of inter-modular hubs based on PC of each region for different age groups was very consistent (**Figure 8**). They were mostly located within posterior cortex, including bilateral PCUN, SPL, and MOG. Large frontal hubs such as bilateral SFGdor appeared in late childhood and remained significant afterwards.

To ensure the change of modular organization didn't result from the different sparsity of the five age groups, an additional analysis was performed where we normalized all group networks sparsity to the lowest sparsity at 0.1142 (early childhood) as weaker connections were removed from the other networks with higher sparsity (late childhood, 0.1262, adolescence, 0.1312, young adult, 0.1289, adult, 0.1329) and re-examined the modularity of the networks. We found extremely consistent modularity (0.56 previous vs. 0.56 with normalized sparsity in late childhood, 0.56 vs. 0.56 in adolescence, 0.55 vs. 0.56 in young adults, and 0.58 vs. 0.56 in adults) and modular organization compared with our original networks. Thus, we could presume that changes in the backbone network and its modular organization were not due to different matrix densities.

## **DISCUSSION**

The present study utilized DTI tractography and network theory to characterize changes to the global structural WM network with age from early childhood to adulthood. Our main results are demonstrations of (1) a non-linear age effect on most network topological properties of brain WM network in development where most changes happen at late childhood stage (10–13 years), such as increased global network efficiency and decreased modularity, suggesting a shift of organization toward a more randomized configuration, (2) consistent hubs involving several major functional systems across all age groups and significant nodal changes only happening from early childhood to adolescence, (3) anatomically localized modules in the development of brain WM network, and (4) key connector hubs during development of the WM network.

First, using graph theoretical analysis, small-world network architecture was demonstrated in the WM networks of all age groups. During the last decade, graph theoretical analysis has been widely applied to both the functional (Stam, 2004; Bassett et al., 2006; Achard and Bullmore, 2007) and anatomical (He et al., 2007; Hagmann et al., 2008; Gong et al., 2009) brain networks and one common finding is the existence of "small-worldness" in all types of network, as defined by high CC and low characteristic path length. Recent structural brain network studies have also

#### **Table 3 | Cortical regions in each module of developing white matter network.**

*(Continued)*

#### **Table 3 | Continued**

revealed that small-world topology and modular organization are established during early brain development (*<*2 years) to support rapid synchronization and information transfer with minimal rewiring cost (Fan et al., 2011; Yap et al., 2011). Thus, our results are in agreement with previous findings that the WM network maintains small-world efficiency at all stages of development.

Total network weight shows increases in three of the four developing stages, although the greatest change is between the two youngest groups pre-adolescence. Our finding is consistent with a previous WM network study that reported a significant increase in the average node strength in a group of subjects aged from 18 months to 18 years where it was suggested that increased network weight indicates increased nodal strength and greater physiological efficacy, particularly of long pathways (Hagmann et al., 2010). A functional network study has also reported increased functional integration due to a decrease of average path length during the same period and suggested it was related to increased axonal diameter and myelin thickness of long association fiber tracts (Supekar et al., 2009). We also found significant age-related decreases in the shortest path length and modularity and increase in the global efficiency of the developing WM network from early childhood to late childhood indicating greater integration among distant brain regions and a shift of topological organization to a more randomized configuration. Previous WM network (Hagmann et al., 2010) and cortical thickness network (Khundrakpam et al., 2013) studies of brain development also demonstrated a similar pattern of network metrics evolution between age 2 and 18 years and between 5 and 18 years, respectively. However, the WM network study had applied a linear fit for all the network metrics vs. age even though network metrics such as efficiency and clustering seemed to have leveled off after late childhood in their paper (Hagmann et al., 2010). Using a similar approach to ours, (Khundrakpam et al., 2013) demonstrated a leveling off of various cortical thickness network metrics after the early adolescence stage.

Consistent global hub regions, indicated by higher regional efficiency, are observed across all age groups. Hub regions are predominately association cortices that receive convergent inputs from multiple cortical regions. Regions such as SFG and PCUN have been constantly identified as the hub regions in both structural (He et al., 2007; Gong et al., 2009) and functional brain networks (Achard and Bullmore, 2007). A recent structural brain network study also identified them as the hub regions from age 2 years suggesting that they are established at a very early age (Hagmann et al., 2010). We also found that the regions with the most age-related increases in efficiency are in the default-mode system, including bilateral PCUN and left DCG. A functional brain network study has reported a less well-developed default mode network connectivity in early childhood compared with adults, especially within posterior regions such as PCUN (Fair et al., 2008). However, evidence from structural covariance network analysis has demonstrated significant pruning in the default mode system from early childhood to late childhood (Zielinski et al., 2010). Thus, we could speculate that nodal efficiency of default mode regions might plateau by late childhood.

In this study, a stable and functionally/anatomically related modular organization was demonstrated in the developing WM network. Six modules comprising regions with known functions or connections were identified in the developing WM network. Modules 1, 2, 5 and 6 were mostly composed of orbitofrontal, occipital, frontal, and posterior parietal regions that could correspond to sensory integration, visual, executive function, and default mode network, respectively, (Duncan and Owen, 2000; Raichle et al., 2001; Kringelbach, 2005). Modular network analysis has provided rich quantitative insights into the organization of complex brain networks.  Six **modules (1–6)** were detected in all age groups and are represented by red, green, purple, yellow, pink, and blue colors. See **Table 3** for a detailed list of modular regions. Module 1 was mostly composed of bilateral orbitofrontal regions (ORBsup, ORBsupmed, ORBmid, REC) in early and late childhood and becomes more lateralized from adolescence onwards. Module 2 consists of mostly occipital regions (SOG, CAL, CUN) bilaterally. Lateralized modules 3 (left hemisphere) and 4 (right hemisphere) consist of regions mostly across frontal, parietal and temporal lobes within each hemisphere. Modules 5 and 6 are two of the most consistent modules during development that include mostly bilateral frontal (SFG, MFG) and posterior parietal (PCUN, DCG, SMA) regions, respectively.

Zhou et al., 2006). Previous neuroimaging studies have also demonstrated anatomically- and functionally-related modules in the human brain structural network using diffusion spectrum imaging (Hagmann et al., 2008) and the functional network using resting-state functional MRI (Salvador et al., 2005; Ferrarini et al., 2009; He et al., 2009; Meunier et al., 2009; Valencia et al., 2009). Also, network modules identified by cortical thickness network analysis are comprised of brain regions known to subserve distinct brain functions such as executive function, vision, and default mode network (Chen et al., 2008, 2011).  Two lateralized modules (3 and 4) that correspond to the frontal-parietal network were also observed in the developing WM network. The adult human brain exhibits distinct hemispheric asymmetries in both structure and function. These asymmetries are thought to originate from evolutionary, developmental, hereditary, experiential, and pathological factors (Toga and Thompson, 2003). Thus, we could speculate that the lateralized network modules might result from the functional and structural hemispheric asymmetries.

Taken together, our results suggest an efficient modular organization in the WM network from early childhood and are consistent with modular behavior reported in previous structural and functional brain network studies and more importantly, a lateralized developmental pattern in some of the modules. The inter-modular hubs are the main connectors between modules and their existence in frontal and posterior cortex in the developing brain are consistent with previous WM network (Yap et al., 2011) and cortical thickness network (Khundrakpam et al., 2013) analysis. Resting state functional networks have also reported a high density of strong functional connections in posterior cortex (Achard et al., 2006). Thus, we could speculate that the inter-modular hubs uncovered in this study are well-established at childhood and are responsible for the connections between different functional systems of the developing brain.

A few methodological issues need to be addressed. Two drawbacks of our study include the acquisition of DTI data with six diffusion directions at low *b* values of 1000 s/mm2 and the use of deterministic tractography which will give errors in such an unsupervised tractography method given abrupt terminations at low FA crossing fiber regions or erroneous connections due to errors in the primary eigenvector direction. Multiple gradient directions can reduce the uncertainty of the primary eigenvector direction and limit potential bias as a function of tract orientation, both concerns for deterministic tractography of WM tracts (Landman et al., 2007). However, a recent study from our group has demonstrated six-direction data can also provide average diffusion measures like FA over a specific tract with comparable robustness to 30- or 60-direction data and yield appropriate parameter values for many major WM tracts (Lebel et al., 2012), which is encouraging as our edge weights were calculated based on the average FA of all voxels over the whole tract connecting two nodes. However, this does not overcome potential false positive connections or missed connections from the deterministic tractography algorithm. We attempted to minimize the former by invoking a minimum FN between regions but an incorrect connection that is consistent among the subjects within a group would still be included in the network analysis. DTI data with more than six directions also permit other advantages such as alternative analysis methods (e.g., probabilistic tractography) (Dennis et al., 2013). Higher *b* values than typically acquired are also advantageous for resolving crossing fibers and increasing the accuracy of tractography derived connections (Tournier et al., 2008). Another limitation of the study is that the age ranges of the groups covered a 3.1 to 4.7 year age range for the four youngest groups. In this study, a general linear model was applied to remove those age effects within all groups before performing the between-group comparison. In future study, smaller age ranges within groups may provide more specific indices of timing for the WM network maturation. Third, a FN threshold of 10 was applied to minimize the inclusion of random connections between two cortical regions. Currently, there are no standard approaches in determining the threshold value for the number of connecting fibers between regions as small thresholds such as 3 streamlines (Shu et al., 2011) produced networks with large sparsity with many spurious connections. Thus, our choice of higher threshold reduces, but does not eliminate, the risk of false-positive connections due to noise or the limitations in deterministic tractography. Recently, a threshold of 10 connecting streamlines or more was also applied in a brain network study (van den Heuvel et al., 2012) in which they considered that edges comprising fewer than 10 streamlines were potentially spurious and were deleted from the connection matrix. To examine the influence of the threshold, we tested a range of thresholds from 3 to 10 fibers and results including all network parameters are shown in **Table A1**. Although the network efficiency decreased as the sparsity decreased, the small worldness of the network remained intact. Most importantly, the group differences among adjacent age groups also remained consistent across all applied thresholds which indicates that the network comparison results are not sensitive to the threshold choices. Cortical regions in our study are defined by an a priori volumetric template (AAL) that was employed to automatically parcellate the entire cerebral cortex into different regions. Different templates used in various studies might cause discrepancy in the specific results, though the main trend of the network properties is expected to remain intact.

Various weighting functions for cortical-cortical connections have been applied in previous brain network analyses of brain development including 1/mean diffusivity (Hagmann et al., 2010) and proportional FN (Dennis et al., 2013), whereas we used the product of tract FA (known to increase exponentially with development over this age range but at unique rates per tract— Lebel et al., 2008) and AAL regional volume-normalized FN that has been used by others in studies of Alzheimer's Disease and aging (Lo et al., 2010; Brown et al., 2011). Other diffusion indices such as mean diffusivity (MD), axial diffusivity (AD) or radial diffusivity (RD) could have been examined instead of FA as a basis of "weighting" the network connections. However, to our best knowledge, while a few studies have applied MD as an edge weighting function (Hagmann et al., 2010; Li et al., 2012), none have used AD or RD. While changes in MD and FA for the WM typically occur together during maturation, with *MD* values decreasing and *FA* values increasing, the processes by which the two parameters change are theoretically different (Schmithorst et al., 2002; Huppi and Dubois, 2006) and they do not change at the same rate (Lebel et al., 2008). Axial and RD, under certain circumstances, may be more specific to underlying biological processes, such as myelin and axonal changes (Song et al., 2002). A recent study has demonstrated changes of FA in corticospinal tract and anterior corona radiata during development (2 to 40 years) that were attributed to the different rate changes in AD and RD (Faria et al., 2010). Thus, one would expect different WM network organization if using different weighting functions. Therefore, future studies could consider using multiple diffusion tensor measures such as FA, MD, AD and RD.

In conclusion, a graph theoretical approach was used to demonstrate age-related alterations in the large scale network properties of the developing WM network from early childhood (6 years) to adulthood (30 years). It was shown that increased network weight signifies a reshaping of the WM network from early childhood to late childhood with increased integration and decreased segregation. These findings are compatible with the notion that structural and functional brain networks become stable after late childhood. Our results also have implications for understanding how the modular organizational alterations in the large-scale structural brain networks underlie maturation of cognitive function in brain development. This study may pave the way for developing novel methods for analyzing disrupted brain connectivity in neurodevelopmental disorders.

## **ACKNOWLEDGMENTS**

Operating funding from the Canadian Institutes of Health Research, and salary awards from the China Scholarship Council (Min Liu) and Alberta Innovates—Health Solutions (Christian Beaulieu).

## **REFERENCES**

Dennis, E. L., Jahanshad, N., McMahon, K. L., de Zubicaray, G. I., Martin, N. G., Hickie, I. B., et al. (2013). Development of brain structural connectivity between ages 12 and 30: a 4-Tesla diffusion imaging study in 439 adolescents and adults. *Neuroimage* 64, 671–684. doi: 10.1016/j.neuroimage.2012.09.004

matter microstructure during adolescence. *Neuroimage* 39, 52–61. doi: 10.1016/j.neuroimage.2007.07.043

MRI. *Cereb. Cortex* 17, 2407–2419. doi: 10.1093/cercor/bhl149

*Sci. U.S.A.* 98, 676–682. doi: 10.1073/pnas.98.2.676

*Neurosci. Lett.* 355, 25–28. doi: 10.1016/j.neulet.2003.10.063

*Neuropsychol. Rev.* 20, 362–375. doi: 10.1007/s11065-010-9145-7

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 01 June 2013; accepted: 09 October 2013; published online: 01 November 2013.*

*Citation: Chen Z, Liu M, Gross DW and Beaulieu C (2013) Graph theoretical analysis of developmental patterns of the white matter network. Front. Hum. Neurosci. 7:716. doi: 10.3389/fnhum. 2013.00716*

*This article was submitted to the journal Frontiers in Human Neuroscience.*

*Copyright © 2013 Chen, Liu, Gross and Beaulieu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## **APPENDIX**

**Table A1 | White matter network parameters derived from different fiber-number node-to-node connection thresholds in all age groups.**

*Adjacent groups show significant differences or trend in network properties are highlighted in shaded cells. FN: fiber number, Eglob: global efficiency, Eloc: local efficiency, Q: modularity, σ: small worldness.*

## Age-related changes in brain structural covariance networks

#### *Xinwei Li 1,2, Fang Pu2, Yubo Fan2, Haijun Niu1,2, Shuyu Li 1,2\* and Deyu Li 2 \**

*1 State Key Laboratory of Software Development Environment, Beihang University, Beijing, China*

*2 Key Laboratory for Biomechanics and Mechanobiology of Ministry of Education, Department of Biomedical Engineering, School of Biological Science and Medical Engineering, Beihang University, Beijing, China.*

#### *Edited by:*

*Yong He, Beijing Normal University, China*

#### *Reviewed by:*

*Xiaobo Li, Albert Einstein College of Medicine, USA Zhang Chen, University of Alberta, Canada*

#### *\*Correspondence:*

*Shuyu Li and Deyu Li, Key Laboratory for Biomechanics and Mechanobiology of Ministry of Education, School of Biological Science and Medical Engineering, Beihang University, No. 37 Xueyuan Road, Haidian District, Beijing 100191, China. e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
Previous neuroimaging studies have suggested that cerebral changes over normal aging are not simply characterized by regional alterations, but rather by the reorganization of cortical connectivity patterns. The investigation of structural covariance networks (SCNs) using voxel-based morphometry is an advanced approach to examining the pattern of covariance in gray matter (GM) volumes among different regions of the human cortex. To date, how the organization of critical SCNs change during normal aging remains largely unknown. In this study, we used an SCN mapping approach to investigate eight large-scale networks in 240 healthy participants aged 18–89 years. These participants were subdivided into young (18–23 years), middle aged (30–58 years), and older (61–89 years) subjects. Eight seed regions were chosen from widely reported functional intrinsic connectivity networks. The voxels showing significant positive associations with these seed regions were used to describe the topological organization of an SCN. All of these networks exhibited non-linear patterns in their spatial extent that were associated with normal aging. These networks, except the primary motor network, had a distributed topology in young participants, a sharply localized topology in middle aged participants, and were relatively stable in older participants. The structural covariance derived using the primary motor cortex was limited to the ipsilateral motor regions in the young and older participants, but included contralateral homologous regions in the middle aged participants. In addition, there were significant between-group differences in the structural networks associated with language-related speech and semantics processing, executive control, and the default-mode network (DMN). Taken together, the results of this study demonstrate age-related changes in the topological organization of SCNs, and provide insights into normal aging of the human brain.

**Keywords: connectivity, structural covariance network, normal aging, neuroimaging, sensorimotor, neurocognition**

## **INTRODUCTION**

The majority of neuroimaging studies of aging have reported a consistent pattern of gray matter (GM) volumetric reductions in the human cortex, involving mainly prefrontal regions, parietal, and temporal association cortices, and the insula and cingulum (Resnick et al., 2003; Sowell et al., 2003; Raz et al., 2005; Du et al., 2006; Terribilli et al., 2011). However, more and more studies have demonstrated that cerebral changes with normal aging are not simply characterized by regional alterations but rather by the reorganization of cortical connectivity patterns (O'Sullivan et al., 2001; Koch et al., 2010; Wu et al., 2011b; Zhu et al., 2012). Using diffusion tensor imaging (DTI), several studies have consistently reported a diffuse loss of axonal integrity in senior populations (Salat et al., 2005; Pagani et al., 2008; Madden et al., 2009), which allowed for inferences regarding changes of structural connectivity in older people compared with younger adults. Building on complex network analysis methods, Gong et al. (2009) reported a reduction in overall cortical connectivity, decreased local efficiency, and a shift in regional efficiency from

parietal and occipital to frontal and temporal neocortex in older brains.

The investigation of structural covariance networks (SCNs) using structural magnetic resonance imaging (sMRI) is another useful method to explore structural brain networks. This approach mainly characterizes the pattern of structural covariance in GM morphology (e.g., volume, thickness and surface area) between brain regions using a general linear model (GLM) framework (Mechelli et al., 2005; Lerch et al., 2006; Nosarti et al., 2010; Zielinski et al., 2010; Montembeault et al., 2012; Soriano-Mas et al., 2012). Many studies have demonstrated underlying relationships among brain areas using structural correlation by sMRI, anatomical connectivity by DTI, and functional correlation by resting-state functional MRI. For example, He and colleagues found that structural networks based on cortical thickness measurements were compatible with known functional networks (He et al., 2007). Greicius et al. measured functional connectivity using independent component analysis and anatomical connectivity using DTI and found there existed white matter tract structural connections between functionally connected regions (Greicius et al., 2009). A recent study reported agreement in the correlations in GM thickness and underlying fiber connections across brain areas, but more information was included for the thickness network than the fiber network (Gong et al., 2012). In addition, Seeley and colleagues demonstrated that SCNs using voxel-based morphometry (VBM) were able to recapitulate the functional connectivity network topologies (Seeley et al., 2009). The covariance of different cortical regions in their GM volumes was considered to be the result of mutual trophic influences (Ferrer et al., 1995) or common experience-related plasticity (Draganski et al., 2004; Mechelli et al., 2004). The consistency among these three networks provides substantial support for SCNs serving a measure of network integrity in cross-sectional studies.

Of note, the analysis of SCNs using VBM has been successfully applied to map the eight SCNs and explore how neural systems build large-scale structural covariance during development (Zielinski et al., 2010). One previous aging study using structural covariance approaches compared two populations (younger vs. older subjects) and reported reduced structural associations in order adults, specifically in high-order cognitive networks (Montembeault et al., 2012). However, because subjects in only two age categories were involved in that experiment, it could not determine the aging trajectories of these sensorimotor and high-order cognitive networks. Zielinski et al. (2010) found that there were non-linear trajectories of primary visual, auditory, and sensorimotor networks during development. Moreover, these networks were provisionally established by early childhood, but underwent significant expansion in early adolescence before contraction or pruning in late adolescence. Many studies have shown that age-related atrophy of some neural regions follows variable, non-linear patterns (Allen et al., 2005; Kennedy et al., 2009; Terribilli et al., 2011). However, it remains unknown whether and how age-related changes of primary sensorimotor and high-order cognitive SCNs exhibit non-linear trajectories during normal aging.

Here, we used an SCN mapping approach to investigate eight large-scale networks in 240 healthy participants aged 18–89. These participants were subdivided into young (18–23 years), middle aged (30–58 years), and older (61–89 years) subjects. Given that SCNs subserving language, social–emotional, and executive control functions have shown gradual deterioration through normal aging, we expected to observe age-related trajectories in the deterioration of these large-scale networks. To test our hypothesis, we first investigated the aging trajectories of three sensorimotor and five high-order cognitive SCNs in three groups composed of an equal number of subjects at different ages, and then compared these SCNs differences between the groups.

#### **MATERIALS AND METHODS**

#### **PARTICIPANTS**

Three hundred sixteen right-handed, healthy subjects were selected from the Open Access Series of Imaging Studies (OASIS) cross-sectional database (http://www.oasis-brains.org) (Marcus et al., 2007). Data from seven subjects were excluded from further analysis because of poor image quality or image preprocessing. Because of relatively few subjects aged 30–60 years, we first selected 80 subjects from this age range as the middle-aged group (30 males and 50 females). Then, we selected 80 separate subjects around 75 years of age as the old group. Finally, we selected the youngest 80 subjects matched for gender as the young group. The names and characteristics of the groups are shown in **Table 1**. All subjects were evaluated using the Mini-Mental State Examination (MMSE) (Folstein et al., 1975) and Clinical Dementia Rating (CDR) scales (Morris, 1993; Morris et al., 2001). MMSE scores were higher than 29, and CDR scores were all zero. For demographic data on all subjects, see Marcus et al. (2007). This dataset has been used in several previous studies (Bakkour et al., 2009; Fjell et al., 2009; Salat et al., 2009; Li et al., 2011).

## **IMAGE ACQUISITION**

For each subject, three to four individual T1-weighted magnetization-prepared rapid gradient echo (MP-RAGE) images were acquired on a 1.5T Vision scanner (Siemens, Erlangen, Germany) within a single session. Head movement was minimized with cushioning and a thermoplastic facemask. Images were motion corrected and averaged to create a single image with a high contrast-to-noise ratio. MP-RAGE parameters were empirically optimized for gray/white contrast: *TR* = 9*.*7 ms; *TE* = 4 ms; flip angle = 10◦; slice number = 128; resolution = 256 × 256 (1 × 1 mm); thickness = 1*.*25 mm.

#### **IMAGE PREPROCESSING**

Structural MR images were processed using a technical computing software program (MATLAB 2010; The MathWorks Inc., Natick, Mass) and Statistical Parametric Mapping software (SPM 8; The Wellcome Department of Imaging Neuroscience, London, UK). Following the inspection of image artifacts, image preprocessing was performed with the VBM8 toolbox (http://dbm.neuro.uni-jena.de/vbm/). Briefly, all native-space MRIs were segmented to extract GM based on an adaptive maximum *a posteriori* technique (Rajapakse et al., 1997) and partial volume estimation method (Tohka et al., 2004) without the need for *a priori* tissue probability information. In addition, a spatially adaptive non-local denoising filter (Manjon et al., 2010) and a hidden Markov random field model (Rajapakse et al., 1997) were applied to minimize the level of noise in the resulting GM segments. Subsequently, the high-dimensional DARTEL (diffeomorphic anatomical registration using exponentiated lie algebra) approach provided non-linear deformation to normalize the images to the DARTEL template in Montreal Neurological Institute (MNI) space, which was derived from 550 healthy control subjects (Ashburner, 2007). Non-brain tissue was also

#### **Table 1 | Characteristics of the subjects in this study.**

*Abbreviations: F, female; M, male.*

removed. Additionally, the Jacobian determinants derived from the spatial normalization were used to modulate the GM value for each voxel to preserve the total amount of GM from the original images (Good et al., 2001). We used non-linear components only, which allowed us to analyze relative (i.e., corrected for individual brain size) differences in regional GM volume. Finally, the resulting modulated and normalized images were smoothed with a 12 mm full width at half maximum isotropic Gaussian kernel.

## **EXTRACTION OF SEED VOLUMES**

To assess the structural covariance pattern of each large-scale network, we extracted individual GM volumes from eight seed regions of interest (ROIs). These ROIs included the primary visual, auditory, and motor cortex, as well as language-related speech and semantic areas, areas related to salience and executive control, and the default-mode network (DMN). We based our ROIs on previous studies (Zielinski et al., 2010; Montembeault et al., 2012). For each region, seeds were defined with the MarsBar ROI toolbox (http://marsbar.sourceforge.net/) as 4-mm radial spheres centered at the following MNI coordinates: right calcarine sulcus (9, −81, 7), right Heschl's gyrus (46, −18, 10), right precentral gyrus (28, −16, 66), left inferior frontal gyrus, pars opercularis (IFGo) (−50, 18, 7), left temporal pole (−38 10, −28), right frontoinsular cortex (38, 26, −10), right dorsolateral prefrontal cortex (DLPFC) (44, 36, 20), and right angular gyrus (46, −59, 23). In addition, we selected contralateral seeds by changing the sign on each seed's x coordinate to perform similar SCNs analyses (**Table A1**, and **Figure A1**).

## **STATISTICAL ANALYSIS**

A voxel-wise statistical analysis was performed on the GM images using the GLM as applied in SPM8. Eight multiple regression models were used to test the strength of the structural covariance between each seeds and all other regions across whole-brain GM for each age group. In each regression model, the extracted mean GM volume from each ROI was entered as a covariate of interest, and gender as a confounding covariate. Because of the unequal sample size across the genders, we removed the effects of gender on the SCN patterns in the correlation analysis. These statistical analyses identified voxels that had a positive covariance with each *a priori* selected ROI in each group. The resulting correlation maps were using height and extent thresholds at *P <* 0*.*05 with family-wise error (FWE) correction. They were displayed on the MNI template to allow qualitative comparisons between the age groups using BrainNet Viewer software (http://www.nitrc.org/projects/bnv/). To quantify differences in normal aging trajectories across networks, we calculated the total number of significant positive ipsilateral, contralateral, and whole-brain voxels and plotted these across the age groups.

To further examine the effects of age on specific regional covariance, we performed between-group difference analyses of SCN patterns for any two groups according to the scientific literature (Lerch et al., 2006). For any pair of voxels in two groups, their structural correlation may have different slopes, and the difference in slope may represent the difference in their structural association. The difference in slope was tested using a classic interaction linear model:

$$V\_i = \beta \alpha + \beta\_1 V\_j + \beta\_2 \text{Group} + \beta\_3 (V\_j \times \text{Group}) + \varepsilon\_i$$

where *Vi* and *Vj* represented the GM volumes of a pair of voxels in two groups. The Group component was modeled using treatment contrasts, and significance tested using the Student's *t* statistic. Specific *t*-value contrasts were established to map the voxels that expressed a significantly different structural association between any two groups. The threshold for the resulting statistical parametric maps was given at *P <* 0*.*05 for the height and extent thresholds, with FWE multiple comparisons correction.

## **RESULTS**

For each network, the age-related SCN trajectory was identified by the spatial extent of each SCN. These results were demonstrated by functional domain, as described in the following sections. We also explored regions that significantly differed in their structural associations between any two groups.

## **PRIMARY SENSORY AND MOTOR NETWORKS**

Seeds within primary visual cortex (right calcarine sulcus) produced SCNs with a relatively preserved pattern (**Figure 1**). However, there were small changes in the distributions of the covariance maps. In the young and middle-aged groups, the covariance regions mainly included the bilateral calcarine sulcus, lingual gyrus, cuneus and right lateral occipital gyrus, whereas only included bilateral calcarine sulcus, lingual gyrus, and cuneus in the old group. Primary auditory cortex (i.e., the right Heschl's gyrus) covaried with the bilateral insula and precuneus, the right posterior cingulate cortex (PCC), parahippocampus, and inferior frontal gyrus, and the left orbital-frontal cortex in the young participants, and underwent significant contraction in the middle-aged and old groups to include only the bilateral insula regions. There was a flat transition in the covariance maps during aging. Primary motor cortex (right precentral gyrus) correlated with the ipsilateral precentral regions in the young participants, which progressed at middle age to include the contralateral precentral and supplementary motor areas. In old participants, these structural associations were similar to the young participants.

In addition, within the primary sensory and motor networks, there were no significant differences observed in the structural association between any two groups when we compared the slopes of the structural associations on a voxel-by-voxel basis.

## **LANGUAGE-RELATED SPEECH AND SEMANTIC NETWORKS**

The two language-related speech and semantic networks followed similar variation patterns (**Figure 2**). From the young participants to the middle group, the covariance regions showed abrupt decreases. They then showed smaller changes at the following ages. For the language-related speech network, the covariance map involved the bilateral anterior insula, medial frontal, cingulate cortex, and temporal regions in the young participants, and then contracted to include seed autocorrelation in the older participants. The semantic regions correlated with the bilateral

**covariance networks. (A)** Statistical maps of regions significantly correlated with the seed region in each group. The results are presented as correlation coefficient values (*P <* 0*.*05, FWE corrected). **(B)** The plots of voxel counts by group indicate small covariance changes throughout the life-span. The y-axis represents the voxel number (×104). Abbreviations: Calcs, calcarine sulcus; HG, Heschl's gyrus; PreCG, precentral gyrus; R, right; L, left; Y, young group; M, middle-aged group; O, old group; CC, correlation coefficient.

temporal cortices, cingulate gyrus, insula and frontal regions in the young participants, and shrank sharply to include only the anterior temporal cortices during normal aging.

Comparing the other SCNs, language-related speech and semantic networks showed more age-related changes. Specifically, the left supplementary motor and superior temporal areas showed significant positive associations with the speech seed region in young adults, whereas this covariance disappeared in the middle-aged group (**Table 2** and **Figure 3**). Decreased positive associations between the left superior temporal region and left temporal pole were found in the middle-aged group compared with the young group (**Table 2** and **Figure 4**). There were

**Table 2 | Contrast analysis of structural covariance network trajectories.**

*The regions listed showed significant between-group differences (P < 0.05, FWE corrected). x, y, z coordinates are reported in standard MNI space. Abbreviations: DMN, default-mode network; ANG, angular gyrus; IFGo, inferior frontal gyrus, pars opercularis; TPole, temporal pole; DLPFC, dorsolateral prefrontal cortex; PFC, prefrontal cortex; PCC, posterior cingulate cortex; ACC, anterior cingulate cortex; SMA, supplementary motor area; STC, superior temporal cortex; SFC, superior frontal cortex; OFC, orbitofrontal cortex; R, right; L, left; Y, young group; M, middle-aged group; O, old group; BA, Brodmann area.*

significant positive associations between the left anterior and the PCC and left IFG in the young participants, and this covariance became a negative correlation in the old group (**Table 2** and **Figure 4**). When compared with the old group, there was a significant positive association between the right PCC and the semantic seed region in the young participants (**Table 2** and **Figure 5**).

#### **SALIENCE, EXECUTIVE CONTROL, AND DEFAULT-MODE NETWORKS**

The three networks associated with social–emotional and cognitive function showed similar patterns to the language-related speech and semantic networks, with more distributed structural covariance in the young participants and limited covariance patterns throughout the later stages (**Figure 5**). Specifically, the right frontoinsular cortex anchored covariance maps included extensive areas of the bilateral lateral and medial frontal cortex, temporal cortices, and cingulate regions. This SCN underwent significant shrinkage in the middle-aged group to include only bilateral insular regions. In the older subjects, this network somewhat extended to include the bilateral medial prefrontal regions. The right DLPFC seed covaried with the bilateral frontal, temporal, and cingulate cortices in the young participants, but shrank into a more focal distribution in the middle-aged and older subjects. A seed in right angular gyrus produced an aging SCN, including the bilateral angular gyrus, middle temporal, cingulum, and prefrontal regions in the young participants, and contracting to include the bilateral PCC and angular gyrus in the middle-aged group, and only the bilateral angular gyrus in the old group.

There were no significant differences observed in the structural associations between any two groups within the salience network. The left PCC showed a significant positive association with the executive control seed region and right angular region (DMN seed) in the young adults, whereas this covariance disappeared in the old group (**Table 2** and **Figure 5**). Furthermore, a less robust positive association between the right prefrontal region and DMN seed region was found in the middle-aged group compared with the young participants (**Table 2** and **Figure 4**).

## **DISCUSSION**

In this study, we investigated the age-related trajectories of eight large-scale networks using an SCN mapping approach in three distinct age groups. All networks exhibited non-linear patterns across normal aging in terms of the spatial extent of the network. Except for the primary motor network, these networks showed a more distributed topology in young participants, which shrank sharply to a more localized topology in the middle-aged group, and maintained this localized topology in the old group. Primary sensory and motor networks showed fewer age-related changes compared with high-order cognitive networks. Moreover, there were significant between-group differences in language-related speech and semantic networks, the executive control network, and the DMN. Taken together, our results provide evidence of variations in the topological organization of SCNs during normal aging.

**FIGURE 3 | Age-related changes in salience, executive control, and default-mode structural covariance networks. (A)** Statistical maps of regions significantly correlated with the seed region in each group. The results are presented as correlation coefficient values (*P <* 0*.*05, FWE corrected). **(B)** The plots of voxel counts by group indicate abrupt contraction in the middle-aged group and mild changes in the old group. The y-axis represents the voxel number (×104). Abbreviations: DMN, default-mode network; DLPFC, dorsolateral prefrontal cortex; FI, frontoinsular cortex; ANG, angular gyrus; R, right; L, left; Y, young group; M, middle-aged group; O, old group; CC, correlation coefficient.

**FIGURE 4 | Group differences for the young vs. the middle-aged subjects.** Significant between-group differences within the DMN **(A)**, speech **(B)**, and semantic **(C)** structural covariance networks were found. For each network, the region of interest (upper) and region showing the most significant structural association (lower) are presented on the right, and a plot of slop differences between the seed region and a 4-mm radius

sphere centered on the peak voxel are presented in the left. Voxels with *P*FWE *<* 0*.*05 are displayed. Abbreviations: DMN, default-mode network; ANG, angular gyrus; IFGo, inferior frontal gyrus, pars opercularis; TPole, temporal pole; PFC, prefrontal cortex; SMA, supplementary motor area; STC, superior temporal cortex; R, right; L, left; Y, young group; M, middle-aged group.

## **COMPARISON OF STRUCTURAL COVARIANCE NETWORKS AND OTHER NETWORKS**

The cerebral cortex is organized into networks of functionally complementary areas. New advances in modern neuroimaging techniques and quantitative analysis of complex networks have made the investigation of brain network topological organization possible. DTI is a useful tool for non-invasively mapping cortico-cortical anatomical connections by examining axonal integrity. Insight into the structural covariance of GM morphology (e.g., volume, thickness, and surface area) between brain regions is provided by sMRI. Resting-state fMRI has allowed for assessments of the strength of functional connections within a network by quantifying correlated activity [i.e., spontaneous, low-frequency fluctuations in the blood oxygen level-dependent (BOLD) signal] between brain regions at rest.

The current studies have demonstrated the underlying relationships among structural correlations, anatomical connectivity, and functional correlations. For example, structural networks based on cortical thickness measurements were consistent with known functional networks (He et al., 2007). A recent study

reported agreement in the correlations in GM thickness and underlying fiber connections across brain areas, but more information was included for the thickness correlation network than the fiber connection network (Gong et al., 2012). In addition, SCNs determine using VBM recapitulated the canonical intrinsic connectivity networks topologies (Seeley et al., 2009). The consistency among the SCNs and other networks provides substantial support for the use of SCNs as a measure of network integrity for cross-sectional studies.

## **AGE-RELATED CHANGING TRAJECTORY OF SCNs**

In this study, we employed three age groups to map the trajectory of SCN changes over age. Except the primary motor network, these SCNs appeared the similar non-linear pattern that had a distributed topology in young participants, a sharply localized topology in middle aged participants, and were relatively stable in older participants. This trajectory was similar to agerelated changes of integrated local efficiency of brain network (Wu et al., 2011b). They employed the graph theory analysis method and reported that the local efficiency in the young group was significantly larger than those of the middle and old groups, whereas no significant difference was found between the middle and old groups. The shrinkage of SCNs as well as the reduction of local efficiency might be explained by the regionally distributed pattern of GM atrophy (Bergfield et al., 2010). As for primary motor network, we found inverted V-curve tendency among three age groups in the SCN of primary motor cortex. The previous study reported increased inter-regional correlations between bilateral primary motor cortex in the aging brain (Chen et al., 2011). In our study, we found the structural covariance of primary motor cortex between two hemispheres in the middle group. Alternatively, we found that the volume of right precentral gyrus reduced with age increasing (**Figure A2**), which might reflect a compensation mechanism that a network may need to work harder by becoming overactive due to regional volumetric atrophy with the network (Reuter-Lorenz and Cappell, 2008).

## **INTRAGROUP PATTERNS AND BETWEEN-GROUP DIFFERENCES OF SCNs**

## *Primary sensory and motor networks*

In this study, we found that seeds within the primary visual, auditory, and motor cortices produced SCNs with smaller spatial extents (i.e., total voxel count exhibiting a significant correlation at the corrected threshold) compared with other networks. This is consistent with previous work (Lerch et al., 2006) with mapped anatomical correlations across cerebral cortex (MACACC) using cortical thickness. Similarly, they used the number of cortical region showing significant correlations with seed region to describe the strength of MACACC. They found the primary motor, sensorimotor, and visual areas had the lowest strengths of correlation, whereas the association cortices had the highest strength. This could be explained by the functions of association cortices, because these cortices receive and integrate inputs from multiple cortical and non-cortical sources, and distribute information to multiple areas.

In addition, primary sensory and motor networks showed fewer age-related changes compared with high-order cognitive networks. There were no significant differences observed in the structural association between any two groups when the slopes of the structural associations were contrasted on a voxel-by-voxel basis within the primary sensory and motor networks. This finding is in line with a recent study (Montembeault et al., 2012) in which between-group differences were only observed in high-order cognitive SCNs between young and old subjects, rather than primary sensory and motor networks.

#### *Language-related speech and semantic networks*

The language-related speech and semantic networks showed significant age-related changes, with a more localized topology with increasing age. These results could account for the decline of language abilities across in normal aging. Federmeier and colleagues observed age-related changes in the timing with which messagelevel information impacted semantic analysis using event-related potentials, which indicated that aging affects high-order language processes (Federmeier et al., 2003). Wierenga and colleagues also reported that healthy older adults had increased difficulty in word retrieval with unchanged semantic knowledge (Wierenga et al., 2008).

Within the language-related speech SCN, a reduced structural association was observed the between inferior frontal gyrus and the supplementary motor area in the middle-aged group compared with the young participants. The speech SCN linked the language and motor systems that enable speech fluency (Seeley et al., 2009). This reduced structural covariance may explain the decline of motor speech skills with the transference of language to a speech code in older adults. Similarly, there was a significant decrease in the covariance between the seed region and the anterior cingulate cortex (ACC) in the old group compared with the young participants. This could be related to the functional degradation of these areas. Takahashi et al. reported that significant age-related reductions in regional cerebral blood flow in the left IFG and the bilateral medial frontal gyri and ACC using single-photon emission tomography (Takahashi et al., 2005). In addition, as the ACC plays a critical role in the control of speech responses (Paus et al., 1993), this reduced covariance may explain some speech difficulties in senior populations.

In semantic SCN analysis, a reduced structural association between the left temporal pole and left superior temporal cortex was found in the middle-aged group compared with the young participants. Many studies have shown that the superior temporal cortex is associated with understanding spoken words (Demonet et al., 1992; Chao et al., 2002; Okada and Hickok, 2006). Thus, this reduced structural association may be related to a decline in auditory word comprehension with increasing age. In the comparisons of the young and old participants, the main differences we found were in the structural associations between the left temporal pole and PCC. The results of functional neuroimaging studies in young, healthy adults provide compelling evidence for the involvement of the PCC in memory retrieval. For example, PCC activation was elicited during recognition of thematic narrative information learned during training sessions (Maguire et al., 1999). The reduced structural association between these two regions may explain the decline of semantic memory in senior populations.

### *Salience, executive control, and default-mode networks*

The three high-order cognitive SCNs showed a similar pattern to the language SCNs, showing a distributed topology in young participants that shrank sharply to a localized topology in the middle-aged group, and maintained a localized topology in the old group. Cepeda and colleagues found agerelated decline of executive control processes by examining taskswitching performance (Cepeda et al., 2001). Kelly and colleagues also demonstrated the stability of executive control decreased in older subjects (West et al., 2002). Several studies have shown that the DMN is altered in old subjects (Andrews-Hanna et al., 2007; Koch et al., 2010; Sambataro et al., 2010), especially in anterior regions (Damoiseaux et al., 2008). Similarly, Chen et al. found the decreased inter-regional correlations with the DMN module (Chen et al., 2011). Cognitive decline with age, such as cognitive control (Persson et al., 2007) and working memory (Sambataro et al., 2010), is associated with decreased DMN connectivity.

In the executive-control SCN analysis, the main differences we found were in the structural association between the right DLPFC and PCC, which exhibited a positive correlation in the young participants and no association in the old group. The DLPFC shows increased activity when experimental stimuli are presented, and is thought to support on-task processing in attention tasks. The PCC, in contrast, shows decreased activity during stimulus presentation, and is thought to support offtask processing in attention tasks. Some studies have reported that activity in the DLPFC shows a negative correlation with activity in the PCC during the resting state (Greicius et al., 2003; Fransson, 2005), and the negative coupling between these regions is lower in older adults (Sambataro et al., 2010). Our findings are inconsistent with functional connectivity between these areas, which should be further investigated in future studies.

Finally, for the DMN, we found a reduction in the structural association between the right angular gyrus and right prefrontal cortex from the young to the middle-aged participants. In the comparison of the young and old group, we found age-related differences in the structural association between the right angular gyrus and PCC. Similarly, Wu and colleagues reported that the left angular gyrus has a significantly reduced correlation with the PCC in older compared with younger participants in a resting-state fMRI study (Wu et al., 2011a). As the PCC is a prominent region within the DMN, the reduced structural association between the seed regions and the PCC indicates that the DMN SCN shrinks with age.

#### **LATERALITY OF THE SCNs**

We found similar SCNs across both hemispheres (**Figure A1**), consistent with previous SCN analyses of the developing brain (Zielinski et al., 2010). However, the auditory and speech networks showed a notable exception. The SCN derived from the left seed of Heschl's gyrus was distinctly smaller than the one derived from the right seed. In contrast, the SCN derived from the left seed of the IFGo was distinctly larger than the one derived from the right seed, which may explain why speech networks are left-dominated (Powell et al., 2006; Xiang et al., 2010).

## **FURTHER CONSIDERATIONS**

To build upon this study, several issues need to be addressed. First, we assigned subjects in their 30s to the middle-aged group because there were fewer subjects in their 40s and 50s in our sample. This resulted in wide age range in the middle-aged group. In future, this age range should be refined by including more subjects. Second, a previous study has reported the brain functional connectivity pattern could be affected by the sex factor (Biswal et al., 2010). In future, it is interesting to explore the sexual differences when the SCNs change with age. Third, it should be highlighted that this study used ROIs based on previous studies (Zielinski et al., 2010; Montembeault et al., 2012). Future studies

## **REFERENCES**

organization of structural cortical network by using cortical thickness from MRI. *Neuroimage* 56, 235–245.

using ROIs extracted from ICNs of the same subjects may obtain more accurate definitions of seed regions.

## **CONCLUSIONS**

In this study, we used an SCN mapping approach to investigate eight large-scale networks in a large cohort of healthy participants. Our results show age-related changes in the topological organization of SCNs and provide insights into the normal aging process of the human brain.

## **ACKNOWLEDGMENTS**

This work was supported by the National Science Foundation of China (Nos. 81171403, 81071212, and 30700182) and the State Key Laboratory of Software Development Environment (No. SKLSDE-2011ZX-10).

of patients for the clinician. *J. Psychiatr. Res.* 12, 189–198.

a comparison of optimized voxelbased morphometry to manual volumetry. *Neurobiol. Aging* 30, 1657–1676.

a positron emission tomography study. *J. Neurophysiol.* 70, 453–469.

impact on working memory performance. *Neurobiol. Aging* 31, 839.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 26 December 2012; accepted: 08 March 2013; published online: 26 March 2013.*

*Citation: Li X, Pu F, Fan Y, Niu H, Li S and Li D (2013) Age-related changes in brain structural covariance networks. Front. Hum. Neurosci. 7:98. doi: 10.3389/fnhum.2013.00098*

*Copyright © 2013 Li, Pu, Fan, Niu, Li and Li. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.*

## **APPENDIX**

#### **Table A1 | Contrast analysis of structural covariance network trajectories derived from contralateral seeds.**

*The regions listed showed significant between-group differences (P < 0.05, FWE corrected). x, y, z coordinates are reported in standard MNI space. Abbreviations: DMN, default-mode network; ANG, angular gyrus; IFGo, inferior frontal gyrus, pars opercularis; TPole, temporal pole; FI, frontoinsular cortex; HG, Heschl's gyrus; IFC, inferior frontal cortex; PCC, posterior cingulate cortex; ITC, inferior temporal cortex; FG, fusiform gyrus; OFC, orbitofrontal cortex; ACC, anterior cingulate cortex; R, right; L, left; Y, young group; M, middle-aged group; O, old group; BA, Brodmann area.*

**FIGURE A1 | Trajectories of the structural covariance networks derived from contralateral seeds.** The plots of voxel counts by group (*P <* 0*.*05, FWE) indicate a similar pattern to the ipsilateral seed, except for the auditory and speech networks, which show differences in their spatial extents. The y-axis represents the voxel number (×104).

Abbreviations: Calcs, calcarine sulcus; HG, Heschl's gyrus; PreCG, precentral gyrus; IFGo, inferior frontal gyrus, pars opercularis; TPole, temporal pole; ANG, angular gyrus; DLPFC, dorsolateral prefrontal cortex; FI, frontoinsular cortex; R, right; L, left; Y, young group; M, middle-aged group; O, old group.

# **HUMAN NEUROSCIENCE**

## A longitudinal study of structural brain network changes with normal aging

## *Kai Wu1,2\*, Yasuyuki Taki 1,3,4, Kazunori Sato1, Haochen Qi 2, Ryuta Kawashima4,5,6 and Hiroshi Fukuda1*

*1 Department of Nuclear Medicine and Radiology, Institute of Development, Aging and Cancer, Tohoku University, Sendai, Japan*

*2 Department of Biomedical Engineering, School of Materials Science and Engineering, South China University of Technology, Guangzhou, China*

*3 Division of Medical Image Analysis, Department of Community Medical Supports, Tohoku Medical Megabank Organization, Tohoku University, Sendai, Japan*