Computational Modeling of the Metabolic States Regulated by the Kinase Akt

Signal transduction and gene regulation determine a major reorganization of metabolic activities in order to support cell proliferation. Protein Kinase B (PKB), also known as Akt, participates in the PI3K/Akt/mTOR pathway, a master regulator of aerobic glycolysis and cellular biosynthesis, two activities shown by both normal and cancer proliferating cells. Not surprisingly considering its relevance for cellular metabolism, Akt/PKB is often found hyperactive in cancer cells. In the last decade, many efforts have been made to improve the understanding of the control of glucose metabolism and the identification of a therapeutic window between proliferating cancer cells and proliferating normal cells. In this context, we have modeled the link between the PI3K/Akt/mTOR pathway, glycolysis, lactic acid production, and nucleotide biosynthesis. We used a computational model to compare two metabolic states generated by two different levels of signaling through the PI3K/Akt/mTOR pathway: one of the two states represents the metabolism of a growing cancer cell characterized by aerobic glycolysis and cellular biosynthesis, while the other state represents the same metabolic network with a reduced glycolytic rate and a higher mitochondrial pyruvate metabolism. Biochemical reactions that link glycolysis and pentose phosphate pathway revealed their importance for controlling the dynamics of cancer glucose metabolism.

Introduction

Growth factors and nutrients are required for cell growth and proliferation in multicellular organisms. As a consequence of growth factors withdrawal normal cells undergo apoptosis, while most transformed cells escape the regulatory mechanisms and acquire the ability to proliferate even in the absence of growth signals (Edinger, 2005). In both normal and cancer cells the onset of proliferation induces important changes in cellular metabolism. Therefore metabolic activities in proliferating cells are fundamentally different from those in non-proliferating cells (DeBerardinis et al., 2008; Lunt and Vander Heiden, 2011).

The correlation between signal transduction pathways and cellular metabolism is mediated by some key components of the growth factor-induced cascades; typically these elements are protein kinases at the core of physiology and disease. Several growth factor-induced signal transduction pathways have been characterized so far and, in particular, the phosphoinositide 3-kinase (PI3K) is a key component downstream of the receptor tyrosine kinases (RTKs; Cantley, 2002). The PI3K is responsible for the production of 3-phosphoinositide lipid second messengers including phosphoinositol trisphosphate (PIP3) at the cell membrane. PIP3, in turn, contributes to the recruitment and activation of a wide range of downstream targets, among which the serine-threonine protein kinase Akt, also known as protein kinase B (PKB; Nicholson and Anderson, 2002; Gonzalez and McGraw, 2009). Akt/PKB is phosphorylated at two sites, one within the T-loop of the catalytic domain by the phosphoinositide-dependent kinase 1 (PDK1) and the other within the carboxyl terminal hydrophobic domain by the mammalian target of rapamycin complex 2 (mTORC2; Alessi et al., 1997; Sarbassov et al., 2005).

Fully activated Akt/PKB translocates from the cell membrane to the cytosol and nucleus where it phosphorylates its substrates (Manning and Cantley, 2007) to regulate multiple functions including cellular metabolism (Figure 1A). One of the chief mechanisms of Akt/PKB promoting cell growth and proliferation is through the activation of mTOR complex 1 (mTORC1), which is regulated by both nutrients and growth factor signaling (Wullschleger et al., 2006; Zoncu et al., 2011). Moreover, mTORC1 directly enhances the transcriptional activity of hypoxia-inducible factor 1α (HIF-1α; Land and Tee, 2007). HIF-1α is known to control the expression of several genes involved in energy metabolism, apoptosis, angiogenesis, and metastasis (Carmeliet et al., 1998; Pugh and Ratcliffe, 2003; Marín-Hernández et al., 2009).

FIGURE 1

Figure 1. The PI3K/Akt/mTOR pathway regulates central carbon metabolism. (A) PI3K/Akt/mTOR pathway. Signaling through the PI3K/Akt/mTOR pathway begins with the activation of RTKs in response to growth factors, leading to auto-phosphorylation on tyrosine residues and trans-phosphorylation of adaptor proteins. The PI3K is responsible for the production of 3-phosphoinositide lipid second messengers, including PIP3, which contributes to the activation of many downstream targets, such as PDK1 and mTORC2. Both PDK1 and mTORC2 activate, through phosphorylation in different sites, the serine-threonine protein kinase Akt. Akt regulates multiple functions including cellular metabolism, by promoting cell growth and proliferation through the activation of mTORC1, which also enhances the transcriptional activity of HIF-1α. Dashed lines represent the negative regulation of the PI3K/Akt/mTOR pathway by the action of mTORC1 feedback mechanism. (B) The metabolic network with the main reactions of glucose metabolism. A schematic representation of the glucose metabolism network, using the SBGN notation (LeNovere:2009]), is presented. The main pathways involved in the glucose metabolism are considered: glycolysis, PPP, the glycogen synthesis and degradation, lactate, and MPM branches. The metabolic targets regulated by PI3K/Akt/mTOR pathway are represented on the network: the PI3K/Akt/mTOR direct regulation is presented in yellow; the PI3K/Akt/mTOR indirect regulation (via HIF-1α) is presented in pink; the PI3K/Akt/mTOR direct and indirect regulation is presented in orange. All the PI3K/Akt/mTOR direct and indirect targets considered here are positively regulated, with the only exception of the MPM. Allosteric regulators (modifiers), activators (+), or inhibitors (−), are depicted in red. Metabolites – ADP, adenosine diphosphate; AMP, Adenosine Monophosphate; ATP, Adenosine Triphosphate; BPG, 1,3-bisphosphoglycerate; DHAP, dihydroxyacetone phosphate; E4P, erythrose-4-phosphate; F6P, fructose-6-phosphate; F16P, fructose-1,6-bisphosphate; G1P, glucose-1-phosphate; G6P, glucose-6-phosphate; GAP, glyceraldehyde-3-phosphate; GLC_c, cytoplasmic glucose; GLC_e, extracellular glucose; GLY, glycogen; LAC, lactate; NAD, nicotinamide adenine dinucleotides; NADH, nicotinamide adenine dinucleotides; NADP, nicotinamide adenine dinucleotide phosphate; NADPH, nicotinamide adenine dinucleotide phosphate; Pi, inorganic phosphate; PEP, phosphoenolpyruvate; PG2, 2-phosphoglycerate; PG3, 3-phosphoglycerate; PGN, 6-phosphogluconolactone; PRPP, phosphoribosylpyrophosphate; PYR, pyruvate; R5P, Ribose-5-phosphate; RU5P, ribulose-5-phosphate; X5P, xylulose-5-phosphate; S7P, sedoheptulose-7-phosphate; reactions – AK, adenylate Kinase; ATPase, ATP hydrolysis; DHase, NADH oxidation; DPHase, NADPH oxidation; ENO, Enolase; G6PDH, glucose-6-P dehydrogenase; GAPDH, glyceraldehyde dehydrogenase; GLUT, glucose transporter; GPa, glycogen phosphorylase A; GPb, glycogen phosphorylase b; GS, Glycogen synthase; FBA, fructose-6-P aldolase; HK, Hexokinase; LDH, lactate dehydrogenase; MPM, mitocondrial pyruvate metabolism; PFK, phosphofructo-kinase; PGDH, phoshogluconolactone dehydrogenase; PGI, phosphoglucoisomerase; PGK, phosphoglycerate kinase; PGLM, phosphoglucomutase; PGYM, 3-phosphoglycerate mutase; PK, pyruvate kinase; PRPPS, phosphoribosylpyrophosphate synthetase; R5PI, ribose-5-P isomerase; RUPE, ribulose-phosphate-3 epimerase; TAL, transaldolase; TKL, transketolase, reaction I; TKL2, transketolase, reaction II; TPI, triose-phosphate isomerase.

Negative regulation of the PI3K/Akt/PKB pathway is primarily accomplished through the action of the PTEN tumor suppressor protein, a lipid and protein phosphatase whose main lipid substrate is PIP3 (Song et al., 2012). Recently, a crucial mTORC1-dependent feedback mechanism has been elucidated (Howell and Manning, 2011). According to the current knowledge mTORC1 can exert a negative regulation of its upstream signaling molecules as depicted in Figure 1A.

Overall, the role of the PI3K/Akt/PKB signaling pathway in oncogenesis has been extensively investigated and altered expression or mutations of many components of this pathway have been implicated in human cancer (Vivanco and Sawyers, 2002; Carnero, 2010). Indeed, expression of constitutively active forms of Akt/PKB can prevent cell death upon growth factor withdrawal. PI3K/Akt/mTOR-mediated survival relies on a profound metabolic adaptation, including aerobic glycolysis (known as Warburg effect).

In the last decade, many efforts have been made to improve the understanding of the control of glucose metabolism. High-rate glycolysis is supported by a number of molecular alterations, which are either unique to cancer in an otherwise healthy organism or druggable with controllable toxicities (Porporato et al., 2011). In general, the goal is to find a therapeutic window between proliferating cancer cells and proliferating normal cells (Vander Heiden, 2011) in the development of successful cancer therapies targeting cellular metabolism. To be an attractive candidate for cancer therapy, there must be a significant difference in enzyme activity between cancer cells and normal proliferating cells (Hamanaka and Chandel, 2012).

In this context, we have modeled the link between the PI3K/Akt/mTOR pathway and metabolism, specifically the glycolytic rate, lactic acid production, and nucleotide biosynthesis. Computational models are a tool to understand, control, and predict the complex behavior of biological systems; for example, computer-based simulations can be a powerful alternative to wet experiments for exploring the impact of perturbations, such as drug treatments, on a molecular network (Kreeger and Lauffenburger, 2010). Indeed, dynamical modeling of complex biological systems has been widely used in the last decade. In particular, Ordinary Differential Equations (ODEs) have been exploited as a standard numerical method in many successful examples (Novak and Tyson, 1993; Tyson et al., 2001; Bazil et al., 2010). Before becoming a predictive tool, a model must be defined and corroborated considering the available knowledge and experimental data. However, current experimental techniques limit our ability to obtain in vivo information about some aspects of metabolic pathways. In fact, while metabolic fluxes can be measured in vivo using metabolically active molecules labeled with stable isotopes, it is hard or even impossible to obtain accurate in vivo measures of enzyme properties. Therefore no existing experimental dataset contains all the information for the biochemical processes and molecular species under analysis. Nevertheless, several experimental datasets, collected using a variety of tissue sources, techniques and experimental conditions (e.g., in vitro), are available for metabolic pathways. Therefore, these data can be carefully integrated to fill the lack of information, so that models can be optimized by means of various simulation and analysis tools (Ghosh et al., 2011). With the addition of new experimental datasets, the model is refined through possible modifications of its structure or parameter values.

Several models have been developed to describe the glucose metabolism (Mosca et al., 2012). To study the accelerated glycolysis of human cancer cells, Marín-Hernández et al. (2011) constructed a model supported by experimental data for enzyme kinetics, steady state pathway metabolite concentrations, and metabolic fluxes. This model was used as the base model for our work due to the type of cells used (a human cancer cell line) and the amount of experimental data collected in the same condition. To accomplish the task of modeling the metabolic effects of the PI3K/Akt/mTOR pathway, we extended the network configuration and replaced simplified rate equations for glycolytic branches with thermodynamically consistent and detailed representations. Subsequently, some parameter values were optimized to obtain a quantitative reproduction of the metabolic state (fluxes and metabolites) of a proliferating cancer cell line. In the resulting extended model, we modulated the activity of PI3K/Akt/mTOR in order to simulate the effects of the kinase levels over the rates of the processes controlled by its (direct and indirect) targets. Then, we identified the reactions exerting a major control over the metabolic network at two different rates of glucose metabolism, corresponding to high and low PI3K/Akt/mTOR activity. Enzymes that have a different relevance in the two metabolic states represent potential targets for a more selective control of the system.

Results

Model Development

The model represents the dynamics of the metabolic network composed by the main metabolites and biochemical processes (transports and reactions) of the glucose metabolism: the glucose transport, the glycolytic reactions, the lactic acid production, the glycogen synthesis/degradation, and the pentose phosphate pathway (PPP; Figure 1B). The rates of the biochemical processes considered in the model are defined by taking into account the details of the respective kinetic mechanisms. Therefore, the main activators and inhibitors forming feedback and feedforward regulations are also considered. To maintain thermodynamic consistency and since equation simplification by irreversibility limits the possible steady states allowed in the system, biochemical processes are represented as thermodynamically balanced and reversible. This approach was also important to enable the reproduction of the wide range of glycolytic rates observed in relation to the regulatory activities operated by oncogenes (Levine and Puzio-Kuter, 2010). The only exception is the ATPase reaction, which is used here as elsewhere (e.g., Lambeth and Kushmerick, 2002; Bazil et al., 2010; Marín-Hernández et al., 2011) to simulate, with a unique flux, all the ATP consuming reactions that are not detailed in the model.

The model is a 28 state system of differential algebraic equations (DAEs) that consists in 25 ODEs and 3 algebraic equations to calculate ADP, NADH, and NADPH concentrations (see Appendix). As previously mentioned, this model is primarily based on a recent cancer glycolysis model of human cervical cancer HeLa cells (Marín-Hernández et al., 2011). The model of HeLa cells glycolysis was developed to reproduce a specific metabolic steady state, making use of detailed rate equations for glycolysis and simplified irreversible reactions for the PPP, the glycogen synthesis and degradation, and the mitochondrial pyruvate metabolism (MPM) branches. Since these glycolytic branches also play an important role in relation to the metabolic states determined by PI3K/Akt/mTOR signaling, we extended the network of reactions substituting simplified branches with thermodynamically balanced and reversible kinetics. In particular, rate equations for PPP, glycogen synthesis, and degradation were defined as in the model of Holzhütter (2004); the kinetics of phosphoglucomutase (PGLM) and glycogen synthase (GS) were derived from Li et al. (2010b); the glycogen phosphorylase (GP) kinetics from Lambeth and Kushmerick (2002); the MPM was defined using a common modular rate law (Liebermeister et al., 2010). The model is available in BioModels database (Li et al., 2010a) with identifier MODEL1210150000.

Model Fitting

The model was optimized assembling a dataset that represents the typical traits of a cancer cell growing in condition of non-limiting glucose: increased glycolytic rate, active synthesis of NADPH and nucleotides, predominant metabolization of glucose to lactate rather than to CO₂ and H₂O through the tricarboxylic acid cycle (Levine and Puzio-Kuter, 2010).

As in several other recent computational models of metabolic pathways (e.g., Bazil et al., 2010; Li et al., 2010b), we adjusted only a small subset of the parameters appearing in the model equations and, more precisely, we focused on the maximum rates of forward reactions (Table 1). Biologically, this can be interpreted as the correction of enzyme concentrations or enzyme activities in order to explain, with the equations used in the model, the metabolic values (fluxes and concentrations) measured experimentally. Crucially, by constructing the model taking into account the Haldane constraint, which relates the equilibrium constant of a biochemical reaction with (forward and reverse) maximum rates and reactant constants (Appendix), the adjustment/optimization of the parameter values will not violate the thermodynamics of any reaction or of the entire network (Lambeth and Kushmerick, 2002), but only physiologically appropriate kinetic parameter values will result in reasonable fluxes and accurate dynamics.

TABLE 1

Table 1. Maximum reaction rates in conditions H and L.

Metabolite concentrations and metabolic fluxes were primarily collected from two experimental studies on human cervix cancer HeLa cells (Reitzer et al., 1980; Marín-Hernández et al., 2011). Reitzer et al. (1980) measured the flux of sugar carbon per unit of protein synthesized in the major pathways of metabolism of HeLa cells growing on 10 mM glucose. Under this condition, they observed that the 10% of the glucose converted to glucose-6-P, which was quantified as 11.3 nmol/min/mg, is diverted toward the oxidative arm of the pentose cycle, 5% to purine metabolism and nucleotide biosynthesis, 1% to glycogen, and 80% to lactate (Figure 2A). The optimized model provided a satisfactory representation of the experimental data considered, in terms of fluxes (Figure 2A) and metabolite concentrations (Figure 2B).

FIGURE 2

Figure 2. Metabolic fluxes and metabolite concentrations. (A) A Simplified illustration of the metabolic network, where the width of the arrows is proportional to the predicted flux values, which are reported close to the respective arrows; experimental values are shown in parentheses; fluxes are in nmol/min/mg unit. (B) Predicted (white) and experimental (black) metabolite concentrations, reported as the log₁₀ of nM values. The full list of predicted and experimental fluxes and concentration values in conditions H and L is available in the Appendix.

In silico Modulation of PI3K/Akt/mTOR Activity

We considered the most established metabolic targets regulated by PI3K/Akt/mTOR and by HIF-1α, as indicated in Figure 1B. The PI3K/Akt/mTOR pathway activation increases the synthesis of GLUT, the main glucose transporter in most cell types (Kohn et al., 1996), and enhances its transcription and its translocation from the cytosol to the plasma membrane, increasing glucose uptake (Barthel et al., 1999; Rathmell et al., 2003). Upon PI3K/Akt/mTOR stimulation, the activity of the hexokinase (HK) is enhanced (Elstrom et al., 2004). PI3K/Akt/mTOR has also effects on other steps of glycolysis; in fact, it has been shown that increasing GLUT and HK expression does not enhance the glycolytic flux to the observed levels with constitutive activation of PI3K/Akt/mTOR (Rathmell et al., 2003). Glycolysis downstream targets of PI3K/Akt/mTOR include PFK2; phosphorylation and activation of PFK2 lead to allosteric activation of PFK1 (Deprez et al., 1997). Moreover, PI3K/Akt/mTOR inhibits GSK, the GS kinase-3beta (Yoeli-Lerner et al., 2009), which inhibits the GS; as a consequence, PI3K/Akt/mTOR has indirect positive effects on the glycogen synthesis. Also the ATP-citrate lyase, the primary enzyme responsible for the synthesis of cytosolic acetyl-CoA, is a substrate for PI3K/Akt/mTOR (Berwick et al., 2002). An important target of PI3K/Akt/mTOR is HIF-1α, which can be overexpressed also under non-hypoxic conditions through PI3K/Akt/mTOR (Lee et al., 2008); in turn, HIF-1α is responsible for the positive regulation of several enzymes of the central metabolism: almost all the glycolytic enzymes (Semenza et al., 1994), G6PDH, PGDH (Guo et al., 2009), TKL (Zhao et al., 2010), and LDH (Firth et al., 1995); HIF-1α also stimulates the inhibition of PDH (Biswas et al., 2010).

To model the metabolic effects of the PI3K/Akt/mTOR signaling pathway, we modified the rate equations for the biochemical processes regulated by the targets of PI3K/Akt/mTOR. Rate equations were modified on the basis of the specific effect PI3K/Akt/mTOR exerts over each of its targets. If PI3K/Akt/mTOR positively (negatively) regulates the protein concentration, we increased (decreased) the value of the maximum rate of the forward reaction (V_f) of the corresponding process, because this quantity is directly proportional to the total protein concentration (Sauro, 2011). This is the case of GLUT, GS, glycolytic enzymes, and PPP enzymes (Figure 2). The enhancement of PFK activity mediated by PI3K/Akt/mTOR is also exerted through F26P, one of PFK allosteric activators, and therefore we increased the F26P concentration. The positive relation between PI3K/Akt/mTOR and fatty acid synthesis was taken into account increasing the rate of NADPH consuming reactions (fatty acid synthesis requires NADPH), while the inhibition of PDH due to PI3K/Akt/mTOR was simulated reducing the rate of the MPM process.

Starting from the metabolic state representing the typical traits of a cancer cell growing in condition of non-limiting glucose (Figure 2) and modulating PI3K/Akt/mTOR activity, we obtained several metabolic steady states (Figure 3). Coherently with what is observed experimentally (DeBerardinis et al., 2008), the flux through glycolysis, glycogen synthesis, and the portion of pyruvate metabolized toward lactate increased as the PI3K/Akt/mTOR activity was increased and vice versa. We considered two metabolic steady states as representatives of two different glycolytic rates, corresponding to high or low PI3K/Akt/mTOR activity: the first, representing a cancer cell line where PI3K/Akt/mTOR promotes a high-rate of glucose metabolism (Figures 2 and 3, condition H); the second characterized by a lower glycolytic rate due to a reduced PI3K/Akt/mTOR signal (Figure 3, condition L).

FIGURE 3

Figure 3. Lower metabolic fluxes are obtained reducing PI3K/Akt/mTORs activity. Time variation of some representative fluxes (nmol/min/mg) before (condition H) and after the reduction of PI3K/Akt/mTOR signaling (at t ( 7dd). The thick lines indicate the dynamics that lead to the representative condition L (right).

In relation to condition H, condition L was selected for its reasonable reproduction of the available experimental data and for conserving the direction of the metabolic fluxes. In fact, on the one hand the variation of glycolysis and lactic acid production rate between state H and L was respectively of ∼1.7- and ∼2.0-fold (Figure 3; Appendix), close to the respective values of ∼1.5 and ∼3.0 recently measured by Fan et al. (2010) after activation of Akt. On the other hand, by conserving the same flux direction between conditions H and L we were able to exclude effects due to a different configuration of the flow of mass within the system, enabling a meaningful comparison focused on metabolic rates between two states that reproduce the metabolic variations measured in relation to PI3K/Akt/mTOR activity (for instance, other possible configurations could be: the non-oxidative arm of the PPP that extracts part of the glycolytic flux and redirects it toward nucleotide biosynthesis; glycogen degradation higher than glycogen synthesis).

Local Sensitivity Analysis

A first-order sensitivity analysis (SA) on the model was performed to determine how robust the steady states under analysis were in relation to local perturbations in the parameter values. These sensitivity indexes represent a measure of the degree to which the steady states H and L are sensitive to the value of individual parameters. A high sensitivity value indicates that changing a given parameter determines significant changes to the considered dynamics. The averages (median) of the sensitivity indexes, of all the 29 adjustable parameters, are respectively 3.198e-03 and 7.320e-04 for conditions H and L (Table 1). Therefore, on average a small perturbation of a given parameter determines a change in the steady state fluxes that is approximately three orders of magnitude less than the perturbation.

Comparison of the Enzyme Control Over the System at Different Rates of Glucose Metabolism

To study the global sensitivity of conditions H and L in relation to each biochemical process of the metabolic network (listed in Table 1), we calculated, for each possible pairs of biochemical processes, the Derivative Based Global Sensitivity Measures (DGSMs; Kucherenko et al., 2009) that summarize the (total) effect exerted by a reaction on the steady state flux of another reaction around the metabolic steady states H and L. Conversely form the Metabolic Control Analysis (MCA; Fell, 1997), DGSMs are a method for SA that copes with model non-linearities and interactions between the quantities under analysis. In fact, this method does not require a linearization of the system and considers the effect of perturbing more than one quantity at the same time, but it is computationally more intensive than MCA (see Materials and Methods).

In state H, steady fluxes are more sensitive to enzymes of the PPP (Figure 4A): G6PDH, first enzyme of the PPP and also target of p53 (Jiang et al., 2011); PRPPS, which catalyzes the phosphoribosylation of ribose-5-phosphate to 5-phosphoribosyl-1-pyrophosphate, a metabolite that is necessary for purine metabolism and nucleotide biosynthesis; TKL, a thiamine-dependent enzyme that channels sugar phosphates between glycolysis and the PPP. GLUT, the source of mass of the system, also exerts a major control. With the exceptions of HK, FBA, TPI, and ENO, the variation of other glycolytic enzyme activity appears to exert a similar effect on the steady state metabolic fluxes in the two conditions. In condition L, steady state fluxes are still sensitive to GLUT, G6PDH, and PRPPS, even if the amount of control differs; in contrast, TKL loses its role and PFK, HK, and ATPase emerge as important players for system dynamics.

FIGURE 4

Figure 4. Heatmap of the steady state flux sensitivity to reactions. (A) Sensitivities of the steady state fluxes (rows) in relation to the reactions (columns) in condition H; red: high sensitivity; green: low sensitivity. (B) Differential ranking of steady state flux sensitivities (rows) to reactions (columns) when comparing conditions H and L; red: high sensitivity in condition H; green: high sensitivity in condition L.

To identify the reactions exerting a different control over the system when using the two conditions H and L, we calculated, for each reaction, the normalized difference in the ranking of all the reactions based on their relevance in state H and L. When comparing the ranking of the sensitivity indexes between the two metabolic steady states, we observed that the PPP enzymes G6PDH and TKL showed a prominent role in condition H; ENO, responsible of the conversion between PG2 and PEP, also had a relatively higher control in H in comparison to L (Figure 4B). Conversely, in condition L the system is more sensitive to GLUT, PFK, and to some of the enzymes of the glycolytic “pay off” phase. This trend is also evident when considering the top ranked enzymes for the control of specific steps (Figure 5).

FIGURE 5

Figure 5. Reactions exerting the most of the control over the metabolic network. Ranking of the top reactions having the highest influence on the flux of GLUT (A), G6PDH (B), PRPPS (C), and LDH (D) in H (black) and L (white). A value of 1 corresponds to the best rank.

Another interesting question that arises with the type of SA we have performed concerns the reactions having a high control over a specific flux or set of fluxes. To this aim we calculated, for each flux, the relative log sensitivity (see Materials and Methods). This analysis underlined that glycolysis is specifically sensitive mainly to GLUT, HK, and PGI (Figure 6).

FIGURE 6

Figure 6. Heatmap of the relative log sensitivities of steady state fluxes (rows) in relation to reactions (columns) in condition H. Red: high sensitivity; green: low sensitivity.

Importantly, the structure of the metabolic network emerges from the clustering of the quantities illustrated in Figures 4 and 6: SA indexes referring to the variation of steady state fluxes of neighboring reactions are similar. This expected result is the consequence of having performed the SA at steady state, where mass conservation constrains the steady state fluxes of neighboring reactions according to their stoichiometry.

Discussion

Several computational models have been proposed to study the control of the glycolytic flux. In Table 2 we summarized the main results obtained by those works that we based the definition of our model on. Marín-Hernández et al. (2006) performed an elasticity-based control analysis of the glycolytic flux in rat hepatoma cells (AS-30D); the authors observed that the major control of the glucose metabolism is due to the upstream part of glycolysis (GLUT, HK, GP, PGI, PFK). In a subsequent analysis based on a kinetic model, the same group calculated that glycolysis is mainly controlled by HK, HPI, and GLUT in AS-30D cells while glycogen degradation, GLUT and HK are key steps for HeLa cells. Using MCA, Schuster and Holzhutter (1995) calculated that the glycolytic flux in erythrocytes is mainly controlled by ATP utilization. Several computational studies have been done on the glycolysis of skeletal muscle cells. For example, Lambeth and Kushmerick (2002) identified in GP and PK the main steps of glycolytic control at resting flux; PFK and PGLM become important when the ATP utilization increases while GP activity gains importance in the network as the glycolytic flux increases.

TABLE 2

Table 2. Major controller of glycolysis reported in other computational and experimental works.

Comparison of the results obtained by different works should be done taking into account the structure of the metabolic network, the rate equations and the model configuration used. For example, a relevant control of the glycogen degradation step over the glycolytic flux was reported when glycogen was one of the sources of glucose for the system, as in Marín-Hernández et al. (2006) and Lambeth and Kushmerick, 2002 (Table 2). This was not the scenario of condition H, but, despite the differences in the metabolic network and metabolic state that we considered in comparison to Marín-Hernández et al. (2011), we also found that GLUT and HK exert a relevant control for glycolysis (Figure 6).

We compared two metabolic states generated by the specific variation of the fluxes regulated by the activity of the PI3K/Akt/mTOR pathway. One state represented the metabolism of a growing cancer cell characterized by aerobic glycolysis and cellular biosynthesis (condition H), while the other (condition L) represented the same metabolic network with a reduced glycolytic rate, a reduced lactic acid production, but a higher MPM, as reported in literature in relation to a lower activity of PI3K/Akt/mTOR (DeBerardinis et al., 2008). Some steps of the metabolic network that link glycolysis and PPP, namely those catalyzed by the G6PDH and TKL enzymes, revealed their importance for the cancer metabolic state. Results from our model may provide insight and assist in the selection of drug targets in anticancer treatments.

Results gained with current models in biology are still restricted by the assumptions made in relation to the limits of the current technologies, the lack of detailed information collected in the same experimental conditions and the need to simplify biological complexity (aspects as enzymes isoforms, molecular crowding and the temporal variation of enzyme concentration are usually not accounted for). Nevertheless, considering that model definition is related to the questions the modeler wants to answer, current models can be a useful tool for directing experimental efforts on a subset of all possible hypotheses. In relation to the crucial steps that emerged in this work for the dynamics of the metabolic network, it is worth to mention that, recently, several authors have recognized G6PDH and TKL as potential anticancer targets (e.g., Tennant et al., 2010; Vander Heiden, 2011). In general, there are several leading therapeutic compounds targeting glucose metabolism in preclinical and clinical phases for many tumor types, such as solid tumors (lung, breast, prostate, gastric), metastatic melanoma and renal cell carcinoma (Tennant et al., 2010; Porporato et al., 2011; Vander Heiden, 2011). Other principal targets that are currently enrolled in clinical and/or preclinical studies are: GLUT, HK, PK, LDH (preclinical and clinical; Mohanti et al., 1996; Singh et al., 2005; Christofk et al., 2008; Mathupala et al., 2009; Le et al., 2010; Wolf et al., 2011); PFK, PGYM, TKL, G6PDH (preclinical only; Kuo et al., 2000; Evans et al., 2005; Clem et al., 2008; Furuta et al., 2010).

It would be interesting to extend our analysis in order to compare other representative steady states, such as, for example, those considered by Sengupta et al. (2011), who evaluated different glucose metabolization patterns between glycolysis and PPP in growing and non-growing conditions. Moreover, our current model could be linked to a model of mitochondrion (e.g., Bazil et al., 2010) in order to include detailed representations of the tricarboxylic acid cycle and other relevant metabolic reactions that take place in the mitochondrion.

Materials and Methods

Numerical Solutions

The DAE system representing the metabolic network was numerically integrated using MATLAB (2008b) and the stiff ode solver ode15s with absolute and relative tolerances of 10⁻⁹ and 10⁻⁶ respectively. Steady states were identified using the MATLAB function fsolve with default options. Model optimization and sensitivity analyses were done on HP(R) workstations equipped with two 2.50 GHz INTEL(R) Quad-core Xeon(R) E5420 processors and 10 GB RAM. The results obtained were displayed using MATLAB.

Model Optimization

Recently, it has been observed that multi-objective optimization have significant benefits compared to single objective approaches (Handl et al., 2007). Model fitting was formulated as a multi-objective optimization problem aiming at the simultaneous minimization of the difference between model predictions and experimentally determined concentrations, enzyme activities, and steady state fluxes. In detail, two objectives [f₁(x), f₂(x)] were defined as

$f_{1,2}\left(\mathbf{\text{x}}\right)=1∕N{\Sigma }_{i=1,\dots ,N}{log}_{10}\left(x_i∕x_i^{*}\right)$

subject to

$\begin{array}{llll}& {\mathbf{\text{J}}}^{*}=\mathbf{\text{J}}\left(\mathbf{\text{x}}\right)& & \\ & \mathbf{\text{x}}>0& & \end{array}$

where $x_i^{*}$ is the experimental value for the concentration of a metabolite (in the case of f₁) or enzyme V_mf (for f₂), x_i is the corresponding value used in the model, N is the number of elements (metabolites or enzymes), J* is the vector of experimental values of enzyme fluxes and J(x) are the respective model predictions obtained using x. The multi-objective optimization problem was solved using the Non-Dominated Sorting Genetic Algorithm II (Deb et al., 2002), which is one of the most popular methods in the field of multi-objective optimization. The NSGA-II algorithm was run in parallel using several populations of 1000 individuals for 100 generations. Every five generations, the worst five individuals of a population were replaced by a random selection of the best five individuals from the other populations.

Modeling the Metabolic Effects of PI3K/Akt/mTOR

The metabolic effects of PI3K/Akt/mTOR were modeled according to the mechanism of interaction with its targets. Parameter values used to create condition L were obtained from condition H, multiplying a specific quantity appearing in the rate equation for the biochemical process regulated by a target of PI3K/Akt/mTOR by quantity α in order to reduce or increase the target activity. In detail, for GLUT, HK, PGI, GS, G6PDH, PGDH, TAL, TKL, TKL2, FBA, TPI, GAPDH, PGK, ENO, PK, LDH, and DPHase, we multiplied the respective V_f by α = 0.56, while for MPM we multiplied the respective V_f by α = 1.16; for PFK, we also multiplied the concentration of its allosteric activator F26P by α = 0.56. The V_f values used to obtain steady states H and L are listed in Table 1. Rate equations are listed in Appendix.

Sensitivity Analysis

Sensitivity analysis can be defined as the study of how uncertainty in the output of a model can be apportioned to different sources of uncertainty in the model input (Saltelli et al., 2000). In most of the current systems biology literature, sensitivities indexes are estimated calculating derivatives of a model output in a specific state of the system (local approach) corresponding to a particular model parameterization; moreover, only the variation of one parameter at a time is considered. For example, control coefficients estimated in the context of MCA are scaled partial derivatives calculated on the model linearized around a steady state; thus, MCA quantifies how a model output is influenced by infinitesimal changes in a parameter. As a consequence, results of MCA are restricted to infinitesimal parameter changes and do not account for interactions between parameters.

To overcome this issue, especially considering that the biological quantities represented with parameters (e.g., the enzyme concentration) are subject to significant variations, in our study we carried out the SA calculating the DGSMs (Kucherenko et al., 2009). The method is based on averaging local derivatives using (quasi) Monte Carlo (MC) sampling.

To compute the DGSMs we generated, for each of the conditions H and L, a set of parameterizations p _k = (p_k1, p_k2, … , p_kn), n = 29, k = (1, 2, … , 500), p _k ∈ [p − εp, p + εp], where each p _k represents a possible perturbations of all the 29 V_f reference values (p) reported in Table 1. This sampling was carried out using the Sobol Sequences, a class of quasi-random low discrepancy sequences (Sobol, 1967). The parameterizations p _k were used to calculate the elementary effects

$\begin{array}{llll}{e}^{ijk}& =\left|\right.\left[J_i\left(p_{k1},\dots ,p_{kj}+\delta ,\dots ,p_{kn}\right)-J_i\left(p_{k1},\dots ,p_{kj},\dots ,p_{kn}\right)\right]\left|\right./& & \\ & \delta \cdot \left(p_{kj}∕\left|\right.J_i\left(p_{k1},\dots ,p_{kj},\dots ,p_{kn}\right)\left|\right.\right)& & \end{array}$

a measure of the impact of having introduced the perturbation δ in the j-th element of p _k on the steady state flux J_i, i = (1, 2, … , 29). Hence, for each of the two conditions H and L, we obtained a total of 500 29-by-29 matrices of elementary effects that were used to calculate the 29-by-29 matrices of scaled sensitivity indexes

$g_{ij}=\left(m_{ij}^2+s_{ij}^2\right)/{\displaystyle {\sum }_{j=1,\dots ,29}\left(m_{ij}^2+s_{ij}^2\right)}$

where m_ij and s_ij are, respectively, the mean and sample standard deviation of the elementary effects e_ijk.

The normalized differential ranking of the reactions was calculated as

$R_{ij}=\left(r_{ij}^L-r_{ij}^H\right)/\left(r_{ij}^L+r_{ij}^H\right)$

where r_ij = (1, 2, … , 29) is the rank of the scaled sensitivity index g_ij when ordering the elements (g_i1, g_i2, … , g_i29) from the highest (strong influence of parameter j on flux i) to the lowest, and the superscript L or H indicates the condition. Hence, R_ij will be positive (negative) for V_f parameters (and the corresponding reaction) exerting a higher (lower) control in condition H in relation to condition L; R_ij will be higher for reactions having a higher control in both conditions.

The relative log sensitivity was calculated as

$RL{S}_{ij}=lo{g}_{10}\left(g_{ij}/m_i^{\prime }\right)$

and indicates the degree of variation of a scaled sensitivity index g_ij in relation to the median m_i’ of all the scaled sensitivity indexes related to the same V_f parameter (g_1j, g_2j, … , g_29j).

For each of the two conditions H and L, the 29-by-29 matrix of the first-order local sensitivity indexes s_ij for steady state flux i to perturbation of V_f of reaction j was obtained by calculating the elementary effects e_ijk using the respective values V_f listed in Table 1, that is p _k = p.

We used ε = 5 · 10⁻⁶ and δ = 10⁻⁶; steady state fluxes J_i were obtained by numerical integration of the DAE system within the interval [0, 10⁴], followed by the identification of the steady state (see above the numerical solutions paragraph). For the calculation of the DGSMs, 2N(n + 1) numerical solutions of the model were required, where N = 500 represents the number of parameterizations, n = 29 is the number of the parameters, and the factor 2 is due to the two conditions H and L.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The work was carried out under the HPC-EUROPA2 (228398) projects with the support of the European Commission – Capacities Area – Research Infrastructures, and under the Italian FIRB-MIUR projects ITALBIONET (RBPR05ZK2Z), BIOPOPGEN (RBIN064YAT), HIRMA (RBAP11YS7K), and INTEROMICS Flagship Initiative. We thank Dr. I. Merelli of Institute for Biomedical Technologies (ITB-CNR) for his technical support with high performance computing and thank John Hatton (ITB-CNR) for proofreading the manuscript.

References

Alessi, D. R., James, S. R., Downes, C. P., Holmes, A. B., Gaffney, P. R., Reese, C. B., et al. (1997). Characterization of a 3-phosphoinositide-dependent protein kinase which phosphorylates and activates protein kinase Balpha. Curr. Biol. 7, 261–269.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Barthel, A., Okino, S. T., Liao, J., Nakatani, K., Li, J., Whitlock, J. P. Jr., et al. (1999). Regulation of GLUT1 gene transcription by the serine/threonine kinase Akt/PKB1. J. Biol. Chem. 274, 20281–20286.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bazil, J. N., Buzzard, G. T., and Rundell, A. E. (2010). Modeling mitochondrial bioenergetics with integrated volume dynamics. PLoS Comput. Biol. 6, e1000632. doi:10.1371/journal.pcbi.1000632

CrossRef Full Text

Berwick, D. C., Hers, I., Heesom, K. J., Moule, S. K., and Tavare, J. M. (2002). The identification of ATP-citrate lyase as a protein kinase B (Akt/PKB) substrate in primary adipocytes. J. Biol. Chem. 277, 33895–33900.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Biswas, S., Troy, H., Leek, R., Chung, Y. L., Li, J. L., Raval, R. R., et al. (2010). Effects of HIF-1alpha and HIF2alpha on growth and metabolism of clear-cell renal cell carcinoma 786-0 xenografts. J. Oncol. 2010, 757908.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cantley, L. C. (2002). The phosphoinositide 3-kinase pathway. Science 296, 1655–1657.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Carmeliet, P., Dor, Y., Herbert, J. M., Fukumura, D., Brusselmans, K., Dewerchin, M., et al. (1998). Role of HIF-1alpha in hypoxia-mediated apoptosis, cell proliferation and tumour angiogenesis. Nature 394, 485–490.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Carnero, A. (2010). The PKB/Akt/PKB pathway in cancer. Curr. Pharm. Des. 16, 34–44.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Casazza, J. P., and Veech, R. L. (1986). The interdependence of glycolytic and pentose cycle intermediates in ad libitum fed rats. J. Biol. Chem. 26, 690–698.

Christofk, H. R., Vander Heiden, M. G., Harris, M. H., Ramanathan, A., Gerszten, R. E., Wei, R., et al. (2008). The M2 splice isoform of pyruvate kinase is important for cancer metabolism and tumour growth. Nature 452, 230–233.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Clem, B., Telang, S., Clem, A., Yalcin, A., Meier, J., Simmons, A., et al. (2008). Small-molecule inhibition of phosphofructo-2-kinase activity suppresses glycolytic flux and tumor growth. Mol. Cancer Ther. 7, 110–120.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Deb, K., Pratap, A., Agarwal, S., and Meyarivan, T. (2002). A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Trans. Evol. Comp. 6, 182–197.

CrossRef Full Text

DeBerardinis, R. J., Lum, J. J., Hatzivassiliou, G., and Thompson, C. B. (2008). The biology of cancer: metabolic reprogramming fuels cell growth and proliferation. Cell Metab. 7, 11–20.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Deprez, J., Vertommen, D., Alessi, D. R., Hue, L., and Rider, M. H. (1997). Phosphorylation and activation of heart 6-phosphofructo-2-kinase by protein kinase B and other protein kinases of the insulin signalling cascades. J. Biol. Chem. 272, 17269–17275.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Edinger, A. L. (2005). Growth factors regulate cell survival by controlling nutrient transporter expression. Biochem. Soc. Trans. 33, 225–227.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Elstrom, R. L., Bauer, D. E., Buzzai, M., Karnauskas, R., Harris, M. H., Plas, D. R., et al. (2004). Akt stimulates aerobic glycolysis in cancer cells. Cancer Res. 64, 3892–3899.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Evans, M. J., Saghatelian, A., Sorensen, E. J., and Cravatt, B. F. (2005). Target discovery in small-molecule cell-based screens by in situ proteome reactivity profiling. Nat. Biotechnol. 23, 1303–1307.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Fan, Y., Dickman, K. G., and Zong, W. X. (2010). Akt and c-Myc differentially activate cellular metabolic programs and prime cells to bioenergetic inhibition. J Biol Chem. 285, 7324–7333.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Fell, D. (1997). Understanding the Control of Metabolism. London: Portland Press.

Firth, J. D., Ebert, B. L., and Ratcliffe, P. J. (1995). Hypoxic regulation of lactate dehydrogenase A. Interaction between hypoxia-inducible factor 1 and cAMP response elements. J. Biol. Chem. 270, 21021–21027.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Furuta, E., Okuda, H., Kobayashi, A., and Watabe, K. (2010). Metabolic genes in cancer: their roles in tumor progression and clinical implications. Biochim. Biophys. Acta 1805, 141–152.

Pubmed Abstract | Pubmed Full Text

Gao, H., and Leary, J. A. (2004). Kinetic measurements of phosphoglucomutase by direct analysis of glucose-1-phosphate and glucose-6-phosphate using ion/molecule reactions and Fourier transform ion cyclotron resonance mass spectrometry. Anal. Biochem. 329, 269–275.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ghosh, S., Matsuoka, Y., Asai, Y., Hsin, K. Y., and Kitano, H. (2011). Software for systems biology: from tools to integrated platforms. Nat. Rev. Genet. 12, 821–832.

Pubmed Abstract | Pubmed Full Text

Gonzalez, E., and McGraw, T. E. (2009). The Akt/PKB kinases: isoform specificity in metabolism and cancer. Cell Cycle 8, 2502–2508.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Guo, S., Miyake, M., Liu, K. J., and Shi, H. (2009). Specific inhibition of hypoxia inducible factor 1 exaggerates cell injury induced by in vitro ischemia through deteriorating cellular redox environment. J. Neurochem. 108, 1309–1321.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hamanaka, R. B., and Chandel, N. S. (2012). Targeting glucose metabolism for cancer therapy. J. Exp. Med. 209, 211–215.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Handl, J., Kell, D. B., and Knowles, J. (2007). Multiobjective optimization in bioinformatics and computational biology. IEEE/ACM Trans. Comput. Biol. Bioinform. 4, 279–292.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Holzhütter, H. G. (2004). The principle of flux minimization and its application to estimate stationary fluxes in metabolic networks. Eur. J. Biochem. 271, 2905–2922.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Howell, J. J., and Manning, B. D. (2011). mTOR couples cellular nutrient sensing to organismal metabolic homeostasis. Trends Endocrinol. Metab. 22, 94–102.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Jiang, P., Du, W., Wang, X., Mancuso, A., Gao, X., Wu, M., et al. (2011). p53 regulates biosynthesis through direct inactivation of glucose-6-phosphate dehydrogenase. Nat. Cell Biol. 13, 310–316.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kohn, A. D., Summers, S. A., Birnbaum, M. J., and Roth, R. A. (1996). Expression of a constitutively active Akt/PKB Ser/Thr kinase in 3T3-L1 adipocytes stimulates glucose uptake and glucose transporter 4 translocation. J. Biol. Chem. 271, 31372–31378.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kreeger, P. K., and Lauffenburger, D. A. (2010). Cancer systems biology: a network modeling perspective. Carcinogenesis 31, 2–8.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kucherenko, S., Rodriguez-Fernandez, M., Pantelides, C., and Shah, N. (2009). Monte Carlo evaluation of derivative-based global sensitivity measures. Reliab. Eng. Syst. Safety 94, 1135–1148.

CrossRef Full Text

Kuo, W., Lin, J., and Tang, T. K. (2000). Human glucose-6-phosphate dehydrogenase (G6PD) gene transforms NIH 3T3 cells and induces tumors in nude mice. Int. J. Cancer 85, 857–864.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lambeth, M. J., and Kushmerick, M. J. (2002). A computational model for glycogenolysis in skeletal muscle. Ann. Biomed. Eng. 30, 808–827.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Land, S. C., and Tee, A. R. (2007). Hypoxia-inducible factor 1alpha is regulated by the mammalian target of rapamycin (mTOR) via an mTOR signaling motif. J. Biol. Chem. 282, 20534–20543.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Le, A., Cooper, C. R., Gouw, A. M., Dinavahi, R., Maitra, A., Deck, L. M., et al. (2010). Inhibition of lactate dehydrogenase A induces oxidative stress and inhibits tumor progression. Proc. Natl. Acad. Sci. U.S.A. 107, 2037–2042.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Le Novère, N., Hucka, M., Mi, H., Moodie, S., Schreiber, F., Sorokin, A., et al. (2009). The systems biology graphical notation. Nat. Biotechnol. 27, 735–741.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lee, B. L., Kim, W. H., Jung, J., Cho, S. J., Park, J. W., Kim, J., et al. (2008). A hypoxia-independent up-regulation of hypoxia-inducible factor-1 by Akt/PKB contributes to angiogenesis in human gastric cancer. Carcinogenesis 29, 44–51.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Levine, A. J., and Puzio-Kuter, A. M. (2010). The control of the metabolic switch in cancers by oncogenes and tumor suppressor genes. Science 330, 1340.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Li, C., Donizelli, M., Rodriguez, N., Dharuri, H., Endler, L., Chelliah, V., et al. (2010a) BioModels Database: an enhanced, curated and annotated resource for published quantitative kinetic models. BMC Syst. Biol. 4, 92. doi:10.1186/1752-0509-4-92

CrossRef Full Text

Li, Y., Solomon, T. P. J., Haus, J. M., Saidel, G. M., Cabrera, M. E., and Kirwan, J. P. (2010b). Computational model of cellular metabolic dynamics: effect of insulin on glucose disposal in human skeletal muscle. Am. J. Physiol. Endocrinol. Metab. 298, 1198–1209.

CrossRef Full Text

Liebermeister, W., Uhlendorf, J., and Klipp, E. (2010). Modular rate laws for enzymatic reactions: thermodynamics, elasticities and implementation. Bioinformatics 26, 1528–1534.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lunt, S. Y., and Vander Heiden, M. G. (2011). Aerobic glycolysis: meeting the metabolic requirements of cell proliferation. Annu. Rev. Cell Dev. Biol. 27, 441–464.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Manning, B. D., and Cantley, L. C. (2007). Akt/PKB signaling: navigating downstream. Cell 129, 1261–1274.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Marín-Hernández, A., Gallardo-Pérez, J. C., Ralph, S. J., Rodríguez-Enríquez, S., and Moreno-Sánchez, R. (2009). HIF-1alpha modulates energy metabolism in cancer cells by inducing over-expression of specific glycolytic isoforms. Mini Rev. Med. Chem. 9, 1084–1101.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Marín-Hernández, A., Gallardo-Pérez, J. C., Rodríguez-Enríquez, S., Encalada, R., Moreno-Sánchez, R., and Saavedra, E. (2011). Modeling cancer glycolysis. Biochim. Biophys. Acta 1807, 755–767.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Marín-Hernández, A., Rodríguez-Enríquez, S., Vital-González, P. A., Flores-Rodríguez, F. L., Macías-Silva, M., Sosa-Garrocho, M., et al. (2006). Determining and understanding the control of glycolysis in fast-growth tumor cells. Flux control by an over-expressed but strongly product-inhibited hexokinase. FEBS J. 273, 1975–1988.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mathupala, S. P., Ko, Y. H., and Pedersen, P. L. (2009). Hexokinase-2 bound to mitochondria: cancer’s stygian link to the “Warburg effect” and a pivotal target for effective therapy. Semin. Cancer Biol. 19, 17–24.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mohanti, B. K., Rath, G. K., Anantha, N., Kannan, V., Das, B. S., Chandramouli, B. A., et al. (1996). Improving cancer radiotherapy with 2-deoxy-D-glucose: phase I/II clinical trials on human cerebral gliomas. Int. J. Radiat. Oncol. Biol. Phys. 35, 103–111.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mosca, E., Barcella, M., Alfieri, R., Bevilacqua, A., Canti, G., and Milanesi, L. (2012). Systems biology of the metabolic network regulated by the Akt pathway. Biotechnol. Adv. 30, 131–141.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Nicholson, K. M., and Anderson, N. G. (2002). The protein kinase B/Akt/PKB signalling pathway in human malignancy. Cell. Signal. 14, 381–395.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Novak, B., and Tyson, J. J. (1993). Numerical analysis of a comprehensive model of M-phase control in Xenopus oocyte extracts and intact embryos. J. Cell Sci. 106, 1153–1168.

Pubmed Abstract | Pubmed Full Text

Porporato, P. E., Dhup, S., Dadhich, R. K., Copetti, T., and Sonveaux, P. (2011). Anticancer targets in the glycolytic metabolism of tumors: a comprehensive review. Front. Pharmacol. 2:49. doi:10.3389/fphar.2011.00049

CrossRef Full Text

Pugh, C. W., and Ratcliffe, P. J. (2003). Regulation of angiogenesis by hypoxia: role of the HIF system. Nat. Med. 9, 677–684.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Rathmell, J. C., Fox, C. J., Plas, D. R., Hammerman, P. S., Cinalli, R. M., and Thompson, C. B. (2003). Akt/PKB-directed glucose metabolism can prevent Bax conformation change and promote growth factor-independent survival. Mol. Cell. Biol. 23, 7315–7328.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Reitzer, L. J., Wice, B. M., and Kennel, D. (1980). The pentose cycle. Control and essential function in HeLa cell nucleic acid synthesis. J. Biol. Chem. 255, 5616–5626.

Pubmed Abstract | Pubmed Full Text

Saltelli, A., Chan, K., and Scott, E. (2000). Sensitivity Analysis. Chichester: Wiley.

Sarbassov, D. D., Guertin, D. A., Ali, S. M., and Sabatini, D. M. (2005). Phosphorylation and regulation of Akt/PKB by the rictor-mTOR complex. Science 307, 1098–1101.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Sauro, H. (2011). Enzyme Kinetics for Systems Biology. Seattle: Ambrosius Publishing and Future Skill Software.

Schuster, R., and Holzhutter, H. G. (1995). Use of mathematical models for predicting the metabolic effect of large-scale enzyme activity alterations. Application to enzyme deficiencies of red blood cells. Eur. J. Biochem. 229, 403–418.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Semenza, G. L., Roth, P. H., Fang, H. M., and Wang, G. L. (1994). Transcriptional regulation of genes encoding glycolytic enzymes by hypoxia-inducible factor 1. J. Biol. Chem. 269, 23757–23763.

Pubmed Abstract | Pubmed Full Text

Sengupta, N., Rose, S. T., and Morgan, J. A. (2011). Metabolic flux analysis of CHO cell metabolism in the late non-growth phase. Biotechnol. Bioeng. 108, 82–92.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Singh, D., Banerji, A. K., Dwarakanath, B. S., Tripathi, R. P., Gupta, J. P., Mathew, T. L., et al. (2005). Optimizing cancer radiotherapy with 2-deoxy-D-glucose dose escalation studies in patients with glioblastoma multiforme. Strahlenther. Onkol. 181, 507–514.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Sobol, I. M. (1967). On the distribution of points in a cube and approximate evaluation of integrals. USSR Comp. Math. Math. Phys. 7, 86–112.

CrossRef Full Text

Song, M. S., Salmena, L., and Pandolfi, P. P. (2012). The functions and regulation of the PTEN tumour suppressor. Nat. Rev. Mol. Cell Biol. 13, 283–296.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Tennant, D. A., Durán, R. V., and Gottlieb, E. (2010). Targeting metabolic transformation for cancer therapy. Nat. Rev. Cancer 10, 267–277.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Tyson, J. J., Chen, K., and Novak, B. (2001). Network dynamics and cell physiology. Nat. Rev. Mol. Cell Biol. 2, 908–916.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Vander Heiden, M. G. (2011). Targeting cancer metabolism: a therapeutic window opens. Nat. Rev. Drug Discov. 10, 671–684.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Vivanco, I., and Sawyers, C. L. (2002). The phosphatidylinositol 3-Kinase Akt/PKB pathway in human cancer. Nat. Rev. Cancer 2, 489–501.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wolf, A., Agnihotri, S., Micallef, J., Mukherjee, J., Sabha, N., Cairns, R., et al. (2011). Hexokinase 2 is a key mediator of aerobic glycolysis and promotes tumor growth in human glioblastoma multiforme. J. Exp. Med. 208, 313–326.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wullschleger, S., Loewith, R., and Hall, M. N. (2006). TOR signaling in growth and metabolism. Cell 124, 471–484.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Yoeli-Lerner, M., Chin, Y. R., Hansen, C. K., and Toker, A. (2009). Akt/PKB/protein kinase b and glycogen synthase kinase-3beta signaling pathway regulates cell migration through the NFAT1 transcription factor. Mol. Cancer Res. 7, 425–432.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Zhao, F., Mancuso, A., Bui, T. V., Tong, X., Gruber, J. J., Swider, C. R., et al. (2010). Imatinib resistance associated with BCR-ABL upregulation is dependent on HIF-1alpha-induced metabolic reprograming. Oncogene 29, 2962–2972.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Zoncu, R., Efeyan, A., and Sabatini, D. M. (2011). mTOR: from growth signal integration to cancer, diabetes and ageing. Nat. Rev. Mol. Cell Biol. 12, 21–35.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Appendix

Differential Algebraic Equations

Here we report the list of equations that define the DAE system. J_x denotes the flux of reaction or, more generally, biochemical process x; V is the volume of the homogeneous compartment in which all biochemical processes take place.

$\begin{array}{llll}& \text{d AMP}/\text{dt}=0\hfill & & \\ & \text{d ATP}/\text{dt}=\left(-J_{\text{HK}}+J_{\text{PK}}-J_{\text{ACL}}+J_{\text{cAK}}-J_{\text{ATPase}}-J_{\text{PFK}}\right.\left.+J_{\text{PGK}}-J_{\text{PRPPS}}+J_{\text{MPM}}-J_{\text{GS}}\right)∕V\hfill & & \\ & \text{d BPG}/\text{dt}=\left(J_{\text{GAPDH}}-J_{\text{PGK}}\right)/V\hfill & & \\ & \text{d CIT}/\text{dt}=0\hfill & & \\ & \text{d DHAP}/\text{dt}=\left(J_{\text{FBA}}+J_{\text{TPI}}\right)/V\hfill & & \\ & \text{d E4P}/\text{dt}=\left(J_{\text{TAL}}-J_{\text{TKL2}}\right)/V\hfill & & \\ & \text{d F6P}/\text{dt}=\left(J_{\text{PGI}}-J_{\text{PFK}}+J_{\text{TKL2}}+J_{\text{TAL}}\right)/V\hfill & & \\ & \text{d F6P}/\text{dt}=\left(J_{\text{PGI}}-J_{\text{PFK}}+J_{\text{TKL2}}+J_{\text{TAL}}\right)/V\hfill & & \\ & \text{d F16P}/\text{dt}=\left(J_{\text{PFK}}-J_{\text{FBA}}\right)/V\hfill & & \\ & \text{d G1P}/\text{dt}=\left(J_{\text{GP}}-J_{\text{PGLM}}-J_{\text{GS}}\right)/V\hfill & & \\ & \text{d G6P}/\text{dt}=\left(J_{\text{HK}}-J_{\text{PGI}}-J_{\text{G6PDH}}+J_{\text{PGLM}}\right)/V\hfill & & \\ & \text{d GAP}/\text{dt}=\left(J_{\text{FBA}}-J_{\text{TPI}}-J_{\text{GAPDH}}+J_{\text{TKL}}+J_{\text{TKL2}}-J_{\text{TAL}}\right)/V\hfill & & \\ & \text{d GLC}/\text{dt}=\left(J_{\text{GLUT}}-J_{\text{HK}}\right)/V\hfill & & \\ & \text{d NAD}/\text{dt}=\left(J_{\text{LDH}}-J_{\text{GAPDH}}+J_{\text{DHase}}\right)/V\hfill & & \\ & \text{d NADP}/\text{dt}=\left(-J_{\text{G6PDH}}-J_{\text{PGDH}}+J_{\text{DPHase}}\right)/V\hfill & & \\ & \text{d Pi}/\text{dt}=0\hfill & & \\ & \text{d PEP}/\text{dt}=\left(J_{\text{ENO}}-J_{\text{PK}}\right)/V\hfill & & \\ & \text{d PG2}/\text{dt}=\left(J_{\text{PGYM}}-J_{\text{ENO}}\right)/V\hfill & & \\ & \text{d PG3}/\text{dt}=\left(J_{\text{PGK}}-J_{\text{PGYM}}\right)/V\hfill & & \\ & \text{d PGN}/\text{dt}=\left(J_{\text{G6PDH}}-J_{\text{PGDH}}\right)/V\hfill & & \\ & \text{d PYR}/\text{dt}=\left(J_{\text{PK}}-J_{\text{LDH}}-\text{1}∕y{J}_{\text{MPM}}\right)/V\hfill & & \\ & \text{d LAC}/\text{dt}=0\hfill & & \\ & \text{d R5P}/\text{dt}=\left(J_{\text{R5PI}}-J_{\text{TKL}}-J_{\text{PRPPS}}\right)/V\hfill & & \\ & \text{d RU5P}/\text{dt}=\left(J_{\text{PGDH}}-J_{\text{R5PI}}-J_{\text{RUPE}}\right)/V\hfill & & \\ & \text{d X5P}/\text{dt}=\left(J_{\text{RUPE}}-J_{\text{TKL}}-J_{\text{TKL2}}\right)/V\hfill & & \\ & \text{d S7P}/\text{dt}=\left(J_{\text{TKL}}-J_{\text{TAL}}\right)/V\hfill & & \end{array}$

Flux Expressions

For each reaction, the corresponding flux expression and the parameter values are listed below with the references to the works they were taken from. When both the maximum rates of the forward (V_f) and reverse (V_b) reactions appear explicitly in the flux expression (e.g., see PGI), V_b was adjusted in agreement with the constraint imposed by the Haldane relationship (Sauro, 2011)

$V_{\text{b}}=V_{\text{f}}\cdot K_{\text{mb}}/\left(K_{\text{eq}}\cdot K_{\text{mf}}\right)$

on the basis of the optimized value of V_f (Table 1), the reactant constants for the forward (K_mf) and backward (K_mb) reaction, and the equilibrium constant (K_eq).

Adenylate kinase

$\text{2ADP}=\text{ATP}+\text{AMP}$

Common modular rate law (Liebermeister et al., 2010).

$\begin{array}{llll}\text{J}\text{\_}\text{AK}& =V_{\text{f}}\cdot \text{AD}{\text{P}}^{\text{2}}\cdot \left(\text{1}-\text{ATP}\cdot \text{AMP}/K_{\text{eq}}\right)/\left[{\left(\text{1}+\text{ADP}\right)}^{\text{2}}+\left(\text{1}+\text{ATP}\right)\cdot \left(\text{1}+\text{AMP}\right)-\text{1}\right]\hfill & & \end{array}$

Consumption of ATP (ATPase)

$\text{ATP}\to \text{ADP}+\text{Pi}$

Common modular rate law (Liebermeister et al., 2010).

$\text{J}\text{\_}\text{ATPase}=V\text{f}\cdot \text{ATP}$

Oxidation of NADH (DHase)

$\text{NADH}=\text{NAD}+\text{H}$

Common modular rate law (Liebermeister et al., 2010).

$\begin{array}{llll}\text{J}\text{\_}\text{DHase}& =V_{\text{f}}\cdot \text{NADH}\cdot \left[\text{1}-\text{NAD}/\left(\text{NADH}\cdot K_{\text{eq}}\right)\right]∕\left[\left(\text{1}+\text{NADH}\right)+\left(\text{1}+\text{NAD}\right)-\text{1}\right]& & \end{array}$

Oxidation of NADPH (DPHase)

$\text{NADPH}=\text{NADP}+\text{H}$

Common modular rate law (Liebermeister et al., 2010).

$\begin{array}{llll}\text{J}\text{\_}\text{DPHase}& =V_{\text{f}}\cdot \text{NADPH}\cdot \left[\text{1}-\text{NADP}/\left(\text{NADPH}\cdot K_{\text{eq}}\right)\right]/\left[\left(\text{1}+\text{NADPH}\right)+\left(\text{1}+\text{NADP}\right)-\text{1}\right]& & \end{array}$

Enolase (ENO)

$\text{PG2}=\text{PEP}$

The rate equation is from Marín-Hernández et al. (2011).

$V_{\text{mr}}=V_{\text{mf}}/K_{\text{eq}}\cdot K_{\text{mp}}/K_{\text{ms}}$

Fructose bisphosphate aldolase

$\text{F16P}=\text{GAP}+\text{DHAP}$

Flux expression from Marín-Hernández et al. (2011).

$\begin{array}{ll}{V}_{\text{mr}}\hfill & =V_{\text{mf}}/K_{\text{eq}}\cdot \left(K_{\text{dhap}}\cdot K_{\text{g3p}}/K_{\text{fbp}}\right)\hfill \\ A\hfill & =\text{F16P};P=\text{DHAP};Q=\text{GAP};\hfill \\ {\text{J}}_{\text{FBA}}\hfill & =\left(V_{\text{mf}}\cdot A/K_{\text{fbp}}-V_{\text{mr}}\cdot P\cdot Q/\left(K_{\text{dhap}}\cdot K_{\text{g3p}}\right)\right)/\left(\text{1}+A/K_{{}_{\text{fbp}}}+P/K_{\text{dhap}}+Q/K_{\text{g3p}}+P\cdot Q/\left(K_{\text{dhap}}\cdot K_{\text{g3p}}\right)\right)\hfill \end{array}$

Glyceraldehyde-3-phosphate dehydrogenase

$\text{GAP}+\text{NAD}+\text{Pi}=\text{BPG}+\text{NADH}+\text{H}$

Flux expression from Marín-Hernández et al. (2011).

$\begin{array}{ll}\hfill A& =\text{NAD};B=\text{GAP};C=\text{Pi};P=\text{BPG};Q=\text{NADH};\hfill \\ \hfill V_{\text{mr}}& =V_{\text{mf}}/K_{\text{eq}}\cdot K_{\text{dpg}}\cdot K_{\text{nadh}}/\left(K_{\text{g3p}}\cdot K_{\text{nad}}\cdot K_{\text{p}}\right);\hfill \\ \text{J}\text{\_}\text{GAPDH}\hfill & =\left(V_{\text{mf}}\cdot A\cdot B\cdot C/\left(K_{\text{nad}}\cdot K_{\text{g3p}}\cdot K_{\text{p}}\right)\right.\left.-V_{\text{mr}}\cdot P\cdot Q/\left(K_{\text{dpg}}\cdot K_{\text{nadh}}\right)\right)/\left(\text{1}+A/K_{\text{nad}}+A\cdot B/\right.\left(K_{\text{nad}}\cdot K_{\text{g3p}}\right)\hfill \\ \hfill & +A\cdot B\cdot C/\left(K_{\text{nad}}\cdot K_{\text{g3p}}\cdot K_{\text{p}}\right)+P\cdot Q/\left.\left(K_{\text{dpg}}\cdot K_{\text{nadh}}\right)+Q/K_{\text{nadh}}\right)\hfill \end{array}$

Glucose-6-phosphate dehydrogenase

As in Holzhütter (2004) we used the sum of the consecutive reactions G6P + NADP + PGLT + NADPH + H and PGLT + PGN

$\text{G6P}+\text{NADP}=\text{PGN}+\text{NADPH}+\text{H}$

Flux expression from Holzhütter (2004).

$\begin{array}{l}\text{J\_G6PDH}=V_{\max }/K_{\text{G6P}}/K_{\text{NADP}}\cdot \left(\text{G6P}\cdot \text{NADP}-\text{PGN}\cdot \text{NADPH}/K_{\text{app}}\right)/(1+\text{NADP}\cdot \left(1+\text{G6P}/K_{\text{G6P}}\right)\\ /K_{\text{NADP}}+\text{ATP}/K_{\text{ATP}}+\text{NADPH}/K_{\text{NADPH}}+\text{BPG}/K_{\text{PGA23}})\end{array}$

Glucose transport

$\text{GLCe}=\text{GLC}$

Monosubstrate reversible Michaelis–Menten equation (Marín-Hernández et al., 2011).

$\begin{array}{llll}& J\text{\_}\text{glut}=V_{maxf}\cdot \left(\text{GLC}\text{\_}\text{e}-\text{GLC}∕k_{\text{eq}}\right)∕\left(K\text{Glc}\text{\_}\text{e}\cdot \left(\text{1}+\text{GLC}∕K_{\text{Glc}}\right)+\text{GLC}\text{\_}\text{e}\right)& & \end{array}$

Glycogen Phosphorylase A

$\text{GLYn}+\text{Pi}=\text{GLYm}+\text{G1P}$

Flux expression from Lambeth and Kushmerick (2002), where GLYn = GLYm = GLY.

$\begin{array}{ll}{V}_{maxr}\hfill & =V_{maxf}\cdot K_{\text{iGLYb}}\cdot K_{\text{iG1P}}∕\left(K_{\text{iGLY}f}\cdot K_{\text{Pi}}\cdot K_{\text{eq}}\right)\hfill \\ J\text{\_}\text{GPa}\hfill & =\left(V_{\text{max }f}\cdot \left(\text{GLY}\cdot \text{Pi}∕\left(K{\text{i}}_{\text{GLYf}}\cdot K\text{Pi}\right)\right)\right.\left.-V_{\text{max}r}\cdot \left(\text{GLY}\cdot \text{G1P}∕\left(K_{\text{GLYb}}\cdot K{\text{i}}_{\text{G1P}}\right)\right)\right)\hfill \\ & ∕\left(\text{1}+\text{GLY}∕K{\text{i}}_{\text{GLY}f}+\text{Pi}∕K\text{iPi}+\text{GLY}∕K{\text{i}}_{\text{GLYb}}+\text{G1P}∕K{\text{i}}_{\text{G1P}}\right)\left.+\text{GLY}\cdot \text{Pi}∕\left(K{\text{i}}_{\text{GLYf}}\cdot K\text{iPi}\right)+\left(\text{GLY}\cdot \text{G1P}\right)∕\left(K{\text{i}}_{\text{GLYb}}\cdot K{\text{i}}_{\text{G1P}}\right)\right)\hfill \end{array}$

Glycogen Phosphorylase B

$\text{GLYn}+\text{Pi}=\text{GLYm}+\text{G1P}$

Flux expression from Lambeth and Kushmerick (2002), where $\text{GLYn}=\text{GLYm}=\text{GLY}$.

$\begin{array}{l}J\_\text{GPb}=\left(V_{\max f}\cdot \left(\text{GLY}\cdot \text{Pi}/K{\text{i}}_{\text{GLYf}}\cdot K\text{Pi}\right)\right)-V_{\max r}\left(\text{GLY}\cdot \text{G1P}/\left(K{\text{i}}_{\text{GLYb}}\cdot K_{\text{G1P}}\right)\right))\\ /\left(1+\text{GLY}/K{\text{i}}_{\text{GLYf}}+\text{Pi}/K\text{iPi}+\text{GLY}/K{\text{i}}_{\text{GLYb}}+\text{G1P}/K{\text{i}}_{\text{G1P}}+\left(\text{GLY}\cdot \text{Pi}\right)/\left(K{\text{i}}_{\text{GLYf}}\cdot K\text{Pi}\right)+\left(\text{GLY}\cdot \text{G1P}\right)/\left(K{\text{i}}_{\text{GLYb}}\cdot K_{\text{G1P}}\right)\right)\\ \cdot \left({\text{AMP}}^{\text{nH}}/K_{\text{amp}}\right)/\left({\text{1+AMP}}^{\text{nH}}/K_{\text{amp}}\right);\end{array}$

Glycogen synthase

We used the sum of three reactions $\text{G1P}+\text{UTP}=\text{UDP}-\text{Glc}+\text{2Pi},\text{UDP}-\text{Glc}+\text{GLYn}=\text{UDP}+\text{GLYm and UDP}+\text{ATP}=\text{UTP}+\text{ ADP}$

$\text{GLY}+\text{G1P}+\text{ATP}=\text{GLY}+\text{ADP}+\text{2Pi}$

Flux expression from Li et al. (2010b).

$\begin{array}{ll}\text{GLYn}\hfill & =\text{GLYm}=\text{GLY}\hfill \\ J\text{\_}\text{GS}\hfill & =V_{\text{maxf}}∕K_{\text{f}}\cdot \text{G1P}\cdot \text{ATP}\cdot \text{GLY}\cdot \left(\text{1}-\text{P}{\text{i}}^{\text{2}}\cdot \text{ADP/(G1P}\cdot \text{ATP}\cdot K_{\text{eq}}\text{)}\right)∕\text{(1 + G1P}\cdot \text{ATP}\cdot \text{GLY/}{K}_{\text{f}}\text{ + GLY}\cdot \text{P}{\text{i}}^{\text{2}}\cdot \text{ADP/}{K}_{\text{r}}\text{)}\hfill \end{array}$

Hexokinase

$\text{GLC}+\text{ATP}=\text{ADP}+\text{G6P}$

Random bi-substrate Michaelis–Menten (Marín-Hernández et al., 2011).

$\begin{array}{ll}\hfill A& =\text{GLC};B=\text{ATP};P=\text{G6P};Q=\text{ADP}\hfill \\ J\text{\_}\text{HK}\hfill & =\left(V_{\text{mf}}∕\left(K_{\text{a}}\cdot K_{\text{b}}\right)\cdot \left(A\cdot B-P\cdot Q∕K_{\text{app}}\right)\right)∕\left(\text{1}+A∕K_{\text{a}}+B∕K_{\text{b}}+A\cdot B∕\left(K_{\text{a}}\cdot K_{\text{b}}\right)+P∕K_{\text{p}}+Q∕K_{\text{q}}+P\cdot Q\right.\hfill \\ \hfill & \left.∕\left(K_{\text{p}}\cdot K_{\text{q}}\right)+A\cdot Q∕\left(K_{\text{a}}\cdot K_{\text{q}}\right)+P\cdot B∕\left(K_{\text{p}}\cdot K_{\text{b}}\right)\right)\hfill \end{array}$

Lactate dehydrogenase

$\text{PYR}+\text{NADH}=\text{LAC}+\text{NAD}$

The rate equation is from Marín-Hernández et al. (2011).

$\begin{array}{ll}\hfill A& =\text{NADH},B=\text{PYR},P=\text{LAC},Q=\text{NAD}\hfill \\ \hfill V_{\text{mr}}& =V_{\text{mf}}∕K_{\text{eq}}\cdot K_{\text{p}}\cdot {\text{K}}_{\text{q}}∕\left(K_{\text{a}}\cdot K_{\text{b}}\right)\hfill \\ J_{\text{LDH}}\hfill & =\left(V_{\text{mf}}\cdot A\cdot B∕\left(\text{alfa}\cdot K_{\text{a}}\cdot K_{\text{b}}\right)-V_{\text{mr}}\cdot P\cdot Q∕\left(\text{beta}\cdot K_{\text{p}}\cdot K_{\text{q}}\right)\right)\hfill \\ \hfill & ∕\left(\text{1}+A∕K_{\text{a}}+B∕K_{\text{b}}+A\cdot B∕\left(\text{alfa}\cdot K_{\text{a}}\cdot K_{\text{b}}\right)+P\cdot Q∕\left(\text{beta}\cdot K_{\text{p}}\cdot K_{\text{q}}\right)+P∕K_{\text{p}}+Q∕K_{\text{q}}\right);\hfill \end{array}$

Mitochondrial pyruvate metabolism

Common modular rate law (Liebermeister et al., 2010).

$\text{1}∕y\text{ PYR}+\text{5}∕\left(\text{2}y\right)\text{O}2+\text{Pi}+\text{ADP}\to \text{3}∕y\text{ CO2}+\text{5}∕y{\text{H}}_{\text{2}}\text{O}+\text{ATP}$

$\begin{array}{l}{J}_{\text{M}}\text{PM}=V_{\text{f}}\cdot ({\text{PYR}}^{1/\text{y}}\cdot \text{Pi}\cdot \text{ADP}\cdot O_2^{5/\left(2\cdot y\right)}\cdot \left[1-\left(\text{ATP}\cdot {\text{CO}}_2^{\text{3}/y}\right)/\left({\text{PYR}}^{1/y}\cdot {\text{O}}_2^{5/\left(2\cdot y\right)}\cdot \text{Pi}\cdot \text{ADP}\cdot K\text{eq}\right)\right]\\ /\left[{\left(1+\text{PYR}\right)}^{1/y}\cdot {\left(1+{\text{O}}_2\right)}^{5/\left(2\cdot y\right)}\cdot \left(1+\text{Pi}\right)\cdot \left(1+\text{ADP}\right)+\left(1+\text{ADP}\right)\cdot {\left(1+{\text{CO}}_{\text{2}}\right)}^{3/y}-1\right]\end{array}$

Phosphofructokinase

$\text{F6P}+\text{ATP}=\text{F16P}+\text{ADP}$

Flux expression from Marín-Hernández et al. (2011).

$A=\text{F6P};B=\text{ATP};P=\text{F16P};Q=\text{ADP}.$

${\text{J}}_{\text{PFK}}=\left(\begin{array}{c}{V}_{\text{m}}\cdot \left(B/K_{\text{atp}}\right)/\left(1+\left(B/K_{\text{atp}}\right)\right)\cdot \left(1+\left(\text{beta}\cdot \text{F26P}\right)/\left(\text{alfa}\cdot \text{Kf26bp}\right)\right)/\left(1+\left(\text{F26P}/\left(\text{alfa}\cdot \text{Kf26bp}\right)\right)\right)\cdot \\ \left(A\cdot \left(1+\left(\text{F26P}/\left(\text{alfa}\cdot \text{Kf26bp}\right)\right)\right)\right)/\left(\text{Kf6p}\cdot \left(\text{1+}\left(\text{F26P}/\text{Kf26bp}\right)\right)\right)\cdot \left(1+\left(A\cdot \left(1+\left(\text{F26P}/\left(\text{alfa}\cdot \text{Kf26bp}\right)\right)\right)\right)\right)/\\ {\left(\text{Kf6p}\cdot \left(1+\left(\text{F26P}/\text{Kf26bp}\right)\right)\right)}^3/\left(L\cdot {\left(1+\text{CIT}/K_{cit}\right)}^4\cdot {\left(1+B/K_{\text{iatp}}\right)}^4\right)/{\left(1+\left(\text{F26P}/\text{Kf26bp}\right)\right)}^4\\ +\left(1+\left(A\cdot \left(1+\left(\text{F26P}/\left(\text{alfa}\cdot \text{Kf26bp}\right)\right)\right)\right)\right)/{\left(\text{Kf6p}\cdot \left(1+\left(\text{F26P}/\text{Kf26bp}\right)\right)\right)}^4-(\left(Q\cdot P\right)/\left(K_{\text{adp}}\cdot K_{\text{fbp}}\cdot K_{\text{app}}\right)/\\ \left(\text{Q}/K_{\text{adp}}\right)+\left(P/K_{\text{fbp}}\right)+\left(\left(Q\cdot P\right)/\left(K_{\text{adp}}\cdot K_{\text{fbp}}\right)+1\right)\end{array}\right)$

Phosphogluconate dehydrogenase

$\text{PGN}+\text{NADP}=\text{RU5P}+\text{CO2}+\text{NADPH}+\text{2H}$

Flux expression from Holzhütter (2004).

$\begin{array}{l}{J}_{\text{PGDH}}=V_{\max }/K_{\text{6PGI}}/K_{\text{NADP}}\cdot \left(\text{PGN}\cdot \text{NADP}-\text{RU5P}\cdot \text{NADPH}/K_{\text{app}}\right)/(\left(1+\text{NADP}/K_{\text{NADP}}\right)\\ \cdot \left(1+\text{PGN}/K_{\text{6PG1}}+\text{BPG}/K_{\text{PGA23}}\right)+\text{ATP}/K_{\text{ATP}}+\text{NADPH}\cdot \left(\text{1+}\text{PGN}/K_{\text{6PG2}}\right)/K_{\text{NADPH}})\end{array}$

Phosphoglucoisomerase

Monoreactant reversible equation with competitive inhibition by E4P, 6PG, and FBP (Marín-Hernández et al., 2011).

$\begin{array}{ll}\text{G6P}\hfill & =\text{F6P}\hfill \\ V_{\text{mr}}\hfill & =V_{\text{mf}}\cdot K_{\text{f6p}}∕\left(K_{\text{eq}}\cdot K_{\text{g6p}}\right)\hfill \\ J_{\text{PGI}}\hfill & =\left(V_{\text{mf}}\cdot \text{G6P}∕K_{\text{g6p}}-V_{\text{mr}}\cdot \text{F6P}∕K_{\text{f6p}}\right)∕\left(\text{1}+A∕K_{\text{g6p}}+P∕K_{\text{f6p}}+\text{E4P}∕K_{\text{ery4p}}+\text{F16P}∕K_{\text{fbp}}+\text{PGN}∕K_{\text{pg}}\right)\hfill \end{array}$

Phosphoglucomutase

$\text{G1P}=\text{G6P}$

Flux expression from Lambeth and Kushmerick (2002).

$\begin{array}{ll}{V}_{\text{maxr}}\hfill & =V_{\text{maxf}}\cdot K_{\text{G6P}}∕\left(K_{\text{G1P}}\cdot K_{\text{eq}}\right)\hfill \\ J_{\text{PGLM}}\hfill & =\left(V_{\text{maxf}}\cdot \text{G1P}∕K_{\text{G1P}}-V_{\text{maxr}}\cdot \text{G6P}∕K_{\text{G6P}}\right)∕\left(\text{ 1}+\text{G1P}∕K_{\text{G1P}}+\text{G6P}∕K_{\text{G6P}}\right)\hfill \end{array}$

Phosphoglycerate kinase

$\text{BPG}+\text{ADP}=\text{PG3}+\text{ATP}$

Flux expression from Marín-Hernández et al. (2011).

$\begin{array}{ll}\hfill A& =\text{BPG},B=\text{ADP},P=\text{PG3},Q=\text{ATP}\hfill \\ \hfill V_{\text{mr}}& =V_{\text{mf}}∕K_{\text{eq}}\cdot K_{\text{p}}\cdot K_{\text{q}}∕\left(K_{\text{a}}\cdot K_{\text{b}}\right)\hfill \\ J_{\text{PGK}}\hfill & =\left(V_{\text{mf}}\cdot A\cdot B∕\left(\text{alfa}\cdot K_{\text{a}}\cdot K_{\text{b}}\right)-V_{\text{mr}}\cdot P\cdot Q∕\left(\text{beta}\cdot K_{\text{p}}\cdot K_{\text{q}}\right)\right)\hfill \\ \hfill & ∕\left(\text{1}+A∕K_{\text{a}}+B∕K_{\text{b}}+A\cdot B∕\left(\text{alfa}\cdot K_{\text{a}}\cdot K_{\text{b}}\right)+P\cdot Q∕\left(\text{beta}\cdot K_{\text{p}}\cdot K_{\text{q}}\right)+P∕K_{\text{p}}+Q∕K_{\text{q}}\right)\hfill \end{array}$

Pyruvate kinase

$\text{PEP}+\text{ADP}=\text{PYR}+\text{ATP}$

The rate equation is from Marín-Hernández et al. (2011).

$\begin{array}{l}{\text{J}}_{\text{PK}}=V_{\text{m}}\cdot \left(B/K_{\text{adp}}\right)/\left(1+\left(B/K_{\text{adp}}\right)\right)\cdot \left(\left(A/K_{\text{pep}}\right)\cdot {\left(1+\left(A/K_{\text{pep}}\right)\right)}^3\right)/(\left(\left(L\cdot {\left(1+\left(Q/K_{\text{iatp}}\right)\right)}^4\right)/{\left(1+\left(\text{F16P}/K_{\text{fbp}}\right)\right)}^4\right)\\ +{\left(1+\left(A/K_{\text{pep}}\right)\right)}^4)-\left(\left(Q\cdot P\right)/\left(K_{\text{atp}}\cdot K_{\text{pyr}}\cdot K_{\text{app}}\right)\right)/\left(\left(\left(Q/K_{\text{atp}}\right)+\left(P/K_{\text{pyr}}\right)\right)+\left(\left(Q\cdot P\right)/\left(K_{\text{atp}}\cdot K_{\text{pyr}}\right)\right)+1\right));\end{array}$

Phosphoribosylpyrophosphate synthetase

$\text{R5P}+\text{ATP}=\text{AMP}+\text{PRPP}$

Flux expression from Holzhütter (2004).

$\begin{array}{c}{J}_{\text{PRPPS}}=V_{\text{max}}\cdot \left(\text{R5P}\cdot \text{ATP}-\text{PRPP}\cdot \text{AMP}∕K_{\text{app}}\right)∕\left(\left(K_{\text{ATP}}+\text{ATP}\right)\cdot \left(K_{\text{R5P}}+\text{R5P}\right)\right)\end{array}$

Phosphoglycerate mutase

$\text{PG3}=\text{PG2}$

The rate equation is from Marín-Hernández et al. (2011).

$\begin{array}{ll}\hfill V_{\text{mr}}& =V_{\text{mf}}∕K_{\text{eq}}\cdot K_{\text{mp}}∕K_{\text{ms}}\hfill \\ J_{\text{PGYM}}\hfill & =\left(V_{\text{mf}}\cdot \text{PG3}∕K_{\text{ms}}-V_{\text{mr}}\cdot \text{PG2}∕K_{\text{mp}}\right)∕\left(\text{1}+\text{PG3}∕K_{\text{ms}}+\text{PG2}∕K_{\text{mp}}\right);\hfill \end{array}$

Ribose-5-phosphate isomerase

$\text{RU5P}=\text{R5P}$

Flux expression from Schuster and Holzhutter (1995).

$\begin{array}{c}{J}_{\text{R5PI}}=V_{\text{max}}\cdot \left(\text{RU5P}-\text{R5P}∕{K_{\text{eq}}}_{\text{R5PI}}\right)∕\left(\text{RU5P}+K_{\text{Ru5P}}\cdot \left(\text{1}+\text{R5P}∕K_{\text{R5P}}\right)\right)\end{array}$

Ribulose-phosphate epimerase

$\text{RU5P}=\text{X5P}$

Flux expression from Holzhütter (2004).

$\begin{array}{c}{J}_{\text{RUPE}}=V_{\text{max}}\cdot \left(\text{RU5P}-\text{X5P}∕{K_{\text{eq}}}_{\text{RUPE}}\right)∕\left(\text{RU5P}+K_{\text{Ru5P}}\cdot \left(\text{1}+\text{X5P}∕K_{\text{X5P}}\right)\right)\end{array}$

Transaldolase

$\text{S7P}+\text{GAP}=\text{F6P}+\text{E4P}$

Flux expression from Schuster and Holzhutter (1995).

$\begin{array}{l}{J}_{\text{TAL}}=V_{\max }\cdot \left(\text{S7P}\cdot \text{GAP}-\text{E4P}\cdot \text{F6P}/K_{\text{eqTAL}}\right)/(\left(K_1+\text{GAP}\right)\cdot \text{S7P}+\left(K_2+K_6\cdot \text{F6P}\right)\cdot \text{GAP}+\left(K_3+K_5\cdot \text{F6P}\right)\cdot \text{E4P}\\ +K_4\cdot \text{F6P}+K_7\cdot \text{S7P}\cdot \text{E4P})\end{array}$

Transketolase

$\text{R5P}+\text{X5P}=\text{GAP}+\text{S7P}$

Flux expression from Schuster and Holzhutter (1995).

$\begin{array}{l}{J}_{\text{TKL}}=V_{\max }\cdot \left(\text{R5P}\cdot \text{X5P}-\text{GAP}\cdot \text{S7P}/K_{\text{eqTKL}}\right)/(\left(K_1+\text{R5P}\right)\cdot X5\text{P}+\left(K_2+K_6\cdot \text{S}7\text{P}\right)\cdot \text{R5P}+\left(K_3+K_5\cdot \text{S}7\text{P}\right)\cdot \text{GAP}\\ +K_4\cdot \text{S7P}+K_7\cdot \text{X5P}\cdot \text{GAP})\end{array}$

Transketolase 2

$\text{X5P}+\text{E4P}=\text{GAP}+\text{F6P}$

Flux expression from Schuster and Holzhutter (1995).

$\begin{array}{l}{J}_{\text{TKL2}}=V_{\max }\cdot \left(\text{E4P}\cdot \text{X5P}-\text{GAP}\cdot \text{F6P}/K_{\text{eqTKL2}}\right)/(\left(K_1+\text{E4P}\right)\cdot \text{X}5\text{P}+\left(K_2+K_6\cdot \text{F}6\text{P}\right)\cdot \text{E4P}+\left(K_3+K_5\cdot \text{F}6\text{P}\right)\cdot \text{GAP}\\ +K_4\cdot \text{F6P}+K_7\cdot \text{X5P}\cdot \text{GAP})\end{array}$

Triose-phosphate isomerase

$\text{GAP}=\text{DHAP}$

Flux expression from Marín-Hernández et al. (2011).

$\begin{array}{ll}{V}_r\hfill & =V_f/K_{\text{eq}}\cdot K_{\text{mp}}/K_{\text{ms}};\hfill \\ {\text{J}}_{\text{TPI}}\hfill & =\left(V_f\cdot \text{GAP}/K_{\text{ms}}-V_{\text{r}}\cdot \text{DHAP}/K_{\text{mp}}\right)/\left(1+\text{GAP}/K_{\text{ms}}+\text{DHAP}/K_{\text{mp}}\right)\hfill \end{array}$

General Model Parameters

Metabolite Concentrations at Steady States H and L

Metabolic Fluxes in Condition H and L