Cardiac age prediction using graph neural networks

Study Overview

The UKB is a prospective longitudinal study that has recruited approximately 500,000 participants from the United Kingdom with deep genetic and phenotypic data.⁷ A sub-study of UKB invited participants for CMR to assess cardiac volumes and function acquired using a standard protocol.⁸ We first partitioned this dataset by selecting 5064 “healthy” participants and 1330 participants diagnosed with hypertension. Baseline characteristics of both groups are shown in Table 1. We then performed image segmentation and motion tracking to represent the heart as a moving three dimensional mesh. We next built a graph representation of cardiac motion for training a GNN to predict age. The trained model was subsequently used to predict age in the 1330 participants with hypertension. Using the age prediction, we then calculate the “age-delta” in participants with hypertension, quantifying the difference between true chronological age and predicted biological age. An outline of the main steps for the analysis are shown in Figure 1.

• Download figure
• Open in new tab

Figure 1. Study flowchart.

A summary of the main steps in our analysis of cardiac motion meshes to prediction of age-delta and analysis. MRI, Magnetic Resonance Imaging

Image phenotyping

We assessed left ventricular volumes and mass, as well as measures of systolic and diastolic function, using segmentation of 319,650 individual cine images by a fully convolutional network (FCN) (Figures 2 and 3A). In total 104 conventional image-derived phenotypes were collected for each participant (selected parameters in Table 1 and Figure 4), which were subsequently used as tabular data for age prediction in a comparative model using gradient boosting. Compared to healthy controls the hypertensive group had higher body mass index (25.5 kg/m² vs 24.3 kg/m², P = 2.7 × 10^-49), cardiac output (5.6 L/min vs 5.2 L/min, P = 8.5 × 10^-19), and left ventricular mass (86.8 g vs 79.5 g, P = 1.1 × 10^-28), but decreased peak strain rate measurements (peak longitudinal strain rate: 1.5s^-1 vs 1.8s^-1, P = 2.9 × 10^-36; and peak radial strain rate: -5.5s^-1 vs -6.1s^-1, P = 1.3 × 10^-14) (Supplementary Table 1). As inputs to a second comparative deep learning model, we derived time-resolved 3D motion models of the heart using non-rigid image registration between cine images. Aligning data to an atlas template allowed motion traits to be aggregated in a standard coordinate space, which demonstrated complex spatio-temporal differences in motion fields between health states (Figure 3B).

• Download figure
• Open in new tab

Figure 2. Overview of image processing framework.

Short axis cine sequences of the heart from cardiac magnetic resonance imaging (3 sections from the ≈10 acquired slices shown) are segmented and shape refined using a multitask network (red - blood pool, green - myocardium). Mesh fitting of the left ventricle is performed using marching cubes and an atlas template registered to establish anatomical correspondence. A deformation field is obtained from the cine images and motion encoded in the mesh models. The output are left ventricular meshes in 5064 participants with 50 temporal phases where coordinates of each vertex are in a standard reference space. (ED, end diastole; ES, end systole; T, temporal phase 1 - 50.)

• Download figure
• Open in new tab

Figure 3. Segmentation and motion results.

A, Example of deep learning left ventricular segmentation of short-axis cine images from the apical to the basal slices, across 50 temporal phases shown here at three time points. Here endocardial and epicardial surfaces are defined from which volumetric and functional parameters were calculated (value for regional radial strain shown). The motion fields and label maps were shape refined and registered to a cardiac atlas to produce a three dimensional mesh. B, Trajectory of left ventricular contraction and relaxation plotted as looped pathlines for a sub-sample of 100 points on a surface-shaded model of the heart at end systole. An outline of the right ventricle is shown for anatomic reference. The colour of the loops projecting from the surface represent the difference in velocity between healthy (n=5064) vs hypertensive (n=1330) participants. Individual-level dense motion coordinates were used as an input to the cardiac-age prediction networks. ED, end diastole; ES, end systole; LV, left ventricle; RV, right ventricle; T, temporal phase 1 - 50.

• Download figure
• Open in new tab

Figure 4. Population distribution of phenotypes.

Ridge plots summarising the distribution densities for selected anthropometric, haemodymanic and cardiac measurements in healthy (n=5064) and hypertensive (n=1330) participants. LV, left ventricle; EDV, end diastolic volume; ESV, end systolic volume; PC, principal component.

Cardiac motion graph representation

We modelled cardiac motion as a graph representation encoding spatio-temporal information. Each temporal phase is represented by an individual graph derived from the atlas-registered cardiac meshes consisting of a set of nodes with feature vectors capturing local geometry and positional information capturing temporal motion. The graph edges connecting nodes define spatial neighbourhood. The graphs were the input to a GNN trained to predict cardiac age. The GNN was a multi-graph network with convolutional submodels using GCNConv for graph convolutional layers.^9,10 In Figure 5, we present an overview of the GNN architecture, with additional information in Supplementary Materials.

• Download figure
• Open in new tab

Figure 5. Overview of the prediction framework.

High resolution time-resolved meshes of the left ventricle produced from the image analysis framework are first down-sampled using geometric decimation. Each frame is consumed by a graph convolutional submodel and transformed into a feature embedding using global average pooling. Cardiac age is then predicted from the concatenation of time-frame embeddings using a set of dense layers. (ED, end diastole; ES end systole; T, temporal phase 1 - 50; GCN, graph convolutional network.)

Age prediction from cardiac motion

The performance comparison for age prediction using CatBoost, dense neural network and GNN models is shown in Table 3. All three models were trained on a normative cohort with healthy participants (n=4020), and evaluated on a holdout test set of healthy participants (n=1044). The CatBoost model comprised pre-defined imaging parameters and had the lowest performance for age prediction. Improved performance was achieved by a dense neural network which took time-resolved cardiac mesh coordinates as training data which were flattened in the input layer. Both models were outperformed by the GNN model, which was trained on graph representations capturing spatio-temporal relationships of cardiac motion. The GNN yielded a coefficient of determination (R²) of 0.35, a Pearson correlation coefficient (|r|) between predicted age and chronological age of 0.59, (P < 1 × 10^-16) and a mean absolute error (MAE) of 4.74 years (Supplementary Materials).

Effect of training data characteristics

We assessed a range of training sample sizes and degrees of mesh decimation (down-sampling) on the performance of the GNN (Table 4) in the same holdout set used for model comparison in Table 3. For these experiments, 100 hyperparameter searches were performed. A higher level of decimation decreases the performance (compared with the results in Table 3, where 90% decimation rate was used). Whilst reducing training set size significantly reduces predictive performance, a GNN model using only half of the dataset is sufficient to obtain a performance that outperforms CatBoost and dense neural network models using the whole dataset.

Age delta estimation

Age predictions for the healthy holdout set are shown in Figure 6A. Consistent with age prediction studies in other domains,^6,11 a bias-correction step was performed to remove the effects of regression to the mean where age predictions are overestimated for younger participants and underestimated for older participants. Following bias-correction there was no residual correlation between cardiovascular age-delta and chronological age (|r| = 1.38 × 10^-8, P ≈ 1 for hypertensive participants, Figure 6B). Subsequently, we performed a 1:1 propensity matching procedure between healthy and hypertensive groups, matching on age and sex.

• Download figure
• Open in new tab

Figure 6. Predictions by the graph neural network.

A, Predicted cardiac age and chronological age in healthy participants (n=1044). Ages jittered, linear model with 95% confidence intervals. B, Cardiac age-delta and chronological age, in healthy participants (n=1044). C, The gap between chronological age and predicted cardiac age (cardiac age-delta) is increased in subjects with hypertension. The grey shaded density plot represent age- and sex-matched controls with subjects with hypertension. The test samples comprised n=383 healthy males, n=612 hypertensive males; n=670 healthy females, n=718 hypertensive females.

Using linear regression to quantify the effect of hypertension on cardiac age-delta, we showed that hypertension was associated with an increase in predicted age in both males (+1.62 years, P < 1 × 10^-16) and females (+4.19 years, P < 1 × 10^-16) (Figure 6C and Supplementary Materials). We performed an additional linear regression including chronological age as an additional covariate, and found chronological age had a non-statistically significant coefficient (coefficient = 0.051, P = 0.176).

Explainability of graph predictions

The ability to understand the predictions from the graph representations is valuable for physiological interpretation of the learned mapping from cardiac motion to ageing. Here, we employ the GNNExplainer algorithm to explain the predictions of the GNN model.¹² Explanations are represented by node masks assigning an importance weighting to individual mesh nodes where the weights indicate a node’s influence for prediction.

To obtain a population-level explanation, we applied the GNNExplainer algorithm when processing an average temporal cardiac mesh sequence determined over all healthy participants. Extracted node masks at end diastole and end systole are shown in Figure 7. The node masks provide a compact visualization of age-related motion traits showing how these vary both anatomically and spatially in the heart. We also present alternative explainability visualizations based on saliency maps in Supplementary Materials.

• Download figure
• Open in new tab

Figure 7. Explainability node masks.

Visualization of the output of the GNNExplainer algorithm when processing an average temporal cardiac mesh sequence determined over all of healthy participants (n=5064) where A is at end diastole and B is at end systole. Red regions indicate higher node importance levels. An outline of the right ventricle is shown for anatomic reference. The models show how age prediction importance varies from base to apex across temporal frames.

Datasets

The data used in this study were drawn from UKB which comprises approximately 500,000 community-dwelling participants aged 40–69 years who were recruited across the United Kingdom between 2006 and 2010.⁷ All participants provided written informed consent for participation in the study, which was approved by the National Research Ethics Service (11/NW/0382). Our study was conducted under terms of access approval number 40616. The UKB dataset was downloaded from https://biobank.ctsu.ox.ac.uk/. From 2014 a subset of participants was invited for multimodal imaging which included comprehensive CMR imaging performed to a standardised protocol.^8,23 From 39,559 participants with CMR available we identified 5064 healthy participants. In order to define this population, we identified participants with a “good” or “excellent” self-reported level of health (UKB data-field 2178), and who were never-smokers (UKB data-field 22506), not taking regular medications (UKB data-field 20003), and did not have a diagnosis of cardiovascular, chronic respiratory and metabolic disease (e.g. hypercholesterolaemia, diabetes and obesity).⁶ This group, who would be expected to show a pattern of normal ageing, was used to train and evaluate our models.

In order to obtain a sample of participants with hypertension, disease definitions were treated as binary traits, and were defined within the UKB based on hospital records and self-reported data. We considered the participant to have hypertension if their records contained any of the ICD9 codes (4010, 4011, 4019, 4039), ICD10 codes (I100, I110, I119, I120, I129, I130, I131, I132, I139, I150, I51, I152, I159) or self-reported codes (1065, 1072, 1073) relevant to hypertension. Details of both healthy and hypertensive groups are shown in Table 1 and Figure 4.

The age of each participant is computed as the number of years between year of birth (UKB data-field 34) and date of CMR (UKB data-field 53). The sex of each participant is determined using UKB data-field 31.

Image processing

Graph neural network training and evaluation

Graph neural networks (GNNs) are machine learning algorithms specifically designed for the processing of graph-like data structures such as 3D meshes. Here, the cardiac meshes generated through our image pre-processing pipeline, consist of a set of nodes and edges between nodes naturally form a graph structure. Each cardiac time phase yields a graph, and each graph contains the same number of corresponding nodes as a result of the template mesh registration. The graph edges define the spatial neighborhood of nodes, and this information is used within the graph convolutional layers for message passing.²⁹ The graph nodes can carry arbitrary feature vectors capturing relevant information for the prediction task at hand.

For our cardiac motion graph representation, a natural choice of features are the (x, y, z) mesh node coordinates. The relative change in node position captures temporal cardiac motion traits (Figure 3B). Additionally, we use geometric features to capture the local mesh geometry encoding both spatial (within mesh) and temporal (across meshes from time phases) anatomical information. Here, we employ the Fast Point Feature Histogram (FPFH) descriptors, originally proposed for mesh registration and recently shown to be powerful feature descriptors for shape classification.^10,30 To determine the FPFH features on a mesh, first a point feature histogram must be generated. For each query point p_r, all nearby points inside a 3D sphere with radius r and centering on point p_r are chosen (k-neighbor points); then, for every pair p_r and p_k in the k-neighbourhood points of p_r with their normals being estimated as n_r and n_k. The point with the smaller angle between the line joining the pair of points and the estimated normals is chosen to be p_r. Finally, a Darboux frame is defined as (u = n_r, v = (p_k − p_r) × u, w = u × v) and the angular variations of n_r and n_k are computed:

Where the SP F H(p_r) (simple point feature histogram) is obtained by calculating the point features of each neighboring point p_k.³¹ The FPFH feature vectors have length 34, which combined with the 3D positional information in terms of mesh node coordinates, yields node feature vectors of length 37. The input to our GNN model from a single subject is thus a set of 50 graphs with 37-dimensional node feature vectors.

We implemented our GNN model using the PyTorch Geometric Python package. The GNN architecture consisted of 50 submodels (one for each cardiac time phase). Each submodel consisted of two convolutional layers for which we employed the GCNConv operator.⁹ Submodels learned time phase specific feature embeddings using the two graph convolutions followed by a global average pooling across nodes. The resulting embeddings from each time phase were then concatenated and flattened into global feature representations which were processed by a prediction head with two dense layers. We used exponential linear units as activation functions after the first GCNConv layer in each submodel and the first dense layer of the prediction head. The GNN input data was preprocessed using the decimate method of PyVista Python package, with 90% decimation rate. In Figure 5, we present an overview of the GNN architecture, with additional details being presented in the Supplementary Methods.

We randomly split the healthy group into separate training (80%, n = 4020) and test (20%, n = 1044) sets. From the training set, a 10% validation holdout set (n = 403) was used for hyperparameter search, using the Python package Optuna with tree-structured Parzen estimator (TPE) search algorithm,³² where 300 models were trained. Additionally, from the remaining training set, a 10% additional validation holdout set was randomly taken to be used for early stopping evaluation of the GNN also using 300 for training.

In Table 5, we present details of the hyperparameter search space where the learning rate and weight decay hyperparameters are optimized in a log base 10 scale.

Dense neural networks

As baseline for comparison, we trained dense neural networks, using the same decimated base meshes as used in the GNNs. An initial quadratic decimation was done on a template mesh, which as before, generates new points in space. Following this step, we use MIRTK registering and matching to find the nearest points in the original template set. We used the indices of those points to decimate all the participants’ meshes. Additional details of the architecture are in the Supplementary Methods.

The data setup for early stopping and hyperparameter search was identical to the GNN and, similarly, 300 models were trained and the model with the lowest mean squared error on the holdout dataset was chosen. In Table 6, we present details of the hyperparameter search space where the learning rate and weight decay hyperparameters are optimized in a log base 10 scale.

CatBoost

As another baseline for comparison, we trained a gradient boosting machine learning algorithm (CatBoost) on pre-defined tabular data. Gradient boosting algorithms are ensemble learning methods that successively train weak estimators (in the case of CatBoost, symmetrical decision trees) on the residuals of the weighted average of the weak estimators. The process of training extra weak estimators can continue until a fixed number of estimators are trained or until early stopping criteria are reached given that the loss on the training set should always decrease (or at least not increase), but the loss on the validation may eventually start increasing.

CatBoost was trained using 104 conventional CMR parameters as summarised in Table 1. These parameters are measures of left ventricular volume, function, mass and strain in longitudinal and radial directions. These low-dimensional features are implicit in the higher-dimensional time-resolved meshes used in our GNN model. CatBoost was trained using the default hyperparameter set (with automatic choice of learning rate and number of decision trees), except for early stopping rounds (patience) of 100, in this case, we used the reserved 10% hyperparameter search dataset as the evaluation for early stopping criteria.

Cardiac age delta calculation

Tabular data have been used to compute “cardiovascular age-delta” in UKB using image-derived features of vascular and cardiac function.⁶ We follow a conceptually similar approach here using high-dimensional cardiac motion mesh data as inputs to the prediction network. We first predict age in the healthy group, and subtract this from the participants’ chronological age to obtain a cardiac age-delta, subsequently apply a bias-correction procedure,¹¹ (Supplementary Materials) and then use this model to predict age in a group of participants with hypertension. Finally, we fit a linear regression with agedelta as the response variable and presence of disease as covariate. We confirmed there was no residual dependence of cardiac age delta on chronological age using linear regression.³³

Graph neural network interpretability

Saliency maps and the GNNExplainer algorithm were used to explore different explainability technqiues of our results.^12,34 Saliency maps were originally proposed for image classification as a visualisation technique applied to one specific image and class at a time but can be adapted to regression and GNNs by leveraging the model’s ability to assign weights for each node of the input graph. We used the versatile implementation in the Captum Python package. GNNExplainer, on the other hand, was proposed to work directly on graph neural network models. The algorithm computes a node mask via iterative optimization identifying a node’s individual influence (or weight) on the final prediction. While saliency maps correspond to importance of specific graph edges, the GNNExplainer’s ability to highlight important nodes leads to better interpretability in our context where cardiac mesh nodes correspond to specific anatomical locations.