- Open Access
Genome-scale network model of metabolism and histone acetylation reveals metabolic dependencies of histone deacetylase inhibitors
© The Author(s). 2019
- Received: 8 November 2018
- Accepted: 21 February 2019
- Published: 1 March 2019
Histone acetylation plays a central role in gene regulation and is sensitive to the levels of metabolic intermediates. However, predicting the impact of metabolic alterations on acetylation in pathological conditions is a significant challenge. Here, we present a genome-scale network model that predicts the impact of nutritional environment and genetic alterations on histone acetylation. It identifies cell types that are sensitive to histone deacetylase inhibitors based on their metabolic state, and we validate metabolites that alter drug sensitivity. Our model provides a mechanistic framework for predicting how metabolic perturbations contribute to epigenetic changes and sensitivity to deacetylase inhibitors.
- Genome-scale metabolic modeling
- Metabolic network
- Protein acetylation
- Metabolic regulation
- Protein deacetylases (KDACs)
- KDAC inhibitors
Protein acetylation is a highly conserved mechanism to regulate many cellular processes including transcription and metabolism. The substrate for acetylation, acetyl-CoA, is a central metabolic intermediate at the crossroads of anabolic and catabolic pathways . Protein acetylation is thus sensitive to the metabolic state of the cell . Biochemical assays in Saccharomyces cerevisiae and human cell lines have shown that levels of acetyl-CoA directly affect protein acetylation [3–5]. Hundreds of proteins, including metabolic enzymes, are regulated by acetylation [6, 7]. Acetylation can also influence gene expression through post-translational modification of histones. Cells rely on histone acetylation to increase chromatin accessibility and influence gene expression [2, 8].
Given its pervasive regulatory role, altered acetylation is believed to play a part in a variety of diseases including cancer and metabolic disorders such as diabetes, obesity, dyslipidemia, and hypertension [5, 9–11]. Since metabolic alterations and dysregulation of protein acetylation are important cancer hallmarks, understanding the interplay between these processes can reveal novel therapeutic targets against cancer. However, predicting the interplay between these two processes is challenging due to acetyl-CoA’s pervasive role in metabolism, and due to the highly interconnected nature of the metabolic network. No theoretical approach exists to predict the impact of the change in cellular metabolism on protein acetylation.
Creating a model of metabolism and protein acetylation can enable the prediction of the impact of nutrient shifts or mutations in metabolic enzymes on the epigenome. This can shed light on metabolic and chromatin dysregulation during tumorigenesis [12, 13]. Compounds that disrupt acetylation machinery such as deacetylase inhibitors are increasingly used for treating cancers and metabolic and immune disorders . Predicting the interplay between metabolism and acetylation can identify cancer cells that are sensitive to deacetylase inhibitors based on their metabolic state.
To address this challenge, here we develop a computational model of metabolism and protein acetylation using constraint-based modeling (CBM). CBM makes use of metabolic network reconstructions that represent the mechanistic relationships between genes, proteins, and metabolites within a biological system. CBM has been successfully used to predict the metabolic state of various mammalian systems, including cancer cells and stem cells [14–17]. We hypothesized that protein acetylation dynamics can be inferred from the metabolic network topology and stoichiometry.
We demonstrate that our metabolic model can explain known acetylation changes associated with nutrient excess and starvation based on the availability of carbon units. We then apply our acetylation model to predict and validate the impact of cellular metabolic state on sensitivity to drugs that disrupt acetylation, specifically protein deacetylase inhibitors that are currently used in the clinic for anticancer therapy. Our approach allowed us to predict the variation in sensitivity between deacetylase inhibitors based on their unique impact on cellular metabolism.
Simulating the effect of the metabolic state on acetylation
To simulate the influence of metabolism on acetylation, a nuclear protein acetylation reaction (protein + acetyl-CoA → acetyl-protein + CoA) was incorporated into the human metabolic network reconstruction by Duarte et al., which contains 3747 reactions, 1496 ORFs, 2004 proteins, and 2766 metabolites . A nuclear ATP citrate lyase reaction and nuclear transport of citrate and oxaloacetate were also included to enable synthesis of acetyl-CoA in the nucleus based on recent biochemical evidence . Since acetyl-CoA and substrates for acetyl-CoA synthesis can diffuse between the cytosol and nucleus through the nuclear pore, the flux through the protein acetylation reaction is representative of acetylation changes in both cytosol and nuclear proteins.
The metabolic model recapitulated known acetylation changes occurring due to a shift in the metabolic state. Switching to high or low glucose significantly affected acetylation levels, which are known to change with glucose flux [4, 19] (Fig. 2a). Bulk histone acetylation in yeast cells has also been shown to be responsive to glucose availability . Addition of acetate, hypoxia, and starvation of amino acids also increased acetylation flux, consistent with experimental observations [19, 22–26]. To further validate the impact of media components on acetylation predicted by the model, we compared predictions with data from recent biochemical studies [19, 27–29] (Fig. 2c). Specifically, these studies measured the impact of addition or removal of amino acids, glucose, acetate, minerals, and vitamins on histone acetylation in cancer cell lines. Our model correctly predicted the impact of 19 out of the 20 nutrient changes. These results suggest that acetylation levels can be predicted based on the carbon flux and energetic state of the cell.
The ability of the model to reproduce these diverse observations provided new insights on the interaction between acetylation and metabolism. Our analysis suggests that bulk acetylation reflects excess carbon flux that is not used for biomass synthesis. Hence, the model predicted acetylation even in the absence of growth, which was observed experimentally in conditions with a carbon source but without a nitrogen source (i.e., depletion of amino acids). In contrast, growth conditions without excess carbon, such as those with amino acids alone, did not display histone acetylation (Fig. 2c).
The inconsistent observation between the model and experiment arose in a media condition lacking the carbon source—glucose. The model did not predict any acetylation, while it was observed experimentally. Acetylation in this condition was only possible in the model if there was an additional carbon source such as acetate or fatty acids. A recent study has found significant levels of acetate in DMEM media, which is usually not accounted for in biochemical studies . This acetate can replenish cellular acetyl-CoA pools , and our modeling suggests that this can potentially support histone acetylation. We then re-evaluated the model predictions after including all nutrients that can be present in trace amounts in serum—amino acids from uptake or degradation of serum proteins, trace amounts of fatty acids, and acetate. These components are sufficient to support acetylation in this condition that contained all major nutrients except glucose. However, the trace nutrients were insufficient to promote acetylation if no major nutrients were present (Additional file 1: Figure S1). Thus, this inconsistent prediction points to the potential role of additional components in culture media that can impact acetylation.
Metabolic modeling predicts bulk histone acetylation levels in diverse cancer cell lines
Since FBA correctly identified the impact of distinct metabolic conditions on acetylation levels, we next tested if FBA can predict more nuanced differences in the basal metabolic state and acetylation levels between different cancer cell lines. We used bulk histone acetylation levels from LeRoy et al.’s study, which measured acetylation levels using quantitative proteomics across a panel of diverse cell lines representing lung, liver, colon, and brain tumors .
To predict cell-type-specific changes in acetylation levels due to their underlying metabolism, we inferred their metabolic state from gene expression data . Metabolic network models have been applied successfully to predict metabolic behaviors of human tissues and cancer cells using transcriptomic data [33–36]. We obtained gene expression data for the cell lines in LeRoy et al. study from the Cancer Cell Line Encyclopedia (CCLE) project  and integrated it with the human metabolic network model using the iMAT approach  (Methods).
Predicting metabolic enzymes that affect protein acetylation
Functional significance of the 24 metabolic genes predicted by the model to impact acetylation
2. Interacts with acetylase?
3. Interacts with deacetylase?
4. Sirtuin target?
ATP citrate lyase
Aconitase 2, mitochondrial
Aldehyde dehydrogenase 1 family, member L1
Fumarylacetoacetate hydrolase (fumarylacetoacetase)
Glutamic-oxaloacetic transaminase 1, soluble (aspartate aminotransferase 1)
Glutathione transferase zeta 1
Isocitrate dehydrogenase 2 (NADP+), mitochondrial
Succinate dehydrogenase complex, subunit A, flavoprotein
Succinate dehydrogenase complex, subunit B, iron sulfur (Ip)
Succinate dehydrogenase complex, subunit C
Succinate dehydrogenase complex, subunit D
Serine hydroxymethyltransferase 1 (soluble)
Solute carrier family 25 (mitochondrial carrier; citrate transporter)
Our analysis uncovered several novel enzymes in addition to those that have been previously reported to impact acetylation. To assess the biological significance of this prediction, we compared the 24 acetylation-impacting metabolic genes with data from the protein acetylation database , which catalogs all known acetylation targets in humans. Proteomic studies have identified several metabolic enzymes as targets of acetylation [6, 7, 39]. Proteins encoded by 18 of the 24 acetylation-impacting genes were themselves targets of protein acetylation. In addition, a significant fraction of the enzymes predicted to affect acetylation by the model were found to be direct targets of Sirtuins, a major class of deacetylases. The fraction of enzymes that were found to be acetylated and deacetylated among those predicted by the metabolic model were significantly higher than would be expected by chance (p value = 8 × 10− 15 for acetylation and p value = 9 × 10− 6 for deacetylation respectively; Table 1). Intriguingly, none of the 39 acetylation-increasing genes were targets of Sirtuins, suggesting that Sirtuins regulate reactions involved in production but not the consumption of acetyl units.
To further corroborate our predictions of enzymes that impact acetylation, physical interaction data from literature were mined for associations between metabolic enzymes and acetylation machinery. Proteins involved in the same cellular process are more likely to physically interact . We hypothesized that the metabolic enzymes identified by our model to affect acetylation levels will have strong physical interactions with acetylation or deacetylated-related enzymes.
Using data from the Biogrid database , which contains physical interaction data curated from literature, interactions between metabolic and acetylation-related enzymes in human cells were identified. Several metabolic enzymes predicted by our model to influence acetylation interacted physically with either a protein acetylase or deacetylase enzyme (Additional file 1: Table S2, Table S3). This overlap is significantly higher than expected by randomly choosing a metabolic enzyme (p value = 0.01 and 0.002 for acetylation-impacting and acetylation-enhancing enzymes respectively, hypergeometric test). Further, we compared the extent of overlap with acetylases and deacetylases separately. We observed a significant overlap between acetylation-enhancing genes with protein deacetylases but not with acetylases (p value = 0.0004 and 0.07 respectively). In contrast, acetylation-impacting genes have a significant number of interactions with both protein acetylases and deacetylases (p value = 0.003 and 0.005 respectively). Thus, acetylases predominantly interact with enzymes that lead to the synthesis of acetyl units. For example, the acetylation-impacting enzyme ACLY had physical interactions with several acetylases. ACLY is known to be transported to the nucleus to facilitate acetylation [42, 43]. The physical interaction suggests that the product of the ACLY metabolic reaction (acetyl-CoA) can be directly channeled as substrate for acetylases. We also found physical interactions between the TCA cycle enzymes citrate synthase (CS) and succinate dehydrogenase (SDH) with deacetylase enzymes HDAC5 and SIRT7 respectively. Overall, our results suggest that enzymes identified by our model to impact acetylation show strong functional association with protein acetylation pathways, thus corroborating our predictions.
Altering sensitivity of cancer cells to protein deacetylase inhibitors by changing their metabolic state
Our analysis so far suggests that genome-scale modeling can predict the impact of the metabolic state on global protein acetylation. We next predicted the impact of the cellular metabolic state on susceptibility to drugs that disrupt cellular protein acetylation, specifically protein deacetylase inhibitors. The human genome encodes 18 different lysine deacetylases (KDACs) that regulate numerous biochemical pathways . KDACs are considered attractive therapeutic targets for treating not just cancers, but also viral infections, inflammation, neurodegenerative diseases, and metabolic disorders . Four deacetylase inhibitors (vorinostat, romidepsin, belinostat, panobinostat) are approved clinically for cancer treatment . However, identifying tumors that are sensitive to these drugs is a significant clinical challenge [45, 46].
Deacetylase enzymes catalyze the removal of acetylation marks from histones and other proteins. They are sensitive to, and directly impact, the metabolic state of the cell . Deacetylase inhibitors block the deacetylation of proteins and are predicted to cause cell death by increasing histone acetylation and triggering apoptosis [45, 46, 48]. Hence, we hypothesized that conditions that favor increased acetylation should enhance sensitivity to these drugs.
Since genome-scale modeling allowed us to predict metabolic conditions that increased protein acetylation, we next used it to predict conditions that enhance sensitivity to deacetylase inhibitors. Specifically, we explored the use of metabolic network modeling to identify cell types that are sensitive to KDAC inhibition based on their metabolic and acetylation state. This can enable in silico identification of tumors most sensitive to these drugs.
To validate our hypothesis, we experimentally measured the sensitivity of the ovarian cancer cell line (HeLa) to vorinostat (also known as suberoylanilide hydroxamic acid (SAHA)). Vorinostat was the first approved histone deacetylase inhibitor. HeLa cells were cultured in 92 different metabolic conditions using Biolog phenotype microarrays and exposed to vorinostat or control (DMSO). The culture conditions in the Biolog arrays included a wide range of carbon sources that result in distinct patterns of metabolic activity in the cell . This broad range of metabolic activity allowed us to assess the impact of the metabolic state of the cell on vorinostat sensitivity.
While the predicted sensitivity showed high concordance with experiments across most substrates, the model incorrectly over-predicted acetylation flux and corresponding sensitivity in conditions with maltose and galactose as carbon sources; galactose is metabolized differently despite similarity to glucose . Intriguingly, vorinostat enhanced the growth of cells under some conditions. KDAC inhibitors are known to reduce glucose metabolism and increase mitochondrial oxidative metabolism . This metabolic impact may be beneficial in conditions that lack glucose or other fermentable sugars. We observed a strong correlation between the extent of growth enhancement and acetylation flux predicted in those conditions (Additional file 1: Figure S2). Hence, KDAC inhibitors may enhance growth in conditions with low acetylation while inhibiting growth in conditions with high acetylation flux. This dependence on the cellular metabolic state may explain the growth-promoting effect of KDAC inhibitor treatment on stem cells and neurons [52, 53].
Predicting cellular sensitivity to protein deacetylase inhibitors using the metabolic state
Since distinct metabolic conditions showed differences in sensitivity to vorinostat, we next tested if our model can predict more nuanced differences in metabolic activity between different cancer cell lines and their corresponding impact on sensitivity to KDAC inhibitors. Large-scale therapeutic screening data from Seashore-Ludlow et al. was used to identify variation in sensitivity to KDAC inhibitors among a panel of 860 cancer cell lines encompassing various tissue and tumor types . The transcriptomic state of these cell lines has been characterized using microarrays as part of the Cancer Cell Line Encyclopedia (CCLE) project .
We compared the sensitivity of the cell lines with high acetylation flux state with those with low acetylation flux. Cell lines with high acetylation flux were found to be significantly more sensitive to vorinostat treatment than those with low acetylation flux (Fig. 5c; p value = 4 × 10− 6, t-test). This is consistent with the observation across metabolic conditions using Biolog phenotype arrays. The low throughput Biolog assays revealed a strong quantitative relationship with the model predictions. However, these were observed with a single cell line. In contrast, the high-throughput drug-screening assays from Seashore-Ludlow et al. involve hundreds of cell lines. Our results across this large dataset suggest that this is a statistically robust phenomenon.
To assess the generality of this observation to other KDAC inhibitors, we compared the pattern of sensitivity between high and low acetylation flux cell lines using data from three other clinically used KDAC inhibitors that were part of the Seashore-Ludlow et al. dataset—panobinostat, entinostat, and belinostat. Entinostat is currently in phase III clinical trials, and the other two drugs are FDA approved. All these drugs were significantly more sensitive to high acetylation flux state cell lines (p value < 10− 6, t-test; Fig. 5c, Additional file 1: Figure S3). Overall, these results suggest that the metabolic and acetylation state is predictive of sensitivity to KDAC inhibitors.
Finally, to assess the specificity of the acetylation flux metric in predicting sensitivity to KDAC inhibitors, we compared the differential sensitivity of all 481 compounds that were part of the Seashore-Ludlow et al. study between the high and low acetylation flux cell lines. KDAC inhibitors were significantly more strongly selective between the two groups than other drug classes (p value = 0.003, t-test; Additional file 1: Figure S4). This trend applies to both experimental and clinically used KDAC inhibitors (Additional file 1: Figure S4). Thus, the observed trend is not due to an indirect correlation with other cellular processes such as the growth rate. Indeed, we observed a stronger trend using the Seashore-Ludlow et al. data after correcting for doubling time of each cell type (p value = 7 × 10− 9, Additional file 1: Figure S5).
Predicting variability between protein deacetylase inhibitors
So far, we have predicted the impact of the cellular metabolic state on sensitivity to protein deacetylase inhibitors. We next incorporated the direct impact of these drugs on cellular metabolism. Vorinostat treatment has been previously shown to reduce glucose and fatty acid metabolism in myeloma cell lines . Incorporating the impact of these drugs on metabolism would allow us to predict differences in sensitivity resulting from both basal metabolic state and the unique impact of each KDAC inhibitor. KDAC inhibitors display variability in their sensitivity to different cell types . KDAC inhibitors block distinct histone deacetylases resulting in transcriptomic changes that influence the levels of metabolic enzymes [48, 53].
To model the unique effect of each drug, we interpreted the pattern of differential gene expression from cells exposed to these drugs with metabolic network models. We used data measured in one representative cell line and used it to predict the metabolic impact of each KDAC inhibitor on all CCLE cell lines. We analyzed transcriptomic signatures of four different KDAC inhibitors (vorinostat, panobinostat, entinostat, and belinostat) before and after treatment from UKN1 cells . Several metabolic genes were uniquely differentially expressed in response to each KDAC inhibitor, which may lead to variation in sensitivity (Additional file 1: Table S6). This list of downregulated and upregulated genes was overlaid onto all the CCLE cell line metabolic models to predict variation in drug sensitivities among the CCLE cell lines. The net acetylation flux due to both basal transcriptomic state and drug-induced transcriptomic changes were determined for all cell lines for the four KDAC inhibitors using the iMAT approach.
In this study, we built a genome-scale biochemical model to predict the impact of cellular metabolic state and availability of nutrients on acetylation. A key observation from our study is that histone acetylation dynamics reflects the flux of excess carbon that is not used for biomass synthesis (overflow metabolism). Overflow metabolism has also been proposed to drive acetylation in prokaryotes [56, 57]. Conditions that lack amino acids or oxygen can increase acetylation as there is excess carbon available for acetylation from glucose. Strikingly, we observed acetylation even in certain nutrient starvation conditions that do not support growth, which is also observed experimentally (Fig. 2c).
Our analysis also identified known and novel metabolic enzymes that can impact acetylation. These predictions were supported by numerous physical interactions between the acetylation-impacting metabolic enzymes predicted by our model with proteins in the acetylation machinery. This suggests that the activity of these metabolic enzymes can impact protein acetylation and are in turn regulated by protein acetylation. This reciprocal regulation to maintain homeostasis has also been observed between metabolic pathways and signaling pathways [58–60] .
Since acetylation plays a central role in signaling and chromatin regulation, mutations in the acetylation-impacting metabolic genes can significantly affect cellular homeostasis. The acetylation-impacting genes, IDH2, SDHA/B/C, and FH, are annotated as being frequently mutated in cancers in the COSMIC database . Aberrant expression of these genes may also lead to altered acetylation in tumors. The genes ACLY, CS, RPE, and IDH2 were found to be among the top 10% of the genes that were frequently overexpressed across 19 cancer types from a meta-analysis study of 1981 tumors (FDR p value < 0.05) . These genes predicted by our model are potential drug targets for disorders with dysregulated epigenome.
We then demonstrated that our model can predict the impact of the metabolic state on sensitivity to deacetylase inhibitors. We hypothesized that conditions that favor increased acetylation should enhance sensitivity to these drugs, which cause cell death by blocking deacetylation. Consistent with this hypothesis, we were able to predict and validate both metabolic conditions and cell types that show increased sensitivity to these drugs. Our experimental screen across diverse metabolic conditions revealed that sensitivity to vorinostat changed significantly with the metabolic state. Notably, our biochemical model accurately predicted the sensitivity to vorinostat by predicting the acetylation flux in diverse nutrient conditions. Nutrients that were modeled in this study such as glucose, acetate, amino acids, and lactate are also metabolized by tumors in vivo [63, 64]. While glucose is the common carbon source, some tumors also use lactate . Our analysis suggests that while glucose condition increases sensitivity, growth in lactate may offer protection against vorinostat. Hence, this drug may be effective only against some tumor types. This is potentially relevant for targeting cancers with distinct metabolic states . While the tumor microenvironment is complex and dynamic, our ability to predict the impact of drug efficacy in simple defined environments is the first step towards addressing the in vivo complexity. Our study paves the way for dynamic models that can simulate the impact of changing the metabolic conditions.
Analysis of KDAC inhibitor sensitivity across 860 cancer cell lines revealed that differences in basal metabolic state among cells can lead to variation in sensitivity. Using metabolic modeling, we were able to group cell lines based on their acetylation state and predict those that show increased sensitivity to deacetylase inhibition. Identifying tumors that are sensitive to specific KDAC inhibitors is a significant clinical challenge . Our approach can help predict efficacy of deacetylase inhibitors against different cancer types using transcriptomic data. Transcriptomic profiling is increasingly done for tumors for treatment stratification and prognosis . A very small tissue sample is needed for transcriptomics compared to the requirements for cell culture and drug screening.
Furthermore, unique differences in the transcriptomic pattern induced by each KDAC inhibitor correlate with variation in sensitivity between these inhibitors. By analyzing sensitivity to four clinically used KDAC inhibitors across all CCLE cell lines, we found that cell lines that showed higher acetylation flux reduction after treatment with a KDAC inhibitor relative to others were more sensitive to the corresponding drug. Our approach complements machine-learning and high-throughput screening approaches  by providing a rational mechanistic strategy to predict and tune sensitivity of KDAC inhibitors based on acetylation flux. Genome-scale modeling has been shown to accurately predict the metabolic state of cancer cells from transcriptomic data. Future analysis using proteomics and metabolomics data could lead to tumor stratification with greater accuracy.
While we have focused on the impact of metabolism and nutrition on acetylation in this study, the activity of acetylases and metabolic enzymes can also be regulated by signaling pathways which can subsequently impact bulk acetylation levels . Further, while our model predicts bulk acetylation levels, it does not provide insights on acetylation changes in specific protein sites. In future, knowledge of specific targets, binding affinity, and expression levels of acetylation and deacetylation enzymes can enable prediction of acetylation dynamics at specific sites in the proteome. Similarly, incorporating the role of redox factors, FADH and NADH, that can influence activity of deacetylases can potentially resolve inconsistencies observed in nutrient conditions such as galactose media that impact redox state . Similar to acetylation, which is sensitive to the flux of acetyl-CoA, recent studies have revealed that histone methylation is sensitive to the flux of the substrate S-Adenosyl Methionine (SAM) [17, 69]. Our study can pave the way for modeling other histone modifications like acylation, succinylation, and methylation that also depend on metabolic intermediates.
Histone acetylation links gene regulation with metabolism, and aberrant acetylation is associated with cancers and metabolic disorders. Here we developed a genome-scale biochemical network model to predict the impact of metabolism on acetylation. Our model revealed that histone acetylation levels change predictably depending on cellular carbon and nutrient availability. We identified several novel metabolic genes that can impact acetylation when deleted. Our approach provided two key insights related to cellular metabolism and sensitivity to deacetylase inhibition. Firstly, our experimental screen across diverse metabolic conditions and analysis across 860 cancer cell lines revealed that sensitivity to deacetylase inhibitors can be predicted from the cellular metabolic state. Secondly, unique rewiring of metabolic pathways induced by each KDAC inhibitor correlates with variability between these inhibitors. This is potentially relevant for tumor treatment stratification based on the metabolic state. Overall, our study demonstrates the strong and predictable interconnection between metabolism and acetylation.
Flux balance analysis (FBA)
In FBA, an optimal metabolic flux state is determined that maximizes an objective (usually the biomass). We used the biomass composition from Shlomi et al.  for the growth objective function. The biomass objective function represents the steady-state consumption of 42 essential metabolites, such as nucleotides, amino acids, and lipids required for cellular proliferation. In addition to maximizing the objective, the following constraints are satisfied—(1) the system is at steady state with mass balanced for all species (i.e., no accumulation or depletion of any intracellular metabolites in the network) and (2) the specified upper and lower bounds for each reaction in the network are not violated. FBA was used for simulating the impact of gene deletions and metabolic environment on protein acetylation. A total of 50 nutrients commonly present in rich media culture were used as starting points for nutrient addition in excess and nutrient removal.
FBA was performed in different nutrient environments with the protein acetylation reaction as the secondary objective (Optimization objective: maximize Biomass_reaction + epsilon × acetylation_reaction). Acetylation flux is optimized after biomass synthesis is maximized. We used a parameter (“epsilon”) to incorporate the relative weights for optimizing biomass flux relative to acetylation. Changing the value of epsilon (i.e., relative weights) over a 10,000-fold range from 0.01 to 10− 6 does not influence the value of acetylation flux and the accuracy of our approach (Additional file 1: Figure S7, Additional file 1: Table S7). The two objectives optimized here can be conflicting or potentially independent depending on nutrient availability. Since acetylation flux is a secondary objective for our optimization, any value of epsilon that is relatively small compared to biomass function would lead to the same value of acetylation flux. The mathematical approach used here with separate weights for different objectives is the most common way of solving such multi-objective problems . This approach is scalable to any number of objectives.
Total carbon flux necessary for acetylation compared to other biosynthetic processes
Using data from total potential acetylation sites in the human genome , we find that diverting a small fraction of the cytoplasmic flux of acetyl-CoA would be sufficient to acetylate all the sites in a cell. The observed flux through cytosolic acetyl-CoA synthesis reaction in human iBMK cells is 2 nmol/h/μlcells , which is equivalent to 109 molecules per hour in each cell. Assuming a total of 109 acetylation sites in a mammalian cell , redirecting 1/24th of the cytosolic acetyl-CoA production flux towards the nucleus in a period of 24 h can be sufficient to saturate all the histone acetylation sites, assuming a relatively low rate of deacetylation . The mitochondrial carbon flux towards acetyl-CoA is 10-fold higher than cytosolic flux . This suggests that redirection of carbon flux towards acetylation will not significantly impact biomass synthesis.
Transcriptomic data integration
To integrate transcriptomic data with the metabolic network model, we used a linear version of the iMAT approach  implemented in the PROM algorithm . In the iMAT approach, the number of active reactions associated with genes that are upregulated or highly expressed is maximized, while flux through reactions associated with downregulated reactions are minimized. The linear optimization version of the iMAT approach allows for continuous restriction of fluxes instead of setting reactions to be completely on or off . The constraints are flexible, i.e., they allow for violation of individual constraints in order to maximize the overall consistency with the data. These inconsistencies may represent sites of regulation by other regulatory mechanisms . We used this approach to infer a cell-type-specific metabolic state using transcriptomic data for CCLE cell lines and to incorporate transcriptomic evidence from KDAC inhibitor profiling data. The optimization problem was solved using the Gurobi mathematical programming solver. The entire series of steps for predicting acetylation flux from transcriptomic data is outlined in Additional file 1: Figure S8 and Additional file 1: Figure S9.
To assess the robustness of our results to the method of transcriptome integration, we repeated all our analyses with an alternative approach for transcriptomic data integration called GIMME . This approach differs from the iMAT approach in that it assumes a direct linear association between extent of downregulation of a gene’s expression and the downregulation of flux through the corresponding enzyme encoded by it. We found that our results from this approach are consistent with our results using iMAT (Additional file 1: Table S8).
Small molecule therapeutic screening data
Area under the growth curve (AUC) values for the small molecule therapeutic screening data for 481 compounds against 860 cell lines were obtained from Seashore-Ludlow et al. along with media used for each of the 860 cell lines. The media composition and transcriptomic data for these cell lines from CCLE was used to infer the metabolic network state for these cell lines as described above. In order to estimate the impact of the cellular doubling rate on drug sensitivity, we also used data from Hafner et al. that controls for this effect . The growth-controlled efficacy metric was found to be superior to existing metrics for assessing the efficacy of drugs in dividing cells. Cell lines with missing growth inhibition values for the KDAC inhibitors were removed from the statistical analysis.
Biolog phenotype microarray (PM) measurement
The phenotype microarray (PM) plates from Biolog Inc. (Hayward, CA, USA) contain diverse energy sources to profile the metabolic capabilities of a cell. Each well in a PM plate contains a single chemical as the sole energy source and production of NADH per well is monitored using a colorimetric redox dye chemistry. As the cells metabolize the energy source, tetrazolium dye in the media is reduced, producing a purple color according to the amount of NADH generated. The PM-M1 plate used in this study contains 92 different carbon energy sources along with negative control wells.
Biolog PM assays were performed as described in Boccuto et al. . Viable cells were counted utilizing a TC20™ Automated Cell Counter. PM-M1 plates were incubated with 20,000 viable cells per well in a volume of 50 μL. The cells were incubated at 37 °C in 5% CO2, using the modified Biolog IF-M1 medium . HeLa cells were exposed to vorinostat (2 μM) or DMSO (control) at the time they were added to the Biolog PM plates, for a 72-h exposure. The experiment was conducted in four replicates using the PM-M1 plate. During the final 24 h of exposure, the optical density of each well was measured every 15 min. The optical readings were normalized using quadruplicate readings from an empty plate (plates run with no cells, just media and dye). The normalized optical density readings were used for statistical analysis.
To control for the impact of slow or no growth, we used only those conditions that supported any growth above negative control (Background). Using various thresholds for growth did not impact the strong correlation between the predicted acetylation flux and growth inhibition (Additional file 1: Figure S2).
The significance of overlap between gene lists was estimated using the hypergeometric test in MATLAB. Biolog data analysis was performed in R.
We are grateful to members of the Chandrasekaran laboratory for feedback and suggestions.
This work was supported by faculty start-up funds from the University of Michigan to SC. The funding body had no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.
Availability of data and materials
The MATLAB implementation of the metabolism-histone acetylation model, the source code for simulating the model and associated data sets are publicly available at the Synapse software repository under the GNU General Public License (http://synapse.org/MetabolismEpigenomeModel) (DOI: https://doi.org/10.7303/syn17114770) . Instructions for running the model and reproducing the data in the manuscript are also provided.
The Biolog phenotype microarray data generated in this study is provided in Additional file 1: Table S4. The small molecule therapeutic screening data is available as Supplementary material (Supplementary Table S3) from Seashore-Ludlow et al. . Bulk histone acetylation levels were obtained from the supplementary information of the LeRoy et al. study (Additional file 1) . Gene expression data for 14 of the 24 cell lines in LeRoy study and 860 cell lines from Seashore-Ludlow et al. were obtained from the Cancer Cell Line Encyclopedia project (CCLE) . Transcriptomics signatures before and after treatment of four different KDAC inhibitors in UKN1 cells were obtained from the Rempel et al. study (Supplementary material 2) . The Synapse website also contains these above-mentioned data sets that were used for building and testing the model.
FZ performed the computational analysis. LB, SS, and RP performed the experimental analysis. SC conceived the study, designed research, supervised, and performed research. SC wrote the manuscript with inputs from other authors. All authors read and approved the final manuscript.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Pietrocola F, Galluzzi L, Bravo-San Pedro JM, Madeo F, Kroemer G. Acetyl coenzyme A: a central metabolite and second messenger. Cell Metab. 2015;21:805–21.View ArticleGoogle Scholar
- Kaelin WG, McKnight SL. Influence of metabolism on epigenetics and disease. Cell. 2013;153:56–69.View ArticleGoogle Scholar
- Cai L, Sutter BM, Li B, Tu BP. Acetyl-CoA induces cell growth and proliferation by promoting the acetylation of histones at growth genes. Mol Cell. 2011;42:426–37.View ArticleGoogle Scholar
- Cluntun AA, Huang H, Dai L, Liu X, Zhao Y, Locasale JW. The rate of glycolysis quantitatively mediates specific histone acetylation sites. Cancer Metab [Internet]. 2015;3:10. Available from: http://www.cancerandmetabolism.com/content/3/1/10
- Lee JV, Shah SA, Wellen KE. Obesity, cancer and acetyl-CoA metabolism. Drug Discov Today Dis Mech. 2013;10:e55–61.View ArticleGoogle Scholar
- Choudhary C, Weinert BT, Nishida Y, Verdin E, Mann M. The growing landscape of lysine acetylation links metabolism and cell signalling. Nat Rev Mol cell Biol. 2014;15:536–50.View ArticleGoogle Scholar
- Drazic A, Myklebust LM, Ree R, Arnesen T. The world of protein acetylation. Biochim Biophys Acta - Proteins Proteomics. 2016;1864(10):1372–401.Google Scholar
- Lu C, Thompson CB. Metabolic regulation of epigenetics. Cell Metab. 2012;16:9–17.View ArticleGoogle Scholar
- Kaypee S, Sudarshan D, Shanmugam MK, Mukherjee D, Sethi G, Kundu TK. Aberrant lysine acetylation in tumorigenesis: implications in the development of therapeutics. Pharmacol Ther. 2016;162:98–119.View ArticleGoogle Scholar
- Kinnaird A, Zhao S, Wellen KE, Michelakis ED. Metabolic control of epigenetics in cancer. Nat Rev Cancer. 2016;16(11):694.Google Scholar
- Wagner GR, Payne RM. Mitochondrial acetylation and diseases of aging. J Aging Res. 2011;2011:13. Article ID 234875. https://doi.org/10.4061/2011/234875.
- Wong CC, Qian Y, Yu J. Interplay between epigenetics and metabolism in oncogenesis: mechanisms and therapeutic approaches. Oncogene. 2017;36(24):3359.Google Scholar
- Martinez-Pastor B, Cosentino C, Mostoslavsky R. A tale of metabolites: the cross-talk between chromatin and energy metabolism. Cancer Discov. 2013;3:497–501.View ArticleGoogle Scholar
- Bordbar A, Monk JM, King ZA, Palsson BO. Constraint-based models predict metabolic and associated cellular functions. Nat Rev Genet. 2014;15:107–20.View ArticleGoogle Scholar
- Nilsson A, Nielsen J. Genome scale metabolic modeling of cancer. Metab Eng. 2017;43:103–12.Google Scholar
- Yizhak K, Chaneton B, Gottlieb E, Ruppin E. Modeling cancer metabolism on a genome scale. Mol Syst Biol. 2015;11:817.View ArticleGoogle Scholar
- Chandrasekaran S, Zhang J, Ross C, Huang Y-C, Asara J, Li H, et al. Comprehensive mapping of pluripotent stem cell metabolism using dynamic genome-scale network modeling. Cell Rep. 2017;21(10):2965–77.Google Scholar
- Duarte NC, Becker SA, Jamshidi N, Thiele I, Mo ML, Vo TD, et al. Global reconstruction of the human metabolic network based on genomic and bibliomic data. Proc Natl Acad Sci U S A. 2007;104:1777–82.View ArticleGoogle Scholar
- Wellen KE, Hatzivassiliou G, Sachdeva UM, Bui TV, Cross JR, Thompson CB. ATP-citrate lyase links cellular metabolism to histone acetylation. Science (80- ). 2009;324:1076–80.View ArticleGoogle Scholar
- Orth JD, Thiele I, Palsson BØ. What is flux balance analysis? Nat Biotechnol Nature Research. 2010;28:245–8.View ArticleGoogle Scholar
- Fan J, Kamphorst JJ, Mathew R, Chung MK, White E, Shlomi T, et al. Glutamine-driven oxidative phosphorylation is a major ATP source in transformed mammalian cells in both normoxia and hypoxia. Mol Syst Biol. 2013;9(1):712.Google Scholar
- Galdieri L, Vancura A. Acetyl-CoA carboxylase regulates global histone acetylation. J Biol Chem. 2012;287:23865–76.View ArticleGoogle Scholar
- Galdieri L, Zhang T, Rogerson D, Lleshi R, Vancura A. Protein acetylation and acetyl coenzyme a metabolism in budding yeast. Eukaryot Cell. 2014;13:1472–83.View ArticleGoogle Scholar
- Takahashi H, McCaffery JM, Irizarry RA, Boeke JD. Nucleocytosolic acetyl-coenzyme a synthetase is required for histone acetylation and global transcription. Mol Cell. 2006;23:207–17.View ArticleGoogle Scholar
- Bánréti Á, Sass M, Graba Y. The emerging role of acetylation in the regulation of autophagy. Autophagy. 2013;9(6):819–29.Google Scholar
- Chaveroux C, Jousse C, Cherasse Y, Maurin A, Parry L, Carraro V, et al. Identification of a novel amino acid response pathway triggering ATF2 phosphorylation in mammals. Mol Cell Biol. 2009;29(24 ):6515–26.Google Scholar
- McBrian MA, Behbahan IS, Ferrari R, Su T, Huang TW, Li K, et al. Histone acetylation regulates intracellular pH. Mol Cell. 2013;49(2):310–21.Google Scholar
- Gao X, Lin SH, Ren F, Li JT, Chen JJ, Yao CB, et al. Acetate functions as an epigenetic metabolite to promote lipid synthesis under hypoxia. Nat Commun. 2016;7:11960.Google Scholar
- McDonnell E, Crown SB, Fox DB, Kitir B, Ilkayeva OR, Olsen CA, et al. Lipids reprogram metabolism to become a major carbon source for histone acetylation. Cell Rep. 2016;17(6):1463–72.Google Scholar
- Kamphorst JJ, Chung MK, Fan J, Rabinowitz JD. Quantitative analysis of acetyl-CoA production in hypoxic cancer cells reveals substantial contribution from acetate. Cancer Metab. 2014;2(1):23.Google Scholar
- LeRoy G, DiMaggio PA, Chan EY, Zee BM, Blanco M, Bryant B, et al.. A quantitative atlas of histone modification signatures from human cancer cells. Epigenetics Chromatin [Internet]. 2013;6:20. Available from: http://epigeneticsandchromatin.biomedcentral.com/articles/10.1186/1756-8935-6-20
- Shlomi T, Cabili MN, Herrgård MJ, Palsson BØ, Ruppin E. Network-based prediction of human tissue-specific metabolism. Nat Biotechnol. 2008;26:1003–10.View ArticleGoogle Scholar
- Frezza C, Zheng L, Folger O, Rajagopalan KN, MacKenzie ED, Jerby L, et al. Haem oxygenase is synthetically lethal with the tumour suppressor fumarate hydratase. Nature. 2011;477:225–8.View ArticleGoogle Scholar
- Folger O, Jerby L, Frezza C, Gottlieb E, Ruppin E, Shlomi T. Predicting selective drug targets in cancer through metabolic networks. Mol Syst Biol EMBO Press; 2014;7:501.Google Scholar
- O’Brien EJ, Monk JM, Palsson BO. Using genome-scale models to predict biological capabilities. Cell Elsevier. 2015;161:971–87.Google Scholar
- Uhlén M, Fagerberg L, Hallström BM, Lindskog C, Oksvold P, Mardinoglu A, et al. Proteomics. Tissue-based map of the human proteome. Science. 2015;347:1260419.View ArticleGoogle Scholar
- Barretina J, Caponigro G, Stransky N, Venkatesan K, Margolin AA, Kim S, et al. The Cancer Cell Line Encyclopedia enables predictive modelling of anticancer drug sensitivity. Nature. 2012;483:603–7.View ArticleGoogle Scholar
- Liu Z, Wang Y, Gao T, Pan Z, Cheng H, Yang Q, et al. CPLM: a database of protein lysine modifications. Nucleic Acids Res. 2014;42:D531–6.View ArticleGoogle Scholar
- Choudhary C, Kumar C, Gnad F, Nielsen ML, Rehman M, Walther TC, et al. Lysine acetylation targets protein complexes and co-regulates major cellular functions. Science (80- ). 2009;325:834–40.View ArticleGoogle Scholar
- Rhodes DR, Tomlins SA, Varambally S, Mahavisno V, Barrette T, Kalyana-Sundaram S, et al. Probabilistic model of the human protein-protein interaction network. Nat Biotechnol. 2005;23(8):951.Google Scholar
- Chatr-Aryamontri A, Oughtred R, Boucher L, Rust J, Chang C, Kolas NK, et al. The BioGRID interaction database: 2017 update. Nucleic Acids Res. 2017;45(D1):D369–D379.Google Scholar
- Sivanand S, Viney I, Wellen KE. Spatiotemporal control of acetyl-CoA metabolism in chromatin regulation. Trends Biochem Sci. 2018;43(1):61–74.Google Scholar
- Sivanand S, Rhoades S, Jiang Q, Lee JV, Benci J, Zhang J, et al. Nuclear acetyl-CoA production by ACLY promotes homologous recombination. Mol Cell. 2017;67(2):252–65.Google Scholar
- Roche J, Bertrand P. Inside HDACs with more selective HDAC inhibitors. Eur J Med Chem. 2016;121:451–83.Google Scholar
- Guha M. HDAC inhibitors still need a home run, despite recent approval. Nat Rev Drug Discov [Internet]. Nature Publishing Group, a division of Macmillan Publishers Limited. All Rights Reserved.; 2015;14:225. Available from: https://doi.org/10.1038/nrd4583
- Li Y, Seto E. HDACs and HDAC inhibitors in cancer development and therapy. Cold Spring Harb Perspect Med. 2016;6(10):a026831.Google Scholar
- Su X, Wellen KE, Rabinowitz JD. Metabolic control of methylation and acetylation. Curr Opin Chem Biol. 2016;30:52–60.Google Scholar
- Van Lint C, Emiliani S, Verdin E. The expression of a small fraction of cellular genes is changed in response to histone hyperacetylation. Gene Expr. 1996;5:245–53.PubMedGoogle Scholar
- Bochner BR, Siri M, Huang RH, Noble S, Lei XH, Clemons PA, et al. Assay of the multiple energy-producing pathways of mammalian cells. PLoS One. 2011;6(3):e18147.Google Scholar
- Gohil VM, Sheth SA, Nilsson R, Wojtovich AP, Lee JH, Perocchi F, et al. Nutrient-sensitized screening for drugs that shift energy metabolism from mitochondrial respiration to glycolysis. Nat Biotechnol. 2010;28:249–55.View ArticleGoogle Scholar
- Wardell SE, Ilkayeva OR, Wieman HL, Frigo DE, Rathmell JC, Newgard CB, et al. Glucose metabolism as a target of histone deacetylase inhibitors. Mol Endocrinol [Internet]. 2009;23:388–401 Available from: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=2654518&tool=pmcentrez&rendertype=abstract.View ArticleGoogle Scholar
- Bug G, Gül H, Schwarz K, Pfeifer H, Kampfmann M, Zheng X, et al. Valproic acid stimulates proliferation and self-renewal of hematopoietic stem cells. Cancer Res. 2005;65(7):2537–41.Google Scholar
- Eckschlager T, Plch J, Stiborova M, Hrabeta J. Histone deacetylase inhibitors as anticancer drugs. Int J Mol Sci. 2017.Google Scholar
- Seashore-Ludlow B, Rees MG, Cheah JH, Coko M, Price EV, Coletti ME, et al. Harnessing connectivity in a large-scale small-molecule sensitivity dataset. Cancer Discov. 2015;5:1210–23.View ArticleGoogle Scholar
- Rempel E, Hoelting L, Waldmann T, Balmer NV, Schildknecht S, Grinberg M, et al. A transcriptome-based classifier to identify developmental toxicants by stem cell testing: design, validation and optimization for histone deacetylase inhibitors. Arch Toxicol. 2015;89:1599–618.View ArticleGoogle Scholar
- Nakayasu ES, Burnet MC, Walukiewicz HE, Wilkins CS, Shukla AK, Brooks S, et al. Ancient regulatory role of lysine acetylation in central metabolism. MBio. 2017;8(6):e01894–17.Google Scholar
- Schilling B, Christensen D, Davis R, Sahu AK, Hu LI, Walker-Peddakotla A, et al. Protein acetylation dynamics in response to carbon overflow in Escherichia coli. Mol Microbiol. 2015;98(5):847–63.Google Scholar
- Wellen KE, Thompson CB. A two-way street: reciprocal regulation of metabolism and signalling. Nat Rev Mol Cell Biol. 2012.Google Scholar
- Jiang P, Du W, Mancuso A, Wellen KE, Yang X. Reciprocal regulation of p53 and malic enzymes modulates metabolism and senescence. Nature. 2013.Google Scholar
- Gut P. Verdin E. Nature: The nexus of chromatin regulation and intermediary metabolism; 2013.Google Scholar
- Forbes SA, Beare D, Gunasekaran P, Leung K, Bindal N, Boutselakis H, et al. COSMIC: exploring the world’s knowledge of somatic mutations in human cancer. Nucleic Acids Res. 2015;43:D805–11.View ArticleGoogle Scholar
- Nilsson R, Jain M, Madhusudhan N, Sheppard NG, Strittmatter L, Kampf C, et al. Metabolic enzyme expression highlights a key role for MTHFD2 and the mitochondrial folate pathway in cancer. Nat Commun. 2014;5:3128.Google Scholar
- Liu X, Cooper DE, Cluntun AA, Warmoes MO, Zhao S, Reid MA, et al. Acetate production from glucose and coupling to mitochondrial metabolism in mammals. Cell. 2018;175(2):502–13.Google Scholar
- Faubert B, Li KY, Cai L, Hensley CT, Kim J, Zacharias LG, et al. Lactate metabolism in human lung tumors. Cell. 2017;171(2):358–71.Google Scholar
- Heiden MG Vander, Cantley LC, Thompson CB. Understanding the warburg effect: the metabolic requirements of cell proliferation. Science. 2009;324(5930):1029–33.Google Scholar
- Roychowdhury S, Chinnaiyan AM. Translating cancer genomes and transcriptomes for precision oncology. CA Cancer J Clin. 2016;66(1):75–88.Google Scholar
- Geeleher P, Loboda A, Lenkala D, Wang F, LaCroix B, Karovic S, et al. Predicting response to histone deacetylase inhibitors using high-throughput genomics. J Natl Cancer Inst. 2015;107(11):djv247.Google Scholar
- Marroquin LD, Hynes J, Dykens JA, Jamieson JD, Will Y. Circumventing the Crabtree effect: replacing media glucose with galactose increases susceptibility of hepG2 cells to mitochondrial toxicants. Toxicol Sci. 2007;97:539–47.View ArticleGoogle Scholar
- Mentch SJ, Mehrmohamadi M, Huang L, Liu X, Gupta D, Mattocks D, et al. Histone methylation dynamics and gene regulation occur through the sensing of one-carbon metabolism. Cell Metab. 2015;22:861–73.View ArticleGoogle Scholar
- Shlomi T, Benyamini T, Gottlieb E, Sharan R, Ruppin E. Genome-scale metabolic modeling elucidates the role of proliferative adaptation in causing the Warburg effect. PLoS Comput biol [internet]. 2011/03/23. 2011;7:e1002018. Available from: http://www.ncbi.nlm.nih.gov/pubmed/21423717.
- Deb K. Multi-objective optimization using evolutionary algorithms. Ser: Syst. Optim; 2001.Google Scholar
- Chandrasekaran S, Price ND. Probabilistic integrative modeling of genome-scale metabolic and regulatory networks in Escherichia coli and Mycobacterium tuberculosis. Proc Natl Acad Sci. 2010;107:17845–50.View ArticleGoogle Scholar
- Becker SA, Palsson BO. Context-specific metabolic networks are consistent with experiments. PLoS Comput Biol [Internet]. 2008/05/17. 2008;4:e1000082. Available from: http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?cmd=Retrieve&db=PubMed&dopt=Citation&list_uids=18483554
- Hafner M, Niepel M, Chung M, Sorger PK. Growth rate inhibition metrics correct for confounders in measuring sensitivity to cancer drugs. Nat Methods. 2016;13:521–7.View ArticleGoogle Scholar
- Boccuto L, Chen C-F, Pittman AR, Skinner CD, McCartney HJ, Jones K, et al. Decreased tryptophan metabolism in patients with autism spectrum disorders. Mol Autism [Internet]. 2013;4:16 Available from: http://molecularautism.biomedcentral.com/articles/10.1186/2040-2392-4-16.View ArticleGoogle Scholar
- Shen F, Chandrasekaran S. Genome-scale network model of metabolism and histone acetylation [internet]. Synapse. 2019. https://doi.org/10.7303/syn17114770 Available from: http://synapse.org/MetabolismEpigenomeModel.
- Broad Institute Cancer cell line encyclopaedia (CCLE) [Internet]. Available from: www.broadinstitute.org/ccle