Predicting physiological aging rates from a range of quantitative traits using machine learning.

It is widely thought that individuals age at different rates. A method that measures "physiological age" or physiological aging rate independent of chronological age could therefore help elucidate mechanisms of aging and inform an individual's risk of morbidity and mortality. Here we present machine learning frameworks for inferring individual physiological age from a broad range of biochemical and physiological traits including blood phenotypes (e.g., high-density lipoprotein), cardiovascular functions (e.g., pulse wave velocity) and psychological traits (e.g., neuroticism) as main groups in two population cohorts SardiNIA (~6,100 participants) and InCHIANTI (~1,400 participants). The inferred physiological age was highly correlated with chronological age (R2 > 0.8). We further defined an individual's physiological aging rate (PAR) as the ratio of the predicted physiological age to the chronological age. Notably, PAR was a significant predictor of survival, indicating an effect of aging rate on mortality. Our trait-based PAR was correlated with DNA methylation-based epigenetic aging score (r = 0.6), suggesting that both scores capture a common aging process. PAR was also substantially heritable (h2~0.3), and a subsequent genome-wide association study of PAR identified significant associations with two genetic loci, one of which is implicated in telomerase activity. Our findings support PAR as a proxy for an underlying whole-body aging mechanism. PAR may thus be useful to evaluate the efficacy of treatments that target aging-related deficits and controllable epidemiological factors.

INTRODUCTION

The study of possible intrinsic aging rates is itself venerable and continues to be motivated by Peter Medawar’s pioneering question: does aging occur by a fundamental process independent of—though affected by—overt disease? [1]. To answer this question, the investigation of human aging rates must confront the experimental hurdle presented by our long lifespans. As a result, longitudinal studies of aging in human cohorts typically require significant resources invested over many years. A reproducible framework for estimating physiological age as a measurement of intrinsic age progression would permit the analysis and evaluation of treatments that target aging-related debilitation without a need for expensive long-term studies [2-4]. To be useful, biological aging rate measurements should be relatively stable across time in longitudinal studies [3, 5, 6], account for aging at different biological levels [7-9], and ideally be associated with disease risk and mortality [3, 8]. Early attempts to estimate “biological age” using individual biomarkers of aging (e.g., telomere length, pulse wave velocity, grip strength, etc) produced only modest correlations with chronological age [10-17]. Since there is little variation in human age at natural death [18], true biological age is unlikely to deviate markedly from chronological age—hence the need for a more tightly correlated measure of biological age.

Measurements of biological age based on DNA methylation patterns have provided significantly higher correlation (R2 > 0.8) with chronological age [19, 20] and can also be associated with disease and death [21-27]. However, the role of DNA methylation in the aging process is still largely unclear [28]; it does not account for higher biological levels such as physiology and phenotype, which are presumably affected by other aging-related mechanisms.

Given that aging is a complex process that manifests in virtually all tissues and organs [29-31], it seems likely that composite indices that summarize many molecular and physiological traits to estimate biological age may provide a route to greater understanding [5, 9]. Many efforts have turned to quantitative models for measuring biological age using sets of biomarkers; these models have had success in measuring variations associated with physiological decline, but may require certain assumptions about biological aging or carefully curated sets of traits [32-39]. An increasingly popular approach for robust measurement of biological age is machine learning, which involves building mathematical models on a set of training data to predict a target variable [40]. Machine learning has been used to predict biological age from DNA methylation patterns [19, 20], blood serum markers [41], image data [42-44], and transcriptomic and proteomic signatures [45-47]. Despite these advances, few biological age measurements integrate traits from distinct hierarchical levels (e.g., proteome, metabolome, and phenome) and across the multiple organ systems in which aging is manifest.

Here we use a machine learning approach with a broad range of biochemical and physiological traits including blood phenotypes (e.g., high-density lipoprotein), cardiovascular functions (e.g., pulse wave velocity) and psychological traits (e.g., neuroticism) as main groups from the SardiNIA longitudinal study of aging [48, 49] to estimate human physiological age, a metric for phenotypic and functional age progression [7]. From the estimated physiological age, we derive a physiological aging rate (PAR) as a measure of individual age progression relative to chronological age. This framework is used to determine the most age-responsive, well-correlated, and relevant traits with respect to physiological aging. Independently, we develop a model for estimating physiological age in the InCHIANTI longitudinal study [50, 51], for which we find results similar to those in the SardiNIA study and show positive correlation between PAR and the epigenetic aging rate (EAR) based on DNA methylation. We determine that PAR is a significant predictor of survival and is appreciably heritable. Additionally, a genome-wide association study (GWAS) reveals single-nucleotide polymorphisms (SNPs) that are significantly associated with PAR. These findings may help to elucidate the physiological and genetic underpinnings of the human aging rate.

MATERIALS AND METHODS

Our primary goal was to predict an individual’s physiological aging rate (PAR) based on his/her quantitative traits by using machine learning methods. We investigated several machine learning models, trait selection and dimensionality reduction methods, and data cleaning and sampling strategies to develop an optimized framework. Each machine learning model was trained on the entire trait set using chronological age as the target variable. The predictive frameworks were evaluated using the coefficient of determination (R2) between the predicted age and chronological age. The most promising model was then used to generate physiological ages for individuals in SardiNIA and separately trained to generate physiological ages for individuals in InCHIANTI.

The general pipeline for estimating individual physiological aging rate is listed below (also see Supplementary Figure 1; we discuss the details in each sub-section of Materials and Methods):

The datasets were cleaned to remove missing values. The strategy for data cleaning was heuristically determined to preserve the maximum number of traits and samples (see Supplementary Materials).

Individuals were binned by chronological age (see Supplementary Figure 2 and Supplementary Materials for details).

Samples were split randomly into the training (90%) and testing (10%) sets in each of the chronological age bins.

Traits were selected by the Fisher score computed on the training set and trait values were transformed with linear discriminant analysis. The number of traits selected was predetermined as the number of Fisher-rank-ordered traits corresponding to the maximum average R2 value between predicted and actual age across 10,000 training-testing splits.

Machine learning models were trained to predict chronological age on the training data set and validated using the testing data set.

Steps 3–5 were reiterated 10,000 times to create enough training-testing splits to ensure wide coverage of all participants. The physiological age of an individual was the mean predicted age across these 10,000 training-testing splits.

The physiological aging rate (PAR) was calculated as the ratio between physiological age and chronological age for a given individual.

For the SardiNIA study, after computing PARs, we estimated the heritability of PAR from pedigree information. We also treated PAR as a quantitative trait to perform a genome-wide association study (GWAS) and to identify genetic variants (i.e., single nucleotide polymorphisms, SNPs) significantly associated with PAR. We examined the reproducibility of PAR measurements across time using the follow-up studies from each dataset. We determined significant differences in PAR measurements between deceased and living subjects. For the InCHIANTI study, we compared our PARs with corresponding DNA methylation age measurements [19, 20].

Datasets

SardiNIA longitudinal study

Funded by the National Institute on Aging in 2001, the SardiNIA Project (age range 14.0 to 101.3 years, with a mean of 43.7 years; 57% female) is a longitudinal study of human aging on the island of Sardinia, which is notable for its long-lived population [48, 49]. Several characteristics made this population ideal for modeling physiological age. First, the study cohort is relatively homogenous, which reduced genetic noise. Second, environmental factors are more controlled in an island population, which reduced possible noise due to environmental heterogeneity. Third, the SardiNIA dataset included sequencing information for most of its cohort, making the study suitable for GWAS [52]. The size of sequence coverage and availability of family relationships permitted heritability estimates from pedigree information. The longitudinal nature of the SardiNIA study (i.e., three follow-up waves in addition to the baseline study) allowed us to study individual aging rate dynamics across time [52].

The dataset included over 6100 participants and more than 200 physiological and cognitive traits in total. After cleaning missing entries in a trait-blind, iterative manner, we retained over 4000 subjects and roughly 140 traits with complete data; they were used for model development.

InCHIANTI longitudinal study

The InCHIANTI study was designed to understand the physiological factors that affect mobility and included measurements from many different physiological systems [50, 51]. The studied population spanned two sites: Greve in Chianti (11,709 inhabitants >65 years) and Bagno a Ripoli (Village of Antella, 4704 inhabitants >65 years). The study started in 1998 with 1,463 subjects (age range 21.0 to 102.0 years, with a mean of 69.0 years; 56% female) and included four follow-up waves with the fourth wave retaining 687 subjects. The majority (1,020) of the subjects were elderly (>65 years old). The study included approximately 2,500 trait variants spanning six main physiological subsystems: central nervous system, peripheral nervous system, perceptual system, muscles, bones/joints, and energy production/delivery [51]. The InCHIANTI study provided similar advantages as the SardiNIA study in genetic and environmental homogeneity and its longitudinal nature. In addition, InCHIANTI included data on DNA methylation for many subjects, which were used to calculate DNA methylation age estimates [50].

For the SardiNIA and InCHIANTI cohorts, we chose to apply the machine learning models to the full set of traits of each cohort. The main reason for this choice is that the vast majority of current machine learning methods require that individuals have all trait values available (i.e., no missing data), and consequently, if we were to try to use the same traits for both cohorts, we could use only the smaller list of overlapping traits. Another advantage of the current approach is that the same method can be applied to data from any new study without requiring that it provides exactly the same set of traits.

Machine learning algorithms

We explored both classification and regression models including the random forest classifier (RFC) [53], elastic net (ElNet) regression [54], k-nearest neighbors regressor (kNNR) [55], LASSO regression [56], multiple linear regression, and support vector machine (SVM) [57] (see Supplementary Materials for model specifications). Classification models required bin labels, which were computed as the center chronological age of each bin. Regression algorithms utilized the individual chronological ages but required binning to uniformly sample subjects from the age bins. The machine learning algorithms were adapted from the open-source Scikit-Learn machine learning library [58].

Random forest classifier

The random forest classifier from the Scikit-Learn library [58] was selected for its optimal predictive performance. Our model used 30 estimators, the Gini impurity criterion, minimum sample splits of 2, and minimum sample leaf sizes of 1. These model parameter choices were determined heuristically. The primary advantage afforded by the random forest model is its natural ability to include ordinal and binary features that are present in both SardiNIA and InCHIANTI studies (e.g., yes/no answers to questionnaires, categorized qualitative ratings of high, moderate, and low, etc).

Trait selection and dimensionality reduction

Fisher trait scoring and selection

We utilized the Fisher scoring method for trait selection [59]. For the training data, the trait values and the corresponding chronological age bin labels were used to assign a Fisher score to each trait. The number of Fisher rank-ordered traits that maximized the average R2 between predicted and actual ages over 10,000 training-testing splits was used for model training. After trait selection, the RFC model included 79 traits while the common-trait RFC model included 52 traits in the baseline SardiNIA study. Further dimensionality reduction was performed using linear discriminant analysis (LDA), which slightly improved the predictive accuracy of the models (see Supplementary Materials, Supplementary Figure 3).

Linear discriminant analysis

Linear Discriminant Analysis (LDA) was used for dimensionality reduction before machine learning. LDA is a generalization of the Fisher discriminant and transforms a set of traits into N-1 linear discriminants where N is the number of discrete classes (i.e., number of age bins). We implemented the default “LinearDiscriminantAnalysis” function from the Scikit-Learn library [58].

Statistical metrics

To measure the predictive performance of the different predictive models, we employed several statistical metrics to determine the accuracy of model-predicted physiological age with respect to chronological age.

Pearson correlation ® and coefficient of determination (R2)

The primary measurement of correlation was the Pearson correlation r and the closely related coefficient of determination R2, which represented the proportion of the variance in the target variable that could be explained by the traits used in the model. The R2 value was computed as the square of the Pearson correlation coefficient r, (R2 = r2):

Mean absolute error (MAE)

The mean absolute error represented the magnitude of difference between two continuous variables. The MAE was calculated for N measurements of two variables (x, y) as:

The MAE measured the average difference between chronological age and physiological age in a cohort of N individuals and is represented in units of years.

Measurements of age progression

Physiological aging rate

The physiological aging rate (PAR) represented the characteristic aging rate for an individual from birth to the time of measurement and was computed as the ratio of predicted physiological age to chronological age:

We chose to use the ratio rather than the difference between physiological and chronological age because the ratio accounts for the age of an individual such that the same difference in ages would result in a larger deviation in the PAR of a younger subject than that of an older subject.

Interpretation of PAR

When PAR < 1, the individual’s physiological age is less than their chronological age, they are aging slower than expected for their age cohort; when PAR > 1, the individual’s physiological age is greater than their chronological age, so they are aging faster than expected.

Testing the association between PAR and mortality

We used two methods to study the effect of PAR on survival (or mortality) using the SardiNIA cohort in which we were able to collect information on mortality. We recorded 263 deaths among the 4,415 SardiNIA participants that had their PAR estimated using all traits. The deceased participants had ages from 16.3 to 76.7 years. We first applied the widely used survival analysis model (Cox proportional hazards model) to test whether PAR is a predictor for survival while adjusting for age and sex. The PAR measurement was base-2 logarithm transformed [i.e., log2(PAR)] such that the effect size of PAR on survival could be conveniently interpreted as the hazard ratio with a two-fold change of PAR. Second, we performed a comparative analysis between the 263 deceased participants and the remaining participants on their aging rates. To remove the potential confounding of chronological age, we randomly paired an age-matched living participant to each deceased participant. The differences in the PAR measurements of the pairs were then calculated as ΔPAR = PARdeceased – PARliving, and a corresponding p-value was obtained from a one-sided (ΔPAR > 0) one-sample t-test. The age-matched comparison was performed 100,000 times and ΔPAR and p-values were calculated for each comparison.

Heritability

We used POLY software (http://people.virginia.edu/~wc9c/poly) to estimate the heritability of the physiological aging rate (PAR) based on the known family structure in the SardiNIA cohort.

GWAS

Genome-wide association studies (GWAS) have been a powerful tool to infer genetic associations with diseases and quantitative traits [60-62]. We treated PAR as a quantitative trait and employed GWAS to determine genetic loci significantly associated with PARs for the SardiNIA study.

Trait ranking

We used three different methods to rank the traits in SardiNIA and InCHIANTI. Each method was designed to assess the relative importance of a trait to a different measure of physiological aging. First, traits were scored for their association to the physiological aging rate (PAR) using the p-value from a two-tailed student’s t-test comparing the mean trait values of individuals in the top PAR quartile and that of individuals in the bottom PAR quartile. The fold change of the trait values between the two quartiles was used to quantify the relation between significant traits and the PAR. Secondly, traits were scored for their association with the physiological age using the Pearson correlation r between the trait values and the measured physiological age of individuals. Finally, traits were scored for their contribution to the predictive performance (R2) of the model. Each trait was independently removed from the model, and the average loss in model performance over 500 training-testing splits was calculated with respect to a baseline performance of R2 = 0.858 (SardiNIA) and R2 = 0.702 (InCHIANTI). We refer to the difference in model performance after removing a given trait as the added value (R2) of the trait.

RESULTS

Intuition about physiological aging rate (PAR)

We used machine learning to infer an individual’s PAR (see Materials and Methods for details), the ratio between physiological age and chronological age, from their physiological profile (i.e., set of quantitative trait values) by comparing that profile to the trait values of other participants in the study (“training data”). For example, a 50-year-old man might have a physiological profile (values of body mass index, pulse wave velocity, NEOPIR assertiveness, etc.) that predicts a physiological age of 48 years. As a result, he would have PAR = 48/50 = 0.96 − i.e., aging slower than expected for his chronological age peers.

Comparison of machine learning models

We applied various machine learning models to predict age from the SardiNIA trait data. Each machine learning model was cross-validated for 10,000 different training-testing splits to ensure convergent values for individuals. Physiological age was then computed as the mean predicted age for an individual across all splits for which the individual was in the testing set. The physiological aging rate (PAR) was calculated as the ratio of physiological age to chronological age. The average PAR and median PAR are close to one (SardiNIA mean PAR = 1.04, median PAR = 1.00). Slowed and accelerated aging correspond to PAR < 1 and PAR > 1 respectively.

We used the coefficient of determination (R2) between the actual and predicted ages to evaluate the performance of each model. Higher R2 values indicated better model performance, with most individuals having concordant chronological and physiological age. All machine learning models produced well-correlated physiological age estimates (R2 > 0.8; see Figure 1A). The physiological age and PAR measurements were also highly correlated between any two machine learning models (R2 > 0.9 for physiological ages; R2 > 0.7 for PARs; see Figure 1B). These similarities suggest the presence of a model-agnostic, intrinsic physiological aging signal. The top-performing model, the random forest classifier (RFC) (R2 = 0.86 for the SardiNIA baseline study; see Figure 1C), reached performance saturation with 37 rank-ordered traits when all individuals were involved in training/testing, or with 75 individuals when all traits were involved in training/testing (see Supplementary Figure 3).

Figure 1

Predictive performance of machine learning models. (A) Comparison of the predictive performance measured by the coefficient of determination (R2) and mean absolute error (MAE) for all machine learning models investigated in the study. The random forest classifier (RFC) model was among the top-performing models. (B) Physiological aging rates were highly correlated between different models. Shown is a tileplot of R2 between PARs obtained from different models where darker green corresponds to higher values. (C) Physiological ages predicted by the RFC model were well-correlated with chronological ages of individuals in the SardiNIA study. (D) Physiological aging rate (PAR) of individuals obtained from the RFC model was weakly correlated with chronological age. All figures shown are for the baseline (W1) SardiNIA study. Similar results as in (C) and (D) were observed in follow-up waves of the SardiNIA study, for the elastic net regression model, and in the InCHIANTI study (see Supplementary Materials).

Results were consistent between the baseline study and all follow-up waves of the SardiNIA study and were also consistent between all study visits in the InCHIANTI study (see Supplementary Figure 4 and Supplementary Figure 5 respectively). In the baseline study, the standard deviation in predicted age was 2.6 years for SardiNIA and 3.3 years for InCHIANTI, and standard deviation in predicted age was uncorrelated with subject age for both studies. In SardiNIA, there was no significant difference between the PAR measures of male and female study participants (see Supplementary Figure 6). The distribution of PARs obtained from the RFC model was also generally independent of age (Figure 1D) and highly similar between the SardiNIA and InCHIANTI studies despite notable differences in the chronological age range of the two cohorts (see Supplementary Materials, Supplementary Figure 7). However, the very youngest individuals were distinctly more likely to have PAR > 1, and the very oldest individuals were more likely to have PAR < 1. This trend may be due to several factors: selection bias, with individuals who survived to advanced age being more likely to have a slower intrinsic aging rate; different underlying biology for the youngest individuals who are undergoing development; or artefacts in the classification procedure caused by the absence of data for individuals younger than the youngest or older than the oldest in the cohort (see “Linear rescaling and trimming to correct for the imbalanced distribution of PARs across age” in Supplementary Methods for details on linear rescaling and/or trimming of predicted ages to reduce the PAR imbalance). Removing the youngest and oldest age bins resulted in a more balanced distribution of PARs across chronological age but did not significantly change the correlation of physiological age with chronological age in either SardiNIA or InCHIANTI (see Supplementary Figure 8).

Using common traits

Models for estimating physiological age are especially useful if they can work with commonly measured traits. Therefore, we trained the RFC model on a reduced set of SardiNIA traits that are part of routine clinical screenings. The physiological age obtained from this model was again strongly correlated with chronological age (R2 = 0.78 for RFC). The addition of ten cardiovascular biomarkers that are commonly measured in cardiology specialty clinics resulted in a total of 65 traits (52 after Fisher trait selection, see Methods; full list in Supplementary Materials) and further increased model performance to a level comparable to the full-trait model (R2 = 0.85 for random forest on the baseline SardiNIA data, see Supplementary Materials). The physiological age and PAR measurements obtained from this “common-trait model” correlated well with their counterparts from the full-trait model (R2 = 0.95 for physiological age; R2 = 0.8 for PAR, see Supplementary Figure 9). These results are consistent with a clinically relevant intrinsic process of aging (see Discussion). Results for the common-trait model were also consistent for follow-up waves of SardiNIA (see Supplementary Figure 10).

Physiological aging rate is associated with mortality

Next, we assessed the clinical relevance of the new aging rate measurement. Specifically, we tested whether PAR is associated with mortality using two methods. First, a Cox proportional hazards model showed that independent of chronological age and sex, PAR was a significant predictor for survival (263 deaths, ~51,000 person-years). As shown in Table 1, a two-fold higher aging rate was significantly associated with a 3.59 fold increase of the death hazard (p-value = 0.010). As age and sex were included as covariates in the model, we also observed that independent of PAR, a one-year age increase would increase the hazard by 1.11 fold (p-value <2 × 10−16) and that males have a hazard 1.67 times higher than females (p-value = 6 × 10−5). As a point of comparison, we used the well-known Klemera-Doubal method (KDM) to compute biological age using six traits: C-reactive protein, creatinine, glycated hemoglobin, total cholesterol, urea nitrogen, and systolic blood pressure [37]. Unlike PAR, the KDM biological aging rate was not significantly associated with mortality (p = 0.53, Cox regression with age, sex, log2 (KDM aging rate); 275 deaths, ~26,800 person-years).

Table 1

The effect of physiological aging rate (and age and sex) on survival from a cox proportional hazards model.

	Coefficient	Exp (Coefficient) -- hazard ratio	z-statistic	p-value
Age	0.11	1.11	15.69	<2e-16
Sex	0.51	1.67	4.02	5.9e-5
log2(PAR)	1.28	3.59	2.57	0.010

To provide more evidence for the association between PAR and mortality, we then randomly age-matched deceased individuals with living participants and compared their aging rates. We performed 100,000 randomizations and found that 77.4% of the age-matched comparisons had significant differences (p < 0.05) in the mean PAR values of the deceased and living subjects with a mean ΔPAR of 0.016 (Figure 2A). Nearly all random age-matched comparisons reported ΔPAR > 0 (Figure 2A). In contrast, the randomized control showed, as expected under null, few (5.2%) p-values < 0.05 and a mean ΔPAR = 0.00 (Figure 2B). Furthermore, lifespans of these individuals were negatively correlated with the corresponding PARs in both the full-trait model (r = −0.491, Figure 2C) and common-trait model (r = −0.469, see Supplementary Figure 11). These results all suggest that individuals with lower PARs outlived their counterparts with higher PARs.

Figure 2

Physiological aging rates are associated with mortality. (A) Age-matched mortality analysis: 265 deceased participants were randomly paired with age- Matched ± 0.5 years) living participants in the baseline SardiNIA study. We calculated the difference in the mean PAR measurements of the two groups, ΔPAR = PARdeceased – PARliving and the corresponding p-value from a one-sided, one-sample t-test for ΔPAR > 0.The age-matched grouping was performed 10000 times and ΔPAR and p-values were calculated for each of the 10000 comparisons. 77.4% of the age-matched comparisons produced significantly greater than zero ΔPAR values (p < 0.05) and the mean ΔPAR across all comparisons was 0.016. Nearly all comparisons (>99%) had ΔPAR > 0, which indicated that PARdeceased > PARliving on average. (B) Randomized, age-matched control comparisons produced a 5.2% frequency of significantly greater than zero (p < 0.05) ΔPAR values and the mean ΔPAR was 0.00. Consistent with random assignment, 50.8% of the ΔPAR values were greater than zero. (C) Lifespans for individuals were negatively correlated with PARs (r = −0.491).

Correlation of physiological aging rate (PAR) with epigenetic aging rate (EAR)

DNA methylation (DNAm) ages were measured for individuals in the baseline study and the most recent follow-up study of InCHIANTI with the Horvath model [20]. Using the same method for calculating PAR, we calculated a DNA methylation-based epigenetic aging rate (EAR) as the ratio between the DNAm age and the chronological age for an individual. The mean EAR was correlated with the mean PAR (R2 = 0.36, r = 0.6, Figure 3A) despite the two modalities representing measures of aging from separate biological levels (i.e., methylome for EAR; proteome, metabolome, and phenome for PAR). The positive correlation between EAR and PAR persisted after eliminating the effect of chronological age and sex as confounding variables (see Supplementary Figure 12). Both PAR and EAR measurements were largely uncorrelated with chronological age. In addition to DNA methylation-based aging rates from [20], the PAR was also positively correlated with DNA methylation-based aging rates [19], GrimAge aging rates [24], and PhenoAge aging rates [23] (see Figure 3B). We also computed the KDM biological age [37] using eight markers from InCHIANTI, which were C-reactive protein, creatinine, glycated hemoglobin, albumin, total cholesterol, alkaline phosphatase, urea nitrogen, and systolic blood pressure. The PAR was better correlated than KDM biological aging rate with all DNA methylation-based aging rates except PhenoAge, which suggests that the PAR is a relatively good proxy for DNA methylation-based aging measures.

Figure 3

Correlation between epigenetic aging rate and physiological aging rate. Correlation between PAR and other aging rate measures. (A) Mean physiological aging rates (PARs) obtained from physiological age measurements were correlated (R2 = 0.362, r = 0.601) with the mean epigenetic aging rates (EARs) calculated for the same individuals across the baseline and latest follow-up InCHIANTI studies using the Horvath DNAm age. (B) Pearson correlation values between different pairings of aging rate measures including PAR, Horvath DNAm age, the Hannum DNAm age, GrimAge, rescaled PhenoAge, and the corrected KDM biological age with eight covariates.

Heritability and GWAS of physiological aging rate

We estimated the extent to which PAR is heritable using pedigrees. PARs from the full-trait RFC model observed moderate heritability (h2 = 0.33), and PARs obtained using the common-trait model were similarly heritable (h2 = 0.30). These values were comparable to previous estimates for the heritability of aging and longevity (Karasik et al., 2005). Furthermore, modifications of the PAR measurements such as linear rescaling of physiological ages increased heritability up to h2~0.45 (see Supplementary Materials and Supplementary Figure 13; see Supplementary Figure 14 for linearly rescaled physiological age results). Like typical traits, the physiological aging rate, a proxy for the intrinsic aging process, appears to have a significant genetic component.

The heritability of PAR encouraged us to perform genome-wide association studies (GWAS) in SardiNIA using PAR as a quantitative trait. GWAS identified two genome-wide significant loci [16 significant (p < 5.00E-8) single-nucleotide polymorphisms] associated with PAR. Six SNPs fell in the CFI/GAR1 locus located on chromosome 4 (Figure 4A, most significant CFI/GAR1 SNP with p = 5.55E-10), which has already been linked to aging in earlier studies (see Discussion). CFI/GAR1 also harbored the most significant SNP associated with physiological age acceleration (PAA), another measure of age progression (p = 4.07E-9) (see Supplementary Materials), and was among the top SNPs associated with the PAR from the common-trait model (p = 1.731E-7) and with other machine learning models (e.g., elastic net regression model, p = 2.95E-7, see Supplementary Materials).

Figure 4

Significant genetic loci obtained from genome-wide association study. (A) CFI/GAR1 was significantly associated with the physiological aging rate (PAR). CFI is a complement factor and has been linked to age-related macular degeneration and other age-related disorders. GAR1 is an accessory protein for the active telomerase complex and is an eQTL target of the top CFI/GAR1 SNPs. (B) LINC00202 was significantly associated with the PAR and corresponded to a long non-coding RNA that has been indirectly linked to age-related disease. Plots were made using LocusZoom (Pruim et al., 2010).

Three significant SNPs and seven significant proximally located SNPs (p < 5.00E-8) were observed at the second locus, LINC00202, also referred to as FAM238C (Figure 4B, most significant LINC00202 SNP with p = 3.94E-10). The LINC00202 locus also contained the top SNPs for PAR measurements from the common-trait model (p = 1.100E-8) and from other machine learning models (e.g., elastic net regression, p = 2.153E-10). LINC00202 is a long intergenic non-coding RNA with largely uncharacterized function (see Discussion). Other genes with suggestive p-values for possible associations with PAR (p < 5.00E-7) include APLF, ARHGAP15, SCD5, SOGA2, HIVEP1, ANKRD26, DYNC2LI1, ZNF518B, TMEM45A, CLDN10, and CSMD1 (see Supplementary Materials for a discussion of these genes).

Trait ranking and biomarkers of aging

We ranked the traits in the SardiNIA study and separately for the InCHIANTI study using three distinct methods (see Methods). These included measuring the difference in trait values between the top and bottom PAR quartiles (Figure 5A), correlation of traits with physiological age (Figure 5B), and the additional contribution of each trait to model performance (Figure 5C), which is similar to variable importance for random forest models [63]. Although, as outlined in the Methods, the sets of traits used in SardiNIA and InCHIANTI were not identical, the top ranked traits showed general concordance in the two studies. For example, a set of cardiovascular traits were ranked as top traits in both cohorts. In SardiNIA, CCA intima media thickness, pulse wave velocity, and diastolic CCA diameter were highly ranked across all three methods (Figure 5). Several other cardiovascular traits also received high scores, including peak systolic velocity, integral time velocity, end diastolic velocity, and mean supine systolic blood pressure. Most of these traits were also highly ranked among the set of common and cardiovascular SardiNIA traits (Spearman correlation of ρ = 0.92, see Supplementary Figure 15), corroborating the comparable predictive performance between the common-trait and the full-trait models. The prevalence of high-scoring cardiovascular traits was likewise seen in the InCHIANTI study cohort, including pulse wave velocity, systolic blood pressure, and repolarization phase among the top 5% of ranked traits.

Figure 5

Top traits determined by three independent methods. (A) Volcano plot of the top traits in the full-trait model, which included pulse wave velocity (pwv), CCA intima media thickness (vasIMT), peak systolic velocity (vasPSV), diastolic CCA diameter (vasDiaDiam), waist circumference (exmWaist), and body mass index (exmBMI). Significant differences between the mean trait values of subjects in the top and bottom PAR quartiles were determined using a two-tailed students t-test on each trait. The dotted line corresponds to a Bonferroni-corrected threshold of p = 3.33 × 10−4 calculated from single-test threshold of p = 5.00 × 10−2. Many top traits were also highly ranked in the common-trait model. (B) Traits rank-ordered by Pearson correlation (r) with physiological age measured using the full-trait RFC model. (C) Traits rank-ordered by approximate added value (R2) for the full-trait RFC model in SardiNIA (R2 = 0.858 with all traits used in model). Added value (R2) was averaged over 500 training-testing iterates for each trait. There was significant overlap in the highest ranked traits across all three scoring methods and between the common-trait and full-trait models (see Supplementary Materials). Trait names are as they appear in the SardiNIA study; descriptions of all traits are available in the Supplementary Materials.

Some NEOPIR (Revised NEO Personality Inventory) personality traits [64] were highly ranked in SardiNIA, including fantasy, compliance, excitement seeking, and assistance with personality test. Other highly ranked SardiNIA traits across multiple methods included general measures of obesity such as waist circumference and body mass index, and blood molecular markers such as cholesterol, nitrogen, glucose, sodium, uric acid, and alanine aminotransferase levels. Many of these traits were also present among the highest-ranking InCHIANTI traits (see Supplementary Figure 16: in the top 5%: blood nitrogen, creatinine, sodium intake, hypertension; in the top 20%: cholesterol, IL-6, fibrinogen; and in the top 30%: uric acid, transferrin, waist circumference). In particular, creatinine level was highly ranked by added value (R2) to model performance and was also the highest-ranked trait in the InCHIANTI study. InCHIANTI trait rankings (see Methods) using the epigenetic age and EAR resulted in similar sets of high-ranked features (see Supplementary Figure 17). Most of the top traits from each scoring method have been suggested earlier as biomarkers of aging, and others may become candidate biomarkers (see Supplementary Materials for full rankings of traits).

DISCUSSION

We present machine learning as a promising framework for measuring physiological age from broad-ranging physiological, cognitive, and molecular traits.

Physiological aging rate model

The measured physiological age is well-correlated with chronological age (R2 > 0.8), and different sets of traits give very comparable estimates of physiological aging rate (PAR). Furthermore, PAR calculated from the traits is correlated with the epigenetic aging rate (EAR) (R2 = 0.36, r = 0.6), even though the two aging rates are based on separate biological frameworks. Comparisons with aging rates derived from newer epigenetic age models, such as GrimAge [24] and phenotypic epigenetic age [22], could provide further insight into this relationship.

The physiological age estimates obtained from the model are well-correlated with chronological age and comparable to epigenetic age measures. The approach described here joins other recent attempts to quantify “aging” from functional parameters. It is encouraging that our results correlate with other approaches that are also predictive of risk of mortality. In particular, the PAR, as compared to the KDM biological aging rate, is better correlated with GrimAge, the epigenetic age score with the strongest mortality association [65]. Unlike many earlier studies on aging metrics, we show that PAR is longitudinally stable for an individual (see Supplementary Figure 18), with temporal correlation comparable to that of epigenetic age acceleration indices [66] and higher for shorter time differences, which is a necessary and desirable property for biological age measurements [3, 5, 6]. Also, because PAR predictions are derived from multiple physiological, cognitive, and molecular inputs that span the proteome, metabolome, and phenome, they potentially represent a larger fraction of the manifestations of aging than earlier models. The model maintains high predictive performance even after reducing the trait set to 65 commonly accessible markers. These traits can be efficiently and inexpensively measured in many clinics already equipped to collect them—making the model framework, like some other current models [37], easily adaptable for clinical studies and workflows.

We focused primarily on the random forest classifier model in this investigation, but all machine learning models detected a strong physiological aging signal (R2 > 0.80). The strength of the signal corroborates the presence of a central aging effect that exists independently of model type. This is an idea that has been previously explored in theoretical models [67, 68] and studies where temporal scaling of a central aging rate was observed [69]. The top traits we found in both SardiNIA and InCHIANTI studies are manifest in many tissues and organs. The findings are thus consistent with the likelihood that PAR represents the rate of a pathology-independent aging process—an answer to Medawar’s famous question [1].

Top traits identified through various methods recapitulate previously identified associations between physiological aging and/or age-related disease and cardiovascular function [70-72], molecular markers, physiological measurements, and personality features. A few novel exceptions are promising candidates for new biomarkers of aging (e.g., NEOPIR compliance and agreeableness, step-climbing endurance, and hand tapping coordination).

PAR as a predictor of survival

Improved prediction of chronological age is not necessarily associated with increased prediction of mortality [73]. We used two methods to test directly for association between PAR and mortality/survival. As shown by a Cox proportional hazards model, a two-fold increase of the aging rate would increase the risk of dying by 3.59 fold, or equivalently, a 10% increase of aging rate would increase the risk by 1.19 fold. Our second method showed that in 100,000 randomizations, deceased individuals had significantly higher PARs (p < 0.05) than age-matched living subjects 77.4% of the time. These results clearly show the effect of aging rate on survival and thus a prospective benefit of slowing aging rates. In the same Cox survival model, we also observed the effects of age and sex on survival (e.g., a one-year age increase would increase the risk of dying by 1.11 fold on average, and males have a risk 1.67 times higher than females of the same-age). These results are consistent with numbers from other resources. For example, based on the 2017 Actuarial Life Table published by US Social Security Administration, a 50-year-old man has a death probability 1.09 times higher than a 49-year-old man, and a 50-year-old man has a death probability 1.60 times higher than a 50-year-old woman. These consistent results support the validity of the apparent effect of PAR on survival. In comparison, aging rates computed using the Klemera-Doubal method, another non-epigenetic biological age model, were not significantly associated with mortality.

Heritability and GWAS signals for PAR

The heritability of PAR (~30%) implies a significant genetic contribution to an individual’s aging rate. Genome-wide association studies identified the CFI/GAR1 and LINC00202 loci to be significantly associated with the physiological aging rate. Both loci have also been previously associated with aging and age-related disease.

Complement factor I (CFI) is involved in eliminating pathogens, triggering inflammation, and removing debris from cells and tissues [74-76]. Variants in the CFI gene have been linked to age-related macular degeneration (AMD) [75, 77–80], chronic infection [81], and broader aging processes [77, 82].

But intriguingly, GAR1, encoded by the other candidate gene at the locus, is one of four accessory proteins associated with the telomerase ribonucleoprotein (RNP) complex [83]. GAR1 polymorphisms have been associated with differences in acute heart rate response to exercise, which is a predictor of all-cause mortality and cardiovascular mortality [84]. Other components and effectors of telomerase have also been linked to aging [60, 62, 85, 86].

We retrieved information on significant single-tissue expression quantitative trait loci (eQTLs) for CFI and GAR1 mRNAs from the GTEx database [87]. All significant CFI/GAR1 SNPs corresponding to significant eQTLs were associated with the expression level of GAR1 mRNA, but not with that of CFI. Further, GAR1 mRNA was the most significantly associated (p = 6.5E-13) for the eQTL corresponding to the top CFI/GAR1 SNP (rs11940869) across all tissues. Significant eQTLs corresponding to rs11940869 were found in several tissue types, including brain (frontal cortex, caudate, putamen) and esophagus mucosa (see Supplementary Materials). Furthermore, we analyzed RNA-seq gene expression profiles of tissues from six diverse GTEx tissues (heart, lung, liver, thyroid, colon, cerebellum). GAR1 mRNA levels were reduced with age in heart (r = −0.135), lung (r = −0.313), and colon (r = −0.231) tissues, whereas CFI mRNA expression levels did not correlate substantially with age in most human organs (see Supplementary Figure 19). These results suggest that it may not be coincidental that telomerase, which is associated with GAR1, was determined to be the critical factor that prevents replicative senescence in the fibroblast model for aging: the supply of active telomerase to aging cells promotes their indefinite growth and division, without activation of aging-related SASP (senescence-associated secretory phenotype) proteins. Our data suggest that effects on telomerase activity and possibly related effects on chromatin are not only markers for aging but may be part of the intrinsic mechanism, implicating cell senescence.

As for the other genome-wide significant signal in LINC00202, expression quantitative trait loci affecting LINC00202 expression have been significantly associated with age-related endophenotypes in earlier genome-wide association studies [88]. Although the role of most lncRNAs remains unknown, a number of lncRNAs have been shown to influence molecular processes that contribute to age-related phenotypes [89], and are known to influence age-related cardiovascular disease [90], oxidative stress responses and cellular senescence [91], and many age-related debilitations [92]. However, an exact role of LINC00202 in the aging process remains to be determined.

Possible next steps

We used data from two independent longitudinal studies (SardiNIA and InCHIANTI) with relatively homogenous populations. Whether the significantly PAR-associated loci and highest-scoring traits reported in this investigation are reproducible in more heterogenous populations must be investigated. Further hyperparameter tuning and recent developments in machine learning for handling sparse data [93] may also be of value in maximizing the utility of trait data for age prediction. Although PAR is measured from molecular, physiological, and cognitive traits, inclusion of traits from other data modalities may improve the metric. Several modalities that are aging-associated and potentially useful for PAR measurement include human transcriptome profiles [45, 94, 95], proteome profiles [46, 47], and image data [42, 43]. One natural extension of our method, which is based on cross-sectional trait data, is to use longitudinal data to predict a physiological age at a subsequent visit (e.g., to predict the aging rate at a subsequent visit based on the change of traits over time or based directly on the aging rates predicted at previous visits). This involves new method development and deserves efforts from a separate study.

Although we have emphasized the genetic contribution to aging rates, the majority of PAR is still attributable to environmental factors, including the effect of pathology. Further investigation of correlations between PAR and smoking or age-related diseases such as diabetes, neurodegenerative disease, cancer, and cardiovascular disease could further aid in characterizing the highest-scoring traits and delineating potential aging-related pathways. In turn, these analyses may offer additional clinical utility for PAR in assessing the personalized risk and natural history of patients for each disease.

The efficacy of treatments aimed at slowing the aging process has traditionally been evaluated using individual biomarkers or limited collections of related biomarkers. Our current study has shown that PAR is a significant predictor for survival and correlated with epigenetic aging rate, providing evidence for a good measurement of “aging”. Therefore, we see two potential clinical applications of this intuitive metric—the physiological aging rate: 1) physicians might evaluate the efficacy of treatments and interventions on the aging process by the change of PAR before and after the treatments; 2) based on the values of traits that contribute to the aging rate, physicians can make recommendations to individuals aimed at slowing the aging process (e.g., making specific diet changes, doing specific exercises). Notably, the recommendations would be much the same that are currently recommended by physicians, but the additional link of PAR to risk of mortality might motivate patients to overcome the typical reluctance to modify their regimen.