Mapping individual voxel-wise morphological connectivity using wavelet transform of voxel-based morphology

Mapping individual brain networks has drawn significant research interest in recent years. Most individual brain networks developed to date have been based on fMRI or diffusion MRI. Given recent concerns regarding confounding artifacts, various preprocessing steps are generally included in functional or structural brain networks. Notably, voxel-based morphometry (VBM) derived from anatomical MRI exhibits high signal-to-noise ratios and has been applied to individual interregional morphological networks. To the best of our knowledge, individual voxel-wise morphological networks remain unexplored. The goal of this research is twofold: to build novel metrics for individual voxel-wise morphological networks and to test the reliability of the proposed morphological connectivity. To this end, anatomical scans of a cohort of healthy subjects were obtained from a public database. The anatomical datasets were preprocessed and normalized to the standard brain space. For each individual, wavelet-transform was applied on the VBM measures to obtain voxel-wise hierarchical features. The voxel-wise morphological connectivity was computed based on the wavelet features. Reliable brain hubs were detected by the z-scores of degree centrality. High reliability was discovered by test-retest analysis. No effects of wavelet scale, network threshold or network type were found on hubs of group-level degree centrality. However, significant effects of wavelet scale, network threshold and network type were found on individual-level degree centrality. Significant effects of network threshold and network type were found on reliability of degree centrality. The results suggested that the voxel-wise morphological connectivity was reliable and exhibited a hub structure. Moreover, the voxel-wise morphological connectivity could reflect individual differences. In summary, individual voxel-wise wavelet-based features can probe morphological connectivity and may be beneficial for investigating the brain morphological connectomes.

Introduction Voxel-based morphometry (VBM) is a neuroimaging tool that can probe local differences of brain anatomy and has been applied in statistical analysis [1,2] and data mining [3][4][5]. The brain is a complex network and contains high-order cortical information, although the local features of the brain morphology do not reflect the brain topology. Individual brain topology has been widely investigated by diffusion MRI and functional MRI [6,7]. Meaningful network topology properties (i.e., hubs, scale-free, small-worldness) have been found in the human brain [8][9][10]. Given recent concerns regarding confounding artifacts, various preprocessing steps are generally implemented in functional or structural brain networks [11,12]. Notably, VBM features that exhibit high signal-to-noise ratios can be derived from anatomical MRI with well-established preprocessing procedures. Hence, investigating the individual morphological relationships based on VBM features may be beneficial for studies of the brain connectome.
Mapping individual morphological relationships is a challenging task based on local features of an anatomical volume. Interregional image density distribution-based Kullback-Leibler (KL) divergence has been recently applied to measure individual morphological connectivity [13]. The interregional KL similarity was shown to be reliable [14] and related to brain maturation [15]. Hence, it is feasible to map individual morphological connectivity using interregional measures. Although kernel density estimation exhibited good performance across networks of different resolutions [14], the different probability distribution functions might influence the morphological connectivity [15]. Thus far, biological explanations of density distribution-based similarity remain unclear [13]. Moreover, the aforementioned studies relied on regional statistics, and therefore, the results cannot be extended to investigations of individual voxel-wise morphological connectivity.
Mapping inter-voxel morphological relationships at the single-subject level remains largely underexplored, as there is a lack of local features in a single voxel of an anatomical volume. Notably, voxel patch-based similarity has been applied to extract individual brain networks from structural MRI scans [16], which has shed light on voxel-wise anatomical networks. To obtain morphological relationships, the brain was divided into thousands of cubes, which served as network nodes. The inter-cube morphological relationship was defined as the correlation between the two cubes of 27 voxels [16]. However, the cubes could not represent the complex structure of the brain, which exhibited variability across subjects. Moreover, the voxel patch could only capture the local features within a cube and did not reflect the global features of functionally or structurally homogeneous brain regions [16]. Thus, the cubes of voxel patch features should be replaced by a vector of hierarchical features that contains both local and global features in a single voxel.
Wavelet transform has been applied in brain morphometry studies [17]. The wavelet transform is a multiscale analysis method that can transform the energy of a signal into a multi-resolution hierarchical organization. Hence, wavelet transform can capture both local and global features of an anatomical MRI dataset. Wavelet-based morphometry (WBM) is a novel method that can probe structural morphometric differences [17]. Moreover, the wavelet transform is a multivariate method that has been applied in diagnostic models of brain disorders [18]. Diagnostic models based on wavelet transform exhibit higher sensitivity and specificity than conventional methods [18]. Therefore, wavelet transform could be beneficial for the analysis of anatomical MRI. Considering the hierarchical organization of wavelet features, investigating the voxel-wise relationships using wavelet transform is a straightforward approach.
In this paper, we aimed to map the individual voxel-wise morphological connectivity based on anatomical MRI and to test the reliability and variability of the voxel-wise morphological connectivity. To this end, a cohort of adult participants was obtained from a publicly available MRI database. The method section describes the application of the wavelet transform to the VBM data to obtain multi-resolution features for a single volume. Wavelet-based similarity was applied as the voxel-wise morphological connectivity. In the result section, the reliability and variability of the voxel-wise morphological connectivity were evaluated using statistical analysis. In the discussion section, we summarize the reliability and variability of the voxelwise morphological networks.

Subjects and MRI protocols
Twenty-one healthy subjects (11 males and 10 females, 22-61 years old) without neurological history were recruited at Vanderbilt University, USA [19]. For each participant, one structural MRI dataset was obtained using a 3T Philips Achieva MR scanner. This study was approved by the institutional review boards of Vanderbilt University, USA. Signed informed consent was provided by each participant. Each subject was scanned twice with a one-hour break between scans. All scans were completed within 2 weeks. The structural MRI dataset was acquired with a magnetization-prepared rapid acquisition gradient echo (MPRAGE) sequence (repetition time = 6.7 ms, echo time = 3.1 ms, inversion time = 842 ms, 1 mm×1 mm×1.2 mm, 6 minutes). The datasets are publicly available without restrictions for academic usage (https://www.nitrc. org/projects/multimodal/).

Data preprocessing
The structural MRI datasets were first reoriented to the standard space, automatic cropped, bias field corrected, and skull-stripped to exclude non-brain tissues. The brain images were then segmented into gray matter, white matter, and cerebrospinal fluid. Third, the original and segmented images were spatially normalized to the standard brain space. Fourth, the spatial normalized gray matter image was multiplied by the Jacobian of the warp field. Finally, the modulated gray matter image was obtained as the VBM features. The spatial normalization procedure was carried out using the fsl_anat command. The VBM procedure was carried out using the fslvbm command.

Wavelet transform
The wavelet transform was applied to extract local and global features for the VBM volumes, resulting in hierarchical features for each single voxel. Thus, the wavelet transform could provide richer information than the conventional VBM method. A wavelet can be defined as a small wave-like oscillation that can be used to decompose a time-domain or spatial-domain signal into different scales [18]. The wavelet transform first decomposes the original 3D images into sub-bands at different spatial scales and then adaptively concentrates the original 3D spatial information into a sets of wavelet coefficients [17]. Here, the multi-resolution analysis was carried out through a discrete orthogonal wavelet transform, which adaptively decomposed the images at different spatial-scale multi-frequency bands based on a small series of wavelet coefficients [17]. The discrete orthogonal wavelet transform was successfully applied in wavelet-based morphometry. For each subject, discrete orthogonal wavelet transform-based feature extraction was conducted using the following steps: 1) The wavedec3 function in MATLAB is used for multilevel 3D wavelet decomposition of the VBM volume; 2) The coefficients of the wavelet components are obtained and saved; 3) The waverec3 function in MATLAB is used to obtain the 3D wavelet reconstruction of approximations and details; 4) Each wavelet component is normalized according to the z-scores; 5) The approximations and details are concentrated into a 4D volume. Here, the approximations represent low-pass components of the wavelet transform, while the details indicate high-pass components of the wavelet transform. Specifically, the approximations denote global features that change over the long space, while the details denote local features that change over the short space.

Voxel-wise morphological connectivity
The nodes and edges of brain networks should be defined before computing voxel-wise morphological relationships. The locations of nodes were defined using the first 90 brain regions in the automated anatomical labeling (AAL) atlas, resulting in 159841 nodes in the 2 mm×2 mm×2 mm 3D space. Given a voxel-wise morphological network, the morphological connectivity (edge) between two voxels (nodes) can be defined by the Pearson correlation coefficient of pairwise wavelet feature vectors in brain voxels.
Here, w denotes the wavelet feature vector within a voxel, cov denotes covariance, and δ denotes the standard deviation of the wavelet feature vector. For a wavelet scale of n, the length of the wavelet feature vector is 2 Ã n.
For a binarized graph, the links are defined as: ( For a weighted graph, the links are defined as: Here, i and j denote the two voxels, and r 0 denotes the threshold for a graph [20]. The thresholds r 0 for morphological connectivity were set to 0.5, 0.6, 0.7, 0.8 and 0.9 for sparsity. Given a wavelet-transformed 4D VBM volume for each individual, voxel-wise morphological relationships were computed by the following steps: 1) The voxel-wise wavelet features are obtained; 2) The voxel-wise features are masked using the AAL atlas; 3) the Pearson correlation coefficients of pairwise wavelet features between the two voxels are derived; 4) The voxelwise morphological connectivity is obtained; 5) Thresholds (i.e.,0.5, 0.6, 0.7, 0.8, 0.9) are applied to the voxel-wise morphological connectivity; 6) The voxel-wise nodal degree of morphological connectivity is obtained. The pipeline to obtain individual voxel-wise morphological connectivity is illustrated in Fig 1. The scripts can be obtained from the Neuroimaging Informatics Tools and Resources Clearinghouse (NITRC) website (https://www.nitrc.org/ projects/mcwt/).

Degree of centrality
Centralities are the most common graph measures for nodal estimators (i.e., degree centrality, betweenness centrality, eigenvector centrality) [7,[20][21][22]. Degree centrality is one of the most common measures of centrality and is defined as the number of links (neighbors) connected to a node [7]. In this paper, degree centrality was applied to investigate the nodal properties of the morphological networks. Degree centrality was calculated using a publicly available script in the Data Processing Assistant for Resting-State fMRI (DPARSF) toolbox [23].
Both binarized and weighted graphs were analyzed through thresholding the symmetrical matrix of morphological connectivity [20]. Degree centrality for binarized networks was derived by computing the sum of the number of edges within the binarized graph. Degree centrality for weighted networks was derived by computing the sum of the values of the weighted edges within the weighted graph. Before statistical analysis, the individual maps of degree centrality were spatially smoothed using the fslmaths function with a Gaussian kernel (sigma = 3 mm). To analyze the effects of wavelet scale (i.e., 3, 4, 5), network threshold (i.e., 0.5, 0.6, 0.7, 0.8, 0.9), and network type (i.e., binarized or weighted), the raw scores of degree centrality were transformed into z-scores using the following formula: where D i indicates the ith voxel-wise degree centrality, D indicates the whole-brain voxel-wise degree centrality, and std indicates the standard deviation.

Hub detection
Hub voxels contain more connections than other voxels within a voxel-wise network. To detect the hubs of the voxel-wise networks, the raw scores of degree centrality were first transformed into z-scores. Then, group-level z-scores of degree centrality were obtained by averaging the individual degree centrality with different effects (i.e., wavelet scales, network thresholds, types of network). Finally, for each effect factor, the hub voxels were indicated by those with a z-score >1, according to previous studies [24][25][26].

Test-retest analysis
A test-retest analysis was applied to investigate the reliability of the proposed voxel-wise morphological network. The test-retest reliability was analyzed via the intraclass correlation coefficients (ICCs) [27,28]. In this paper, the ICCs were computed based on a random effects model and were given by the ratio of the intersubject variance to the total variance. According to previous studies [29], the ICCs were partitioned into five levels: excellent reliability (ICC>0.8), high reliability (0.6<ICC<0.79), moderate reliability (0.4<ICC<0.59), fair reliability (0.2<ICC<0.39), and poor reliability (ICC<0.2).

Statistical analysis
To analyze the effects of wavelet scale (i.e., 3, 4, 5), network threshold (i.e., 0.5, 0.6, 0.7, 0.8, 0.9), and network type (i.e., binarized or weighted) on the network hubs and network reliability, three-way Analysis of Variance (ANOVA) was performed on the group-level degree centrality or reliability maps using the following formula: Here, Y denotes group-level degree centrality or reliability maps, and A, B, and C denote wavelet scales, network thresholds, and network types, respectively.
To analyze the effects of wavelet scale, network threshold, and network type on individual differences in degree centrality, repeated-measures ANOVA was performed on the individuallevel degree centrality maps using the following formula: Here, DC denotes the individual-level degree centrality, and A, B, and C denote wavelet scales, network thresholds, and network types, respectively. The results were corrected by the false discovery rate (FDR). Table 1 shows the network sparsity based on different network thresholds. Table 2 and Table 3 present the results of the three-way ANOVA on group-level degree centrality among the different kinds of networks. There are no effects of wavelet scale, network threshold or network type on group-level degree centrality. Fig 2 shows the correlations of group-level degree centrality among different kinds of networks. The maps of degree centrality are significantly correlated with each other, and most of the correlation coefficients are close to 1. The correlation coefficients are also reliable with different sessions. The group-level degree centrality exhibits reliable spatial patterns across different wavelet scales, network thresholds and network types. The hub brain regions in the voxel-wise networks Fig 3 shows the brain hubs (i.e., precuneus, cingulate gyrus, medial frontal gyrus, superior frontal gyrus, superior temporal gyrus, middle temporal gyrus, middle occipital gyrus, cuneus, insula, hippocampus, and amygdala) of the voxel-wise networks. The brain hub regions in the voxel-wise networks exhibit reliable spatial patterns across the two scan sessions. Fig 3A shows the brain hubs based on the average networks of session 1, and Fig 3B shows the brain hubs based on the average networks of session 2. The brain hubs are indicated by a z-score >1, which represents the cores of connections. The average percentages of high degree centrality (i.e., z-score >1) for session 1 and session 2 are 14.16% and 14.25%, respectively. The 10 brain hub regions of session 1 are listed in Table 4, and the 11 brain hub regions of session 2 are listed in Table 5. The sizes and locations of the hub regions are comparable between the two sessions. Table 6 show the results of three-way ANOVA for reliability of degree centrality. There are no effects of wavelet scale on the reliability of voxel-wise degree centrality. However, the network threshold and network type might affect the reliability of voxel-wise degree centrality.  Table 7 shows the mean and standard deviation of reliability. Most of the mean values of ICCs exceed 0.7 with standard deviations of less than 0.14. High-to-excellent reliability of degree centrality can be observed in Fig 5 and  Table 7. Fig 6 shows the maps of reliability with a network threshold of 0.8. The voxel-wise degree centrality exhibits high test-retest reliability across different networks.

Individual-level degree centrality of the voxel-wise networks
The repeated-measures ANOVA analysis revealed significant (p<0.05, FDR corrected) main effects and interaction effects of the three different factors (i.e., threshold, scale, type) on individual-level degree centrality. A significant main effect of network threshold was found in 91.62% of the whole-brain voxels. A significant main effect of wavelet scale was found in 97.16% of the whole-brain voxels. A significant main effect of network type was found in 81.84% of the whole-brain voxels. A significant threshold:scale interaction was found in 99.1% of the whole-brain voxels. A significant threshold:type interaction was found in 88.24% of the whole-brain voxels. A significant scale:type interaction was found in 96.12% of the whole brain voxels. A significant threshold:scale:type interaction was found in 98.55% of the wholebrain voxels. Individual voxel-wise morphological connectivity using VBM

Discussion
In this study, we sought to map voxel-wise morphological connectivity at the single-subject level based on anatomical MRI. To this end, wavelet transform was applied on the VBM features. The morphological connectivity was obtained via the inter-voxel hierarchical wavelet features. The voxel-wise morphological networks exhibited good test-retest reliability across different wavelet scales, different network thresholds, and different network types. In addition, the voxel-wise morphological networks exhibited hub structures, most of which were located in the default mode network (DMN) area. Overall, the inter-voxel morphological relationships were successfully obtained for individuals through wavelet transform. The voxel-wise morphological connectivity might provide additional measures for studies of the human connectome.  Wavelet transform can capture the local and global features of individual MRI datasets in a multi-resolution manner. Thus, individual brain networks can be built through wavelet-based hierarchical features in a straightforward manner. Wavelet transform has been employed as an efficient tool for feature extraction and signal representation [17]. The components of wavelet transform might reflect the biological meanings of the original images. The low-frequency information decoded by wavelet transform might be beneficial in assisting in the diagnostic process of brain disorders [18]. Therefore, wavelet transform could provide rich information for disease classification and advance the analysis of neuroimaging data. Moreover, wavelet features can indicate inherent hierarchical structures, which could be used to probe multiscale information of brain imaging datasets. In contrast to previous studies, this study for the first time mapped the voxel-wise morphological networks, which correspond to multivariate measures. Furthermore, the results revealed that the wavelet scales had no effects on group-level degree centrality or its reliability, suggesting that wavelet transform could capture reliable features that could be beneficial to the brain connectome. However, the biological meanings of the morphological network remain unclear. Possible explanations can be drawn from the axon tension theory [13] or that morphologically connected brain regions may share similar structures for information transmission.
The hub regions, which are connected to more brain nodes than other regions, may reflect the functional or structural cores within the brain networks. In this study, we found certain brain regions that exhibited hub organization (i.e., the precuneus, the cingulate gyrus, the medial frontal gyrus, the superior frontal gyrus, and the superior temporal gyrus), most of Individual voxel-wise morphological connectivity using VBM which were located in the DMN area. DMN-related morphological hubs might reflect structural or functional connections. Pervious evidence has indicated structural cores located in brain regions of the posterior components of the DMN [21]. High global connectivity has been found in the DMN regions and cognitive control regions, suggesting the important roles of hub regions in information processing [30]. Moreover, the low percentage of network nodes with high degree centrality found in this study suggested that the morphological connectivity might exhibit scale-free organization. The human brain functional network has been shown to be regulated by highly connected cortical cores and exhibits power-law degree distribution [31]. Therefore, the brain hubs detected by morphological connectivity could be supplementary to conventional functional and structural cores. In summary, the morphological hubs might possess functional and structural meanings and might be related to brain disorders [32] and development [33]. The test-retest reliability is a key feature for evaluating novel neural metrics [34]. The active patterns (i.e., Amplitude of Low Frequency Fluctuations (ALFF), Regional Homogeneity (ReHo), entropy) of resting-state blood-oxygenation level-dependent (BOLD) signals were reliable at the voxel level [24,27,35]. Reliable voxel-wise measures of VBM were found using multicenter datasets [36]. In this paper, the voxel-wise degree centrality maps were derived from the VBM features. Thus, the proposed voxel-wise morphological networks should exhibit certain test-retest reliability. Reliable hub structures of voxel-wise morphological networks at the group level were found across the two scan sessions. High reliability of the voxel-wise morphological networks was confirmed by the test-retest analysis. The reproducible hubs suggested that wavelet-based morphological connectivity might be a reliable metric to probe voxel-wise brain connectivity. However, the network thresholds and network types might Binarized scale 3 0.71 ± 0.14 0.72 ± 0.14 0.73 ± 0.14 0.74 ± 0.14 0.74 ± 0.14 Binarized scale 4 0.71 ± 0.14 0.72 ± 0.14 0.73 ± 0.14 0.74 ± 0.14 0.74 ± 0.14 Binarized scale 5 0.71 ± 0.14 0.72 ± 0.14 0.73 ± 0.14 0.74 ± 0.14 0.74 ± 0.14 Weighted scale 3 0.72 ± 0.14 0.73 ± 0.14 0.73 ± 0.14 0.74 ± 0.14 0.74 ± 0.14 Weighted scale 4 0.72 ± 0.14 0.73 ± 0.14 0.73 ± 0.14 0.74 ± 0.14 0.74 ± 0.14 Weighted scale 5 0.72 ± 0.14 0.73 ± 0.14 0.73 ± 0.14 0.74 ± 0.14 0.74 ± 0.14 https://doi.org/10.1371/journal.pone.0201243.t007 Individual voxel-wise morphological connectivity using VBM impact the test-retest reliability of voxel-wise networks. The results suggested that the network thresholds and network types should be carefully selected to verify the test-retest reliability of morphological networks. Nevertheless, high ICC values were obtained across different wavelet scales, network thresholds, and network types, implying that the morphological connectivity is a reliable measure for clinical applications. Individual differences are a fundamental property of a novel neural metric. Previous individual connectivity was based on diffusion MRI or functional MRI, which are both sensitive to confounding artifacts [13]. Given the well-established scanning and preprocessing procedures, anatomical MRI could be beneficial for understanding individual variability. Notably, previous structural covariance networks based on group-level measures could not represent individual differences [37,38]. Source-based morphometry (SBM) is a novel method based on grouplevel independent component analysis (ICA) of anatomical MRI data [39]. The method proposed in this study could be extended to individual-level SBM, as a wavelet-based 4D volume was obtained for each subject. Thus, individual differences of spatial sources based on SBM could be investigated using wavelet features. In this study, repeated-measures ANOVA revealed significant individual differences in the morphological connectivity. Significant main and interaction effects of the three factors (i.e., network threshold, wavelet scale, network type) were found throughout the whole brain, suggesting that the network factors could influence individual differences. Therefore, network factors should be taken into account when studying individual differences in the proposed morphological connectivity. To the best of our knowledge, this is the first study to map individual voxel-wise morphological connectivity based on anatomical MRI. The proposed morphological connectivity could advance the statistical analysis of individual differences and machine learning-based diagnostic models.
This study is characterized by several advantages. First, the voxel-wise morphological network was investigated for single subjects in a wavelet-based framework. Evaluating the voxelwise morphological connectivity for single subjects might be beneficial for investigating individual differences. Second, the wavelet transform could yield hierarchical features (i.e., global features, local features). Thus, the wavelet-based morphological connectivity could capture multi-resolution information. Third, the voxel-wise morphological connectivity detected reliable hub structures. The reliable hub structures identified by the framework might reflect the functional or structural basis of voxel-wise morphological connectivity. Fourth, the voxel-wise morphological estimators exhibited good test-retest reliability. Reliability is an essential feature for novel neural metrics. The reliable voxel-wise morphological estimators implied their potential values in clinical applications. In summary, the proposed morphological connectivity method was reliable across different network levels and could advance the study of individual variability.
One limitation of this study is that the biological meanings of the proposed voxel-wise morphological connectivity feature could not be elucidated. Although reliable voxel-wise morphological connectivity was discovered by wavelet-based features, supplemental experiments are still required to test the hypothesis that morphological connectivity can reflect structural and functional connectivity. Moreover, there are numerous similarity and dissimilarity metrics (i.e., KL divergence, correlative features) that can be applied to measure the morphological connectivity in additional studies. Another limitation of this study is the lack of physiological signals for references. The morphological similarity-based results could be interpreted with supplementary physiological measures (i.e., electroencephalography (EEG), magnetoencephalography (MEG)) in future studies. The third limitation of this study is the relatively small sample size of participants. The results should be verified using different scanners and different populations in subsequent studies. Together, there are still many challenges that must be overcome in the study of voxel-wise morphological networks. Additional physiological measures might help interpret the novel individual morphological networks.

Conclusion
This paper proposed a voxel-wise wavelet-based similarity measure to evaluate morphological connectivity at the single-subject level based on structural MRI. The morphological networks were reliable across different wavelet scales, network thresholds, and network types. The voxel-wise morphological connectivity exhibited excellent reliability and reflected individual differences. In summary, the voxel-wise wavelet-based similarity could probe individual morphological connectivity and could be beneficial for investigating the brain morphological connectome.