Quantitative determinants of aerobic glycolysis identify flux through the enzyme GAPDH as a limiting step

Abstract
eLife digest
Introduction
Results
Discussion
Materials and methods
References
Article and author information
Metrics

Abstract

Aerobic glycolysis or the Warburg Effect (WE) is characterized by the increased metabolism of glucose to lactate. It remains unknown what quantitative changes to the activity of metabolism are necessary and sufficient for this phenotype. We developed a computational model of glycolysis and an integrated analysis using metabolic control analysis (MCA), metabolomics data, and statistical simulations. We identified and confirmed a novel mode of regulation specific to aerobic glycolysis where flux through GAPDH, the enzyme separating lower and upper glycolysis, is the rate-limiting step in the pathway and the levels of fructose (1,6) bisphosphate (FBP), are predictive of the rate and control points in glycolysis. Strikingly, negative flux control was found and confirmed for several steps thought to be rate-limiting in glycolysis. Together, these findings enumerate the biochemical determinants of the WE and suggest strategies for identifying the contexts in which agents that target glycolysis might be most effective.

https://doi.org/10.7554/eLife.03342.001

eLife digest

Cells generate energy from a sugar called glucose via a process called glycolysis. This process involves many enzymes that catalyze 10 different chemical reactions, and it essentially converts glucose step-by-step into a simpler chemical called pyruvate.

Pyruvate is then normally transported into structures within the cell called mitochondria, where it is further broken down using oxygen to release more energy. However, in cells that are rapidly dividing, pyruvate is converted into another chemical called lactate—which releases energy more quickly, but releases less energy overall. Cancer cells often convert most of their glucose into lactate, rather than breaking down pyruvate in their mitochondria: an observation known as the ‘Warburg effect’. And while many factors affect how a cell releases energy from pyruvate, it remains unclear what regulates which of these biochemical processes is most common in a living cell.

In this study, Shestov et al. have developed a computational model for the process of glycolysis and used this to investigate the causes of the Warburg Effect. The model was based on the known characteristics of the enzymes and chemical reactions involved at each step. It predicted that the activity of the enzyme called GAPDH, which carries out the sixth step in glycolysis, in many cases affects how much lactate is produced. This suggests that this enzyme represents a bottleneck in the pathway.

Next, Shestov et al. performed experiments where they used drugs to block different stages of the glycolysis pathway, and confirmed that the GAPDH enzyme is important for regulating this pathway in living cancer cells too. In these treated cells, the levels of a chemical called fructose-1,6-biphosphate (which is made in a step in the pathway between glucose and pyruvate) were either very high or very low. Shestov et al. proposed that the flow of chemicals through the glycolysis pathway is controlled by the GAPDH enzyme when the chemicals used by the enzymes upstream of GAPDH in the pathway (which includes fructose-1,6-biphosphate) are plentiful. However, if these chemicals are limited, other enzymes that are involved in earlier steps of the pathway regulate the process instead.

The findings of Shestov et al. reveal that the regulation of glycolysis is more complex than previously thought, and is also very different when cells are undergoing the Warburg Effect. In the future, these findings might help to identify the types of cancer that could be effectively treated using drugs that target the glycolysis process, which are currently being tested in pre-clinical studies.

https://doi.org/10.7554/eLife.03342.002

Introduction

Proliferating cells increase their glucose consumption and secrete lactate as opposed to completely oxidizing the glucose in the mitochondria (Warburg et al., 1927) and is known as aerobic glycolysis or the Warburg Effect. Currently, this altered metabolism is exploited for diagnostics and is subjected to multiple drug development efforts (Koppenol et al., 2011; Vander Heiden, 2011; Hamanaka and Chandel, 2012). Numerous studies have identified genes such as KRAS, PIK3CA, and cMYC and microenvironments such as hypoxia and hypoglycemia that promote aerobic glycolysis but a complete understanding of the necessary and sufficient biochemical alterations associated with this phenotype is unknown. Furthermore successful translation for biomedical applications is limited by understanding the contexts in which therapies that target glycolysis might be effective.

Computational modeling has a successful history in the study of metabolism (Rapoport et al., 1976; Fell, 1992; Schilling et al., 1999; Cascante et al., 2002). Genome-scale stoichiometric models of metabolism have been developed to study the effects of drug targets in human metabolism and have had success in predicting the WE (Molenaar et al., 2009; Vazquez et al., 2010; Folger et al., 2011; Shlomi et al., 2011). However, a comprehensive quantitative understanding of the WE requires knowledge of enzyme activities and metabolic control.

Therefore, we collected and integrated multiple forms of data into a modeling framework involving flux balances of glycolysis, detailed chemical kinetics based on reaction mechanisms and parameters measured, physico-chemical constraints from thermodynamics and mass conservation, metabolic control analysis, and Monte Carlo sampling of parameter space. We next use mass spectrometry and isotope tracing to probe concentrations and fluxes through the pathway and their responses to several perturbations. Together, we elaborate the determinants of aerobic glycolysis and identify and confirm novel points of regulation in glycolysis that have remained unidentified for over 50 years since the discovery of the pathway.

Results

Biochemical kinetic model of aerobic glycolysis

We investigated the kinetics of the glycolytic pathway from glucose uptake to oxidation of pyruvate in the mitochondria or export of lactate out of the cytosol. We modeled each step of the pathway according to enzymatic mechanism and known modes of allosteric control resulting in a set of differential equations (Figure 1A, ‘Materials and methods’, Supplementary file 1). While it is not possible to model every possible interaction explicitly, the aim is to capture enough of the pathway so that a large range of experimentally realized measurements can be obtained and relationships between variables can be observed.

Figure 1

Download asset Open asset

A quantitative model and statistical simulation method captures the diversity of metabolic states observed in tumor and proliferating cells.

(A) Schematic of the glycolysis model with chemical reactions and allosteric points of regulation described. Abbreviations: GLC—glucose, G6P—glucose-6-phosphate, F6P—fructose-6-phosphate, FBP—fructose-1,6,-bisphosphate, F26BP—fructose-2,6,-bisphosphate, GAP—glcyceraldehyde-3-phosphate, DHAP—dihydroxyacetone phosphate, BPG—1,3 bisphosphoglycerate, 3PG—3-phosphoglycerate, 2PG—2-phosphoglycerate, PEP—phosphoenolpyruvate, PYR—pyruvate, SER—Serine, GLY—glycine, Lac—lactate, MAL—malate, ASP—aspartate, Pi—inorganic phosphate, CI—creatine, PCI—phosphophocreatine, GTR—glucose transporter, HK—hexokinase, PGI—phosphoglucoisomerase, PFK—phosphofructokinase, ALD—aldolase, TPI—triosephosphoisomerase, GAPDH—glyceraldehyde-phosphate dehydrogenase, PGK—phosphoglycerate kinase, PGM—phosphoglycerate mutase, ENO—enolase, PK—pyruvate kinase, LDH—lactate dehydrogenase, MCT—monocarboxylate transporter, PDH—pyruvate dehydrogenase, CK—creatine kinase. (B) Overview of the algorithm and simulation method. (C) Measured values of the NADH/NAD⁺ ratio across a population of MCF10A breast epithelial cells. Three values of glucose concentration are considered (0.5 mM blue, 5.5 mM green, and 25 mM red). (D) Calculated fluxes (mM/hr) for glycolysis rate (Glycolysis) are defined as the rate of glucose to pyruvate (per molecule of pyruvate), pyruvate to lactate flux (LDH), rate of oxygen consumption (OxPhos), rate of NADH turnover (NADH), and ATP turnover (ATPase). (E) Calculated probability density function (PDF) of NAD⁺ concentrations. (F) Calculated probability density function (PDF) of NADH/NAD⁺ ratio. (G) Calculated probability density function (PDF) of ATP concentrations. (H) Calculated probability density function (PDF) of ATP/ADP ratio. (I) Box plots showing the distribution of concentrations computed from the simulation for each intermediate in glycolysis.

https://doi.org/10.7554/eLife.03342.003

Since glycolysis is the most extensively studied biochemical pathway, there is a wealth of information on the kinetic parameters and enzyme expression that govern the equations. Nevertheless, it is also not possible to capture cellular physiology in any biochemical model with single values of kinetic parameters (Daniels et al., 2008). This difficulty arises from the tremendous amount of heterogeneity within cells at multiple levels. The origins of this heterogeneity vary from genetic variation observed across cancer types, tumor types, differences in signaling mechanisms that affect post-translational modifications in each cell, and the differences in microenvironmental pressures (e.g., the oxygen availability) that each cell within a given tumor experiences, as well as the inherent cell to cell variation common to all cells (Marusyk et al., 2012). Therefore, we developed an integrated algorithm to evaluate the statistics of the kinetics of glycolysis that accounts for the possible variation within metabolism (Figure 1B, ‘Materials and methods’).

First the model is constrained using mass conservation constraints that conserve the balance of glucose, redox state, and energy status. Next, thermodynamic constraints are used to constrain fluxes according to the free energy of the reactions determined by Haldane relationships. These physical constraints are combined with the kinetic mechanisms that define each step of glycolysis and the chemical reactions involving the redox-associated metabolites NAD⁺, and NADH and energy-associated metabolites ATP, ADP, and AMP. Next, the model is constrained to expression data so that protein concentrations are subsumed in the Vmax values and chosen from typical concentrations in cancer cells, measured kinetic parameters, and measured concentrations of nutrients such as glucose, oxygen, and total intracellular adeno-nucleotide concentration.

At this stage, the model is subjected to a thorough statistical analysis. A Monte Carlo simulation is conducted for which the parameters within the model are randomized and resulting differential equations are solved. The distributions for each parameter are chosen to capture observed ranges of variation. After each simulation, numerical stability and thermodynamics are assessed and all simulations that are unstable or appear thermodynamically infeasible (i.e., a positive net flux through glycolysis is required) are rejected. In each simulation, concentrations, fluxes, metabolic control coefficients, and thermodynamic quantities are computed and recorded. This statistical analysis explores the space of glycolysis in the context of the Warburg Effect. By assessing the statistics, inferences can be made on the determinants of the Warburg Effects and its context dependence to pharmacological intervention and nutrient environment. For experimental confirmation of the model, we first utilized a recently developed a NADH/NAD⁺ fluorescent reporter (Hung et al., 2011) and measured the ratio of NADH/NAD⁺ across a population of cells where some experience hypoxia and others glucose deprivation (Figure 1C). This variation in nutrient availability across individual cells occurs in epithelial cell cultures due to differences in diffusion and available cell surface area (Sheta et al., 2001). The distribution is plotted for three media concentrations ranging from extreme hypoglycemia to the hyperglycemic conditions typically used in cell culture (0.5 mM, 5.5 mM, and 25 mM).

Next we noted that the distribution of resulting fluxes (Figure 1D) was first found consistent with known measurements. An analysis of the simulation revealed a distribution of NAD⁺ levels and NADH/NAD⁺ redox potential consistent with those observed experimentally (Figure 1E,F). ATP levels peaked around 3 mM with little variation (Figure 1G) and the ATP/ADP energy state is observed to be bimodal (Figure 1H). Concentrations along the glycolytic pathway varied with means similar to those measured in normal tissues with fructose-1,6-bisphosphate (FBP) being most variable (Figure 1I). Together, these findings indicate that these simulations capture the range of measured cellular concentrations and provide confidence for further assessment of the behavior within the model.

Evaluation of the Warburg Effect and its relationship to metabolic variables

Having developed a model of glycolysis and its regulation, we assessed the relationship of aerobic glycolysis to other characteristics of glycolysis. The distribution of values (Figure 2A) of the Warburg Effect (the ratio of flux to lactate over that entering the mitochondria) ranged from less than one (primarily oxidative metabolism) to over 95% of glucose being converted to lactate. The dynamic range over which the Warburg Effect was observed in our model allowed us to investigate to other variables in metabolism. We correlated the calculated value of the Warburg Effect with metabolite concentrations of intermediates in glycolysis. From an analysis of these correlations, a pattern within glycolysis emerges. The levels of intermediates in the beginning steps in glycolysis positively correlate with the Warburg Effect, the levels leading up to the oxidation of glyceraldehyde-3-phopshate (GAP) by GAPDH negatively correlate with the Warburg Effect, and those following GAPDH positively correlate with the Warburg Effect (Figure 2B) suggesting together that a bottleneck exists in the pathway that determines the extent of fermentation. An analysis of the correlation of enzyme expression with the Warburg Effect revealed that the enzyme expression for any single step within glycolysis did not completely correlate with the Warburg Effect. However, across glycolysis, GAPDH, the enzyme that carries out oxidative phosphorylation of glyceraldehyde-3-phosphate to yield NADH and 1,3-diphosphoglycerate most strongly correlated with aerobic glycolysis (Figure 2C). Notably many of the correlations although significant, appear not very strong indicating that the expression of individual enzymes is not sufficient for induction of aerobic glycolysis. The model also captures many of the activities that are known to correlate with aerobic glycolysis. These include glucose transport, pyruvate kinase, and lactate dehydrogenase activities. Interestingly, the expression of enzymes such as hexokinase and phosphofructokinase were either uncorrelated or negatively correlated with the extent of aerobic glycolysis likely indicating their role in creating bottlenecks at other points along the pathway. Together, these findings identify enzyme expression patterns that determine the extent of the Warburg Effect.

Figure 2

Download asset Open asset

Evaluation of the statistics of the Warburg Effect and relationships to other variables in metabolism.

(A) Probability density function (PDF) of the Warburg Effect (WE) defined as the ratio of flux through LDH to that of flux into the mitochondria. (B) Pearson correlations of intermediate metabolite levels in glycolysis with the extent of the Warburg Effect (WE). (C) Pearson correlations of the expression levels of glycolytic enzymes with the extent of the Warburg Effect (WE). (D) Pearson correlations of coupled metabolic parameters with the extent of the Warburg Effect (WE).

https://doi.org/10.7554/eLife.03342.004

We next investigated the relationship between the Warburg Effect and other physiological variables including the NADH/NAD⁺ redox status, the energy state defined as the ATP/ADP ratio, lactate, oxygen levels, and phosphocreatine levels (Figure 2D). Oxygen concentration, ATP, and phosphocreatine levels positively correlated with the Warburg Effect and NADH/NAD⁺ redox status and NADH levels in the cytosol negatively correlated with the Warburg Effect suggesting that positive and negative feedback inherent to the circuitry of glycolysis contributes to buffering the Warburg Effect. Together, these results identify multiple relationships between the extent of aerobic glycolysis and measurable variables in metabolism.

Flux control in glycolysis

After assessing how metabolic parameters and concentrations within glycolysis determine the flux to lactate, we next investigated how each node within glycolysis exerts its control on the Warburg Effect. Metabolic control analysis (MCA) provides a mathematical framework for evaluating the extent that a change in metabolic activity affects a given flux. In MCA, a pure rate-limiting step occurs when the flux control coefficient (FCC) (Supplementary Information) is one at that step and zero at all other steps. In most cases, values of FCCs are distributed with gradual values across a pathway. We carried out a metabolic control analysis using MCA within our statistical algorithm to investigate the influence of each node in glycolysis on lactate flux across the ensemble of statistical realizations of glycolysis (Figure 3A). We first computed the distribution of FCCs for lactate production for each step in the pathway (Figure 3B). It was found that for each step within glycolysis, the average control exerted was near zero. Steps early in glycolysis involving enzymes Hexokinase (HK), Phosphoglucoisomerase (PGI), and Phosphofructokinase (PFK) exhibit both positive and negative control, and GAPDH on average exhibits the most positive control on the flux through the pathway. In addition to steps within glycolysis, ATPase activity exerts the most influence over lactate production with oxygen consumption having less of an affect. Together, these results both confirm the long-standing hypothesis that ATPase activity is the most prominent rate-determining step in glycolysis but also suggest that GAPDH can often exert a large control on the flux to lactate. In addition, an analysis of the statistics indicates that depending on the context, any point along the pathway can exert large flux control on lactate production with some steps (e.g., pyruvate kinase) less likely to exert control over the Warburg Effect than others. Although the finding that pyruvate kinase typically does not exert substantial control over glycolysis may appear surprising, this finding is consistent with studies that have observed only modest changes in glycolytic flux due to changes in pyruvate kinase activity (Christofk et al., 2008; Israelsen et al., 2013).

Figure 3

Download asset Open asset

Metabolic control analysis and its relationship to metabolic variables.

(A) Schematic of workflow for global sensitivity analysis. After the model is constructed and feasible solutions obtained, each realization of glycolysis is subjected to metabolic control analysis (MCA). The resulting analysis is then subject to a statistical evaluation. (B) (left) Box plots of flux control coefficient (FCC) for lactate production for each enzymatic step in glycolysis (FCC = dlnJ_lac/dln E_i) where J_lac is the rate of pyruvate conversion to lactate, and E_i is the i^th enzyme in glycolysis for each step of glycolysis. (right) Box plots of flux control coefficient (FCC) for lactate production for Oxygen consumption (OxPhos) and ATP consumption (ATP). (C) Pearson correlations between lactate FCC values for each step in glycolysis. Heat map is colored ranging from the minimum value (green) to the maximum value (purple). (D) Pearson correlations between metabolite concentrations in glycolysis and lactate FCC values for each step in glycolysis. Heat map is colored ranging from the minimum value (green) to the maximum value (purple). (E) Pearson correlations between metabolic parameters and lactate FCC values for each step in glycolysis. Heat map is colored ranging from the minimum value (green) to the maximum value (purple). (F) Pearson correlations between ratios and lactate FCC values for each step in glycolysis.

https://doi.org/10.7554/eLife.03342.005

To further investigate the contexts in which steps along glycolysis can be limiting, we correlated the FCCs with one another and carried out a hierarchical clustering that revealed modules of co-occurring flux control (Figure 3C). It was found that enzymes residing in proximal regions of upper glycolysis tend to have positively correlated flux control. However, for lower glycolysis the flux control is uncorrelated suggesting that when a step is limiting in lower glycolysis, it can be more readily disrupted. This behavior is in contrast with that of upper glycolysis where limiting steps co-occur. Furthermore, an analysis of ATPase and Oxygen consumption activities reveals their control to be related to the control of enzymes clustered at different points in glycolysis (Figure 3C) indicating that flux control exerted by these ancillary fluxes are tied to different regions of glycolysis. Since these simulations yield for each realization of glycolysis, a set of concentrations, fluxes and flux control coefficients, relationships between flux control coefficients and measurable concentrations can be obtained. An analysis of these relationships using hierarchical clustering (Figure 3D) revealed a bi-modal relationship in glycolysis where increases in metabolites in upper glycolysis led to steps in lower glycolysis exerting flux control and increases in lower glycolysis resulted in enzymes in upper glycolysis exerting flux control. Further analysis of FCCs and metabolic parameters (Figure 3E,F) revealed additional relationships. Together, these findings yield a comprehensive map of the flux control in glycolysis and its connections to metabolic variables.

Experimental flux control in glycolysis

Surprisingly, in several instances, negative flux control is observed for several enzymes that have been previously thought to be the rate-limiting steps. Negative flux control implies that inhibiting the enzyme would actually increase the rate of flux of glucose to lactate. Given the surprising findings on flux control in glycolysis, we sought to experimentally investigate these predictions. We therefore considered systematic perturbations to glycolysis in cancer cells.

We first devised a combined flux profiling and metabolite profiling method using ¹³C isotope tracing from glucose (Figure 4A). Cells were incubated with U–¹³C glucose and media were extracted at different time points and subjected to high resolution LC-MS analysis (Liu et al., 2014). The linearity of the time course allows for an exact measurement of the flux from glucose to lactate. We next carried out direct measurements of glycolysis state and flux control by perturbing glycolysis acutely with pharmacological agents and subsequently measuring metabolite levels and the flux from glucose to lactate. Genetic manipulations such as RNA interference were not possible since the measurement of dose-dependent acute effects was necessary. In each case, the compounds considered have been reported to exhibit direct inhibition of the enzyme in question and to our knowledge do not directly inhibit other enzymes in glycolysis. Nevertheless, the general specificity of these compounds is not established and off target effects that exist are likely reduced by considering acute treatments.

Figure 4

Download asset Open asset

Experimental flux control coefficients.

(A) Schematic of experimental flux control analysis. Cells are pre-incubated with ¹³C glucose and treated with differing concentrations of inhibitors that target glycolysis at different points in the pathway. Media and intracellular metabolites are collected, subjected to (liquid chromatography high resolution mass spectrometry) LC-HRMS, and subjected to flux analysis. (B) Changes in metabolite levels observed from treatment with 3PO an inhibitor of PFK2. (C) Changes in metabolite levels observed from treatment with IA an inhibitor of GAPDH. (D) Changes in metabolite levels observed from treatment with FX11 and inhibitor of LDH. For B–D, the logarithm (log₂) of the fold change of treated to vehicle across intermediates in glycolysis is shown for each concentration of compound denoted in the figure legend. Abbreviations are the same as described in Figure 1 except that HP denotes all hexose phosphates that were measured and not distinguished in the current mass spectrometry method. (E) Lactate flux from glucose as a function PFK2 inhibition. (F) Lactate flux from glucose as a function GAPDH inhibition. (G) Lactate flux from glucose as a function LDH inhibition. For E–G, the plot on the left shows the measured glucose to lactate flux as a function of the estimated fraction of enzyme inhibited (left) inhibitor concentration (right).

https://doi.org/10.7554/eLife.03342.006

We considered inhibiting glycolysis at three separate points along the pathway that are predicted to have widely variable flux control of the pathway. In the beginning of glycolysis we used 3PO, a compound that targets PFK2 thus inhibiting the phosphofructokinase step (Schoors et al., 2014). Next, we considered iodoacetate (IA), a compound that targets GADPH (Campbell-Burk et al., 1987). Finally FX11, a compound that targets LDH was considered (Granchi et al., 2011).

In each case (Figure 4B–D) differential, complex, nonlinear responses of metabolite levels to inhibiting glycolysis were observed. The exception was treatment with IA (Figure 4C) which exhibited an expected accumulation of intermediates upstream of GAPDH and depletion of intermediates downstream of the target. We next measured directly the flux control coefficients for each compound. An analysis of the flux control (Figure 4E–G) revealed several interesting responses. Strikingly, as predicted by the model, a negative flux control was observed for targeting PFK implying that inhibiting this step in these circumstances increases glycolytic flux and Warburg Effect (Figure 4E). The largest flux control was observed with inhibition of GAPDH (Figure 4F). Inhibition of LDH as predicted also resulted in little flux control (Figure 4G). Together, these findings confirm the mechanisms of flux control of glycolysis predicted by the model and demonstrate the novel regulation of glycolytic flux that can be differentially perturbed by pharmacological compounds.

FBP status and signatures of glycolysis state

We have thus far observed experimentally the dynamic behavior of flux control in glycolysis including both positive and negative flux control in the pathway and a high variability of flux control depending on the point of inhibition in the pathway. Having established the complex relationships between glycolytic flux and susceptibility to specific targeted inhibition of the pathway, we sought to investigate whether there was any predictive capacity and new mechanisms of biochemical regulation related to these findings.

We noticed that Fructose-(1,6)-bisphosphate (FBP) levels exhibited highly dynamic and counterintuitive behavior with each drug perturbation. An analysis of metabolite levels across a series of 14 conditions in triplicate involving pharmacological perturbations of PFK2, LDH, and GAPDH at different concentrations and two separate vehicle treatments (Figure 5A) revealed large magnitude, dynamic responses in FBP levels. We next investigated the extent that FBP levels could characterize the metabolic state of glycolysis. The simulated PDF of FBP levels exhibited a bimodal distribution (Figure 5B) consisting of a state of low FBP where the concentration was in the high micromolar range. In addition, there existed a state of FBP where the concentration was several orders of magnitude high in the millimolar range. When demarcating the experiments into two groups (high and low FBP) we correlated the levels with lactate flux experimentally (Figure 5C,E) and found the results to match those observed computationally (Figure 5D,F). In addition, correlations in these two states with the remainder of the concentrations of glycolytic intermediates agreed well with experiments. Together these findings suggest in addition to the model being able to capture a diverse set of experimental metabolic conditions, the FBP status in cells was able to determine the state of glycolysis and its rate-limiting steps.

Figure 5

Download asset Open asset

FBP levels predict distinct mechanisms in glycolysis.

(A) Variation of metabolite levels across glycolysis over 14 conditions in triplicate resulting in 42 independent experiments involving cells growing in basal conditions and those with differing extents of inhibition of glycolysis from results in Figure 4. (B) Simulated distribution of FBP levels in glycolysis. (C) Correlation of lactate flux with measured glycolytic intermediates for low FBP levels. The left panel shows data for FBP and right panel reports the values of the Spearman correlation coefficients for each metabolite. (D) Simulated correlation of lactate flux with metabolite levels of glycolytic intermediates in conditions of low FBP levels. (E) Correlation of lactate flux with measured glycolytic intermediates for high FBP levels. The left panel shows data for FBP and right panel reports the values of the Spearman correlation coefficients for each metabolite. (F) Simulated correlation of lactate flux with metabolite levels of glycolytic intermediates in conditions of low FBP levels.

https://doi.org/10.7554/eLife.03342.007

A unified model of aerobic glycolysis

Together, the combined computational model, metabolite profiling, and flux analysis points to a different picture of glycolysis (Figure 6). When the levels of FBP are low, metabolite levels of glycolytic intermediates tend to be more evenly distributed across the pathway. In this circumstance, flux through the pathway is controlled largely through the initial steps in glycolysis involving hexokinase and phosphofructokinase. Bottlenecks downstream of these canonical rate-limiting steps are not affecting flux through the pathway. Under these conditions, FBP is also not allosterically activating pyruvate kinase. In contrast, when the levels of FBP are high, there is a disconnect in the relative concentration of glycolytic intermediates that is marked by a separation between upper and lower glycolysis at the GAPDH step. In this case, there is an accumulation of intermediates in upper glycolysis, most notably FBP and a depletion of intermediates downstream of GAPDH. Under these circumstances, GAPDH exerts control as the most rate-determining step in the pathway since increased activity through GAPDH will serve then to create a balance along the pathway by pulling the metabolites from upper glycolysis into lower glycolysis. Under these conditions, the beginning steps of glycolysis can exert negative control on the pathway since inhibiting them will result in even greater increases in the intermediates in upper glycolysis. This state has analogies with a recently observed state of glycolysis observed in yeast where an accumulation of FBP leads to cellular toxicity (van Heerden et al., 2014). Notably in this case, in order to balance the fluxes in higher and lower glycolysis, an imbalance in the concentrations of metabolites in upper and lower glycolysis results.

Figure 6

Download asset Open asset

A unified model of aerobic glycolysis.

A unified picture of flux control in aerobic glycolysis. (left) Under conditions where there is an accumulation of intermediates in upper glycolysis and depletion of intermediates in lower glycolysis a bottleneck exists at the step involving GAPDH. This bottleneck is due to the status of energy and redox metabolism and the thermodynamics of the pathway that together mediate the flux through GAPDH. As a result, inhibiting flux through glycolysis is most sensitive to a perturbation in GAPDH activity. (right) Under conditions where the metabolites in glycolysis are distributed more evenly with levels together being either high or low, no such bottleneck exists. Instead flux through glycolysis leading to lactate production is most determined by the canonical pacemaking steps in glycolysis involving PFK and HK. The relative levels of glycolytic intermediates are denoted by the size of the text.

https://doi.org/10.7554/eLife.03342.008

Discussion

Together, our analysis yields a comprehensive, quantitative framework for understanding glycolysis and its regulation in the context of the Warburg Effect. Historically, glycolysis is thought to have a rate-limiting step at several points in the pathway (Chance and Hess, 1956; Wu, 1965). These points correspond to positions in the pathway where large free energy differences arise including the ATP-coupled enzymes Hexokinase, Phosphofructokinase, and Pyruvate Kinase (Rose and Warms, 1966; Rapoport et al., 1976; Hue and Rider, 1987). Surprisingly, we identified a strong context-dependence with both positive and negative control in glycolysis at each of these steps. In the case of inhibiting flux through PFK, the observed negative control in certain conditions implies that inhibiting this point in the pathway can lead to increased rates of fermentation. This finding provides a possible explanation for why its efficacy may be more prevalent in stromal cells (Schoors et al., 2014).

Unexpectedly, GAPDH was found to be a recurrent rate-controlling step in aerobic glycolysis. This finding, first documented to our knowledge in parasitic bacteria feeding on high glucose (Bakker et al., 1999), is in contrast to the longstanding notion that GAPDH is not a rate-determining enzyme in glycolysis with the activity of enzymes such as hexokinase, phosphofructokinase, and pyruvate kinase thought to be more controlling. In the case of GAPDH, the bottleneck occurs due to its unique placement in the pathway where it can be regulated by ATP, NAD⁺, and the levels of glucose-derived intermediates in the pathway that affects the thermodynamics of glycolysis. There are multiple mechanisms that lead to this finding. ATP consumption has previously been reported to control the rate of glycolysis and this effect likely occurs to some extent through GAPDH (Racker, 1976; Locasale and Cantley, 2011; Lunt and Vander Heiden, 2011). NAD⁺ regeneration that is mediated by the malate–aspartate shuttle and lactate dehydrogenase also affects flux through glycolysis (DeBerardinis et al., 2008; Locasale and Cantley, 2011). In addition, the activity of the pathway upstream and downstream of GAPDH changes the balance of the levels of intermediates in glycolysis and results in driving the thermodynamics of the reactions out of equilibrium also can result in greater flux control (Noor et al., 2014). Each of these mechanisms separately or together acts to allow for GAPDH to exert flux control over the glycolytic pathway.

Enzymes along glycolysis that are believed to control flux have many documented regulatory mechanisms. For example pyruvate kinase and phosphofructokinase have numerous small molecule effectors and post-translational modifications that affect their activities (Mor et al., 2011; Chaneton and Gottlieb, 2012; Keller et al., 2012; Yi, et al., 2012). It is notable that GAPDH also is subject to multiple forms of regulation including post-translational modifications such as nitrosylation and reactive oxygen species (ROS) that interacts with the catalytic cysteine in GAPDH to inhibit its activity (Gaupels et al., 2011; Tristan et al., 2011; Moellering and Cravatt, 2013). In the context of ROS, it is tempting to speculate that alterations in ROS could lead to selective modulation of glycolytic flux as has been suggested to occur with pyruvate kinase (Anastasiou et al., 2011). While ROS-mediated inhibition is unlikely in conditions of exponential growth, it may be more apparent in physiological conditions of hypoxia and glucose deprivation with higher concentrations of ROS. Another critical aspect of the flux control that GAPDH exerts over the glycolytic pathway is high expression of GAPDH in cells undergoing aerobic glycolysis. Indeed, reports of quantitative protein abundance in mammalian cells have identified enzymes in the pathway glycolysis as the most highly expressed collective of proteins in cells (Moghaddas Gholami et al., 2013). Interestingly, it was found that within glycolysis, GAPDH is often the mostly highly concentrated protein in glycolysis suggesting that the role of this high expression in these cases is to support the increased amount of glycolytic flux in these cells. Nevertheless, although mechanisms that regulate phosphofructokinase and pyruvate kinase for example have been shown to mediate cell growth and proliferation, whether regulatory mechanisms of GAPDH have functional roles in cell growth and proliferation related to aerobic glycolysis is not known.

Two states of glycolysis were observed with different extents of flux control, tendencies for aerobic glycolysis, and concentration patterns along glycolysis. Notably, the FBP levels have implications on the activity of pyruvate kinase that is also allosterically activated by FBP. At low FBP concentration, pyruvate kinase is not activated by FBP but this occurs only in the high FBP state. This finding could have implications in understanding the contexts in which pharmacologically activating pyruvate kinase may have efficacy. Furthermore, it remains to be seen if any signature of these metabolite states could manifest in the alterations of peripheral metabolism involving pathways whose fluxes emanate from glycolysis. If this were the case, then a tempting possibility would be that these states of glycolysis could be predicted from measurements of peripheral metabolites that could be excreted into circulation allowing for the possibility of developing serum biomarkers for the status of glycolysis in tumors beyond what can be resolved with positron emission tomography using radioactive glucose.

Finally, the surge of interest in metabolism and its contribution to pathogenesis has created an expectation that therapeutics that target glucose metabolism will be clinically successful. Targeting glycolysis in metabolism has raised interest but is limited by the development of biomarkers that could determine the contexts in which targeting glucose metabolism in malignancy would be efficacious (Vander Heiden, 2011; Galluzzi et al., 2013; Vander Heiden, 2013). From our model analysis, the consequences of inhibiting glycolysis appear enormously complex that limit biomarker development due to the nonlinear mechanisms that determine the response of the pathway to a drug perturbation. Nevertheless, with the development of these predictive models that can capture diverse behaviors of glycolysis, it is worth considering whether they may have predictive capacity in pre-clinical and clinical settings.

Materials and methods

Kinetic model of glycolysis

Request a detailed protocol

The model includes a compartment involving enzymes in glycolysis and additional compartments involving with reactions that are coupled to glycolysis. The following compartments are considered in addition to the enzymes in glycolysis are considered:

Glucose uptake through the glucose transporter,
Lactate dehydrogenase (LDH) activity and lactate transport through the monocarboxylate transporter (MCT),
Oxygen consumption, transport and activity of oxidative phosphorylation (OxPhos),
ATP-buffering mechanisms involving adenylate Kinase, creatine Kinase, and ATPase activity,
NADH/NAD⁺-mediated mitochondrial shuttles including the Malate–Aspartate Shuttle and the Glycerol-3-Phosphate Shuttle.

More generally glycogen metabolism, pentose phosphate pathway and alanine biosynthesis could also be included but are omitted since the conclusions drawn in this study do not depend on their activities. A schematic representation of the kinetic model and resulting network is shown in Figure 2A.

We first introduce notation and other conventions. Each symbol indicates the respective concentration. A subscript ‘0’ refers to the steady-state value. Reaction rates are expressed in millimoles per hour per unit intracellular volume. Together with initial state vectors $C_{0} \in R^{N}, V_{o} \in R^{M}$ , of metabolites and fluxes respectively.

The dynamics of the system are therefore formulated as initial value problem for ordinary differential equations (ODEs). Starting parameter values based on published data are defined in Supplementary file 1. Reaction rates are defined through a consideration of the respective enzyme mechanisms with additional feed forward and feed back regulation giving rise to allosteric activation and inhibition.

Therefore, the model is designed to describe: (a) steady-state and dynamic behavior of energy metabolism, including the Warburg effect ( $W = \frac{J_{L a c}}{J_{O x}}$ ), energy state (ATP/ADP ratio), redox state (NADH/NAD⁺ ratio) together with glucose, lactate, and oxygen supply through exchange fluxes from the intercellular and extracellular compartment, (b) steady states and time courses of all variables of the metabolic network, and (c) effects of metabolic parameter perturbations on the overall system output. Cellular energy homeostasis through ATP is supported by several mechanisms with different relaxation times and regulation mechanisms: (a) directly by glycolysis that converts glucose into intracellular pyruvate, (b) mitochondrial respiration through consumption of pyruvate and oxygen via the TCA cycle, and (c) the buffering effect of creatine kinase that facilitates the reversible reaction of phosphocreatine (PCr) with ADP to produce creatine (Cr) and ATP, and (d) adenylate kinase activity that catalyzes interconversion of adenine nucleotides: 2ADP ⇔ ATP + AMP, and (e) the additional NADH pool produced in cytosol and transported to mitochondria by different shuttle mechanisms. For oxygen consumption, we assumed that mitochondrial respiration (O₂ consumption) depends on cellular pyruvate and oxygen concentrations and that respiratory chain activity is activated by ADP concentration.

The state of the system is represented by the state vector of time-dependent metabolite concentrations C_n, (n = 1, …, N) and includes 27 state variables whose names, balance equations, and steady-state values are presented in Supplementary file 1. The model also includes three mass conservations laws, that reduce the number of state variables. The temporal profile of the system is governed by the set of ODE based on the model network that distinguishes media and cellular domains. For the metabolite i in the intracellular compartment, the general form is:

\frac{d C_{i}}{d t} = J_{i}^{t r} + \sum_{j} ν_{i}^{p j} V_{i}^{p j} - \sum_{k} ν_{i}^{u k} V_{i}^{u k},

where the intracellular concentration of species i is C_i; $J_{i}^{t r}$ is the net transport flux for boundary species i between the media and cell; $V_{i}^{p j}$ , $V_{i}^{u k}$ are normalized reaction fluxes that produce (j) or utilize (k) cellular species i, $ν_{i}^{p j}$ , and $ν_{i}^{u k}$ are corresponding stoichiometric coefficients. For extracellular boundary metabolites in the media, the dynamic mass balance has the following general form:

\frac{d C e}{d t} = r_{i e} \cdot J_{i}^{t r},

where the concentration of extracellular (e) species is C_e and r_ie is the ratio of cell volume to media volume. To ensure that steady states are obtained during perturbations of parameters without loss of generality, we assume that the media boundary metabolites have constant concentrations. The transport of boundary species is taken to be either facilitated (glucose, lactate, serine, glycine) or passive (oxygen). For passive diffusion of oxygen the equation for the net transport flux is:

J_{O_{2}}^{t r} = k_{O_{2}} (C_{O 2 e} - C_{O 2 i}),

where $k_{O_{2}}$ is the effective rate constant for passive O₂ diffusion, and $C_{e O 2}$ and $C_{i O 2}$ are medium and intracellular oxygen concentrations. For facilitated transport, the equation for the net transport flux is:

J_{i}^{t r} = \frac{J_{t r}^{m a x} (\frac{C_{e i}}{K_{m i}^{t r}} - \frac{C_{i}}{K_{m i}^{t r}})}{1 + \frac{C_{e i}}{K_{m i}^{t r}} + \frac{C_{i}}{K_{m i}^{t r}}},

where C_ei is the extracellular concentration for species i, $J_{t r}^{m a x}$ is the maximal transport rate, and $K_{m i}^{t r}$ is the Michaelis–Menten constant for transport. NADH transport from the cytosol to the mitochondria is mediated by the malate–aspartate and glycerol phosphate shuttles. We modeled the NAD⁺-mediated mitochondrial shuttle flux from NADH to NAD⁺ by assuming that the total shuttle flux is balanced with the NADH-generation reactions involving LDH and PHGDH.

V_{m a s} = V_{o x} + 2 V_{p h g d h}

When possible, in order to minimize the number of parameters, we utilized a ‘one-step binding enzyme mechanism’ with a rate law of the form (Segel, 1975):

V = (\frac{V_{f}^{m a x} \frac{\prod_{i} C_{s i}^{ν_{i}}}{K m f} - V_{r}^{m a x} \frac{\prod_{j} C_{p j}^{ν_{j}}}{K m r}}{1 + \frac{\prod_{i} C_{s i}^{ν_{i}}}{K m f} + \frac{\prod_{j} C_{p j}^{ν_{j}}}{K m r}}) χ (α, ι),

where $χ (α, ι)$ is a control function that accounts for the effects of activation $α$ , or inhibition $ι$ . A general form for the control function is (Liebermeister and Klipp, 2006):

χ (α, ι) = {\begin{matrix} \frac{[A]}{([A] + K_{A})} i n c a s e o f a c t i v a t i o n, α \\ \frac{K_{I}}{([I] + K_{I})} i n c a s e o f i n h i b i t i o n, ι \end{matrix},

where $[A]$ and $[I]$ are the concentrations of the activator or inhibitor, K_A and K_I are the corresponding activation and inhibition constants. With this convention, there are 4 (or 5 including the control function) independent parameters for each metabolic flux with a one-step binding mechanism. For fluxes with no allosteric regulation, $χ (α, ι) = 1$ .

For all metabolic fluxes involving NAD⁺/NADH oxidoreductase activity, we utilized a more complex rate law than that of a one step binding reaction (Cornish-Bowden, 2012). We chose to consider a more complicated rate-law since our initial observations indicated that a one-step binding mechanism produced numerical instabilities in the solutions. We therefore considered a general rate law corresponding to a random bi–bi mechanism (for LDH and PHGDH) and a random ter-bi reaction for GAPDH. For example, in its general form the rate law for LDH is:

V_{l d h} = \frac{V_{f l d h}^{m a x} \frac{P Y R * N A D H}{K m l d h f} - V_{r l d h}^{m a x} \frac{L A C_{i} * N A D^{+}}{K m l d h r}}{1 + \frac{P Y R}{K_{m p y r}^{l d h}} + \frac{N A D H}{K_{m n a d h}^{l d h}} + \frac{P Y R * N A D H}{K m l d h f} + \frac{L A C_{i}}{K_{m l a c}^{l d h}} + \frac{N A D^{+}}{K_{m n a d}^{l d h}} + \frac{L A C_{i} * N A D^{+}}{K m l d h r}} .

However, it has also been suggested that the enzymatic reaction for oxidoreductase activity proceeds in an ordered fashion resulting in a rate law of the form:

V_{l d h} = \frac{V_{f l d h}^{m a x} \frac{P Y R * N A D H}{K m l d h f} - V_{r l d h}^{m a x} \frac{L A C_{i} * N A D^{+}}{K m l d h r}}{1 + \frac{N A D H}{K_{m n a d h}^{l d h}} + \frac{P Y R * N A D H}{K m l d h f} + \frac{N A D^{+}}{K_{m n a d}^{l d h}} + \frac{L A C_{i} * N A D^{+}}{K m l d h r}}

When considering this mechanism, it was found that statistical evaluation of the model was robust to changes in the choice of rate law.

The rate equations for each metabolic reaction and metabolic parameters values are in Supplementary file 1. Taking into account additional thermodynamic constraints by invoking the Haldane equation that relates the forward and reverse fluxes and their Vmax values allowed for a further reduction in the number of independent parameters to three parameters per reaction flux. The Haldane equation has the form:

V_{r}^{m a x} = \frac{V_{f}^{max} K m r}{K m f K_{e q}},

where K_eq is the reaction apparent equilibrium constant: $K_{e q} = \exp (- Δ G^{°} / R T)$ , where $Δ G^{°}$ is standard reaction Gibbs free energy. $Δ G^{°}$ values for glycolytic reactions are listed in Supplementary file 1.

We first used two independent sets of experimental data from a compendium of steady-state metabolite concentrations of liver cells and steady-state fluxes of DB1 melanoma cells evaluated using a cultured bioreactor (Konig et al., 2012; Shestov et al., 2013). The model was further validated by a comparison of measurements of redox status, assessment of ATP, ADP, and AMP concentrations. Together, those data were able to generate a thermodynamically feasible and experimentally observable model of glycolysis that was then used as the starting point for exhaustive Monte Carlo simulations and coupled metabolic control analysis.

Monte Carlo simulations for global sensitivity analysis

Request a detailed protocol

The specific aim of this work is to enumerate the variables within and coupled to glycolysis that determine the extent of aerobic glycolysis. We therefore developed a novel algorithm to assess these features. The algorithm involves a global sensitivity analysis based on a Monte Carlo method. The Monte Carlo analysis allows for sampling of parameter space over a broad range of simultaneous variations of parameters (enzyme expression levels) followed by statistical assessment of the resulting solution. The idea of this method is to inject uncertainty of the parameters in the model by randomly selecting parameter values from uniform probability distributions. This was achieved by a Monte Carlo method using a randomly drawn set of parameters. The range of parameters chosen with uniform sampling was within two orders of magnitudes of each Vmax (which is proportional to enzyme activity level and was in the range of 10 times less and 10 times greater than the initial Vmax value). Such a sampling was chosen to be large enough to cover all feasible solutions of aerobic glycolysis. The sampling was carried out using a Latin Hypercube sampling method (Oguz et al., 2013). For each chosen set of random parameter vectors the model was simulated between 2000 and 5000 times and its output was calculated. Briefly, we divided the range of the i-th normalized parameter into n (n = 2000–5000) subintervals of equal size. Then, we randomly sample n values (e.g., p_i, i = 1, …, n), one from each subinterval, for the i-th parameter. We next randomly permuted the n values for each parameter to get the parameter vector. We then evaluated the following: the Warburg Effect value, energy and redox states, metabolite steady-state concentrations and fluxes and flux control coefficients (FCCs). The flux control coefficient that an enzyme E_i exerts on the lactate flux J_Lac is defined as (Heinrich and Rapoport, 1974; Fell, 1992; Shestov et al., 2013):

F C C_{i} = \frac{\partial \ln J_{L a c}}{\partial \ln E_{i}} .

To compute the FCC values numerically, we considered a perturbation 0.01 times of the enzyme value and then evaluated the change in flux. The model was implemented in Matlab as a system of 27 ordinary differential equations (ODEs) with an ODE solver designed for stiff ODE systems.

Cell culture and metabolite extraction

Request a detailed protocol

HCT116 cells were cultured with a growth medium, which contains RPMI 1640, 10% Fetal Bovine Serum (FBS), 100 U/ml penicillin, and 100 μg/ml streptomycin. The cells were obtained as a gift from Lewis Cantley's laboratory. The cells were cultured in 37°C with 5% CO₂. For treatments, the cells were seeded in 6-well plate at a density of 2 × 10⁵ to 5 × 10⁵ cells per well. After overnight incubation, full growth media were removed, and cells were washed with 2 ml PBS before the addition of RPMI media (without glucose), supplemented with 10% dialyzed and heat inactivated FBS, 100 U/ml penicillin, 100 μg/ml streptomycin and 5 mM 13C-U-glucose (Cambridge Isotope Laboratory, Tewksbury, MA) (fresh growth media). For drug treatments and flux measurements, full growth media were replaced with fresh growth media, containing either 0.1% DMSO, or the concentration of the indicated drug. Media were then collected at 30, 60, 90, and 120 min after incubation. From 20 µl media, 80 µl of ice cold H₂O was added, together with 400 µl ice cold methanol (Fisher, Optima LC/MS grade). After vigorous vortexing, the solution was then centrifuged at 20,000×g at 4°C for 10 min and then the supernatant was dried under vacuum. At 120 min, the media were removed as completely as possible, and then 6-well plates were immediately placed on dry ice, followed by the addition of 1 ml extraction solvent, 80% MeOH/H₂O (Fisher, Optima LC/MS Grade), which was pre-cooled in −80°C freezer for at least 1 hr. The dishes were then transferred to the −80°C freezer. The plates were left for 15 min and then cells were scraped into the solvent on dry ice and then transferred to two 1.7-ml eppendorf tubes, and centrifuged with the speed of 20,000×g at 4°C for 10 min. Metabolite extracts are prepared from three separate wells to make three replicate samples. The supernatant is then transferred to new eppendorf tubes, and dried under vacuum. For analysis, each extract was then re-constituted into water (15 μl for cell extract and 50 μl for medium extract) and then diluted with an equal volume of 50% Methanol/Acetonitrile. Finally, 5 μl was injected into the column for analysis.

Flux analysis

Request a detailed protocol

Absolute concentrations of ¹³C lactate were first measured by mixing a standard of unlabeled lactate with a defined concentration into the medium. For each calculated flux, an estimation of the flux from glucose to lactate was obtained from the slope of the time course of ¹³C lactate production using the time points of media collection described above. Slopes were computed using Graphpad Prism with time courses determined to be linear from the goodness of fit (R² > 0.98 with few exceptions). For corresponding plots for different agents, fraction of inhibition was computed from K_i values that were taken based on previous reports (Foxall et al., 1984; Le et al., 2010; Telang et al., 2012).

Mass spectrometry

Request a detailed protocol

The Q Exactive Mass Spectrometer (QE-MS) is equipped with a heated electrospray ionization probe (HESI), and the relevant parameters are as listed: heater temperature, 120°C; sheath gas, 30; auxiliary gas, 10; sweep gas, 3; spray voltage, 3.6 kV for positive mode, and 2.5 kV for negative mode. Capillary temperature was set at 320°C, and S-lens was 55. A full scan range from 60 to 900 (m/z) was used. The resolution was set at 70,000. The maximum injection time was 200 ms with typical injection times around 50 ms. These settings resulted in a duty cycle of around 550 ms to carry out scans in both positive and negative mode. Automated gain control (AGC) was targeted at 3,000,000 ions.

Liquid chromatography

Request a detailed protocol

Liquid chromatography (Ultimate 3000 UHPLC) is coupled to the QE-MS for metabolite separation and detection. An Xbridge amide column (100 × 2.1 mm i.d., 3.5 μm; Waters) is employed for compound separation. The mobile phase A is 20 mM ammonium acetate and 15 mM ammonium hydroxide in water with 3% acetonitrile, pH 9.0, as described above, and mobile phase B is acetonitrile. A linear gradient was used as follows: 0 min, 85% B; 1.5 min, 85% B; 5.5 min, 35% B; 10 min, 35% B; 10.5 min, 35% B; 14.5 min, 35% B; 15 min, 85% B; and 20 min, 85% B. The flow rate was 0.15 ml/min from 0 to 10 min and 15 to 20 min, and 0.3 ml/min from 10.5 to 14.5 min. All solvents are LC-MS Optima grade and purchased from Fisher Scientific.

NADH/NAD⁺ imaging

Request a detailed protocol

Human mammary epithelial MCF-10A cells (CRL-10317; ATCC) stably expressing Peredox-NLS were generated and cultured as previously described (Debnath et al., 2003; Hung et al., 2011). 2–4 days prior to imaging, ∼500 cells were plated onto the center of each well of a 96-well plate. On the day of the experiment, cells were placed in custom DMEM/F12 (Gibco) containing no glucose and supplemented to the levels mentioned. Fluorescence images were acquired using a Nikon inverted Eclipse Ti microscope, equipped with a Nikon 20×/0.75 Plan Apo objective, three different regions of interest chosen, images sequentially acquired every 8 min with 50–100 ms exposure and 2 × 2 binning. Using a custom MATLAB algorithm previously developed, we subtracted background, set a threshold for cell segmentation, and analyzed the data as previously described (Hung et al., 2011). The concentration of glucose is maintained by considering a culture system in which cells are seeded at low density and media are present in vast excess. An estimate of glucose maintenance in the media in the low glucose condition can be considered using values the glucose uptake rate of the cells (∼100 fmol/cell/hr), the cell number (∼2000), the volume of media used (400 μl), and concentration of glucose in the media (750 μM, 500 μM from supplementation + 250 μM from the Horse Serum). Thus in 24 hr, we estimate that about 5 nmols of glucose are consumed. The media contain roughly 200 nmols of glucose which is in excess of the amount consumed by the cells.

Data analysis

Request a detailed protocol

Raw data collected from the QE-MS is processed using Sieve 2.0 (Thermo Scientific). Peak alignment and detection are performed according to the protocol described by Thermo Scientific. For a targeted metabolite analysis, the method ‘peak alignment and frame extraction’ is applied. An input file of theoretical m/z and detected retention time known metabolites is used for targeted metabolite analysis with data collected in both positive and negative mode. m/z width is set at 10 ppm. The output file including detected m/z and relative intensity in different samples is obtained after data processing. For resulting simulation data, hierarchical clustering was carried out using spearman ranked correlations and the Gene-e software package (Broad Institute). Box-plots (25%/75% percentile, mean and median) was calculated and made with the Graphpad Prism software package.

References

1. Anastasiou D
2. Poulogiannis G
3. Asara JM
4. Boxer MB
5. Jiang JK
6. Shen M
7. Bellinger G
8. Sasaki AT
9. Locasale JW
10. Auld DS
11. Thomas CJ
12. Vander Heiden MG
13. Cantley LC
(2011) Inhibition of pyruvate kinase M2 by reactive oxygen species contributes to cellular antioxidant responses
Science 334:1278–1283.

https://doi.org/10.1126/science.1211485
- Google Scholar
(1999) What controls glycolysis in bloodstream form Trypanosoma brucei?
The Journal of Biological Chemistry 274:14551–14559.

https://doi.org/10.1074/jbc.274.21.14551
- Google Scholar
(1987) 31P NMR saturation-transfer measurements in Saccharomyces cerevisiae: characterization of phosphate exchange reactions by iodoacetate and antimycin A inhibition
Biochemistry 26:7483–7492.

https://doi.org/10.1021/bi00397a043
- Google Scholar
(2002) Metabolic control analysis in drug discovery and disease
Nature Biotechnology 20:243–249.

https://doi.org/10.1038/nbt0302-243
- Google Scholar
1. Chance B
2. Hess B
(1956) On the control of metabolism in ascites tumor cell suspensions
Annals of the New York Academy of Sciences 63:1008–1016.

https://doi.org/10.1111/j.1749-6632.1956.tb50908.x
- Google Scholar
1. Chaneton B
2. Gottlieb E
(2012) Rocking cell metabolism: revised functions of the key glycolytic regulator PKM2 in cancer
Trends in Biochemical Sciences 37:309–316.

https://doi.org/10.1016/j.tibs.2012.04.003
- Google Scholar
(2008) The M2 splice isoform of pyruvate kinase is important for cancer metabolism and tumour growth
Nature 452:230–233.

https://doi.org/10.1038/nature06734
- Google Scholar
Book
1. Cornish-Bowden A
(2012)
Fundamentals of enzyme kinetics

Weinheim: Wiley-Blackwell.
- Google Scholar
(2008) Sloppiness, robustness, and evolvability in systems biology
Current Opinion in Biotechnology 19:389–395.

https://doi.org/10.1016/j.copbio.2008.06.008
- Google Scholar
(2008) The biology of cancer: metabolic reprogramming fuels cell growth and proliferation
Cell Metabolism 7:11–20.

https://doi.org/10.1016/j.cmet.2007.10.002
- Google Scholar
(2003) Morphogenesis and oncogenesis of MCF-10A mammary epithelial acini grown in three-dimensional basement membrane cultures
Methods 30:256–268.

https://doi.org/10.1016/S1046-2023(03)00032-X
- Google Scholar
1. Fell DA
(1992)
Metabolic control analysis: a survey of its theoretical and experimental development

The Biochemical Journal 286:313–330.
- Google Scholar
1. Folger O
2. Jerby L
3. Frezza C
4. Gottlieb E
5. Ruppin E
6. Shlomi T
(2011) Predicting selective drug targets in cancer through metabolic networks
Molecular Systems Biology 7:501.

https://doi.org/10.1038/msb.2011.35
- Google Scholar
(1984) The inhibition of erythrocyte glyceraldehyde-3-phosphate dehydrogenase. In situ PMR studies
Biochimica et Biophysica Acta 804:209–215.

https://doi.org/10.1016/0167-4889(84)90151-4
- Google Scholar
(2013) Metabolic targets for cancer therapy
Nature Reviews Drug Discovery 12:829–846.

https://doi.org/10.1038/nrd4145
- Google Scholar
(2011) Upstream and downstream signals of nitric oxide in pathogen defence
Current Opinion in Plant Biology 14:707–714.

https://doi.org/10.1016/j.pbi.2011.07.005
- Google Scholar
1. Granchi C
2. Roy S
3. Giacomelli C
4. Macchia M
5. Tuccinardi T
6. Martinelli A
7. Lanza M
8. Betti L
9. Giannaccini G
10. Lucacchini A
11. Funel N
12. León LG
13. Giovannetti E
14. Peters GJ
15. Palchaudhuri R
16. Calvaresi EC
17. Hergenrother PJ
18. Minutolo F
(2011) Discovery of N-hydroxyindole-based inhibitors of human lactate dehydrogenase isoform A (LDH-A) as starvation agents against cancer cells
Journal of Medicinal Chemistry 54:1599–1612.

https://doi.org/10.1021/jm101007q
- Google Scholar
1. Hamanaka RB
2. Chandel NS
(2012) Targeting glucose metabolism for cancer therapy
The Journal of Experimental Medicine 209:211–215.

https://doi.org/10.1084/jem.20120162
- Google Scholar
1. Heinrich R
2. Rapoport TA
(1974) A linear steady-state treatment of enzymatic chains. General properties, control and effector strength
European Journal of Biochemistry 42:89–95.

https://doi.org/10.1111/j.1432-1033.1974.tb03318.x
- Google Scholar
1. Hue L
2. Rider MH
(1987)
Role of fructose 2,6-bisphosphate in the control of glycolysis in mammalian tissues

The Biochemical Journal 245:313–324.
- Google Scholar
1. Hung YP
2. Albeck JG
3. Tantama M
4. Yellen G
(2011) Imaging cytosolic NADH-NAD(+) redox state with a genetically encoded fluorescent biosensor
Cell Metabolism 14:545–554.

https://doi.org/10.1016/j.cmet.2011.08.012
- Google Scholar
1. Israelsen WJ
2. Dayton TL
3. Davidson SM
4. Fiske BP
5. Hosios AM
6. Bellinger G
7. Li J
8. Yu Y
9. Sasaki M
10. Horner JW
11. Burga LN
12. Xie J
13. Jurczak MJ
14. DePinho RA
15. Clish CB
16. Jacks T
17. Kibbey RG
18. Wulf GM
19. Di Vizio D
20. Mills GB
21. Cantley LC
22. Vander Heiden MG
(2013) PKM2 isoform-specific deletion reveals a differential requirement for pyruvate kinase in tumor cells
Cell 155:397–409.

https://doi.org/10.1016/j.cell.2013.09.025
- Google Scholar
1. Keller KE
2. Tan IS
3. Lee YS
(2012) SAICAR stimulates pyruvate kinase isoform M2 and promotes cancer cell survival in glucose-limited conditions
Science 338:1069–1072.

https://doi.org/10.1126/science.1224409
- Google Scholar
(2012) Quantifying the contribution of the liver to glucose homeostasis: a detailed kinetic model of human hepatic glucose metabolism
PLOS Computational Biology 8:e1002577.

https://doi.org/10.1371/journal.pcbi.1002577
- Google Scholar
(2011) Otto Warburg's contributions to current concepts of cancer metabolism
Nature Reviews Cancer 11:325–337.

https://doi.org/10.1038/nrc3038
- Google Scholar
1. Le A
2. Cooper CR
3. Gouw AM
4. Dinavahi R
5. Maitra A
6. Deck LM
7. Royer RE
8. Vander Jagt DL
9. Semenza GL
10. Dang CV
(2010) Inhibition of lactate dehydrogenase A induces oxidative stress and inhibits tumor progression
Proceedings of the National Academy of Sciences of the United States of America 107:2037–2042.

https://doi.org/10.1073/pnas.0914433107
- Google Scholar
1. Liebermeister W
2. Klipp E
(2006) Bringing metabolic networks to life: integration of kinetic, metabolic, and proteomic data
Theoretical Biology & Medical Modelling 3:42.

https://doi.org/10.1186/1742-4682-3-42
- Google Scholar
1. Liu X
2. Ser Z
3. Locasale JW
(2014) Development and quantitative evaluation of a high-resolution metabolomics technology
Analytical Chemistry 86:2175–2184.

https://doi.org/10.1021/ac403845u
- Google Scholar
1. Locasale JW
2. Cantley LC
(2011) Metabolic flux and the regulation of mammalian cell growth
Cell Metabolism 14:443–451.

https://doi.org/10.1016/j.cmet.2011.07.014
- Google Scholar
1. Lunt SY
2. Vander Heiden MG
(2011) Aerobic glycolysis: meeting the metabolic requirements of cell proliferation
Annual Review of Cell and Developmental Biology 27:441–464.

https://doi.org/10.1146/annurev-cellbio-092910-154237
- Google Scholar
(2012) Intra-tumour heterogeneity: a looking glass for cancer?
Nature Reviews Cancer 12:323–334.

https://doi.org/10.1038/nrc3261
- Google Scholar
1. Moellering RE
2. Cravatt BF
(2013) Functional lysine modification by an intrinsically reactive primary glycolytic metabolite
Science 341:549–553.

https://doi.org/10.1126/science.1238327
- Google Scholar
1. Moghaddas Gholami A
2. Hahne H
3. Wu Z
4. Auer FJ
5. Meng C
6. Wilhelm M
7. Kuster B
(2013) Global proteome analysis of the NCI-60 cell line panel
Cell Reports 4:609–620.

https://doi.org/10.1016/j.celrep.2013.07.018
- Google Scholar
(2009) Shifts in growth strategies reflect tradeoffs in cellular economics
Molecular Systems Biology 5:323.

https://doi.org/10.1038/msb.2009.82
- Google Scholar
(2011) Control of glycolysis through regulation of PFK1: old friends and recent additions
Cold Spring Harbor Symposia on Quantitative Biology 76:211–216.

https://doi.org/10.1101/sqb.2011.76.010868
- Google Scholar
1. Noor E
2. Bar-Even A
3. Flamholz A
4. Reznik E
5. Liebermeister W
6. Milo R
(2014) Pathway thermodynamics highlights kinetic obstacles in central metabolism
PLOS Computational Biology 10:e1003483.

https://doi.org/10.1371/journal.pcbi.1003483
- Google Scholar
1. Oguz C
2. Laomettachit T
3. Chen KC
4. Watson LT
5. Baumann WT
6. Tyson JJ
(2013) Optimization and model reduction in the high dimensional parameter space of a budding yeast cell cycle model
BMC Systems Biology 7:53.

https://doi.org/10.1186/1752-0509-7-53
- Google Scholar
1. Racker E
(1976) Why do tumor cells have a high aerobic glycolysis?
Journal of Cellular Physiology 89:697–700.

https://doi.org/10.1002/jcp.1040890429
- Google Scholar
(1976)
The regulatory principles of glycolysis in erythrocytes in vivo and in vitro. A minimal comprehensive model describing steady states, quasi-steady states and time-dependent processes

The Biochemical Journal 154:449–469.
- Google Scholar
1. Rose IA
2. Warms JV
(1966)
Control of glycolysis in the human red blood cell

The Journal of Biological Chemistry 241:4848–4854.
- Google Scholar
(1999) Metabolic pathway analysis: basic concepts and scientific applications in the post-genomic era
Biotechnology Progress 15:296–303.

https://doi.org/10.1021/bp990048k
- Google Scholar
1. Schoors S
2. De Bock K
3. Cantelmo AR
4. Georgiadou M
5. Ghesquière B
6. Cauwenberghs S
7. Kuchnio A
8. Wong BW
9. Quaegebeur A
10. Goveia J
11. Bifari F
12. Wang X
13. Blanco R
14. Tembuyser B
15. Cornelissen I
16. Bouché A
17. Vinckier S
18. Diaz-Moralli S
19. Gerhardt H
20. Telang S
21. Cascante M
22. Chesney J
23. Dewerchin M
24. Carmeliet P
(2014) Partial and transient reduction of glycolysis by PFKFB3 blockade reduces pathological angiogenesis
Cell Metabolism 19:37–48.

https://doi.org/10.1016/j.cmet.2013.11.008
- Google Scholar
Book
1. Segel IH
(1975)
Enzyme kinetics: behavior and analysis of rapid equilibrium and steady state enzyme systems

New York: Wiley.
- Google Scholar
1. Shestov AA
2. Barker B
3. Gu Z
4. Locasale JW
(2013) Computational approaches for understanding energy metabolism
Wiley Interdisciplinary Reviews System Biology and Medicie.

https://doi.org/10.1002/wsbm.1238
- Google Scholar
(2013) Metabolic network analysis of DB1 melanoma cells: how much energy is derived from aerobic glycolysis?
Advances in Experimental Medicine and Biology 765:265–271.

https://doi.org/10.1007/978-1-4614-4989-8_37
- Google Scholar
(2001) Cell density mediated pericellular hypoxia leads to induction of HIF-1alpha via nitric oxide and Ras/MAP kinase mediated signaling pathways
Oncogene 20:7624–7634.

https://doi.org/10.1038/sj.onc.1204972
- Google Scholar
1. Shlomi T
2. Benyamini T
3. Gottlieb E
4. Sharan R
5. Ruppin E
(2011) Genome-scale metabolic modeling elucidates the role of proliferative adaptation in causing the Warburg effect
PLOS Computational Biology 7:e1002018.

https://doi.org/10.1371/journal.pcbi.1002018
- Google Scholar
1. Telang S
2. Clem BF
3. Klarer AC
4. Clem AL
5. Trent JO
6. Bucala R
7. Chesney J
(2012) Small molecule inhibition of 6-phosphofructo-2-kinase suppresses t cell activation
Journal of Translational Medicine 10:95.

https://doi.org/10.1186/1479-5876-10-95
- Google Scholar
1. Tristan C
2. Shahani N
3. Sedlak TW
4. Sawa A
(2011) The diverse functions of GAPDH: views from different subcellular compartments
Cellular Signalling 23:317–323.

https://doi.org/10.1016/j.cellsig.2010.08.003
- Google Scholar
1. Vander Heiden MG
(2011) Targeting cancer metabolism: a therapeutic window opens
Nature Reviews Drug Discovery 10:671–684.

https://doi.org/10.1038/nrd3504
- Google Scholar
1. Vander Heiden MG
(2013) Exploiting tumor metabolism: challenges for clinical translation
The Journal of Clinical Investigation 123:3648–3651.

https://doi.org/10.1172/JCI72391
- Google Scholar
1. van Heerden JH
2. Wortel MT
3. Bruggeman FJ
4. Heijnen JJ
5. Bollen YJ
6. Planqué R
7. Hulshof J
8. O'Toole TG
9. Wahl SA
10. Teusink B
(2014) Lost in transition: start-up of glycolysis yields subpopulations of nongrowing cells
Science 343:1245114.

https://doi.org/10.1126/science.1245114
- Google Scholar
1. Vazquez A
2. Liu J
3. Zhou Y
4. Oltvai ZN
(2010) Catabolic efficiency of aerobic glycolysis: the Warburg effect revisited
BMC Systems Biology 4:58.

https://doi.org/10.1186/1752-0509-4-58
- Google Scholar
(1927) The metabolism of tumors in the Body
The Journal of General Physiology 8:519–530.

https://doi.org/10.1085/jgp.8.6.519
- Google Scholar
1. Wu R
(1965)
Rate-Limiting factors in glycolysis and inorganic orthophosphate transport in rat liver and kidney slices

The Journal of Biological Chemistry 240:2373–2381.
- Google Scholar
1. Yi W
2. Clark PM
3. Mason DE
4. Keenan MC
5. Hill C
6. Goddard WA
7. Peters EC
8. Driggers EM
9. Hsieh-Wilson LC
(2012) Phosphofructokinase 1 glycosylation regulates cell growth and metabolism
Science 337:975–980.

https://doi.org/10.1126/science.1222278
- Google Scholar

Article and author information

Author details

Alexander A Shestov

Division of Nutritional Sciences, Cornell University, Ithaca, United States

Contribution
AAS, Conception and design, Acquisition of data, Analysis and interpretation of data

Contributed equally with
Xiaojing Liu

Competing interests
No competing interests declared.
Xiaojing Liu

Division of Nutritional Sciences, Cornell University, Ithaca, United States

Contribution
XL, Conception and design, Acquisition of data, Analysis and interpretation of data

Contributed equally with
Alexander A Shestov

Competing interests
No competing interests declared.
Zheng Ser

Division of Nutritional Sciences, Cornell University, Ithaca, United States

Contribution
ZS, Acquisition of data, Analysis and interpretation of data

Competing interests
No competing interests declared.
Ahmad A Cluntun

Field of Biochemistry and Molecular Cell Biology, Department of Molecular Biology and Genetics, Cornell University, Ithaca, United States

Contribution
AAC, Acquisition of data, Analysis and interpretation of data

Competing interests
No competing interests declared.
Yin P Hung

Department of Neurobiology, Harvard Medical School, Boston, United States

Contribution
YPH, Acquisition of data, Analysis and interpretation of data

Competing interests
No competing interests declared.
Lei Huang

Field of Computational Biology, Department of Biological Statistics and Computational Biology, Cornell University, Ithaca, United States

Contribution
LH, Acquisition of data, Analysis and interpretation of data

Competing interests
No competing interests declared.
Dongsung Kim

Field of Biochemistry and Molecular Cell Biology, Department of Molecular Biology and Genetics, Cornell University, Ithaca, United States

Contribution
DK, Acquisition of data, Analysis and interpretation of data

Competing interests
No competing interests declared.
Anne Le

Department of Pathology and Oncology, Johns Hopkins University School of Medicine, Baltimore, United States

Contribution
AL, Contributed to the analysis of the experiments with FX11, Contributed unpublished essential data or reagents

Competing interests
No competing interests declared.
Gary Yellen

Department of Neurobiology, Harvard Medical School, Boston, United States

Contribution
GY, Contributed to the NADH/NAD+ measurements, Contributed unpublished essential data or reagents

Competing interests
No competing interests declared.
John G Albeck

Department of Cell Biology, Harvard Medical School, Boston, United States

Contribution
JGA, Acquisition of data, Analysis and interpretation of data

Competing interests
No competing interests declared.
Jason W Locasale
1. Division of Nutritional Sciences, Cornell University, Ithaca, United States
2. Field of Biochemistry and Molecular Cell Biology, Department of Molecular Biology and Genetics, Cornell University, Ithaca, United States
3. Field of Computational Biology, Department of Biological Statistics and Computational Biology, Cornell University, Ithaca, United States
Contribution
JWL, Conception and design, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article

For correspondence
locasale@cornell.edu

Competing interests
JWL: A patent related to this work has been filed. US Provisional Patent Appln. No. 61/908,953.

"This ORCID iD identifies the author of this article:" 0000-0002-7766-3502

Funding

National Institutes of Health (R00CA168997)

Jason W Locasale

National Institutes of Health (R01AI110613)

Jason W Locasale

International Life Sciences Institute Research Foundation (Future Leader Award)

Jason W Locasale

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We are grateful to Chi V Dang, Navdeep Chandel, Dimitris Anastasiou, and members of the Locasale lab for helpful comments.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.