System‐level metabolic modeling facilitates unveiling metabolic signature in exceptional longevity

Abstract Although it is well known that metabolic control plays a crucial role in regulating the health span and life span of various organisms, little is known for the systems metabolic profile of centenarians, the paradigm of human healthy aging and longevity. Meanwhile, how to well characterize the system‐level metabolic states in an organism of interest remains to be a major challenge in systems metabolism research. To address this challenge and better understand the metabolic mechanisms of healthy aging, we developed a method of genome‐wide precision metabolic modeling (GPMM) which is able to quantitatively integrate transcriptome, proteome and kinetome data in predictive modeling of metabolic networks. Benchmarking analysis showed that GPMM successfully characterized metabolic reprogramming in the NCI‐60 cancer cell lines; it dramatically improved the performance of the modeling with an R 2 of 0.86 between the predicted and experimental measurements over the performance of existing methods. Using this approach, we examined the metabolic networks of a Chinese centenarian cohort and identified the elevated fatty acid oxidation (FAO) as the most significant metabolic feature in these long‐lived individuals. Evidence from serum metabolomics supports this observation. Given that FAO declines with normal aging and is impaired in many age‐related diseases, our study suggests that the elevated FAO has potential to be a novel signature of healthy aging of humans.


| INTRODUC TI ON
Population aging is an increasingly urgent issue confronting many countries worldwide (Chang et al., 2019). As most disabilities and fatal human diseases are age-related (DALYs & Collaborators, 2018), understanding the mechanisms of aging will help with the development of therapeutics for aging-related diseases. Metabolic control plays a crucial role in regulating the health span and life span of various organisms, for example, worms (Leiser et al., 2015) and primates (Mattison et al., 2012). Dysregulated metabolism often leads to premature aging and certain diseases in humans (López-Otín et al., 2016). In contrast, long-lived people, such as centenarians (CENs), may have "healthy" metabolic profiles that support them in resisting age-and metabolic-related diseases, although the exact mechanism remains elusive (López-Otín et al., 2016).
Human metabolism is a complex network that contains thousands of reactions and metabolites, and systematic identification of metabolic changes in health and diseases remains challenging (Brunk et al., 2018). Constraint-based reconstruction and analysis (COBRA) is based on flux balance analysis theory (Orth et al., 2010) and uses different types of constraints, including metabolite availability, nutrient limits, and the most widely available data, gene expression from either microarray or RNA-seq, to build tissue-specific metabolic models (Blazier & Papin, 2012;O'Brien et al., 2015;Vlassis et al., 2014;Wang et al., 2012) and perform metabolic modeling Le Novere, 2015;Mardinoglu et al., 2014aMardinoglu et al., ,2014bYizhak et al., 2015). Many COBRA methods have been developed to perform the metabolic modeling, and most of them were merged into the COBRA toolbox, a desktop software suite of interoperable COBRA methods (Becker et al., 2007;Heirendt et al., 2019;Schellenberger et al., 2011). COBRA methods have been widely used for modeling cellular metabolism (Bintener et al., 2020;Heirendt et al., 2019;Nam et al., 2012), and discovering disease mechanisms (Lewis et al., 2010;Mardinoglu et al., 2018), targets (Larsson et al., 2020;, and drug candidates Bintener et al., 2020).
A major challenge in previous metabolic modeling studies is the "low accuracy" in predicting metabolic fluxes in human cells, largely due to the fact that they considered merely qualitative data, rather than quantitative information. In the most commonly used methods (Agren et al., 2012;Blazier & Papin, 2012;Pacheco et al., 2019;Shlomi et al., 2008;Vlassis et al., 2014;Wang et al., 2012), quantitative gene expression or proteomics data need to be translated into qualitative values (O'Brien et al., 2015). Such kind of doing inevitably leads to the loss of most of the quantitative information and thus introduces biases in predicting metabolic fluxes in human cells.
In this work, we present a systems biology approach to quantitatively integrate omics (i.e., transcriptome and proteome) data and kinetome information into genome-wide precision metabolic modeling (GPMM), aiming to accurately identify metabolic changes in human health and diseases. To benchmark its performance, we applied GPMM and other methods commonly used for metabolic modeling on the same transcriptome data from the NCI-60 cell lines (Reinhold et al., 2012) to compare the predicted metabolic fluxes with the experimentally measured values. GPMM robustly and reliably predicted the experimentally measured fluxes and significantly outperformed the existing methods. We then applied GPMM to study the metabolism of a Chinese centenarian cohort to understand why CENs can delay or avoid many serious age-related diseases. We found that elevated fatty acid oxidation (FAO) is the most significant metabolic feature in the CENs. Further serum metabolomic data showed that the decreased serum fatty acid concentration was the most significant feature in the CENs, supporting our observations from metabolic modeling results. Our study suggested a new signature in exceptional longevity.

| Developing genome-wide precision metabolic modeling method
In the present study, we developed a genome-wide precision metabolic modeling (GPMM) method to address the "low accuracy" challenge. The method quantitatively integrates the enzyme kinetics information from knowledge bases and the enzyme levels from transcriptome and proteome data into metabolic models ( Figure 1 and Figures S1-S3). Specifically, we first curated the generic human metabolic model (Recon 3D;Brunk et al., 2018) (Table S1) and set the upper bounds for the main nutrient uptake rates in blood using information from the literature. To reduce noise from reactions without enzyme kinetics information (the turnover number, Kcats), we then constructed a reduced Recon 3D model to maximize the number of reactions with Kcats and minimize the number of reactions without Kcats (Table S1). We next quantitatively integrated the enzyme kinetic parameters and gene expression levels to constrain the upper and lower bounds of each reaction (Figures S1-S3). According to Michaelis-Menten kinetics, the upper bound of a reaction is the product of the concentration and turnover number (Kcat) of its enzyme. Finally, we used flux variability analysis (FVA) to obtain individual models by removing the reactions with zero flux and then performed Markov Chain Monte Carlo (MCMC) sampling to identify metabolic changes and key regulators.
Notably, GPMM enabled several in silico analyses that were not included in the state-of-art COBRA methods toolbox (COBRA toolbox v3.0), and hence can be broadly applicable in metabolic Research Project, Grant/Award Number: 202101AS070058, 202101AT070299 and 2019FB094 K E Y W O R D S aging, GPMM, longevity, metabolic modeling, omics integration, systems biology engineering and therapy ( Figure 2a and Table S2). These include not only genome-wide single and combinatorial knock-in and knock-out analysis, but also quantitative inhibition and activation analysis to identify key regulators for target discovery ( Figure 2a and Table S2). In addition, GPMM can be applied to conduct personalized metabolic modeling for precision medicine ( Figure 2a and Table S2).

| GPMM robustly and precisely captured the experimentally measured fluxes
Since the input transcriptome may carry out noise from the experimental or mapping procedures, we analyzed the robustness of GPMM to demonstrate whether GPMM has the ability to tolerate gene expression noise. We first constructed a noise-induced transcriptome by adding random values (viz. artificial noise) into the original expression data (viz. the genuine expression values) of each gene. Specifically, we induced 1%, 5%, and 10% noise into the gene expression data of NCI-60 cell lines to construct noise-induced transcriptomes. Then, we performed the metabolic modeling on the genuine and the noise-induced expression datasets, and then compared the obtained flux results between both datasets. If a method is robust, this method should be able to tolerate a certain extent of noise on quantified gene expression; then, the correlation (measured by R 2 or R-squared) between the genuine and noise-induced samples should approach 1. The results showed that the R-squared in each cell line is larger than 0.98 under either 1%, 5%, or 10% noise ( Figure 2b). To further investigate the robustness of GPMM, we next performed multiple random sampling to investigate whether GPMM is still robust at different sampling times (Figure 2c). We induced 5% gene expression noise to the H460 (one of NCI-60 cell lines) for 100 times and performed metabolic modeling using GPMM. We obtained an average R 2 of 0.984 (Figure 2c) between the genuine and the noise-induced samples. Some important fluxes in cancer cells, such as ATP production and lactate secretion, were also consistent among these 100 simulations (Figure 2d). These results indicated that GPMM is a robust method. F I G U R E 1 Flowchart of genome-wide precision metabolic modeling. A generic human metabolism model (Recon 3D) was first curated from the literature, and transcriptome data were then used to estimate enzyme abundance using a steady-state mathematical model. Next, a reducing model was constructed, and the upper bound of each reaction was calculated using the product of the concentration and turnover number (Kcat) of its enzyme. Finally, flux variability analysis (FVA) was performed to reconstruct individual models, and Markov Chain Monte Carlo (MCMC) sampling was used to detect metabolic differences and key regulators To benchmark the performance of GPMM, we chose to utilize the transcriptome and metabolic flux (uptake and secretion) data of the NCI-60 cell lines (Table S3; Jain et al., 2012). Four other methods commonly used for metabolic modeling, GIMME (Blazier & Papin, 2012), Fastcore (Vlassis et al., 2014), rFASTCORMICS (Pacheco et al., 2019), and ecModel (Robinson et al., 2020;Sánchez et al., 2017), were chosen as comparisons. We applied each of the four methods to the transcriptome data of the NCI-60 cells and compared the computationally calculated metabolic fluxes with the reported experimental measurements. The results showed that the GPMM method had a much higher correlation between the predicted and experimental values ( Figure S4a, R 2 = 0.72, p = 2.3e-106) than GIMME ( Figure   S4b, R 2 = 0.011), Fastcore ( Figure S4c The Warburg effect, indicated by an increase in lactate secretion, is one of the most important cancer hallmarks in NCI-60 cells (Jain et al., 2012). We thus compared the predicted lactate secretion fluxes with the experimental values, and found that GPMM Comparisons between predicted metabolic fluxes and experimentally measured lactate fluxes in NCI-60 cells using GPMM (e), GIMME (f), Fastcore (g), rFASTCORMICS (h) and ecModel (i). GPMM, Fastcore, and rFASTCORMICS had R 2 values of 0.86, 0.088 and 0.33, respectively, whereas GIMME failed to predict lactate secretion. ecModel has the ability to predict lactate secretion, but the magnitude of the predicted fluxes is different from the experimental values. Note: the ecModel reconstruction and the flux detection were derived from Zenodo (https://doi.org/10.5281/zenodo.3577466), and only 11 ecModels are available well-predicted the secretion of lactate in NCI-60 cells (R 2 = 0.86, p = 2.2e-24) ( Figure 2e). For other four methods, neither GIMME nor Fastcore could predict lactate secretion (Figure 2f,g). One of the upto-date methods, rFASTCORMICS, returned reasonable results with the R 2 of 0.33 ( Figure 2h). However, the other up-to-date method, ecModel, fails to quantitatively predict lactate secretion (Figure 2i), although its overall prediction performance is reasonable ( Figure   S4e). These results showed that GPMM can precisely capture the experimentally measured fluxes and significantly outperformed the existing methods.

| Metabolic modeling reveals elevated fatty acid oxidation as the most significant metabolic feature in centenarians
To shed light on the metabolic characteristics of the CENs to better understand why these individuals are able to delay or avoid many serious age-related diseases that afflict the normal population (Evert et al., 2003), we aimed to study the metabolism of longevity in a CEN cohort sampled from Hainan Province, China. The cohort included 76 CENs, 54 centenarian-children (F1s), and 41 spouses of centenarian-children (F1SPs ; Table 1), whose RNA-sequencing data were reported in our previous study (Xiao et al., 2018).
We next applied GPMM to study the metabolic features of longevity in this cohort. In total, we developed 171 individual GPMM the CENs and younger controls (viz. F1SPs) and adjusted for age and gender effect, we obtained 343 upregulated and 90 downregulated fluxes. We observed that the overall CEN flux signature was slightly negatively correlated with the aging effect (r = −0.15, p = 7.1e-12) ( Figure S5a,b), suggesting that the CENs contain some signatures that are different from the ones associated with age.
The most striking signature in all four metabolic processes consistently indicated that long-chain fatty acid beta-oxidation ( Consistent with these observations, total cellular ATP production capacity was also significantly enhanced (p = 0.032) in the CENs ( Figure S6b).

| Serum metabolomics supports the metabolic modeling observations
Because a higher systemic FAO leads to higher uptake and con- are likely heritable. Intriguingly, among the upregulated metabolites, the most significant ones were bile acids, a group of metabolites for fatty acid absorption (Figure 4c). These results suggest that the decreased serum fatty acid concentration was the most significant feature in our centenarian metabolomics data. This observation also explains our previous epidemiological survey, which found that total cholesterol is decreased in the CENs compared with F1SPs (He et al., 2014). A similar result was also obtained by analyzing the clinical data of the same longevity cohort studied here ( Figure S7).

| DISCUSS ION
Identifying metabolic signatures in centenarians is important for healthy aging. Constraint-based reconstruction and analysis (COBRA) is a promising method that can capture metabolic signatures in health and diseases, but has been a long-standing challenge to quantitatively predict molecular phenotypes (Lewis et al., 2010;O'Brien et al., 2015). In this study, we present a solution to this criti-  Figure S4). As most parameters and datasets are precalculated, GPMM is easy to use, and the only required input is the transcriptome (RNA-sequencing or microarray). Therefore, GPMM is able to systematically analyze cellular metabolic profiles using only the transcriptome and enables broad computational studies on discovering disease mechanisms.
Previous studies indicated that the existing continued methods, such as PRIME (Yizhak et al., 2014) and RegrEX (Robaina Estévez & Nikoloski, 2015), are less robust than the methods that utilize discretization workflows, such as rFASTCORMICS (Pacheco et al., 2016(Pacheco et al., , 2019. However, the results showed that our developed quantitative method can robustly and well predict experimentally measured fluxes. The reason may be because we not only translated the transcripts to proteome data using a simple but efficient model, but also used enzymatic parameters to restrain the maximum rate for each reaction. In addition, to avoid the effect of the unconstrained reactions on the metabolic simulation, we also reduced the generic models (i.e., Recon3D) to maximize the number of reactions with Kcats and minimize the number of reactions without Kcats. Similar to rFASTCORMICS, GPMM can also use the secretion information to improve the model performance by adding the secretion reaction information in the exchange input file. These improvements thus largely overcome the performance and robustness issues of the existing methods. Therefore, the dramatic improvement of GPMM will enable many computational studies to discover biomedical mechanisms.
Utilizing this method, we studied the metabolic profiles of CENs and identified the elevated fatty acid oxidation as the most significant metabolic feature in the CENs. As the input of GPMM is the transcriptome, we investigated the main beta-oxidation genes from Recon 3D (Brunk et al., 2018) and found that 11 of 14 (78.6%) FAO-related genes were slightly upregulated in the CENs ( Figure   S8a), including four essential peroxisomal beta-oxidation genes (EHHADH, HSD17B4, ACAA1, and ACOX1) and five key mitochondrial beta-oxidation genes (HADHB, ACAA2, ECHS1, ACADL, and ACADVL) ( Figure S8a). We also observed that 9 of 14 (64.3%) FAO-related genes were slightly downregulated with aging in F1SP samples ( Figure S8b). These results support our findings that the elevated fatty acid oxidation as a metabolic signature in the CENs.  Figure   S9b). Consistently, serum metabolism results also showed that most of down-regulated metabolites in serum are FALs (12 of 20, 60%); and most of FALs with significant differences between the two groups are downregulated (12 of 19, 63%) in F1s ( Figure S9c).
Interestingly, we also observed that F1s has lower free fatty acid levels (i.e., trans-vaccenic and trans-vaccenic levels) than F1SP ( Figure 4e,f). Given that F1s might have a higher probability of long lifespan than F1SPs, these results added further support to our conclusion that the elevated fatty acid oxidation is a signature involved to healthy human aging and longevity.
There are many studies in both model organisms and humans showing associations between FAO decline and aging (Gong et al., 2017;Levadoux et al., 2001;Short et al., 2005). However, whether FAO also declines in the CENs, the paradigm of healthy human aging, remains unclear. Interestingly, multiple lines of evidence from our study argues for an enhanced FAO in the CENs, which well explains the previous observations that decreased fatty acid levels, especially the PEs, were found in centenarians' serum (Pradas et al., 2019). In addition, impaired FAO is frequently observed in many age-related diseases, including atherosclerosis (Freigang et al., 2013) and diabetes (Wei et al., 2016). Importantly, elevated FAO is reported to be causally associated with metformin-induced longevity in Caenorhabditis elegans (Pryor et al., 2019). Elevated fatty acid beta-oxidation related genes extend the lifespan of worms . Collectively, these results suggest that the elevated long-chain FAO function in the CENs, at least in female CENs, represents a "healthy" metabolic profile of longevity, which may convey survival advantages to long-lived individuals by reducing lipid accumulation and lowering the risks of common age-related diseases, especially those involved in lipid metabolic disorders.
In summary, we have developed a novel systems biology approach to effectively integrate omics data in the modeling of metabolic mechanisms in human health and disease. This approach dramatically improved the performance over the existing methods.
Our method thus immediately enables many computational studies on discovering disease mechanisms and candidate drug targets, as well as further developments of the algorithms. We applied this method to investigate the metabolic profiles of CENs, and suggested the enhanced fatty acid oxidation as a novel metabolic signature of healthy aging in exceptional longevity.
Nevertheless, there are some limitations that need to be over-

| Knowledgebase curations
To reconstruct genome-wide metabolic models (GEMs) by the GPMM approach, we collected several relevant knowledge bases.
Next, we manually curated the global human metabolic network of Recon 3D using thermodynamic analysis  and the precured Recon 2 model . Because adenosine monophosphate (AMP) cannot be directly changed into adenosine triphosphate (ATP) in any reaction, some reversible reactions, such as FACOAL150 and RE1514M, were curated as irreversible.
The curated Recon 3D model is shown in Table S1.

| Setting quantitative upper and lower bounds of biochemical reactions
For each biochemical reaction, the flux of any reaction has the following equation: where V is the flux of a reaction, V max is the maximum reaction rate according to Michaelis-Menten kinetics, [E] is the enzyme concentration, and Kcat is the turnover number of the enzyme. The Kcat values of human enzymes were obtained from the BRENDA database (Placzek et al., 2017). If an enzyme had multiple Kcat records, their median was used. Where experimental data were missing, we used Kcats from other species.
We obtained 2602 Kcat records in the 4352 reactions with an EC number in Recon 3D (Table S1). Although 42% reactions with an EC number lacked Kcat records, the enzyme abundance percentage was smaller than 10% ( Figure S1).
A previously published method, named GIMME (Blazier & Papin, 2012), was used to reduce the Recon 3D model to maximize the number of reactions with Kcats and minimize the number of reactions without Kcats. The GIMME objective functions were set to ATP production and biomass reaction, as described in previous studies (Blazier & Papin, 2012;Nam et al., 2014

| Predicting enzyme abundance using gene expression data
Recently, a simple but efficient mathematical model was proposed to predict protein abundance using gene expression data (Wilhelm et al., 2014). Changes in enzyme abundance can be determined by the number of proteins synthesized from mRNA minus the number of proteins degraded. In the steady-state, we have following equation: where E and M are the enzyme and corresponding mRNA abundances, respectively; α is the enzyme synthesis rate from mRNA; and γ is the enzyme degradation rate.
Thus, in the steady-state, we can predict enzyme abundance as follows: To estimate the ∕ ratio, we downloaded microarray data of 12 normal tissues with GSE7307 (http://www.ncbi.nlm.nih.gov/ geo/query/ acc.cgi?acc=GSE7307) and RNA-seq data of 15 normal tissues from the Human Protein Atlas Dataset (Uhlen et al., 2015).
We also obtained the corresponding protein abundance data from the MOPED database (Montague et al., 2014) with the unit of nmol/L. MOPED uses the human body map dataset and estimates protein concentration from protein abundance (Montague et al., 2014).
Thus, the ∕ ratio was estimated using the median ratio of protein/ mRNA across multiple tissues.
The correlation between the transcriptome and proteome is usually quite low (Pearson correlation of 0.4-0.5, Figure S2). Notably, using the steady-state kinetic method, the Pearson correlation between the predicted proteome and experimental measurements reached 0.8-0.9 ( Figure S3), indicating that protein abundance could be correctly estimated using transcriptome data.

| Steps of the metabolic modeling
First, enzyme abundance was predicted using the above men-

| Predicting fluxes of NCI-60 cells using GPMM
The gene expression dataset (RNA-seq) was downloaded from the CellMiner website (https://disco ver.nci.nih.gov/cellm iner/loadD ownlo ad.do), and the uptake rate for each cell was obtained from the above experimental fluxes (Jain et al., 2012). Using the GPMM toolbox mentioned above, and setting ATP production and biomass reaction as the optimized functions, genome-wide precision metabolic modeling of NCI-60 cells was performed to obtain the flux matrix.
The predicted secretion fluxes were then compared with the experimental dataset to evaluate the overall performance of the GPMM.

| Predicting fluxes of NCI-60 cells using GIMME
The quantitative gene expression of NCI-60 was the first transformed into qualitative present/absent logical values using a Fragments Per Kilobase of transcript per Million mapped reads (FPKM) cutoff = 3.0. Metabolic models were then reconstructed using the "GIMME" function in the Cobra toolbox 3.0 (Heirendt et al., 2019). Next, MCMC sampling was conducted to obtain the distribution of fluxes, and the average flux for each reaction was then calculated to obtain the GIMME-based flux matrix.

| Predicting fluxes of NCI-60 cells using Fastcore
The consistent Recon 3D model was first constructed using the "fastcc" function in the Fastcore toolbox (Vlassis et al., 2014) with an epsilon of 1e-4 using the linear solver of cplex (https://www.ibm. com/analy tics/cplex -optim izer). The metabolic models were then reconstructed using the "fastcore" function in the Fastcore toolbox, and MCMC sampling was conducted using the "ACHRSampler" function in the Cobra toolbox 3.0 (Heirendt et al., 2019).

| Predicting fluxes of NCI-60 cells using rFASTCORMICS
rFASTCORMICS is an updated version of Fastcore that uses discretization workflows instead of the heuristic thresholds method (Pacheco et al., 2019). The consistent Recon 3D model was also first constructed using the "fastcc" function in the Fastcore toolbox (Vlassis et al., 2014) with an epsilon of 1e-4 using the linear solver of cplex (https://www.ibm.com/analy tics/cplex -optim izer). The metabolic models were then reconstructed using rFASTCORMICS (Pacheco et al., 2019), and MCMC sampling was conducted using the "ACHRSampler" function in the Cobra toolbox 3.0 (Heirendt et al., 2019).

| Predicting fluxes of NCI-60 cells using ecModel
The model reconstruction and flux detection of 11 NCI-60 cell lines were derived from Zenodo (https://doi.org/10.5281/zenodo.3577466). Note: as each constructed ecModel has over 20,000 reactions and has a different model framework from the COBRA toolbox, performing MCMC sampling is difficult. Therefore, the fluxes for each ecModel were estimated using the suggested method "minProSimulation" from Human 1 (Robinson et al., 2020).

| Comparisons of the predicted and experimentally measured fluxes
The predicted secretion fluxes using different methods, that is, GPMM, GIMME, and Fastcore, were compared with the experimental flux datasets as mentioned above. To avoid linear optimization precision error, we removed the fluxes with absolute values smaller than 1e-6 mmol/L/min. Pearson correlations between predicted fluxes and experimental measurements were calculated to compare the performance of the different methods.

| Robustness analysis of GPMM
We first constructed noise-induced transcriptomes by adding random numbers (viz. artificial noise) to the original expression data (or the genuine transcriptome) of each gene. Specifically, 1%, 5%, and 10% noise was induced in each NCI-60 cell line transcriptome. In these processes, the noise inducing is produced from a uniform distribution of [0.99, noise. For example, in the 5% noise translation procedure, if a gene in a cell line has a gene expression (e.g., fpkm) of 1.0, the adding random number of this gene is ranged from −0.05 to 0.05, such as 0.03; then, the noised-induced gene expression is a number ranged from 0.95 to 1.05, such as 1.03. Second, we performed the metabolic modeling on the genuine and the noise-induced transcriptomes and compared the obtained flux results to evaluate GPMM robustness. In addition, to further test whether multiple sampling affects robustness, we also induced 5% noise to the gene expression of H460 cell line 100 times and performed metabolic modeling to determine the stability of GPMM. As shown in Table 1 and Table S4, 96% of CENs lived independently (e.g., eating, walking, and talking). Compared with F1SPs, CENs had significantly higher diastolic blood pressure (146.0 vs. 137.9, p = 0.03), similar systolic blood pressure (83.2 vs. 86.1 p = 0.19), slightly lower blood glucose (5.98 vs. 6.70, p = 0.14, t-test), lower total cholesterol (4.68 vs. 5.43 p = 0.009, t-test), and lower low-density lipoprotein cholesterol (2.45 vs. 2.97, p = 0.043, t-test). These results are also consistent with our previous studies, where levels of risk factors for cardiovascular diseases, including blood glucose, triglyceride, and total cholesterol, were significantly lower in the CENs than those of the general older population from the same province, and the diagnoses of type 2 diabetes mellitus, hypertriglyceridemia, and hypertension were lower in the CENs than Chinese national levels (He et al., 2014).

| The Chinese centenarians study
The relatively healthy status of CENs suggests that they can serve as a good model for healthy aging studies.
For transcriptome analysis, peripheral blood samples were treated with red blood cell lysis buffer (Tiangen Biotech) and then centrifuged at 1800 g for 10 min to isolate white blood cells. For metabolomics and proteomics measurements, peripheral blood samples were allowed to clot at room temperature for 30 min and then centrifuged for 10 min at 1500 g to extract the serum.

| Metabolic modeling
Gene expression (FPKM) levels in 170 individuals from the Hainan longevity cohort, including 76 CENs, 52 centenarian-children, and 42 spouses of centenarian-children (F1SPs), were derived from our previously published data (Xiao et al., 2018). The upper bounds of the white blood cell metabolite uptake rates were separated into three categories: nutrient uptake for energy production, cofactors, and iron/oxygen uptake (Table S5). The nutrient uptake rates, including those of glucose, l-glutamine, and fatty acids, were derived from published literature (Table S5). The essential amino acid uptake rates were set to a small number, whereas those of cofactors, iron, oxygen, and primers for glycogen synthesis were set to unlimited (Table   S5). After preparing the gene expression and nutrient uptake rates, genome-wide precision metabolic modeling was conducted using our developed GPMM toolbox.

| Identification of metabolic flux profiles of centenarians
As the CENs (aged 98-108) were older than the controls (45-75), we used two linear models to distinguish centenarian signatures and aging effects. Model l was applied to CEN and F1SP samples to determine unfiltered centenarian signature: Model 2 was applied to F1SP samples to determine aging effects: After the unfiltered centenarian signatures and aging effect were determined, we obtained the actual centenarian signatures by excluding the overlapping fluxes between the unfiltered centenarian signature and the age effect. Therefore, upregulated fluxes in CENs were defined as fluxes with p < 0.05 and beta >0 in model 1 but not significant or beta <0 in model 2. Downregulated fluxes were defined vice versa.

| Identifying significant metabolic subsystems using differential abundance scores
The differential abundance (DA) score was calculated using previously published methods (Hakimi et al., 2016). For each metabolic subsystem, the ith DA score (DA i ) was calculated as follows: To obtain the significance of DA scores, we used a "bootstrap without replacement" method to calculate p-values. Briefly, we first randomly shuffled the sample labels 1000 times. Second, for each randomly shuffled label, the corresponding random DA scores were calculated using the above formula. We thus obtained 1000 random (4) Model 1: flux ∼ lm (centenarians + sex) (5) Model 2: flux ∼ lm (age + sex) where DA i is the centenarian DAscore of the ith subsystem and RandomDAs i is the random DAscore of the ith subsystem. The adjusted p-value was calculated using the false discovery rate (FDR).

| Sample preparation
The collected serum samples were thawed on ice. Samples (100 μl) were extracted with 750 μl of methanol/acetonitrile/water solution (V methanol :V acetonitrile :V water = 2:2:1), with 30 μl of 1 mg/ml l-2chlorophenylalanineas then added as an internal standard, followed by vortexing for 10 s and sonicating for 10 min on ice. After that, the extract was incubated for 1 h at −20°C. Following centrifugation . Ion source gas 1 was 60, ion source gas 2 was 60, curtain gas was set to 30 L/h, source temperature was set to 550°C, and ion spray voltage floating (ISVF) was set to 5500 V or −4500 V in positive or negative modes, respectively. 4.5.3 | Data preprocessing and annotation UPLC-MS raw data (.wiff) were converted to mzXML, with ProteoWizard Peak exaction, identification, integration, alignment, and retention time correction processed with XCMS (R package, v3.2). The preprocessing results generated a data matrix that consisted of retention time (RT), mass-to-charge ratio (m/z) values, and peak intensity. The R package CAMERA was used for peak annotation after XCMS data processing. The in-house MS2 database was applied for metabolite identification, and only the metabolites with MS2 >0.8 remained.

| Identification of metabolomic profiles of centenarians
Using equations to those shown in flux analysis, we also used two linear models to distinguish centenarian signatures and aging effects in the metabolomic data. Model l was applied to CEN and F1SP samples to determine unfiltered centenarian signatures: Model 2 was applied to F1SP samples to determine aging effects: Then, the actual centenarian signature was calculated by excluding the overlapping metabolites between the unfiltered centenarian signature and the age effect. Therefore, upregulated metabolites in CENs were defined as the log2 transformed metabolic abundance with p < 0.05 and beta >0 in model 1 but not significant or beta <0 in model 2. Downregulated metabolites were defined vice versa.

| Identifying significant metabolic classes in metabolomic data
Similar to the flux subsystem analysis in Equations (6) and (7). The significance of the metabolic class was also analyzed using the differential abundance (DA) score (for Figure 4c). For each metabolic class, the DA score was calculated as the number of upregulated metabolites minus the downregulated metabolites and then divided by the total number of metabolites in the given class. Similarly, as presented in flux subsystem analysis, we used a "bootstrap without replacement" method to calculate p-values and then adjusted these p-values using the false discovery rate (FDR).

ACK N OWLED G EM ENTS
This work was supported jointly and equally by grants from the