In this research, we used a dataset composed of 432 brain T1 MRI images, relative to 126 normal control (NC), 214 MCI and 92 AD subjects from the Alzheimer's Disease Neuroimaging Initiative (ADNI). ADNI is a longitudinal multicenter study designed to obtain early diagnosis and monitoring of AD through the study and development of clinical, imaging, biochemical, and genetic biomarkers (https://adni.loni.usc.edu/about/). ADNI images were normalized using the MNI152 brain template with 197 × 233 × 189 mm³ size and 1 × 1 × 1 mm³ resolution; from now onward voxel and 1 mm³ will be interchangeably used. In the following Table 1 the number of instances, gender information, age, years of education, and Mini Mental State Examination (MMSE) score are enlisted. This cohort of MRI scans passed quality control (QC) using the Laboratory of Neuro Imaging (LONI) QC System (Kim et al., 2019) so that only high quality images were considered without signal alterations.

The proposed approach involves two main stages: the construction of the network model and the learning and explainability phase (see Figure 1). First, a graph theory-based model is employed to study structural brain connectivity and identify early signs of Alzheimer's disease. Accordingly, each patient's brain is modeled as a complex network. Brain MRIs are parceled into rectangular boxes called “patches” of fixed dimensions that represent the network's nodes; the links are obtained in terms of nodes' pairwise similarity measured by absolute Pearson's correlation. As a difference with previous works where multiplex networks were adopted, here we consider single-subject networks to achieve a simpler description and, therefore, a more interpretable model. In fact, using single-subject networks, network features can be directly related to single nodes and then to specific brain regions.

Image processing was the first algorithmic step. Using Oxford FMRIB software library (FSL) (Jenkinson et al., 2012), images were skull-stripped and spatially normalized as well in intensity to mitigate data heterogeneity, an aspect of fundamental importance for studying a multi-center database such as the one provided by ADNI. First bias field correction and skull stripping was performed using FSL Brain Extraction Tool (BET) (Smith, 2002). Thereafter, spatial normalization was performed to ensure co-registration to the MNI152 template using FSL Linear Registration Tool (FLIRT) (Jenkinson et al., 2002); in particular, an affine registration was adopted with default parameters.

After registration, by using the medial longitudinal fissure, normalized brains were divided into the two hemispheres; each hemisphere was then covered by an equal number of patches of fixed dimensions l₁×l₂×l₃ (patches overlapping the template with < 10% of voxels were neglected). For the patch dimension, we considered the 10 × 15 × 20 mm³ configuration, corresponding to a total volume of 3, 000 mm³, for a total of 549 patches (see Figure 2).

Accordingly, we built an undirected weighted network for each MRI scan with 549 nodes. Pearson's correlation coefficient was chosen for links for three main reasons: (i) it has an affordable computational cost; (ii) it is easy to implement, and (iii) it easy interpret in terms of brain atrophy. Absolute values were considered to take into account left/right symmetry of the brain; finally, a threshold was used to remove weak correlations (< 0.3) that could raise noisy connections.

Once structural connectivity had been modeled in terms of a complex network, we exploited such characterization to evaluate some network metrics and outline the effects of cerebral atrophy on brain connectivity. Furthermore, since network metrics can be easily divided in local and global ones, to acquire a detailed description of how AD differently affects distinct anatomical districts, nodal centrality measures were preferred. In particular, three graph metrics were used: strength, betweenness, and eigenvector centrality. The strength of a node is a local centrality measure defined by the sum of edge weights of each node; as a difference, a global centrality measure such as eigenvector centrality takes into account also the influence of a node in a network based on the amount of nodes with high number of connections which it is connected to. Finally, a dynamic centrality measure is considered: the betweenness; betweenness is the ratio between the shortest paths connecting two nodes passing through a specific node and all available shortest paths. Accordingly, it is a measure capturing the information flow within a network (Amoroso et al., 2021; Bellantuono et al., 2021; Sheng et al., 2021).

After calculating these three metrics for each node of the adjacency matrices, a matrix representation M×n of the data was obtained, where M is the number of enrolled patients and n is the total number of features (number of metrics used times number of nodes). Here, the resulting matrix had dimensions of 432 × 1, 647.

We have used the previously defined matrix representation to train a machine learning model. The main objectives we set during this learning phase were to achieve a sound representation of AD and verify the robustness and reliability of the connectivity network model. To this aim, a three-label classification was performed for AD, MCI, and NC subjects. All these analyses were carried out within a 10-fold cross-validation framework; the procedure was iterated 50 times to estimate performances and uncertainties; besides, to ensure balanced cross-validation splits between training and test, stratified sampling was used.

To ensure that the observed performance was due to the informative content provided by the network features independently from the classification model used, we compared the performance of several classifiers: Random Forest (RF) (Breiman, 2001), Support Vector machine (SVM) (Cortes and Vapnik, 1995), eXtreme Gradient Boosting (XGBoost) (Chen et al., 2015), Naive Bayes (NB) (Rish et al., 2001), and Logistic Regression (LR) (Ng and Jordan, 2001) were used.

After classification performance accuracy and reliability had been assessed, explainability analyses were carried out by means of SHapley Additive exPlanations (SHAP) (Lundberg and Lee, 2017). Accordingly, it is possible to evaluate to which extent each feature affects the model's predictions. By averaging the impact of a feature on the whole dataset, a feature importance map can be achieved; at the same time, this explanation allows for the identification of the factors which determine the classifier's decisions for each subject. Also, this analysis provides additional information compared with standard feature importance approaches in that it does not determine the importance of one feature but also it explains how its value, combined with the other features' values, leads to the decision taken by the algorithm. Moreover, it is possible to compare whether the impact of features on a single decision is coherent with the overall feature importance (the one related to the three-class classification model) or not. Finally, an important aspect to remark is that SHAP can be adopted with different machine learning algorithms (Strumbelj and Kononenko, 2010), thus making the proposed approach algorithm-independent.

A fundamental aspect of SHAP concerns their computational burden, as the computational load increases exponentially with the number of features. To overcome this issue, a specific experimental design is proposed. Explainability analyses were nested in the repeated 10-fold cross-validation. Using a RF classifier, for each cross-validation round a subset of important features was selected by means of the combined use of Mean Decrease Accuracy (MDA) and statistical significance.

This procedure can be described in the following two steps: (i) firstly, the features exceeding the 95th percentile of the MDA distribution were selected within each cross-validation round; (ii) then, these newly selected features were further reduced considering statistical significance: only the features with the 5% significance using the one-sided binomial test were kept. On average, the number of features selected during each cross-validation round was ~50. During each cross-validation round, this final set of important features was used to train a second classification model on the same training set with reduced number of features, in order to evaluate its informative content and to carry out the explainability analyses. Shapley values were calculated for each of the test subjects of each cross-validation round. The same analysis was carried out for each available class by considering three one-vs.-all classifiers to outline the presence of class-specific behaviors (Whitwell et al., 2007; Byun et al., 2015; Cabral et al., 2015). The mentioned analyses were carried out with R, with the DALEX and shapviz packages (Biecek, 2018; Biecek and Burzykowski, 2020; Michael Mayer, 2022). Finally, to validate the results from a clinical perspective, we determined which specific brain regions were associated with the selected important features (as each feature was directly related to a patch). To this aim, the Talairach atlas (Lancaster et al., 1997, 2000) was used.

Firstly, we investigated the validity and reliability of the brain connectivity model. To this aim, a 10-fold cross-validation analysis was repeated 50 times, results are shown in Figure 3.

These findings demonstrate that, except for LR, the adopted models are consistent and the informative power of the network features can reach satisfactory values: in terms of accuracy, RF resulted the best model with a mean three-class accuracy of 0.66 ± 0.07.

To obtain class-specific evaluations, ROC curves and their area under the curve (AUC) were taken into account (see Figure 4).

The results show an accurate classification for all classes: 0.81 ± 0.02 for AD, 0.77 ± 0.01 for MCI and 0.83 ± 0.01 for NC. The contingency table allows to appreciate how MCI and AD include most of misclassifications, 50 and 46, respectively. We also evaluated the agreement of the predictions of the best performing classifiers: RF and SVM by means of Pearson's correlation. We found a 0.79 correlation for AD, 0.70 for MCI and 0.78 for NC. Hence, for subsequent analyses, only the RF model was considered.

Once demonstrated the reliability of the base of knowledge, we investigated whether the features driving classification remained the same despite class-specific differences. In Figure 5, it is possible to observe the overall (related to all three classes) ranking of features; only the top 20 features are shown for display purposes.

The smooth trend demonstrates that it is difficult to find few features dominating over the others. Interestingly, the top three positions were occupied by betweenness; however, the significant presence of both strength and eigenvector centrality features suggested that no metrics prevailed over the others.

Using the three class-specific models, three importance rankings were also constructed; besides, to outline class-related patterns, explainability analyses were carried out using only correctly classified subjects. These class-rankings were compared with the overall ranking by means of the Spearman's correlation coefficient ρ. The results are presented in Table 2.

While MCI and NC rankings were significantly correlated to the overall ranking, the AD ranking showed a very weak correlation. These findings allow us to highlight how much brain connectivity is influenced by brain atrophy; particularly in AD patients, in which the high heterogeneity is more impactful than in subjects of the other two clinical classes.

Lastly, we analyzed the correlation among the three feature importance rankings restricted to only correctly classified subjects (see Table 3).

The degree of correlation between these rankings allows us to understand to what extent the correct classification of subjects of different classes may depend on the order of importance of the features. Using Spearman's correlation, it is possible to determine the strength and direction of the monotonic relationship between two ranks. Although the observed correlations are definitely lower than the previous ones, at 5% significance anti-correlation between AD vs. MCI and AD vs. NC is detected. These negative Spearman correlation coefficients highlight a decreasing monotonic trend among the ranks of the features that better correctly classify the AD class compared to those of the other two clinical classes. These correlations, along with the previous ones, suggest the presence of coherent and distinguishable patterns for NC and MCI, whereas a more heterogeneous and elusive condition characterizes the AD class.

SHAP were also used to examine the effects of features on predictions: how and how much each feature affected the prediction score.

This analysis provided additional information compared to the simple rankings showed in the previous section: the features are ranked from top to bottom, from the most important to the least important, and they are distributed based on both their impact on the prediction as well as their value, for high values, colors tend toward yellow, and for low values, colors tend toward purple. A positive Shapley value represents a positive impact on the correct predicted outcome, and a negative value implies a negative impact (see Figure 6).

Observing the beeswarm plots, the features in common between the three cases are of particular interest. In particular, for AD and MCI the betweenness of patches 364, 368, and 380 were found along with the strength of patches 374 and 380. Analogously, for AD and NC, strength of patches 27 and 61, eigenvector centrality of patches 27, 301, and 304 and betweenness of patch 293 were highlighted. Finally, MCI and NC showed in common: betweenness of patches 188, 255, 298, 313, 384, and 420; eigenvector centrality of patches 247, 255, 305, 417, and 523; strength of patches 42, 247, 255, 308, and 417. Notably, in some cases, the same patch was detected from different metrics.

Let us consider some specific cases. For AD and MCI classes, it is possible to observe the opposite behaviors of betweenness of patches 364 and 380. Thus, these features are important for the correct classification of AD subjects if their values are high; conversely, low values of these features are significant for the correct classification of MCI subjects. Similar findings can be obtained for other classes, for example, NC subjects are characterized by low values of eigenvector centrality of patch 27 while AD subjects tend to show are higher values. Interestingly, the importance of a feature for classification can differ from one class to another, for example the betweenness of patch 293 ranks third for NC subjects, but it is ranked last for AD patients. Understanding this findings without the help of anatomical interpretation can be extremely complicated, this is why in the following section the XAI visualization is presented. Another aspect which deserves to be outlined is the possibility to carry out personalized XAI analyses, which explain for each patient how classification was determined (see Figure 7).

The Figure 7 shows three correctly classified subjects belonging to AD, MCI, and NC, respectively, which correspond to three example cases to show possible differences with global XAI analysis. These graphs are waterfall plots and provide an analysis tool to visualize how each feature impacts on the average of the model output over the training data (baseline value E[f(x)]) for a specific subject. The bars represent the Shapley values of the features that most influentially drive the individual prediction. As can be seen from the three plots, positive values indicate a greater impact on the prediction, increasing the baseline predicted score. Longer bars indicate a greater influence of the feature and adding them together we obtain the final prediction f(x); the Figure 7 shows how this value has increased the initial baseline value. It is apparent from these graphs that some features extremely important from a global perspective (e.g., eigenvector centrality of patch 27 or betweenness of patch 61) and that are essential to the global analysis of classes, do not appear in the class-specific analyses. While global XAI explains on average which features are important for classification, subject-level XAI is of fundamental importance to clinically validate what happens for a specific patient; in anticipation of future developments toward personalized diagnostics.

Finally, to clinically validate the methodology we associated an anatomical district to each important feature (see Figure 8).

The goal was to show the model consistency by demonstrating that the outlined regions are clinically related to AD. For what concerns the overall ranking, the brain regions found were putamen, middle and superior temporal gyrus, anterior cingulate, precentral and postcentral gyrus, insula, sub-gyral, thalamus, culmen, lingual gyrus, cuneus, middle occipital gyrus, and brain stem. Analogously, we found the most relevant brain regions for the class-specific rankings.

There are a number of brain regions that were included in the best 30 for AD, including the amygdala, parahippocampal gyrus, fusiform gyrus, precuneus, posterior cingulate, as well as other regions previously ranked in global rankings (insula and lingual gyrus) that represent additional features nodes. For the MCI class, the following brain regions have been included in the top 30: posterior cingulate, precuneus and others already found in the global ranking, such as middle temporal gyrus, sub-gyral, and insula. The following regions are ranked in the top 30 for the NC class: uncus, parahippocampal gyrus, lateral ventricle, caudate, precuneus, posterior cingulate, and other regions already found in the global ranking (culmen, thalamus).