Single-subject morphological brain networks: connectivity mapping, topological characterization and test-retest reliability.

INTRODUCTION: Structural MRI has long been used to characterize local morphological features of the human brain. Coordination patterns of the local morphological features among regions, however, are not well understood. Here, we constructed individual-level morphological brain networks and systematically examined their topological organization and long-term test-retest reliability under different analytical schemes of spatial smoothing, brain parcellation, and network type.
METHODS: This study included 57 healthy participants and all participants completed two MRI scan sessions. Individual morphological brain networks were constructed by estimating interregional similarity in the distribution of regional gray matter volume in terms of the Kullback-Leibler divergence measure. Graph-based global and nodal network measures were then calculated, followed by the statistical comparison and intra-class correlation analysis.
RESULTS: The morphological brain networks were highly reproducible between sessions with significantly larger similarities for interhemispheric connections linking bilaterally homotopic regions. Further graph-based analyses revealed that the morphological brain networks exhibited nonrandom topological organization of small-worldness, high parallel efficiency and modular architecture regardless of the analytical choices of spatial smoothing, brain parcellation and network type. Moreover, several paralimbic and association regions were consistently revealed to be potential hubs. Nonetheless, the three studied factors particularly spatial smoothing significantly affected quantitative characterization of morphological brain networks. Further examination of long-term reliability revealed that all the examined network topological properties showed fair to excellent reliability irrespective of the analytical strategies, but performing spatial smoothing significantly improved reliability. Interestingly, nodal centralities were positively correlated with their reliabilities, and nodal degree and efficiency outperformed nodal betweenness with respect to reliability.
CONCLUSIONS: Our findings support single-subject morphological network analysis as a meaningful and reliable method to characterize structural organization of the human brain; this method thus opens a new avenue toward understanding the substrate of intersubject variability in behavior and function and establishing morphological network biomarkers in brain disorders.

Introduction

The human brain is universally appreciated as a highly integrative system wherein ongoing signaling, information exchange and processing occur in response to a wide variety of endogenous neurophysiological processes and external cognitive demands. There is a growing body of evidence supporting that the brain function is not solely attributable to the dynamics of individual regions but rather integrative processes and dynamic interactions across multiple distributed systems (Bullmore and Sporns 2009; Barch 2013; Sporns 2013, 2014; Cole et al. 2014; Stam 2014; Vertes and Bullmore 2015). This nature makes the brain particularly amenable to study with complex brain network analysis, a powerful tool to investigate how the entire assemblage of brain regions adaptively reorganizes in the face of various cognitive demands and brain disorders (Bullmore and Bassett 2011; Kelly et al. 2012; Borsboom and Cramer 2013; Craddock et al. 2013; Fornito and Bullmore 2015).

Human brain networks (i.e., the connectome) (Sporns et al. 2005; Biswal et al. 2010) can be constructed using multimodal neuroimaging techniques in vivo. Currently, functional MRI (fMRI) and diffusion tensor imaging (DTI) are the two most commonly used approaches to construct individual brain networks by estimating interregional functional connectivity (Biswal et al. 1995; Salvador et al. 2005) or axonal pathways (Hagmann et al. 2007; Iturria‐Medina et al. 2007), respectively. Besides fMRI and DTI, structural MRI (sMRI) has attracted increasing attention recently in delineating whole‐brain morphological connectivity patterns by calculating interregional morphological correlations across a cohort of participants (He et al. 2007; Bassett et al. 2008). Compared with fMRI/DTI, sMRI has distinct advantages in its easy access, high signal‐to‐noise ratio, and relative insensitivity to artifacts (e.g., head motion). Thus, an sMRI‐based network approach is promising to serve as another canonical tool in characterizing network‐level brain organization under both healthy and pathological conditions (Alexander‐Bloch et al. 2013a; Evans 2013). Nevertheless, it should be noted that this methodology can obtain only one network for a group of participants, thereby ignoring interindividual variability and making the examination of brain–behavior relationships and health‐disease classification impossible.

Accordingly, several new methods have been developed to construct individual brain networks based on sMRI data (Raj et al. 2010; Zhou et al. 2011; Tijms et al. 2012). Specifically, in terms of an axon tension theory where axon‐connectivity of cortical areas is believed to have an influence on morphology (Van Essen 1997; Hilgetag and Barbas 2005), Tijms and colleagues proposed an intuitive, test–retest (TRT) reliable method to construct individual morphological networks at a cube (i.e., 3 × 3 × 3 voxels) resolution. Although this cube‐based method keeps the 3D structure of the cortex intact, it ignores remarkable variability of geometry (e.g., shape and size) among different brain regions. Moreover, the rigid extraction of the cubes might not match well with the convolutions of the brain (Tijms et al. 2012). To overcome these limitations, Kong and colleagues introduce a new method to estimate interregional morphological connectivity in terms of the Kullback–Leibler (KL) divergence (Kullback and Leibler 1951) and demonstrate the success of this method in revealing longitudinal changes in morphological connectivity profiles of the thalamus after long‐term sleep deprivation.

Here, we extended the work of Kong et al. (2014) to construct whole‐brain morphological networks at individual‐level and further characterized their topological organizations at both global and nodal levels. Moreover, we systematically evaluated the influences of several analytical factors on the network topology including spatial smoothing, brain parcellation and network type. Finally, we examined the long‐term TRT reliability of this method in mapping morphological connectivity patterns and capturing their underlying topological architecture. Figure 1 illustrates the schematic representation of the main analytical process used in the current study. Our results showed that the derived morphological brain networks were specifically organized, analytical scheme‐dependent and long‐term TRT reliable, which suggest the current method as a potentially promising framework to associate morphological network organization with interindividual behavior differences in both typical and atypical populations. Notably, a similar study (Kong et al. 2015) was published recently while this work was under peer review. In comparison with what was done by Kong and colleagues, this work has several distinct features, such as the systematic evaluation of effects of different analytical strategies on morphological brain networks (we provide a more detailed comparison of these two studies and results in the discussion).

Figure 1

A flowchart illustrating the main analytical process in the current study. Briefly, individual structural images were first segmented into gray matter, white matter and cerebrospinal fluid (A). The gray matter maps (smoothed and nonsmoothed) were then divided into different numbers of regions according to prior brain atlases (AAL and HOA) (B). For each region, the gray matter volume values within it were extracted and used to estimate the probability distribution function (C). Subsequently, the KL divergence‐based similarity was calculated between any pair of regions in their probability distribution functions, resulting in a similarity matrix (D). The resultant similarity matrix was further thresholded into both binary and weighted networks (E), which could be visualized as graphs (F). Finally, several graph‐based network measures were employed to topologically characterize the graphs at both global and nodal levels (G). AAL, Anatomical Automatic Labeling atlas; HOA, Harvard‐Oxford atlas; PCUN, precuneus; TPOmid, temporal pole: middle temporal gyrus; L, left; R, right.

Methods

Participants

A publicly available TRT dataset (1000 Functional Connectomes, RRID: SCR_005361, http://fcon_1000.projects.nitrc.org/indi/CoRR/html/bnu_1.html) was used in the current study, which is a subset of the Connectivity‐based Brain Imaging Research Database at Beijing Normal University. The dataset contained a total of 57 right‐handed participants (male/female: 30/27; age: 19–30 years, mean = 23.09 ± 2.36 years) with no history of neurological and psychiatric disorders. All participants completed two MRI scan sessions (session 1 and session 2) within an interval of approximate 6‐weeks (40.94 ± 4.51 days).

MRI data acquisition

All MRI scans were performed on a 3T Siemens Tim Trio MRI scanner at the Imaging Center for Brain Research, Beijing Normal University. The T1‐weighted images were acquired using a magnetization prepared rapid gradient echo sequence with the following imaging parameters: repetition time = 2530 ms; echo time = 3.39 ms; inversion time = 1100 ms; slice thickness = 1.33 mm; flip angle = 7°; no interslice gap; 144 sagittal slices covering the whole brain; matrix size = 256 × 256; field of view = 256 × 256 mm2. Other modalities were not used in the current study and therefore were not described here.

Gray matter volume calculation

All data preprocessing (session 1 and session 2) were carried out with the VBM8 toolbox (http://dbm.neuro.uni-jena.de/vbm8) based on Statistical Parametric Mapping 8 (SPM8, RRID: SCR_007037, http://www.fil.ion.ucl.ac.uk/spm/). First, the raw MRI data were checked manually to ensure no obvious artifacts. Then, individual structural MRI images were segmented into gray matter (GM), white matter and cerebrospinal fluid using an adaptive Maximum A Posterior technique. The resultant GM images were subsequently normalized to the MNI space using a high‐dimensional “Diffeomorphic Anatomical Registration Through Exponential Lie Algebra” approach (Ashburner 2007) and further nonlinearly modulated to compensate for spatial normalization effects. The nonlinear modulation essentially corrected for individual differences in brain size. After these steps, a GM volume map was obtained for each participant (1.5 mm isotropic voxels).

Spatial smoothing

Spatial smoothing is a typically used step for voxel‐based morphology analysis that increases the signal‐to‐noise ratio and improves intersubject anatomical correspondence of sulci and gyri in the brain. However, on the other hand, this step may introduce spurious local anatomical connectivity (see below for the definition of anatomical connectivity in the current study) due to the fusion of signals for spatially adjacent regions. Considering that little is known about the effects of spatial smoothing on the subsequent connectivity mapping and topological characterization, we performed all the following analyses separately for GM volume maps with or without spatial smoothing (Gaussian kernel with 6‐mm full width at half maximum).

Construction of individual morphological brain networks

In the current study, we constructed large‐scale morphological brain networks for each participant based on their GM volume images. A brain network is comprised of a collection of nodes and edges interconnecting the nodes, wherein nodes represent brain regions and edges represent interregional similarity in the distributions of regional GM volume here.

Definition of network nodes

To define network nodes or brain regions, we parcellated the brain into different regions of interest (ROIs) in terms of prior brain atlases. A convergent finding from different neuroimaging modalities has shown that different brain parcellation schemes are associated with different topological organizations for the resultant brain networks (Wang et al. 2009; Sanabria‐Diaz et al. 2010; Zalesky et al. 2010). Thus, we employed two widely used structural brain templates of the Anatomical Automatic Labeling atlas (AAL, RRID: SCR_003550) (Tzourio‐Mazoyer et al. 2002) and the Harvard‐Oxford atlas (HOA, RRID: SCR_001476) (Kennedy et al. 1998; Makris et al. 1999) to divide the brain into different numbers of ROIs (90 for the AAL in Table S1 and 112 for the HOA in Table S2). This allows us to estimate the robustness of our findings against different parcellation schemes.

Definition of network edges

To estimate internodal or interregional network edges, we utilized a KL divergence‐based similarity (KLS) measure to quantify morphological connectivity between two regions (Kong et al. 2014). For each participant, we first extracted GM volume values for all the voxels within each ROI under each brain parcellation scheme. The probability density function of these values was then estimated using the kernel density estimation (KDE) (Rosenblatt 1956; Parzen 1962) with bandwidths chosen automatically (Botev et al. 2010). This analysis was performed using public Matlab code posed by Botev (function: kde, http://www.mathworks.com/matlabcentral/fileexchange/14034-kernel-density-estimator). Further, the probability distribution function (PDF) was calculated for the obtained probability density function. Subsequently, the KL divergence was calculated between any pair of ROIs in their PDFs. KL divergence is an index from probability theory to measure the difference between two probability distributions or the information lost when a probability distribution is used to approximate another from the perspective of information theory (Burnham and Anderson 2002). Formally, the KL divergence from distribution Q to P is defined as: where P and Q are two PDFs, n is the number of sample points (see “Results” for the selection of n in the current study). It is worth noting that is not equal to . To derive a symmetric measure, we calculated a variation of the KL divergence as follows:

Finally, the KLS is computed as: where e is nature exponential. Through this transformation, the KLS ranges from 0 to 1, where 1 is for two identical distributions. After all these analyses, four KLS‐based morphological connectivity matrices were generated for each participant (self‐connections were set to 0): one for AAL‐based parcellation on smoothed GM maps, one for AAL‐based parcellation on unsmoothed GM maps, one for HOA‐based parcellation on smoothed GM maps and the other for HOA‐based parcellation on unsmoothed GM maps.

Determination of the number of sampling points during KDE

Before application of KDE to human brain data, it is important to investigate how many sampling points are needed to generate stable estimation. Therefore, we first examined the effects of different numbers of sampling points on the estimation of regional probability density function. Specifically, we estimated the probability density function for each region of each participant under each combination between spatial smoothing and brain parcellation atlas when different numbers of sampling points were used (24, 25, 26, 27, 28, 29 and 210). Then, the Fréchet distance (Alt and Godau 1995) was calculated for each region between any pair of probability density functions that were estimated using adjacent numbers of sampling points (i.e., 24 and 25, 25 and 26, 26 and 27, 27 and 28, 28 and 29, and 29 and 210). The Fréchet distance provides a measure to quantify the similarity/dissimilarity between two curves that not only takes into account the location and ordering of the points along the curves but also can deal with curves with different lengths. After this procedure, two 57 (participants) × 90 (regions) × 6 (adjacent pairs in sampling number) Fréchet distance matrices were obtained when the AAL‐based parcellation was used: one for smoothed data, the other for nonsmoothed data. Similarly, two 57 (participants) × 112 (regions) × 6 (adjacent pairs in sampling number) Fréchet distance matrices were obtained when the HOA‐based parcellation was used. Subsequently, the four matrices were separately averaged along the first dimension (i.e., participants) to derive four group‐level, regional mean Fréchet distance matrices (two 90 × 6, the other two 112 × 6). Finally, any pair of nearby columns of each the resultant matrix was compared by using a nonparametric permutation test (10,000 iterations). We found that the differences in the regional mean Fréchet distances between nearby columns became nonsignificant from the comparison between column 3 and column 4 for all the combinations between spatial smoothing and brain parcellation (Table 1). This means that the similarities/dissimilarities between regional probability density functions estimated with 26 and 27 sampling points (i.e., column 3) were comparable with those derived from 27 and 28 sampling points (i.e., column 4), and so on. Accordingly, 27 sampling points was conservatively chosen in the current study to provide a trade‐off between stable estimation of regional probability density functions and computational complexity.

Table 1

The P values of pairwise comparisons of the Fréchet distances between curves estimated with nearby sample points

Categories	A versus B	B versus C	C versus D	D versus E	E versus F
Nonsmo‐AAL	0.016	0.051	0.215	0.451	0.071
Nonsmo‐HOA	0.477	0.432	0.439	0.390	0.272
Smo‐AAL	0.011	0.133	0.135	0.155	0.174
Smo‐HOA	0.005	0.012	0.093	0.135	0.211

Nonsmo, no spatial smoothing; Smo, spatial smoothing; AAL, Anatomical Automatic Labeling atlas; HOA, Harvard‐Oxford atlas. A – the Fréchet distances between curves estimated with 24 and 25 sample points; B – the Fréchet distances between curves estimated with 25 and 26 sample points; C – the Fréchet distances between curves estimated with 26 and 27 sample points; D – the Fréchet distances between curves estimated with 27 and 28 sample points; E – the Fréchet distances between curves estimated with 28 and 29 sample points; F – the Fréchet distances between curves estimated with 29 and 210 sample points.

Threshold selection

Prior to topological characterization of the derived morphological connectivity matrices, a thresholding procedure is typically used to exclude noisy elements. Here, we employed a sparsity threshold, S (defined as the ratio of the number of actual edges divided by the maximum possible number of edges in a network) to convert each matrix C  = [c ] into a binary and a weighted network by applying a subject‐specific KLS threshold: and a weighted network

This thresholding method ensures the same number of nodes and edges for the resultant networks across participants. Due to the lack of a conclusive method to select a single threshold, a consecutive sparsity range of 0.05 < S < 0.4 (interval = 0.02) was chosen, where the resultant networks have sparse properties (Achard et al. 2006; He et al. 2007; Wang et al. 2009) and are estimable for the small‐world attributes (Watts and Strogatz 1998). All the following network analyses were performed at each of the threshold level in this range, therefore resulting in functions or curves of sparsity for the topological measures listed below.

Network analysis

Based on the analyses above, we obtained 8 = 2 (spatial smoothing: yes vs. no) × 2 (brain parcellation: AAL vs. HOA) × 2 (network type: binary vs. weighted) morphological brain networks at each sparsity level for each participant. For each of these networks, we calculated both global (clustering coefficient, C p, characteristic path length, L p, local efficiency, E loc, global efficiency, E glob and modularity, Q) and nodal (nodal degree, k i, nodal efficiency, e i and nodal betweenness, b i) metrics with the GRETNA toolbox (Wang et al. 2015a). Detailed formulas, usages and explanations of these metrics in the brain networks can be found in our previous study (Wang et al. 2011) and in an excellent methodological review (Rubinov and Sporns 2010).

To determine whether the morphological brain networks were nonrandomly organized, all the global network measures were separately normalized by the corresponding mean of 100 matched random networks. The random networks were generated using a topological rewiring algorithm (Maslov and Sneppen 2002) which preserved the same number of nodes and edges and the same degree distribution as real brain networks. Typically, an efficient, small‐world and modular network should fulfill the following conditions: normalized E loc > 1 and normalized E glob ~ 1, normalized C p > 1 and normalized L p ~ 1 and normalized Q > 1. Notably, to simplify subsequent TRT reliability and statistical analyses, we also calculated the area under curve (AUC, i.e., the integral over sparsity range) to provide a threshold‐independent summary scalar for each global and nodal network metric of each participant under each analytical combination

TRT reliability

We utilized a common index of intra‐class correlation (ICC) (Shrout and Fleiss 1979) to investigate the TRT reliability of the current single‐subject method in mapping morphological connectivity patterns and characterizing their topological organization. Specifically, for each connectivity element (i.e., KLS value) or network metric (global or nodal) derived from each appropriate combination of the three factors mentioned above, individual values were first merged into a 57 × 2 matrix with rows corresponding to participants and columns corresponding to sessions. Using a one‐way analysis of variance (ANOVA), we then split the total sum of the squares into between‐subject (MSb) and within‐subject (MSw, i.e., residual error) sum of squares. The ICC was then calculated as: where k is the number of repeated observations per participant (2 here). ICC is close to 1 for reliable measures that show low within‐subject variance relative to between‐subject variance and 0 (negative) otherwise. Consistent with our previous study (Wang et al. 2011), the reliability was categorized into poor (0 < ICC < 0.25), low (0.25 < ICC < 0.4), fair (0.4 < ICC < 0.6), good (0.6 < ICC < 0.75) and excellent (0.75 < ICC < 1.0).

Statistical analysis

To determine whether different analytical factors of spatial smoothing, brain parcellation and network type significantly affect topological descriptions of the morphological networks, a three‐way repeated measures ANOVA was performed for each global network parameter (clustering coefficient, characteristic path length, local efficiency, global efficiency and modularity and their corresponding normalization versions). For nodal centralities (nodal degree, efficiency and betweenness), a two‐way repeated measures ANOVA was performed for each nodal metric under the AAL and HOA schemes, respectively, due to the different numbers of nodes. Analogously, the three‐way and two‐way repeated measures ANOVA were separately conducted to examine the effects of different analytical strategies on the TRT reliability (ICC values) of global and nodal network measures.

Results

KLS‐based morphological brain networks

Patterns of morphological similarity matrices

Figure 2A and B show the mean morphological similarity matrices under all combinations between spatial smoothing and brain parcellation scheme for data session 1 and session 2, respectively. Visual inspection found that the connectivity patterns were complex but specifically organized with strong interhemispheric morphological similarity between bilaterally homologous regions. This was further validated by statistical analyses showing that the similarity values between homotopic regions were significantly greater than the others in the matrix, a robust finding against the analytical choices of spatial smoothing, brain parcellation scheme and data session (t‐test, all P < 10−3). Moreover, quantitative spatial correlation analyses showed that the overall patterns of morphological similarity matrices were highly similar between session 1 and session 2 at both group (r > 0.99) (Fig. 2C) and individual levels (mean r > 0.91) regardless of the analytical strategies used.

Figure 2

Spatial pattern, intersession similarity and TRT reliability of KLS‐based morphological connectivity matrices. The mean KLS matrices (across participants) derived from session 1 (A) and session 2 (B) were highly correlated with each other, a robust finding against different choices of spatial smoothing and brain parcellation (C). TRT reliability analysis revealed high reliability for most elements in the matrices (D). KLS, KL divergence‐based similarity; ICC, intra‐class correlation; Nosmo, no smoothing; Smo, smoothing; AAL, Anatomical Automatic Labeling atlas; HOA, Harvard‐Oxford atlas.

TRT reliability of morphological similarity matrices

Intra‐class correlation‐based TRT reliability analysis on individual elements in the morphological similarity matrices demonstrated high reliability for interregional KLS values (Nosmo‐AAL: 0.816 ± 0.098; Nosmo‐HOA: 0.816 ± 0.140; Smo‐AAL: 0.905 ± 0.050; Smo‐HOA: 0.899 ± 0.056) (Fig. 2D). Specifically, under each brain parcellation scheme, more than 85% elements derived from nonsmoothed data and more than 97% elements derived from smoothed data exhibited excellent reliability. However, in contrast to the abovementioned finding of higher similarities for interhemispheric connections between bilaterally homologous regions, lower TRT reliability were found for these connections than the others in the matrices regardless of the analytical choices of spatial smoothing and brain parcellation scheme (t‐test, all P < 10−3). This implies that the similarities between bilaterally homologous regions are related to higher within‐subject or lower between‐subject variance.

Global organization of morphological brain networks

Small‐worldness, efficiency and modularity

Compared with random networks, the morphological brain networks exhibited larger values in the clustering coefficient, local efficiency and modularity but approximately equal values in the characteristic path length and global efficiency under each analytical combination of spatial smoothing, network type and brain parcellation. This resulted in a pattern of normalized clustering coefficient, local efficiency and modularity > 1 and normalized characteristic path length and global efficiency ~ 1 (Fig. 3). These findings suggest conserved global network organization of high‐efficient, small‐world and modular architectures for morphological brain networks. In addition, we also presented the modular structure derived from the group‐level mean morphological brain network based on spatially smoothed data under the AAL parcellation scheme (Fig. 4). Seven modules were found (Q = 0.628, z‐score = 27.030) that corresponded well with known neuroanatomical systems, and the modules were similar to those derived from diffusion and functional MRI data, such as the visual module of (medial) occipital system and subcortical module (Chen et al. 2008; Hagmann et al. 2008; He et al. 2009b; Meunier et al. 2009).

Figure 3

The normalized global network measures as a function of sparsity for morphological brain networks from nonsmoothed (A) and smoothed (B) data. All the normalized global network measures significantly deviated from 1 (dotted lines) regardless of different analytical strategies, indicating obviously different organizations of morphological brain networks from matched random networks. The results were from data session 1. Bin, binary network; Wei, weighted network; AAL, Anatomical Automatic Labeling atlas; HOA, Harvard‐Oxford atlas.

Figure 4

Matrix and brain surface presentations of modular structure for morphological brain networks. A group‐level backbone matrix (left) was derived from the mean morphological brain network across all participants based on spatially smoothed data under the AAL parcellation scheme. The backbone matrix was further reordered according to the optimal modular partition as shown in the brain surface (right). Seven modules were identified that corresponded well with known neuroanatomical systems, such as the medial occipital module and subcortical module. KLS, Kullback–Leibler divergence‐based similarity. AAL, Anatomical Automatic Labeling atlas.

Effects of spatial smoothing, brain parcellation and network type on global network measures

Generally, the results showed that the three factors significantly influenced numerical values of all the global network measures in a complex interactive pattern (Fig. 5). This suggests that quantitative characterization of morphological brain networks depends on the analytical strategies. No further post hoc analyses were done because of the complex patterns and more importantly comparisons of numerical values make no sense regarding the selection of optimal analytical strategy.

Figure 5

Effects of spatial smoothing, brain parcellation and network type on global network measures. The statistical analysis was done for each global network measure based on their areas under curve. Generally, the three factors significantly influenced quantitative descriptions of all the studied global network measures in a complex pattern. These results were from data session 1. The results were from data session 1. AUC, the area under curve; Bin, binary network; Wei, weighted network; AAL, Anatomical Automatic Labeling atlas; HOA, Harvard‐Oxford atlas. a, significant main effect of spatial smoothing; b, significant main effect of brain parcellation; c, significant main effect of network type. Of note, all the interactive effects were also significant.

TRT reliability of global network measures

Figure 6 shows the TRT reliability for all the global network measures as a function of sparsity threshold. All the global measures exhibited fair to good reliability (i.e., ICC ranged from 0.40 to 0.75) over most of the sparsity range studied (mean ICC over sparsity threshold and across network measures ranged from 0.523 to 0.648). Further analysis based on the AUCs revealed even higher reliability for the global network measures (mean ICC across network measures ranged from 0.613 to 0.781) (Fig. 7). In addition, we compared the reliability between original or first‐order global network measures (clustering coefficient, characteristic path length, local efficiency, global efficiency and modularity) and their normalized or second‐order versions and found no significant difference (P > 0.05).

Figure 6

The test–retest (TRT) reliability of global network measures as a function of sparsity for morphological brain networks from nonsmoothed (A) and smoothed (B) data. Generally, despite of different analytical strategies, the global network measures showed fair to excellent reliabilities (i.e., 0.4 < ICC < 1) over almost the whole sparsity range studied. Bin, binary network; Wei, weighted network; AAL, Anatomical Automatic Labeling atlas; HOA, Harvard‐Oxford atlas; ICC, intra‐class correlation.

Figure 7

The test–retest (TRT) reliability of global network measures based on the area under curve. Fair to excellent reliabilities (i.e., 0.4 < ICC < 1) were found for all global network measures irrespective of different choices of spatial smoothing, brain parcellation and network type. ICC, intra‐class correlation; Bin, binary network; Wei, weighted network; AAL, Anatomical Automatic Labeling atlas; HOA, Harvard‐Oxford atlas.

Effects of spatial smoothing, brain parcellation and network type on TRT reliability of global network topology

To examine whether the ICC of global network topology depends on different analytical strategies, a three‐way repeated ANOVA measures was conducted across different network measures. The results showed that only spatial smoothing significantly affected the reliability of global network measures (F 1,9 = 17.543, P = 0.002). No effects were observed for brain parcellation, network type or any interaction (P > 0.05). Post hoc comparisons showed that performing spatial smoothing significantly improved the reliability of global network measures (t 78 = 3.848, P < 10−3).

Local nodal characteristics of morphological brain networks

Hubs

All nodal results were based on the AUCs. Therefore, a total of 12 = 2 (spatial smoothing: yes and no) × 2 (network type: binary and weighted) × 3 (nodal centrality metrics: degree, efficiency and betweenness) nodal centrality maps were obtained for each participant under each brain parcellation scheme (AAL and HOA). We first calculated the Pearson correlations between any pair of these nodal centrality maps after averaging them across participants. This resulted in a 12 × 12 correlation matrix under the AAL and HOA parcellation scheme, respectively (Figs. 8A, 9A). We found that the correlation coefficients were significantly larger for nonsmoothed (AAL: 0.786 ± 0.179; HOA: 0.791 ± 0.170) and smoothed (AAL: 0.742 ± 0.213; HOA: 0.748 ± 0.208) data than those between nonsmoothed and smoothed data (AAL: 0.299 ± 0.206; HOA: 0.473 ± 0.188) (permutation test, 10,000 iterations, P < 0.001 for all comparisons). These findings indicate that spatial smoothing significantly affects the distribution of nodal centralities in the brain for morphological brain networks. Therefore, we separately presented nodal centralities for nonsmoothed (Fig. 8B for AAL and Fig. 9B for HOA) and smoothed (Fig. 8C for AAL and Fig. 9C for HOA) data and highlighted the top 10% regions (9 for AAL and 11 for HOA) with the highest centrality values, which may serve as potential hub regions.

Figure 8

Nodal centralities under the AAL parcellation scheme. Pairwise correlation analyses revealed that spatial smoothing significantly modulated spatial patterns of nodal centralities (A). Therefore, (B) and (C) were used to illustrate specific patterns of nodal centralities (weighted network analysis) for nonsmoothed and smoothed data, respectively. Regions with the highest centralities (top 10%) were highlighted in orange. A hub score was further calculated to identify regions that consistently showed high centralities over different network types and centrality metrics (D; left for no n‐smoothed data and right for smoothed data; also see Table 2). The results were from data session 1. All regional abbreviations can be found in Table S1. AAL, Anatomical Automatic Labeling atlas.

Figure 9

Nodal centralities under the HOA parcellation scheme. The results were presented in a similar manner to that in Fig. 8. All regional abbreviations can be found in Table S2. HOA, Harvard‐Oxford atlas.

Table 2

Regions consistently showing high centralities (i.e., hubs)

Name	Score	Category	Name	Score	Category
Nonsmo‐AAL			Smo‐AAL
Superior frontal gyrus, medial orbital (R)	6	Paralimbic	Superior frontal gyrus, orbital part (R)	6	Paralimbic
Anterior cingulate and paracingulate gyri (L)	6	Paralimbic	Inferior temporal gyrus (R)	5	Association
Anterior cingulate and paracingulate gyri (R)	6	Paralimbic	Inferior frontal gyrus, orbital part (R)	4	Paralimbic
Cuneus (L)	6	Association	Median cingulate and paracingulate gyri (L)	4	Paralimbic
Cuneus (R)	6	Association	Parahippocampal gyrus (L)	4	Paralimbic
Calcarine fissure and surrounding cortex (R)	5	Primary	Parahippocampal gyrus (R)	4	Paralimbic
Inferior occipital gyrus (R)	4	Association	Lingual gyrus (R)	4	Association
Angular gyrus (L)	4	Association	Supramarginal gyrus (L)	4	Association
Precuneus (R)	4	Association
Nonsmo‐HOA			Smo‐HOA
Cingulate gyrus, posterior division (R)	6	Paralimbic	Inferior temporal gyrus, anterior division (L)	6	Association
Frontal orbital cortex (R)	6	Paralimbic	Parahippocampal gyrus, posterior division (L)	6	Paralimbic
Intracalcarine cortex (R)	6	Primary	Cingulate gyrus, posterior division (R)	5	Paralimbic
Inferior temporal gyrus, anterior division (L)	4	Association	Parahippocampal gyrus, anterior division (R)	5	Paralimbic
Temporal fusiform cortex, posterior division (R)	4	Association	Temporal fusiform cortex, posterior division (R)	5	Association
Frontal medial cortex (R)	4	Paralimbic	Temporal fusiform cortex, posterior division (L)	4	Association
Subcallosal cortex (L)	4	Limbic	Amygdala (L)	4	Subcortical
Subcallosal cortex (R)	4	Limbic	Paracingulate gyrus (L)	4	Paralimbic
Paracingulate gyrus (L)	4	Paralimbic	Precuneous cortex (R)	4	Association
Paracingulate gyrus (R)	4	Paralimbic	Parahippocampal gyrus, anterior division (L)	4	Paralimbic
Cingulate gyrus, anterior division (R)	4	Paralimbic

Nonsmo, no smoothing; Smo, smoothing; AAL, Anatomical Automatic Labeling atlas; HOA, Harvard‐Oxford atlas; L, left; R, right.

To further locate brain regions that were consistently identified as hubs, we assigned each node a hub‐score (ranged from 0 to 6) indicating the times that the node fell within the top 10% nodes across nodal centrality metrics and network types. For example, if a region was identified as a hub by all the nodal centrality metrics (degree, efficiency and betweenness) of both binary and weighted networks, its hub‐score was 6. Figures 8D, 9D show the spatial distributions of nodal hub‐scores in the brain under the brain parcellation schemes of AAL and HOA, respectively. The most consistent hubs (scored > 3), as summarized in Table 2, were mainly limbic/paralimbic and association cortices but with a slight overlap between nonsmoothed and smoothed data irrespective of the brain parcellation schemes.

Effects of spatial smoothing and network type on nodal centralities

A two‐way repeated ANOVA was separately performed to test the effects of spatial smoothing and network type on the mean (across participants) values of nodal degree, efficiency and betweenness. The results showed that no matter which parcellation scheme was used, network type significantly affected nodal degree and efficiency (P < 0.05, corrected by false discovery rate procedure across centrality metrics) but not nodal betweenness, and the effects were further modulated by the other factor of spatial smoothing. Again, no further post hoc analyses were needed here.

TRT reliability of nodal centralities

The mean ICC values of nodal centralities under different analytical categories were visualized in Fig. 10. Generally, high TRT reliability was found regardless of nodal centrality metrics and analytical schemes employed (0.741 ± 0.076 overall all factors). To understand nodal reliability more deeply, we also calculated their pairwise spatial correlations across regions and obtained a 12 × 12 correlation matrix under each brain parcellation scheme (Fig. 11A for AAL and Fig. 12A for HOA). Again, we found that the correlation values among nodal reliability maps derived from nonsmoothed data (AAL: 0.688 ± 0.268; HOA: 0.517 ± 0.394) and from smoothed data (AAL: 0.703 ± 0.234; HOA: 0.670 ± 0.271) were significantly larger than those calculated between nonsmoothed and smoothed data (AAL: 0.235 ± 0.113; HOA: 0.157 ± 0.138) (permutation test, 10,000 iterations, P < 0.001 for all comparisons). These findings indicate that, in addition to nodal centralities, spatial smoothing remarkably affects spatial distribution of nodal reliability for morphological brain networks (Fig. 11B vs. Fig. 11C for AAL and Fig. 12B vs. Fig. 12C for HOA). Thus, we used a similar procedure to hub detection to identify regions that consistently exhibited high reliability (top 10%) over nodal centrality metrics and network types but for nonsmoothed and smoothed data (Fig. 11D for AAL and Fig. 12D for HOA). Overall, regions with high reliability were mainly association cortices of temporal, occipital and parietal regions depending on the spatial smoothing and brain parcellation scheme. Table 3 further summarizes the most consistent regions showing high reliability (scored > 3). Notably, the right superior occipital gyrus and middle occipital gyrus under the AAL parcellation and the left posterior division of inferior temporal gyrus and right superior division of lateral occipital cortex under the HOA parcellation were identified for both nonsmoothed and smoothed data. In addition, we compared nodal ICC values among the three centrality metrics (one‐way ANOVA) and found significant main effects (all P < 0.001) driven by higher reliability of nodal degree and efficiency than nodal betweenness. This was independent on choices of spatial smoothing, brain parcellation and network type. Finally, we also compared nodal ICC values between cortical and subcortical regions and consistently found no significant differences regardless of nodal centrality metrics and analytical strategies (all P > 0.05).

Figure 10

The test–retest (TRT) reliability of nodal centrality metrics under the Anatomical Automatic Labeling atlas (AAL) (A) and Harvard‐Oxford atlas (HOA) (B) parcellation schemes. Fair to excellent reliability (i.e., 0.4 < ICC < 1) was observed for almost all nodal centrality metrics under all analytical combinations of spatial smoothing, brain parcellation and network type. Notably, nodal degree and efficiency exhibited obviously higher reliability than nodal betwenness regardless of the parcellation schemes. ICC, intra‐class correlation; Bin, binary network; Wei, weighted network; Nosmo, no smoothing; Smo, smoothing.

Figure 11

The test–retest (TRT) reliability of nodal centralities under the Anatomical Automatic Labeling atlas (AAL) parcellation scheme. Pairwise correlation analyses revealed that spatial smoothing significantly modulated spatial patterns of nodal reliabilities (A). Therefore, (B) and (C) were used to illustrate specific patterns of nodal reliabilities (weighted network analysis) for nonsmoothed and smoothed data, respectively. Regions with the highest reliabilities (top 10%) were highlighted in orange. A reliability‐based hub score was further calculated to identify regions that consistently showed high reliability over different network types and centrality metrics (D; left for nonsmoothed data and right for smoothed data; also see Table 3). The results were from data session 1. All regional abbreviations can be found in Table S1.

Figure 12

The test–retest (TRT) reliability of nodal centralities under the Harvard‐Oxford atlas (HOA) parcellation scheme. The results were presented in a similar manner to that in Fig. 11. All regional abbreviations can be found in Table S2.

Table 3

Regions consistently showing high reliabilities

Name	Score	Category	Name	Score	Category
Nonsmo‐AAL			Smo‐AAL
Precentral gyrus (L)	4	Primary	Superior occipital gyrus (R)	6	Association
Middle frontal gyrus (R)	4	Association	Inferior parietal, but supramarginal and angular gyri (R)	5	Association
Superior occipital gyrus (R)	4	Association	Cuneus (R)	4	Association
Middle occipital gyrus (R)	4	Association	Middle occipital gyrus (R)	4	Association
Angular gyrus (R)	4	Association	Inferior parietal, but supramarginal and angular gyri (L)	4	Association
Middle temporal gyrus (R)	4	Association	Precuneus (R)	4	Association
			Superior temporal gyrus (R)	4	Association
Nonsmo‐HOA			Smo‐HOA
Angular gyrus (L)	5	Association	Middle temporal gyrus, posterior division (L)	4	Association
Frontal orbital cortex (L)	5	Paralimbic	Middle temporal gyrus, temporooccipital part (L)	4	Association
Lingual gyrus (R)	5	Association	Inferior temporal gyrus, posterior division (L)	4	Association
Parietal operculum cortex (R)	4	Association	Lateral occipital cortex, superior division (R)	4	Association
Precentral gyrus (R)	4	Primary	Parahippocampal gyrus, posterior division (R)	4	Paralimbic
Inferior temporal gyrus, posterior division (L)	4	Association
Inferior temporal gyrus, temporooccipital part (L)	4	Association
Lateral occipital cortex, superior division (R)	4	Association
Lateral occipital cortex, inferior division (L)	4	Association

Nonsmo, no smoothing; Smo, smoothing; AAL, Anatomical Automatic Labeling atlas; HOA, Harvard‐Oxford atlas; L, left; R, right.

Effects of spatial smoothing and network type on TRT reliability of nodal centralities

A two‐way repeated ANOVA was performed for each nodal centrality metric (degree, efficiency and betweenness) under each parcellation scheme (AAL and HOA). The results showed that spatial smoothing and network type significantly modulated nodal reliability in an interactive manner regardless of the brain parcellation schemes and nodal centrality metrics used (all P < 0.05). Post‐hoc comparisons revealed that performing spatial smoothing significantly improved nodal reliability of all the three centrality metrics for both binary and weighted networks but to different extents (all P < 0.001). In addition, binary network analysis outperformed weighted network analysis with respect to reliability of nodal betweenness for both smoothed (P = 0.037 for AAL and <0.001 for HOA) and nonsmoothed data (P < 0.001 for AAL and HOA).

Relationship between nodal centrality and ICC

To determine whether nodal ICCs were related to their centralities in morphological brain networks, we calculated the Pearson correlation coefficients across regions between each nodal centrality metric (averaged across participants and data sessions) and their corresponding ICC values under each analytical strategy. Generally, significantly positive correlations were observed that were largely reproducible regardless of the analytical strategies and nodal centrality metrics used (Fig. 13). These results indicate higher reliability for more central nodes in morphological brain networks.

Figure 13

Relationship between nodal centralities and their reliabilities. In most cases, significantly positive correlations (solid circles) were observed between nodal centralities and their corresponding ICC values for each nodal centrality metric in particular for nodal efficiency. ICC, intra‐class correlation; Bin, binary network; Wei, weighted network.

Discussion

In this study, we constructed KLS‐based, individual‐level, whole‐brain morphological brain networks using structural MRI data and systematically investigated their topological organization under different analytical strategies. We found that the morphological brain networks were specifically organized in a high‐efficient, small‐world and modular manner with several highly connected hubs. Moreover, we demonstrated that these configurations of morphological brain networks were dependent on different choices of spatial smoothing, brain parcellation and network type. Further examination of long‐term TRT reliability showed that all the topological properties had fair to excellent reliability but spatial smoothing significantly improved reliability and nodal degree and efficiency outperformed nodal betweenness. Interestingly, nodal centralities were positively correlated with their ICC values, suggesting higher reliability for more central nodes. Taken together, these findings suggest that KLS‐based, single‐subject morphological brain network is a meaningful and reliable method in characterizing the brain organization and thus opens a new avenue toward understanding structural substrate of intersubject variability in behavior and function.

Currently, morphological brain networks are mainly derived by estimating interregional correlations in morphological features (e.g., cortical thickness or GM volume) at either group (He et al. 2007; Bassett et al. 2008) or individual (Tijms et al. 2012; Batalle et al. 2013) level. These correlation‐based methods require the normal distribution of morphological features across subjects or spatial locations or voxels. In contrast, the KLS‐based method used here has no such restriction, therefore may be more suitable for studying the human brain particularly given its complex folding structure. Moreover, based on empirical data, several recent studies have demonstrated the effectiveness of the KL/KLS method in studying the human brain in different populations (Velazquez and Galan 2013; Kong et al. 2014, 2015). Additionally, compared with a recently proposed single‐subject morphological network method (Tijms et al. 2012), the KLS‐based method has distinct advantages in network node definition which allows flexible choices of brain parcellation without any restriction on regional size and shape. This allows researchers to freely generate study‐specific, customized brain atlas according to their study objectives. All these features suggest that KLS‐based morphological network analysis could serve as a promising method to study morphological organization of the human brain at an individual level.

Specifically organized, analytical strategy‐dependent and long‐term reliable morphological brain networks

We showed that the KLS‐based, single‐subject morphological brain networks globally exhibited high‐efficient, small‐world and modular architecture. These findings are consistent with previous morphological brain network studies (He et al. 2007; Bassett et al. 2008; Chen et al. 2008; Tijms et al. 2012). Currently, the human brain is universally believed to have evolved to support both specialized or modular processing in local regions and distributed or integrated processing over the entire brain (Sporns et al. 2004; Bullmore and Sporns 2009, 2012; Meunier et al. 2009; He and Evans 2010). Here, our findings provide further empirical evidence to support the theory that the human cortical morphology has evolved into a complex but efficient neuronal architecture which confers an optimal balance between local specialization and global integration to maximize the power of parallel information processing. Locally, nodal centrality analysis of morphological brain networks revealed a heterogeneous distribution over the brain with several association and paralimbic regions (e.g., the precuneus, angular gyrus, cingulate gyrus, and parahippocampal gyrus) topologically holding central positions. These regions are largely comparable with the putative hubs reported in previous morphological, structural and functional brain networks (Achard et al. 2006; He et al. 2007, 2009b; Hagmann et al. 2008; Buckner et al. 2009; Gong et al. 2009; Tomasi and Volkow 2010). Recently, identifying brain hubs and further studying their vulnerability in various brain disorders are becoming a hot research topic (van den Heuvel and Sporns 2013; Crossley et al. 2014). Therefore, the current work provides an alternative method for researchers to explore hubs of the brain under both healthy and pathological conditions. Overall, these findings suggest that KLS‐based, single‐subject morphological brain networks are wired in an organized manner rather than randomly connected. Particularly, our global and local findings are largely consistent with those reported in a recent study which also investigated the topological organization of individual‐level morphological brain networks based on the KLS‐based method (Kong et al. 2015). These findings suggest a stable intrinsic architecture of KLS‐based morphological brain networks.

Nonetheless, compared with what was done by Kong and colleagues, the current work addressed several key questions which are important before the practical application of the KLS‐based approach. We first examined long‐term TRT reliability of KLS‐based morphological brain networks based on datasets collected with an interval of approximate 6 weeks. This is in contrast with what was done by Kong et al. (2015) who used datasets acquired on the same day to examine short‐term reliability. Long‐term reliability analysis is particularly important for newly developed methods before their clinical applications to identify stable biomarkers. We showed that most global and nodal network measures exhibited fair to excellent long‐term TRT reliability. Together with the high short‐term TRT reliability as reported in Kong et al. (2015), we propose the KLS‐based single‐subject morphological brain network analysis as a stable approach for future research to study effects of interest, such as topological reorganization in various brain disorders. Furthermore, we evaluated the effects of several factors on topological organization of KLS‐based morphological brain networks, including network type, brain parcellation and spatial smoothing. Insights into these issues are vital for researchers to determine their analytical strategies when using this approach. We showed that quantitative characterization of KLS‐based morphological brain networks was significantly affected by different analytical strategies. For example, spatial smoothing remarkably modulated spatial distribution of nodal centrality. The dependence of morphological brain networks on choices of preprocessing and network construction methods is consistent with previous brain network studies based on other neuroimaging modalities (Wang et al. 2009; Cole et al. 2010; Zalesky et al. 2010; Liang et al. 2012). Thus, the current results together with previous findings collectively highlight that for imaging connectomics studies, researchers should carefully choose their analytical pipelines and explain the findings when results are compared across studies with different processing methods. Apart from numerical values, the TRT reliability analysis further showed that different analytical strategies also affected the stability of KLS‐based morphological brain networks, consistent with previous brain network studies (Bassett et al. 2011; Wang et al. 2011; Liang et al. 2012; Buchanan et al. 2014; Zhao et al. 2015). Particularly, performing spatial smoothing significantly improved the reliability of both global and nodal network properties. This could be due to the elevated correspondence of brain structures among participants. To test this interpretation, we computed Pearson correlations in the KLS matrices across matrix elements for any pair of participants before and after spatial smoothing and found significantly higher intersubject similarities when spatial smoothing was conducted (P < 0.001 for both AAL and HOA brain parcellations). In addition, we found that nodal degree and efficiency outperformed nodal betweenness with respect to reliability. The higher reliability of nodal degree is consistent with our previous functional brain network study (Wang et al. 2011), possibly due to its conciseness and directness in definition. These findings provide important guidance on how to choose reliable analytical strategies and network metrics for KLS‐based, single‐subject morphological brain networks. Finally, we noted that several association and paralimbic regions were consistently identified as hubs (e.g., the precuneus, angular gyrus, parahippocampal gyrus, temporal and lateral occipital cortex) and consistently exhibited high reliability. Further correlation analysis revealed positive correlations between nodal centralities and their ICC values, suggesting higher reliability for more central regions of morphological brain networks. This is consistent with our previous finding of functional brain networks (Wang et al. 2011). These findings collectively imply that hub architecture is a stable organizational principle for the human brain networks (Buckner et al. 2009). Given the fact that brain hubs are generally implicated in various brain disorders (Crossley et al. 2014; Stam 2014), the high reliability of hubs makes them potential candidates to serve as reliable markers for disease diagnosis and prognosis, an interesting research topic in the future.

Taken together, the KLS‐based, single‐subject morphological brain networks are specifically organized, analytical strategy sensitive and TRT reliable, therefore opening a new avenue in linking morphological network variability and interindividual differences in behavior and cognition, although more studies are needed to elucidate the underlying biological significance.

Possible biological interpretation of KLS‐based morphological brain networks

In the current study, we calculated interregional similarity (as quantified by the KLS) in their distributions of GM volume to define morphological connectivity. However, the biological meaning underlying the similarity is not clear (Kong et al. 2014). Nevertheless, it should be noted that interregional covariance in morphological features has been observed as early as 1997 for several components of the human visual system (Andrews et al. 1997). Not limited to the visual system, the morphological coordination is demonstrated to expand to the whole brain, forming a morphological covariance network (He et al. 2007; Bassett et al. 2008). Moreover, an increasing number of studies have shown that the morphological covariance networks exhibit adaptive reorganization during normal development (Zielinski et al. 2010; Fan et al. 2011; Alexander‐Bloch et al. 2013b) and aging (Chen et al. 2011; Wu et al. 2012; Zhu et al. 2012) and in various brain disorders (He et al. 2008, 2009a; Seeley et al. 2009; Zhang et al. 2012). These studies jointly suggest that morphological covariance networks are biological meaningful in capturing potential mechanisms involved in these processes. Although the biological significance of morphological covariance is still not fully understood, accumulating evidence indicates that heredity, experience‐related plasticity, mutually trophic influences or coordinated neurodevelopment and aging trajectories play important roles in the formation morphological brain networks [for recent reviews, see (Alexander‐Bloch et al. 2013a; Evans 2013)]. In addition, there is another possible explanation coming from the axonal extension theory (Van Essen 1997), which suggests that connected areas tend to be pulled together by a tension from the axons between them. Therefore, it is plausible to speculate that these factors may also contribute to KLS‐based, individual‐level morphological brain networks studied here. Insights into this speculation could benefit from future studies by linking morphological brain organization with behavior performances or cognitive abilities for the same cohort of participants, or examining genetic, developmental, training‐induced, or brain disorder‐related changes of KLS‐based single‐subject morphological brain networks. Of note, the KLS‐based approach has been demonstrated to be capable of capturing age‐related changes in morphological brain networks (Kong et al. 2015).

Limitations and Future Directions

First, this study only examined several most prevalent topological attributes for single‐subject morphological brain networks. Besides these features, there are many other topological properties that are consistently observed in human brain networks, such as the “rich‐club” organization (van den Heuvel and Sporns 2011), hierarchy (Bassett et al. 2008) and heavy‐tailed degree distribution (Newman 2003; Avena‐Koenigsberger et al. 2015; Roberts et al. 2015). Future studies are needed to determine whether these configurations hold in KLS‐based, single‐subject morphological brain networks. Second, in the current study, we employed the AAL and HOA, two widely used atlases in previous brain network studies, to examine the effects of different brain parcellation schemes on topological organization of KLS‐based individual‐level morphological brain networks. Nevertheless, it should be noted that how to divide the brain into different ROIs to define network nodes is still an open question (de Reus and van den Heuvel 2013). Apart from the two atlases studied here, there are many other brain atlases available, such as the Automated Non‐linear Image Matching and Anatomical Labeling algorithm (Collins et al. 1995) and the LONI Probabilistic Brain Atlas (LPBA40) (Shattuck et al. 2008). Particularly, the Freesurfer provides a classification technique for automatically labeling individual brains into different regions that is robust to intersubject anatomical variability (Fischl et al. 2004; Desikan et al. 2006). In the future, it is important to compare these atlases to provide more comprehensive insights into how different parcellation schemes affect topological organization of individual‐level morphological brain networks. In addition, most of the current brain atlases are comparable with respect to the number of brain regions (i.e., similar spatial scale). Future studies are also needed to investigate how individual‐level morphological brain networks topologically organize at different scales by employing atlases that span several orders of magnitude in the number of regions or using an iterative random parcellation method (Fornito et al. 2010; Zalesky et al. 2010). Third, based on data from healthy participants we demonstrated that the KLS‐based method was reproducible and reliable in characterizing single‐subject morphological brain networks. The next step is to examine whether this method could reveal sensitive and reliable biomarkers associated with different states and various brain disorders. Third, single‐subject morphological brain networks were constructed by calculating interregional similarities in regional GM volume in the current study. Straightforwardly, this method could be extended to other morphological features (e.g., cortical thickness) or even images of other modalities (e.g., positron emission tomography). Furthermore, in addition to the KLS measure used here, there are also other measures to quantify the similarity of two curves, such as the Bhattacharyya distance (Zhou and Chellappa 2006) or related measures (De Maesschalck et al. 2000; Comaniciu et al. 2003). Future studies are required to examine unique insights into morphological brain networks from different choices of these factors. Finally, an increasing number of studies have examined the relationship between anatomical and functional brain networks and find that functional connectivity profiles are largely shaped but not fully determined by structural pathways (Wang et al. 2015b). In this regard, it is an interesting topic for future research to determine the similarities and differences between single‐subject morphological brain networks and those derived from other neuroimaging modalities in mapping and characterizing the human connectome.

Conclusion

This study demonstrates that KLS‐based, single‐subject morphological brain networks are specifically organized with several nontrivial topological features, which are TRT reliable but depend on choices of analytical strategies. This method therefore could complement the current methodology of neuroimaging connectomics and open a new avenue toward understanding structural substrate of intersubject variability in behavior and function from a network perspective.

Conflict of Interest

The authors declare that there are no conflicts of interest.

Table S1. Regions of interest from the AAL atlas.

Table S2. Regions of interest from the HOA atlas.

Click here for additional data file.