A Comparative Study of Population-Graph Construction Methods and Graph
  Neural Networks for Brain Age Regression
The difference between the chronological and biological brain age of a
subject can be an important biomarker for neurodegenerative diseases, thus
brain age estimation can be crucial in clinical settings. One way to
incorporate multimodal information into this estimation is through population
graphs, which combine various types of imaging data and capture the
associations among individuals within a population. In medical imaging,
population graphs have demonstrated promising results, mostly for
classification tasks. In most cases, the graph structure is pre-defined and
remains static during training. However, extracting population graphs is a
non-trivial task and can significantly impact the performance of Graph Neural
Networks (GNNs), which are sensitive to the graph structure. In this work, we
highlight the importance of a meaningful graph construction and experiment with
different population-graph construction methods and their effect on GNN
performance on brain age estimation. We use the homophily metric and graph
visualizations to gain valuable quantitative and qualitative insights on the
extracted graph structures. For the experimental evaluation, we leverage the UK
Biobank dataset, which offers many imaging and non-imaging phenotypes. Our
results indicate that architectures highly sensitive to the graph structure,
such as Graph Convolutional Network (GCN) and Graph Attention Network (GAT),
struggle with low homophily graphs, while other architectures, such as
GraphSage and Chebyshev, are more robust across different homophily ratios. We
conclude that static graph construction approaches are potentially insufficient
for the task of brain age estimation and make recommendations for alternative
research directions.
Kyriaki-Margarita Bintsi, Tamara T. Mueller, Sophie Starck, Vasileios Baltatzis, Alexander Hammers, Daniel Rueckert
A Comparative Study of Population-Graph Construction Methods and Graph Neural Networks for Brain Age Regression

###### Abstract

The difference between the chronological and biological brain age of a subject can be an important biomarker for neurodegenerative diseases, thus brain age estimation can be crucial in clinical settings. One way to incorporate multimodal information into this estimation is through population graphs, which combine various types of imaging data and capture the associations among individuals within a population. In medical imaging, population graphs have demonstrated promising results, mostly for classification tasks. In most cases, the graph structure is pre-defined and remains static during training. However, extracting population graphs is a non-trivial task and can significantly impact the performance of Graph Neural Networks (GNNs), which are sensitive to the graph structure. In this work, we highlight the importance of a meaningful graph construction and experiment with different population-graph construction methods and their effect on GNN performance on brain age estimation. We use the homophily metric and graph visualizations to gain valuable quantitative and qualitative insights on the extracted graph structures. For the experimental evaluation, we leverage the UK Biobank dataset, which offers many imaging and non-imaging phenotypes. Our results indicate that architectures highly sensitive to the graph structure, such as Graph Convolutional Network (GCN) and Graph Attention Network (GAT), struggle with low homophily graphs, while other architectures, such as GraphSage and Chebyshev, are more robust across different homophily ratios. We conclude that static graph construction approaches are potentially insufficient for the task of brain age estimation and make recommendations for alternative research directions.

Keywords:Brain age regression Population graphs Graph Neural Networks

## 1 Introduction

Alzheimer's disease [9], Parkinson's disease [24], and schizophrenia [18], among other neurodegenerative diseases, cause an atypically accelerated aging process in
the brain. Consequently, the difference between an individual's biological brain age and their chronological age can act as an indicator of deviation from the normal aging trajectory [2], and potentially serve as a significant biomarker for neurodegenerative diseases [7, 13].

Graph Neural Networks (GNNs) have been recently explored for medical tasks, since graphs can provide an inherent way of combining multi-modal information, and have demonstrated improved performance in comparison to graph-agnostic deep learning models [23, 1]. In most cases, the whole population is represented in a graph, namely population-graph, and the structure of the graph is pre-decided and remains static throughout the training. However, since the structure is not given but has to be inferred from the data, there are multiple ways that a graph can be constructed, which will possibly lead to different levels of performance. Every subject of the cohort is a node of the graph, with the imaging information usually allocated as node features. The interactions and relationships among the subjects of the cohort are represented by the edges. Two nodes are chosen to be connected based on a distance measure, or similarity score, usually based on the non-imaging phenotypes. The population graph is used as input to a GNN for the task of node prediction, most commonly node classification.

Brain MRI data have been widely used in population graphs in combination with GNNs. Parisot et al. were the first to propose a way to construct a static population graph using a similarity score that takes into account both imaging and non-imaging data [23]. Many works were based on this work and extended the method to other applications [16, 28]. When it comes to brain age regression, which is an inherently more complicated task than classification, there is little ongoing work using population graphs. To our knowledge, only Stankeviciute et al. [25] worked on the construction of the static population graph using the non-imaging information, but the predictive performance was relatively low.

In all of the cases mentioned above, the graph is chosen based on some criterion before training and remains static, thus it cannot be changed. However, since there are multiple ways that it can be constructed, there is no way to ensure that the final structure will be optimized for the task. This problem has been identified in the literature and multiple metrics have been described in order to evaluate the final graph structure [29, 21], with homophily being the one most commonly used. A graph is considered homophilic, when nodes are connected to nodes with the same class label, and hence similar node features. Else, it is characterized as heterophilic. It has been found that some GNNs, such as Graph Convolutional Network (GCN) [17], and Graph Attention Network (GAT) [26], are very sensitive to the graph structure, and a meaningless graph along with these models can perform worse than a graph-agnostic model. On the other hand, other models, such as GraphSage [14] and Chebyshev [10], are more resilient to the graph structure, and their performance is not affected as much by a heterophilic graph [30].

In this work, we implement and evaluate the performance of different static graph construction methods for the task of brain age regression. We test their
performance on the UK Biobank (UKBB), which offers a variety of both imaging and non-imaging phenotypes. The extracted population graphs are used along with the most popular GNNs, and more specifically GCN [17] and GAT [26], which are architectures highly sensitive to the graph structure, as well as GraphSage [14] and Chebyshev [10], which have been found to be resilient to the population-graph structure. The quantitative results and the visualization of the graphs allow us to draw conclusions and make suggestions regarding the use of static graphs for brain age regression. The code is available on GitHub at: [https://github.com/bintsi/brain-age-population-graphs](https://github.com/bintsi/brain-age-population-graphs).

## 2 Methods

We have a dataset consisting of a set of \(N\) subjects, each described by \(M\) features. This dataset can be represented as \(\mathbf{X}=[\mathbf{x}_{1},\ldots,\mathbf{x}_{N}]\in\mathbb{R}^{N\times M}\), where \(\mathbf{x}_{i}\) represents the feature vector for the \(i\)-th subject. Additionally, each subject has a label denoted by \(\mathbf{y}\in\mathbb{R}^{N}\). Every subject \(i\) is also characterized by a set of \(K\) non-imaging phenotypes \(\mathbf{q}_{i}\in\mathbb{R}^{K}\).

To establish the relationship between subjects, we introduce a population graph, denoted as \(\mathcal{G}=\{\mathcal{V},\mathcal{E}\}\). The graph consists of two components: \(\mathcal{V}\) represents the set of nodes, where each subject corresponds to a unique node, and \(\mathcal{E}\) represents the set of edges that define the connectivity between nodes.

In this paper, we explore four different graph construction approaches, which are described in detail below. In all of the graph construction methods, the node features consist of the imaging features described in the Dataset section 3.1.

#### 2.0.1 No edges

As the baseline, we consider a graph with no edges among the nodes. This is the equivalent of traditional machine learning, where the node features act as the features used as input in the model. The machine learning model in this case is a Multi-Layer Perceptron (MLP). We refer to this approach as _"No edges"_.

#### 2.0.2 Random Graph

With the next approach, we want to explore whether the way of constructing the graph plays an important role on the performance of the GNN. Thus, we build a random Erdos-Renyi graph [11], where a random number of nodes is chosen as neighbors for every node.

#### 2.0.3 Clinical similarity score (Stankeviciunte et al.)

An approach of creating a population graph specifically for the task of brain age regression was proposed by [25], and it does not include the imaging features at all in the extraction of the edges. Instead, the edges are decided only from the non-imaging information. More specifically, the similarity function for two subjects \(i\) and \(j\) is given by:

\[sim(i,j)=\frac{1}{K}\sum_{k=1}^{K}\mathbf{1}[q_{ik}=q_{jk}], \tag{1}\]
where \(q_{ik}\) is the value of the \(k\)-th non-imaging phenotype for the \(k\)-th subject, and \(\mathbf{1}\) is the Kronecker delta function. Intuitively, the Kronecker delta function will only return 1 if the values of a particular non-imaging phenotype of two subjects match. Two nodes \(i\) and \(j\) are connected if \(sim(i,j)\geq\mu\), with \(\mu\) being a similarity threshold decided empirically.

#### 3.1.1 Similarity score (Parisot et al.)

In most of the related works in the literature, the way of creating the adjacency matrix \(W\) was originally suggested by [23] and it is calculated by:

\[W(i,j)=Sim(\mathbf{x}_{i},\mathbf{x}_{j})\sum_{k=1}^{K}\gamma(q_{ik},q_{jk}), \tag{2}\]

where \(Sim(\mathbf{x}_{i},\mathbf{x}_{j})\) is a similarity measure between the node features \((\mathbf{x}_{i},\mathbf{x}_{j})\) of the subjects \(i\) and \(j\), in our case cosine similarity, \(\gamma(\cdot,\cdot)\) is the distance of the non-imaging phenotypes between the nodes, and \(q_{ik}\) is the value of the \(k\)-th non-imaging phenotype for the \(i\)-th subject.

The two terms in Eq. 2 indicate that both imaging, and non-imaging information are taken into account for the extraction of the edges. The second term is similar to the term \(sim(i,j)\) of Eq. 1.

The computation of \(\gamma(\cdot,\cdot)\) is different for continuous and categorical features. For categorical data, \(\gamma(\cdot,\cdot)\) is defined as the Kronecker delta function \(\mathbf{1}\), as before. For continuous data, \(\gamma(\cdot,\cdot)\) is defined as a unit-step function with respect to a specific threshold \(\theta\). Intuitively, this means that the output of the \(\gamma\) function will be 1 in case the continuous phenotypes of the two nodes are similar enough.

#### 3.1.2 kNN graph

The last two graph construction approaches are based on the Nearest Neighbors (NN) algorithm. We connect the edges based on a distance function, in this case, cosine similarity. Each node is connected to its 5 closest neighbors.

In the first approach, we use the neuroimaging information for the node features, and we also use the node features in order to estimate the distances of the nodes. We refer to this approach as _"kNN (imaging)"_.

Similarly to before, in the second approach, we use the cosine similarity of a set of features in order to extract the graph structure. The node features incorporate the imaging information as before. The difference here is that we estimate the cosine similarity of the non-imaging phenotypes of the subjects with the purpose of finding the 5 closest neighbors. We refer to this approach as _"kNN (non-imaging)"_.

Finally, we use all the available phenotypes, i.e. both the imaging, and non-imaging information, to connect the nodes, again using cosine similarity. We refer to this approach as _"kNN (all phenotypes)"_.
## 3 Experiments

### Dataset

For the experiments of the comparative study we use the UK Biobank (UKBB) [3], which not only offers an extensive range of vital organ images, including brain scans, but also contains a diverse collection of non-imaging information such as demographics, biomedical data, lifestyle factors, and cognitive performance measurements. Consequently, it is exceptionally well-suited for brain age estimation tasks that require integrating both imaging and non-imaging information.

To identify the most important factors influencing brain age in the UKBB, we leverage the work of [6] and select 68 neuroimaging phenotypes and 20 non-imaging phenotypes that were found to be the most relevant to brain age in UKBB. The neuroimaging features are obtained directly from the UKBB and include measurements derived from both structural MRI and diffusion-weighted MRI. All phenotypes are standardized to a normalized range between 0 and 1.

The study focuses on individuals aged 47 to 81 years. We include only those subjects who have the necessary phenotypes available, resulting in a group of approximately 6500 subjects. The dataset is split into three parts: 75% for training, 5% for validation, and 20% for testing.

\begin{table}
\begin{tabular}{l c c c c c} \hline \hline
**Graph** & **GCN** & **GraphSAGE** & **GAT** & **Chebyshev** & **Homophily** \\ \cline{2-5}
**Construction** & \multicolumn{5}{c}{**MAE** (years)} \\ \cline{2-5} Random graph & 5.19 & 3.72 & 5.38 & 3.83 & 0.7495 \\ Parisot [23] & 4.21 & 3.74 & 4.35 & 3.77 & 0.7899 \\ Stankeviciütte [25] & 4.61 & 3.73 & 4.90 & 3.77 & 0.7743 \\ kNN (imaging) & **3.89** & 3.77 & 4.09 & 3.75 & 0.8259 \\ kNN (non-imaging) & 4.76 & **3.68** & 4.98 & **3.72** & 0.7796 \\ kNN (all phenotypes) & 3.93 & 3.76 & **4.07** & 3.73 & 0.8191 \\ \cline{2-5}  & \multicolumn{5}{c}{\(R^{2}\) **score**} \\ \cline{2-5} Random graph & 0.26 & 0.59 & 0.2 & 0.54 & \\ Parisot [23] & 0.49 & 0.59 & 0.47 & 0.55 & \\ Stankeviciütte [25] & 0.4 & 0.59 & 0.3 & 0.58 & \\ kNN (imaging) & **0.56** & 0.58 & 0.51 & 0.59 & \\ kNN (non-imaging) & 0.38 & **0.59** & 0.3 & 0.55 & \\ kNN (all phenotypes) & 0.56 & 0.58 & **0.53** & **0.6** & \\ \hline \hline \end{tabular}
\end{table}
Table 1: Performance of the different graph construction methods along with various GNNs on the test set. For every population graph, the homophily ratio is estimated. The baseline MLP (_No edges_) has a performance of MAE=3.73 years and and \(R^{2}\) score of 0.56. Best performance for each GNN model is highlighted in bold.
### Results

The graph construction methods described in the previous section are leveraged by four GNN architectures and evaluated on the task of brain age regression. For every experiment, both spectral and spatial methods, and more specifically, GCN [17], GAT [26], GraphSage [14], and Chebyshev [10], are used. The results for all combinations of graph construction methods and GNN architectures are shown in Table 1.

As a baseline, we use an MLP (_No edges_) that achieves a MAE of 3.73 years and a \(R^{2}\) score of 0.56. For the case of GCN and GAT, we notice that even though all of the graph costruction methods manage to extract a graph that outperforms the random graph, none of them manages to outperform the simple MLP. The methods that used only the imaging information, i.e. the _kNN graph (imaging)_, or all of the phenotypes, such as _Parisot et. al [23]_, extract graphs that when used as input to the GCN or the GAT models, they perform better (MAE of around 4 years) compared to the ones that only use the non-imaging features for the extraction of the graph, such as _kNN graph (non-imaging)_ (MAE of 4.5-5 years). The GraphSAGE and Chebyshev models perform on similar levels regardless of the graph construction methods, with a MAE similar to the MAE given by the MLP. The graph construction method that performs the best though, is the _kNN graph (non-imaging)_ along with the GraphSAGE, with a MAE of 3.68 years and a \(R^{2}\) score of 0.59.

Finally, we estimate the homophily ratio of the population graphs extracted from the different construction methods based on [22] (Table 1). The most homophilic graphs tend to be the ones that use only imaging or a combination of imaging and non-imaging information for the extraction of the edges.

### Visualizations

Visualizing the extracted static graphs, can provide more insights about the performance of the different GNN architectures. Therefore, we visualize the different population graphs colored based on the age of the subjects in Figure 1. The various graph construction methods result in very different graph structures. More specifically, the graph construction methods that take into account the imaging phenotypes, either alone or in a combination with the non-imaging ones, provide graphs that are more meaningful, as subjects of similar ages are closer to each other. On the contrary, methods such as Stankeviciute et al. [25], create graphs that are similar to the random graph, with the neighborhoods not being very informative regarding the age of the subjects.

### Implementation Details

The static graphs extracted are chosen to be sparse for complexity reasons. The hyperparameters of the different construction methods were selected in such a way that all the extracted graphs have about 40000 to 50000 edges. Regarding the model hyperparameters, every GNN contains a graph convolutional layer
consisting of 512 units, followed by a fully connected layer with 128 units, prior to the regression layer. ReLU activation is chosen. The number of layers and their dimensions are determined by conducting a hyperparameter search based on validation performance. During training, the networks are optimized using the AdamW optimizer [19] with a learning rate of 0.001 for 150 epochs, with the best model being saved. For the similarity threshold we use \(\mu=18\) and for the unit-step function threshold, we use \(\theta=0.1\). Both hyperparameters are selected based on validation performance and sparsity requirements. The implementation utilizes PyTorch Geometric [12] and a Titan RTX GPU.

## 4 Discussion and Conclusions

In this work, we implement and evaluate static population graphs that are commonly used in the literature for other medical tasks, for brain age regression. We use the extracted graphs along with a number of popular GNN models, namely GCN, GAT, GraphSAGE, Chebyshev in order to get insights about the behavior of both the extraction methods, as well as the performance of the different models. By visualizing the graphs and estimating their homophily, we can provide further intuition in why the different static graphs do not work as expected

Figure 1: Visualizations of the population graphs extracted from the different graph construction methods. The nodes are colored based on the age of the subject. Colder colors (i.e. blue) indicate subjects with an older age, while warmer colors (i.e. red) correspond to younger subjects.
and we highlight the problem that extracting static graphs from the data is not straightforward and possibly not suitable for brain age regression.

The reported results in Table 1 indicate that the GCN and GAT are highly sensitive to the graph structure, which is in agreement with the relevant literature [30, 20]. It is clear that the graph construction methods that provide graphs of higher homophily, lead to better performance for GCN and GAT compared to the ones that provide a random-like population graph.

On the contrary, GraphSAGE and Chebyshev are not negatively affected by a graph structure with low homophily. This is because GraphSAGE encodes separately the node's encodings and the neighbors' encodings, which in our case are dissimilar, and hence it is affected less by the graph structure. When it comes to Chebyshev, the model is able to aggregate information from k hops in one layer, while the other GNNs achieve this through multiple layers. Being able to get information from higher order neighborhoods allows the model to find more relevant features, which would not be possible in the 1-hop neighborhood as this is highly heterophilic. But even these GNNs, that are more resilient to the graph structure, cannot distinctly outperform the MLP. This behavior is expected as, according to the relevant literature, these models perform similarly to a MLP, and they can outperform it only after a specific homophily threshold [30]. What is also observed is that the _kNN (non-imaging)_ works slightly better than the MLP both for the GraphSAGE and the Chebyshev, probably because the models were able to capture the information added in the edges.

To further explore the behavior of the graph construction methods and the population-graph they produce, we calculate the homophily ratio and we visualize the graphs. The graphs that are more homophilic perform better along with the models that are more sensitive to the graph structure, such as the GCN and the GAT. We note here that all of the homophily ratios are higher than expected compared to the homophily reported in classification tasks, since we would expect that the random graph would have a homophily of 0.5. This is possibly because of the implementation of homophily for regression, as well as due to the imbalanced nature of the dataset. However, the trend is very clear and in agreement with the graph visualizations and the performance of the models.

All in all, the extraction of static population graphs for brain age regression, and in general for medical tasks for which the graph is not given, does not look very promising. In our opinion, there are multiple directions that should be explored. Firstly, one approach could be to learn the edges of the graph along with the training of the GNN, which allows the extraction of an optimized graph for the specific task at-hand. There is some ongoing work on this adaptive graph learning [27, 15, 8, 5], but more focus should be given. In addition, the creation of GNN models for graphs with high heterophily [29, 30], or the exploitation of graph rewiring techniques that could make the existing GNNs work better on heterophilic graphs [4], have proved to be useful for classification tasks. In the case of medical tasks, it might also be beneficial to incorporate the existing medical insights along with the above. It is also important to make the models and the graphs interpretable, as interpretability can be vital in healthcare, and
it is something that is not currently widely explored in heterophilic graphs [29]. Last but not least, introducing metrics that evaluate population graphs is of high importance, since this can help us understand better the structure of the graph. Even though there is some ongoing work when it comes to node classification [20, 29], metrics regarding node regression have only recently started to being explored [22].

#### Acknowledgements

KMB would like to acknowledge funding from the EPSRC Centre for Doctoral Training in Medical Imaging (EP/L015226/1).
