A Mathematical Model of HIF-1 Regulated Cellular Energy Metabolism

In this study, we formulated a mathematical model of hypoxia-inducible factor 1 (HIF-1) mediated regulation of cellular energy metabolism describing the reprogramming of cell metabolic processes from oxidative phosphorylation to glycolysis under reduced oxygen levels. The model considers the dynamics of fifteen biochemical species and the proton concentration with the underlying reaction processes localized in three intracellular compartments, i.e. the cytoplasm, mitochondrion and nucleus. More than sixty parameters of the model were calibrated using both the published data and the system steady-state based identification procedure. The model was validated by generating predictions which could be compared to empirical observations. The model behaviors representing the cell metabolism switching over in response to transitioning from a normoxic to hypoxic environment are consistent with the current views of the role of HIF-1 in hypoxia.


Introduction
Energy metabolism drives the functioning of organs and individual cells. There is a growing need in controlling specific metabolic processes that are skewed in many pathologies, such as cancer, infectious diseases, and sepsis [9,19,25]. Hypoxia (low oxygen tension) is an essential hallmark of the pathological processes. The metabolic reprogramming of cells that occurs in hypoxia is regulated by hypoxia-inducible factor 1 (HIF-1) [17,[21][22][23]. The activity of HIF-1 is determined by the oxygen level. In turn, HIF-1 regulates the balance between oxidative and glycolytic metabolism [21]. With oxygen as the substrate, the respiratory chain paves the way for the efficient generation of ATP from citric acid cycle products. Lowered oxygen levels lead to increased production of reactive oxygen species by respiratory complexes and without further protection, this may lead to cellular dysfunction or cell death. HIF acts as a master regulator adapting the metabolism to such hypoxic conditions by initiating several mechanisms which reduce mitochondrial activity and amplify anaerobic ATP production [21]. Although ATP levels cannot be maintained at normoxic values, the cell benefits from this switch, because ROS production in the mitochondria is prevented [32].
Understanding the complex regulation of cellular energy metabolisms under normal and hypoxic settings requires the development of mechanistic mathematical models of cellular energetics. Some mathematical models of intracellular respiration have already been proposed [1,2,10,15,28]. These models mostly describe the biophysical processes of oxidative phosphorylation but do not touch upon HIF-mediated regulation. Recently, two studies addressed the transcriptional activity of the HIF-1 signaling network [16] and the interleukin-15 dependent regulation of HIF-1a in Natural Killer cells [6]. As a mathematical model of HIF-1 regulated reprogramming of cellular metabolism does not yet exist, we aim to formulate and calibrate such an integrative model.
In this study, we formulate and calibrate the mathematical model of HIF-1 mediated regulation of cellular energy metabolism describing the reprogramming of cell metabolic processes from oxidative phosphorylation to glycolysis under reduced oxygen levels. The model considers biochemical and biophysical processes underlying the dynamics of sixteen chemical species and the proton concentration localized in the cytoplasm, mitochondrion and nucleus compartments. We identify unknown parameters such that the system exhibits a stable stationary state under normoxic conditions (see Section 4.1). The stability of the steady state is shown by computer experiments (simulation of the evolution of the metabolites for perturbed (±10%) steady states), see also Fig. 2, and by a linear stability analysis. The existence and uniqueness of a global solution to the model equations is shown in Section 3. Comparing numerical simulations under hypoxic conditions with and without HIF-mediated adaption confirms the ability of HIF to indirectly limit respiratory chain activity (by inhibiting the pyruvate dehydrogenase) and to increase glycolytic ATP production (by activating glycolytic enzymes), see Figs. 4, 5, and 6 in Section 4.

The Model
The mathematical model is given by a system of ordinary differential equations (ODEs) and consists of the main parts of the cellular energy metabolism, namely glycolysis and lactic acid fermentation in the cytoplasm as well as citric acid cycle and oxydative phosphorylation in the mitochondria, as summarized in Fig. 1. In addition, a HIF-mediated activation of the glycolysis and inactivation of the citric acid cycle at low oxygen concentrations is implemented. Reaction rates based on the law of mass action are preferably used to model the biochemical reactions. Oxydative phosphorylation is modeled using phenomenological rate equations. To keep the extent of the model to a minimum, properties of enzymatically driven reactions such as inhibition, saturation or other regulatory mechanisms are not considered for large parts of the model and subsequent reactions are often interpreted as one entity. Irreversibility of according reaction rates is often justified by the irreversibility of at least one of the underlying reactions. To make the paper more accessible for readers who Fig. 1 Biochemical scheme of the mathematical model of the cellular energy metabolism describing the species, metabolic reaction networks and corresponding regulatory mechanisms. Number i at reaction arrows refers to the corresponding rate equation v i . A simplified fatty acid catabolism is outlined in pale gray, but not yet considered in the model are not familiar with the biomedical application, we described all species considered in our model together with the corresponding acronym in Table 1. All parameters of the model are listed with their physiological meaning in Table 2.

The Cytoplasm
Glycolysis. Glucose uptake from the extracellular medium is modeled by a constant rate, assuming that the cell is always provided with sufficient extracellular glucose and intracellular processes do not affect its transport across the membrane: Next, glucose is further processed at the expense of two molecules of ATP, in order to prepare it for the following degradation process. Repeated phosphorylation and isomerization on the way from glucose to fructose 1,6-bisphosphate is summarized as one mechanism in the second rate equation: The subsequent cleavage of fructose 1,6-bisphosphate into two three-carbon units and their conversion into pyruvate is considered as the second and last step of glycolysis in the model. It is accompanied by the formation of four molecules ATP and two molecules NADH, making glycolysis an overall ATP-producing pathway.
F ructose 1,6-bP To obtain this rate equation, two more assumptions are made: Firstly, the amount of adenosine bodies is a conserved quantity in the model and adenosine only exists as ATP or ADP (A tot cyt = [AT P ] cyt + [ADP ] cyt ) and secondly, there is always a sufficient and constant amount of inorganic phosphate such that its concentration does not affect the velocity of the reaction.
Lactic Acid Fermentation. Analogously to the adenosine bodies, the conservation N tot cyt = [NADH ] cyt + [NAD] cyt is considered. While almost all energy consuming processes in the cell use ATP as an energy source and therefore regenerate ADP, the energy of NADH cannot be used directly and the cell would run out of NAD without further metabolic pathways. Therefore, lactic acid fermentation enables NAD to be regained in the cytoplasm, as pyruvate and NADH are converted into lactate and NAD.
. Reversibility of the reaction is not taken into account here and so the only possible fate of lactate is to be transported out of the cell: ATP Consumption. All energy consuming processes in each compartment of the cell are represented by the following rate equation, which was adopted from a model of Koenig et al. [13]: It consists of Michaelis-Menten type kinetics, accounting for basal life-sustaining processes and a linear term, which respects the ability of the cell to adjust its energy demand according to the current ATP supply.

Metabolite Transport Across the Inner Mitochondrial Membrane
Pyruvate Transport An alternative fate of cytosolic pyruvate is to be transported into mitochondria, where it can fuel the citric acid cycle. The transport is assumed to be reversible in the model, because mitochondrial pyruvate could otherwise accumulate arbitrarily when the citric acid cycle is inhibited. The direction of the transport is prescribed by the difference in particle numbers directly at the membrane. As the model does not resolve the distribution of pyruvate spatially, the respective metabolite concentration is divided by the volume of the membrane and then compared. Depending on the compartment, the resulting value is converted into cytosolic or mitochondrial volume units: wherek ANT := k ANT V 2 m . K eq,ANT is introduced such that ATP is transported out of the mitochondria under normal conditions.

The Mitochondrion
The model for the mitochondrial energy metabolism is adopted from Nazaret et al. [15] with extensions regarding the oxygen dependence of the oxydative phosphorylation. Citric Acid Cycle An additional preliminary step of the citric acid cycle is the oxydative decarboxylation of mitochondrial pyruvate into acetyl-CoA. It is the first occurrence of mitochondrial NADH production, one of the main objectives of the mitochondrial energy metabolism, and also an important regulatory point.
The concentration of Coenzyme A is assumed to be constant and without impact on the reaction rate. A round of the citric acid cycle starts, when the acetyl group of acetyl-CoA is transferred into oxaloacetate, forming citrate.
Oxaloacetate + Acetyl-CoA The reaction from citrate to α-ketoglutarate via isocitrate with NADH production is simplified to: α-ketoglutarate undergoes further transformations and finally, oxaloacetate is formed to complete the circle. During this reaction sequence, two molecules NADH, one molecule FADH 2 and one molecule GTP are gained. For the model, only NADH is considered and GTP is replaced by ATP as they are energetically equivalent.
The model respects some additional features of the TCA cycle. α-ketoglutarate and oxaloacetate are also linked via a reversible reaction within the scope of the amino acid metabolism. For the sake of simplicity, additional reactants are omitted: An output of the TCA cycle is modeled by the concentration dependent degradation of oxaloacetate in order that the system can reach a steady state.
Finally, the ATP dependent transformation of pyruvate into oxaloacetate complements the citric acid cycle by an replenishing reaction: Oxydative Phosphorylation Energy is now intermediately stored as NADH and the mitochondrial regeneration of NAD is the task of the respiratory chain. It comprises four enzyme complexes which are located in the impermeable inner mitochondrial membrane. The high energy electrons of NADH are transferred sequentially onto oxygen and the energy that is released during this process is used to translocate ten protons per NADH molecule from the mitochondrial matrix into the intermembrane space. When the matrix and the intermembrane space, separated by the inner mitochondrial membrane, are interpreted as a capacitor, the resulting proton difference can be expressed as a potential. In accordance with Nazaret et al., the proton concentration outside the mitochondrial matrix is assumed to be constant and therefore the temporal change of the membrane potential is linearly linked to the change of the intramitochondrial proton concentration: There are several processes which affect this concentration, firstly the respiratory chain and the ATPase. The latter couples the re-entrance of protons from the intermembrane space into the matrix to the synthesis of ATP. Instead of modeling the complex mechanisms at the respiratory chain in detail, a phenomenological approach is chosen. The four enzyme complexes are assembled into one rate equation, which exhibits a Michaelis-Menten type dependence in the NADH concentration as well as a complete inhibition once the membrane potential becomes too high to translocate more protons. For details, see [15]. For our purposes, this reaction rate is extended by a Michaelis-Menten type term for the oxygen dependence of the respiratory chain: .
The formation of ATP from ADP and P i , together with the backflow of protons into the mitochondrial matrix, is reversible: To gain a phenomenological, reversible rate equation, Nazaret et al. consider the overall Gibb's free energy of the coupled reaction: [ADP ] eq [P i ] eq . From the ansatz G = 0, they derive a formula for the critical ATP concentration in the equilibrium, dependent on the membrane potential. The difference between the actual ATP concentration and the critical ATP concentration then prescribes the direction of the reaction in the following phenomenological rate: For more details about the derivation, see again [15].

Proton Leak
Since the inner mitochondrial membrane is not completely impermeable for protons, some of them flow back into the matrix without contributing to the synthesis of ATP. This process is taken into account by a linear dependence on the membrane potential: Glycerol Phosphate Shuttle The glycerol phosphate shuttle for the feeding of high energy electrons from cytosolic NADH is assumed to have the same phenomenological properties as the respiratory chain: Michaelis-Menten like dependencies on the oxygen and cytosolic NADH concentration as well as an inhibition at high membrane potentials.
. We do not model the transport of oxygen from the blood into the mitochondria but interpret [O 2 ] as the mitochondrial oxygen concentration and set it manually to a constant value for the simulations.

Hypoxia Inducible Factor 1
The HIF Pathway. HIF1, assembled from the two subunits HIF-1α and HIF-1β, acts as a master regulator for the adaption of cellular metabolism to low oxygen levels. It is a transcription factor and initiates the enhanced transcription of a large number of different genes. For this model, the increased transcription of glycolysis-related genes is of interest, as well as the HIF-mediated transcription of the pyruvate dehydrogenase kinase. The latter can phosphorylate and thus inactivate the enzyme which catalyzes the oxydative decarboxylation in the mitochondria, resulting in a lowered citric acid cycle activity.
While HIF-1β is an abundant protein, HIF-1α is not present in relevant amounts during normoxia. Instead, the HIF-1α-gene is constantly transcribed and translated, but subsequently degraded in an oxygen-dependent manner. Only when oxygen concentrations are sufficiently low, HIF-1α can accumulate and unite with HIF-1β to affect the metabolism.
The mathematical model for the HIF accumulation during hypoxia is partly based on the a dynamic model of the HIF-1α network by Nguyen et al. [16]. A first simplification is to merge HIF-1α and HIF-1β into one entity, HIF1. The permanent and oxygen-independent formation of HIF1 is modeled by a constant rate: HIF1, like all proteins, is then subject to an unspecific protein turnover: HIF1-prolyle-hydroxylases (HIF-PHDs), the primary oxygen sensors of the cell, mark the HIF1 protein for degradation by the proteasomal machinery, when the oxygen concentration is sufficiently high: A negative feedback loop for HIF is incorporated via the modeling of the prolyl hydroxylase concentration. As the PHD-gene is also a target gene of HIF, the corresponding enzyme concentration is expected to rise during hypoxia: This guarantees the dampening of the HIF concentration after some time and also promotes a fast return to the HIF-independent metabolism once the cell is re-oxygenated.

Updated Rate Equations For HIF1-Dependent Reactions
We model the effect of high HIF levels on the glycolytic rates and on the rate of oxadative decarboxylation by imposing HIF-dependent activation/inactivation terms in the following way: For high HIF levels, the additional term approaches the value prescribed by parameter P i .

The ODE System
For the ODE system, we use the following notation:

Existence And Uniqueness Of Global Solutions
Proposition 1 We consider the initial value problem (IVP) consisting of (18)-(33) together with initial values x i (0) = x i,0 , i = 1, . . . , 16 which are non-negative except for component 14 (which represents the mitochondrial membrane potential) and satisfy The solution remains non-negative except for component x 14 and for all t ∈ (0, ∞) it holds that Proof Existence and uniqueness of a local solution on the maximal time interval follows due to the local Lipschitz-continuity of the right-hand side of (18)- (33). Next, we prove that the solution exists globally to the right. Let us consider the set First, we show that D is positive invariant for the IVP corresponding to (18)-(33), i.e., for all x 0 ∈ D, the solution x to the IVP satisfies x(t) ∈ D, for all t ≥ 0 belonging to the maximal existence interval. This follows by using e.g., Appendix 1, Theorem in [24] or Theorem 5.6.1 in [14] (see also [20]) and the following properties of the right-hand side of the system (18) To prove that the remaining species x i , i = 1, . . . , 16, i = 5, 6, 12, 13, 14 are bounded from above on the maximal existence interval, we set Thus, Gronwall's lemma implies that S cannot be infinite in finite time. Finally, we have to study the behavior of x 14 . We have The estimates

Numerical Results
Computer experiments allow simulating processes which are not (yet) accessible by biological experimental methods and may help to deepen the understanding-or to predict the behavior-of a system under certain conditions. Prior to that, the underlying model has to be calibrated, i.e. the parameters of the model have to be set to physically meaningful values. This is often an inverse problem, where parameters have to be determined such that the system reproduces some observed data. Unfortunately, appropriate data on the time course of metabolites, in particular under hypoxic conditions, is, to our knowledge, barely available. Therefore, we are restricted to stationary information and have to collect the metabolite concentrations from the literature with the aim of representing a metabolic equilibrium state (see Table 1). The parameters are then determined such that the model is able to reproduce this prescribed equilibrium.
Parameters which are not rate constants (such as for instance parameters of HIF-regulated activation/inactivation or K M -values) are either obtained from the literature and previous modeling studies or fixed manually to biologically reasonable values (see Table 2). Rate constants k AT P and k resp are adopted from [15].

Parameter Identification
The identification of the unknown rate constant vector k = (k 1 , . . . , k 14 ,k ANT , k GP S , k leak ) T ∈ R 17 of the system defined by (18)-(31) is based on the assumption that the metabolite concentrations are constant over time in the cell under normal conditions. Hence we look for a k such that the model reproduces the biological equilibrium which is prescribed by the steady state metabolite concentrations in Table 1. Thus, k has to be chosen in such a way that the righthand side of the ODE system (18)-(31) vanishes when the metabolite concentrations are at their equilibrium values. This leads to the following constrained optimization problem:  Membrane potential Δ Mitochondrion 150 mV [15] Hypoxia-inducible factor HIF 5 · 10 −6 mM [16] Prolyl-hydroxylase PHD 1 · 10 −5 mM [16] wherex = (x 1 (t), . . . , x 14 (t)) T andf is the right-hand side of the ODE system (18)- (31). Note that the parameters k 15 − k 19 of the HIF-subsystem are determined separately and that the HIF concentration is for now set constant to its basal value in Table 1. Constraint (35) guarantees non-negativity of the rate constants and constraints (36)-(37) lead to reaction rates, which are all within one order of magnitude. We use the default interiorpoint-algorithm of the MATLAB-function fmincon to solve the COP numerically and obtain the rate constants in Table 2. It remains to determine the parameters of the HIF-subsystem, which is composed of (32) and (33). The K M -values for v 17 are adopted directly from Nguyen et al. ( [16]). The rate constants of v 15 -v 19 are chosen manually (cf. Table 2) such that the HIF-subsystem exhibits the following characteristic properties: 1) The system is in an equilibrium for the basal concentrations of HIF1 and PHD (see Table 1) during normoxia. 2) HIF1 protein accumulation during hypoxia takes place in the time range of hours ( [11,27]). 3) Hypoxia-induced increase of HIF1-levels is attenuated after some time and HIF1-levels eventually even decrease again during prolonged hypoxia ( [11,27]). 4) A negative feedback loop between HIF1 and PHD promotes rapid return to the basal HIF1-concentration after reoxygenation ( [7,16]).

Numerical Simulations Under Normoxic Conditions
In the following, all numerical simulations of the initial value problem for the ODE-system are performed with the MATLAB-function ode13s (https://www.mathworks.com/). Using the parameters identified in Section 4.1, see Table 2, we parameterize the system (18)- (33) and obtain the stationary state given in Table 1. To examine the local stability of this steady state, the Jacobian matrix of the system is computed by using symbolic differentiation in   [16] MATLAB (https://www.mathworks.com/). The eigenvalues of the Jacobian evaluated in the steady state are determined numerically and without exception exhibit a negative real part. Thus we conclude exponential stability. This stability property is visualized in Fig. 2, where the initial values of all species are chosen ±10% away from the respective steady state values. The system eventually moves back to the equilibrium after some damped oscillations. Lactate takes the longest time to approach its steady state value, but as it does not interact with other metabolites, this does not influence the rest of the system. The system is robust to a severe drop in the ATP concentration. Indeed, when the ATP levels in the Cytosol and Mitochondrion compartments are reduced to 50% of their equilibrium values as shown in Fig. 3, the steady state of the metabolic system is also rapidly recovered. Mitochondrial ATP levels rise on the expense of the membrane potential and since ATP is transported from mitochondria into the cytosol during the simulations, the cytosolic ATP levels reach their steady state value even faster than the mitochondrial ATP levels. This is accompanied by some damped oscillations in the mitochondria and in the cytosol. Note that, if not indicated otherwise, the concentration curves in all plots are divided by their respective steady-state concentration (cf.

Numerical Simulations Under Hypoxic Conditions
The following computer experiments simulate hypoxic conditions by imposing an oxygen concentration of 0.05μM to the system (the normoxic cellular oxygen concentration is set to be 100μM).

Computer Experiments Without HIF-Mediated Regulatory Mechanisms
In the simulation presented in Fig. 4, the HIF mechanism is not considered, i.e., meaning that there is no additional activation of the glycolysis or inhibition of oxydative decarboxylation due to HIF accumulation. After a transient phase, the system approaches a new steady-state (apparent for longer simulation times, not entirely shown) with decreased but still relatively high ATP values (47% of the equilibrium value in the cytosol, 15% of the equilibrium value in the mitochondria). Maintenance of such high ATP concentrations is possible, because the reaction rates of the respiratory chain and of the ATPase are still at 50% and 40% of their equilibrium activity. Although the Michaelis-Menten-like term for the oxygen dependence of the respiratory chain (see (16)) reduces from approximately 1 to 0.05, the only slightly decreased membrane potential and the increased NADH levels as seen in Fig. 4b) enable the relatively high value of the respiratory rate. This guarantees a sufficient membrane potential for the perpetual synthesis of ATP within the scope of oxydative phosphorylation. Accumulation of glucose can be observed in the cytosol (cf. Fig. 4a)), because its transformation to fructose 1,6-bP is limited due to the low ATP level. Increased cytosolic NADH concentration (follows from NAD decrease) amplifies the lactic acid fermentation, resulting in lactate accumulation and cytosolic pyruvate depletion. The latter in turn effects the time-delayed decline of mitochondrial pyruvate, after its concentration initially rises due to the lack of reaction partners for further breakdown.

Computer Experiments With HIF-mediated Regulatory Mechanisms
When the same experiment is repeated in Fig. 5, but this time with an active HIF-pathway, changes in the metabolite concentrations can be observed after 2-6 hours of simulated time. The HIF-mediated inhibition of the oxydative decarboxylation and amplification of lactic acid fermentation together re-establish high NAD concentrations in both compartments. Therefore the respiratory chain remains at less than 1% of its equilibrium activity and subsequently, the ATPase, which was modeled as a reversible mechanism in (17), now actually catalyses the backreaction from ATP to ADP. This leads to a ATP concentration close to zero in the mitochondria (cf. Fig. 5a)). In contrast, the cytosolic ATP levels only drop to approximately 15% of their equilibrium value due to the increased glycolytic flux. HIF-mediated increase of glycolytic reaction rates apparently leads to the accumulation of glycolytic metabolites (see Fig. 5b)) and eventually, the system reaches a new steady-state (for longer simulation times, not shown). These results are consistent with experiments from [32], where wildtype mouse embryo fibroblasts (MEFs) possessed lower cellular ATP levels than HIF1α −/− knock-out MEFs under hypoxic conditions.
Simulation results presented in Fig. 6 outline the advantage of HIF-driven adaptation to hypoxic conditions. The cytosolic ATP levels attained by the system after 24 hours of adaptation under different oxygen concentrations have been computed, together with the corresponding respiration rate v resp (see (16)) at that time. When HIF is not considered (blue curve), both, ATP levels (Fig. 6a)) and respiratory rate (Fig. 6b)) remain high even for very low oxygen concentrations and only drop once it is very close to zero. Such behavior is characteristic for HIF KO cells and associated with high ROS formation [21]. Inhibition of the oxydative decarboxylation in the simulations (orange curve) results in lowered ATP levels already for moderately low oxygen levels (Fig. 6a)), but comes with the advantage of decreased respiratory chain activity (Fig. 6)). Additional amplification of the glycolytic rates (yellow and purple curves) leads to comparably low respiratory chain activity at least for low oxygen concentrations and helps to maintain higher ATP levels in the cytosol.

Sensitivity Analysis Under Normoxic And Hypoxic Conditions
We use the direct differential method to gain insights into the local sensitivity of the presented model with respect to 34 selected system parameters [26,33]. In general, the system of sensitivity equations which is associated with some parameter of interest p, reads: where we define The solution of the sensitivity equations with respect to k 1 in relative terms under normoxic conditions is shown in Fig. 7. There, the original system defined by (18)-(33) and parameters from Table 2 starts in the steady state given in Table 1 and remains there during the simulation while the sensitivity equations associated with parameter k 1 are solved simultaneously.
Analogous calculations are performed for 34 selected system parameters. To facilitate the evaluation and comparison of the results, we define a measure for the sensitivity of the system with respect to a certain parameter p as:  Fig. 7 Solution to the Sensitivity Equations for Parameter k 1 . The sensitivities of the metabolite concentrations with respect to the glucose uptake parameter k 1 are solved numerically. Using the same numeration as in Section 2.5, the left panel shows the relative sensitivities for metabolites 1-8 and the right panel shows the relative sensitivities for metabolites [9][10][11][12][13][14][15][16] where · 2 is the euclidean norm. By t end we denote the endpoint of the simulated time interval, which is t end = 6 · 10 3 s in the case of the sensitivity analysis under normoxic conditions. The outcome is compared to sensitivity results under hypoxic conditions in Fig. 7. The latter are obtained as follows: Again the system starts with the initial values from Table 1 and with parameters from Table 2, but this time the oxygen concentration is set to [O 2 ] = 0.05μM, as already in the simulations for Fig. 5, and the simulated time interval is extended to t end = 24h. The sensitivity equations are solved simultaneously for all selected parameters and respective sensitivity measures are computed and depicted in Fig. 7. These findings suggest that the system is dominated by the rate constant of the glucose uptake k 1 and of the ATP consumption k 6 during normoxia. The high influence of these two parameters is slightly dampened under hypoxic conditions and especially the HIF accumulation rate k 15 as well as the HIF-dependent amplification factor of the glucose uptake P 1 become more important.
Finally, the dependence of the steady-state concentration of cytosolic ATP on the most sensitive parameters from Fig. 8, k 1 , k 6 , k 15 , and P 1 , is investigated during normoxia and hypoxia. In Fig. 9, the cytosolic ATP concentration after 24h of simulated time is recorded, while varying the respective parameter over the interval from 50% to 200% of its baseline value from Table 2.
As expected, the ATP concentration does not change with the HIF-associated parameters k 15 and P 1 during normoxia (cf. Fig. 9a) and therefore, the respective curves are constant and coincide. However, the ATP concentration increases with decreasing ATP consumption rate constant k 6 . Diminished glucose uptake, represented by the interval where k 1 is below 100%, results in decreased ATP content, whereas increasing k 1 beyond 100% has only a small effect, since the ATP concentration in the cytosol is bounded by the total concentration of adenosine, A tot cyt . In the case of simulated hypoxia (Fig. 9b), the trend of the ATP curves with respect to variations of parameters k 1 and k 6 is similar to the normoxic case. In contrast with Fig. 9a, the ATP curve maintains an almost constant slope and does not come close to its limit when increasing k 1 beyond 100%, since the ATP level is generally lower in hypoxia. Not surprisingly, the ATP curve which is obtained by varying P 1 , is almost identical to the k 1 -dependent black curve. During hypoxia, the small constant c in the rate equation for the glucose uptake v * 1 can be neglected and hence the term scales linearly with both, k 1 and P 1 . Apparently, the ATP concentration in the cytosol is not at all sensitive to the imposed changes of the HIF-accumulation rate k 15 during hypoxia. In fact, the high sensitivity measure for k 15 in Fig. 8 appears, since the HIF concentration itself is very sensitive to k 15 and can reach very high levels for slightly increased k 15 , leading to a allegedly high sensitivity of the model to k 15 . But due to the saturation structure of the HIF-dependent activation/inactivation terms, a higher HIF concentration does not necessarily lead to a stronger response of the rest of the model.

Conclusions
In this study, we formulated the mathematical model of HIF-1 mediated regulation of cellular energy metabolism describing the reprogramming of cell metabolic processes from oxidative phosphorylation to glycolysis under reduced oxygen levels. The model considers the dynamics of sixteen biochemical species and the proton concentration with the underlying reaction processes localized in three intracellular compartments, i.e. the cytoplasm, mitochondrion and nucleus. More than sixty parameters of the model were calibrated using both the published data and the system steady-state based identification procedure. The model was validated by generating predictions which could be compared to empirical observations. The model behavior representing the cell metabolism switching over in response to transitioning from a normoxic to hypoxic environment is consistent with the current perception of the role of HIF-1 in hypoxia as detailed below: -Firstly, HIF-1 is known to prevent the persistence of potentially harmful ROS levels through the induction of pyruvate dehydrogenase kinase 1 (PDK1) [21,22]. -Secondly, glycolysis generates 2 molecules of ATP from one molecule of glucose whereas oxidative phosphorylation produces 36 ATP molecules per glucose molecule [17]. This difference in efficiency suggests that low oxygen tension will severely reduce cell energy availability [23].
The calibrated model presents a first step, i.e., a basic module for oxidative phosphorylation and glucolysis regulation, in the development an integrative multiscale model of cellular energetic metabolism. Such a model should consider the whole spectrum of metabolic fuels including the fatty acid metabolism, glutaminolysis. As glutamine is also required for the biosynthesis of proteins, it will be critical to identify an optimal physiological balance between its use for energy production and as a building block of proteins required for immune system function. Importantly, the immune cells metabolism is also regulated by multiple signaling networks, initiated by the formation of receptor synapses, determining the decisions on the cell fate under normal, hypoxic and inflammatory conditions. These also need to be considered on the way towards a comprehensive model of the cell.
There is growing need in developing pharmacological strategies for controlling specific metabolic processes that are skewed in many pathologies, such as autoimmune diseases, tumors, sepsis. Overall, we hope that our modeling approach should pave the way towards a systems level predictive model to be used for targeting the specific processes in the pathologies related to a skewed cellular metabolism [9].