Introduction

Physical inactivity is a worldwide health problem, and is ranked as the fourth leading behavioral risk factor for global mortality1. The imperative to maintain body activity, physically and metabolically, is on the rise. The concept of active aging, inspired by Robert Havighurst’s activity theory2, suggests that maintaining an active lifestyle is crucial for the well-being of older individuals. The thought of active aging emerged and began to develop in the 1990s, placing strong emphasis on the link between activity and health3.This focus became particularly pertinent due to the worldwide aging population, leading to concerns about inactivity permeating various social domains4. Within the transition into the 2020 s, there has been an escalating emphasis on harnessing technology to foster healthy aging5,6,7. Beyond longevity, active aging encompasses regular physical activity, better management of chronic diseases, and improved quality of life8.

Conventionally, several studies have examined various physical aspects of active aging, such as sleep, sedentary time, muscle strengthening activities, and movement behaviors9,10. In this study, we focus on a specific aspect of body activity — physical performance (e.g., muscle endurance and strength)- and explore its potential association with individual’s metabolomics profiling. As the new technology to diagnose diseases, metabolomics involves the comprehensive analysis of small-molecule metabolites (<10 kDa) present in a biological sample, including metabolic intermediates, hormones, signaling molecules, and secondary metabolites11,12. Functioning as the culmination of all biological processes in the body, metabolites play a pivotal role in energy generation, signal transmission, and carrying essential information about the body’s status and ongoing functions. Consequently, important metabolites possess the potential to serve as aging biomarkers or as integral components of the metabolic signature. This signature mirrors the active state of the organism as it traverses the aging process13. The development of metabolomics empowers us to scrutinize health issues at the molecular level14. Notably, amid the COVID-19 pandemic, metabolomics has demonstrated its potency as a diagnostic, prognostic, and drug intervention tool15. As expected, COVID-19 has been extensively investigated using metabolomics methodologies, contributing to biomarker studies16,17 and evaluations of drug impacts18,19. Beyond specific disease diagnosis, metabolomics can also illuminate our comprehension of bodily activities (active aging). Recent endeavors have delved into metabolic profiling within aging studies13,20, providing us with overarching insights into metabolic changes during the aging process. Physical performance (e.g., walking test, muscle strength) is a critical aspect of active aging. In this work, we specifically focus on the metabolomic profiling of older adults and aim to detect key biomarkers and important metabolic interactions linked to physical performance and active aging.

The emerging large scale datasets from OMICS (metabolomics, proteomics, transcriptomic and genomics) measurements empower us to scrutinize any question in biology from a systemic perspective21,22. In the field of systems biology, a central goal is to identify biomarkers and infer biochemical regulations from large-scale metabolomics data23. Statistical methods, especially when combined with machine learning techniques, have shown power on constructing accurate classifiers capable of distinguishing between diverse sample groups and revealing underlying biomarkers24,25,26,27. However, statistical methods offer limited insights into how information is transferred within a biochemical network, the critical regulatory steps involved, and how regulatory mechanisms change under different conditions11,21,23,28. Several studies have emphasized the necessity of analyzing the dynamic behavior of metabolism to understand the evolution and maintenance of stable metabolic homeostasis under varying environmental conditions29,30,31,32,33.

Mathematically, kinetic models can provide systemic insights into metabolic networks, but constructing these models and estimating parameters poses challenges, particularly for large-scale models34. In recent years, several studies have focused on the steady-state Jacobian investigation of metabolomics data35,36,37,38,39 integrating with fluxomics or time-series measurements. In addition, with only large sampled metabolomics measurements, recent studies have developed inverse differential Jacobian algorithms, which provide a convenient way to infer differences in the dynamics of metabolic networks between different conditions21,30,31,32,33,40,41,42,43,44. Among them, the most recent study has developed a novel method and workflow termed COVRECON for analyzing key biochemical regulations through the solution of a differential Jacobian problem21,31,32,33. The COVRECON approach integrates the covariance matrix of metabolomics data with automatic metabolic network modeling based on genome-scale metabolic reconstructions and biochemical reaction databases32.

Figure 1 illustrates the workflow of this study, which consists of three main steps. In step 1, we aim to cluster the original samples into different groups based on physical and functional measurements, where each group represents different body activity conditions calculated from individual’s physical performance measurements. Step 2 involves building machine learning-based classifiers to identify these different groups using metabolomics data, thus enabling the identification of key metabolites as biomarkers. Finally, in step 3, we employ the inverse Jacobian analysis and the COVRECON workflow to uncover the most important biochemical regulations associated with the identified body activity conditions. By conducting this approach, we intend to contribute to the understanding of active aging from a metabolomics perspective and shed light on the key biochemical regulations underlying different fitness conditions.

Fig. 1: Work scheme of the proposed approach.
figure 1

This figure illustrates the key steps of our proposed methodology. In the first stage, we employ the Canonical Correlation Algorithm (CCA) to derive a body activity index that is highly linked to metabolomics data. Subsequently, we perform sample clustering, wherein samples are grouped into 2, 4, and 6 clusters based on this index. To address the non-linear effects of metabolites, we employed a set of automated machine-learning classifiers through repeated double CV approach to predict the groups from metabolomics data, allowing us to identify crucial variables. Finally, we employ the COVRECON workflow32 to construct a topological metabolic interaction model for measured metabolites and to solve the differential Jacobian problem of the system.

Full size image

Results

This study was performed in 5 retirement homes in Vienna managed by Curatorship of Viennese Retirement Homes, with the main aim to assess the impact of resistance training and protein-vitamin supplementation or a cognitive training on physical performance45. The subjects were randomly assigned to three groups (see Figure S1 in Supplementary Data 1): resistance training (RT), resistance training and supplements (RTS) and cognitive training, acting as a control group (CT). Blood samples were collected at the baseline (T1), after three months (T2) and after six months (T3). Oesen et al. have shown resistance training and supplementation can improve the physical performance of the older adults. In this study we analyze all samples from these three groups together, focusing on the relationship between physical performance and the metabolomics profiles of the older adults.

In this secondary analysis, we conducted plasma metabolomics measurement. Plasma plays a crucial role in maintaining body health by transporting nutrients, hormones, and waste products throughout the body. Changes in plasma composition can serve as indicators of various diseases, making it a valuable medium for health monitoring and diagnosis. Numerous studies have demonstrated associations between plasma metabolomic profiles and different health conditions46,47,48. Here we focus on the plasma metabolomics changes and the identification of potential biomarkers and biochemical processes related to physical performance in older adults. Thus, treatment-specific improvements in physical strength are not central to our objective. Instead, we analyze all samples collectively to investigate the relationship between active aging and metabolomics patterns. The cohort of older adults with an average age above life expectancy consisted of 117 participants at baseline and altogether we measured 263 plasma metabolomics samples.

In order to establish the relationship between body physical performance and metabolomics profiles, we initially investigated the physical data measurements, which consisted of two types: body physical performance and body shape. Moreover, the body physical performance measurements can be further divided into resistance exercise and endurance exercise types. Table 1 shows the group differences of the physical measurements across the three groups. As expected, compared to the Cognitive Training (CT), the resistance training groups (RT and RTS) exhibit better resistance measurements. Nevertheless, there was no influence on endurance measurements (e.g. walking distance). Notably, endurance exercise has been reported more related to body aging conditions than resistance measurements49,50. This is also consistent with the experimental design, where the old adults were randomly assigned to the three groups regardless of their fitness.

Table 1 Differences in the physical fitness of the three treatment groups
Full size table

Canonical Correlation (CCA) based clustering to assess physical fitness in a cohort of older institutionalized adults

Since our aim was to investigate the relationship between metabolomics and body physical performance, we first employed Canonical Correlation Analysis (CCA) to generate a body activity index (BAI) based on the physical performance measurements. Subsequently, we clustered the old adults and samples into two, four, and six groups based on this body activity index.

As demonstrated in Fig. 2a, the generated body activity index has a high correlation to the metabolomics index (Pearson Coeff = 0.8471, \(p=1.5* {10}^{-19}\)), where the CCA loadings of the body activity index is listed aside. Among all physical performance indexes, walking distance showed the most dominant effect within the body activity index. This observation is biologically reasonable since walking distance directly reflects an individual’s endurance condition, which is directly related to the aging process49,50.

Fig. 2: Body activity index and metabolomics index from canonical correlation analysis (CCA) and cluster analysis.
figure 2

a The scatter plot of the generated Body Activity Index and Metabolomics Index and loadings of Body Activity Index. b The scatter plots of physical (left) and metabolomics (right) profiles, where all old adults are clustered into two groups based on Body Activity Index. c The scatter plots of physical (left) and metabolomics (right) profiles with all samples clustered into 2, 4, and 6 groups.

Full size image

Considering the potential non-linear relationship between the generated body activity index and the metabolomics index, we constructed an automated machine-learning classifier using the XGBoosting algorithm as described in the method part. The automatic classifier was trained with 50 maximum models, over 30 random datasets separation, the averaged AUCs calculated on the hold-out test sets through a repeated double CV approach were determined to be 91.50%, 82.36%, and 62.17% for the two, four, and six-group clusters, respectively (Fig. 2c). This indicates that the CCA generated body activity index and metabolomics index exhibit a strong correlation. Meanwhile, we group all the old adults into two groups for the inverse Jacobian analysis using the mean body activity index as shown in Fig. 2b.

For comparison, we also performed CCA analysis between the metabolomics data and body shape features, such as gender, height, and age. The biplots of the CCA from the metabolomics and body shape analysis are presented in Fig. S2 in Supplementary Data 1. The highest Pearson’s correlation coefficient obtained was only 0.4963 (\(p=1.5* {10}^{-19}\)) for the age index. Meanwhile, this correlation can partially originate from the physical performance difference generated in this body shape index (e.g., it is evident that age and BMI may influence individual’s physical performance). Additionally, we conducted a further CCA analysis considering the metabolomics data along with both body functionality and body shape data. However, the Pearson correlation coefficient increased only marginally from 0.8471 to 0.8574. This indicates that the metabolomics data are primarily influenced by the body strength/ functionality aspects. Consequently, this validates both the body activity index and the metabolomics index that we developed.

The results of the CCA-based cluster analysis highlight the strong relationship between our derived body activity index and the metabolomics index. The dominance of walking distance as a key factor indicates its significance as a reflection of an individual’s health condition and metabolic activity. In the following analysis, we will focus on the two old adult groups clustered based on the body activity index, labeled as active group and less-active group.

Machine learning based classifiers and variables importance reveals strong association of metabolites and fitness

In this following section, we developed several machine learning based classifiers to predict the active/less-active groups from the metabolomics dataset. This approach can provide us valuable insights on the nonlinear influence between the metabolomics index and body activity index. As described in the Methods section, we evaluated the predictive performance of five machine learning algorithms: XGBoosting (XGB), DRF, GLM, GBM, and DeepLearning (DL) algorithms. For each algorithm, the optimal model was selected using an automated machine learning (AutoML) framework.

Classifier construction followed a repeated double cross-validation scheme based on the area under the receiver operating characteristic curve (AUC). The dataset was divided into four parts. In each iteration, one part was used as the test set, while the remaining three parts were used for model selection via five-fold cross-validation. The selected model was then evaluated on the held-out test set. The final classifier performance was reported as the average of the evaluation metrics across the four folds. This entire double cross-validation procedure was repeated 30 times with different random splits of the dataset.

As shown in Fig. 3, we compared the average performance of the classifiers trained on the full dataset and on three-fourths of the data across the five algorithms. To control potential bias due to random sampling of group labels, we conducted a permutation test with 1500 repeats. All five machine learning models achieved statistically significant results, and all permutation repeats for the five algorithms show negative results, resulting in permutation test p values < 1/1500 ~ 0.00067. The details of permutation test are presented in “Methods” section.

Fig. 3: Different machine learning classifiers results.
figure 3

a Average AUC on 30 hold out test sets through repeated double CV of the five machine learning algorithms for the prediction of two body activity index groups from metabolomics data: XGBoosting algorithm (XGB) (0.915), Distributed Random Forest and Extremely Randomized Trees (DRF) (0.894), Generalized Linear Model with regularization (GLM) (0.903), H2O GBM (GBM) (0.913) and DeepLearning models (DL) (0.910). For each algorithm, we assess the effect of sample size by building a separated classifier through repeated double CV approach with 1/4 training set removed. b The AUC metrics plots of all five algorithms for all the 30 hold out test sets. c The paired t test p value between each two algorithms, it is evident that the two boosting algorithms show statistical significance compared to the other three algorithms. d Comparison of all evaluation statistic metrics for the five algorithms.

Full size image

Figure 3a illustrates the averaged AUC values calculated on the hold-out test datasets through the repeated double CV approach, averaged across 30 random data splits for each algorithm. The detailed AUC results are provided in Supplementary Data 2. Among the five methods, the two boosting methods—XGBoosting (XGB) and Gradient Boosting Machine (GBM)—achieved the highest predictive performance with average AUC values of 0.9150 and 0.913, respectively. As shown in the Fig. 3c, these results are statistically significant (Paired T test P < 0.01) compared to DeepLearning algorithm, which itself outperformed both the Generalized Linear Model (GLM) and Distributed Random Forest (DRF) algorithms (Paired T test P < 0.01). Additional performance metrics comparing these algorithms are also presented in Fig. 3d, with detailed descriptions available in “Methods” section and corresponding statistical values reported in Supplementary Data 2. The superior performance of boosting-based methods suggests the presence of nonlinear patterns in the data, potentially arising from complex interactions within the metabolic network. To evaluate the influence of sample size on classifier performance, we randomly removed one quarter of the training data and reassessed all five algorithms. As expected, AUC values for each method declined slightly, yet remained high overall. This indicates that classifier performance had not yet plateaued and may improve further with larger datasets.

In order to assess the importance of metabolites directly related to the two body activity groups, we ranked the metabolites extracted from the five algorithms based on the testing dataset. We identified the top 10 metabolites for each algorithm by calculating the average variable importance across the 25 repeats. The algorithm-metabolite bipartite graph is shown in Fig. 4a, where Aspartate, Proline, Fructose, Pyruvate and Malic Acid were consistently identified as the top metabolites across almost all classifiers. The detailed metabolite importance values of each algorithm are presented in Supplementary Data 2.

Fig. 4: Variable importance results of the machine learning based classifiers.
figure 4

a Bipartite graph of the top metabolites extracted from the five machine-learning algorithms. For each algorithm, we keep the metabolites if they are identified over 15 times in the top ten metabolites over 30 hold out test sets. b SHAP summary plot of the XGBoosting algorithm in one hold out test set. It shows the contribution of the features for each instance (row of data). c The variable importance for the XGBoosting algorithm in one hold-out test set.

Full size image

For a better understanding of the variable importance results, we applied a multi-algorithm auto-machine learning approach, including all five algorithms with a maximum of 100 models, using the ‘automl’ function in the H2o.py package. XGB achieved the best performance, as shown in Table S1 in Supplementary Data 1. The Pareto front plot in Figure S3 in Supplementary Data 1 determined the optimal subset classifier, which included XGBoosting and GBM classifiers, highlighting the superiority of boosting methods for this task. Figure 4b, c present the variable importance and SHAP summary plot for the leading XGBoosting classifier on the test set. The analysis revealed that Aspartate was the most important metabolite, accounting for over 90% of the importance. This highlights the direct influence of the metabolomics aspect on the body activity index. The Pearson’s correlation heatmap shown in Fig. 5 further supports this observation, with Aspartate exhibiting the most significant correlation with body strength data. Although other metabolites, such as Proline, Malic Acid, and Pyruvate, had lower importance values, they consistently appeared among the top 10 metabolites across different classifiers. In Fig. 5, we also did the t test for all metabolites between the two groups, where the differences with significance are plotted. Interestingly, they didn’t fully match the classifier results, e.g. Pyruvate is identified as key metabolites by all classifiers but didn’t show significance. This may suggest that the effect of Pyruvate is non-linear between the two groups. In addition, as shown in Fig. 4c, the SHAP plot of the classifier top metabolites still shows good separation between two groups, albeit with less pronounced distinctions compared to Aspartate. This further indicates that they play a role in reflecting non-linear metabolic effects on the body activity index.

Fig. 5: Statistic results between the two groups.
figure 5

a Heatmap of the Pearson correlation matrix between the measured physical data and metabolomics data and b metabolites’ concentrations comparison between two groups.

Full size image

We choose the eight most important metabolites: Aspartate, Proline, Fructose, Malic Acid, Pyruvate, Valine, Citrate and Ornithine, and map them to the KEGG pathways as shown in Fig. S4 in Supplementary Data 1. We can see aside from a few large comprehensive pathways, the top metabolites identified in the classifier results are most related to Central carbon metabolism in cancer and 2-Oxocarboxylic acid metabolism. However, it merely revealed a surface-level connection between active aging and these pathways, which falls short of providing a comprehensive understanding of the underlying biochemical regulations of the active aging dynamics.

Predictive inverse metabolic interaction modeling using the COVRECON platform

While the machine learning and classifier results provide insights into the variable importance between the measured metabolites and the body activity index, this does not explain the mechanistic change between the two groups. Since for each old adult, metabolomics analysis was done three times, first time point, after 3 months and after 6 months, we plotted the correlation heatmap of the change of all body features and metabolomics measurement changes within two the time intervals in Fig. 6. It is evident that the correlation patterns within the metabolomics measurement changes show high similarity. This reflects the internal dynamics of the metabolic networks. Nevertheless, when we check the highly correlated metabolites, we may find no biochemical reactions between the two metabolites from any database. This situation frequently happens, e.g. in Fig. 6, Threonine, Tyrosine and Valine show a high correlation, yet no direct biochemical reactions occur among them. This is because the high correlations originate from the network dynamics. Thus, finding the causal interactions among the metabolites is crucial.

Fig. 6: The correlations of metabolomics changes.
figure 6

Heatmaps of the Pearson correlation matrix of all body features and metabolomics measurement changes within two the time intervals (left: first 3 months, right: another 3 months).

Full size image

In recent years, inverse differential Jacobian algorithms have been developed, providing a convenient way to infer causal dynamics of metabolic networks from metabolomics data21,30,31,32,40. Besides the metabolomics measurements, metabolic reconstruction is used as complementary information to build a topological model for metabolic interaction network. Based on this, we have developed the COVRECON toolbox (available at: https://bitbucket.org/mosys-univie/covrecon/src/main/)32. As shown in the method part, we applied the COVRECON workflow to the two group datasets. The COVRECON workflow consists of two steps: building the metabolic interaction network and the inverse Jacobian calculation.

As described in COVRECON32, we used a default setting in the Sim-Network part to generate a metabolic super-pathway network of the measured metabolites. Each edge in the network represents a feasible pathway between two nodes (metabolites) and reflects a non-zero component in the system Jacobian matrix. The default setting assigns a fixed weight of one to each reaction, and the reverse reaction weight is based on the log value of its delta Gibbs free energy. Additionally, a pathway-step limitation of 4 is set. Detailed information about reactions, enzymes and genes of the resulting metabolic interaction network can be found in the Supplementary Data 3. By integrating the covariance of the metabolomics data from both groups and the Jacobian structure matrix, we can perform the inverse Jacobian analysis in the second part of COVRECON toolbox. The COVRECON workflow and toolbox address the ill-conditioned matrix problem associated with the inverse Jacobian approach through a regression loss-based algorithm, significantly improving its stability and feasibility32,33,42. However, given that the inverse Jacobian approach is based on the Jacobian structure and is more reliable in smaller-sized models, we selected a tailored core part of the whole model containing 10–20 metabolites based on the classifier variable importance results as described in method part. The same network reduction strategy as in Sim-Network was employed, with additional indirect connections added to the reduced model. For example, an additional connection from Proline to Aspartate was added to account for the indirect effects through the connections from Proline to Asparagine and from Asparagine to Aspartate (Fig. 7). In Supplementary Data 1, Figure S8 presents 12 typical results in the repeated calculation. All the repeated results are available in Supplementary Data 5. It is evident that even though the local results are different due to the influence from the Jacobian structure information, the Inverse Jacobian approach shows stability on several highlighted metabolic interactions. For example, the interactions Proline->Aspartate, Ornithine->Aspartate, Citrate->Aspartate and Glutamate->2-oxo glutaric acid are high valued in the resulted differential metabolic interaction network of many repeats. To present the overall metabolic interaction importance, we integrated all the 200 local results into the full differential Jacobian (DJ) by calculating the average value of each metabolic interaction within the repeats. The final R* matrix and the differential interaction network are presented in Fig. 7a, b, respectively. In Fig. 7b, we plot only the highlighted metabolic interactions with calculated value (scaled to 0–1) above 0.5. Here we note, the result showed robustness, with similar overall R* using 100, 200 and 500 repeats. Further results are using 200 repeats.

Fig. 7: The overall inverse Jacobian results integrating all 200 local calculation repeats.
figure 7

Panel a is the average inverse Jacobian matrix; b only plotted highlighted metabolic interactions with inverse Jacobian calculation value above 0.5 (scaled to 0–1); c listed the most relevant enzymes, where Aspartate Transaminase showed 11 times out of the 15 highlighted interactions, with the accumulated value 6.2681.

Full size image

Through this COVRECON approach, we are able to find several important perturbed metabolic interactions between the two body activity index clustered groups. The highlighted interactions and the detailed reactions, enzymes and gene information are presented in Supplementary Data 4. These findings provide valuable insights into the regulatory interactions and dynamics of the metabolic network related to Aspartate, further supporting its importance as the dominant biomarker in the classifiers results. As shown in Fig. 7c, several reactions are consistently identified in several highlighted metabolic interactions. Among these, enzyme aspartate transaminase (AST, EC number 2.6.1.1) is identified in 11 out of the 15 highlighted interactions and shown in all the largest valued interactions: Proline->Aspartate, Valine->Aspartate, Citrate->Aspartate and Glutamate->2-oxo glutaric acid. The enzyme Glutamic-Pyruvic Transaminase (ALT, EC number: 2.6.1.2) is also highlighted. Notably, both AST and ALT are important enzymes in amino acid metabolism, and recently there is indication of their involvement in health-related issues of older adults51,52,53. Furthermore, enzyme asparagine synthetase B (EC number: 6.3.5.4) was identified in 8 out of the 15 highlighted interactions. This enzyme is less studied for health issues of elderly peoples. However, asparagine synthetase (ASNS) deficiency was recently discovered as a metabolic disorder of non-essential amino acids54. Moreover, it is evident that most identified enzymes in Fig. 6c belong to enzyme class of transaminases (EC:2.6.1.-). The transaminase enzymes are important in the production of various amino acids, and measuring the concentrations of various transaminases in the blood is important in diagnosing and tracking of many diseases55.

For a further analysis of the enzymes, we conducted routine blood tests measurements of the old adults across the three time points. Four metabolic enzymes were measured: AST, ALT, Gamma-glutamyltransferase (GGT) and Creatine Kinase (CK). The data measurements are presented in Supplementary Data 2. As shown in Fig. 8 and Figure S6 in Supplementary Data 1, we compared the enzyme measurements between the two groups (active/less active). The results suggested significant differences in AST and ALT, while GGT and CK did not exhibit such significant variations. This observation validates the inverse Jacobian results in Fig. 7. Furthermore, we compared the AST and ALT changes within the two 3-months’ time intervals. As demonstrated in Fig. 8, both AST and ALT showed significant changes in the “active group”, while the changes were not significant in the “less active group” during both 3-months intervals. Notably, the changes also exhibited significant differences between the two groups. Specifically, in the “active group”, AST and ALT demonstrated a significant larger decrease during the first 3 months, followed by a significant larger increase in the subsequent 3-months interval. This suggests a larger plasticity of enzymatic liver and muscle systems in individuals with a high level of body activity. Interestingly, a few studies have revealed similar observations while investigating enzyme variations. In a long-term study of 29 routine laboratory measurements of 30 athletes, AST and ALT exhibited significantly larger variations over an 11-months period compared to those reported for general population56,57. Moreover, various studies have evidenced the enzyme fluctuations within healthy individuals’ blood samples from physical activity and exercises10,58,59,60,61,62.

Fig. 8: The enzyme measurements for the two enzymes identified in inverse Jacobian results.
figure 8

a aspartate Transaminase (AST) and b glutamic-Pyruvic Transaminase (ALT). The enzymes are measured in three time points: first measurement, 3 months later and 6 months later. The enzymes’ concentrations are compared between two groups (physical active and less active groups) and different time points, where significances are highlighted.

Full size image

As a comparison to the COVRECON results, we performed disease and pathway enrichment analyses using MetaboAnalyst. The results are presented in Fig. 9. While the pathway enrichment analysis yielded limited informative results, the disease enrichment analysis interestingly highlighted cirrhosis—a condition commonly observed in older adults. Cirrhosis refers to advanced liver scarring resulting from various causes, including hepatitis and alcohol use disorder. Although our dataset does not include clinical information indicating whether participants were diagnosed with cirrhosis, the enrichment of cirrhosis-related metabolites suggests a degree of heterogeneity among participants in terms of liver health. This observation aligns with the findings from our COVRECON network analysis, though MetaboAnalyst provides less detailed insight into specific network-level changes.

Fig. 9: The metabolites enrichment analysis results generated by MetaboAnalyst.
figure 9

a Disease Enrichment results and b Pathway Enrichment results. Disease Enrichment results (a) and pathway enrichment results b generated form MetaboAnalyst.

Full size image

We also conducted a metabolite network analysis using MetaboAnalyst, with the results shown in Figure S9 in Supplementary Data 1. This network was constructed based on a precision matrix, which does not incorporate structural information from biochemical databases such as KEGG. In contrast, COVRECON explicitly focuses on differential networks between conditions and leverages prior biochemical knowledge to highlight detailed and biologically meaningful connections.

Discussion

In this article, we measured 263 plasma metabolomics samples to study active aging and fitness in a cohort of very old adults close to or above the average life expectancy. Using a CCA approach, we clustered all old adults and samples into two groups based on a body activity index. Then we identified several key biomarkers between these two groups through machine- and deep learning analysis. The identified metabolites are Aspartate, Proline, Fructose, Malic Acid, Pyruvate, Valine, Citrate and Ornithine, where Aspartate showed dominant effects. XGboosting showed the best performance. In a further analysis, we applied the COVRECON (Li et al., 2023) approach to the two groups of old adults. Through this method, we identified several key metabolic interaction changes between the active versus the less active group. Many of these interactions are related to aspartate, which is also consistent with the machine learning results. By checking the detailed enzyme information of the highlighted metabolic interactions, we identified several important enzyme regulations. The enzyme AST showed a relation to most highlighted interactions in the COVRECON analysis. Blood measurements of all individuals across the three time points validated the results. Existing studies also showed that AST and ALT are highly related to health issues and dementia of older adults63.

Since the study is conducted at three time points with different treatments, we also analyze the metabolic context for resistance training. As shown in Table S2 in Supplementary Data 1, we conducted a group difference t-test for the metabolomics measurements with alpha Tocopherol showing significant differences between the nutritional supplement intake group (E) and the other two groups, as it is a part of the supplement FortiFit. The metabolites Linoleic acid, Methionine, Palmitic acid, Succinate and Tyrosine showed a significant difference between the control group (K) and the resistance training groups (T & E). Interestingly, this divergence contrasts with the results obtained from the body activity classifiers, suggesting distinct metabolic mechanisms for resistance exercise and endurance exercise. This mechanistic difference between endurance and resistance exercise has been previously explored64, where the metabolite changes induced by endurance or resistance exercise are identified in two different modes.

Moreover, several studies have reported that endurance exercise but not resistance exercise has a high relevance to aging related questions. In a study by Cao Dinh et al., 2019, it is reported that among 100 old women (aged over 65 years) strength endurance training significantly reduced senescence-prone T cells, which is widely recognized as age-related65, while intensive training showed no significant influence. In another study, Weiner et al., 2019 concluded that endurance but not resistance training has anti-aging effects while examining a total of 124 healthy previously inactive individuals49. These studies provide additional support for our body activity index and metabolic network analysis.

Furthermore, we generated scatter plots and calculated correlations between changes in our CCA-derived metabolic index and changes in the body activity index for each individual across two time intervals: from baseline to 3 months, and from 3 to 6 months. The results are presented in Fig. 10. We can see the changes in metabolic index and body activity index also show high correlations in each individual group. Most large changes are in the same direction for metabolic index and body activity index. This further highlights the correspondence between the metabolic index and the body activity index.

Fig. 10: Comparison of metabolic index changes and body activity index changes.
figure 10

Scatter plots and calculated correlations between changes in our CCA-derived metabolic index and changes in the body activity index for each individual treatment group across two time intervals: from baseline to 3 months, and from 3 to 6 months. Correlations were computed after outlier adjustment using Winsorization.

Full size image

As shown in the results section, Aspartate is a dominant blood biomarker for body activity and one of the 22 protein-generic amino acids. It is involved in the malate-aspartate shuttle, which facilitates the transfer of electrons and energy between the cytoplasm and mitochondria, ultimately contributing to the production of ATP and the efficient functioning of cellular energy metabolism66. Thus, it is particularly important in tissues with high-energy demands, such as muscle, liver and the heart. This may account for the larger aspartate metabolism in the “active group”. From this point, several groups have evidenced the effect of aspartate as an important supplement for attenuation of exercise-induced hyperammonemia and an increase in exercise endurance67,68. On the other hand, aspartate is involved in the removal of ammonia from the body through the urea cycle69. Performing exercise can lead to ammonia production as a byproduct of energy metabolism. Aspartate may be used to help detoxify ammonia, potentially altering its levels.

Another result generated from COVRECON approach is that old adults with better body activity index have larger plasticity of enzymatic liver and muscle system. AST and ALT are two of the routine blood test enzymes highly related to individual’s liver but also muscle and heart health70, where elevated levels of AST and ALT enzymes beyond a specified threshold may indicate medical condition like hepatitis, liver disease or myonecrosis. The ratio AST/ALT is a significant sign of liver disease. We plotted the AST/ALT ratio changes over the three time points for the two groups in Figure S7 in Supplementary Data 1. The results showed no significant changes across the time points and groups. This suggests that AST and ALT variations originate from non-disease related factors. Furthermore, investigations have provided evidence that physical exercise and improved fitness levels can also lead to a transient elevation of these enzyme levels within a healthy range for individuals without underlying liver issues58,59,60. This exercise-induced transaminase elevation is a well-documented phenomenon, commonly observed in response to vigorous physical activity. It is essential to recognize that these exercise-related increases in AST and ALT levels are typically temporary and return to baseline levels shortly after physical exertion. This indicates larger AST and ALT variations for individuals with better body functionality/activity, as observed in Fig. 8. This viewpoint is also suggested in a long-term study of 29 routine laboratory measurements of 30 athletes, where AST and ALT exhibited significantly larger variations over an 11-months period compared to those reported for general population56,57.

Furthermore, since natural dementia is a hallmark of aging, we explored its associations with our primary findings—despite physical performance being the central focus of this study. Emerging evidence indicates that in Alzheimer’s disease, brain levels of D-aspartate are dysregulated and neuronal N-acetyl-L-aspartate (NAA) is reduced, reflecting impaired neurotransmission and neuronal integrity71; aspartate serves as an agonist at N-methyl-D-aspartate (NMDA) receptors, and its metabolic imbalance contributes to synaptic dysfunction and cognitive decline72. Peripherally, lower mid-life alanine aminotransferase (ALT) and aspartate aminotransferase (AST) levels are associated with increased long-term dementia risk73, while an elevated AST/ALT ratio correlates with poorer cognitive performance and hippocampal atrophy in older adults74. Collectively, these findings underscore regular physical activity as a potent, modifiable factor in reducing natural dementia risk75. Together, central aspartate metabolism and peripheral liver enzyme alterations underscore a liver–brain axis in dementia pathogenesis, suggesting novel biomarkers and therapeutic targets.

On the other hand, regular physical activity has been robustly linked to a lower risk of natural (age-related) dementia, with meta-analyses demonstrating that higher activity levels reduce all-cause dementia incidence by approximately 28% (HR 0.72, 95% CI 0.65–0.80)75,76. Even minimal exercise—such as walking over 6000 steps per day—has been shown to prevent incident dementia in older adults77. Mechanistic studies highlight that endurance exercise elevates neurogenesis and brain-derived neurotrophic factor (BDNF), protecting against neurodegeneration78. And physical activity programs in dementia patients slow cognitive decline and improve walking quality79.

Taken together, this suggests a tightly linked triad of physical activity, metabolomic alterations, and natural dementia. While our work concentrates on the nexus between physical activity and metabolomics, the relationships among the other two pairs are well studied elsewhere. We acknowledge that, without direct investigation in older adult or patient cohorts, our conclusions regarding dementia remain hypothetical and propose that this compelling topic warrants dedicated future research.

In conclusion, in this study we integrate machine learning statistical analysis and COVRECON inverse Jacobian analysis. In metabolomics analysis, machine learning based statistical methods aid us to find the key metabolites. As for the dynamical analysis, aside from kinetic modeling which needs many parameters fitting processes, we showed the predictive metabolic interaction modeling using the inverse differential Jacobian approach. This approach provides a powerful tool to find important dynamic causal molecular regulations between two conditions. By integrating the machine learning results, we showed a robust approach for the inverse differential Jacobian calculation. By robustly identifying aspartate as a biomarker for active aging in the combined ML/COVRECON approach we provide novel insights into fitness parameter such as the body activity index and their linkage to metabolic processes but also a potential link to dementia. In future we propose to integrate metabolomics studies and COVRECON into a unified study linking physical activity, metabolomic alterations, and cognitive decline; exploring this triadic relationship of natural dementia in more detail. Further developments of COVRECON will include weighing of enzyme-level regulation or reaction kinetics—factors that can introduce more precision into the inferred network structure33. Addressing these gaps will be key objectives in our next phase of research.

Methods

Experimental design

This study was performed in 5 retirement homes in Vienna managed by Curatorship of Viennese Retirement Homes. The aim of this study was to assess the impact of strength training, strength training and protein-vitamin supplement or cognitive training on very old, institutionalized adults. This study was conducted in a randomized, controlled, observer-blind design. The subjects were randomly assigned to three groups: resistance training (RT), resistance training and supplements (RTS) and cognitive training, acting as a control group (CT). The details are presented in Supplementary Data 1. Blood samples were collected at the baseline (T1), after three months (T2) and after six months (T3).

One hundred and seventeen subjects were recruited from five senior residences (Fig. S1 in Supplementary Data 1). The exclusion criteria consisted of physical fitness (Short Physical Performance Battery >4) and mental performance (Mini Mental State Examination ≥23). Moreover, they were free of severe diseases such as diabetic retinopathy, CVDs and regular use of cortisone-containing drugs. Before starting the intervention the health and nutritional status was assessed by specialists in internal medicine and gerontology45. All subjects signed informed consent before inclusion in accordance with the Declaration of Helsinki. The study was approved by the ethics committee of the City of Vienna (EK-11-151-0811) and registered at ClinicalTrials.gov, NCT0177511145.

Subject characteristics

The sex distribution (87.6% women; 12.4% men) among participants was representative for the population living in nursing homes. The mean age of the study population was 82.9 ± 6.0 years for women and 84.9 ± 6.7 years for men. The participants had a BMI of 29.27 kg/m2 ± 5.00 kg/m2 45.

Treatment

Resistance training (RT)

The participants performed resistance training twice a week, supervised by a sport scientist. They were able to use elastic bands, chairs and their own body weights. The session consisted of a 10 min warm-up, 30–40 min of strength training that consisted of ten exercises for the main muscle groups (shoulders, arms legs, back, abdomen and chest) and ended with a 10 minutes cool down (Oesen et al. 2015). The exertion was adjusted to the participants’ individual fitness level by adapting the resistance of the elastic band. 15 repetitions were performed and as soon as the exercise could be easily performed by the subjects, the resistance was increased to perform a more difficult version of the exercise and thus obtain a higher training effect (Oesen, et al., 2015).

Resistance training and supplementation (RTS)

The participants performed the same resistance training as the resistance training group. In addition, they received supplements every day and after each training session. The intake of this supplement was controlled. Nutritional Supplement FortiFit, produced by NUTRICIA GmbH, Vienna, Austria, contained 20.7 g protein (56 energy (En)%, 19.7 g whey protein, 3.0 g leucine, >10 g essential amino acids), 9.3 g carbohydrates (25 En%, 0.8 BE), 3.0 g fat (18 En%), 1.2 g roughage (2 En%), 800IU (20 μg) of vitamin D, 250 mg calcium, vitamins C, E, B6 and B12, folic acid and magnesium (Oesen et al. 2015).

Cognitive training (CT)

The participants performed twice a week memory training and finger dexterity exercises in sitting position. Therefore, minimizing the “bias” being alone and not being a part in group activities (socialization factor) (Oesen et al., 2015). Participants of all groups were instructed to maintain their regular food intake.

Blood plasma metabolite extraction and analysis

Several studies addressed the choice of blood sample, revealing that Heparin plasma produces a smaller side effect in the chromatogram spectrum80,81. Concordant with these findings, Heparin was used as an anticoagulant, while blood plasma was separated from fresh blood samples and kept in –80 °C for further clinical analysis. Metabolite profiles of obtained human plasma samples were measured using a gas chromatograph coupled to mass spectrometer82,83. The samples were thawed on ice for 45 min and were vigorously vortexed for 10 s. The extraction consisted of two steps. First, 100 µl plasma were transferred into 1.5 ml Eppendorf tubes, followed by the addition of 600 µl ice cooled MeOH, immediately vortexed for 10 s and left one ice for 15 min for incubation. In order to remove proteins, the samples were centrifuged at 14,000 × g for 4 min at 4 °C. The supernatant was transferred into new tubes and dried down in a SpeedVac. Afterwards the dried pellets were stored at –20 °C.

The second step consisted of extraction with CHCl3. 300 µl of CHCl3 were added to pellets. The further procedure was a repetition of the first step. The supernatant was transferred into new Eppendorf tubes and dried down in SpeedVac. Metabolite extractions were performed in batches of 30 samples of randomly selected subjects.

Quality control-mix

A quality control (QC) and calibration mixture consisted of specific metabolites, including organic acids, amino acids, mono- and disaccharides and substrates of the TCA cycle. The table of metabolites for the QC-Mix is attached to Supplementary Data 2. A calibration curve was prepared with concentrations of 2 µl, 5 µl, 10 µl, 20 µl, 40 µl, 80 µl and 100 µl.

Internal standard (10 µl Pinitol and 10 µl Sorbitol) were added to each sample and to each QC just the day before GC-MS analysis. Afterwards, the samples were dried in a SpeedVac.

Derivatization

First, addition of 20 µl of 40 mg mL–1 of methoxyamine hydrochloride (MeOX) dissolved in pyridine were added to each sample in order to dissolve MeOX in pyridine appropriately, the solution was vigorously vortexed several times and tube was put into hot water. After that, samples were vortexed until pellets were completely dissolved, followed by agitation at 30 °C for 90 min at 750 rpm with a thermoshaker.

N-Methyl-N-(trimethylsilyl) trifluoroacetamide (MSTFA) flasks of 1 ml content was spiked with 30 µl retention index marker solution of alkanes from C10- C40 in hexane. After addition of 80 µl of prepared MSTFA, samples were incubated at 37 °C for 30 min at 750 rpm, followed by centrifugation at 14,000 × g for 2 min at room temperature (24 °C). Immediately after this step, 70 µl of the supernatant were transferred to GC-vials with micro inserts and closed with crimp caps.

GC- MS analysis

Finally, samples were analysed using GC-MS (LECO Pegasus® 4D GCxGC-TOF-MS, Mönchengladbach, Germany) according to Weckwerth et al. 2004 and Leitner et al., 201782,83. Immediately after derivatization, 1 µl of sample were injected utilizing a split ratio of 1:5. The split/splitless injector was kept at a constant temperature of 230 °C equipped with a single-tapered liner with deactivated wool. The GC-MS consisted of an Agilent 6890 (Agilent Technologies, Glostrup, Denmark) using helium as carrier gas at a flow rate of 1 mL min–1. Gas separation was performed on the HP-5MS column (30 m 3 0.25 mm 3 0.25 mm, Agilent Technologies).

The initial temperature of the GC oven was set to 70 °C isothermal for 1 min, followed by a heating ramp of 9 °C \({\min }^{-1}\) to reach 330 °C and hold for 7 min.

Transfer line temperature was 250 °C, and ion source temperature was set to 200 °C. The MS detector was switched off during the first 260 s. Mass spectra were acquired with an acquisition rate of 20 spectra \({{\rm{s}}}^{-1}\) and were recorded in the range of 40–600 m/z, utilizing a detector voltage of 1550 V and electron impact ionization of 70 eV. The metabolite assessment required an exchange of the liner every 70 injections, thus every 2 batches in a row.

The whole data acquisition was performed within 14 batches. Each batch was measured in the same chronological order. At the beginning and at the end of each batch an alkane mix containing C10-C40 and different concentrations of the QC mix were measured to allow for external calibration and check for instrument performance. In order to be able to estimate the carry over effect, and to keep the machine clean from the most abundant metabolites, blank samples that contained only dried extraction reagents and derivatization solvents were injected each 5 or 7 samples. Each batch consisted of plasma samples from 20-30 subjects and was analyzed within 24–32 h. One pooled sample was measured for each batch, in order to assess instrument stability. At the end of every batch, the same QC was measured again to monitor instrumental performance over time. To minimize systematic bias induced by preparation order, samples were randomly distributed into 14 batches. Additionally after each batch (around 60–70 total injections), the liner was exchanged. Each batch consisted of a representative cross section of total samples and was comparable to the total experimental population. For batch effect analysis we performed a PCA analysis84. There is no pronounced batch effect on every PCA component. The top 8 PCA components biplots are presented in Fig. S10 in Supplementary Data 1.

Metabolite identification, peak integration and alignment

After GC-MS analysis the raw data consisted of ion peaks and were preprocessed using LECO Chroma-TOF. The ion fragmentation spectra were matched to fragmentation spectra in NIST library and inhouse libraries and scored with a match probability, taking into account only metabolites with at least 700 similarity score. Analytes were identified by comparison of ion fragments to a reference library of chemical standards and by calibration based on calibration curves generated with the QC mixtures (see Supplementary Data 2 and the raw data is uploaded in MetaboLights). Alkanes measured at the beginning of each batch provided retention indices that were assigned to all ion peaks. Peak integration and alignment was performed with Chroma-TOF software from LECO. The data is presented in Supplementary Data 2.

Blood test and enzyme determination

Venous blood samples were collected from participants following an overnight fast and processed within two hours. Serum levels of aspartate aminotransferase (AST), alanine aminotransferase (ALT), gamma-glutamyltransferase (GGT), and creatine kinase (CK) were measured using standard automated enzymatic colorimetric methods on a clinical chemistry analyzer (e.g., Roche Cobas 8000 or equivalent). Quality control procedures were performed daily according to manufacturer protocols. All assays were conducted in the hospital’s certified clinical laboratory, and results were expressed in units per liter (U/L). These enzymes serve as routine biomarkers of hepatic and muscular function.

Data processing

The data processing steps involved several procedures. Initially, missing values in the metabolomics measurements were imputed using the K-Nearest Neighbors (KNN) method. Following, normalization was performed to reduce heteroscedasticity and adjust for the offset between high and low intensity features, where the log transformation of each metabolite by centering it around its mean (x̅) and scaling it by its standard deviation (s): \({\hat{x}}_{{ij}}=\left(\frac{{\log }_{2}({x}_{ij})-\overline{{\log }_{2}({x}_{i})}}{s}\right)\)

Data clustering

To identify biomarkers and perform the inverse Jacobian analysis, the samples were firstly clustered into distinct groups. The clustering process comprised the following steps. Firstly, based on the information provided in Supplementary Data 2, it was observed that physical measurements could be categorized into two types: “body-shape” data (e.g., gender and height) and “body-functional” data (e.g., walking distance and left standing time). In order to generate a body activity index that reflects body functionality while minimizing the influence of body-shape differences, Canonical Correlation Analysis (CCA)85 was applied. The loadings of this body activity index are presented in Fig. 2a, where it can be observed that walking distance exhibits the strongest effects. The metabolomics-related body activity index generated through CCA was then used to cluster the samples using the k-means method, grouping them based on this body activity index.

Machine learning based classifiers

While the CCA-based clustering approach analyzes the relationship between the body activity index and the metabolic index as a linear method, it may not fully capture the dynamic nature of the metabolic mechanism, which inherently exhibits predominantly non-linear behavior. To capture this non-linear influence and achieve higher accuracy with the identification of important variables, several machine learning based classifiers were employed within an automated machine learning framework, implemented using the H2o package in Python. The classifiers are build to predict body activity groups from the metabolomics data. The features are all the metabolites measurements, thus, the feature dimension is 35. The methods utilized are as follows:

1, Generalized Linear Models (GLM): GLM implements regularized linear models with stochastic gradient descent (SGD) learning. The model is updated iteratively using a decreasing strength schedule, estimating the loss gradient for each sample at a time. This method offers a baseline for the linear effects.

2, Random Forest Classifier (DRF): A random forest is an ensemble meta-estimator that fits multiple decision tree classifiers on different sub-samples of the dataset, utilizing averaging to improve predictive accuracy and mitigate overfitting.

3–4, Boosting Methods: Boosting is an ensemble meta-algorithm that reduces bias and variance in supervised learning. It integrates a family of machine learning algorithms that convert weak learners to strong ones86. The main variation between many boosting algorithms is their method of weighting training data points and hypotheses. We employed two common boosting methods, LGBMClassifier and XGBClassifier. LGBMClassifier (GBM) is a distributed gradient-boosting framework based on decision tree algorithms, originally developed by Microsoft87, while XGBClassifier (XGB) is an open-source library for regularizing gradient boosting88.

5, Autoencoder + deep learning: Deep learning (DL), also known as deep neural networks, is a powerful machine learning method extensively used in pattern recognition, image processing, and bioinformatics89. Prior to training the model, we employed an autoencoder to pre-train it, using the entire unlabeled data, improving model performance, preventing random weight initialization.

In our approach, each of these machine learning methods was integrated into an automated framework that encompasses hyper-parameter optimization. Hyper-parameter optimization entails the selection of ideal parameter values that govern the learning process, aiming to enhance model performance90. In Supplementary Data 1, the Figure S5 provides an overview of the scope of hyper-parameters associated with each machine learning method.

Repeated double cross validation

To optimize hyperparameters and more comprehensively evaluate model performance, we employed a repeated double cross-validation (rdCV) strategy91,92. This advanced validation technique is particularly well-suited for small datasets, as it helps to optimize model complexity while yielding robust and realistic estimates of predictive performance.

The rdCV procedure comprises three main components:

1, Outer Loop (Test Set Evaluation): The dataset is randomly divided into four segments. In each iteration, one segment is held out as the test set, while the remaining three serve as the calibration set. This outer loop is used to evaluate the model’s predictive performance on unseen data, providing an unbiased estimate of generalization ability.

2, Inner Loop (Model Optimization): Within each calibration set from the outer loop, an additional cross-validation is performed to optimize hyperparameters. Model selection is based on the area under the ROC curve (AUC), ensuring a balance between model complexity and predictive accuracy, and reducing the risk of overfitting.

3, Repetition Loop: The entire double cross-validation process is repeated 30 times with different random data partitions. This repetition improves the stability and reliability of the performance metrics and offers insights into the variability of model complexity and selection across different data splits.

Model performance and permutation test

To assess the performance of each classifier, we employed a comprehensive set of evaluation metrics:

1, AUC (Area Under the ROC Curve): AUC measures a model’s ability to distinguish between positive and negative classes across all possible classification thresholds. An AUC of 1.0 indicates perfect classification, whereas an AUC of 0.5 reflects performance equivalent to random guessing.

2, AUCPR (Area Under the Precision-Recall Curve): AUCPR evaluates the trade-off between precision and recall across various thresholds. It is particularly informative when dealing with imbalanced datasets. AUCPR is computed as the weighted average of precision over recall, where the weights correspond to the probability distribution over thresholds.

3 Gini Index: The Gini index is derived from the AUC and quantifies the discriminatory power of a classifier. Gini=2×AUC − 1. A Gini index of 1 indicates perfect discrimination, while a value of 0 indicates no discriminative ability.

4, F1 Score: The F1 score is the harmonic mean of precision and recall, providing a single metric to evaluate classification effectiveness on the positive class. An F1 score of 1 signifies that both precision and recall are perfect; a lower F1 score indicates a trade-off between false positives and false negatives.

5, Logarithmic Loss (LogLoss): LogLoss measures the accuracy of predicted probabilities rather than hard classification labels. It penalizes both overconfident incorrect predictions and underconfident correct predictions. Lower LogLoss values indicate better probabilistic predictions.

6, Mean Squared Error (MSE): MSE quantifies the average squared difference between predicted probabilities and actual class labels. Although more common in regression, MSE is also informative for probabilistic classifiers.

7, Root Mean Squared Error (RMSE): RMSE is the square root of MSE and provides an interpretable scale for average prediction error, in the same units as the predicted values.

8, Mean Per-Class Error: This metric computes the average misclassification rate across all classes. For binary classification, it reflects the average of the error rates in each class, offering a balanced perspective on performance, especially for imbalanced datasets.

To further validate each classifier, we conducted permutation tests to assess the statistical significance of the model performance92,93. Although high evaluation metrics (e.g., AUC) may suggest strong discriminative power, such values can occasionally be achieved by chance due to favorable random splits of training, validation, and test sets. Therefore, it is important to determine whether the observed performance truly reflects meaningful classification ability.

In a permutation test, the class labels of the samples are randomly shuffled, and the same classifier is rebuilt using the permuted labels. This process is repeated multiple times—in our case, 1500 iterations—to generate a null distribution of the performance metric under the hypothesis of no association between predictors and labels. If a classifier trained on permuted labels outperforms the original model, the observed classification performance is likely due to chance, rendering the original model unreliable. In case of AUC metric, the upper threshold P for the P-value of each algorithm is calculated as \(P=\frac{1+\#({{AUC}}_{{Perm}}\ge {{AUC}}_{{obs}})}{N}\), where \({{AUC}}_{{Perm}}\) refers to the \({AUC}\) values obtained from the permuted datasets, \({{AUC}}_{{obs}}\) is the observed \({AUC}\) from the actual data, and \(N\)=1500 is the number of permutation tests. A similar formulation is applied to other performance metrics, replacing the inequality with the appropriate comparison. In the case where none of the permuted models achieve a better performance than the original (i.e., zero out of 1500 permutations exceed the observed value), the resulting p-value is bounded by: \(P < \frac{1}{1500}\approx 0.00067\).

The machine learning classifiers and permutation test code is available in Supplementary Data 6 in Supplementary Data 1.

Feature importance

Feature importance was estimated using a model-based approach, considering a feature to be important if it significantly contributed to the model’s performance. Here, the ‘varimp’ function within the H2o.py package was utilized to rank the important metabolites of each classifier. The importance value is averaged over the 25 training-test separations, and we choose the top 10 metabolites for each machine-learning method.

Predictive metabolic modeling using an inverse Jacobian approach

Statistical and machine learning methods face inherent limitations in elucidating biochemical network dynamics, identifying critical regulatory steps, and capturing condition-specific regulatory changes24. To address this, inverse differential Jacobian algorithms have recently been developed as a powerful approach to infer dynamic regulation of metabolic networks from metabolomics data21,30,31,32,33,40,41,42,94.

In previous studies, we introduced the COVRECON workflow and Matlab toolbox as the standard inverse Jacobian workflow31,32. COVRECON combines the covariance matrix of metabolomics data with automated network modeling based on genome-scale metabolic reconstructions and biochemical reaction databases.

Consider a metabolic network with n metabolites \(\{{X}_{i}{\}}_{i=1\ldots n}\), modeled by a system of ODEs:

$$\frac{d{\boldsymbol{M}}}{{dt}}={\boldsymbol{F}}\left({\boldsymbol{M}}\right),\,{\boldsymbol{M}}=\left\{{M}_{i}\right\}=\left\{\left|{X}_{i}\right|\right\}$$
(1)

where \({\boldsymbol{M}}\) are the concentrations of the n metabolites, and \({\boldsymbol{F}}=\left\{{f}_{i}\left({\boldsymbol{M}}\right)\right\}\) denotes their reaction rates (e.g., mass action or Michaelis–Menten kinetics).

The steady-state Jacobian matrix is defined as: \({{\boldsymbol{J}}}_{{ij}}={\frac{\partial {f}_{i}}{\partial {M}_{j}}|}_{{steady}}\)

$${\boldsymbol{J}}={\frac{\partial {\boldsymbol{F}}}{\partial {\boldsymbol{M}}}}_{{steady}}=\left[\begin{array}{ccc}\begin{array}{cc}\frac{\partial {f}_{1}}{\partial {M}_{1}} & \frac{\partial {f}_{1}}{\partial {M}_{2}}\\ \frac{\partial {f}_{2}}{\partial {M}_{1}} & \frac{\partial {f}_{2}}{\partial {M}_{2}}\end{array} & \cdots & \begin{array}{c}\frac{\partial {f}_{1}}{\partial {M}_{n}}\\ \frac{\partial {f}_{2}}{\partial {M}_{n}}\end{array}\\ \vdots & \ddots & \vdots \\ \begin{array}{cc}\frac{\partial {f}_{n}}{\partial {M}_{1}} & \frac{\partial {f}_{n}}{\partial {M}_{2}}\end{array} & \cdots & \frac{\partial {f}_{n}}{\partial {M}_{n}}\end{array}\right]steady$$
(2)

It captures first-order interactions among metabolites and encodes dynamic regulatory relationships. Steuer et al. 94 derived the following Lyapunov equation linking the covariance matrix \(C\) of metabolite concentrations to the Jacobian \({\boldsymbol{J}}\):

$${\boldsymbol{J}}* C+C* {{\boldsymbol{J}}}^{T}=-2{\boldsymbol{D}}$$
(3)

where the fluctuation matrix D represents the covariance of noise sources acting on the system.

The differences between two conditions can be quantified by the differential Jacobian \(D{\boldsymbol{J}}\), which is calculated from the Jacobians of the two groups:

$${D{\boldsymbol{J}}}_{{ij}}=\left\{\begin{array}{c}\max \left(\left|\frac{{({{\boldsymbol{J}}}_{{\boldsymbol{d}}})}_{{\boldsymbol{ij}}}}{{{({\boldsymbol{J}}}_{{\boldsymbol{h}}})}_{{\boldsymbol{ij}}}}\right|,\left|\frac{{({{\boldsymbol{J}}}_{{\boldsymbol{h}}})}_{{\boldsymbol{ij}}}}{{{({\boldsymbol{J}}}_{{\boldsymbol{d}}})}_{{\boldsymbol{ij}}}}\right|\right)\\ 1,{\,{{\boldsymbol{if}}({\boldsymbol{J}}}_{{\boldsymbol{h}}})}_{{\boldsymbol{ij}}}={\bf{0}}.\end{array}\right.$$
(4)

To compare two conditions, the differential Jacobian \(D{\boldsymbol{J}}\) is defined elementwise from the corresponding Jacobians. It highlights condition-specific regulatory changes. Solving the inverse problem (inferring \(D{\boldsymbol{J}}\) from metabolomics data) requires both the structural information of \({\boldsymbol{J}}\) and an optimization strategy.

COVRECON addresses this by integrating automated metabolic network assembly with an inverse Jacobian algorithm. The method reformulates the Lyapunov equation as a regression problem, where variation in regression loss is more robust than variation in regression solutions. Based on this property, a regression loss matrix R* is constructed to approximate the relative importance of elements in \(D{\boldsymbol{J}}\). Larger R* values indicate stronger regulatory differences between two conditions. For robustness, R* is computed across multiple random realizations of the fluctuation matrix \(D\) with the final score normalized to [0, 1]32.

In the COVRECON workflow and its Matlab toolbox32, we integrate automated metabolic network reconstruction with inverse differential Jacobian analysis via a regression-loss-based algorithm. The method first assembles a metabolic interaction network that encodes the Jacobian structure, then estimates the differential Jacobian by calculating a regression loss matrix R*. The results are visualized in Matlab figures, where interaction pathways can be interactively examined.

In this approach, the Lyapunov Eq. (3) is reformulated as a set of linear equations:

$$\begin{array}{c}{A}_{h}\,{{\rm{q}}}_{h}={{\rm{b}}}_{h}\\ {A}_{d}\,{{\rm{q}}}_{d}={{\rm{b}}}_{d}\end{array}$$
(5)

Where \(A,{\rm{q}},{\rm{b}}\) are generated from corresponding \(C,J,D\) respectively. Li, et al. 32 verified that under numerical variations in \({\rm{b}}\) the variation of the regression solution \({\rm{q}}\) is much larger compared to the variation in the regression loss \(r\). Based on this property, we construct a “regression loss matrix” \({{\boldsymbol{R}}}^{* }\) to capture the relative importance of Jacobian elements rather than directly estimating their absolute values. Specifically, for each element \({J}_{{ij}}\)

$${q}_{s}^{* }={{({A}_{c}}^{T}{A}_{c})}^{-1}{{A}_{c}}^{T}{b}_{s}{R}_{{ij}}^{* }={\min }_{{b}_{s}}{||}{b}_{s}-{A}_{c}{q}_{s}^{* }{||}$$
(6)

Where \({A}_{c}\) is calculated by combining \({A}_{h}\) and \({A}_{d}\) in Eq. (5) with additional constraint that only that single element \({J}_{{ij}}\) is the same between the Jacobians and \({b}_{c}=[{b}_{h};{b}_{d}]\). If this assumption fails, an additional regression loss arises, reflected in larger \({R}_{{ij}}^{* }\)32.

Because only the structure (not the values) of \({D}_{h}\) and \({D}_{d}\) is known, multiple realizations of D are sampled according to its nonzero structure. In practice, 1000 samples are used, and the final \({R}_{{ij}}^{* }\) is taken as the minimum loss across all samples. Larger \({R}_{{ij}}^{* }\) indicates stronger regulatory differences between conditions. For interpretability, \({R}^{* }{\rm{is\; normalized\; to}}\) [0,1]32.

Unlike correlation-based methods, which often include indirect or spurious associations due to the absence of biochemical priors, COVRECON incorporates curated reaction databases to reconstruct direct interactions. This enables the identification of potentially causal mechanisms underlying metabolic regulation. Through this framework, we aim to reveal key components and regulatory interactions embedded in the differential Jacobian.

Integrate classifier biomarkers and group differential Jacobian analysis

Since we have clustered the samples into two groups in the data clustering part, we are now able to do the inverse Jacobian analysis for the two groups. As discussed in Supplementary Data 1, similar to the general approach of most kinetic models, we consider the dynamics within each group is simulated in a group model, thus the steady state dynamics can be represented as a group Jacobian. Consequently, the inverse Jacobian algorithm can offer valuable information of the regulated dynamics between the two groups.

The results from the inverse Jacobian analysis are closely linked to the structural information of the Jacobian obtained from the automatically generated super-pathway metabolic interaction networks. It is essential to highlight that we combine the significance of classifier variables in the context of inverse Jacobian analysis. Simply put, we retain the pivotal biomarkers and introduce a controlled mix of randomly chosen additional metabolites. The augmented networks, encompassing 10–20 metabolites, are subsequently subjected to the COVRECON workflow. Notably, in COVRECON results, large values serve as indicators of the dynamics difference between the two distinct groups. We are able to identify the important reactions or enzymes involved in the active aging context by checking the detailed information behind these large values32.