Extended multimodal whole-brain anatomical covariance analysis: detection of disrupted correlation networks related to amyloid deposition.

BACKGROUND: An anatomical covariance analysis (ACA) enables to elucidate inter-regional connections on a group basis, but little is known about the connections among white matter structures or among gray and white matter structures. Effect of including multiple magnetic resonance imaging (MRI) modalities into ACA framework in detecting white-to-white or gray-to-white connections is yet to be investigated. NEW
METHOD: Proposed extended anatomical covariance analysis (eACA), analyzes correlations among gray and white matter structures (multi-structural) in various types of imaging modalities (T1-weighted images, T2 maps obtained from dual-echo sequences, and diffusion tensor images (DTI)). To demonstrate the capability to detect a disruption of the correlation network affected by pathology, we applied the eACA to two groups of cognitively-normal elderly individuals, one with (PiB+) and one without (PiB-) amyloid deposition in their brains.
RESULTS: The volume of each anatomical structure was symmetric and functionally related structures formed a cluster. The pseudo-T2 value was highly homogeneous across the entire cortex in the PiB- group, while a number of physiological correlations were altered in the PiB + group. The DTI demonstrated unique correlation network among structures within the same phylogenetic portions of the brain that were altered in the PiB + group. COMPARISON WITH EXISTING
METHOD: The proposed eACA expands the concept of existing ACA to the connections among the white matter structures. The extension to other image modalities expands the way in which connectivity may be detected.
CONCLUSION: The eACA has potential to evaluate alterations of the anatomical network related to pathological processes.

Introduction

The brain is an organ that consists of basic functional elements, wiring that connects the units to form purposeful and efficient functional modules, and systems to support and maintain the working environment for cells in the brain. The gray matter forms the functional units and the white matter fibers, as wires, connect these units. Advancements in neuroimaging technologies and neurophysiology allow the in vivo observation of the connections among gray matter structures, which has evolved into an emerging concept of the connectome. Among various methods to infer connections among brain regions, anatomical covariance analysis (ACA) (Alexander-Bloch et al., 2013; Evans, 2013; Lerch et al., 2006) has been emerged to infer inter-regional connections on a group basis. The ACA is based on an assumption that anatomical structures that form a functional unit are tightly connected with axons and synapses, and therefore, their biological properties are similar. For example, if region A and B are connected, a person with high gray matter density in the region A also has the region B with high gray matter density (and vice versa). Such relationships, termed correlation networks, might be genetically defined, or could be acquired later in life through experience, training, or degeneration caused by various diseases. One of the well-described disease-related anatomical correlations is diaschisis, which was originally defined as a dysfunction of areas that are connected to a damaged brain area, but can also result in atrophy of the areas connected to the damaged area (e.g., crossed cerebellar atrophy (Jimenez-Caballero, 2012; Menon, 2012; Olsen et al., 1963)). For the past decade, efforts have been made to systematically and comprehensively investigate the correlations of cortical thickness, volume, density, or myelination (Ma and Zhang, 2017) among entire cortical areas, which was summarized in (Evans, 2013). The pattern of anatomical correlation is similar to that observed in resting-state fMRI analysis, but varies depending on the stage of brain development (Raznahan et al., 2011) and aging (Brickman et al., 2007). Disease-related anatomical covariance has been elucidated in various neurological diseases with brain atrophy, including Alzheimer's disease (He et al., 2008), frontotemporal dementia (Hafkemeijer et al., 2016), semantic dementia (Seeley et al., 2009), progressive non-fluent aphasia (Seeley et al., 2009), corticobasal syndrome (Seeley et al., 2009), multiple sclerosis (He et al., 2009), and temporal lobe epilepsy (Bernhardt et al., 2008), as well as disorders with less morphological abnormalities, such as schizophrenia (Wible et al., 1995, 2001) or fibromyalgia (Kim et al., 2015).

To systemically and comprehensively apply the ACA strategy to capture the anatomical features of the brain, several attempts have been made to analyze correlation network among white matter as well as gray matter structures. Such methods utilized the white matter parcellation based on a projection of cortical areas (Mitelman et al., 2005) or white matter regions-of-interest derived from voxel-based group comparison (Hulshoff Pol et al., 2004). However, existing approaches cover limited amount of preselected white matter areas or ignore boundary of the white matter structures that should be defined by orientation of fiber bundles. Therefore, structure-by-structure correlation network of the white matter area still need to be elucidated.

In this study, we extended the ACA method in two directions: multi-structural and multi-modality. Our approach, the extended anatomical covariance analysis (eACA) method, analyzed correlations among whole brain areas, including gray and white matter structures (multi-structural) in various types of imaging modalities (multi-modality). To parcellate both gray and white matter structures, we applied the multi-atlas label fusion method (Wu et al., 2016). Following this multi-dimensional extension in the image domain, the biological rationale for the “connectome” seen in the eACA needs to be re-defined. Connectome originally implies the connection among or within gray matter structures where neurons exist. The spatial extension to the white matter structures expands the concept to the connections between the white matter structures and between the gray and the white matter structures, in addition to the conventional gray-to-gray connections. The extension to multiple image modalities expands the way in which connectivity may be detected, which was originally based only on the regional property of the cortical gray morphology measured from T1-weighted images. We included diffusion MRI-based measures, fractional anisotropy (FA), and the sum of the three eigenvalues (trace) as markers of white matter integrity or the rotationally invariant magnitude of water diffusion within brain tissue, both of which are known to be affected by brain development (Dubois et al., 2008; Faria et al., 2010; Huppi and Dubois, 2006; Schmithorst et al., 2002, 2005), aging (Sexton et al., 2014; Xie et al., 2016) or diseases (for reviews (Dong et al., 2004; Lerner et al., 2014; Rollins, 2007; Sundgren et al., 2004). An estimate of T2 relaxation time (pseudo-T2 or pT2, as referred to in (Nakamura et al., 2014)), was included as a measure by which to evaluate the status of myelination and macromolecules (Du et al., 2012; Miot-Noirault et al., 1997; Mosher and Dardzinski, 2004).

We hypothesized that the re-defined connectome in the eACA could be applied to capture the alterations in the brain correlation network, especially related to brain pathologies seen in both gray and white matter structures. To test this hypothesis, we investigated whether the eACA could detect correlation networks related to amyloid deposition in cognitively normal individuals. Increasing evidence indicates that amyloid deposition is seen in 20–40% of cognitively normal individuals over 65 years of age (Aizenstein et al., 2008; Bouwman et al., 2009; Mintun et al., 2006; Peskind et al., 2006; Shaw et al., 2009). Amyloid pathology is known to alter both gray (Chang et al., 2015; Ten Kate et al., 2018; Tijms et al., 2016) and white matter structures (Kantarci et al., 2014; Molinuevo et al., 2014) through the neuronal toxicity (Bloom, 2014; Iadanza et al., 2018) and vascular changes (Dallaire-Theroux et al., 2017; Lee et al., 2018). Recent studies suggest that whether subtle anatomical alterations related to amyloid deposition are MRI-detectable or not depends on the demographics of the target population and methods applied to analyze MRIs. The previous findings include atrophy located in the cingulate cortex (Storandt et al., 2009), the temporal neocortex (Chetelat et al., 2012), the hippocampal areas (Hsu et al., 2015), the inferior frontal cortex (Oh et al., 2014), and the basal forebrain (Grothe et al., 2014), or no atrophy (Bourgeat et al., 2010). Based on our hypothesis that the eACA could provide unique anatomical features that are difficult to capture through conventional regional comparison, we investigated the correlation network of the amyloid-positive and amyloid-negative normal cognitive groups to determine whether the eACA could detect differences between the two groups.

Methods

Participants

All participants enrolled in this study were selected from the BIOCARD study (Biomarkers of Cognitive Decline Among Normal Individuals (http://www.biocard-se.org) (Albert et al., 2014), which aims to identify biomarkers associated with the progression from normal cognitive status to cognitive impairment or dementia, with a particular focus on Alzheimer's disease. These analyses are based on a subset of participants who have been scanned with 3T neuroimaging protocols since 2015, including high-resolution anatomical MRI and amyloid positron emission tomography (PET). In these analyses, 95 participants were scanned with multi-modal MR imaging (T1-and T2-weighted images, and diffusion tensor images), and the Pittsburgh compound, B (PiB-PET)). Eight subjects were excluded because of a voxel size coding error or multi-modal mis-coregistration that caused miscalculation of the distribution volume ratio (DVR) map of PiB-PET scanning. The remaining 87 subjects were selected for the subsequent processing. Table 1 summarizes the demographic information of the participants involved in this study. The research protocols of the BIOCARD study were approved by local institutional review boards. All subjects signed written, informed consent for their participation.

Table 1

Characteristics of subjects.

Characteristric	PiB- (n = 63)	PiB+ (n = 24)	Statistic	p value
Age (y)	66.5 ± 10.7	70.4 ± 5.4	2.23†	0.029*
Sex (female)	42 (67%)	15 (63%)	0.013§	0.910
MMSE1	29.5 ± 0.7	29.2 ± 2.1	0.713†	0.483
Total CVLT2	65.3 ± 9.2	61.4 ± 14.2	1.26†	0.218
cDVR3	0.95 ± 0.04	1.23 ± 0.16	8.373†	<0.001***

* p < 0.05 and *** p < 0.001.

MMSE = Mini-Mental State Examination.

CVLT = California Verbal Learning Test.

cDVR = mean cortical distribution volume ratio.

statistic in Welch's two-sample t-test.

χ2 in chi-squared test.

Scan protocol

MRI scans were obtained on a Philips Achieva 3.0T scanner. T1-weighted structural images were obtained from the three-dimensional (3D) magnetization-prepared rapid gradient-echo (MPRAGE) sequence (repetition time [TR] = 6.8 ms, echo time [TE] = 3.1 ms, inversion time [TI] = 843 ms, field of view [FOV] = 256 × 256, thickness/gap = 1.2/1.2 mm, flip angle = 8, sagittal plane). Proton density and T2-weighted images were obtained using a double-echo fast spin echo sequence (TR = 4.2 s, TE1/TE2 = 28/80 ms, FOV = 212 × 212, thickness/gap = 2.2/2.2 mm, flip angle = 90, axial plane), from which T2 maps were generated. To obtain diffusion tensor images, b0 and diffusion weighted images were acquired from a spin echo sequence (TR = 7.5 s, TE = 75 ms, FOV = 260 × 260, thickness/gap = 2.2/2.2 mm, flip angle = 90, b-value = 700, number of gradients = 33, axial plane).

Dynamic 11C-PiB PET studies were performed on a GE scanner in 3D mode. To minimize motion during scanning, participants were fitted with a thermoplastic head-mask. After an IV bolus injection of a mean (SD) 14.7 (0.8) mCi of 11C-PiB, PET data were scanned immediately and acquired according to the protocol: 4×0.25; 8×0.5; 9×1, 2×3; and 10 × 5 min (70 min total, 33 frames). Dynamic images were reconstructed using filtered back-projection with a ramp filter (image size = 128 × 128, pixel size = 2 × 2 mm, slice thickness = 4.25 mm), yielding a spatial resolution of approximately 4.5 mm full width at half maximum (FWHM). Overall quality control review was performed by one of the co-authors (K.O.).

Image processing

The pipeline of the proposed eACA was summarized in Fig. 1. Details of each step will be described explicitly below.

Fig. 1

Schematic flowchart of the proposed extended anatomical covariance analysis. A: parcel-based metrics for each image modality were extracted by image preprocessing, coregistration and brain parcellation, PiB+/PiB- grouping was achieved based on cortical DVR images, B: parcel-by-parcel correlation matrices of PiB+ and PiB- group were generated for each image modality, C: hierarchical clustering was performed on each correlation matrix, D: significant between-group difference for each parcel-pair correlation was evaluated by permutation test. The x-axis is the absolute difference of inter-parcel correlation coefficients between groups.

T1-weighted images

To remove potential confounders of low-frequency-intensity non-uniformity present in the MR images, N4ITK bias correction was applied on the 3D T1-weighted images (Mintun et al., 2006; Tustison et al., 2010).

pT2 maps

The T2 maps (Duncan et al., 1996) were calculated from proton-weighted images (first echo) and T2-weighted images (second echo) obtained from the double-echo fast spin echo sequence. Pixel-by-pixel T2 maps were obtained from the two-echo T2 images byWhere S1 and S2 are the first echo and second echo image intensities, with echo times TE1 and TE2, respectively. Rigid transformation implemented in SPM12 (cost function = Normalized Mutual Information; reslicing with trilinear interpolation; no wrapping) was applied to transform the T2-weighted images (S2) to the corresponding T1-weighted images; then, the transformation matrix was applied to co-register the T2 maps to the corresponding T1-weighted images. We used the term “pT2,” since the dual-echo T2 estimation is regarded as semi-quantitative compared to that based on multi-echo sequences with 16 ≤ echos.

DTI

The tensor calculation function implemented in MRICloud (www.MRICloud.org) (Mori et al., 2016) was used to calculate a tensor field from the raw b0 and diffusion weighted images. The tensor-calculation pipeline adopted a state-of-the-art corrections scheme for simultaneous corrections of motion and distortion by eddy current using a cost function based on tensor-fitting residuals and gradient-specific distortion-modeling (Li et al., 2013). The pipeline also applied voxel-by-voxel outlier rejections for occasional intensity loss due to local elastic motions. The MRICloud generated scalar maps, such as the fractional anisotropy (FA) and the trace maps, from the tensor field. An experienced neurologist performed the second quality check, based on visual inspection of the color-coded orientation maps calculated from the tensor field, to ensure the data had no noticeable motion artifacts. The DTIs were linearly co-registered to the corresponding T1-weighted images using rigid transformation, in two steps: 1) co-registration of the b0 image to the corresponding T2-weighted image; and 2) co-registration of the result of 1) to the corresponding T1-weighted image, by applying the transformation matrix obtained from T2-to T1-weighted image co-registration (2.3.2). This concatenated approach was adopted based on the fact that b0 images share very similar contrast with T2-weighted images.

PET

The time frames of each dynamic PiB-PET scan were aligned to the average of the first 2 min to correct for motion (Jenkinson et al., 2002). For registration purposes, we obtained static images by averaging the first 20 min of each dynamic PET scan. The 3D T1-weighted images were rigidly registered onto their corresponding 20-minute average. We generated a study-specific template using T1-weighted images in a diffeomorphic registration framework (Avants et al., 2010b), and then segmented this template using FreeSurfer (version 5.1, http://surfer.nmr.mgh.harvard.edu) (Desikan et al., 2006). The appropriate MRI template was registered onto each subject's T1-weighted image using diffeomorphic registration (Avants et al., 2008). The FreeSurfer segmentation in template space was transformed accordingly onto the PET scan. Consecutive transformations were concatenated so that FreeSurfer segmentations were interpolated only once. Distribution volume ratio (DVR) images were calculated in the native space of each PET image using the simplified reference tissue model with the cerebellar gray matter as reference tissue (Zhou et al., 2007). The brain parcellation map generated by the MRICloud (see 2.5 for more details) on the T1-weighted MRI space was transformed on the DVR image by applying the transformation matrix obtained from the rigid transformation of the T1-weighted image to the 20-minute average PET image. The cDVR was calculated as the average of the DVR values of the cingulate, frontal, parietal (including precuneus), lateral temporal, and the lateral occipital cortical regions, excluding the sensorimotor strip.

Definition of amyloid-positive and amyloid-negative based on DVR

To separate the 87 participants into individuals with minimal amyloid deposition (PiB-) and individuals with elevated levels of amyloid deposition (PiB-), a mean cDVR threshold of 1.02 was used, estimated by the best intersection of the probability density functions of the PiB+ and PiB- groups based on a Gaussian mixture model (Bilgel et al., 2016).

MRI parcellation and matrix generation

The eACA was performed on an atlas-based matrix, which was created by an anatomical parcellation map (a set of labels that follow anatomical ontology) on each T1-weighted image space. The parcellation map was generated by a multi-atlas label-fusion method, in which multiple brain MRIs with the parcellation map are transformed to a target image, and the transformed parcels are fused (label fusion) to create the anatomical parcellation map specific to the target image. The “T1-multi-atlas segmentation” function implemented in the MRICloud was applied to parcellate each of the T1-weighted images into a parcellation map that contained 285 ontological parcels, including both gray and white matter structures (Wu et al., 2016), which was suitable for this study targeting both gray and white matter structures. The function offers a fully automated multi-atlas image parcellation algorithm, which combines the image transformation algorithm, Large Deformation Diffeomorphic Metric Mapping (LDDMM) (Christensen et al., 1997; Granander and Miller, 1996; Miller et al., 1997; Oishi et al., 2009), and the atlas-fusion algorithm based on the joint-label fusion method (Wang et al., 2013), along with a local search algorithm (Coupe et al., 2011). The multi-atlas label-fusion algorithm relies on a collection of pre-parcellated T1-weighted brain MRIs as teaching files (atlas library version 9B, which includes 26 atlases that contain MRIs from cognitively-normal individuals and individuals with cognitive impairments or dementia; 50–90 years old (Wu et al., 2016)) to achieve high parcellation accuracy against a wide range of atrophic changes seen in patients with memory impairments or dementia (Liang et al., 2015; Zu et al., 2017). The parcels are listed in the Appendix Table 1. The parcellation map generated on the T1-weighted image was applied to a co-registered pT2 map and DTI to perform parcel-by-parcel quantification of the volume, pT2, FA, and trace. The atlas-based image parcellation and quantification method was chosen for two reasons: information integration and dimension reduction. Introducing a common ontological framework among the T1-weighted image, the pT2 map, and DTI enables the integration of multimodality information, and can effectively reduce the dimension of the data space while minimizing the loss of meaningful information, which is essential to avoid the curse of dimensionality (Miller et al., 2013; Mori et al., 2013).

The normalized volume of each parcel was calculated by dividing the volume of each parcel by that of the whole brain. Parcels that were located outside the brain parenchyma, such as the cerebrospinal fluid space, skull, bone marrow, scalp, and the choroid plexus, or that were too small for accurate identification, such as the mammillary body, fimbria, optic tracts, substructures in the brainstem and the subcallosal and rostral white matter around the anterior cingulate cortex (ACC), were excluded from the eACA. The parcellation map was further modified by removing non-brain structures, such as the ventricles and sulci, and applied to other image modalities that were co-registered to the T1-weighted images, to quantify the pT2, FA, and trace values of each parcel (see Fig. 1A). To eliminate the potential effects of partial volume (e.g., contamination of signals from non-brain structures) or that of slight misregistration due to geometric distortion, we excluded pixels with a pT2 value of more than 200 ms, a trace value of more than 0.0045 mm2s-1 (to exclude voxels that potentially include signal from the cerebrospinal fluid), and a trace value of less than 0.0015 mm2s-1 (to exclude voxels that potentially include an artifact related to spatial distortion due to adjacent air cavity) from the quantification. The above-mentioned procedures resulted in a volume matrix (87 individuals × 172 parcels), a pT2 matrix (87 individuals × 142 parcels), a FA matrix (87 individuals × 142 parcels), and a trace matrix (87 individuals × 142 parcels). Overall quality control review was performed by one of the authors (K.O.).

Correction for age and sex

Age and sex are known to affect the local volume of the brain, DTI-derived measures, and functional connectivity (Barnes et al., 2010). An adjustment for age and sex is thus necessary for the eACA. We adopted the method proposed by (Dukart et al., 2011) to adjust for the confounding effects. Essentially, confounding variables (age and sex) were modeled as predictors to determine all parcels on each MR modality (normalized volume, pT2, FA, and trace) based on the PiB- group. The following general linear model for a single image modality was solved for the parameters and ,

In the above model, and are age and sex vectors for PiB- subjects. is a raw image feature vector for each parcel, is the image features without effects of age and sex. To obtain image features for the PiB + group, the estimated parameters and were applied to general linear model in terms of .

Given that the PiB- group was significantly younger than the PiB + group (Table 1), we also rebuilt the above linear model by adjusting confounding variables based on both PiB+ and PiB- groups.

Group-wise anatomical correlation matrix

The anatomical correlation was defined as statistical associations in volume or image-intensity between brain structures defined by the parcellation map (see Section 4). Pearson correlation coefficients across subjects within one group (PiB- or PiB+) were calculated to form the anatomical correlation matrix (N×N, where N represents the number of parcels) for each modality, such as the normalized volume matrices (172×172) and other intensity matrices (pT2, FA, trace, 142×142), as shown in Fig. 1B. Since the obtained interregional correlation matrices were symmetric, only the lower triangle part of each matrix, including N×(N-1)/2 connections served our interests.

The agglomerative hierarchical cluster analysis (clustergram)

To explicitly describe the pattern of each interregional correlation matrix, we further applied agglomerative hierarchical cluster analysis for each correlation matrix. Clustering was generated by grouping parcels with similar correlation vectors (1×N) according to minimum linkage based on correlation distance between every pair, as described in (Liu et al., 2012). The hierarchical tree was cut by a prefix threshold (60% of maximum linkage across all parcels) for all modalities. The threshold was defined based on the volume matrix, in an attempt to decrease the threshold of maximum linkage to increase number of clusters, while setting a constraint of a minimum number of parcels (two parcels) included in each cluster.

The eACA and statistical comparisons between PiB+ and PiB-

To investigate the effects of brain amyloid deposition on the anatomical correlation, we compared the anatomical correlation matrices between the PiB- and PiB + groups. Based on previous research that suggested strong correlations among functionally connected anatomical areas, as observed in healthy cohorts (Chen et al., 2008; Lerch et al., 2006), and among areas vulnerable to specific diseases, as observed in cohorts of diseases (Bernhardt et al., 2008; Hafkemeijer et al., 2016; He et al., 2008, 2009; Seeley et al., 2009), we hypothesized that the anatomical correlations in normal subjects would be disrupted by amyloid deposition. To extract this pairwise correlation from all parcel-pairs (N×(N-1)/2), a permutation test was applied to assess the significance of each inter-parcel correlation coefficient between PiB- and PiB+ (Bassett et al., 2008). Permutation was repeated for 5000 times to construct sampling distribution of all correlation coefficients under the null hypothesis that observed correlation differences were not determined by true group factor. A p < 0.05 was considered significant difference after the false discovery rate (FDR) correction (Benjamini and Hochberg, 1995). The results were represented in a Manhattan plot for each image modality, where parcel pairs were depicted along the horizontal axis, and the negative logarithm of the p-value for each parcel pair was depicted on the vertical axis. Under a limited number of permutations (5,000 in our study), the estimated p value could be zero when the effect size is large. On such occasions, we used p = 0.01 instead for the visualization.

Results

The eACA

The correlation matrices derived from multiple image modalities are demonstrated in Fig. 2, in which the correlation coefficients of each combination of parcels were color-coded from blue (-1), to green (0), to red (1). The findings for each modality are summarized in Sections 3.3–3.6, below.

Fig. 2

Correlation matrices for various image modalities. A: PiB + group, B: PiB- group. For each correlation matrix, the structural labels along each axis for volume, pseudo-T2, FA, and trace are presented in the color bar above the matrix.

The hierarchical cluster analysis

The hierarchical clustering results for each correlation matrix were mapped to a brain template, as shown in Fig. 3. Parcels with a similar set of covaried parcels were clustered and represented by colors that were arbitrarily assigned. The details of cluster components across every image modality are listed in Appendix Table 2 – 9. The findings for each modality are summarized in Sections 3.3–3.6, below.

Fig. 3

The agglomerative hierarchical cluster analysis for each correlation matrix. A: volume, B: pseudo-T2, C: FA, D: trace. Clustering was generated by grouping parcels with similar correlation vectors according to the minimum linkage between each pair. The clusters are represented by colors that were arbitrarily assigned.

Normalized volume

The hierarchical clustering identified 23 clusters for the PiB- group, which are visualized in (Fig. 3A – PiB-), and which had a left-right symmetric appearance. The limbic structures (hippocampus, amygdala, cingulum, fornix, stria terminalis, parahippocampal gyrus, and entorhinal cortex) formed a single cluster. In contrast, 18 clusters were identified for the PiB + group, which are visualized in Fig. 3A – PiB+, and which showed a more asymmetric appearance than the PiB- group, especially in the orbitofrontal, parietal, and occipital areas.

No significant difference in the correlation coefficients between the PiB- and the PiB + groups was identified among 14706 (172×171/2) parcel-pairs of normalized volume (Fig. 4A).

Fig. 4

Manhattan plot of anatomical correlation with various image modalities. A: volume, B: trace, C: FA, D: pseudo-T2, E: Scatterplot of correlation between the right entorhinal cortex and the right cerebellum with p < 0.001 (PiB- > PiB+). Each dot in the Manhattan plots represents significance for different correlation coefficients between a pair of anatomical structures. The blue horizontal line represents p = 0.05 (FDR-corrected).

pT2 maps

The hierarchical clustering identified seven clusters for the PiB- group, which are visualized in Fig. 3B – PiB-, and show that all the cortical and deep white matter parcels were clustered to form a single component (Fig. 3B – PiB-). In contrast, five clusters were identified for the PiB + group, which are visualized in Fig. 3B– PiB+, and which indicated two major cortical components, the cerebrum and the cerebellum (Fig. 3B– PiB+). The pT2 values of the parcels of the entire cerebrum were highly correlated.

Among 10011 (142×141/2) parcel-pairs of pT2 values, the correlation coefficients of 45 parcel-pairs were significantly higher in the PiB- group than in the PiB + group (Fig. 4E). Most of the significant parcel-pairs were related to the cerebellum or the hippocampus. Among them, one of the most significant correlations was identified the right entorhinal cortex and the right cerebellum (r for PiB + group = -0.242, r for PiB- group = 0.66, corrected p < 0.001), as shown in Fig. 4F. The complete list of parcel-pairs with significantly different correlations between groups is demonstrated in Table 2.

Table 2

Pairs of structures with a significant difference in the correlation coefficient between the PiB+ and the PiB- groups. NA: not applicable (non of the pairs demonstrated sifnificant difference).

Image modality	Pair of structures		Correlation coefficient for PiB+	Correlation coefficient for PiB-	Corrected p value
	Name of parcel1	Name of parcel1
Volume	NA	NA	NA	NA	NA
pT2	ST_L	PLIC_L	-0.312	0.483	0.044
	Hippo_L	IO_R	-0.043	0.680	0.044
	Hippo_L	EC_L	-0.269	0.651	0.044
	Hippo_L	IFO_L	-0.165	0.671	<0.001
	Hippo_L	MF_L	0.003	0.665	<0.001
	Hippo_L	IF_opercularis_L	-0.003	0.648	0.044
	Hippo_L	IF_triangularis_L	-0.185	0.701	<0.001
	Hippo_L	IF_triangularis_R	-0.047	0.630	<0.001
	Caudate_tail_L	EC_L	-0.581	0.450	<0.001
	MO_L	SCR_L	-0.315	0.478	<0.001
	Cerebellum_L	CGC_L	-0.314	0.614	<0.001
	Cerebellum_L	EC_L	-0.312	0.703	0.044
	Cerebellum_L	ACR_L	-0.430	0.609	<0.001
	Cerebellum_L	SCR_L	-0.404	0.562	<0.001
	Cerebellum_L	SF_L	-0.201	0.707	0.044
	Cerebellum_L	SF_PFC_L	-0.282	0.693	<0.001
	Cerebellum_L	CGC_R	-0.341	0.605	0.044
	Cerebellum_L	EC_R	-0.383	0.694	<0.001
	Cerebellum_L	IFO_R	-0.324	0.688	<0.001
	Cerebellum_L	ACR_R	-0.447	0.596	0.044
	Cerebellum_L	PoC_R	-0.024	0.755	<0.001
	Cerebellum_L	SF_PFC_R	-0.286	0.670	<0.001
	PVWa_R	A_L	-0.353	0.455	<0.001
	PVWa_R	MO_L	-0.282	0.457	<0.001
	PVWa_R	A_R	-0.211	0.507	0.044
	PVWa_R	IO_R	-0.356	0.430	0.044
	SCR_R	CerebellumWM_L	-0.127	0.703	0.044
	Caudate_tail_R	IO_R	-0.143	0.630	0.044
	Caudate_tail_R	EC_R	-0.340	0.465	0.044
	Caudate_tail_R	CGC _L	-0.344	0.421	0.044
	Caudate_tail_R	EC_L	-0.512	0.462	<0.001
	Fu_R	PLIC_L	-0.427	0.449	<0.001
	PCC_R	PLIC_L	-0.358	0.432	<0.001
	Cerebellum_R	EC_L	-0.301	0.678	0.044
	Cerebellum_R	ACR_L	-0.342	0.599	<0.001
	Cerebellum_R	SCR_L	-0.437	0.552	<0.001
	Cerebellum_R	CGC_L	-0.179	0.584	0.044
	Cerebellum_R	IF_triangularis_R	-0.038	0.631	0.044
	Cerebellum_R	SF_L	-0.110	0.687	<0.001
	Cerebellum_R	SF_pole_R	-0.052	0.609	0.044
	Cerebellum_R	SF_PFC_L	-0.109	0.682	<0.001
	Cerebellum_R	EC_R	-0.242	0.661	<0.001
	Cerebellum_R	ACR_R	-0.305	0.581	<0.001
	Cerebellum_R	PrC_R	-0.008	0.737	0.044
	Cerebellum_R	SF_PFC_R	-0.155	0.648	<0.001
FA	MO_L	RLIC_L	-0.455	0.407	<0.001
	CGC_R	IT_R	-0.584	0.531	<0.001
	MFO_L	Subgenual_ACC_L	-0.391	0.556	<0.001
	SCC_R	MF_R	-0.079	0.576	<0.001
Trace	NA	NA	NA	NA	NA

Abbreviation of brain structures is explained in Appendix Table 1.

DTI

FA

The hierarchical clustering identified seven clusters for the PiB- group, which are visualized in Fig. 3C – PiB-. The FA values of the parcels within the neocortex (a cluster marked by yellow), within the paleocortex (piriform cortex, periamygdalar area, and olfactory areas; a cluster marked by red), and within the diencephalon-rhombencephalon complex (cerebellum-thalamus; a cluster marked by green), were all highly correlated. In contrast, eight clusters were identified for the PiB + group. The paleocortex cluster seen in the PiB- was fragmented and a new cluster was seen in the deep white matter area. Moreover, most of the cerebrum and cerebellar cortices were highly correlated.

Among 10011 (142×141/2) parcel-pairs of FA values, four parcel-pairs of FA between the right cingulum and the right interior temporal gyrus (r for PiB + group = -0.584, r for PiB- group = 0.531, corrected p < 0.001), the left retrolenticular part of internal capsule and the left middle occipital gyrus (r for PiB + group = -0.455, r for PiB- group = 0.407, corrected p < 0.001), the right splenium of corpus callosum and the right middle frontal gyrus (r for PiB + group = -0.079, r for PiB- group = 0.576, corrected p < 0.001), and the left subgenual anterior cingulate cortex and the left middle fronto-orbital gyrus (r for PiB + group = -0.391, r for PiB- group = 0.556, corrected p < 0.001) were identified as correlation coefficients that were significantly lower in the PiB + group than in the PiB- group, as shown in Fig. 4C and Table 2.

Trace

The hierarchical clustering identified nine clusters for the PiB- group, as visualized in Fig. 3D – PiB-. Among them, two major clusters were observed in the cortex; one related to the frontopolar-orbitofrontal and the right frontotemporal areas (a cluster marked by purple in Fig. 3D), and another related to the parieto-occipital and left frontotemporal areas (a cluster marked by green in Fig. 3D). In contrast, eight clusters were identified for the PiB + group, as visualized in Fig. 3D – PiB+. An additional cluster in the cortical area, the dorsal frontal cluster, which was not seen in the PiB- group, was seen in the PiB + group (marked by brown in Fig. 3D – PiB+).

No significant difference in the correlation coefficients between the PiB- and the PiB + groups was identified among 10011 (142×141/2) parcel-pairs of trace values (Fig. 4B).

Confounding adjustment based on the entire cohort

Even after replacing the covariates from PiB- with those from the entire cohort in the eACA confounding adjustment, the significant pairs of structures remained the same for the volume and DTI (shown in Appendix Table 10) analysis. The only exception was seen in the pT2. Thirty of 46 significant pairs of structures for pT2 were unaffected by adjustment changes. As described in Section 3.4, most of the significant parcel-pairs were still related to the cerebellum or the hippocampus.

Discussion

The eACA method described here provides the ability to visualize a “normal correlation network” and the disruption related to amyloid deposition in a comprehensive way, via the correlation matrix and the clustergram. Anatomical correlation is assumed to be modulated by developmental coordination, synaptic maturation, and other neuropathologies (Alexander-Bloch et al., 2013; Hof and Morrison, 2004; Montembeault et al., 2016). Previous studies about anatomical covariance related to neurological diseases have primarily targeted on cortical thickness (Avants et al., 2010a), surface area (Palaniyappan et al., 2011), and gray (Seeley et al., 2009) and white matter volumes (Hulshoff Pol et al., 2004). Our study is the first, to our knowledge, to systematically explore the whole-brain anatomical covariance associated with amyloid pathology in a cognitively normal population. The primary finding was evidence for a disruption of the physiological intrinsic brain network apparent on pT2 values and FA in the PiB + group.

Cross-modality comparisons in the eACA

Overall, the eACA results demonstrated that each modality had a unique pattern of physiological correlations. The number of clusters and the patterns observed in this study were based on 60% of the maximum linkage in the hierarchical tree (Liu et al., 2012), although the method to define an optimal modality-specific threshold has not been fully investigated. Further studies are needed to investigate the optimal threshold to capture the landscape of anatomical correlations, the molecular and cellular underpinnings of the macroscopic pattern observed in our study, and the mechanisms of covariate patterns specific to each modality.

The left-right symmetric appearance of the physiological correlation would suggest strong connections between bilateral corresponding anatomical areas via commissural fibers (Mechelli et al., 2005). The result of hierarchical clustering indicated that the clusters obtained from the PiB- group represented functional units formed by regionally connected areas, such as a cluster that involves the Papez circuit (the entorhinal cortex, parahippocampal gyrus, hippocampus, fornix, cingulum, thalamus, hypothalamus, and the stria terminalis) and the amygdala, which are crucial in memory modulation (Hamann et al., 1999; Markowitsch, 1997). Although such symmetric and functionally meaningful clusters seen in the PiB- were disrupted in the clustergram obtained from the PiB+, the pairwise correlation analysis demonstrated that the physiological correlation seen in the PiB- group was preserved in the PiB + group.

T2 or pT2 values in the brain are known to be affected by various micro-environmental factors of the brain tissue, especially the existence of myelin and macromolecules (Du et al., 2012; Miot-Noirault et al., 1997) (Mosher and Dardzinski, 2004). In this study, pT2 values of entire brain parcels were well correlated in both the PiB- and the PiB + groups, except for a dissociation between the cerebrum and the cerebellum observed in the PiB + group.

The pairwise correlation analysis demonstrated that the physiological correlation seen in the PiB- group was widely altered in the PiB + group. Among the 71 parcel-pairs with alteration, 47 pairs were related to the cerebellum, indicating stable pT2 values for the cerebellum independent of the pT2 values of the cerebral structures. This finding was congruent with the results of clustering analysis, suggesting that, similar to the evaluation of amyloid-PET (Catafau et al., 2016), the pT2 value of the cerebellum could be a reference for the evaluation of amyloid deposition in the cerebrum. Namely, for the PiB- group, pT2 value among the regions connected to the cerebellum was homogeneous (little intra-individual variation) but the inter-individual variation of the pT2 value was substantial. For the PiB + group, this cerebellum-related correlation network was disrupted. Seven pairs were related to the hippocampus, indicating consistently low pT2 values of the hippocampus, independent of the pT2 values of other structures. The scatterplot (Fig. 4F) clearly demonstrates the strength of the eACA approach, in which the two groups cannot be differentiated based on the pT2 of each structure, but a difference in the regression lines can distinguish the groups. Interestingly, for the cerebellum- and hippocampus-related correlation networks, the PiB + group demonstrated the negative correlation between pair of regions, while PiB- group always demonstrated strong positive correlation. The implication is not straightforward, but this tendency seen in the PiB + group: pT2 value of the cerebellum or hippocampus tends to be lower for the individuals with higher pT2 value in the connected structures, might indicate physiological or histopathological alterations that result in lower pT2 value in these regions. Indeed, both cerebellum and hippocampus are known as structures with altered BOLD response related to amyloid pathology (Kennedy et al., 2018; Quiroz et al., 2010).

The left-right symmetric appearance of the physiological correlation was remarkable in FA, which suggested strong connections between bilateral corresponding anatomical areas via commissural fibers.

The pairwise correlation analysis demonstrated that the physiological correlation between the cingulum bundle and the inferior temporal lobe, and between the superior longitudinal fasciculus and the superior parietal lobule, were altered in the PiB + group. A relationship between amyloid burden and altered DTI measures in the selective white matter bundles, particularly the cingulum has been repeatedly reported in amnestic MCI (Douaud et al., 2013; Mito et al., 2018; Svard et al., 2017) and even in cognitively-normal elderly (Gold et al., 2014; Molinuevo et al., 2014; Racine et al., 2014; Rieckmann et al., 2016). The dissociation between the altered white matter bundles and preserved cortical structure might cause the altered physiological correlation.

There was a striking physiological asymmetry seen in the trace clusters. The left language areas were clustered with the bilateral prefrontal area and the bilateral parieto-occipital areas, while the right side of the corresponding areas were clustered with basal and polar aspects of the bilateral frontal lobe, the pattern of which follows the macroscopic anatomical asymmetries seen in the normal brain (Lancaster et al., 2003; Toga and Thompson, 2003). Although the right dominant cluster was separated into two sub-clusters in the PiB + group, the pairwise correlation analysis demonstrated no differences between the PiB- and PiB + groups, indicating the robustness of the trace value to amyloid deposition.

Limitations

Several limitations should be noted. First, the ground-truth of the physiological intrinsic brain network seen in this study is unknown. This lack of ground-truth is a common problem among connectome- or connectivity-related studies. Validations based on histology-MRI comparative studies are highly anticipated. Second, differences in the eACA between the PiB- and the PiB + groups should be interpreted as features that are related to amyloid deposition in the brain, but which may not be a definitive indication of future Alzheimer's disease, since it is not yet known which of these individuals might progress to mild cognitive impairment or dementia. There is a potential that the results of eACA might vary, depending on other physiopathological factors, such as genotype (Chen et al., 2011; Panizzon et al., 2012) and Tau deposition (see 4.5). Further classification of the study population might be needed in a future study. Third, interpretation of the hierarchical clustering followed by visualization of difference between groups was descriptive; while the permutation test applied to the correlation matrices was intended to identify group differences, the clustering approach was intended to capture spatial distribution of the correlation network. Forth, the generalizability of the proposed eACA approach to images obtained on different scanners or scan protocols needs to be evaluated. In particular, it is known that scalar values derived from DTI are affected by scan parameters. External validation of the results should be performed in the future. Fifth, DTI-to-T1-weighted image co-registration was affected by geometric distortion. Although we applied thresholds (see Section 2.5) to ameliorate the effect of cerebrospinal fluid contamination or signal dropout due to susceptibility artifacts, the effects of distortion and misregistration might have remained in the results. Sixth, each modality has its own noise patterns and variations in measurements related to both biological and methodological effects, as seen in Fig. 4. However, there is no general consensus about the appropriate ways to set modality-specific thresholds to detect pathological changes. As an initial attempt, we applied the same statistical threshold for all modalities, but the biologically optimal method for the modality-specific thresholding has yet to be investigated. Seventh, the small sample size of our study might have resulted in low sensitivity to detect mild anatomical features related to amyloid deposition. Ultimately, to investigate the long-term dynamics of the altered physiological correlation, longitudinal studies are necessary.

Conclusion

In conclusion, we adopted an eACA approach to describe the physiological correlation identified in a comparison of cognitively normal individuals with and without evidence of significant amyloid deposition. The pattern of physiological correlation was modality-specific, and congruent with previous research based on conventional ACA. The eACA approach could detect altered physiological correlation in T2 and FA in cognitively normal elderly with amyloid deposition. The eACA has the potential to evaluate alterations of anatomical status related to pathological processes.

Declarations

Author contribution statement

Chenfei Ye, Kenichi Oishi: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Wrote the paper.

Timothy Brown: Analyzed and interpreted the data.

Murat Bilgel, Johnny Hsu: Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data.

Marilyn Albert, Ting Ma, Michael I. Miller, Susumu Mori: Contributed reagents, materials, analysis tools or data.

Brian Caffo: Conceived and designed the experiments; Analyzed and interpreted the data.

Funding statement

This work was made possible by the Dana Foundation Clinical Neuroscience Research Program, grant P41EB015909 from the National Institutes of Health, the Pilot Project Discovery Program from the Johns Hopkins Individualized Health Initiative (inHealth), the Fakhri Rad BriteStar award from the Department of Radiology Johns Hopkins University School of Medicine, a grant from the Basic Research Foundation of Shenzhen Science and Technology Program (JCYJ20150403161923510), the Basic Research Foundation Key Project Track of Shenzhen Science and Technology Program (JCYJ20170413110656460), and was supported, in part, by the Intramural Research Program of the National Institute on Aging.

Competing interest statement

The authors declare the following conflict of interests: Dr Miller and Dr Mori are the co-founders, and Dr Oishi is a consultant, for Anatomy Works. The terms of this arrangement are being managed by the Johns Hopkins University in accordance with its conflict of interest policies.

Additional information

No additional information is available for this paper.