Leveraging metabolic modeling and machine learning to uncover modulators of quiescence depth

Abstract

Quiescence, a temporary withdrawal from the cell cycle, plays a key role in tissue homeostasis and regeneration. Quiescence is increasingly viewed as a continuum between shallow and deep quiescence, reflecting different potentials to proliferate. The depth of quiescence is altered in a range of diseases and during aging. Here, we leveraged genome-scale metabolic modeling (GEM) to define the metabolic and epigenetic changes that take place with quiescence deepening. We discovered contrasting changes in lipid catabolism and anabolism and diverging trends in histone methylation and acetylation. We then built a multi-cell type machine learning model that accurately predicts quiescence depth in diverse biological contexts. Using both machine learning and genome-scale flux simulations, we performed high-throughput screening of chemical and genetic modulators of quiescence and identified novel small molecule and genetic modulators with relevance to cancer and aging.

Keywords: aging; cancer; genome-scale models; machine learning; quiescence deepening.

Fig. 1.

Graphical overview. In this article, we take a multifaceted approach to identifying modulators of quiescence depth: (i) metabolic modeling and flux simulations that enable in silico genetic knockouts and (ii) a multi-cell type transcriptomics quiescence depth predictor that facilitates screening chemical and genetic perturbations for their effect on quiescence depth. In doing so, we simulate metabolic and epigenetic features of quiescence deepening and explore quiescence depth in embryogenesis and aging. Predictions are validated using a cancer genomics database.

Fig. 2.

Metabolic pathway analysis of quiescence deepening. A) Workflow for simulating metabolic reaction fluxes. Transcriptomics from nine stages of quiescence depth were integrated into a metabolic model (Rat-GEM) to constrain the set of allowable fluxes (using INIT). Then, we used CHRR to randomly sample from the set of potential fluxes for each quiescence depth. B) Top subsystems associated with quiescence depth. For each subsystem, each reaction in the subsystem was correlated with quiescence depth and the average Pearson correlation was computed. Subsystems with the highest average Pearson correlation coefficients are shown here (Q-value < 1E-5, one-sample t test). A positive mean correlation coefficient implies a positive association with quiescence depth while a negative mean correlation coefficient indicates a negative association with quiescence depth. C) Subsystem enrichment analysis of reactions negatively correlating with QDS at a P-value of <0.15. Subsystems (pathways) ranked by fold enrichment are shown here along with an enrichment P-value calculated using the hypergeometric function and corrected with the Benjamini–Hochberg adjustment (Q-value < 0.01 subsystems shown here). D) Subsystem enrichment analysis of reactions positively correlating with QDS at a P-value of <0.15. Subsystems (pathways) ranked by fold enrichment are shown here along with an enrichment P-value calculated using the hypergeometric function and corrected with the Benjamini–Hochberg adjustment (Q-value < 0.01 subsystems shown here). Top 10 subsystems meeting criteria are shown here, ranked by fold enrichment. E) Normalized reaction fluxes across quiescence depths for all reactions negatively correlating with QDS at a P-value of <0.15. F) Normalized reaction fluxes across quiescence depths for all reactions positively correlating with QDS at a P-value of <0.15.

Fig. 3.

Lipid anabolism and lipid catabolism display largely opposing trends during quiescence deepening. A) Mean Pearson correlation coefficient of all reactions with QDS for lipid anabolism subsystems. Only subsystems with a mean correlation coefficient different from zero (Q-value < 0.05, one-sided t test) are shown here. B) Mean Pearson correlation coefficient of all reactions with QDS for lipid catabolism subsystems. Only subsystems with a mean correlation coefficient different from zero (Q-value < 0.05, one-sided t test) are shown here. C) Metabolic pathway plot for example lipid anabolism subsystem: fatty acid biosynthesis (odd-chain). Reactions are colored by the Pearson correlation coefficient with quiescence depth. D) Metabolic pathway plot for example lipid catabolism subsystem: beta oxidation of di-unsaturated fatty acids (N-6) (mitochondrial). Reactions are colored by the Pearson correlation coefficient with quiescence depth. E) Distribution of Pearson correlation coefficients with quiescence depth for transcripts belonging to the GO biological process “lipid biosynthesis.” Only correlation coefficients with Q-value < 0.05 are shown here. F) Distribution of Pearson correlation coefficients with quiescence depth for transcripts belonging to the GO biological process “lipid catabolism.” Only correlation coefficients with Q-value < 0.05 are shown here.

Fig. 4.

Quiescence deepening is associated with diverging trends in epigenetic histone acetylation and methylation. A) Diagram of chromatin highlighting histone post-translational modifications and their relationship to metabolism. SAM is the substrate for histone methyltransferases (HMT) to methylate histone residues and convert SAM to S-adenosylhomocysteine (SAH). Acetyl-CoA is the substrate for histone acetyltransferases (HAT) to acetylate histone residues with CoA as a byproduct. Citrate is converted to acetyl-CoA by ATP citrate lyase (ACLY). These reactions are simulated in the metabolic model. B) Simulated mean nuclear histone acetylation flux vs. quiescence depth reveals a downward trend in histone acetylation as quiescence deepens. Pearson correlation performed between quiescence depth and mean nuclear histone acetylation flux. C) Top subsystems whose knockout reverse the acetylation trend, ranked by Pearson correlation between QDS and histone acetylation upon subsystem knockout. D) Simulated mean nuclear histone methylation flux vs. quiescence depth reveals an upward trend in histone methylation with increasing quiescence depth. Pearson correlation performed between quiescence depth and mean nuclear histone methylation flux. E) Top subsystems whose knockout reverses the methylation trend, ranked by Pearson correlation between QDS and histone methylation upon subsystem knockout.

Fig. 5.

Gene knockout simulations predict a range of quiescence depth modulators in alignment with experimental data. A) Overview of gene knockout simulations. Transcriptomics from shallow and deep quiescence states were used to constrain a metabolic model and perform flux simulations, enabling the simulation of the metabolic state of a cell in shallow and deep quiescence. Then, we used the rMTA to predict gene knockouts that could move a cell to (i) a shallower quiescence state or (ii) a deeper quiescence state. B) Accuracy of gene perturbation simulations for top 50 hits. To validate these results experimentally, we used a cancer genomic database called CancerMine. The CancerMine database validation relied on the following principle: genes whose knockdown lowered quiescence depth was considered validated if they were annotated as tumor suppressors (genes whose downregulation promotes cancer growth), genes whose knockdown raised quiescence were considered validated if they were annotated as oncogenes (genes whose downregulation lowers cancer growth). C) Enrichment analysis. For quiescence deepeners, we calculated the fold enrichment and P-value of tumor suppressors among the top 50 hits compared to the background list of all genes knocked out with rMTA. For quiescence reducers, we calculated the fold enrichment of cancer oncogenes or drivers compared to the background set of genes. P-values are computed using the hypergeometric function and listed above the bar chart bars. **P-value < 0.01. D) Top 20 hits of the gene perturbation simulation for quiescence reducers: The top 20 genes whose knockout is predicted to move the cells toward a deeper quiescence state are shown here (ranked by rMTA robust transformations core). Here, the CancerMine annotation described in B) is shown for each prediction. E) Top 20 hits of the gene perturbation simulation for quiescence deepeners: The top 20 genes whose knockout is predicted to move the cells toward a shallower quiescence state are shown here (ranked by rMTA robust transformations core). As before, the CancerMine annotation described in B) is shown here for each prediction. The gene MPC1 is highlighted. This gene has experimental evidence for its inhibition lowering quiescence depth (as predicted here).

Fig. 6.

Flux-based machine learning model captures metabolic signatures of quiescence deepening for multiple types of input data. A) Diagram of the machine learning model. Reaction fluxes are inputted, and the model is trained to predict QDS (days since serum starvation) using elastic net regression. B) Top subsystems enriched in the model's positive coefficients (up with quiescence deepening) and negative coefficients (down with quiescence deepening). The top 10 subsystems with Q-value < 0.05 for the hypergeometric fold enrichment test are shown here. C) K-fold cross-validation accuracy plot. The model is trained and optimized on 90% of the data and tested in the remaining 10% for 10 random subsets of the data. Note: all parameter are un-optimized in training set except for fold-change threshold, which is set close to one to maximize the impact of genes on fluxes. D) External validation shows that the model predicts deeper quiescence (higher QDS) for cells in a quiescent phase vs. a proliferative phase for multiple types of quiescence and proliferation induction and multiple types of input omics data (transcriptomics, translatomics, and proteomics). Fluxes are generated using genes upregulated in quiescence relative to proliferation vs. fluxes generated using genes upregulated in proliferation relative to quiescence with a fold change threshold of 1.00. (i) Proliferation is induced using epidermal growth factor (EGF), horse serum, insulin, and hydrocortisone. Quiescence is induced with serum starvation, Mek inhibition, CDK46 inhibition, and contact inhibition. Transcriptomics measured. Comparisons were performed with a two-sample t test. Source: GSE122927 (41). (ii) Proliferation is induced using EGF or tissue plasminogen activator (TPA). Quiescence induced with serum starvation. Transcriptomics measured. Comparisons were performed with a two-sample t test. Source: Ref. (42). (iii) QDS for proliferating vs. quiescent cells. Proliferation is induced by IL-3. Quiescence is induced with IL-3 starvation. Proteomics measured. Comparisons were performed with a paired t test. Source: Ref. (43). (iv) QDS for proliferating vs. quiescent cells. Proliferation is induced using EGF. Quiescence induced with serum starvation. Translatomics measured. Comparisons were performed with a paired t test. Source: GSE112295 (44). (v) Proliferation is induced using fetal bovine serum (FBS). Quiescence induced with 7-day contact-inhibition. Transcriptomics measured. Comparisons were performed with a two-sample t test. Source: GSE117444 (45). (vi) Proliferation is induced using FBS. Quiescence induced with serum starvation. Transcriptomics measured. Comparisons were performed with a two-sample t test. Source: GSE74620 (46).

Fig. 7.

Multi-cell type machine learning model accurately predicts QDS in a wide range of contexts. A) Diagram of the machine learning model. Transcriptomics from three datasets measuring quiescence deepening is inputted, and the model is trained to predict QDS (days since serum starvation, normalized to the maximum of each dataset) using elastic net regression. B) GO pathway enrichment analysis for positive and negative coefficients of model. The top eight biological processes for negative coefficients (down in quiescence deepening, Q-value < 0.65) and positive coefficients (up in quiescence deepening, Q-value < 0.05) as ranked by Q-value of enrichment. C) K-fold cross-validation accuracy plot. The model is trained and optimized on 90% of the data and tested in the remaining 10% for 10 random subsets of the data. D) External validation: model predictions for proliferating vs. quiescent cells in various cell types and experimental conditions. (i) Proliferation is induced using FBS. Quiescence induced with 7-day contact-inhibition. Comparisons were performed with a two-sample t test. Source: GSE117444 (45). (ii) Proliferation is induced using FBS. Quiescence induced with serum starvation. Comparisons were performed with a two-sample t test. Source: GSE74620 (46). (iii) Proliferation is induced using EGF, horse serum, insulin, and hydrocortisone. Quiescence is induced with serum starvation, Mek inhibition, CDK46 inhibition, and contact inhibition. Comparisons were performed with a two-sample t test. Source: GSE122927 (41). (iv) QDS for proliferating vs. quiescent cells. Proliferation is induced using EGF. Quiescence induced with serum starvation. Comparisons were performed with a two-sample t test. Source: GSE112295 (44). (v) Proliferation is induced using EGF or TPA. Quiescence induced with serum starvation. Comparisons were performed with a two-sample t test. Source: Ref. (42).

Fig. 8.

Changes in quiescence depth are associated with stem cell aging and embryogenesis. A) HSC aging: (i) Older HSCs have a trend toward higher QDS than young HSCs. Source: GSE109546. (ii) Older LT HSCs and ST HSCs have higher QDS than young cells of the same type. No aging trend was observed in MPPs. Single-cell RNA-seq of aging HSCs. Source: GSE59114. B) Median QDS and age correlate in mouse subventricular zone (SVZ) neural stem cell (NSC) niche (Pearson correlation). Source: GSE104651. C) During embryogenesis, quiescence depth is predicted to fall and then rise between the zygote and 32-cell state. Single cell RNA-seq measured for preimplantation mouse blastomeres. Source: GSE136714. D) GO transcriptomic analysis of proliferation markers. Shown here are the mean z-score-normalized transcript levels for transcripts belonging to GO cell growth and GO stem cell proliferation. While GO cell growth-related genes show no definitive pattern, stem cell proliferation-related genes rise and then fall in a pattern that (except for the zygote stage) aligns with the fall and rise in QDS seen in early embryogenesis.

Fig. 9.

Small molecules and gene knockouts predicted to alter quiescence align strongly with experimental results. A) Overview of perturbation screening: L1000 data measured as part of the library of integrated network-based cellular signatures program for two classes of perturbations (small molecules and gene knockouts) is inputted into the QDS predictor to screen perturbations for their effect on quiescence depth. B) Top small molecules whose knockout is predicted to alter quiescence depth. For each of the 1,713 small molecules tested, conditions are controlled to a 24-h timepoint and 10 μM dose across multiple cell lines. Top 25 predictions with Q-value < 0.05 shown here, as ranked by magnitude of quiescence depth score modulation. Q-value is the Benjamin–Hochberg-corrected P-value of a two-sample t test between predicted QDS for the control condition (drug vehicle) and treatment condition (small molecule compound). C) Top genes whose knockout is predicted to alter quiescence depth (Q-value < 0.05). D) Enrichment of tumor suppressors among quiescence deepeners (genes whose knockout lowers QDS and increases proliferative potential), and enrichment of oncogenes/cancer drivers among quiescence reducers (genes whose knockout raises QDS and decreases proliferative potential). P-value calculated using the hypergeometric function. Results are shown at three Q-value thresholds for selecting top gene knockout hits: 0.05, 0.1, and 0.85. Cancer annotations sourced from the CancerMine database. Genes with a majority of tumor suppressor annotations in CancerMine are classified as tumor suppressors, and genes with a majority of cancer driver or oncogene annotations are classified as cancer drivers.