Fifty-five right-handed healthy subjects (24 male ages ranging from 20 to 29 with mean = 22.92 and standard deviation = 2.89; 31 female ages ranging from 20 to 28 with mean = 21.71 and standard deviation = 2.19) were selected from the Open Access Series of Imaging Studies Database (http://www.oasis-brains.org/). All the subjects were questioned about their medical histories and use of psychoactive drugs before a trained physician began image acquisition. For details on the clinical and demographic information of the subjects, please refer to Marcus et al. (2007).

For each subject, three or four individual T1-weighted magnetization-prepared rapid gradient-echo images were acquired using a single Siemens 1.5 T Vision scanner with the following parameters: repetition time = 9.7 ms, echo time = 4.0 ms, inversion time = 20 ms, flip angle = 10°, sagittal orientation with 128 slices, and resolution = 1 × 1 × 1.25 mm. Multiple T1 images obtained for each subject were motion-corrected and then averaged to achieve an improved signal-to-noise ratio. For additional details on the post-processing information regarding the raw images, please refer to Marcus et al. (2007).

All the images were preprocessed through FreeSurfer (version 5.3.0), which can be freely downloaded online (http://surfer.nmr.mgh.harvard.edu/). Studies have assessed FreeSurfer's performance with images acquired from various MRI scanners or sequences and received reliable results for cortical and subcortical measurements (Khan et al., 2008; Shi et al., 2013; Mulder et al., 2014). Details on FreeSurfer's processing pipeline can be found in Dale et al. (1999) and Fischl et al. (1999). In short, the raw images were resampled into 256 × 256 × 256 isotropic resolutions with a voxel size of 1 × 1 × 1 mm. The intensity bias was subsequently corrected in a volume-based stream, and the skull was stripped, followed by volumetric labeling and white matter segmentation (Figure 1A). In the surface-based stream, the surface of white matter was extracted as the inner surface, and then nudged to the gray–pial interface to generate the outer surface (Figure 1B). The FreeSurfer surface was created as a mesh with a non-uniform grid, also known as a vertex. Hence, the cortical morphometric features can be measured using every vertex between the inner and outer surfaces (Figure 1C). Finally, the built-in Desikan-Killiany cortical atlas (Desikan et al., 2006) was applied to obtain the regional measurements (Figure 1D). This atlas presents the parcellation of 34 regions for each hemisphere on the basis of the structural patterns of the gyrus and sulcus (Supplementary Table 1).

Nine morphometric features were initially employed in this study (Figure 1E): (1) the number of vertices in each region, (2) regional gray matter volume, (3) surface area, (4) mean cortical thickness, (5) standard deviation of cortical thickness, (6) mean curvature, (7) Gauss curvature, (8) curvature index, and (9) fold index. In particular, the first two features are tallied as the sum of the vertices and voxels in a region, respectively. Surface area was computed as the total of the areas of all vertices falling within a region (Panizzon et al., 2009). The cortical thickness is measured as the distance between the inner and outer cortical surfaces at each vertex (Fischl et al., 1999). The curvature is calculated as the reciprocal of the radius of the inscribed circle for each vertex, and the signs of the gyrus and sulcus are opposite (Li et al., 2014). The Gaussian curvature of a surface at a given point is the product of the principal curvatures, which are the eigenvalues of the shape operator at the point (Pienaar et al., 2008). They measure how the surface bends by different amounts in different directions at that point. The folding index is calculated as the number of cortices buried within the sulcal folds, as compared with the number of cortices on the visible outer cortex (Schaer et al., 2008). Notably, all the surface computations were performed in the native space, enabling the above-mentioned features to be measured without deformation.

To form the network, the regional descriptor was defined as a feature vector, which is the concatenation of the nine morphometric features in each region, as shown in Figure 1. Therefore, the interconnected matrix is generated by computing the Pearson correlation coefficient for each pair of feature vectors (Figure 1E). Significant differences exist between the orders of magnitude of each morphometric feature (10⁻² to 10³). Thus, the z-score values were computed for each feature as the standardized values before the correlation calculation. The z-score values were calculated from each morphometric feature by first subtracting the mean of all values from each individual value and then dividing each remainder by the standard deviation of all the values, as realized by the normalization function in the statistical analysis software SPSS v22.0 (SPSS Inc., Chicago IL, USA). The correlation matrix was then binarized to generate the unweighted and undirected networks (Figure 1F), because this matrix captures the underlying anatomical connection patterns, while providing a simple method for investigating the network (He et al., 2007). To maintain both positive and negative high correlations, the absolute values of matrixes were applied for binarization. Lastly, such network properties as small-world configurations and definition of hubs (Figure 1G) were analyzed in accordance with graph theory.

For binarization, sparsity was used to express the extent of thresholding, which is the ratio of existing connections to the total possible ones. However, the appropriate sparsity remains an open question (Zhu et al., 2012). Therefore, the network properties were evaluated as a function of sparsity ranging from 20 to 40% with a step size of 1%. The range was determined using the requirements stating that the minimum sparsity guarantees no isolated node in the network, and the maximum one ensures small-worldness. The other advantage of using a certain range of sparsity is that it allows us to investigate the network properties at the same level for all the subjects, as the same number of existing connections are found in each case. Additionally, the sparsity of 23% was highlighted to enable direct comparison with previous studies.

The network properties were computed using the Brain Connectivity Toolbox (BCT) (Rubinov and Sporns, 2010) and Graph-theoretical Network Analysis (GRETNA) toolkit (Wang et al., 2015). SPSS was used for the single-feature normalization (see individual-based morphological network construction) and all the statistical analyses, while the spatial distribution of hubs was visualized using the BrainNet Viewer toolkit (Xia et al., 2013). The multiple comparison was corrected using the false discovery rate (FDR) proposed by Genovese et al. (2002).

To evaluate the reproducibility of morphometric features and the present method, the intraclass correlation coefficient (ICC) was applied. ICC was first defined by Mcgraw and Wong (1996) as the fraction of the variance of chosen graphic property between subjects to the total variance, which is the summation variance of between and within subjects of that property.

If the measurements of repeated scans are consistent for each subject, the ICC would be close to one. An ICC value above 0.75 is considered excellent, and one ranging from 0.6 to 0.75 is considered good (Bennett and Miller, 2010). ICC was computed using the “Reliability Analysis” function of SPSS.

Initially, apart from the subjects' respective brain networks, C_p and L_p were also derived from the corresponding random networks at each sparsity for each subject. Hence, the comparisons were conducted between the different subjects' brain networks and random networks (Figures 2A,B) or sparsities (Figures 2C–E). The comparison shows that C_p and L_p are profoundly higher than the random ones all over the sparsity range with the maximum p = 1.14 × 10⁻¹¹ (t-values ranged from 11.75 to 27.77 with FDR correction) and p = 6.15 × 10⁻³ (t-values ranged from 2.98 to 5.00 with FDR correction), as revealed using the independent two-sample t-test (Figures 2A,B, respectively). Likewise, γ was larger than one (max = 1.99, min = 1.37) throughout the whole sparsity range, while λ was close to one (max = 1.28, min = 1.02). Hence, as expected σ are found larger than one throughout the entire sparsity range (max = 1.56, min = 1.35) which confirms the existence of small-worldness. Overall, the individual networks exhibit significantly higher C_p than the random network, while remaining similar L_p. Thus, small-worldness is confirmed for each subject at every sparsity. Moreover, the increase and decrease of small-world configurations, as sparsity increased, is illustrated in Figure 2C. Such small-world configuration fluctuations are in accordance with the variation tendency documented in previous reports (He et al., 2007; Kong et al., 2015).

Furthermore, the sparsity of 23% is highlighted for convenient comparison with previous studies, including those based on individual-based morphological brain networks (Tijms et al., 2012; Kong et al., 2015), those involving population-based morphological brain networks with various regional descriptors (He et al., 2007; Sanabria-Diaz et al., 2010; Zhu et al., 2012), and functional brain network studies (Zhang et al., 2011). As listed in Table 1, the population-based morphological networks and functional network exhibit smaller results than the current findings, whereas the individual-based morphological network studies present relatively similar results, particularly that of Kong et al. (2015). This finding suggests that the individual morphological brain networks may demonstrate a stronger integration and segregation because the inter-individual variability is highly reserved (Kanai and Rees, 2011).

In addition to the assessment of small-world configurations, the spatial distribution of hubs was investigated for different subjects and sparsities. First, the effect of sparsities on hub identification was computed. The BC was averaged through all the subjects at each sparsity to investigate the extent of the sparsity range for each hub region (msjBC). A total of ten regions identified as hubs achieved a sparsity range of over 20% (Figure 3), including the bilateral entorhinal cortex (left 100.00%, right 85.71%), superior temporal gyrus (left 100.00%, right 100.00%), lateral occipital gyrus (left 95.24%, right 100.00%), frontal pole (left 71.43%, right 28.57%) and left caudal middle frontal gyrus (100.00%), left rostral anterior cingulate cortex (71.43%), left isthmus cingulate cortex (33.33%), left parahippocampal gyrus (33.33%). In addition, a total of 4 regions were regarded as hubs in a sparsity range of 4.76 to 14.29%, and the rest of the 52 regions were never observed presenting a hub at any sparsity. The specific extent of sparsity range for each node identified as a hub may be found in Supplementary Table 2.

Thereafter, the influence of individual differences on hub determination was investigated. Likewise, the BC was averaged throughout the sparsity range for each subject to determine the number of subjects presenting each region as a hub (mspBC). A total of 14 nodes are hubs for over 20% of the subjects (Figure 4), including the superior temporal gyrus (left 54.55%, right 43.64%), caudal middle frontal (left 38.18%, right 27.27%), entorhinal cortex (left 38.18%, right 25.45%), lateral occipital gyrus (left 32.73%, right 38.18%), parahippocampal gyrus (left 23.64%, right 29.09%), and left isthmus cingulate cortex (27.27%), left rostral anterior cingulate cortex (23.64%), left frontal pole (23.64%), left inferior temporal gyrus (21.82%). A total of 53 regions are considered hubs in 1.82 to 18.18% of the subjects, and just one region comprised no subjects. The specific number of subjects presenting each node as a hub is found in Supplementary Table 3.

The similarity of BC across subjects and sparsities was evaluated through the proximity matrix generated by SPSS (Figure 5). As illustrated in the figure (yellow: similar, blue: dissimilar), the inter-subject and inter-sparsity similarities of BC were observed, and the highly comparable BC was found in adjacent sparsities (close to diagonal) (Figure 5A). The individual differences in subjects significantly influence the BC value (Figure 5B). Clearly, the similarities were significantly high across sparsities. These findings are also in agreement with the above-mentioned results, which state that the hubs can be mostly retained for different sparsities, while hub identification varied profoundly for different subjects.

Furthermore, Figure 6 shows the visualization of hub distribution based on the average BC of all the subjects and the entire sparsity range. The results indicate that ten hubs were identified, comprising six heteromodal or unimodal associative regions and four paralimbic regions (Table 2). All these hub regions have been reported in at least one previous study.

The method's reproducibility was evaluated by measuring the ICC index of network properties for scans with acquisitions of two different time points in the same subjects. The ICC was investigated throughout the entire sparsity range. Only C_p, L_p, and BC were examined in this study because of the instability of random network generation processing.

The results indicated that C_p was highly reproducible (minimum ICC = 0.71, average ICC = 0.83), as plotted in Figure 7A. Moreover, the robustness of L_p (minimum ICC = 0.63, average ICC = 0.81) was fairly stable (Figure 7B), as the similarities to the BC (minimum ICC = 0.629, average ICC = 0.78) are shown in Figure 7C. Most of the results are significant, except for L_p and BC at sparsity of 40% (p = 5.6 × 10⁻² and p = 5.7 × 10⁻², separately). Hence, the reliability of the present method is well-documented in accordance with the standard proposed by Bennett and Miller (2010). The specific value of ICC for the three measurements across the entire sparsity range is shown in Supplementary Table 4.

The interregional morphological similarity has been observed and verified repeatedly in recent studies based on such morphometric features as cortical thickness and regional volume (Lerch et al., 2006; He et al., 2007; Bassett et al., 2008; Sanabria-Diaz et al., 2010). The interpretation of this observation implies that the brain regions exhibit covaring morphometric features that can form the connected structure, which may capture long-term neurobiological effects (Mechelli et al., 2005). However, the underlying mechanism of this covariance pattern remains elusive. Some conjectures have been debated in the scientific literature, including mutually trophic effects (Ferrer et al., 1995; Aid et al., 2007), environment-related plasticity (Maguire et al., 2000; Draganski et al., 2004; Mechelli et al., 2004), genetic influence (Schmitt et al., 2008), and normal development (Raz et al., 2005; Chen et al., 2011). The axonal tension theory has also been raised recently (Tijms et al., 2012; Kong et al., 2015), stating that the interconnected areas are becoming either thicker or thinner as a result of being pulled by a mechanical force (Van Essen, 1997).

Moreover, the relationships among morphological (M), functional (F), and anatomical (A) connectivities were investigated in previous studies. For instance, Reid et al. (2017) have demonstrated that the extent of accordance between F-M modalities varied remarkably across seed regions. Wang et al. (2015) have illustrated a tight coupling of F-M modalities in connectivity strength and network topologic organizations (i.e., modules, rich club, and motifs), further indicating that the neuropsychiatric disorders may considerably break this coupling. In addition, Gong et al. (2012) have found an A-M convergence in 35–40% of morphological connectivities, with the majority of the convergences observed in the positive morphological connectivities.

Further research in neurobiology would significantly assist the construction and exploration of morphological brain networks. In addition, the data on different modalities were typically used or analyzed separately because they were supposed to capture the distinct underlying phenomena. However, the correspondences in topological structures have been recently observed across modalities, suggesting that the combination of multimodal research will become increasingly important in neuroscience.

The spatial distribution of hubs (nodes with higher BC of average subjects and the entire sparsity range) in the present study was found to be strikingly similar to the results of functional studies. The consistent regions are the bilateral superior temporal gyrus, entorhinal cortex, and lateral occipital gyrus, and left isthmus cingulate, frontal pole (Table 2). Such accordance is quite thought-provoking, although the different criteria were employed to identify the hubs (BC, degree, and L_p). Hence, the findings may suggest that the individual-based morphological brain networks may exhibit an improved consistency with functional studies.

In addition, the hub spatial distribution of different subjects (averaged by sparsities) shows striking diversity in the present study (Supplementary Table 3). As illustrated in Figures 4, 5 (right), inter-individual variability is implied as a major influence on hub determination. Indeed, a number of studies have demonstrated unmistakable discrepancies in brain structures across individuals (Kanai and Rees, 2011). For example, genetic differences can lead to brain morphological changes (Thompson et al., 2001; Pol et al., 2006), as well as the variance in occupations (Maguire et al., 2000; Gaser and Schlaug, 2003) and skills (Mechelli et al., 2004). Moreover, cognition is found to be correlated with cerebral structures. For instance, individual differences in degrees of empathy were observed to be associated with gray matter density (Eres et al., 2015) and volume (Banissy et al., 2012). Psychiatric diseases such as epilepsy also involve unpredictable foci locations and numbers (Engel, 2006). Therefore, individual-based investigations are indispensable in clinical and scientific research.

Hence, investigation of the organization of morphological brain networks, along with functional and anatomical network organization on a case-by-case basis is both intriguing and indispensable. The individual-based accordance across modalities may also be an indicator in future cognitive and psychiatric research.

Several methodological issues are found in the present study, which should be addressed and solved in future research.

First, the interregional connections were measured as the Pearson correlation coefficient in the present study, as in many earlier studies (He et al., 2007; Chen et al., 2011; Shi et al., 2013). However, this measurement shows that the observed connections between any two regions are actually the summation of the direct and indirect correlations. A widely used approach for eliminating the influence of other regions is to alternatively use partial correlation (Bassett et al., 2008). The number of variables should be less than a number of samples of each variable in partial correlation computation. However, in the present study, a total of 68 regions (i.e., variables) with nine morphometric features (i.e., samples) are used, which goes against the regular partial correlation. Consequently, in future research, the partial correlation computation must be adjusted to increase the accuracy of results addressing interregional connections.

Second, the networks were binarized with a sparsity range from 20 to 40% for topological analyses, as in previous studies (He et al., 2007; Kong et al., 2015). This binarization leads to the exclusion of 60–80% of connection information simply because such connections were deemed “not important” based on the correlation algorithms. The underlying mechanism of morphological connections remains unclear (see interregional morphological similarity). Therefore, whether or not the “not important” connections are really useless has not been confirmed (Barrat et al., 2004). As such, further exploration of brain networks as the full-connected weighted graph is essential to obtaining additional insights into the human brain.

Third, the nine morphometric features applied in the present study were far from painstakingly selected. However, the anatomical information including cortical thickness, gray matter volume, surface area, and curvature is varied (He et al., 2007; Sanabria-Diaz et al., 2010). In addition, redundant information may have existed in the three curvature-related measurements (mean curvature, Gauss curvature, and curvature index). Therefore, an optimal selection of morphometric features will be further explored.

Fourth, the present study only investigated the individual-based morphological connectivities, whereas the functional and anatomical connections were found to demonstrate a mutual relationship with the morphological connections in the present and previous studies (Sui et al., 2012; Mueller et al., 2013). Hence, in future research, the multimodal images should be employed in individual-based brain network analyses.

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.