Gaussian Graphical Models Reveal Inter-Modal and Inter-Regional Conditional Dependencies of Brain Alterations in Alzheimer's Disease

Alzheimer's disease (AD) is characterized by a sequence of pathological changes, which are commonly assessed in vivo using various brain imaging modalities such as magnetic resonance imaging (MRI) and positron emission tomography (PET). Currently, the most approaches to analyze statistical associations between regions and imaging modalities rely on Pearson correlation or linear regression models. However, these models are prone to spurious correlations arising from uninformative shared variance and multicollinearity. Notably, there are no appropriate multivariate statistical models available that can easily integrate dozens of multicollinear variables derived from such data, being able to utilize the additional information provided from the combination of data sources. Gaussian graphical models (GGMs) can estimate the conditional dependency from given data, which is conceptually expected to closely reflect the underlying causal relationships between various variables. Hence, we applied GGMs to assess multimodal regional brain alterations in AD. We obtained data from N = 972 subjects from the Alzheimer's Disease Neuroimaging Initiative. The mean amyloid load (AV45-PET), glucose metabolism (FDG-PET), and gray matter volume (MRI) were calculated for each of the 108 cortical and subcortical brain regions. GGMs were estimated using a Bayesian framework for the combined multimodal data and the resulted conditional dependency networks were compared to classical covariance networks based on Pearson correlation. Additionally, graph-theoretical network statistics were calculated to determine network alterations associated with disease status. The resulting conditional dependency matrices were much sparser (≈10% density) than Pearson correlation matrices (≈50% density). Within imaging modalities, conditional dependency networks yielded clusters connecting anatomically adjacent regions. For the associations between different modalities, only few region-specific connections were detected. Network measures such as small-world coefficient were significantly altered across diagnostic groups, with a biphasic u-shape trajectory, i.e., increased small-world coefficient in early mild cognitive impairment (MCI), similar values in late MCI, and decreased values in AD dementia patients compared to cognitively normal controls. In conclusion, GGMs removed commonly shared variance among multimodal measures of regional brain alterations in MCI and AD, and yielded sparser matrices compared to correlation networks based on the Pearson coefficient. Therefore, GGMs may be used as alternative to thresholding-approaches typically applied to correlation networks to obtain the most informative relations between variables.


Introduction
Alzheimer's disease (AD) is characterized by a range of pathological brain alterations that can be assessed in vivo using various neuroimaging methods, including MRI and PET. Several studies suggest that information obtained from combining different imaging modalities could provide reliable markers of cerebral reserve capacity and might be used to predict and monitor the evolution of AD and its relative impact on cognitive domains in pre-clinical, prodromal, and dementia stages of AD (see e.g. reviews [Teipel et al., 2015a, Teipel et al., 2016). However, there is still an unmet need for appropriate analysis methods for assessing statistical associations between individual brain regions and between different pathology markers derived from multiple neuroimaging modalities. Up to date, multimodal studies employ one of the following approaches: (i) Correlation of pathology maps on a voxel level [Altmann et al., 2015, Grothe and Teipel, 2016, La Joie et al., 2012; (ii) linear regression analysis with a-priori specified regions-of-interest [Buckner et al., 2005, Seeley et al., 2009, Villain et al., 2010, Kljajevic et al., 2014, Chang et al., 2015, Grothe et al., 2016, Teipel and Grothe, 2016; (iii) stratification of subjects into distinct groups (e.g. amyloid positive/negative) to compare differences in other imaging modalities [Buckner et al., 2005, Kljajevic et al., 2014, Grothe et al., 2016; (iv) comparison of graph-theoretical measures and statistics between modalities [Stam et al., 2006, Buckner et al., 2009, Zhou et al., 2012, Sepulcre et al., 2013, Sepulcre et al., 2017; and (v) estimation of generative models for comparing spreading mechanisms of amyloid-β deposition and its contribution to neurodegeneration [Iturria-Medina et al., 2017, Torok et al., 2018. Commonly employed statistical models, such as linear regression analysis, provide limited ability to assess the interactions between dozens of variables in the same model, as they cannot derive reliable estimates regarding the individual contribution of highly collinear predictors and suffer from variance inflation [Dormann et al., 2013]. Calculation of covariance/connectivity matrices based on the Pearson correlation between each pair of variables has led to practical problems in deriving meaningful results, i.e. these matrices are commonly thresholded to an a-priori defined density and binarized [Buckner et al., 2009, Zhou et al., 2012, Sepulcre et al., 2013. More recently, summary statistics based on graph-theory have been proposed [Watts andStrogatz, 1998, Stam et al., 2006] and are currently widely applied [Buckner et al., 2009, Zhou et al., 2012, Sepulcre et al., 2013, Sepulcre et al., 2017. However, this approach has been criticized, as for instance, group differences in small-worldness of the brain network might be sensitive to the specific density threshold [Hlinka et al., 2017, Mårtensson et al., 2018.
We suggest the application of Gaussian graphical models (GGMs), which are able to estimate the partial correlation between various multicollinear predictors [Hastie et al., 2013, ch. 7.3]. GGMs yield sparse conditional dependency matrices, that are conceptually expected to closer reflect the underlying causal relationships [Bontempi andFlauder, 2015, Koller andFriedman, 2009, chapter 21.7]. This makes GGMs an interesting candidate for studying properties of the brain network; an example is illustrated in Figure 1. The partial correlation derived from GGMs is conceptually similar to the partial correlation obtained from a series of linear regression models, which estimate the statistical association of the dependent and independent variables while controlling for the confounding variables. Additionally, GGMs extend this concept by estimating the partial correlation matrix as a set of coupled regression problems, in contrast to separate regression problems modeled by traditional linear regression [Meinshausen andBühlmann, 2006, Hastie et al., 2013, ch. 7.3]. Technically, GGMs are naively realized by matrix inversion of the covariance matrix. In more robust and efficient approaches, regularization techniques [Meinshausen and Bühlmann, 2006, Ravikumar et al., 2011, Ryali et al., 2012, Cai et al., 2013, Wang et al., 2016 or efficient sampling schemes [Mohammadi andWit, 2015, Mohammadi andWit, 2019] are applied.
In this paper, we tested the applicability and clinical utility of GGMs to reveal the conditional dependency structure between regional pathology measures. For this purpose, we assessed inter-regional statistical associations within and between three main imaging markers of Alzheimer's disease using GGMs based on a whole-cortex parcellation of the brain. The assessed imaging markers included amyloid-β deposition (florbetapir/AV45-PET), glucose metabolism (FDG-PET), and gray matter volume (T 1 -weighted MRI). Based on our previous results with only six representative brain regions [Dyrba et al., 2017], we hypothesized that regional amyloid deposition has low contribution to gray matter atrophy, whereas hypometabolism was expected to be stronger related to atrophy. Further, we expected a few hub-nodes influencing pathology in other regions. For graph-theoretical measures, we expected a linear trajectory of decreasing clustering coefficient and increasing path length with stronger disease severity, as previously reported in the literature for connectivity analyses based on Pearson correlation [He et al., 2008, Yao et al., 2010, Li et al., 2012, Morbelli et al., 2012, Tijms et al., 2013, Pereira et al., 2016, John et al., 2017, Titov et al., 2017.

Study participants
Data were obtained from the Alzheimer's Disease Neuroimaging Initiative (ADNI) database (http://adni.loni.usc.edu). The ADNI was launched in 2003 by the National Institute on Aging, the National Institute of Biomedical Imaging and Bioengineering, the Food and Drug Administration, private pharmaceutical companies, and non-profit organizations, with the primary goal of testing whether neuroimaging, neuropsychological, and other biological measurements can be used as reliable in vivo markers of Alzheimer's disease pathogenesis. A complete description of ADNI and up-to-date information is available at http://www.adni-info.org. For this study, 529 subjects with amnestic mild cognitive impairment (MCI), 189 patients with Alzheimer's dementia (AD), and 254 cognitively healthy control subjects (CN) were selected from the ADNI-GO, ADNI-2 and ADNI-3 extensions of the ADNI project, based on the availability of concurrent structural MRI, FDG-PET, amyloid-sensitive AV45-PET, and neuropsychological assessments. In ADNI, two MCI subgroups exist, which only differ by the less severe impairment of memory function for early MCI (EMCI) compared to late MCI (LMCI) subjects. Detailed inclusion criteria for the diagnostic categories can be found at the ADNI website (http://adni.loni.usc.edu/methods, ADNI2 manual page 27). Demographics and neuropsychological profiles of the different diagnostic groups are summarized in Table 1. 2.2 Imaging data and feature extraction ADNI-GO/2 MRI, FDG-and AV45-PET data were downloaded from the ADNI image archive. ADNI-GO/2 MRI data were acquired on multiple 3T MRI scanners using scanner-specific T1-weighted sagittal 3D MP-RAGE/IR-SPGR sequences. To increase signal uniformity across the multicenter scanner platforms, original T1 acquisitions underwent standardized image preprocessing correction steps (http://adni.loni.usc.edu/methods/mri-tool/mri-pre-processing/). FDG-and AV45-PET data were acquired on multiple instruments of varying resolution and following different platform-specific acquisition protocols. Similar to the MRI data, PET data in ADNI were also subject to standardized image preprocessing correction steps, with the aim of increasing data uniformity across the multicenter acquisitions (http://adni.loni.usc.edu/methods/pet-analysis-method/pet-analysis/). Imaging data were processed by using statistical parametric mapping (SPM8, Wellcome Centre for Human Neuroimaging, University College London) and the VBM8 toolbox (Structural Brain Mapping Group, University of Jena) implemented in MATLAB R2013b (Math-Works, Natick, MA) as previously described in [Grothe et al., 2016, Grothe andTeipel, 2016]. First, MRI T1 scans were segmented into gray matter, white matter and cerebrospinal fluid partitions using the segmentation routine of the VBM8 toolbox. Then, the resulting gray matter and white matter segments were spatially normalized to an aging/AD-specific reference template [Grothe et al., 2013] using the DARTEL algorithm. Additionally, voxel values of the normalized gray matter segments were modulated for volumetric changes introduced by the high-dimensional normalization, such that the total amount of gray matter volume present before warping was preserved. Each subject's FDG-and AV45-PET scans were rigidly coregistered to the corresponding skull-stripped T1 scan. Then, the PET scans were corrected for partial volume effects using a three-compartment model and the MRI-derived tissue segments [Müller-Gärtner et al., 1992, Gonzalez-Escamilla et al., 2017. Corrected PET scans were spatially normalized (without modulation) by applying the deformation fields of the T1-weighted scans. All original data and normalized scans were visually inspected to ensure a high quality of the data. Subsequently, mean gray matter volumes and mean FDG-/AV45-PET uptake values were calculated for 108 cortical and subcortical regions defined by the Harvard-Oxford atlas [Desikan et al., 2006] after projecting the atlas to the aging/AD-specific reference space and removing voxels with a gray matter probability of less than 50% in the aging/AD template. Finally, regional gray matter volumes were proportionally scaled by total intracranial volume (TIV), regional FDG-PET values were proportionally scaled to pons uptake, and regional AV45-PET values were proportionally scaled to whole-cerebellum uptake. To be able to directly compare the different modalities, all regional values were normalized using the congitively normal subjects as reference group [La Joie et al., 2012]. As described previously [Dyrba et al., 2017], we used the so-called W -scores, which are analogous to Z-scores but are adjusted for specific covariates; age, gender, and education in the present case. Like Z-scores, W -scores have a mean value of 0 and a standard deviation of 1 in the control group, and values of +1.65 and −1.65 correspond to the 95 th and 5 th percentiles, respectively. To calculate the W -scores, regression models were estimated for the control group using age, gender, and education as independent variables and the mean value of each region as dependent variable. Then, W -scores were computed using W = (x ij − e ij )/s res,j ; with x ij being the i th subject's raw value for region j; e ij being the value expected for region j in the control group for the i th subject's age, gender, and education; and s res,j being the standard deviation of the residuals for region j in controls.

Statistical modeling
Graphical models provide an effective way for describing statistical patterns in multivariate data and for estimating the conditional dependency between the various brain regions and imaging modalities based on GGMs [Lauritzen, 1996, Mohammadi andWit, 2015]. For data following a multivariate normal distribution, undirected GGMs are commonly used. In these graphical models, the graph structure is directly characterized by the precision matrix, i.e. the inverse of the covariance matrix: non-zero entries in the precision matrix show the edges in the conditional dependency graph. Notably, simple inversion of the covariance matrix usually does not work in real world data sets, as already slight noise in the empirical data causes the precision matrix to contain almost no zero entries. To overcome this problem, regularization techniques or efficient sampling algorithms have been proposed that reduce the effect of noise by additionally employing a sparsity assumption and, thus, only detect the most probable conditional dependencies. For our analyses, we employed a computationally efficient Bayesian framework implemented in the R package BDgraph. More specifically, this framework implements a continuous-time birth-death Markov process for estimating the most probable graph structure and edge weights that correspond to the observed partial correlations [Mohammadi andWit, 2015, Mohammadi andWit, 2019]. For this study, BDgraph was substantially extended by multi-threaded parallel processing and marginal pseudo-likelihood approximation to speed up computations.

Experimental setup
First, we estimated GGMs based on the combined data of EMCI, LMCI and AD patients to study the conditional dependency between brain regions and modalities. Second, we estimated GGMs for each diagnostic group separately to assess alterations of the graph structures. For the combined model, regional W -scores of all MCI and AD patients (N = 718) and all three imaging modalities were entered. Initially, we took all 108 cortical and subcortical regions included in the Harvard-Oxford atlas [Desikan et al., 2006] into consideration, corresponding to P = 3 * 108 = 324 variables. The sampling process included 1,000,000 burn-in iterations 1 , starting from a random estimate for the inverse covariance matrix and converging to estimates with higher posterior probability giving the training data. The burn-in iterations were then discarded, and subsequently 150,000 sampling iterations followed to obtain the estimates for the inverse covariance matrix. Because results were showing a strong left-right hemisphere symmetry, we repeated model estimation including only the 54 regions in the left hemisphere (P = 3 * 54 = 162 variables) to increase model stability. From the final model we set a probability threshold of P avg > 0.5 for selecting the edges, with the notion that a specific edge was considered to be present if it existed in at least half of all model iterations [Madigan et al., 1996]. For the second analysis of group differences, we estimated individual GGMs for each group based on the multimodal data of the left hemisphere. Sampling was again performed with 1,000,000 burn-in iterations followed by 150,000 sampling iterations.
For comparison, these analyses were also repeated (i) using data of the right hemisphere to validate the results and (ii) using the traditional approach of constructing correlation networks based on the Pearson correlation coefficient.

Graph-theoretical analyses
To assess group differences of the estimated graph structure we calculated the three graph-theoretical measures that are most commonly reported in the literature; clustering coefficient, characteristic path length, and their ratio, the small-world coefficient. The path length quantifies the distance of connections between two nodes along the shortest path. The weighted characteristic path length is the average minimum distance between a node i ∈ N and all other nodes, L i = j∈N,j =i d ij /(n − 1), where d ij = auv∈gi↔j ω uv is the shortest weighted path length between i and j, g i↔j defines the shortest path, and ω uv defines the distance between two nodes. Here, the distance matrix was defined as Ω = 1 − abs(Θ), that is one minus the absolute pair-wise partial correlation as derived from the GGMs or the absolute Pearson coefficient, respectively [Rubinov and Sporns, 2010]. The weighted clustering coefficient indicates the inter-connectedness of neighboring nodes C i = 2t i /(k i (k i − 1)), where t i = 0.5 j,h∈N (ω ij ω ih ω jh ) 1/3 is the geometric mean of triangles around node i, and where k i = j∈N a ij is the number of nodes connected to node i [Onnela et al., 2005, Rubinov andSporns, 2010]. k i is often referred to as the degree of the node i, and the link status a ij = 1 if node i is connected to another node j, or a ij = 0 otherwise. The small-world coefficient is defined as the ratio of the clustering coefficient C and characteristic path length L in comparison to a random network, S = (C/C rand )/(L/L rand ), with S 1 in small-world networks [Rubinov and Sporns, 2010]. To simplify calculations, we omitted defining a random network to estimate C rand and L rand , and directly took the ratio S i = C i /L i for group comparisons. Notably, we later report the distribution of graph measures for single regions, as the dependency measures were derive from the whole group of subjects. Graph metrics were compared between diagnostic groups using analysis of variance (ANOVA) and Tukey's honest significant difference tests.

Conditional dependency of Alzheimer's pathology
The conditional dependency matrix obtained using the GGM approach for all region of the left hemisphere is given in Figure 2 (right). For the partial correlation between all pairs of brain regions, we obtained 960 significant associations (7% network density) surviving the posterior probability threshold of P > 0.5 (see Supplementary Figure S1 showing the probability of links). For comparison, the Pearson correlation matrix is given in Figure 2 left). We obtained approximately 6,000 significant Pearson correlations (P < 0.05, Bonferroni corrected), corresponding to a network density of 46% of the total number of possible edges.
For intra-modal associations, i.e. within the same imaging modality, brain regions directly adjacent to each other formed smaller clusters of high partial correlation around the main diagonal (Figures 3,4,and 5). When considering inter-modal associations, i.e. between different imaging modalities, we obtained a consistent pattern of significant positive intra-regional conditional dependency for the pairs amyloid-β deposition and metabolism with a mean partial correlation of ρ = 0.21 for 43 significant associations. These are visible as the higher intensities in the diagonal of Supplementary Figure S2. Between amyloid-β and gray matter volume as well as between metabolism and gray matter volume, only few significant intra-regional associations were found (Supplementary Figures S3 and S4).

Group comparison of the graph structures
When estimating separate models for each diagnostic group based on the multimodal data, graph structures derived from Pearson and partial correlation matrices (Figures 6,7,and 8) both differed in their density, leading to significant alterations of the clustering coefficient, characteristic path length, and small-world coefficient ( Figure 9 and Supplementary Figure S9). We observed a biphasic trajectory of the graph measures. This means that the clustering coefficient and small world coefficient initially increases when comparing early MCI and CN participants ( Figure 9). When Alzheimer's disease progresses, i.e. in the late MCI and dementia groups, both measures decrease again, with late MCI being approximately on the same level as CN participants ( Figure 9). The characteristic path length showed a similar pattern across groups, but with inverted directionality. All blocks showed significant differences in mean between groups, one-way analysis of variance (ANOVA), df = 215, F ≥ 4, p < 0.01, η 2 ≥ 0.055. Detailed results are provided in Supplementary Table S2. P-values for Tukey's honest significant difference tests are provided in Table 2 and Supplementary Table S1. Graph statistics obtained from the right hemisphere data (Supplementary Figure S8) were largely consistent with strongest agreement for the characteristic path length metric.

Conditional dependency between brain regions.
The GGMs estimated the strongest conditional dependencies mainly within imaging modalities. We expected adjacent brain regions to form clusters with high inter-cluster similarity for amyloid-β deposition (Figure 3), as it is known to have low variability in spatial distribution and, therefore, is often used as a dichotomic variable after applying a certain threshold to the global amyloid tracer uptake [Chételat et al., 2013, Landau et al., 2013 or as four-stage variable derived from a linear spreading pattern [Grothe et al., 2017, Sakr et al., 2019. We also found such clustering patterns for metabolism ( Figure 4) and gray matter volume ( Figure 5), matching previous studies on metabolism and gray matter covariance networks based on Pearson correlation [Yao et al., 2010, Carbonell et al., 2016, Pereira et al., 2016 or principal component analysis [Di and Biswal, 2012, Spetsieris et al., 2015, Savio et al., 2017. Clusters of high covariance have been found in the lateral and medial parietal lobe, lateral frontal lobe, and lateral and medial temporal lobe, and had been associated with simultaneous growth during brain development, functional co-activation, and axonal connectivity in the literature [Gong et al., 2012, Alexander-Bloch et al., 2013.
Our analyses yielded only few and relatively weak associations between different modalities ( Figure S2-S4), except for the direct intra-regional dependency between amyloid-β and metabolism as well as between amyloid and gray matter volume (diagonal of Figure S2 and Figure S4), which matched our previous analysis with six selected regions of interest [Dyrba et al., 2017]. The positive dependency between amyloid-β and metabolism was strongest in the early MCI group and matches previous results for partial correlation obtained from linear regression models [Altmann et al., 2015]. This previous study reported a markedly reduced number and strength of negative associations between regional amyloid-β and metabolism when correcting for global amyloid load. They concluded that the negative association between amyloid deposition and metabolism is more related to the global amyloid level than to the distinct regional level. The pattern of intra-regional dependency between amyloid-β and metabolism as well as between amyloid-β and gray matter volume was strongest in the early MCI group, which could refer to the early phase of the disease and, therefore, a high variation in regional amyloid-β deposition and a strong contribution of the amyloid level on both metabolism and volume [Drzezga et al., 2011, Carbonell et al., 2016. Notably, conditional dependencies between metabolism and volume were obtained only for few regions including hippocampus and putamen, but not for other expected regions such as posterior cingulate cortex [Teipel and Grothe, 2016] (Supplementary Figure S3).

Alterations of graph measures
Various studies reported a network disruption of AD in comparison to cognitively healthy controls for gray matter volume [He et al., 2008, Yao et al., 2010, Li et al., 2012, Tijms et al., 2013, John et al., 2017 and glucose metabolism [Morbelli et al., 2012, Titov et al., 2017, and intermediate levels for volume in MCI [Yao et al., 2010, Pereira et al., 2016; which we could replicate in our sample (Supplementary Figure  S9). However, it has to be noted that for Pearson correlation matrices usually high thresholds are applied to obtain sparser graphs. Chung et al. [Chung et al., 2016] and Voevodskaya et al. [Voevodskaya et al., 2017] reported a high influence of the selected graph density threshold on the graph measures, leading to divergent increases and decreases of the global clustering coefficient metric. To circumvent such problems, we used weighted versions of the graph measures [Rubinov and Sporns, 2010] and proposed GGMs to obtain sparse conditional dependency matrices. Our results suggest that graph statistics for regional dependency networks follow a biphasic trajectory in the course of AD, a pattern that was recently also reported for cortical thinning and mean diffusivity [Montal et al., 2017] and resting-state fMRI connectivity [Schultz et al., 2017].
In the current study, the EMCI group displayed the strongest alterations of network structure with an increase of the clustering coefficient, which may relate to the process of amyloid accumulation taking place in several regions simultaneously in this group increasing the intra-cluster correlation. For amyloid-β and volume, LMCI subjects showed a clustering coefficient and small-world coefficient comparable to controls, in contrast to metabolism, where this group showed strongest deviation from the other groups (Table 2). The lowest alterations of graph measures were obtained for the gray matter network.
GGMs were recently applied as clustering algorithm for brain networks in a few other single-modality applications. De Vos et al. [de Vos et al., 2017] found them useful for increasing group separation between AD and controls compared to classical Pearson correlation networks in resting-state functional connectivity. Titov et al. [Titov et al., 2017] compared metabolic networks for the differential diagnosis between AD and frontotemporal lobar degeneration (FTLD). They also proposed an algorithm to estimate if an individual subject shows a more AD or FTLD pattern of regional metabolism. Munilla et al. [Munilla et al., 2017] systematically evaluated the influence of the number of subjects and the regularization strength on the GGM stability and graph structure. They found that the estimated GGM graph structure and small-world coefficient converged to a stable level when including 40 or more subjects in their study sample. For regularization-based approximation of GGMs, they showed that the probability of an edge to exist in the estimated graph structure almost linearly corresponds to the magnitude of their partial correlation. Thus, this finding confirms our initial decision, that sampling-based Bayesian estimation of the graph structure might be more useful for detecting even low associations.

Limitations
It has to be noted that our methodological framework can currently only be applied as a group statistic but not for individual subjects. Therefore, GGMs can be used for exploratory analyses as alternative to Pearson correlation networks, and may aid generating new hypotheses about the interrelation of clinical variables or feature selection. Then, derived hypotheses can be validated using classical statistical methods such as regression or mediation analysis.
Another limitation is the high uncertainty in the statistical model to estimate the partial correlations. This is due to the theoretically hard problem of matrix inversion on the one hand, and due to the high number of possible graph edges in comparison to the sample size on the other hand. Thus, the model might be fragile with respect to the obtained values and requires large training samples to get stable results. Here, we repeated the model estimation on the whole data for ten times to observe the effect on model stability, which was yielding largely consistent results for strong links with high partial correlation, but getting more variable for weaker links with low partial correlation. Replicating the results using the right hemisphere data also yielded largely consistent results with highest agreement for the characteristic path length metric. Apparent deviation in clustering coefficient and consequently in small-world coefficient (=ratio of both) might be explained by the asymmetry of the brain and the lateralization reported for Alzheimer's disease in the literature (e.g. stronger left hippocampus atrophy in ADNI) Teipel, 2016, Weise et al., 2018]. However, our findings still need to be replicated in independent cohorts.
We observed a saturation of the conditional dependency network when adding many variables. This means, the model parameters might strongly change when having only few variables in the model and adding another variable; in contrast to very stable estimates of larger models with dozens of variables, which are hardly altered when adding another variable. Actually, this problem is well-known for linear regression models and related to multicollinearity in the data [O'brien, 2007, Dormann et al., 2013, Teipel et al., 2015b. Recent developments in stochastic block models may help to overcome these limitations, as they try to infer the underlying clustering block structure and separately estimate statistical associations within and between clusters [Sun et al., 2014, Hosseini andLee, 2016].

Conclusion
We applied GGMs to assess inter-modal and inter-regional dependencies of high-dimensional multimodal neuroimaging data of AD-related brain alterations. Our results showed that conditional dependency networks estimated by GGMs provide useful information within imaging modalities and could be used as alternative to Pearson-correlation networks. Nonetheless, GGMs did not detect some expected associations between modalities and, therefore, may have limited applicability for large-scale data with dozens of variables.

Conflict of interest
The authors declare no conflict of interest with respect to this study.

Data availability statement
MRI and PET data being used in this study can be retrieved from ADNI (http://adni.loni.usc.edu/data-samples/access-data/). Processed imaging data and extracted regional mean values are available from the corresponding authors upon request. The R package BDgraph can be downloaded from CRAN (https://cran.r-project.org/web/packages/BDgraph) or GitHub (https://github.com/cran/BDgraph).     The node u is statistically independent from v given the node dis, formally p(u, v|dis) = p(u|dis)p(v|dis). (B) Pearson correlation matrix, showing a "spurious" correlation between nodes u and v. Notably, when considering only u and v alone, the independence assumption does not hold; formally p(u, v) = p(u)p(v). (C) Partial correlation matrix derived from Gaussian graphical models. Using this model, we can approximately recover the underlying dependency structure, with u ⊥ v|dis =⇒ cor(u, v|dis) = 0. Figure 2. Pearson correlation matrix (left) and partial correlation matrix (right) for the three imaging modalities (left hemisphere data only) estimated for the combined data of EMCI, LMCI and AD patients. For better readability, each individual block of the partial correlation matrix is shown in Figures 3-5 and Supplementary Figures S2-S4.
EMCI/LMCI: early and late amnestic mild cognitive impairment, AD: Alzheimer's dementia. Figure 3. Partial correlation matrix for glucose amyloid-β in the left hemisphere estimated for the combined data of EMCI, LMCI and AD patients. Averaged over ten repetitions. Associations of lowest magnitude were not present in all iterations.
EMCI/LMCI: early and late amnestic mild cognitive impairment, AD: Alzheimer's dementia, amy: amyloid-β. Figure 4. Partial correlation matrix for glucose metabolism in the left hemisphere estimated for the combined data of EMCI, LMCI and AD patients. Averaged over ten repetitions. Associations of lowest magnitude were not present in all iterations.
EMCI/LMCI: early and late amnestic mild cognitive impairment, AD: Alzheimer's dementia, metab: glucose metabolism. Figure 5. Partial correlation matrix for gray matter volume in the left hemisphere estimated for the combined data of EMCI, LMCI and AD patients. Averaged over ten repetitions. Associations of lowest magnitude were not present in all iterations.
EMCI/LMCI: early and late amnestic mild cognitive impairment, AD: Alzheimer's dementia, vol: gray matter volume. Figure 6. Partial correlation matrix for amyloid-β in the left hemisphere by group. Averaged over ten repetitions. Associations of lowest magnitude were not present in all iterations.
CN: cognitively healthy elderly controls, EMCI/LMCI: early and late amnestic mild cognitive impairment, AD: Alzheimer's dementia, amy: amyloid-β. Figure 7. Partial correlation matrix for glucose metabolism in the left hemisphere by group. Averaged over ten repetitions. Associations of lowest magnitude were not present in all iterations.
CN: cognitively healthy elderly controls, EMCI/LMCI: early and late amnestic mild cognitive impairment, AD: Alzheimer's dementia, metab: glucose metabolism. Figure 8. Partial correlation matrix for gray matter volume in the left hemisphere by group. Averaged over ten repetitions. Associations of lowest magnitude were not present in all iterations.
CN: cognitively healthy elderly controls, EMCI/LMCI: early and late amnestic mild cognitive impairment, AD: Alzheimer's dementia, vol: gray matter volume. Figure 9. Comparison of graph statistics for the partial correlation matrices of the left hemisphere stratified by diagnostic group and image modality. Estimates based on Gaussian graphical models using multimodal neuroimaging data. The distribution of the weighted clustering coefficient, characteristic weighted path length, and small-world coefficient for individual brain regions is shown. Boxes display median, first and third quartile of the distributions, and whiskers indicate ±1.5×interquartile range. All blocks showed significant differences in mean between groups, one-way analysis of variance (ANOVA), df = 215, F ≥ 4, p < 0.01, η 2 ≥ 0.055. Details are given in Supplementary Tab. S2. P-values for Tukey's honest significant difference tests are given in Tab. 2.
CN: cognitively healthy elderly controls, EMCI/LMCI: early and late amnestic mild cognitive impairment, AD: Alzheimer's dementia, amy: amyloid-β, metab: glucose metabolism, vol: gray matter volume.    Figure 3. Partial correlation matrix for amyloid-β deposition in the left hemisphere estimated for the combined data of EMCI, LMCI and AD patients. Averaged over ten repetitions. Associations of lowest magnitude were not present in all iterations. EMCI/LMCI: early and late amnestic mild cognitive impairment, AD: Alzheimer's dementia, amy: amyloid-β.     Clustering coefficient q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q

Comparison of graph statistics (left hemisphere, 10x repeated)
q q q q q q q q q q q q q q q q q q q q q q q q q q q CN EMCI LMCI AD q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 8 10 12 14 16 18 Characteristic path length q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q  Figure 9. Comparison of graph statistics for the partial correlation matrices of the left hemisphere stratified by diagnostic group and image modality. Estimates based on Gaussian graphical models using multimodal neuroimaging data. The distribution of the weighted clustering coefficient, characteristic weighted path length, and small-world coefficient for individual brain regions is shown. Boxes display median, first and third quartile of the distributions, and whiskers indicate ±1.5×interquartile range. All blocks showed significant differences in mean between groups, one-way analysis of variance (ANOVA), df = 215, F > 4, p < 0.01. P-values for Tukey's honest significant difference tests are given in Tab. 2. CN: cognitively healthy elderly controls, EMCI/LMCI: early and late amnestic mild cognitive impairment, AD: Alzheimer's dementia, amy: amyloid-β, metab: glucose metabolism, vol: gray matter volume.    Figure S2. Partial correlation matrix for amyloid-β deposition and glucose metabolism in the left hemisphere estimated for the combined data of EMCI, LMCI and AD patients. Averaged over ten repetitions. Associations of lowest magnitude were not present in all iterations. EMCI/LMCI: early and late amnestic mild cognitive impairment, AD: Alzheimer's dementia, amy: amyloid-β, metab: glucose metabolism.  Figure S4. Partial correlation matrix for amyloid-β deposition and gray matter volume in the left hemisphere estimated for the combined data of EMCI, LMCI and AD patients. Averaged over ten repetitions. Associations of lowest magnitude were not present in all iterations. EMCI/LMCI: early and late amnestic mild cognitive impairment, AD: Alzheimer's dementia, amy: amyloid-β, vol: gray matter volume.   Figure S5. Comparison of weighted clustering coefficient stratified by brain region, diagnostic group and modality for the partial correlation matrices of the left hemisphere. Averaged over ten repetitions. CN: cognitively healthy elderly controls, EMCI/LMCI: early and late amnestic mild cognitive impairment, AD: Alzheimer's dementia, amy: amyloid-β, metab: glucose metabolism, vol: gray matter volume.   Figure S6. Comparison of characteristic path length stratified by brain region, diagnostic group and modality for the partial correlation matrices of the left hemisphere. Averaged over ten repetitions. CN: cognitively healthy elderly controls, EMCI/LMCI: early and late amnestic mild cognitive impairment, AD: Alzheimer's dementia, amy: amyloid-β, metab: glucose metabolism, vol: gray matter volume.   Figure S7. Comparison of small-world coefficient stratified by brain region, diagnostic group and modality for the partial correlation matrices of the left hemisphere. For better readability, individual values were upscaled by a factor of 1,000. Averaged over ten repetitions. CN: cognitively healthy elderly controls, EMCI/LMCI: early and late amnestic mild cognitive impairment, AD: Alzheimer's dementia, amy: amyloid-β, metab: glucose metabolism, vol: gray matter volume.  Figure S8. Comparison of graph statistics for the partial correlation matrices of the right hemisphere stratified by diagnostic group and image modality. Estimates based on Gaussian graphical models using multimodal neuroimaging data. The distribution of the weighted clustering coefficient, characteristic weighted path length, and small-world coefficient for individual brain regions is shown. Boxes display median, first and third quartile of the distributions, and whiskers indicate ±1.5×interquartile range. CN: cognitively healthy elderly controls, EMCI/LMCI: early and late amnestic mild cognitive impairment, AD: Alzheimer's dementia, amy: amyloid-β, metab: glucose metabolism, vol: gray matter volume.  Figure S9. Comparison of graph statistics for the Pearson correlation matrices of the left hemisphere stratified by diagnostic group and image modality. The distribution of the weighted clustering coefficient, characteristic weighted path length, and small-world coefficient for individual brain regions is shown. Boxes display median, first and third quartile of the distributions, and whiskers indicate ±1.5×interquartile range. Prior to calculating the graph measures, the correlation matrices were thresholded such that correlations with p > 0.05, i.e. approximately r < 0.12, were set to zero. CN: cognitively healthy elderly controls, EMCI/LMCI: early and late amnestic mild cognitive impairment, AD: Alzheimer's dementia, amy: amyloid-β, metab: glucose metabolism, vol: gray matter volume.