Explainable Brain Age Prediction using coVariance Neural Networks
In computational neuroscience, there has been an increased interest in
developing machine learning algorithms that leverage brain imaging data to
provide estimates of "brain age" for an individual. Importantly, the
discordance between brain age and chronological age (referred to as "brain age
gap") can capture accelerated aging due to adverse health conditions and
therefore, can reflect increased vulnerability towards neurological disease or
cognitive impairments. However, widespread adoption of brain age for clinical
decision support has been hindered due to lack of transparency and
methodological justifications in most existing brain age prediction algorithms.
In this paper, we leverage coVariance neural networks (VNN) to propose an
explanation-driven and anatomically interpretable framework for brain age
prediction using cortical thickness features. Specifically, our brain age
prediction framework extends beyond the coarse metric of brain age gap in
Alzheimer's disease (AD) and we make two important observations: (i) VNNs can
assign anatomical interpretability to elevated brain age gap in AD by
identifying contributing brain regions, (ii) the interpretability offered by
VNNs is contingent on their ability to exploit specific eigenvectors of the
anatomical covariance matrix. Together, these observations facilitate an
explainable and anatomically interpretable perspective to the task of brain age
prediction.
Saurabh Sihag, Gonzalo Mateos, Corey McMillan, Alejandro Ribeiro
# Explainable Brain Age Prediction using coVariance Neural Networks

###### Abstract

In computational neuroscience, there has been an increased interest in developing machine learning algorithms that leverage brain imaging data to provide estimates of "brain age" for an individual. Importantly, the discordance between brain age and chronological age (referred to as "brain age gap") can capture accelerated aging due to adverse health conditions and therefore, can reflect increased vulnerability towards neurological disease or cognitive impairments. However, widespread adoption of brain age for clinical decision support has been hindered due to lack of transparency and methodological justifications in most existing brain age prediction algorithms. In this paper, we leverage coVariance neural networks (VNN) to propose an anatomically interpretable framework for brain age prediction using cortical thickness features. Specifically, our brain age prediction framework extends beyond the coarse metric of brain age gap in Alzheimer's disease (AD) and we make two important observations: (i) VNNs can assign anatomical interpretability to elevated brain age gap in AD by identifying contributing brain regions, (ii) the interpretability offered by VNNs is contingent on their ability to exploit specific eigenvectors of the anatomical covariance matrix. Together, these observations facilitate an explainable perspective to the task of brain age prediction.

## 1 Introduction

Aging is characterized by progressive changes in the anatomy and function of the brain [1] that can be captured by different modalities of neuroimaging [2; 3]. Importantly, individuals can age at variable rates, a phenomenon described as "biological aging" [4]. Numerous existing studies based on a large spectrum of machine learning approaches study brain-predicted biological age, also referred to as brain age, which is derived from neuroimaging data [5; 6; 7; 8; 9; 10; 11; 12]. Accelerated aging, i.e., when biological age is elevated as compared to chronological age (time since birth), may predict age-related vulnerabilities like risk for cognitive decline or neurological conditions like Alzheimer's disease (AD) [13; 14]. In this domain, the metric of interest is _brain age gap_, i.e., the difference between brain age and chronological age. We use the notation \(\Delta\)-Age to refer to the brain age gap.

Inferring \(\Delta\)-age from neuroimaging data presents a unique statistical challenge as it is essentially a qualitative metric with no ground truth and is expected to be elevated in individuals with underlying neurodegenerative condition as compared to the healthy population[15; 12]. The existing machine learning approaches often rely on models that have been trained to predict chronological age for a healthy population and then applied to a cohort with neurodegenerative condition. A layman overview of the procedure of inferring \(\Delta\)-Age is included in Appendix B. Several criticisms for brain age
evaluation approaches have been identified that stem from the lack of transparency in the machine learning models. Major criticisms include the coarseness of \(\Delta\)-Age that results in lack of specificity of brain regions contributing to the elevated brain age; and an unexplained reliance on the prediction accuracy for chronological age in the design of these machine learning models [5]. A significant portion of existing literature that studies brain age using deep learning models considers the ability of their models to accurately predict chronological age for healthy controls [16; 17; 18] as a relevant metric for assessing the quality of their methodological approach. Simultaneously, deep learning models that have a relatively moderate fit on the chronological age of healthy controls can also provide better insights into brain age than the ones with a tighter fit [17; 19]. Thus, there is a lack of conceptual clarity in the role of training to predict chronological age of healthy controls in predicting a meaningful \(\Delta\)-age [20].

In this paper, we propose a principled framework for brain age prediction using cortical thickness data that accommodates interpretability by isolating the contributing brain regions to elevated \(\Delta\)-Age in neurodegeneration. Our framework is based on the recently studied coVariance neural networks (VNNs) [21]. VNN is a graph convolutional network (GCN) with sample covariance matrix as the graph. The foundational analyses of VNNs in [21] showed that VNNs draw similarities with principal component analysis (PCA) while overcoming its significant drawbacks concerning reproducibility. Cortical thickness evolves with normal aging [22] and is impacted due to neurodegeneration [23; 24]. Thus, the age-related and disease severity related variations also appear in anatomical covariance matrices evaluated from the correlation among the cortical thickness measures across a population [25]. We focus our analysis on open access OASIS-3 dataset consisting of cortical thickness features from cognitively normal individuals and subjects in various stages of cognitive decline [26]. The utility of VNNs in predicting \(\Delta\)-Age has been explored previously in [27] but no insights were provided regarding their explainability. See Appendix A for other relevant studies.

**Contributions.** Our contributions in this paper are summarized as follows.

1. **VNNs provide anatomically interpretable \(\Delta\)-Age:**\(\Delta\)-Age in individuals with AD diagnosis was elevated as compared to healthy controls and significantly correlated with a clinical marker of dementia severity. Moreover, we could identify contributing brain regions to elevated \(\Delta\)-Age by analyses of the outputs at the final layer of the VNNs. Hence, by exploiting the VNN architecture, we could characterize \(\Delta\)-Age with anatomical interpretability (Fig. 2).
2. **Anatomical interpretability correlate with eigenvectors of the anatomical covariance matrix:** Our experiments demonstrated that certain eigenvectors of the anatomical covariance matrix were strongly correlated with the features that facilitated anatomical interpretability for \(\Delta\)-Age (Fig. 3). Thus, \(\Delta\)-Age was linked to the ability of VNNs to exploit specific eigenvectors of the anatomical covariance matrix.

The aforementioned findings also helped clarify the role of the preliminary step of training the VNNs to predict chronological age of a healthy population in \(\Delta\)-Age prediction. Specifically, this step equipped the VNNs with the ability to exploit the eigenvectors of the anatomical covariance matrix associated with elevated \(\Delta\)-Age or accelerated aging in AD-driven neurodegeneration.

## 2 coVariance Neural Networks

We begin with a brief introduction to VNNs. VNNs inherit the architecture of GCNs and operate on sample covariance matrix as the graph [21]. A dataset consisting of \(n\) random, independent and identically distributed (i.i.d) samples, given by \(\mathbf{x}_{i}\in\mathbb{R}^{m\times 1},\forall i\in\{1,\ldots,n\}\), can be represented in matrix form as \(\mathbf{X}_{n}=[\mathbf{x}_{1},\ldots,\mathbf{x}_{n}]\). Using \(\mathbf{X}_{n}\), the sample covariance matrix is estimated as

\[\mathbf{C}\triangleq\frac{1}{n-1}\sum_{i=1}^{n}(\mathbf{x}_{i}-\bar{\mathbf{x }})(\mathbf{x}_{i}-\bar{\mathbf{x}})^{\mathsf{T}}\;, \tag{1}\]

where \(\bar{\mathbf{x}}\) is the sample mean of samples in \(\mathbf{X}_{n}\). The covariance matrix \(\mathbf{C}\) can be viewed as the adjacency matrix of a graph representing the stochastic structure of the dataset \(\mathbf{X}_{n}\), where the \(m\) dimensions of the data can be thought of as the nodes of an \(m\)-node, undirected graph and its edges represent the pairwise covariances between different dimensions.
### Architecture

Similar to GCNs that rely on convolution operations modeled by _linear-shift-and-sum_ operators [28; 29], the convolution operation in a VNN is modeled by a coVariance filter, given by

\[\mathbf{H}(\mathbf{C})\triangleq\sum_{k=0}^{K}h_{k}\mathbf{C}^{k}\, \tag{2}\]

where scalar parameters \(\{h_{k}\}_{k=0}^{K}\) are referred to as filter taps. The application of coVariance filter \(\mathbf{H}(\mathbf{C})\) on an input \(\mathbf{x}\) translates to combining information across different sized neighborhoods. For \(K>1\), the convolution operation combines information across multi-hop neighborhoods (up to \(K\)-hop) according to the weights \(h_{k}\) to form the output \(\mathbf{z}=\mathbf{H}(\mathbf{C})\mathbf{x}\).

A single layer of VNN is formed by passing the output of the coVariance filter through a non-linear activation function \(\sigma(\cdot)\) (e.g., \(\mathsf{ReLU},\tanh\)) that satisfies \(\sigma(\mathbf{u})=[\sigma(u_{1}),\ldots,\sigma(u_{m})]\) for \(\mathbf{u}=[u_{1},\ldots,u_{m}]\). Hence, the output of a single layer VNN with input \(\mathbf{x}\) is given by \(\mathbf{z}=\sigma(\mathbf{H}(\mathbf{C})\mathbf{x})\). The construction of a multi-layer VNN is formalized next.

**Remark 1** (Multi-layer VNN).: _For an \(L\)-layer VNN, denote the coVariance filter in layer \(\ell\) of the VNN by \(\mathbf{H}_{\ell}(\mathbf{C})\) and its corresponding set of filter taps by \(\mathcal{H}_{\ell}\). Given a pointwise nonlinear activation function \(\sigma(\cdot)\), the relationship between the input \(\mathbf{x}_{\ell-1}\) and the output \(\mathbf{x}_{\ell}\) for the \(\ell\)-th layer is_

\[\mathbf{x}_{\ell}=\sigma(\mathbf{H}_{\ell}(\mathbf{C})\mathbf{x}_{\ell-1}) \quad\text{ for }\quad\ell\in\{1,\ldots,L\}\, \tag{3}\]

_where \(\mathbf{x}_{0}\) is the input \(\mathbf{x}\)._

Furthermore, similar to other deep learning models, sufficient expressive power can be facilitated in the VNN architecture by incorporating multiple input multiple output (MIMO) processing at every layer. Formally, consider a VNN layer \(\ell\) that can process \(F_{\ell-1}\) number of \(m\)-dimensional inputs and outputs \(F_{\ell}\) number of \(m\)-dimensional outputs via \(F_{\ell-1}\times F_{\ell}\) number of filter banks [30]. In this scenario, the input is specified as \(\mathbf{X}_{\text{in}}=[\mathbf{x}_{\text{in}}[1],\ldots,\mathbf{x}_{\text{in }}[F_{\text{in}}]]\), and the output is specified as \(\mathbf{X}_{\text{out}}=[\mathbf{x}_{\text{out}}[1],\ldots,\mathbf{x}_{\text{ out}}[F_{\text{out}}]]\). The relationship between the \(f\)-th output \(\mathbf{x}_{\text{out}}[f]\) and the input \(\mathbf{x}_{\text{in}}\) is given by \(\mathbf{x}_{\text{out}}[f]=\sigma\Big{(}\sum_{g=1}^{F_{\text{in}}}\mathbf{H}_{ fg}(\mathbf{C})\mathbf{x}_{\text{in}}[g]\Big{)}\), where \(\mathbf{H}_{fg}(\mathbf{C})\) is the coVariance filter that processes \(\mathbf{x}_{\text{in}}[g]\). Without loss of generality, we assume that \(F_{\ell}=F,\forall\ell\in\{1,\ldots,L\}\). In this case, the set of all filter taps is given by \(\mathcal{H}=\{\mathcal{H}^{\ell}_{fg}\},\forall f,g\in\{1,\ldots,F\},\ell\in\{ 1,\ldots,L\}\), where \(\mathcal{H}_{fg}=\{h^{\ell}_{fg}[k]\}_{k=0}^{K}\) and \(h^{\ell}_{fg}[k]\) is the \(k\)-th filter tap for filter \(\mathbf{H}_{fg}(\mathbf{C})\). Thus, we can compactly represent a multi-layer VNN architecture capable of MIMO processing via the notation \(\Phi(\mathbf{x};\mathbf{C},\mathcal{H})\), where the set of filter taps \(\mathcal{H}\) captures the full span of its architecture. We also use the notation \(\Phi(\mathbf{x};\mathbf{C},\mathcal{H})\) to denote the output at the final layer of the VNN. Various aspects of the VNN architecture are illustrated in Fig. 6 in Appendix C. The VNN final layer output \(\Phi(\mathbf{x};\mathbf{C},\mathcal{H})\) is succeeded by a readout function that maps it to the desired inference outcome.

### Properties of VNNs

The foundational analyses of VNNs in [21] established the following properties relevant to the data analysis in this paper.

**VNNs and PCA.** Given the eigendecomposition of \(\mathbf{C}=\mathbf{V}\mathbf{\Lambda}\mathbf{V}^{\mathsf{T}}\), the spectral properties of the coVariance matrix are established by studying the projection of the coVariance filter output \(\mathbf{z}=\mathbf{H}(\mathbf{C})\mathbf{x}\) on the eigenvectors \(\mathbf{V}\) (similar to that for a graph filter using graph Fourier transform [31; 32]). Theorem 1 in [21] establishes the equivalence between processing data samples with PCA and processing data samples with a specific polynomial on the covariance matrix \(\mathbf{C}\). Hence, it can be concluded that input data is processed with VNNs, at least in part, by exploiting the eigenvectors of \(\mathbf{C}\).

**Stability to perturbations in C.** VNNs are stable to perturbations in \(\mathbf{C}\)[21; Theorem 3]. This property implies that the inference performance achieved by VNN is likely to be reproducible, and hence, robust to the number of samples \(n\) used to estimate the covariance matrix \(\mathbf{C}\).

**VNNs are scale-free.** The filter taps in the coVariance filter in (2) are independent of the dimension of \(\mathbf{C}\). Hence, the covariance matrix \(\mathbf{C}\) can readily be replaced with another covariance matrix \(\mathbf{C}^{\prime}\) (not necessarily of the same dimensionality as \(\mathbf{C}\)) to process another dataset. By extension, the set of filter taps \(\mathcal{H}\) for a VNN is scale-free and can be used to process a dataset of arbitrary dimensionality. This
observation is particularly relevant for data analysis in neuroimaging, where the number of features is highly variable across datasets, but different datasets contain similar information [33].

We leverage the aforementioned properties of VNNs in the context of brain age to demonstrate the relationships of the inference outcomes with the eigenvectors of \(\mathbf{C}\), robustness of the results to the number of samples used to estimate \(\mathbf{C}\), and cross-validate various findings by leveraging datasets of different dimensionalities.

### VNN Learning

The VNN model is trained for a regression task using a dataset of \(m\) cortical thickness features and the chronological age for a population. Since the VNN architecture has \(F\) number of \(m\)-dimensional outputs in the final layer, \(\Phi(\mathbf{x};\mathbf{C},\mathcal{H})\) is of dimensionality \(m\times F\). The regression output is determined by a readout layer, which evaluates an unweighted mean of all outputs at the final layer of VNN. Therefore, the regression output for an individual with cortical thickness \(\mathbf{x}\) is given by

\[\hat{y}=\frac{1}{Fm}\sum_{j=1}^{m}\sum_{f=1}^{F}[\Phi(\mathbf{x};\mathbf{C}, \mathcal{H})]_{jf}\;. \tag{4}\]

Prediction using unweighted mean at the output implies that the VNN model exhibits permutation-invariance (i.e., the final output is independent of the permutation of the input features and covariance matrix) and VNN retains its scale-free property. Moreover, it allows us to associate a scalar output with each brain region among the \(m\) regions at the final layer. Specifically, we have

\[\mathbf{p}=\frac{1}{F}\sum_{f=1}^{F}[\Phi(\mathbf{x};\mathbf{C},\mathcal{H})] _{f}\;, \tag{5}\]

where \(\mathbf{p}\) is the vector denoting the mean of filter outputs in the final layer's filter bank. Note that the mean of all elements in \(\mathbf{p}\) is the prediction \(\hat{y}\) formed in (4). In the context of cortical thickness datasets, each element of \(\mathbf{p}\) can be associated with a distinct brain region. Therefore, \(\mathbf{p}\) is a vector of "regional contributions" to the output \(\hat{y}\) by the VNN. This observation will be instrumental to establishing the interpretability offered by VNNs in the context of \(\Delta\)-Age prediction in Section 3. For a regression task, the training dataset \(\{\mathbf{x}_{i},y_{i}\}_{i=1}^{n}\) (where \(\mathbf{x}_{i}\) are the cortical thickness data for an individual \(i\) with chronological age \(y_{i}\)) is leveraged to learn the filter taps in \(\mathcal{H}\) for the VNN \(\Phi(\cdot;\mathbf{C},\mathcal{H})\) such that they minimize the training loss, i.e.,

\[\mathcal{H}_{\text{opt}}=\min_{\mathcal{H}}\frac{1}{n}\sum_{i=1}^{n}\ell(\hat {y}_{i},y_{i})\;, \tag{6}\]

where \(\hat{y}_{i}\) is evaluated similarly to (4) and \(\ell(\cdot)\) is the mean-squared error (MSE) loss function.

## 3 Methods Overview for Brain Age Prediction

In this section, we provide an overview of the brain age prediction framework based on VNNs (see Fig. 1 for a summary). Our results primarily focus on the dataset described below.

**OASIS-3 Dataset**. This dataset was derived from publicly available freesurfer estimates of cortical thickness (derived from MRI images collected by 3 Tesla scanners and hosted on central.xnat.

Figure 1: Workflow for brain age and \(\Delta\)-Age prediction using VNNs.
org), as previously reported [26], and comprised of cognitively normal individuals (HC; \(n=611\), age = \(68.38\pm 7.62\) years, \(351\) females) and individuals with AD dementia diagnosis and at various stages of cognitive decline (\(n=194\), age = \(74.72\pm 7.02\) years, \(94\) females). The cortical thickness features were curated according to the Destrieux (DKT) atlas [34](consisting of \(m=148\) cortical regions). For clarity of exposition of the brain age prediction method, any dementia staging to subdivide the group of individuals with AD dementia diagnosis into mild cognitive impairment (MCI) or AD was not performed, and we use the label AD+ to refer to this group. The individuals in AD+ group were significantly older than those in HC group (t-test: \(p\)-value = \(2.46\times 10^{-23}\)). The boxplots for the distributions of chronological age for HC and AD+ groups are included in Fig. 14 in Appendix H. For \(191\) individuals in the AD+ group, the clinical dementia rating (CDR) sum of boxes scores evaluated within one year (365 days) from the MRI scan were available (CDR sum of boxes = \(3.45\pm 1.74\)). CDR sum of boxes scores are commonly used in clinical settings to stage dementia severity [35] and were evaluated according to [36].

The findings obtained via the analyses of the OASIS-3 dataset were validated using a smaller, independent dataset (described in Appendix G).

### Training the VNNs on HC group

We first train the VNN model to predict chronological age using the cortical thickness features for the HC group. This enables the VNN models to capture the information about healthy aging from the cortical thickness and associated anatomical covariance matrix. The hyperparameters for the VNN architecture and learning rate of the optimizer were chosen according to a hyperparameter search procedure using the package Optuna [37]. The VNN model had \(L=2\)-layers with a filter bank, such that we had \(F=5\), and \(6\) filter taps in the first layer and \(10\) filter taps in the second layer. The learning rate for the Adam optimizer was set to \(0.059\). The number of learnable parameters for this VNN was \(290\). The HC group was split into a \(90/10\) training/test split, and the covariance matrix was set to be the anatomical covariance evaluated from the training set. A part of the training set was used as a validation set and the other used for training the VNN model. We trained a set of 100 VNNs on different permutations of the training data. The training process was similar for all VNNs and is described in Appendix D. No further training was performed for the VNN models in the subsequent analyses.

### Analyses of regional residuals in AD+ and HC groups

Next, the VNN models trained to predict the chronological age for the HC group and (5) were adopted to study the effect of neurodegeneration on brain regions for AD+ group. Since the impact of neurodegeneration is expected to appear in the anatomical covariance matrix, we report the results when anatomical covariance \(\mathbf{C}_{\text{HA}}\) from the combined cortical thickness data of HC and AD+ groups was deployed in the trained VNN models. Because of the stability property of VNNs [21, Theorem 3], the inference drawn from VNNs is expected to be stable to the composition of combined HC and AD+ groups used to estimate the anatomical covariance matrix \(\mathbf{C}_{\text{HA}}\).

For every individual in the combined dataset of HC and AD+ groups, we processed their cortical thickness data \(\mathbf{x}\) through the model \(\Phi(\mathbf{x};\mathbf{C}_{\text{HA}},\mathcal{H})\) where parameters \(\mathcal{H}\) were learned in the regression task on the data from HC group as described previously. Hence, the vector of mean of all final layer outputs for cortical thickness input \(\mathbf{x}\) is given by \(\mathbf{p}=\frac{1}{F}\sum_{f=1}^{F}[\Phi(\mathbf{x};\mathbf{C}_{\text{HA}}, \mathcal{H}]_{f}\) and the VNN output is \(\hat{y}=\frac{1}{148}\sum_{j=1}^{148}[\mathbf{p}]_{j}\). Furthermore, we define the residual for feature \(a\) (or brain region represented by feature \(a\) in this case) as

\[[\mathbf{r}]_{a}\triangleq[\mathbf{p}]_{a}-\hat{y}. \tag{7}\]

Thus, (7) allows us to characterize the residuals with respect to the VNN output \(\hat{y}\) at the regional level, where brain regions are indexed by \(a\). Henceforth, we refer to the residuals (7) as "regional residuals". Recall that these are evaluated for an individual with cortical thickness data \(\mathbf{x}\).

In our experiments, for a given VNN model, the residual vector \(\mathbf{r}\) was evaluated for every individual in the OASIS-3 dataset. Also, the population of residual vectors for the HC group is denoted by \(\mathbf{r}_{\text{HC}}\), and that for individuals in the AD+ group by \(\mathbf{r}_{\text{AD+}}\). The length of the residual vectors is the same as the number of cortical thickness features (i.e., \(m=148\)). Since each element of the residual vector is associated with a distinct brain region, ANOVA was used to test for group differences between
individuals in HC and AD+ groups. Also, since elevation in \(\Delta\)-Age is the biomarker of interest in this analysis, we hypothesized that the brain regions that exhibited higher means for regional residuals for AD+ group than HC group would be the most relevant to capturing accelerated aging. Hence, the results are reported only for brain regions that showed elevated regional residual distribution in AD+ group with respect to HC group. Further, the group difference between AD+ and HC groups in the residual vector element for a brain region was deemed significant if it met the following criteria: i) the corrected \(p\)-value (Bonferroni correction) for the clinical diagnosis label in the ANOVA model was smaller than \(0.05\); and ii) the uncorrected \(p\)-value for clinical diagnosis label in ANCOVA model with age and sex as covariates was smaller than \(0.05\). See Appendix F for an example of this analysis.

Recall that 100 distinct VNNs were trained as regression models on different permutations of the training set of cortical thickness features from HC group. We used these trained models to establish the robustness of observed group differences in the distributions of regional residuals.

**Deriving regional profile for \(\Delta\)-Age from the robustness of findings from regional analyses.** We performed the regional analysis described above corresponding to each trained VNN model and tabulated the number of VNN models for which a brain region was deemed to be associated with a significantly elevated regional residual for the AD+ group. A higher number of VNN models isolating a brain region as significant suggested that this region was likely to be a highly robust contributor to accelerated aging in the AD+ group.

### Individual-level Brain Age Prediction

Finally, a scalar estimate for the brain age was obtained from the VNN regression output through a procedure consistent with the existing studies in this domain. Note that \(100\) VNNs provide \(100\) estimates \(\hat{y}\) of the chronological age for each subject. For simplicity, we consider \(\hat{y}\) to be the mean of these estimates. A systemic bias in the gap between \(\hat{y}\) and \(y\) may potentially exist when the correlation between \(\hat{y}\) and \(y\) is smaller than 1. Such a bias can confound the interpretations of brain age due to underestimation for older individuals and overestimation for younger individuals [38]. Therefore, to correct for this age-induced bias in \(\hat{y}-y\), we adopted a linear regression model-based approach [38, 39]. Specifically, the following bias correction steps were applied on the VNN estimated age \(\hat{y}\) to obtain the brain age \(\hat{y}_{\text{B}}\) for an individual with chronological age \(y\):

**Step 1.** Fit a regression model for the HC group to determine scalars \(\alpha\) and \(\beta\) in the following model:

\[\hat{y}-y=\alpha y+\beta. \tag{8}\]

**Step 2.** Evaluate brain age \(\hat{y}_{\text{B}}\) as follows:

\[\hat{y}_{\text{B}}=\hat{y}-(\alpha y+\beta). \tag{9}\]

The gap between \(\hat{y}_{\text{B}}\) and \(y\) is the \(\Delta\)-Age and is defined below. For an individual with cortical thickness \(\mathbf{x}\) and chronological age \(y\), the brain age gap \(\Delta\)-Age is formally defined as

\[\Delta\text{-Age}\triangleq\hat{y}_{\text{B}}-y\, \tag{10}\]

where \(\hat{y}_{\text{B}}\) is determined from the VNN age estimate \(\hat{y}\) and \(y\) according to steps in (8) and (9). The age-bias correction in (8) and (9) was performed only for the HC group to account for bias in the VNN estimates due to healthy aging, and then applied to the AD+ group. Further, the distributions of \(\Delta\)-Age were obtained for all individuals in HC and AD+ groups.

\(\Delta\)-Age for AD+ group was expected to be elevated as compared to HC group as a consequence of elevated regional residuals derived from the VNN model. To elucidate this, consider a toy example with two individuals of the same chronological age \(y\), with one belonging to the AD+ group and another to the HC group. Equation (9) suggests that their corresponding VNN outputs (denoted by \(\hat{y}_{\text{AD+}}\) for individual in the AD+ group and \(\hat{y}_{\text{HC}}\) for individual in the HC group) are corrected for age-bias using the same term \(\alpha y+\beta\). Hence, \(\Delta\)-Age for the individual in the AD+ group will be elevated with respect to that from the HC group only if the VNN prediction \(\hat{y}_{\text{AD+}}\) is elevated with respect to \(\hat{y}_{\text{HC}}\). Since the VNN predictions \(\hat{y}_{\text{AD+}}\) and \(\hat{y}_{\text{HC}}\) are proportional to the unweighted aggregations of the estimates at the regional level [see (4) and (5)], larger \(\hat{y}_{\text{AD+}}\) with respect to \(\hat{y}_{\text{HC}}\) can be a direct consequence of a subset of regional residuals [see (7)] being robustly elevated in AD+ group with respect to HC group across the \(100\) VNN models. When the individuals in this example have different chronological age, the age-bias correction is expected to remove any variance due to chronological age in \(\Delta\)-Age. We also verified that the differences in \(\Delta\)-Age for AD+ and HC group were not driven by age or sex differences via ANCOVA with age and sex as covariates.
Results

### Chronological age prediction for the HC group

The performance achieved by the VNNs on the chronological age prediction task for the HC group has been summarized over the \(100\) nominal VNN models. VNNs achieved a mean absolute error (MAE) of \(5.82\pm 0.13\) years and Pearson's correlation of \(0.51\pm 0.078\) between the chronological age estimates and the ground truth on the test set. Moreover, on the complete dataset, the MAE was \(5.44\pm 0.18\) years and Pearson's correlation was \(0.47\pm 0.074\). Thus, the trained VNNs were not overfit on the training set.

Next, for every individual in the HC group, we evaluated the mean of the inner products (also equivalently referred to as dot product) between the vectors of contributions of every brain region [**p** in (5)] and the eigenvectors of the anatomical covariance matrix for all \(100\) VNN models. The strongest alignment was present between the first eigenvector of the anatomical covariance matrix (i.e., the eigenvector associated with the largest eigenvalue) and the vectors of regional contributions to the VNN output (\(0.987\pm 0.0005\) across the HC group), with relatively smaller associations for second (\(0.051\pm 0.003\)), third (\(0.075\pm 0.004\)), and fourth (\(0.094\pm 0.003\)) eigenvectors. Additional details are included in Appendix E. Thus, the VNNs exploited the eigenvectors of the anatomical covariance matrix to achieve the learning performance in this task. The first eigenvector of the anatomical covariance matrix predominantly included bilateral anatomic brain regions in the parahippocampal gyrus, precuneus, inferior medial temporal gyrus, and precentral gyrus. These findings were cross-validated on the HC group by the VNNs that had been trained on a different dataset (Appendix G).

We remark that several existing studies on brain age prediction have utilized deep learning and other approaches to report better MAE on their respective healthy populations [16; 17; 18; 40]. In contrast, our contribution here is conceptual, where we have explored the properties of VNNs when they are trained to predict the chronological age for the HC group. Subsequently, our primary focus in the context of brain age is on demonstrating the anatomical interpretability offered by VNNs and relevance of eigenvectors of the anatomical covariance matrix. Thus, we further provide the insights that have not been explored (or are infeasible to obtain) for most existing brain age evaluation frameworks based on less transparent deep learning models.

### Analyses of regional residuals derived from VNNs reveal regions characteristic of AD

Figure 1(a) displays the robustness (determined via analyses of 100 VNN models) for various brain regions being associated with significantly larger residual elements for the AD+ group than the HC group. The most significant regions with elevated regional residuals in AD+ with respect to HC were concentrated in bilateral inferior parietal, temporal, entorhinal, parahippocampal, subcallosal, and precentral regions. All these brain regions, except for precentral and subcallosal, mirrored the cortical atrophy (Fig. 15 in supplementary material), and these regions are known to be highly interconnected with hippocampus [41]. Hence, brain regions characteristic of AD had significant differences in regional residual distributions for AD+ group as compared to HC group. Similar findings were recovered by the analysis of the regional residuals derived from the VNNs that were trained on an independent dataset (Appendix G).

Although the results in Fig. 1(a) provided a meaningful regional profile for AD+ group, we further performed exploratory analysis to check whether the regional residuals had any clinical significance. To this end, we evaluated the Pearson's correlations between CDR sum of boxes and the regional residuals derived from final layer VNN outputs for the AD+ group for all 100 VNN models. Interestingly, the brain regions with the largest correlations with the CDR sum of boxes scores in Fig. 1(c) were concentrated in the parahippocampal, medial temporal lobe, and temporal pole regions (also isolated in Fig. 1(a)). This observation provides the evidence that the regional residuals for the AD+ group that led to the result in Fig. 1(a) could predict dementia severity.

### \(\Delta\)-Age is elevated in AD+ group and correlated with CDR

We evaluated the \(\Delta\)-Age for HC and AD+ groups according to the procedure specified in Section 3.3. We also investigated the Pearson's correlation between \(\Delta\)-Age and CDR sum of boxes scores in AD+ group. Figure 1(b) illustrates the distributions for \(\Delta\)-Age for HC and AD+ groups (\(\Delta\)-Age for HC: \(0\pm 2.83\) years, \(\Delta\)-Age for AD+: \(3.54\pm 4.49\) years). The difference in \(\Delta\)-Age for AD+ and
HC groups was significant (Cohen's \(d\) = 0.942, ANCOVA with age and sex as covariates: \(p\)-value \(<10^{-20}\), partial \(\eta^{2}=0.158\)). Also, age and sex were not significant in ANCOVA (\(p\)-value \(>0.4\) for both). Hence, the group difference in \(\Delta\)-Age for the two groups was not driven by the difference in the distributions of their chronological age. Figure 2d plots CDR sum of boxes scores versus \(\Delta\)-Age for the AD+ group. Pearson's correlation between CDR sum of boxes score and \(\Delta\)-Age was \(0.474\) (\(p\)-value \(=2.88\times 10^{-12}\)), thus, implying that the \(\Delta\)-Age evaluated for AD+ group captured information about dementia severity. Hence, as expected, the \(\Delta\)-Age for AD+ was likely to be larger with an increase in CDR sum of boxes scores. For instance, the mean \(\Delta\)-Age for individuals with CDR sum of boxes greater than 4 was \(6.04\) years, and for CDR sum of boxes \(\leq\) 4 was \(2.42\) years.

Given that the age-bias correction procedure is a linear transformation of VNN outputs, it can readily be concluded that the statistical patterns for regional residuals in Fig. 2a and Fig. 2c lead to elevated \(\Delta\)-Age and correlation between \(\Delta\)-Age and CDR sum of boxes scores. Therefore, our framework provides a feasible way to characterize accelerated biological aging in AD+ group with a regional profile. Additional figures and details pertaining to VNN outputs and brain age before and after the age-bias correction was applied are included in Appendix H.

Regional residuals derived from VNNs trained on OASIS-3 are correlated with eigenvectors of the anatomical covariance matrix

We further investigated the relationship between regional residuals derived from VNNs and the eigenvectors of \(\mathbf{C}_{\text{HA}}\) to determine whether any specific eigenvectors (principal components) of \(\mathbf{C}_{\text{HA}}\) were instrumental to recover the findings in Fig. 2b. For this purpose, we evaluated the inner product of normalized residual vectors (norm = 1) obtained from VNNs and the eigenvectors of the covariance matrix \(\mathbf{C}_{\text{HA}}\) for the individuals in AD+ group. The normalized residual vector is denoted by \(\mathbf{\bar{r}}_{\text{AD+}}\). For every individual, the mean of the absolute value of the inner product \(|<\mathbf{\bar{r}}_{\text{AD+}},\mathbf{v}_{i}>|\) (where \(\mathbf{v}_{i}\) is the \(i\)-th eigenvector of \(\mathbf{C}_{\text{HA}}\)) was evaluated for the 100 VNN models that were used to derive the results in Fig. 2. Figure 3 plots the mean inner product for eigenvectors associated with \(50\) largest eigenvalues of \(\mathbf{C}_{\text{HA}}\). The three largest mean correlations with the regional residuals in AD+ group were observed for the third eigenvector of \(\mathbf{C}_{\text{HA}}\) (\(|<\mathbf{\bar{r}}_{\text{AD+}},\mathbf{v}_{3}>|=0.645\pm 0.016\)), fourth eigenvector (\(|<\mathbf{\bar{r}}_{\text{AD+}},\mathbf{v}_{4}>|=0.305\pm 0.02\)), and the first eigenvector (\(|<\mathbf{\bar{r}}_{\text{AD+}},\mathbf{v}_{1}>|=0.299\pm 0.001\)). Similar correlation patterns between the regional residuals and the eigenvectors

Figure 2: **Interpretable \(\Delta\)-Age evaluation in OASIS-3 dataset.** Panel **a** displays the robustness of the significantly elevated regional residuals for AD+ group with respect to HC group for different brain regions. For every VNN model in the set of \(100\) nominal VNN models that were trained on HC group, the analyses of regional residuals helped isolate brain regions that corresponded to significantly elevated regional residuals in AD+ group with respect to HC group. After performing this experiment for \(100\) VNN models, the robustness of the observed significant effects in a brain region was evaluated by calculating the number of times a brain region was identified to have significantly elevated regional residuals in AD+ group with respect to HC group. The number of times a brain region was linked with significantly elevated regional residuals in AD+ group with respect to HC group is projected on the brain template. Panel **b** displays the distribution of \(\Delta\)-Age for HC and AD+ groups. The elevated brain age effect here is characterized by regional profile in Panel **a**. Panel **c** projects the mean Pearson’s correlation between regional residuals (derived for each VNN model in the set of 100 nominal VNN models) and CDR sum of boxes for AD+ group on the brain template. Panel **d** displays the scatter plot for CDR sum of boxes versus \(\Delta\)-Age in AD+ group. The correlation between \(\Delta\)-Age and CDR sum of boxes could be attributed to the observations made in Panel **c**.
of \(\mathbf{C_{HA}}\) for the AD+ group were recovered by the VNNs that had been trained on an independent dataset and used to process OASIS-3 dataset (Appendix G). These eigenvectors are plotted on a brain template in the expanded Fig. 9b in the supplementary material. Inspection of the first, third, and fourth eigenvector of \(\mathbf{C_{HA}}\) suggested that subcallosal, entorhinal, parahippocampal and temporal pole regions could be the most relevant contributors to the elevated \(\Delta\)-Age in AD+ group in Fig. 2.

### Additional Experiments

**Cross-validation:** Using the VNNs trained to predict chronological age for an independent dataset of healthy population (cortical thickness curated according to a different brain atlas), we recovered a regional profile similar to Fig. 2a on the OASIS-3 dataset. Furthermore, the regional residuals derived from this set of VNNs had similar correlation patterns with the eigenvectors of \(\mathbf{C_{HA}}\). This observation and Fig. 3 led to the conclusion that training to predict chronological age of a healthy population was instrumental for the VNNs to gain the ability to exploit eigenvectors of the anatomical covariance matrix that are relevant to brain age prediction. See Appendix G for details.

**Stability to perturbations in \(\mathbf{C_{HA}}\):** As a consequence of the stability of VNNs, we observed that the regional profile for \(\Delta\)-Age in Fig. 2a was stable even when the covariance matrix \(\mathbf{C_{HA}}\) was estimated by a variable composition of individuals from the HC and AD+ group (Appendix J).

**Anatomical covariance matrix and brain age:** Use of anatomical covariance matrix derived from only HC group provides results consistent with Fig. 2, albeit with a slightly smaller group difference between the \(\Delta\)-Age distributions for HC and AD+ groups. See Appendix K for details.

## 5 Discussion

Our study has provided a foundational contribution to the methodology of brain age prediction. Anatomical interpretability of \(\Delta\)-Age was derived from the features extracted at the final layer of the VNNs and these features were correlated with marker of dementia severity and certain eigenvectors of the anatomical covariance matrix. Importantly, training the VNNs to predict chronological age helped fine tune their parameters to exploit the relevant eigenvectors of the anatomical covariance matrix. Thus, the role of the age-bias correction step was restricted to projecting the VNN outputs onto a space where one could infer accelerated biological aging with respect to the chronological age from a layman's perspective. By associating \(\Delta\)-Age with a regional profile, VNNs also provide a feasible tool to distinguish pathologies if the distributions of \(\Delta\)-Age for them are overlapping. Moreover, a near-perfect chronological age prediction for healthy controls by itself may not be a determinant of the quality of brain age prediction in neurodegeneration. A larger focus is needed on principled statistical approaches for brain age prediction that can capture the factors that lead to accelerated aging. Locally interpretable and theoretically grounded deep learning models such as VNNs can provide a feasible, promising future direction to build statistically and conceptually legitimate brain age prediction models in broader contexts. Incorporating other modalities of neuroimaging or alternative metrics of aging other than chronological age (such as DNA methylation aging [42]) in the training of VNNs provide promising future directions that can help improve our understanding of aging.

**Limitations.** We provide findings derived on a single dataset (although corroborated using VNNs trained on an independent dataset) and further investigation on larger and more diverse datasets may improve the confidence in our conclusions. Existing studies, including this paper, fall short at

Figure 3: Bar plots of the mean inner products between the normalized vector of regional residuals (norm = 1) of VNN outputs (VNNs trained on OASIS-3) obtained from AD+ group and the eigenvectors of \(\mathbf{C_{HA}}\) (covariance matrix of combined HC and AD+ group) with respective standard deviations as error bars. Results with coefficient of variation \(>30\%\) across the AD+ group have been excluded.
concretely defining the notion of optimal brain age. From a broader perspective, quantifying biological age even for a healthy population is a complex task due to various factors that can contribute to accelerated aging in the absence of an adverse health condition [43; 44; 45]. Moreover, the impacts of the quality of MRI scans and brain atlases across datasets on \(\Delta\)-Age must also be explored.

## 6 Data and Code Availability

OASIS-3 dataset is publicly available and hosted on central.xnat.org. Access to OASIS-3 dataset may be requested through [https://www.oasis-brains.org/](https://www.oasis-brains.org/). Code for demonstrating brain age evaluation is available at [https://github.com/pennbindlab/VNN_Brain_Age](https://github.com/pennbindlab/VNN_Brain_Age). Requests for details regarding IDs of individuals in OASIS-3 and source data for all figures may be sent to sihags@pennmedicine.upenn.edu.

## 7 Acknowledgements

OASIS-3 data were provided by Longitudinal Multimodal Neuroimaging: Principal Investigators: T. Benzinger, D. Marcus, J. Morris; NIH P30 AG066444, P50 AG00561, P30 NS09857781, P01 AG026276, P01 AG003991, R01 AG043434, UL1 TR000448, R01 EB009352. The MRI data for XYZ dataset was provided by the Penn Frontotemporal Degeneration Center (NIH AG066597) and Penn Institute on Aging. Cortical thickness data were made available by Penn Image Computing and Science Lab at University of Pennsylvania.
