In this section we have shown that the relationship of graph properties between health conditions strongly depends on the methods used for graph construction. For our data we find more significant mean differences for control-AD and control-MCI then for MCI-AD. With respect to clustering and edge definition methods, the largest number of significant differences are found for RGS and Ward clustering, and for model-free edge definitions. These results show that conclusions on how graph properties change due to AD have to be drawn carefully, and ideally validated by other methods, as they can be highly sensitive to the methods used for graph construction.

### 2.3. Evaluation of Graph Construction Methods Based on Negative Surprise

Having examined the consequences of particular choices for clustering and edge definition techniques in the previous sections, we now evaluate their combinations by considering their ability to help a clinician to discriminate among patient groups. This discrimination is achieved by using the graph data within a statistical model, which specifies the likelihood of the graph data. The model is described in section 5.6; the likelihood is a distribution which depends on a set of parameters. In general, the kind of graph data—i.e., their construction method—and the statistical model with its parameters are interdependent: they cannot be freely varied separately. Therefore our evaluations of the predictive power of the various graph construction methods have to be understood with a caveat: they depend on our specific choice of statistical model.

To quantify the discriminating power for each graph construction combination, we use a metric based on the final probabilities for the correct health conditions known as the log-probability, or negative surprise (Bartlett, 1952; Good, 1956, 1957a,b, 1983): a sure event, i.e., with unit probability, has surprise equal to zero; whereas an impossible event, i.e., with zero probability, has surprise equal to infinity, reflecting the fact that its occurrence would be contrary to all our expectations. A high surprise (in absolute value) therefore signals a low predictive power of the data we are using. The expectation or average of the surprises is the Shannon entropy (Shannon, 1948; Bartlett, 1952; McCarthy, 1956; Bernardo, 1979; Jaynes, 2003: section 11.3).

Another possibility, of a more decision-theoretical character, is to consider a metric based on the average utilities obtained with each particular graph-construction method. Given several possible courses of action (e.g., treat or dismiss) and their utilities or costs with respect to each health condition (e.g., treating an Alzheimer patient, dismissing a healthy patient, dismissing an Alzheimer patient, or treating a healthy one), the clinician

should choose the action that maximizes the expected utility, the expectation being calculated from the final probabilities for the possible health conditions (Sox et al., 2013). This kind of metric therefore requires not only the final probabilities—which depend on the graph-construction method—but also a table of utilities.

Numerical tests show that the two kinds of metric yield similar results, at least for utility tables close to the identity (treating an ill patient and dismissing a healthy one have unit utility; the remaining combinations have zero utility). We therefore choose a metric based on the negative surprise, which is simpler and more intuitive than a utility metric.

In order to have an approximate idea of the relative predictive powers of the graph-construction methods we would like to use a statistical method that can be kept the same, as much as possible, across different methods. For this reason we choose a model based on the working hypothesis of sufficiency of mean and correlations of past data, as explained in the Introduction. This model ignores any restricted range of variability of graph quantities (e.g., positive or bounded). As explained in Porta Mana et al. (2018), this choice is non-standard but does not entail contradictions. The model has some free parameters; their values reflect the fact that the units of measure for the graph quantities make the latter of order unity. This choice of a generic, common statistical model allows us to sidestep the demanding problem of tailoring it for the different graph quantities from our 850 graph-construction methods.

**Figure 8** shows the obtained negative surprises for all combinations of graph construction methods except H2D, which is left out due to an inadequacy of the statistical model, resulting in unrealistic values between −1.26 and −0.66 with a mean and standard deviation of −0.94 ± 0.19.

The differences in negative surprise between the different graph construction method are in general small. The best results are obtained for ward4 clustering combined with mutual information (MIT) based edge definition. Across edge definition methods, linear correlation (corr) and mutual information give the best results and transfer entropy (TE) the worst. The rather poor performance of TE edge definition is in line with the small number of significant differences found for this method (compare **Figure 6**). Comparing the different clustering methods, atlas and ward4 clustering give the best results, as long as the edge definition is not TE. These two clustering methods have in common a very small number of graph nodes and (correspondingly) the highest number of voxels per cluster (compare **Figure 2**).

As explained above, the comparison of graph-construction methods can be affected by the statistical model and its parameters, especially for small datasets. As a complementary analysis we compare the negative surprises with the classification performances of a support vector machine (SVM, section 5.7.1) based on the same graph constructions. In a clinical setting, a misclassification between control and AD has more severe consequences than between MCI and AD. To avoid introducing an asymmetric misclassification penalty, we perform the classification between pairs of classes only (control-AD, C-MCI, MCI-AD).

**Figure 9** shows the relationship between the SVM performance (measured as proportion of correct classifications) against the negative surprise. As long as TE edge definition is excluded, the two performance measures are positively correlated. In particular RGS clustering achieves low performance in both negative surprise and SVM classification. Furthermore, atlas clustering achieves a high classification performance across all edge definitions. The exact SVM classification results for each realization of graph construction method are depicted in **Figure S2** (see Supplemental Material).

**Figure 10** demonstrates that thresholding graphs has only a minor effect on the negative surprise for small thresholds up to 0.2. No systematic relationship can be observed for the effect of larger thresholds; for example, increasing the threshold to 0.4 causes a decrease in negative surprise for RGS clustering with linear correlations or mutual information, but an increase for atlas clustering with transfer entropy edge detection. Likewise, the creation of highly connected and rich club sub-graphs typically decreases the negative surprise, but in some cases increases it (e.g., RGS H2U). Overall the highest negative surprise (−0.66) is obtained for ward4 clustering combined with BMITU1 thresholded at wmin = 0.1.

These results suggest that the best combination of graph construction techniques to use for this data set is the atlasbased or ward4 clustering combined with linear correlation methods or mutual information transfer. Thresholding the graph edges, which might reduce experimental noise and does lower computational complexity, has only a minor effect on the predictive power, as long as threshold values are small. Reducing the graphs complexity via larger thresholds or extracting the rich-club of the graph should be done with care, since the results can change in either direction. Although transfer entropy yields lower negative surprises then the model-free functional connectivity measures, we would not conclude that this edge definition performs worse in general, since it achieves high values in SVM classification. It is very likely that our choice of statistical

model is not ideal, and a more tailored choice would improve performance.

#### 3. DISCUSSION

In this article we have compared different techniques for constructing and analyzing graphs. By applying a statistical model, we have demonstrated a principled method for choosing a combination of techniques for a given data set. By examining the varied outcomes of the techniques, we have shown how sensitive the results of graph theoretical analyses, such as significant differences in mean properties, can be to the choice of clustering or edge definition technique.

With regards to the predictive power of the graph construction techniques, measured in terms of negative surprise, we find that Ward and atlas clustering yield the highest performance of the clustering techniques, and region growing and selection clustering (RGS) the lowest. In particular, the variant of Ward clustering that produces large clusters and small numbers of nodes (ward4) achieved the highest performance values. Analogously for the edge detection methods, we find better performance for the model-free methods (linear and nonlinear correlations, mutual information transfer) than for the model-based method of transfer entropy. For this particular data set, a combination of ward4 clustering with mutual information derived edges achieves best results. Therefore, we would recommend this combination as the primary target for a more narrowly focussed investigation based on a larger data set.

The performances we obtain are above chance level but still far away from optimal prediction of the three health conditions. One reason for this sub-optimal prediction might lie in our choice of statistical model and its parameters. With our small data set (26 controls, 16 MCI, 14 AD) the model and its parameters have a high influence on the final probabilities, and thus on the performance (Porta Mana et al., 2018). We avoided tailoring the statistical model for the theoretic and practical reasons explained in section 2.3. Even if the model is not tailored, the results are consistent with the classification performance of support vector machines (see **Figure 9** and **Figure S2**), for the model-free edge definition techniques.

It remains unclear why Ward and atlas clustering are more successful than RGS, especially in combination with model-free edge definition. One possibility is that this is related to the large variability in graph sizes generated by RGS (**Figure 2**). In addition, the variance of weight distributions across subjects, and the variance of the cluster distances, are much larger in RGS then in Ward clustering (**Figure 3**). This could be related to the variance in the number of nodes; however, choosing graphs similar in size causes even higher variances (section 2.1). Therefore we assume that the number and connectivity of the small functional units extracted by RGS are highly variable across subjects. This variance might be even higher across subjects within a health condition than across health conditions, such that changes due to AD cannot be detected. This assumption might at first glance seem to contradict the high number of significant comparisons observed (**Figure 5**). However, we only calculate the significance level for the means of the distributions and not their entire shape. In addition, it is likely that some graph properties correlate with the graph size, and thus that apparent significant differences in graph properties are simply reflecting significant differences in numbers of nodes detected, and do not provide further information useful for classification or understanding the nature of the disease. Further investigation is needed on this matter.

The low negative surprise of transfer entropy (TE) compared with other model-free functional connectivity measures might

FIGURE 10 | Negative surprise for different graph edge thresholds wmin (wmin = 0 for complete graphs, indicated by a vertical dashed line) and rich club graphs (rich) for different edge definitions: corr (First Panel), H2U (Second Panel), MIT (Third Panel), TE (Fourth Panel) and different clustering methods Ward (light gray), RGS (dark gray) and atlas (black). Each dot is the result of averaging across all possible parameters of a general graph construction method (for wmin = 0 the average across all points of a swarm in Figure 8). Since some methods yield small edge weights, some graphs become unconnected for large wmin such that the statistical analysis is not conducted; no values are depicted in this case. Markers are connected for better visual comprehension.

have several reasons. The comparison of the negative surprise with the support vector machine classification suggests that a better choice of a statistical model is possible: the classification results for TE are similar to those of the model-free measures. In TE the data of a certain time interval in the past is used in order to calculate how much the uncertainty of the future is reduced. Here we use the data of the last 15 s. This time period might be poorly chosen, influencing the overall negative surprise. In addition TE is more sensitive to short recording periods than other methods, which may well also result in a reduced performance (Pereda et al., 2005).

With regards to the robustness of the graph theoretical outcomes, we discovered that relationships between mean graph properties, such as closeness centrality, edge weight or clustering coefficient (**Figures 4**–**6**) were sensitive to choice of clustering and edge definition techniques, to parameter choices for a given technique, and to the manner in which sub-graphs were defined (thresholding value and rich club). For most relationships between graph properties X, we could find significant (p < 0.05) differences in both directions, i.e., both XAD > XC and XAD < XC, for specific choices of clustering and edge definition technique. This strongly suggests that a degree of suspicion should be applied to studies reporting such significant differences, especially if these results are argued to give insight into how a disease affects brain properties, unless the significance level is much more compelling or the reported differences can be validated with alternate methods.

We also investigated the sensitivity to method choice of the ordering of subjects according to a graph theoretic metric (**Figure 7**). In this analysis, transfer entropy was the most consistent. Nevertheless, the distributions of the negative surprises is as broad for transfer entropy as for other edge definitions (**Figure 8**). In general, the exact parameter selection within an edge definition method causes only slight changes in the negative surprise, more crucial is the exact realization of the clustering method: ward4 clustering generally achieves a better performance then ward3 clustering. These two variants differ only in the number of predefined clusters (see Supplemental Material **Figure S1**). Applying a lower threshold wmin on the graph's edge weights has little effect on the negative surprise for all methods, as long as only small weights (up to 0.2) are set to zero. Thresholding higher weights or extracting the graph's rich club has unpredictable effects on the results, and so should be used with caution (**Figure 10**). Atlas clustering was least consistent in the subject ordering analysis, suggesting that although it may provide a good basis for a diagnostic tool, care should be taken in reporting discoveries of particular relationships in graph properties between health conditions, as these may well turn out to be critically dependent on the edge definition method used.

Due to the intense computational requirements of the survey performed in this article, we recognize that it would be advantageous to develop heuristics for choosing between graph construction methods without performing the full calculation for each combination. Our results suggest that properties visible at the clustering stage, such as average heterogeneity, may give some indication of predictive performance: graph constructions that result in different degrees of heterogeneity between the health conditions seem to be more discriminable by the later steps of the calculation. More research is needed in this area, which is outside the scope of the current study. In addition, it is tempting to consider t-test results of the mean graph properties as a heuristic. Our results suggest that this approach is largely inadequate. It holds for edge definition via transfer entropy, which gives very few significant results and the negative surprise is rather small compared with the model-free edge definitions. Conversely, region growing clustering yields most significant differences but a generally poor negative surprise. This may be due to graph properties being highly correlated, and so not providing additional information to the statistical model. In addition we used the first four moments (wherever possible) in our statistical model, rather than just the mean, which may also partially account for this apparent contradiction.

In addition to considering the predictive power and robustness of graph construction techniques, we can also evaluate them according to their practicality, i.e., speed of calculation and the extent to which they are easily available in established medical infrastructure and diagnostics. In general, applying graph theoretic measures to fMRI data for improving AD diagnosis makes sense, since MRI scans are already implemented in AD diagnostics for detecting structural changes such as hippocampal dystrophy caused by AD or AD-unrelated pathology (e.g., brain tumors). Softwares such as SPM (Tzourio-Mazoyer et al., 2002) and FSL (Jenkinson et al., 2012) are frequently used in medical research and mainly support clustering that is atlas and independent component analysis based. Ward clustering, which is the fastest of all these clustering methods, is a standard hierarchical clustering method and implemented in all standard programming softwares such as Python and Matlab. The region growing algorithm is not implemented in established softwares and is also computational very demanding. Given that it does not out-perform atlas or Ward clustering, we therefore do not recommend it. For edge definition and graph properties, several software packages are available based on Matlab (Wang et al., 2014; Kruschwitz et al., 2015) or Python1 , which provide a comprehensive range of edge definition and graph analysis methods.

In general we recommend using statistical models and not pure classifiers such as support vector machines as diagnostic tools, since statistical models calculate a probability of a diagnosis rather than assign a classification, i.e., "Given the fMRI scan, person x has a 80% probability of having Alzheimer's disease," rather than "Given the fMRI scan, person x has Alzheimer's disease." Probabilities can be easily combined with other probabilities of other diagnostic tests (Porta Mana et al., 2018) such as cognitive assessment, amyloid beta and tau protein occurrence in cerebrospinal fluid, blood tests, and structural MRI2 (Johnson et al., 2012). This allows the medical doctor to conclude, for example: "Given the results of the cognitive test and cerebrospinal fluid analysis and structural and functional MRI scan, person x has a 95% probability of having Alzheimer's disease." After the estimation of the probability for a disease, she has to decide on a treatment, also taking into consideration such factors as "how harmful would the treatment be for a healthy person," which can be expressed in a utility function (Porta Mana et al., 2018). In addition, the statistical model used in this work allows an estimation of how much the model can be trusted, and therefore evaluate whether the sample size is sufficiently large (Porta Mana et al., 2018).

#### 3.1. Relationship to Previous Studies

Studies focusing on the graph properties extracted from restingstate fMRI in AD and its pre-stages generally have one of two aims. The first aim is to identify significant differences in the graph properties between health conditions, and to use these to gain insight into the effects of AD on the physical brain and its cognitive processes. These studies complement the picture revealed by investigations based on structural MRI and functional changes on the basis of EEG and MEG recordings. Typically a variety of graph properties (e.g., nodal degree, clustering coefficient, averaged shortest path, local efficiency, betweenness centrality, global efficiency, small worldness) are calculated, and used to motivate an account of how diseaserelated modifications to these properties result in a reduced capacity to transfer and process information.

However, such studies reveal entirely contradictory results. For example, the value of the clustering coefficient in AD with respect to controls has been reported to be increased, unchanged, and decreased, respectively (Supekar et al., 2008; Sanz-Arigita et al., 2010; Zhao et al., 2012). Analogous contradictions have been found for the comparative length of the shortest path (Supekar et al., 2008; Sanz-Arigita et al., 2010; Zhao et al., 2012). These contradictions could be caused by methodological differences or by not separating the different states of AD. Our results show ample evidence that the precise choice of graph construction techniques can easily account for contradictory findings, even for atlas based clustering, in which the number and size of clusters is held constant across all subjects (**Figure 5**). Evidence that the separation of different AD stages is relevant was provided by Kim et al. (2015), who demonstrated a nonmonotonic behavior of global efficiency, local efficiency and betweenness centrality across different stages of AD and MCI. In our study, we could reproduce the pattern of increase and decrease of closeness centrality across conditions (**Figure 4**). However, we also demonstrate that the same analysis based on the rich club sub-graph yields a different pattern, and that contradictory (but significant) results can be obtained for the same graph construction techniques with different choices of threshold. We thus conclude that differences in graph properties between health conditions are currently ill-suited to provide an account of disease mechanisms in AD, unless either: (1) a specific method of graph construction can be shown to be more representative of the underlying connectivity than other methods, (2) the differences can be shown to be robust to choice of graph construction, (3) the differences can be validated by another analytical approach, or (4) the significance level is shown to be substantially more persuasive than p < 0.05.

1https://github.com/dpisner453/PyNets.

2https://www.alz.org/research/diagnostic\_criteria/.

The second category of studies use graph theoretical information as input for machine learning algorithms to classify the health conditions of the subjects. Note that for this purpose it is irrelevant if a difference between health conditions is not robust to method choice, as the goal is not to understand the effects of the disease but to robustly distinguish between conditions. Recent studies have reached very high performance: 100% accuracy in discriminating AD and control (Khazaee et al., 2015), and 93% for AD, MCI and control classification (Khazaee et al., 2017). In the latter work they extract more than two dozen local and global graph properties, resulting in roughly 3, 000 features, since each of the local properties is calculated for all brain areas. Only a small subset of features is then used for classification, e.g., in-degree of the left middle temporal gyrus. They found that the classification power of local graph measures is larger than that of the global ones. Local changes in graph properties that do not propagate to global mean values have also been reported for area specific (frontal cortices, parietal and occipital regions) synchronization levels (Sanz-Arigita et al., 2010).

In this work we do not compare node-specific graph properties, because Ward and RGS clustering do not result in the same spatial location of clusters across subjects. Instead, we consider, wherever possible, the first four moments of the entire distributions of graph properties. This is more information than typically used for global measures, where often only the first moment (the mean) of a graph property distribution is taken into consideration. Nevertheless, it is still possible that considering single nodes, of which some may be more damaged by AD than others, could yield a better diagnostic performance. This requires further study in a survey considering only atlas based clustering. Again, this is out of scope of the current study, but we remark that the statistical model methodology we employ here would be equally applicable to such an investigation. The advantage of taking the entire distribution lies in the possibility of using purely data driven clustering algorithms (e.g., Ward clustering) that can be substantially faster than atlas based clustering, since they do not depend on a time and memory consuming registration of the individual brain image to standard space.  Finally, a short recording time might be expected to have a weaker influence on entire graph property distributions then on single nodes. Thus we conclude that global measures are preferable, if a good diagnostic performance can be reached. Although the goal of this work was not classification, we note that we obtain up to (80–90%) correct classification using an offthe-shelf support vector machine on leave-one-out subsets of our data for pairwise (C-AD, C-MCI, AD-MCI) comparisons. Whether global measures can reach the impressive performance shown by Khazaee et al. (2017) can only be investigated on a sufficiently large data set, ideally with several hundred participants.

### 3.2. Limitations of This Study

In each step of the graph construction and analysis pipeline (**Figure 1**) we set limits to the endless space of possible methods and their corresponding parameters. Here we will shortly summarize the reasons motivating the selection of the methods examined here and the exclusion of others, given the constraint of limited computational and temporal resources. As a general principle, we aimed to include the most commonly used method(s) and additional methods that we found to be reasonable, even if they are not currently frequently used.

 The global signal (the average activity across all brain voxels) is assumed to originate partly from vascular and respiratory processes that do not represent neuronal activity. However, there is also evidence that it contains neuronal-signaling based components, since it is negatively correlated with the EEG signal and strongly correlated with the activity of the largest network in the brain (the default mode network, which plays a major role in rest state activity) when noise levels are low (Murphy and Fox, 2017). Without global signal regression, the Pearson correlation distribution derived from the signal of all voxels, or the average activity of clustered voxels, is biased to the right such that negative values are rare and small. The correction for the global signal centers this distribution, such that negative values are much more prominent. This also changes the properties of the graphs extracted from such data, for example an increase in modularity combined with fewer unconnected nodes has been reported (Schwarz and McGonigle, 2011; Hayasaka, 2013).

Speaking against global signal regression is the finding that correction for white matter, CSF and motions yield the most stable graph properties across sessions compared with additional applied global regression (Schwarz and McGonigle, 2011). In diagnostics it is important to have only small variance in the outcome across different sessions if the health condition of a subject is stable, such that small changes that indicate a worsening of the health condition can be rapidly detected. Moreover, we define the edges of our graphs as the absolute values of the functional connectivity values. As the negative part of the correlation distribution is small without global regression, different possible treatment of negative correlations (taking the absolute values or setting them to zero) should have only a small influence on the resulting graph properties, at least when the underlying functional connectivity are based on correlations. Consequently, we elect not to include global signal regression in our pipeline.

In the clustering step, the most commonly used method is to define clusters based on cortical regions defined by a brain atlas. We supplemented this with two data-driven clustering approaches: Ward clustering and RGS clustering. We selected Ward clustering, as it has been shown to perform better than alternative hierarchical clustering methods with respect to reproducibility and accuracy (Thirion et al., 2014). RGS, a method derived from image processing (Lu et al., 2003), was selected because we could adjust the method to produce functionally homogeneous clusters. In this formulation, the only free parameter of the algorithm is the minimal cluster size. For both data-driven methods, we selected parameters such that graphs did not exceed a maximal size of 1,500 nodes, due to computational limitations. We excluded clustering based on independent component analysis, because of its laborious implementation and the requirement for domain expertise to distinguish noise from activity-related components. We also excluded all clustering algorithms that do not take functional consistency into account, e.g., dividing the voxels into cuboid patches, as has been proposed for structural data (Amoroso et al., 2017).

With regards to methods for edge definition, we limit our survey to functional connectivity measures that act in the time domain and not in the frequency domain, thus omitting frequency based wavelet analysis (Supekar et al., 2008), synchronization likelihood (Sanz-Arigita et al., 2010) and coherence (Wang et al., 2014). The most commonly used and simplest functional connectivity measure is the Pearson correlation coefficient (e.g., Zhao et al., 2012), which we name BCorrU in our work. We also test two additional model-free and one model-based method. A further model-based method based on Granger causality was excluded because it is too computationally expensive for larger graphs (Wang et al., 2014).

A thresholding operation is often applied to graphs extracted from fMRI, setting all values below wmin to zero. The aim of this step is to reduce experimental noise, which mainly manifests in the weaker edges, and to make the computation of graph properties computationally less demanding (Bordier et al., 2017). The threshold wmin can be defined in several ways: it can be set arbitrarily, without satisfying a certain demand, or such that certain properties of the graphs are preserved, e.g., average number of edges per vertex (Sanz-Arigita et al., 2010), node density (Zhao et al., 2012), small world behavior (Bassett et al., 2008) or a fixed cluster coefficient. Alternatively, it can be set such that information on the network's community structure is maximized; see, e.g., Bordier et al. (2017). In a variant of the thresholding approach, it has been proposed to transform the edge weights by applying a power law (Schwarz and McGonigle, 2011). In this study, for the sake of simplicity, we examine graph properties as a function of wmin without targeting any specific value of a graph property. Potentially, our results would reveal a different picture if wmin was optimized for each subject to attain, for example, a specific average nodal degree. However, comparison of these two different thresholding mechanisms resulted in no major difference in the relationships of graph properties between the control and AD groups (Sanz-Arigita et al., 2010).

We do not binarize our graphs (setting all values below wmin to zero and those above it to one) as is frequently done (e.g., Zhao et al., 2012), as this leads to a loss of information, and moreover some distributions of graph properties would become discrete (e.g., only ones and zeros for edge weights distributions), such that higher moments would be uninformative. The disadvantage of using weighted graphs lies in the limitation of possible graph properties. Most graph properties are well-defined for binary graphs and have been partly extended to weighted graphs. Here, we calculate the (normalized) weighted degree, shortest path, closeness centrality, clustering coefficient, and the modularity. We only investigate the most commonly used metrics and do not include more complex methods such as the minimal spanning tree (Çiftçi, 2011).

In addition to the restrictions of scope with regards to the examined techniques, a clear limitation of this study is the small data set. As our aim here is primarily to provide a methodology for evaluating and comparing analysis methods, rather than to draw conclusions on the effect of Alzheimer's disease on the graph properties of the cortex, a small data set is less problematic. Indeed, for the explorative survey carried out here, a large data set would have been prohibitively expensive with respect to computational resources. Moreover, many studies applying graph analysis to fMRI data are based on similarly sized data sets, which highlights the importance of raising awareness about the methodological artifacts we have identified.

The results of our survey indicate which combinations of methods are promising in view of Alzheimer diagnosis and should be investigated further in future studies based on larger data sets. Naturally, such studies could yield some quantitatively different results to those reported here, particularly with regard to the classification performance. Nonetheless, we would like to summarize some conclusions of the work that are unlikely to change with a larger data set. First, our results show that different combinations of methods can lead to contradictory findings with regard to significant differences in mean properties (section 2.2). This effect is unlikely to be resolved by a larger sample size. Second, methods showing good robustness with respect to parameter choice for a small sample size (e.g., TE edge definition, see **Figure 7**), are likely to remain robust with increasing sample size. Likewise, there is no reason to assume that methods performing well in all circumstances for the small data set, e.g., Ward clustering combined with corr edge definition (section 2.3), would perform worse for larger data sets. Finally, we assert that thresholding the graphs of a large data set with a small wmin (as shown in section 2.3) would similarly not result in a sudden jump in negative surprise.

### 3.3. Application of Approach to Other Analysis Techniques

We have demonstrated a systematic, quantitative approach for comparing and evaluating sequences of algorithms that result in classification of fMRI data based on the first four moments of simple graph theoretic metrics defined on the whole graph. However, the approach we present is equally well suited for assessing pipelines based on other metrics, as we briefly outline in the following.

One possibility is to consider the graph properties of individual nodes, as these have been shown to be very informative (Xia et al., 2014; Khazaee et al., 2015; Wang et al., 2016; Dillen et al., 2017).

This entails the use of atlas based clustering. We speculate that a global analysis of graph properties would be both faster and more robust to brain abnormalities and short recording times, and so would be the preferable approach if equivalent performance levels can be attained.

A second possibility is to extend our approach to a hierarchical analysis. This could potentially be of great use, as previous studies based on PET imaging have suggested that in Alzheimer's disease, long range connections become weaker but local clustering increases (Pagani et al., 2016, 2017). These alterations would not be observable using the graph analyses so far considered, although we have taken the first step by calculating the modularity, which compares the ideal dissection of the given graph into modules with that of a random graph with similar edge weights.

To capture the graph meta-structures it is necessary to cluster graph nodes into modules, or sub-graphs. Modules can be defined either purely functionally, such that each node (ideally) has the strongest connections to the nodes in its own cluster, and the weakest connections to nodes of other clusters, or based on anatomic structures, such that nodes in a cluster are part of large, anatomo-functionally similar brain areas. Analogous to the variety of methods for spatial clustering and edge definition investigated in this study, there are many techniques used to cluster nodes into modules (e.g., k-clustering, hierarchical clustering and spectral clustering, for a review see Schaeffer, 2007 or anatomo-functional clustering, see Pagani et al., 2016), and likewise multiple options for analysing the characteristics of the resulting modular structure (e.g., module degree or participation coefficient; see Guimerá and Nunes Amaral, 2005). Such a comprehensive study is outside the scope of the current work, but could well provide great insight into health condition related alterations in the global network structure of the brain.

### 4. CONCLUSIONS

In order to achieve a robust and successful Alzheimer's disease diagnosis based on graphs extracted from fMRI data, we recommend clustering that results in rather small graphs with large clusters. Ward clustering, in which the number of clusters can be predefined, is fast, but requires programming knowledge to implement it. Atlas clustering is well established standard fMRI analysis software applications, but it is slow and might be affected by morphologic abnormalities in the brain, such as atrophy which is a common symptom of AD.

Edge weights should be calculated based on correlations or mutually information transfer, especially if a focus of the study is uncovering significant differences in mean graph properties between health conditions. We emphasize that the existence, magnitude and direction of such significant differences can be very sensitive to the methods chosen, and the parameterization of those methods, and so such findings should be reported with care, especially if a biological interpretation of said findings is claimed. Transfer entropy rarely gives significant results, but is more robust toward parameter changes in the algorithm and different clustering algorithms. Finding appropriate statistical models may be an additional challenge for this method.

Weak thresholding may be used for complexity reduction as it has little effect on performance. Applying a higher threshold or extracting the rich club sub-graph (The 10% of nodes with highest degree) causes unsystematic changes in the negative surprise and should therefore be used with caution, and validated against the full graph.

In summary, our quantitative evaluation and comparison of graph construction and analysis methods provides insight into how contradicting results come about in studies of graph properties of fMRI data, and identifies a number of potential methodological artifacts. Moreover, it provides a blueprint for establishing appropriate analysis pipelines, and serves as a well-founded starting point for future research on larger data sets.

### 5. METHODS

#### 5.1. Data Acquisition

The recruitment and neuropsychological assessment of the study participants is given in Dillen et al. (2017). Demographic information is given in **Table 2**.

Anatomical MRI and resting state fMRI (rfMRI) images were obtained from a 3T MR-Brain-PET scanner (Siemens, Erlangen, Germany) in the Memory Clinic Cologne Juelich. The parameters for the single-shot echo planar imaging sequence of the functional (T2\* weighted) image are the following: TR = 3, 000 ms, TE = 30 ms, FA = 90◦ , FOV = 200 × 200 mm2 , matrix = 80 × 80, voxel resolution = 2.5 × 2.5 × 2.8, 50 oblique slices parallel to the infra-supratentorial line, gap = 0.28 mm, interleaved, scan time = 7 min. Parameters of the high-resolution T1-weighted structural image based on a magnetization-prepared rapid gradient echo sequence: TR = 2, 250 ms, TE = 3.03 ms, FA = 9 ◦ , FOV = 256 × 256 mm2 , matrix = 256 × 256, voxel resolution = 1 mm isotropic, 176 sagittal slices, no gap, interleaved, scan time = 314 s. For more detail see Dillen et al. (2017).

#### 5.2. Preprocessing of fMRI-Data and Extraction of Cortical Data

Image preprocessing is accomplished using FMRIB's Software Library tools (FSL; Woolrich et al., 2009; Jenkinson et al., 2012). We carry out the following steps for the structural T1-weighted image: skull-stripping (Smith, 2002) with bias field correction (Keihaninejad et al., 2010; Leung et al., 2011; Popescu et al., 2012) and for the functional T2-weighted image: discarding the first 10 volumes (out of 140 each taken after 3 sec), motion

Average and minimal and maximal values [min, max] are given for age and years of education; female (f), male (m).

correction (Beckmann and Smith, 2004), spatial smoothing using a 4 mm full width at half maximum Gaussian kernel, highpass temporal filtering at 0.02 Hz and a six-parameter, rigidbody linear transformation procedure in MCFLIRT (Jenkinson et al., 2002). More details can be found in Dillen et al. (2017), where the same preprocessing is applied. In addition we carry out white matter and cerebrospinal fluid regression (FSL regfilt, MELODIC) to the functional image in order to reduce noise.

In order to extract only cortical voxels from the entire brain fMRI image, as needed for the data-driven clustering described in the next section, we first register cortical regions (frontal- , occipital-, temporal-, and insular-cortex) defined in the MNI structural atlas (Collins et al., 1995) to the structural and then to the functional space. For this registration we apply the transformation matrix obtained from registering the entire standard brain first to the individual structural brain (linear registration with FSL/FLIRT; Jenkinson and Smith, 2001; Jenkinson et al., 2002) and then to the functional space (nonlinear registration with Advanced Normalization Tools, ANTs; Avants et al., 2011). In order to extract only gray matter tissue, we apply the gray matter image of the structural space (segmentation with FSL-FAST; Zhang et al., 2001) registered to functional space as described above, as a mask to the to the functional image.

#### 5.3. Data-Driven and Atlas Based Clustering of Cortical Voxels

In order to construct graphs we cluster cortical voxels into regions using three different methods. Two of these methods, the Ward clustering and the region growing and selection algorithm (RGS) are data driven, such that only neighboring voxels with similar activity are combined into a single region. For these algorithms the number of regions per brain and the participating voxels in a region can differ for each individual and strongly depend on predefined algorithm-specific parameters. The atlasbased cluster algorithm, in contrast, produces the same number of clusters and a constant number of voxels per region across individuals, because the individual brains are mapped onto a standard brain.

#### 5.3.1. Atlas-Based Clustering

For each subject we linearly register the rfMRI image first to the structural, skull-removed image (image segmentation for skull removing with SPM8, Wellcome Department of Cognitive Neurology, London, UKFSL; linear registration with FSL/FLIRT; Jenkinson and Smith, 2001; Jenkinson et al., 2002) and then, through a non-linear mapping, to the MNI standard brain [nonlinear registration with Advanced Normalization Tools (ANTs; Avants et al., 2011); MNI 152 standard brain, non-linear 6th generation (Grabner et al., 2006)]. Regions of interest (ROIs) of the resulting functional image in standard space are extracted such that they match the 94 regions identified by the Oxford lateral cortical atlas (regions have a probability above 50%) (Desikan et al., 2006). A demonstration of how the brain is clustered according to the brain areas is given in the first panel of **Figure 12**.

TABLE 3 | Parameters used for the different clustering algorithms.

Ward clustering (ward), region growing and selection (RGS), atlas-based clustering (atlas).

#### 5.3.2. Ward Clustering

Ward clustering (Python: sklearn.cluster.AgglomerativeClustering, Pedregosa et al., 2011) is a data-driven clustering algorithm, which is initiated by defining each voxel as a cluster and then, in each iteration step, merging the two neighboring clusters (even of different sizes) that after merging show minimal intra-cluster variance compared with all other possible variations of combining two adjacent clusters. In this way, the number of clusters is reduced by one in each iteration step. In our case the clustering stops after k clusters (**Table 3**) are formed. Afterwards, we discard away all clusters that contain less then p voxels (**Table 3**). An example of the outcome of Ward clustering algorithm is depicted in the second panel of **Figure 12**.

#### 5.3.3. Region Growing and Selection

The region growing and selection algorithm is a modified version of the algorithm described in Lu et al. (2003). Region growing implies that each voxel serves as an initial seed (center) and neighboring voxels are added iteratively if they fulfill a certain growing criteria. (**Figure 11A**) The condition proposed for adding a voxel to a region is based on the Pearson correlation coefficient R between the averaged time-varying signals of the pre-merged region and the signal of the voxel to be tested (Lu et al., 2003). If this correlation is higher then a pre-defined threshold T (**Table 3**), the voxel is merged to the region. We tighten the growth criteria by imposing a second condition that allows the merging of voxels only if, in addition to exceeding the correlation threshold, the resulting cluster is also functionally homogeneous. Here, functional homogeneity means that the time-varying signals of all voxels can be expressed as instances of a single signal with varying levels of noise. The number of independent signals in a cluster can be estimated by the spatial functional heterogeneity h (Marrelec and Fransson, 2011):

$$h = n\_0 + \frac{e\_{n\_0} - b\_{n\_0}}{(e\_{n\_0} - e\_{n\_0 + 1}) - (b\_{n\_0} - b\_{n\_0 + 1})},\tag{1}$$

where en are the eigenvalues of the N**x**N covariance matrix of all N time varing signals in a cluster that exceed the eigenvalues generated by the broken-stick model bn, such that en > bn = PN i=n 1/i . The index n0 accounts for the smallest eigenvalues that fulfill this inequality equation, such that

FIGURE 11 | Region growing and selection algorithm. (A) Region growing, left: each voxel (colored squares) serves as center for a cluster, right: example of a growing region (purple), only adjacent voxels that fulfill the fusion criteria are added to the growing cluster. (B) Region selection. Small regions (pink) with centers overlapping with larger regions (green) get deleted (from left to right) in a iterative manner. Remaining regions can still overlap as long as their centers do not cover other regions. This illustration is in 2D for simplicity, the algorithm used for fMRI data acts in 3D following the same rules.

en0 > bn0 and en0+1 < bn0+1 . A value of h = 1 indicates a homogeneous cluster.

The region selection algorithm iteratively selects the largest region and deletes all clusters that have their centers in that region, excluding the possibility that centers overlap with other regions. However, clusters can still overlap (**Figure 11B**). Applying this framework does not guarantee that clusters remain spatially connected after deleting regions with overlapping centers. Nevertheless, a check for spatial consistency reveals that only a negligible fraction of the clusters are disrupted in that way. Finally, we took only the clusters that comprised a minimum number of voxels p (**Table 3**). The outcome of RGS is illustrated in the last panel of **Figure 12**.

#### 5.4. Edge Definition

A graph consists of nodes (vertices) that are connected through edges, that might be weighted or binary and directed or undirected. We construct individual brain graphs by defining nodes that represent clusters as described in section 5.3, such that the mean activity of a cluster becomes a node attribute. We presume that all graphs are fully connected and edge weights are defined in terms of functional connectivity. Since functional connectivity can be calculated in several ways, we apply a range of different connectivity measures. In Wang et al. (2014) many such methods are evaluated, taking the structural connectivity of a toy model as reference. As a starting point, for each proposed category of functional connectivity, measured in time, we select the analysis measurement that captures structural connectivity best. We follow this strategy for all proposed measurement categories in Wang et al. (2014), leaving out only Granger causality measures, due to limited computational resources. We thus use linear and non-linear correlation (corr and H2) and mutual information transfer (MIT) for the model-free category and transfer entropy (TE) for the model-based category. In all groups the bivariate methods perform better then the partial ones. In conclusion we select for each of the families the bivariate implementation that can be both directed and undirected. For consistency we use the same abbreviations for the different methods as in Wang et al. (2014) and the same Matlab toolbox Mulan3 which they made public. Here we provide only a short description of the applied methods and more details can be inferred from Wang et al. (2014).

Linear correlation (corr) are measured based on the Pearson correlation coefficient (Rodgers and Nicewander, 1988) in a pairwise manner. For directed connectivity (BCorrD) delays of up to 5 time steps (**Table 4**) are considered and the largest connectivity value is selected. We do not take into account time lags for undirected correlation (BCorrU).

Non-linear correlations (H2) are based on piece-wise linear correlations of two time signals on which the non-linear curve is fitted (da Silva et al., 1989). Bivariate directed (BH2D) and bivariate undirected (BH2U) are defined as above for linear correlations.

Mutual information indicates how much information is shared between two time varying signals by means of Shannon entropy (Grassberger et al., 1991). For BMITD1 individual histograms of two time series are contrasted to the joint histogram across different time delays. No delays are taken into account in BMITU.

Transfer entropy (Schreiber, 2000) describes how far in the past the activity of a node can reduce the uncertainty of the future activity of another node for which the past activity is also considered. Bivariate directed (BTED, Chicharro, 2011) and bivariate undirected (BTEU) are defined as above for linear correlations.

All methods were tested for a window size that comprises the whole time range (130 time points/6.5 min) and for a sliding window of 50 time points (2.5 min) with an overlap of 10 time points (0.5 min), see **Table 4**. If the methods revealed negative weights, the absolute value was considered. The resulting graphs are directed or undirected weighted graphs with values between TABLE 4 | Parameters of the different functional connectivity measures.

Bivariate (B), undirected (U), directed (D), linear correlation (corr), non-linear correlation (H2U), mutual information entropy (MIT), transfer entropy (TE).

zero and one for all methods except non-linear correlations, where values can exceed one.

Many studies transfer weighted graphs into binary ones by setting all values below a threshold wmin to zero and above to one e.g., Zhao et al. (2012). Following this strategy we also investigate the effect of setting all weights below wmin to zero but leaving higher weights unchanged. As far as the remaining graphs are still connected (left panels in **Figure 13**) and single nodes are not disconnected from the network (right panels in **Figure 13**) we study the disease diagnosis capacity for wmin ∈ {0.1, 0.2, ...0.7, 0.8}. In addition we extract the rich club of the graphs. The rich club is a subgraph that comprises the nodes that are most strongly connected to the network. In this work we define the rich club as the 10% of nodes with highest degree.

#### 5.5. Graph Properties

This section describes the different graph properties that are either characteristics of single nodes (weighted degree, closeness centrality, cluster coefficient), of pairs of nodes (shortest path) or of the entire network (modularity). In the first two cases we get a range of values for each graph. Since we do not know, which are the important features of the resulting distributions, we take the first four moments for our statistical analysis. Because graphs based on data-driven clustering contain different number of nodes and the calculated graph properties might be dependent on the number of nodes, we also include the number of nodes in the subsequent analysis (section 5.6).

#### 5.5.1. Weighted Degree

The weighted degree degw describes how strongly a node is connected to all other vertices of the network, obeying the

Bachmann et al. Exploring fMRI-Graphs for AD Diagnosis

3https://github.com/HuifangWang/MULAN.

Left) and larger then 0.9 (wmin > 0.9, Right). In the according weight histograms (Lower Panel) the red bars correspond to the edges drawn in the upper graph. Edges corresponding to the black bars are not shown.

equation:

$$\deg\_{\mathcal{W}}(\nu) = \sum\_{\boldsymbol{\nu} \in V \backslash \{\boldsymbol{\nu}\}} \mathcal{W}\_{\boldsymbol{\nu}\boldsymbol{\nu}} \tag{2}$$

where wuv is the weight on the edge between nodes u and v of all nodes V in the graph. This definition implies a high dependency of the weighted degree on the number of nodes in a graph. To address this problem, we normalize the weighted degree

$$\deg\_n(\nu) = \frac{\deg\_n(\nu)}{\deg(\nu) \cdot \nu\_{\max}} \tag{3}$$

with wmax being the maximal weight of the graph. The resulting values are between 0 and 1.

#### 5.5.2. Shortest Path and Closeness Centrality

The shortest path distw(u, v) between a pair of nodes u and v describes the path that minimizes the sum of the weights of its participating edges. A small shortest path should indicate a strong functional connectivity, therefore we consider the inverse of the graph weights for its calculation. Its computation is carried out using Dijkstra's algorithm (Rivest et al., 2000), which requires the weights to be positive.

Based on the shortest paths of a network we calculate closeness centrality Cw(v) - a measure that indicates how strongly a node v participates in all shortest paths of the graph. It is given by:

$$C\_{\mathbf{w}}(\nu) = \frac{n - 1}{\sum\_{\mathbf{u} \in V \backslash \{\nu\}} \text{dist}\_{\mathbf{w}}(\boldsymbol{\mu}, \nu)}\tag{4}$$

Here, n is the number of all nodes V in the graph.

#### 5.5.3. Clustering Coefficient

The clustering coefficient cc(v) describes to what degree the neighbors of a node v are connected among each other and with node v. Since our network is weighted, we use the Zhang-Horvath clustering coefficient (Zhang and Horvath, 2005; Kalna and Higham, 2007), which is an extension to the "standard" algorithm applied to binary graphs:

$$\mathcal{cc}(\boldsymbol{\nu}) = \frac{\sum\_{i \neq \boldsymbol{\nu}} \sum\_{j \neq i, j \neq \boldsymbol{\nu}} \hat{\boldsymbol{\nu}}\_{\boldsymbol{\nu} i} \hat{\boldsymbol{\nu}}\_{\boldsymbol{\nu} i} \hat{\boldsymbol{\nu}}\_{ij} \hat{\boldsymbol{\nu}}\_{j\boldsymbol{\nu}}}{\left(\sum\_{i \neq \boldsymbol{\nu}} \hat{\boldsymbol{\nu}}\_{\boldsymbol{\nu} i}\right) \left(\sum\_{i \neq \boldsymbol{\nu}} \hat{\boldsymbol{\nu}}\_{\boldsymbol{\nu} i}^2\right)} \tag{5}$$

for i, j neighbors of v and wˆ denoting the weights normalized by the highest weight in the network, such that 0 ≤ ˆw ≤ 1.

#### 5.5.4. Modularity

A graph can be partitioned into smaller components. Modularity measures the deviation of the properties of these components as compared to the components of a random graph with the same edge weights. Accordingly, the modularity of a partition p of a network G into communities c is given by Newman (2004):

$$Q(p) = \frac{1}{2m} \sum\_{i,j \in V} \left( w\_{ij} - \frac{\deg\_{\mathcal{W}}(i) \cdot \deg\_{\mathcal{W}}(j)}{2m} \right) \delta\_{\mathfrak{c}\_i \mathfrak{c}\_j} \tag{6}$$

where δcicj is 1, if the community ci of node i is the same as the community cj of node j, and 0 otherwise, and m = 1 2 P i,j∈V wij is the total sum of edge weights in a network. Although there are many different definitions in literature about what a community consists of, we define a community as a group of strongly interconnected nodes that make only weak connections to other communities. In addition, a node can maximally contribute to one community. Hence we want to find the partition that maximizes modularity, which is computationally very demanding, so it is important to use a very effective algorithm. We therefore use the fast algorithm by Blondel et al. (2008), which is implemented in the Python packages community. Unfortunately this implementation is only suitable for undirected graphs, so we investigate modularity only for these type of graphs.

#### 5.6. Statistical Model

The generated graph data is used as input for an exchangeable parametric statistical model. Let us recall that the purpose of the fMRI scan of a patient is to give the clinician a likelihood for the patient's health condition,

P(graph data from fMRI scan | health condition ∧ prior info), (7)

which she combines with the likelihoods from other tests and her initial probability assignment, to obtain via Bayes's theorem a final probability for the health condition (Sox et al., 2013):

final probability z }| { P(health condition| results of all tests ∧ prior info) ∝ likelihoods P(graph data from fMRI scan| health condition ∧ prior info) ×P(results of other tests| health condition ∧ prior info) × · · ·

$$\times \underbrace{\text{P(heath condition|prior info)}}\_{\text{favorable probability}}.\tag{8}$$

The prior information also includes test results from previous patients, so that the prediction becomes more accurate and reliable, the more patients have been previously observed.

The functional dependence of the likelihood on the graph data is determined by the statistical model we use, and may be different for each health condition. The statistical model is determined by additional assumptions or hypotheses. Such hypotheses and the functional form of the likelihood may depend on the particular space of graph data (e.g., real-valued, or positive, or bounded within a finite range, or combinations thereof), and therefore on the graph construction method.

As explained in section 2.3, our purpose is to assess as far as possible the relative predictive power of the different graph construction methods. We therefore would like the functional dependence on the graph data space to be minimal. In the present study we adopt the working hypothesis that only the first and second empirical moments—means and correlations—of the graph data from past patients with the same health condition are relevant to make predictions about a new patient. This hypothesis is adopted for all graph construction methods. We also assume our initial knowledge of the graph data to be approximately invariant under rescalings of their values (Minka, 2001). Finally, we do not take into account the natural range of variability (positive, bounded, etc.) of the graph data; this choice does not seem to impact the predictive power of the model (Porta Mana et al., 2018).

These assumptions almost uniquely determine the statistical model and the likelihood (Porta Mana et al., 2018): it turns out to be a multivariate t distribution (Minka, 2001; Kotz and Nadarajah, 2004; Murphy, 2007). More precisely: select a particular health condition, e.g., Alzheimer's disease. Denote with **f** 0 the d-dimensional vector of graph data obtained from the patient's fMRI scan via a particular graph construction method, and with (**f** i ) the graph data of n previous patients with the selected health condition. Then the likelihood that the present patient has the selected health condition is

$$\begin{aligned} \mathrm{p}[f\_0 \mid (f\_i), \kappa\_0, \delta\_0, \upsilon\_0, \Delta\_0, M] &= \mathrm{p}(f\_0 \mid \kappa, \delta, \upsilon, \Delta, M) \\ &= \mathrm{t} \Big[ f\_0 \mid \upsilon - d + 1, \delta, \frac{\kappa + 1}{\kappa \left( \upsilon - d + 1 \right)} \Delta \Big] \end{aligned} \tag{9}$$

with (10)

$$\begin{aligned} \boldsymbol{\delta} &= \frac{\kappa\_0 \,\delta\_0 + n \overline{f}}{\kappa\_0 + n}, \quad \mathbf{A} = \mathbf{A}\_0 + n \,\mathrm{Cov}(\mathbf{f}) \\ &+ \frac{\kappa\_0 \, n}{\kappa\_0 + n} (\overline{f} - \delta\_0)(\overline{f} - \delta\_0)^\mathsf{T}, \end{aligned}$$

κ = κ0 + n, ν = ν0 + n,

where t is a multivariate t distribution with ν − d + 1 degrees of freedom, mean δ, and scale matrix κ+1 κ (ν−d+1)∆, and

$$\overline{f} := \frac{1}{n} \sum\_{i} f\_{i\text{'}} \qquad \text{Cov(f)} := \frac{1}{n} \sum\_{i} (f\_i - \overline{f})(f\_i - \overline{f})^\mathsf{T} \tag{11}$$

are the empirical mean and covariance matrix of the previous graph data.

The parameters κ0, ν0, δ0, ∆0 should reflect our initial knowledge of the graph parameters. For the reasons explained above and in section 2.3, we fix one set of values identically for all graph construction methods: κ = 1, (δ0)a = 0.5, ∆0 = 2.5**I**, where **I** is the identity matrix. These values yield an initial distribution (before any data from previous patients) centered on positive values of unit order of magnitude, as all the graph data indeed are for each graph construction method.

#### 5.7. Supportive Evaluation Measures of Graph Construction Methods 5.7.1. Significance Test

We measure the significance level of the mean values of a graph property distribution between pairs of the three healthy conditions (control-AD, control-MCI, MCI-AD) based on the Student's t-test, if variances are equal (F-test ), and Welch's t-test otherwise. The underlying null hypothesis is that the means of the two data arrays are assumed to be equal, which is rejected for p-values smaller then 0.05.

#### 5.7.2. Dendrograms of Subject Order

Subjects indexed from 1 to 56 (total number of participants) across all health conditions are ordered according to the mean values of a given graph property distribution. The indices of the ordering (the rank) calculated for each graph construction method is then used in order to construct the dendrogram. In the dendrogram, the Euclidean distance between two indices arrays is indicated by the height of the top of the U-link linking the two arrays. In addition, arrays with a small distance are clustered together.

#### 5.7.3. Support Vector Machines

For all complete graphs constructed by all different graph construction methods, we apply a support vector classification (Python: sklearn.svm.SVC) on each pair of health conditions (control-AD, control-MCI, MCI-AD). Hereby we choose the graph properties such that the performance of the algorithm maximizes. We use the default parameters and do not optimize performance by varying the kernel coefficient or the penalty parameter of the error term.

#### ETHICS STATEMENT

This study was part of a larger study, which was approved by the local ethics committee, in accordance with the declaration of Helsinki and performed after informed written consent of each participant. Healthy participants were reimbursed. AD patients were not reimbursed since imaging was part of their diagnostic procedures. We did, however, pay for and organize their traveling costs and lunch.

#### AUTHOR CONTRIBUTIONS

CB constructed the graphs and calculated and analyzed the graph properties. She also applied the statistical analysis, formulated together with PP, to the data. KD, HJ, NR, BvR, JD, OO, K-JL, GF and JK contributed to the conception of the study design and recruited patients. KD, NR, BvR, and JD organized and performed fMRI scanning. KD and HJ applied primary preprocessing to the fMRI data. The manuscript was written by CB, AM, and PP, with additional editing by HJ and JK.

#### REFERENCES

#### FUNDING

We acknowledge partial support by the Helmholtz Alliance through the Initiative and Networking Fund of the Helmholtz Association and the Helmholtz Portfolio theme Supercomputing and Modeling for the Human Brain and the German Research Foundation (DFG; grant DI 1721/3-1 [KFO219-TP9]). This work was also supported by a DFG individual grant JA 2336/1-1 (HJ) and by a grant of the Marga and Walter Boll Foundation, Kerpen, Germany, to GF and JK.

#### ACKNOWLEDGMENTS

PP thanks Mari & Miri for continuous encouragement, affection, and support; the kind staff at Iris; and Buster Keaton and Saitama for filling life with awe and inspiration. We are grateful to Simone Buttler for her important help with regard to the calculation of graph properties, and Fahad Khalid and Andreas Müller of the SimLab Neuroscience at the Jülich Supercomputing Center for their expertise in graph visualization. We also acknowledge the support and expert advice by Alper Yegenoglu, Paulina Dabrowska, and Dr. Jyotika Bahuguna. We would like to thank Dr. Gabriele Stoffels, Dr. Christian Filss, and Nathalie Judov for their assistance and generous support. We also acknowledge the technical support and advice of Prof. Dr. Hans Herzog, Dr. Elena Rota Kops, Lutz Tellmann, and Dr. Daniel Pflugfelder. Finally, we are grateful to Kornelia Frey, Suzanne Schaden, and Silke Frensch for their important help in data acquisition. Thanks are extended to Prof. Dr. Nadim Jon Shah for support with the MRI.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins. 2018.00528/full#supplementary-material

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Bachmann, Jacobs, Porta Mana, Dillen, Richter, von Reutern, Dronse, Onur, Langen, Fink, Kukolja and Morrison. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Analysis of Alzheimer's Disease Based on the Random Neural Network Cluster in fMRI

Xia-an Bi\*, Qin Jiang, Qi Sun, Qing Shu and Yingchao Liu

College of Information Science and Engineering, Hunan Normal University, Changsha, China

As Alzheimer's disease (AD) is featured with degeneration and irreversibility, the diagnosis of AD at early stage is important. In recent years, some researchers have tried to apply neural network (NN) to classify AD patients from healthy controls (HC) based on functional MRI (fMRI) data. But most study focus on a single NN and the classification accuracy was not high. Therefore, this paper used the random neural network cluster which was composed of multiple NNs to improve classification performance. Sixty one subjects (25 AD and 36 HC) were acquired from the Alzheimer's Disease Neuroimaging Initiative (ADNI) dataset. This method not only could be used in the classification, but also could be used for feature selection. Firstly, we chose Elman NN from five types of NNs as the optimal base classifier of random neural network cluster based on the results of feature selection, and the accuracies of the random Elman neural network cluster could reach to 92.31% which was the highest and stable. Then we used the random Elman neural network cluster to select significant features and these features could be used to find out the abnormal regions. Finally, we found out 23 abnormal regions such as the precentral gyrus, the frontal gyrus and supplementary motor area. These results fully show that the random neural network cluster is worthwhile and meaningful for the diagnosis of AD.

#### Edited by:

Mohammad Amjad Kamal, King Abdulaziz University, Saudi Arabia

#### Reviewed by:

Rifai Chai, University of Technology Sydney, Australia Na Li, Central South University, China

> \*Correspondence: Xia-an Bi<EMAIL_ADDRESS>
Received: 03 February 2018 Accepted: 22 August 2018 Published: 07 September 2018

#### Citation:

Bi X-a, Jiang Q, Sun Q, Shu Q and Liu Y (2018) Analysis of Alzheimer's Disease Based on the Random Neural Network Cluster in fMRI. Front. Neuroinform. 12:60. doi: 10.3389/fninf.2018.00060 Keywords: random neural network cluster, fMRI, classification, Alzheimer's disease, functional connectivity

### INTRODUCTION

Alzheimer's disease (AD) is degenerative and irreversible which results from the cognitive decline. The progress of AD increases with age and the disease easily triggers the other psychiatric diseases, and eventually causing dementia. In 2006, the number of AD patients is 26.6 million all around the world and the number would be quadruple by 2050. Therefore, it is meaningful for clinician to track its progression and diagnose the disease. Currently, there have been many different neuroimaging techniques which can be applied to diagnose AD, such as ElectroEncephaloGram (EEG) (Engels et al., 2015), Single Photon Emission Computed Tomography (SPECT) (Prosser et al., 2015), Positron Emission Tomography (PET) (Pagani et al., 2017), Magnetoencephalographic (MEG) (Engels et al., 2016), and functional magnetic resonance imaging (fMRI) (Griffanti et al., 2015). Among these techniques, fMRI is widely used in the diagnosis of AD, because it adopts a non-invasive way and could be used to find the differences of brain regions between AD and healthy controls (HC) (Challis et al., 2015).

Machine learning is a method of pattern recognition and has been used for the study of AD in recent years (Moradi et al., 2015; Khazaee et al., 2016). Among many machine learning methods,

artificial neural network (ANN) is a useful classification method which evolves from human brain (Er et al., 2016). Several previous studies showed that neural network (NN) was able to be applied to diagnose neurological disease. Suk et al. (2016) proposed a method that combined deep learning and state-space model to classify Mild Cognitive Impairment (MCI) patients from HC, and the accuracy was 72.58%. Gao et al. (2017) employed an advanced convolution neural network (CNN) with 2D and 3D to diagnose AD, and the average accuracy reached to 87.6%. Ortiz et al. (2016) used deep learning architectures to classify AD patients from HC, and the accuracy approximately reached to 90%. Ortiz et al. (2013) used Learning Vector Quantization (LVQ) algorithm to classify AD patients from HC, and the accuracy was close to 90%. Luo et al. (2017) applied CNN to classify AD patients from HC, and the sensitivity and specificity of classification was 1 and 0.93 respectively. Suk et al. (2014) used deep learning to classify AD patients from HC, and the accuracy was close to 93.35%.

In existing studies, a single NN is often used to classify patients with neurological diseases and HC, and the accuracy of classification is considerable which indicates that NN is a powerful classification model (Hirschauer et al., 2015; Anthimopoulos et al., 2016). As the features of neuroimaging data are characterized by high dimension, significant information of original variables would be lost in the process of dimensionality reduction in traditional methods such as principal component analysis, local linear embedding and linear discriminate analysis (Mattioni and Jurs, 2003; McKeown et al., 2007; Mannfolk et al., 2010). In this paper, the method of the random neural network cluster is proposed to classify AD from HC. This method not only could be used in the classification, but also could be used for feature selection. The procedure of this method is as follows. Firstly, we chose Elman NN as the optimal base classifier from five types of neural networks [Back Propagation (BP) NN, Elman NN, PNN, Learning Vector Quantization (LVQ) NN and Competitive NN] based on the results of feature selection, and the accuracies of the random Elman neural network cluster could reach to 92.31% which is the highest and stable. Then we used the random Elman neural network cluster to select significant features and these features could be used to find out the abnormal regions. Finally, we found out 23 abnormal regions, such as the precentral gyrus, the frontal gyrus and supplementary motor area. These results fully show that the random neural network cluster is worthwhile and meaningful for the diagnosis of AD.

#### MATERIALS AND METHODS

#### Subjects

The experimental data was collected from the Alzheimer's Disease Neuroimaging Initiative (ADNI)1 dataset which includes a variety of neuroimaging data. The ADNI study was approved by Institutional Review Board (IRB) of each participating site. All ADNI subjects together with their legal representatives should have written informed consent before collecting clinical, genetic

1http://adni.loni.usc.edu/

and imaging data. The following two criteria need to meet when selecting data. One criterion is that the selected data should be resting-state fMRI data. Another criterion is that the selected data should have mini-mental state examination (MMSE) and clinical dementia rating (CDR) scores, and this criterion ensures that the selected data is homologous. Finally, 61 subjects were selected which consisted of 25 AD patients and 36 HC.

#### Data Collection and Preprocessing

Scanning images were acquired on a Philips Medical Systems 3 Tesla MRI scanner. Acquisition parameters included: pulse sequence = GR, TR = 3,000 ms, TE = 30 ms, matrix = 64∗ 64, slice thickness = 3.3 mm, slice number = 48, flip angle = 80◦ .

To decrease the influence of signal-to-noise ratio of the fMRI images, the selected data should be preprocessed. The data was preprocessed based on the Data Processing Assistant for Resting-State fMRI (DPARSF) software (Chao-Gan and Yu-Feng, 2010; Wang et al., 2013). The process of the data preprocessing included: converting DICOM format into NIFTI format; removing first 10 time points; slicing timing (Sarraf et al., 2016); realigning (Jenkinson et al., 2002); normalizating images into the echo planar imaging (EPI) template (Misaki et al., 2010); temporal smoothing; removing the effect of low-level (<0.01 HZ) and high-level (>0.08 HZ) noise by a high-pass temporal filtering (Challis et al., 2015); removing covariates such as the whole brain signal and cerebrospinal fluid signal.

#### Functional Connectivity of Brain

After the preprocessing steps, the analysis of the functional connectivity of brain was followed. In this paper, we chose functional connectivity as the sample feature. The extraction of the functional connectivity is as follows.  Secondly, the time series of each region were extracted. Thirdly, the Pearson correlation coefficient between two separated brain regions could be defined as the functional connectivity (Friston et al., 1993). Finally, 4,005 (90 × 89/2) functional connectivity of each subject were taken as their features.

#### The Random Neural Network Cluster

The 4,005 functional connectivity belongs to the highdimensional feature which causes the problems of computation difficulty and dimensions of disaster. Moreover, the highdimensional features are likely to result in numerous redundant and irrelevant features which may decrease the classification performance. Therefore, the irrelevant features are needed to be removed by feature selection (Azar and Hassanien, 2015).

There are many methods of feature selection, such as principal component analysis, local linear embedding and linear discriminate analysis. However, in these methods, the process of selecting features may cause the loss of the original information, and the classification performance may be unsatisfactory (Zhou et al., 2015; Jolliffe and Cadima, 2016; Alam and Kwon, 2017).

To solve the problem above, this paper proposes the random neural network cluster by randomly selecting samples and features. The random neural network cluster could be used

to classify AD patients from HC and select features. In addition, the new method could achieve the purpose of reducing dimensionality, avoid losing significant information and improve classification performance.

#### The Design and Classification Accuracy of the Random Neural Network Cluster

The basic thought of the random neural network cluster is ensemble learning whose basic classifier is neural network. The detailed process of establishing the random neural network cluster is as follows. Firstly, the sample dataset D is randomly divided into a training set N1 and a test set N2, where D = N1 + N2. The real label of HC and AD is +1 and −1 respectively. Secondly, we randomly select n samples from the training set and m features from the 4,005 features. Thirdly, the selected samples and features are used to establish a single NN and the process of the second and third step is repeated for k times. Thus, k NNs are obtained which construct the random neural network cluster. It can be seen from this process that the method and the conventional process are essentially not the same in essence. **Figure 1** shows the formation of the random neural network cluster.

When a new sample enters the random neural network cluster, we could obtain k class labels from k NNs and the majority of class label is made as the predicted label of the new sample. Similarly, we could get the label of each sample in the test set N2. Then the predicted label is compared with the real label to judge whether they are consistent, and we assume C as the number of consistent situations. The accuracy of the random neural network cluster equals to C/N2.

### Extracting Features From the Random Neural Network Cluster

As the features are randomly selected, the NNs constructed by these features have different characteristics. In this paper, the accuracy of each NN is used to evaluate the significance of selected features. The features that make significant contributions to the accuracy of the NN are called the "significant features". The process of extracting significant features is as follows.

Firstly, the test samples are used to obtain the accuracy of each NN in the random neural network cluster. If the accuracy of a NN is high, the corresponding features are significant. Next, the features in each NN with high accuracy are extracted to form a feature matrix. Finally, we count the frequency of each feature in the feature matrix, and extract the features with high frequency which are called as "significant features." **Figure 2** shows the process of selecting "significant features." The significant features tremendously contribute to the accuracy of a NN, thus they also make great contributions to the accuracy of the random neural network cluster. In this paper, we use the significant features to find the difference between AD patients and HC.

#### The Abnormal Brain Regions

As mentioned above, the significant features make great contributions to the accuracy of the random neural network cluster, thus we could find the difference between AD and HC

through these significant features. In this paper, the significant features are regarded as the abnormal functional connectivity. As the functional connectivity is defined as the relationship of two brain regions, the extracted significant features could be used to find the abnormal brain regions between AD and HC. In order to estimate the abnormal degree of a brain region, the number of features related to a certain region is considered as the weight. When there is no functional connectivity related to a certain brain region, the weight of the region is 0. When the weight is greater, the abnormal degree of the brain region is higher.

#### Experiment Design

The process of the experiment involves six steps in this paper.

Step 1. Building the random neural network cluster. Firstly, 61 subjects are divided into a training set and a test set according to the proportion of 8:2. Thus, the number of training samples and test samples is 48 and 13 respectively. Then, 45 subjects are randomly selected from 48 training samples and 120 features are randomly selected from 4,005 features to build a NN. Similarly, we build 1,000 NNs to construct the random neural network cluster. The result of the random neural network cluster is calculated by using neural network toolbox.

Step 2. Selecting significant NNs. We select the NNs whose accuracy is >0.6 from the 1,000 NNs and call these NNs as significant NNs.

Step 3. Selecting significant features. The features corresponding to the significant NNs form the feature matrix. Then we count the frequency of each feature, and sort the features with a descending order. Finally, we retain the first 240 features which are regarded as the original significant features.

Step 4. Determining the optimal number of the significant features. We change the number of original significant features from 140 to 240 features with a step of 10 as the new significant features, thus there are 11 types of numbers of new significant features. Then we select 120 features from the new significant features to construct a NN. We select 120 features from the original significant features, and change the number of significant features from 140 to 240 features with a step of 10. Thus, we could obtain 11 random neural network clusters. The number of the significant features corresponding to the random neural network cluster with highest accuracy is the optimal number.

Step 5. We repeat the step1-step 4 by using five types of NNs. They are BP NN, Elman NN, PNN, LVQ NN, and Competitive NN. Then we choose one of them as the best base classifier who has the highest accuracy in step 4.

Step 6. Finding abnormal brain regions through the significant features in step 4 based on the best base classifier.

#### RESULTS

#### The Demographic Information of Participants

In this study, the selected 61 subjects include 25 AD patients and 36 HC. The gender and age difference between the AD group and HC group data are examined by two-sample t-test and chi-square test respectively. The result is shown in **Table 1**. It is referred that the two groups have no statistical significance in the gender and the age.

#### Classification Results

In this study, five different types of NNs are applied to construct five different types of random neural network clusters. The

AD, Alzheimer's disease; HC, Healthy Control.

number of significant features in each random neural network cluster changes from 140 to 240 and the step is 10. Therefore, we can obtain 11 results for each type of random neural network cluster, and their parameters have been made appropriate adjustments to get better results. The accuracies of five types of random neural network clusters are shown in **Figure 3**. It is referred that the accuracies of random Competitive neural network cluster are not stable; the accuracies of random Elman neural network cluster and random Probabilistic neural network cluster are high, and their highest accuracy reaches to 92.31%; the accuracies of random BP neural network cluster and random LVQ neural network cluster are lower than the random Elman neural network cluster. The accuracies of the random neural network cluster changes when the number of features changes. As the highest accuracies of the random Elman neural network cluster are relatively stable and high, this paper chooses Elman neural network as the best base classifier. From **Figure 3** we can learn that when the number of significant features is 180, the accuracy is the highest. Therefore, 180 is the optimal number of significant features.

In order to better compare the performances of the random neural network cluster and a single NN, we display the corresponding 1,000 NNs' accuracies of five types of random neural network clusters in the **Figure 4**. From **Figures 3**, **4** we could learn that the accuracy of random neural network cluster is higher and more stable than a single NN except for the random Competitive neural network cluster. In addition, we also show the training errors, test errors and running time of five types of random neural network clusters in **Table 2**.

#### The Abnormal Brain Regions

The first 180 features constitute the optimal feature set which is used to find the abnormal brain regions between AD patients and HC. **Table 3** shows the abnormal brain regions with higher weight, their abbreviation and volume. **Figure 5** shows the abnormal degree of brain regions by using the Brain-NetViewer2 .

2http://www.nitrc.org/projects/bnv/

FIGURE 3 | The accuracies of five different types of random neural network clusters.

fninf-12-00060 September 6, 2018 Time: 18:23 # 4

The node represents the brain region. The size of a node represents the weight of the brain region and it also indicates the abnormal degree of the brain region.

In this paper, we focus on the brain region whose weight is >17. **Figure 6** shows the functional connectivity between the 23 brain regions in the optimal feature set. **Figure 7** shows the functional connectivity between PreCG and other brain regions. The node in **Figures 6**, **7** also represents the brain region and the size of a node represents the weight of the brain region. Besides, the line represents the functional connectivity between two brain regions.

### DISCUSSION

#### Classification Performance

A number of researchers have tried to classify and diagnose AD patients from HC in the past few years. For instance, Khazaee et al. (2016) applied machine learning methods to classify MCI patients, AD patients and HC, and the accuracy of AD and HC is 72%. Kim et al. (2016) applied DNN to classify schizophrenia (SZ) patients and HC based on fMRI image, and the accuracy is 85.8%. Zhang et al. (2014) proposed a kernel support vector machine decision tree (kSVM-DT) to classify MCI patients, AD patients and HC based on the MRI data, and the accuracy of AD and HC is 96%. However, there were several disadvantages in these studies. For instance, as the features of neuroimaging data are characterized by high dimension, significant information of

TABLE 2 | The errors and running time of five types of random neural network clusters.

original variables would be lost in the process of dimensionality reduction in traditional methods.

In this paper, the new method was proposed to avoid the loss of information and improve classification performance. We used five different types of NNs as the base classifier to build five different types of random neural network clusters. In the five types of random neural network clusters, the classification accuracies of the random Elman neural network cluster and the random PNN cluster could reach to 92.3%. As the resting-state fMRI data dynamically changes in a period of time and the Elman NN is able to deal with the time and spatial domain data, the Elman NN is more suitable for the fMRI data. The accuracies of the random Probabilistic neural network cluster could also reach to 92.3%, but the running time of the PNN is long which makes it not suitable for the base classifier of the random neural network cluster. Thus, we finally choose the random Elman neural network cluster as the best base classifier. To solve the problems caused by the high-dimensional data, the traditional methods suffer from the loss of some information in the process of reducing dimensionality. In this paper, we proposed a random neural network cluster which could be a good solution in dealing with the calculation of large samples and high-dimensional data. Moreover, the highest accuracies of the random Elman neural network cluster could reach to 92.3%. Therefore, this method could process high-dimensional data without information loss and improve classification performance.

The basic classifier used in this paper is the neural network. This classifier is indeed an existing classifier, and it is not the innovation point of this paper. The innovation of this article is the integrated cluster, which is the innovation in structure.

This paper researches on fMRI image, there is no innovation in the classification indicators of this article. But the classification indicators have been used in research. We used the random neural network cluster to classify the subjects and feature selection, and we got good results. That is the innovation in application.

#### The Additional Details of the Random Neural Network Cluster

In the part of method, we have introduced many elementary details in our method and in this part we would introduce some additional details.

In order to make the random neural network cluster less complex, each NN occupies the same weight in the random neural network cluster. When determining the label of a new sample, each NN in the random neural network cluster predicts the label of the sample. The majority of class label is made as the predicted label of the sample which is equivalent to that the weight of each NN is the same. From **Figures 3**, **4** we could learn that the accuracy of a single NN is lower than the accuracy of the random neural network cluster. Moreover, when the number of significant features changes, the accuracies of five types of random neural network clusters are high and stable except for the random Competitive neural network cluster. This fully demonstrates that the robustness of the random neural network cluster is good.

#### TABLE 3 | Abnormal brain regions.

fninf-12-00060 September 6, 2018 Time: 18:23 # 6

In terms of the complexity of random neural network cluster, it is mainly reflected in the following two aspects. One aspect is that the number of base classifiers is large which makes the process of constructing a random neural network cluster complex. Another aspect is that we need to construct multiple random neural network clusters during finding the optimal feature set which also makes the method complex.

There are two types of parameters in our method. One type is the parameters of the random neural network cluster which are decided by the accuracy of the random neural network cluster through hundreds of experiments. Another type is the parameters

of the NN which are decided by the neural network toolbox that automatically selects the optimal parameters.

The random neural network cluster is constructed by the training set and the performance is tested by the test set. The experimental results show that the random neural network cluster not only performs well on the training set, but also performs well on the test set. This fully demonstrates that

the overfitting does not exist in the random neural network cluster.

 They are the PreCG, the OLF, the ORBsup, the IFGtriang, SMA, the SFGdor, the ORBmid, the ORBinf, the SFGmed.L, the IFGoperc, the ROL, and the MFG.

Previous studies have concluded that these abnormal brain regions are associated with AD patients. Khazaee et al. (2015) found that the lingual gyrus, the occipital gyrus and the superior frontal gyrus are abnormal brain regions in AD patients. Agosta et al. (2012) pointed out that the functional connectivity changed in the default mode network (DMN) and frontal networks in AD patients. Wang et al. (2007) discovered the decreased abnormal functional connectivity in the prefrontal and parietal lobes, meanwhile the increased abnormal functional connectivity in the occipital lobe in AD patients. Greicius et al. (2004) found that AD patients showed less deactivation in the anterior frontal, precuneus and posterior cingulate cortex. He et al. (2007) concluded that regional coherence of AD patients significantly decreased in the posterior cingulate cortex/precuneus. Binnewijzend et al. (2012) declared that the abnormal brain regions of AD patients located in the precuneus and posterior cingulate cortex. Golby et al. (2005) concluded that impaired activation changed in the temporal lobe and fusiform regions in AD patients. Grossman et al. (2003) found the activation of AD patients in the left poster lateral temporal and inferior parietal cortex. It is referred that our results are consistent with previous studies. In this paper, we focus on the precentral gyrus and the frontal gyrus which have larger weights.

## Precentral Gyrus (PreCG)

The PreCG has the greatest weight in the abnormal brain regions. It is referred that the PreCG makes a great contribution to classify AD and HC in the random Elman neural network cluster. The PreCG locates in the primary motor area (Hopkins et al., 2017), and the superior part of PreCG is responsible for motor hand function (Yousry et al., 1997; Nebel et al., 2014).

Existing studies have found that the PreCG is abnormal in AD patients. Kang et al. (2013) showed that AD patients performed significant cortical thinning in superior and medial frontal gyrus, left precentral gyrus, postcentral gyrus, paracentral lobule, precuneus and superior parietal lobule. Bi et al. (2018) proposed a method of random support vector machine cluster to diagnose AD, and found out several disorder regions including inferior frontal gyrus, superior frontal gyrus, precentral gyrus and cingulate cortex. Zhang et al. (2015) pointed out that the middle occipital gyrus, the postcentral gyrus, the PreCG and precuneus were important in diagnosing AD. Brickman et al. (2015) discovered that the PreCG and right middle frontal gyrus were abnormal in AD patients based on T2-weighted MRI. Dani et al. (2017) used white matter hypointensity (WMH) volume to diagnose AD based on MRI and PET data, and they found out the difference between the PreCG, and the right medial and anterior part of orbital gyrus.

The abnormal functional connectivity between PreCG and other brain regions may lead to the physical movement dysfunction in AD patients. The above results revealed that PreCG may be a clinical diagnosis of AD in future.

#### Frontal Gyrus (FG)

The FG had a relatively higher weight in the abnormal regions. It is referred that the FG makes a great contribution to classify AD and HC in the random Elman neural network cluster. The left of inferior frontal gyrus is associated with language (Costafreda et al., 2006). The superior frontal gyrus contributes to higher cognitive functions and particularly the learning and working memory (WM) (Boisgueheneuc et al., 2006; Carter et al., 2006). Eliasova et al. (2014) explored that the improvement of attention and psychomotor speed resulted from the abnormity of the right IFG in AD patients.

Existing studies have found that the FG is abnormal in AD patients. Schultz et al. (2015) found out several abnormal brain regions in AD patients including the hippocampus, the posterior cingulate, the anterior cingulate and the middle frontal gyrus. Griffanti et al. (2016) pointed out that there was difference in orbito frontal network between AD patient and HC. Yetkin et al. (2006) evaluated brain activation in AD patients and HC based on fMRI while performing a WM task and found that AD group showed more activation in the right superior frontal gyrus, bilateral middle temporal, middle frontal, anterior cingulate and fusiform gyri. Karas et al. (2007) used voxel-based morphometry (VBM) to examine AD patients, and they found out several disorder regions including the left superior and inferior temporal gyrus, and the left superior frontal gyrus. Our experiment results are consistent with the existing studies.

The abnormal functional connectivity between FG and other brain regions may lead to the memory dysfunction in

AD patients. This abnormal brain region has significant effect on the identification of potentially effective biomarkers for the diagnosis of AD.

In this paper, the random neural network cluster is proposed to classify AD patients and HC and the abnormal brain regions are found out on the basis of fMRI data. Moreover, these findings suggest that the random neural network cluster might be an appropriate approach for diagnosing AD. This new method has some advantages. Firstly, the NNs make different contributions to the random neural network cluster which could make full use of each NN's characteristics, thus it could avoid losing information. Secondly, it is able to effectively deal with a large dataset even when there is many missing data. Finally, it is also able to select significant features from high-dimensional features. The new method presents a good classification performance with accuracy of 92.3% and detects several abnormal brain regions which would have influence in diagnosing AD. However, it has several limitations. Firstly, this paper only used the fMRI data and the future studies could integrate other imaging data to obtain comprehensive brain activity. Secondly, this paper only studied the brain activity, and the future studies could combine the brain and cerebellum activity. Thirdly, this paper only studies the brain difference between AD patients and HC, and the future studies could focus on the brain relationship in AD patients. Finally, the functional connectivity was regarded as feature in this paper, and the future studies could choose other attributions of the brain as feature such as clustering coefficient, degree and shortest path.

#### REFERENCES

#### ETHICS STATEMENT

This study was carried out in accordance with the recommendations of National Institute of Aging-Alzheimer's Association (NIA-AA) workgroup guidelines, Institutional Review Board (IRB). The study was approved by Institutional Review Board (IRB) of each participating site, including the Banner Alzheimer's Institute, and was conducted in accordance with Federal Regulations, the Internal Conference on Harmonization (ICH), and Good Clinical Practices (GCP).

#### AUTHOR CONTRIBUTIONS

X-aB proposed the design of the work and revised it critically for important intellectual content. QSu and QJ carried out the experiment for the work and drafted part of the work. YL and QSh collected, interpreted the data, and drafted part of the work. All the authors approved the final version to be published and agreed to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.

#### FUNDING

This work was supported by the National Science Foundation of China (Nos. 61502167 and 61473259).

classification. IEEE Trans. Geosci. Remote Sens. 53, 1082–1095. doi: 10.1109/ TGRS.2014.2333539

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Bi, Jiang, Sun, Shu and Liu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Decreased Bilateral FDG-PET Uptake and Inter-Hemispheric Connectivity in Multi-Domain Amnestic Mild Cognitive Impairment Patients: A Preliminary Study

Xiao Luo1†‡, Kaicheng Li1† , Qingze Zeng1 , Peiyu Huang1 , Yeerfan Jiaerken1 , Tiantian Qiu1 , Xiaojun Xu1 , Jiong Zhou2 , Jingjing Xu1 and Minming Zhang1 \* for the Alzheimer's Disease Neuroimaging Initiative (ADNI)§

1 Department of Radiology, The Second Affiliated Hospital of Zhejiang University School of Medicine, Hangzhou, China, 2 Department of Neurology, The Second Affiliated Hospital of Zhejiang University School of Medicine, Hangzhou, China

Background: Amnestic mild cognitive impairment (aMCI) is a heterogeneous condition. Based on clinical symptoms, aMCI could be categorized into single-domain aMCI (SDaMCI, only memory deficit) and multi-domain aMCI (MD-aMCI, one or more cognitive domain deficit). As core intrinsic functional architecture, inter-hemispheric connectivity maintains many cognitive abilities. However, few studies investigated whether SD-aMCI and MD-aMCI have different inter-hemispheric connectivity pattern.

Methods: We evaluated inter-hemispheric connection pattern using fluorine-18 positron emission tomography – fluorodeoxyglucose (18F PET-FDG), resting-state functional MRI and structural T1 in 49 controls, 32 SD-aMCI, and 32 MD-aMCI patients. Specifically, we analyzed the 18F PET-FDG (intensity normalized by cerebellar vermis) in a voxel-wise manner. Then, we estimated inter-hemispheric functional and structural connectivity by calculating the voxel-mirrored homotopic connectivity (VMHC) and corpus callosum (CC) subregions volume. Further, we correlated inter-hemispheric indices with the behavioral score and pathological biomarkers.

Results: We found that MD-aMCI exhibited more several inter-hemispheric connectivity damages than SD-aMCI. Specifically, MD-aMCI displayed hypometabolism in the bilateral middle temporal gyrus (MTG), inferior parietal lobe, and left precuneus (PCu) (p < 0.001, corrected). Correspondingly, MD-aMCI showed decreased VMHC in MTG, PCu, calcarine gyrus, and postcentral gyrus, as well as smaller mid-posterior CC than the SD-aMCI and controls (p < 0.05, corrected). Contrary to MD-aMCI, there were no neuroimaging indices with significant differences between SD-aMCI and controls, except reduced hypometabolism in bilateral MTG. Within aMCI patients, hypometabolism and reduced inter-hemispheric connectivity correlated with worse executive ability. Moreover, hypometabolism indices correlated to increased amyloid deposition.

#### Edited by:

Mohammad Amjad Kamal, King Fahad Medical Research Center, King Abdulaziz University, Saudi Arabia

#### Reviewed by:

Arun Bokde, Trinity College Dublin, Ireland Panteleimon Giannakopoulos, Université de Genève, Switzerland

\*Correspondence:

Minming Zhang<EMAIL_ADDRESS>
†These authors have contributed equally to this work.

‡orcid.org/0000-0003-1743-7842

§The data used in the preparation of this article were obtained from the Alzheimer's Disease Neuroimaging Initiative (ADNI) database (adni.loni.usc.edu). As such, the investigators within the ADNI contributed to the design and implementation of ADNI and provided data but did not participate in the analysis or writing of this report.

> Received: 22 February 2018 Accepted: 14 May 2018 Published: 05 June 2018

#### Citation:

Luo X, Li K, Zeng Q, Huang P, Jiaerken Y, Qiu T, Xu X, Zhou J, Xu J and Zhang M for the Alzheimer's Disease Neuroimaging Initiative (ADNI) (2018) Decreased Bilateral FDG-PET Uptake and Inter-Hemispheric Connectivity in Multi-Domain Amnestic Mild Cognitive Impairment Patients: A Preliminary Study. Front. Aging Neurosci. 10:161. doi: 10.3389/fnagi.2018.00161

**64**

Conclusion: In conclusion, patients with MD-aMCI exhibited the more severe deficit in inter-hemispheric communication than SD-aMCI. This long-range connectivity deficit may contribute to cognitive profiles and potentially serve as a biomarker to estimate disease progression of aMCI patients.

Keywords: mild cognitive impairment, cerebral metabolism, corpus callosum, resting-state functional MRI, voxelmirrored homotopic connectivity

### INTRODUCTION

fnagi-10-00161 June 5, 2018 Time: 12:2 # 2

The amnestic mild cognitive impairment (aMCI) represents an intermediate stage between normal aging and Alzheimer's disease (AD) (Petersen et al., 2014). According to impaired cognitive domains number, aMCI patients could be further categorized into two subtypes: single-domain aMCI (SD-aMCI), characterized by relatively selective memory impairment and the multi-domain aMCI (MD-aMCI), indicating extensive deficits involving at least one other domain (Winblad et al., 2004; Busse et al., 2006). Previous epidemiology study shows that MD-aMCI has a higher risk of clinical progression than the SD-aMCI (Golob et al., 2007). Therefore, to investigate the mechanism underlying these aMCI subtypes may efficiently facilitate clinical early intervention and management. Recently, neuroimaging studies pointed out that the MD-aMCI displays more diffuse gray matter atrophy (Haller et al., 2010; Zhang et al., 2012) and lower brain activity than the SD-aMCI, mainly involving default mode network (DMN) and frontoparietal regions (Li et al., 2014). Despite these studies shedding light into the aMCI pathological mechanism to some extent, it remains unclear whether SD-aMCI and MD-aMCI have different inter-hemispheric connection pattern.

Compared to the other neuroimaging indices, insufficient attention paid to the direct inter-hemispheric connectivity of aMCI patients. Anatomically, the connectivity between hemispheres is a core mode of the brain intrinsic functional architecture. Moreover, this bi-hemispheric communication procedure substantially affects many cognitive domains, including executive and memory functions (Salvador et al., 2008; Stark et al., 2008; Saar-Ashkenazy et al., 2016). Until now, only two functional MRI studies directly explored the interhemispheric functional connectivity in aMCI patients. However, two studies drew the entirely different conclusion. Specifically, Wang et al. (2015) demonstrated that MCI patients exhibited an enhanced inter-hemispheric functional connectivity in the sensorimotor cortex to resist the cognitive decline; however, another work reported the negative result between aMCI and controls (Qiu et al., 2016). Given that aMCI is a heterogeneous condition, we thus hypothesized that this inconsistency might attribute to the different damage pattern of inter-hemispheric connectivity between SD-aMCI and MD-aMCI patients.

Notably, some previous findings focusing on corpus callosum (CC) supported this hypothesis. Anatomically, CC acts as the most robust commissural white matter bundle to maintain the functional connectivity between the hemispheres (Roland et al., 2017). Specifically, some studies demonstrated that the AD and MCI patients exhibit the CC shape change, atrophy, and impaired diffusivity indices impairment, especially in the posterior part (Janowsky et al., 1996; Di Paola et al., 2010; Ardekani et al., 2014; Rasero et al., 2017). Subsequently, some studies further demonstrated patients with MD-aMCI, but not SD-aMCI, have reduced mean diffusivity (MD) in the whole CC compared to controls; moreover, decrease of MD in the CC body is associated with decreased general cognition and executive ability (Li et al., 2013). On the other hand, as a valid diagnostic biomarker in dementia studies, fluorine-18 fluorodeoxyglucose ( 18F FDG) positron emission tomography (PET) research also shows more diffuse hypometabolism in MD-aMCI than the SDaMCI. Interestingly, these hypometabolism regions are mostly located in the bilateral homotopic precuneus (PCu) and the temporoparietal cortex (Caffarra et al., 2008; Cerami et al., 2015). Considering the evidence of CC degeneration and bilateral homotopic hypometabolism, we inferred that different interhemispheric connectivity damage pattern might exist between aMCI subgroups.

To test this hypothesis, we analyzed the structural T1 image, resting-state functional MRI (rsfMRI), 18F PET-FDG, neuropsychological scales, and pathological biomarkers in a relatively large aMCI sample from ADNI database. Firstly, we analyzed inter-hemispheric homotopic functional connectivity by using a voxel-wise method, namely voxel-mirrored homotopic connectivity (VMHC) (Zuo et al., 2010; Luo et al., 2018b). The potential inter-hemispheric structural connectivity was evaluated by estimating CC volume. We divided CC into five part due to its different anatomical connection [i.e., genus connects the frontal areas, while body and splenium part connects the temporal and parietal areas (Fischl, 2012)]. Moreover, to test whether regions with inter-hemispheric disruption accompanied by metabolic abnormalities, we also analyzed 18F FDG-PET data in a voxelwise manner. Additionally, we also correlated neuroimaging indices with cognition and pathological biomarkers (reflecting by CSF and amyloid imaging). Our study aims to (i) compare the alteration of the inter-hemispheric connectivity and the metabolism between aMCI subtypes; (ii) explore the possible interactions between different neuroimaging modalities; and (iii) explore the relationships between neuroimaging indices and cognition.

#### MATERIALS AND METHODS

#### Alzheimer's Disease Neuroimaging Initiative

Data used in this study were obtained from the Alzheimer's disease Neuroimaging Initiative (ADNI) database1 . The ADNI was launched in 2003 by the National Institute on Aging (NIA),

1 adni.loni.usc.edu

the National Institute of Biomedical Imaging and Bioengineering (NIBIB), the Food and Drug Administration (FDA), private pharmaceutical companies, and non-profit organizations, as a \$60 million, 5-year public–private partnership. The primary goal of ADNI has been to test whether serial magnetic resonance imaging (MRI), PET, other biological markers, and clinical and neuropsychological assessment can be combined to measure the progression of mild cognitive impairment (MCI) and early Alzheimer's disease (AD). Determination of sensitive and specific markers of very early AD progression is intended to aid researchers and clinicians in developing new treatments and monitor their effectiveness, as well as lessen the time and cost of clinical trials.

#### Subjects

All procedures performed in studies involving human participants were following the ethical standards of the institutional and national research committee and with the 1975 Helsinki Declaration and its later amendments or comparable ethical standards. The ADNI project was approved by the Institutional Review Boards of all participating institutions, and all participants at each site signed informed consent. Based on ADNI GO and ADNI 2 database, we identified 49 healthy righthanded subjects (normal controls, NC) and 64 aMCI patients, who had undergone structural scans, rsfMRI scans, FDG-PET scans, and neuropsychological evaluation. We downloaded the study data from the ADNI publicly available database before April 15, 2017. According to ADNI protocol, to be classified as healthy controls, the Mini-Mental State Examination (MMSE) for the subject should be between 24 and 30 (inclusive), and clinical dementia rating (CDR) score should be 0. Besides, the subject has no signs of depression (Geriatric Depression Scale, GDS <6) or possible dementia. To be defined as aMCI, the subject had an MMSE score between 24 and 30 (inclusive), memory complaint, objective abnormal memory evidence, and a CDR score of 0.5. Besides, aMCI patients' general cognition preserved well (cannot meet AD diagnostic criteria). Also, there were no signs of depression (GDS score <6) in aMCI patients.

Followed by exclusion standards: the subject who has a history of apparent head trauma, other neurological or major psychiatric disorder, and alcohol/drug abuse. Additionally, subjects were also excluded if they exhibited a significant vascular disease risk history, defined as Hachinski Ischemia Scale (HIS) scores higher than 4.

### Neuropsychological Evaluation and MCI Subtypes Diagnosis

All subjects underwent comprehensive neuropsychological tests. Content includes general mental status (MMSE score, Alzheimer's Disease Assessment Scale, ADAS) and other cognitive domain, including memory function (Auditory Verbal Learning Test, AVLT; Wechsler Memory Scale-Logical Memory, WMS-LM, including immediate and delayed score), processing speed (Trail-Making Test, Part A, TMT-A), visuospatial function (Clock-Drawing Test, CDT), executive function (Trail-Making Test, Part B, TMT-B), and language ability (Boston Naming Test, BNT). According to the previous studies, we used the composite scores for executive functioning and memory (Crane et al., 2012; Gibbons et al., 2012).

The aMCI patients were divided into two subtypes by the performance in cognitive domains (Shu et al., 2012). Specifically, SD-aMCI refers to aMCI patients having an impairment in memory alone; MD-aMCI refers to aMCI patients having an impairment in memory and other cognitive domain (at least one). We defined impairment as the presence of a test scoring 1.5 standard deviations (SD) below the average score, from age- and education-matched healthy controls from ADNI. Specifically, the total number of healthy controls from ADNI is 198, mean age is 73.03 ± 6.20, and mean education attainment is 16.47 ± 2.45. No significant difference existed in terms of age and education level (p > 0.05) between ADNI controls and aMCI subjects (mean age: 73.33 ± 5.69; education level: 15.84 ± 2.45). Finally, we enrolled 49 out of 198 ADNI controls who undergone structural scans, rsfMRI scans, and neuropsychological evaluation for subsequent analyses. Meanwhile, 32 MD-aMCI and 32 SD-aMCI patients entered the subsequent analyses.

## CSF Data

We downloaded the CSF dataset from the ADNI database. Based on CSF samples, Aβ1−42, total tau (t-tau) and phosphorylated tau (p-tau181) were measured (Bittner et al., 2016). These CSF biomarkers, including Aβ1−42, t-tau, and p-tau181, are also useful candidates as they are intimately related to amyloid plaques, neuronal death, and accumulation of tangles, respectively (Jovicich et al., 2016). Moreover, the ratio of p-tau/Aβ1−42 was also calculated due to its association with cognition and disease progression (Landau et al., 2010). Notably, not all subjects in the present study had the CSF samples due to the invasive procedure of lumbar puncture. Thus, the final CSF sample for analyses included 26 out of the 49 NC, 30 out of the 32 SD-aMCI, and 28 out of the 32 MD-aMCI patients.

### Amyloid PET Data

Given that amyloid PET shows more potent than CSF markers for MCI prognosis, we further downloaded the results of 18F-florbetapir PET data from the ADNI database, namely UCBERKELEYAV45\_11\_14\_17.csv (Bouallègue et al., 2017). The detailed processing procedure was described previously (Landau et al., 2013). In the following analysis, we used the results of composite SUVR value (intensity-normalized by a whole cerebellum region).

### Image Acquisition

The 3.0-Tesla Philips MRI scanner was used to scan all participants. Based on MPRAGE T1-weighted sequence, structural images were acquired with the following parameters: repetition time (TR) = 2300 ms; echo time (TE) = 2.98 ms; inversion time (TI) = 900 ms; 170 sagittal slices; within plane FOV = 256 × 240 mm2 ; voxel size = 1.1 × 1.1 × 1.2 mm3 ; flip angle = 9 ◦ ; and bandwidth = 240 Hz/pix. Based on echoplanar imaging sequence, rsfMRI images were acquired with the following parameters: time points = 140; TR = 3000 ms; TE = 30 ms; flip angle = 80 ◦ ; the number of slices = 48; slice thickness = 3.3 mm; spatial resolution = 3.31 × 3.31 × 3.31 mm3 ; and matrix = 64 × 64. Meanwhile, PET images were obtained by either a 30-min six frame scan acquired 30–60-min post-injection or a static 30-min single-frame scan acquired 30–60-min post-injection.

### Image Preprocessing

fnagi-10-00161 June 5, 2018 Time: 12:2 # 4

The image data were preprocessed by using the Data Processing Assistant for rsfMRI, DPASF2 , based on the Statistical Parametric Mapping (SPM 12)3 and rsfMRI Data Analysis Toolkit, REST4 . Briefly, we removed the first 10 image volumes of rsfMRI scans were for the signal equilibrium and the subjects' adaptation to the machine noise. The remaining 130 image volumes were corrected for both timing differences and head motion (24 Friston; Friston et al., 1996). We excluded the data with more than 2.0 mm maximum displacement in any directions (X-, Y-, and Z-axis) and 2.0 degrees of any angular motion. Subsequently, the T1-weighted images were co-registered to the mean rsfMRI image based on a thorough rigid-body transformation and normalized to the Montreal Neurological Institute (MNI) space, subsequent resliced into of 3 mm × 3 mm × 3 mm cubic voxels. Besides, we performed linear trends and temporal filter (0.01 Hz < f < 0.08 Hz). To remove any motion residual effects and other non-neuronal factors, we corrected covariates (including six head motion parameters, WM signal, CSF) in the following functional connectivity analysis. After preprocessing, we obtained each subject's 4D residual functional volume in native functional space. These 4D data were registered to the MNI152 space with 2 mm resolution (affine transformation). Due to the dispute of removing the global signal, we omit to regress out the global signal. Given the effect of the micro-motion artifact, we evaluated the framewise displacement (FD) value for each subject. Subjects were screened and excluded if FD value >0.5 mm on more than 35 volumes.

### Voxel-Mirrored Homotopic Connectivity Analysis

For VMHC computation, firstly, the mean T1 image was generated from the average of 113 spatially normalized T1 images. Then, we flipped the left hemisphere along the midline of the x-axis to obtain the symmetric brain template, further, to create the final template. Subsequently, the each subject's T1 image was co-registered nonlinearly to the customized symmetric template. The same deformation field subsequently applied to the rsfMRI. More details regarding VMHC data processing are available in the literature (Anderson et al., 2010; Zuo et al., 2010). Based on Pearson's correlation, the homotopic RSFC between any pair of symmetrical inter-hemispheric voxels was estimated, then transformed to Fisher's Z map. Finally, we defined these resultant values as the VMHC and used it in following analyses.

3http://www.fil.ion.ucl.ac.uk/spm/

## Corpus Callosum Volume

We evaluated the CC volume and its subregions based on FreeSurfer software package (Version 5.15 ) as described previously (Fischl, 2012; Reuter et al., 2012). In detail, FreeSurfer automatically segmented CC into five sections, with each section representing a fifth of the total area. We defined five segments as anterior (CC-A), mid-anterior (CC-MA), central (CC-C), midposterior (CC-MP), and posterior CC part (CC-P).

### FDG-PET Voxel-Wise Analysis

We downloaded the 18F FDG-PET data from the ADNI database in their most processed formats (PET Pre-processing protocol6 ). Specifically, pre-processed scans were generated following co-registration dynamic (to the first acquired frame), averaging frames, spatial re-orientation (AC–PC line), intensity normalization (within subject-specific mask), and smoothing (uniform isotropic spatial resolution 8 mm full width at half maximum kernel). The subsequent processing procedures of FDG-PET were conducted by combining T1 MRI images, which were similar to VMHC processing mentioned above. Firstly, each FDG-PET scan was co-registered to the individual's corresponding T1 MRI image in native space. Then, these structural scans were segmented, and the resulting GM partitions were spatially normalized using the SPM 12. The spatial normalization parameters (deformation fields) that had been estimated from the gray matter segment normalization were subsequently applied to the co-registered FDG-PET images so that both the GM, VMHC, and the FDG-PET images were in MNI space. To remove interindividual nuisance variability in tracer metabolism, we intensitynormalized FDG-PET image via dividing it by the average FDG-PET value of the reference region (cerebellar vermis, defined by the AAL regions within the MNI atlas, manually checked by experienced radiologists, MMZ). Finally, we created the standardized uptake value ratio (SUVR) images for the following analysis.

#### Statistics

All statistical analyses were performed using IBM SPSS statistical software. The TMT-A, TMT-B, and ADAS-cog were logtransformed because of a positively skewed distribution. We used ANOVA to detect group differences in terms of age, education level, cognitive ability, and AD-related pathological results. We used the post hoc pairwise t-tests if ANOVA was significant (p < 0.05, corrected by Bonferroni).

Regarding VMHC results, the ANCOVA was used to explore different brain regions among the three groups, with controlling for age, education level, and gender. For the objective to explore, the threshold set at height p < 0.05 and cluster level p < 0.05 (Gaussian random field, GRF corrected). We performed the post hoc pairwise t-tests within the brain regions identified by the ANCOVA (p < 0.005 Bonferroni corrected). Regarding the FDG-PET SUVR image, we performed ANCOVA to explore metabolic differences, with controlling for age, education level,

2http://www.rfmri.org/DPARSF

4http://www.resting-fmri.sourceforge.net/

5http://surfer.nmr.mgh.harvard.edu/

6http://adni.loni.usc.edu/methods/pet-analysis/pre-processing/

and gender [p < 0.001 at height and p < 0.05 at the cluster level (GRF corrected)]. Similarly, we performed the post hoc pairwise t-tests if ANCOVA was significant (p < 0.005, Bonferroni corrected). Also, we calculated the between-group CC subregions volume difference by using ANCOVA corrected by age, education level, gender, and total intracranial volume (TIV). We performed post hoc pairwise t-test if ANOVA was significant (p < 0.05, corrected by the least significant difference, LSD).

Further, we correlated imaging measures to cognitive abilities and also explored the possible interactions between neuroimaging modalities. To achieve the best visual effect, we overlapped between-group difference results of VMHC and FDG-PET SUVR image by the same statistical standard (ANCOVA, corrected by age, education level, and gender, p < 0.05 at height, p < 0.05 at the cluster level, GRF corrected). Meanwhile, to explore the possible interaction between CC subregions and other modalities, we performed voxel-wise regression analysis (p < 0.05 at height, p < 0.05 at the cluster level, GRF corrected). Notably, we calculated correlations limited to those indices having significant betweengroup differences.

### RESULTS

#### Patient Characteristics

We displayed quantitative variables as the mean and its SD, and the categorical variable as absolute and relative frequency. **Table 1** shows the demographic characteristics, cognitive abilities, and pathological biomarkers for each group. There were no significant differences in terms of age, education level, and gender composition among the three groups. Regarding the general cognitive ability, we found that three groups exhibited the significant differences in the score of MMSE and ADAS-cog, with MD-aMCI displaying the worst performance. Regarding other cognitive domains, group effects were significant, with the best performance in NC and also the worst performance in MD-aMCI patients.

Regarding AD-related pathological biomarkers, MD-aMCI patients exhibited a significantly higher composite PET PiB SUVR value than NC (p < 0.01). Meanwhile, although MDaMCI group presented a trend of the increased level of t-tau and p-tau181 as well as decreased Aβ1−42, no significant differences existed among groups (p > 0.05). Also, we also failed to detect the difference in the p-tau181/Aβ1−42 ratio among groups.

#### Glucose Metabolism Differences

Three groups demonstrated the significant metabolic differences in the following regions: namely bilateral middle temporal gyrus (MTG), bilateral inferior parietal lobe (IPL), and left PCu. Subsequent pairwise group comparisons carried out on these clusters revealed that both aMCI patients exhibited a lower metabolism in bilateral MTG. Also, the MD-aMCI demonstrated a reduced metabolism in left PCu and bilateral IPL compared to the SD-aMCI and NC (**Figure 1** and **Table 2**).

## Inter-Hemispheric Connection Differences

Among the three groups, we established the significant VMHC difference in these regions: namely middle MTG, PCu, postcentral gyrus (PCG), and calcarine gyrus (CG). As expected, subsequent pair-wise group comparisons suggested that MDaMCI exhibited a decreased VMHC in MTG, CG, PCu, and PCG relative to SD-aMCI and NC (p < 0.005, corrected by Bonferroni); however, there were no significant differences of VMHC between SD-aMCI and NC (**Figure 2** and **Table 3**).

We observed the between-group effect on CC-MP subregion (**Figure 1**). Then, the post hoc t-test revealed that MD-aMCI exhibited the smallest CC-MP volume among groups, and no significant volume difference existed between SD-aMCI and NC (p < 0.05, corrected by LSD).

#### Neuroimaging Indices Correlate Cognitive Abilities and Pathological Biomarkers

All the correlation analyses we performed were within aMCI patients (**Figure 3**). We observed that reduced left IPL metabolism (r = 0.27) was related to worse executive ability (p < 0.05). In addition, regarding functional connectivity, we identified that reduced inter-hemispheric connection in PCu (r = 0.39), PCG (r = 0.40), and CG (r = 0.54) were related to worse composite score of executive functioning (p < 0.001). Also, the CC-MP volume was related to the composite score of executive functioning (r = 0.40, p < 0.001).

Regarding AD-related pathological biomarkers, we analyzed the correlation only within the indices having the between-group difference. Our results revealed that hypometabolism of left PCu (r = −0.25, p < 0.05) and right IPL (r = −0.25, p < 0.05) was related to the increased amyloid accumulation reflecting composite 18F-florbetapir PET value.

### Interaction Between Modalities