Modeling dioxygenase enzyme kinetics in familial paraganglioma

Maintaining oxygen homeostasis is vital to many cellular processes. When a cell is deprived of oxygen, ATP production drastically decreases from 36 to 2 ATP molecules per molecule of glucose metabolized. To sense and adapt to a low oxygen environment, mammalian cells rely on two key mediators, HIF1α and HIF2α (hypoxic-inducible factors; referred to as HIFα) (Pugh and Ratcliffe, 2003; Wenger, 2002). HIFs are transcription factors with an oxygen sensitive alpha-subunit (HIFα) and a constitutively expressed beta-subunit (HIFβ). In normoxia, HIFα is constitutively hydroxylated by prolyl 4-hydroxylase domain (PHD) protein, which mediates HIFα interaction with the von Hippel-Lindau (VHL) E3 ubiquitin ligase complex that marks the transcription factor for rapid proteasomal degradation. During hypoxia, the enzymatic rate of PHD decreases allowing HIFα to remain stable and translocate to the nucleus, where it can dimerize with HIFβ to form transcriptionally active HIF complex. The HIF complex activates genes involved in glucose uptake, energy metabolism, angiogenesis, erythropoiesis, cell proliferation and apoptosis (Pugh and Ratcliffe, 2003; Wenger, 2002).

Mammalian cells express three PHD isoforms: PHD1, PHD2, and PHD3 (Berra et al., 2003). PHD1 is found exclusively in the cytosol. PHD2 is located in both the cytosol and nucleus, whereas PHD3 is only in the nucleus. PHD2 is believed to be the main oxygen sensor determining HIFα levels in a number of cell types. In fact, intracellular levels of PHD1 and PHD3 are almost undetectable in normoxia (Appelhoff et al., 2004). PHD enzymes operate in an ordered tri-tri reaction mechanism, requiring iron as a cofactor, α-ketoglutarate (αKG), and oxygen to hydroxylate HIFα, and producing carbon dioxide and succinate as by-products (Fig. 1) (Myllyla et al., 1977). In this mechanism, a 2-histidine-1-carboxylate facial triad (Hegg and Que, 1997) anchors the iron atom in the active site and maintains three coordination sites available for the substrates (Costas et al., 2004), which is consistent with the crystal structure of PHD2 (McDonough et al., 2006). Under conditions of reduced iron, αKG, or oxygen, the rate of PHD catalysis is severely reduced. PHD enzymes are also subject to competitive inhibition by Krebs cycle metabolites such as fumarate (Isaacs et al., 2005; Sudarshan et al., 2009; Xiao et al., 2012) and succinate (Letouze et al., 2013; Selak et al., 2005; Xiao et al., 2012) and by 2-hydroxyglutarate (Zhao et al., 2009) because of the structural similarity to αKG.

To date, a number of diseases have been hypothesized to involve PHD inhibition and HIFα activation (Pollard et al., 2005; Selak et al., 2005; Zhao et al., 2009). Here we focus on loss of succinate dehydrogenase (SDH) as it relates to the development of familial paraganglioma (PGL). In PGL, SDH loss leads to succinate accumulation in the mitochondria. Succinate can diffuse into the cytoplasm and competitively inhibit PHD resulting in the stabilization of HIFα (Selak et al., 2005).

While several mathematical models have been developed since 2004 to describe HIF signaling (Cavadas et al., 2013), most of the previous models are designed around the HIF pathway (Dayan et al., 2009; Kooner et al., 2006; Nguyen et al., 2013; Schmierer et al., 2010; Yucel and Kurnaz, 2007) – including at a minimum HIFα, PHD, VHL, and HIF-responsive elements (HRE) but also incorporating a combination of HIFβ, factor inhibiting HIF (FIH), vascular endothelial growth factor (VEGF), arrest-defective protein (ARD), and the transactivation domains of HIFα (N-TAD and C-TAD) – or the entire hypoxia response control network (Kohn et al., 2004; Yucel and Kurnaz, 2007). At a pathway or network level these models ignore the roles of α-ketoglutarate and succinate. One previous model did explicitly incorporate all of the PHD reactants (Qutub and Popel, 2006) but did not include the product succinate until its subsequent expansion (Qutub and Popel, 2007). In this expanded model, HIFα was completely stabilized at 21% oxygen independent of the presence of accumulated succinate. This prediction is not consistent with our observations in SDH deficient cells (Her et al., 2015). We recently developed a human embryonic kidney HEK293 lentiviral shRNA SDHB knockdown cell culture model where no detectable HIFα stabilization was observed upon succinate accumulation in 21% oxygen (Fig. 2A). Only when oxygen was lowered to 10% or 2% was HIFα stabilized. We also developed an immortalized mouse embryonic fibroblast Sdhc knockout model, which shows minimal HIFα stabilization at 21% oxygen and overall decreased HIFα responsiveness (Fig. 2B). From these cell culture models, we deduced that succinate accumulation and hypoxia synergistically inhibit PHD and concluded that an improved mathematical model was needed.

Can in silico tools be generated that accurately predict HIF levels at both physiologically-relevant as well as atmospheric oxygen conditions? Can this model capture the effects of succinate accumulation on PHD activity? Can this initial framework be expanded to more fully describe the HIF network, i.e. additional isoforms, greater crosstalk, and more compartments? And once developed, will this model make meaningful predictions, such as αKG and hyperoxia levels required for therapeutic treatment of PGL? These questions are the driving force for this work. Using our in vitro data as motivation and founded on PHD enzyme kinetics we developed an ordinary differential equations-based model to provide a quantitative framework for understanding HIF1α hydroxylation mediated by PHD2 and a plausible mechanistic explanation for the synergistic effect of succinate and hypoxia. We show that this model predicts that succinate and hypoxia synergistically inhibit PHD, and this inhibition can be overcome by the addition of αKG. Simulations from our model suggest new testable hypotheses, which can guide future experiments in a rational way, and this work further reinforces the notion that αKG and/or hyperoxia could be therapeutic in PGL.

The challenges in modeling the HIF pathway are as follows: there are three PHD isoforms and each may differentially contribute to HIF regulation. There exists a compartmentalization of these proteins and other factors in the HIF response, and specific localizations within the cell can affect rates and activities. Additionally, three isoforms of HIFα exist and each can regulate a specific set of genes. Thus, mathematical models face enormous obstacles in trying to account for all of the possible signaling crosstalk and to account for the abundance (and potential modified state) of each of the players in all of their potential subcellular locations. Moreover, there exists feedback regulation within the HIF pathway as well as crosstalk and synergism with other pathways, e.g. NFκB and mTOR (Cavadas et al., 2013). Consequently, several assumptions must be imposed due to a lack of experimental data and to simplify the modeling process.

Our model incorporates key molecular interactions in the HIF1α hydroxylation process mediated by PHD2. The molecular components and steps of the hydroxylation reaction are presented in Fig. 1. Detailed discussion of the model including reactions, reaction rates, and parameters are described in Materials and Methods. Under normoxic conditions, HIF1α protein is produced at a steady rate, but the protein is degraded either by non-specific protein turnover or from oxygen-sensitive hydroxylation and the resulting von Hippel-Lindau (VHL)-mediated proteasomal degradation. In this initial stage of model development, we omit any molecular details of the VHL complex and we assume that once hydroxylated, HIF1α is committed towards a proteasomal destruction fate. This is accomplished by setting to zero the rate constant for the one term influenced by hydroxylated HIF1α (see Materials and Methods). Future iterations of the model will incorporate the VHL pathway in greater detail.

Although the PHD isoforms have specific cytoplasmic and/or nuclear localizations, for simplicity in this initial development of the model we have considered only a single PHD entity – PHD2. Furthermore, since the levels of the various reactants and products (and their fluxes) are not specifically known in the mitochondria, the nucleus, and the cytoplasm, we assume a single cellular compartment as opposed to a multi-compartment model. Conceptually, this is equivalent to a well-mixed cytosolic compartment where metabolites like αKG and succinate can freely move in and out of the mitochondria and PHD2 can freely translocate in and out of the nucleus. In such a model, changes to these global levels can reflect different subcellular or environmental contexts. Finally, we have focused our attention on only the HIF1α isoform at present, leaving the incorporation of other isoforms for future iterations.

Our in vitro experiments (Her et al., 2015) showed that HIF1α is stabilized at 2% oxygen in control cells and at both 10% and 2% oxygen upon succinate accumulation due to loss of SDH function (Fig. 2A,B). Stabilization of HIF1α was not observed at 21% oxygen with or without succinate accumulation. We first applied our mathematical model to simulate the synergistic inhibition of PHD by succinate and hypoxia, using the intracellular concentrations of αKG and succinate measured in our SDH-loss cell culture experiments (Table 1). Succinate is a competitive inhibitor of PHD2 in addition to being a product of the reaction (Myllyla et al., 1977). Fig. 3 shows HIF1α concentration as a function of time using measured levels of αKG and succinate determined in untreated cells (Table 1) and selected oxygen amounts ranging from 0.01% (∼anoxia) to 21% (atmospheric). The model predicts that under hypoxic conditions, HIF1α protein is rapidly stabilized, reaching steady-state levels within ∼12 h (the HIF1α concentration achieved being inversely related to oxygen amount). At steady state, PHD2 activity balances the level of HIF1α being synthesized with the level of HIF1α being degraded. Fig. 4 shows the dependence of HIF1α concentration on the combined effects of PHD inhibition by succinate and hypoxia oxygen levels simulated at 48 h post exposure to hypoxia. This time point mimics our cell culture experiments (Fig. 2A,B). At 21% oxygen, both non-specific target (scr.) and shRNA SDHB target (shRNA) lentiviral-treated HEK293 cells did not stabilize HIF1α (Figs 2A and 4A). At 10% oxygen, we observed a dramatic increase in HIF1α levels (shRNA versus scr.), and at 2% oxygen HIF1α is stabilized under both conditions. Our mathematical model qualitatively confirms this observation (Fig. 4A). Similarly in the knockout model, Sdhc fl/− cre− control cells only exhibit appreciable HIF1α stabilization at 2% oxygen, whereas Shhc fl/− cre+ cells have stabilized HIF1α at both 10% and 2% oxygen (Fig. 4B). For comparison, data from HeLa cells (Jiang et al., 1996) are plotted along with the results in Fig. 4B. While cell-type specific differences undoubtedly exist, it is remarkable how well the model agrees with the experimental data from another cell type. However, the HeLa cell data suggest an even stronger dependence on oxygen levels, showing a majority of HIF1α degraded already at 2% oxygen.

The dependence of HIF1α concentration on oxygen levels is more fully appreciated in log-log plots (Fig. 4C,D). There is marked sensitivity in the range of 2% to 10% oxygen for control cells, and this shifts to 10% to 21% upon succinate accumulation in knockdown or knockout cells. HIF1α appears maximally stabilized at 0.2% oxygen (bold arrows, Fig. 4C,D), consistent with experimental data from HeLa cells where knockdown of both PHD2 and PHD3 in these cells did not lead to much further HIF1α stabilization (Fig. 5B of Stiehl et al., 2006). The model can be utilized to make predictions about the therapeutic rescue of hyperoxia. Consider 5% oxygen in control HEK293 cells, where the concentration of HIF1α is 0.01 nM (Fig. 4C). In order to achieve this concentration of HIF1α in SDH-loss cells, the model predicts that the oxygen level would need to be increased to 16% (gray lines, Fig. 4C). Thus, our enzyme kinetic model confirms our previous experimental result (Her et al., 2015) that PHD is synergistically inhibited by succinate accumulation and hypoxia in cell culture SDH-loss models of PGL.

As originally shown by others (MacKenzie et al., 2007; Smith et al., 2007), we also showed in our SDH-loss models of PGL that addition of cell-permeable αKG could overcome succinate inhibition of PHD (Her et al., 2015). This is rationalized by the fact that succinate and αKG compete for the active site of PHD2 (Letouze et al., 2013; MacKenzie et al., 2007; Xiao et al., 2012). In our SDH-loss cell culture experiments, we observed that addition of 250 μM cell-permeable octyl-αKG reversed the synergistic inhibition of PHD by succinate and hypoxia (10% oxygen) in both knockdown and knockout cells (Fig. 2C,D). To model this effect, we mimicked the cell culture experiments and simulated a 48-h incubation at 10% oxygen before adding 250 μM αKG for an additional 12 h. The results demonstrate that addition of 250 μM αKG completely reverses the combined succinate and hypoxia inhibition of PHD2 to produce HIF1α concentrations observed in the control cells (dot-dash line, Fig. 4C). In fact, the model prediction for the knockout cells is an over-rescue at this αKG concentration (dot-dash line, Fig. 4D). This result further validates our mathematical model, showing that it qualitatively confirms the results of our in vitro experiments showing αKG reversal of PHD inhibition by succinate and hypoxia.

After shifting to hypoxic conditions, HIF1α is stabilized and initiates hypoxic signaling. A negative feedback response partially decreases HIF1α concentration through an increase in PHD2. It has been proposed that the feedback control of HIF1α creates a new set point in the system, and allows room for a future increase in HIF1α stabilization in response to a further hypoxic shift (Stiehl et al., 2006). HEK293 cells shifted to 1% oxygen for 72 h responded to a shift to 0.2% oxygen with increased stabilization of HIF1α (Stiehl et al., 2006). This response was also followed by a resetting of stabilized HIF1α levels over the course of 4 h, presumably in preparation for a subsequent insult. Our model nicely illustrates this concept (Fig. 6A). Without feedback regulation (gray curves, Fig. 6A), HIF1α concentrations approach the theoretical maximum. Since oxygen levels are already low, the production of additional enzyme is necessary so that the system is once again poised to respond to any further decrease in oxygen. This results in a resetting of HIF1α levels below a certain oxygen threshold (black curves, Fig. 6A). The observed feedback regulation has important therapeutic implications. When we model our cell-permeable αKG experiment in the context of a feedback system, we see that the model predicts an over-rescue of the cells with 250 μM αKG (Fig. 6B), reminiscent to that predicted in the knockout cells (Fig. 4D). This result suggests that negative feedback regulation of HIF1α will make cells more sensitive to αKG therapy.

The experimental values for intracellular succinate and αKG concentrations shown in Table 1 demonstrate that succinate accumulates in SDH-loss cell culture models, as expected. Interestingly, there is also a concomitant decrease in αKG concentration. These reciprocal changes exacerbate the observed ratio of αKG to succinate. To explore this idea, we used our expanded mathematical model with feedback regulation to simulate the effect of addition of cell-permeable succinate to control HEK293 cells. Experimental data at 10% and 21% oxygen with addition of 20 mM dimethyl succinate show increased HIF1α stabilization only at 10% oxygen (supplementary material Fig. S4C of Her et al., 2015), though it is difficult to determine how much of the added compound enters cells and is de-esterified. Again our simulations mimic the experimental conditions with a 48-h incubation in hypoxia before adding 20 mM succinate for an additional 12 h. The model predicts a dramatic stabilization of HIF1α at 10% oxygen but not at 21% oxygen, even with this high concentration of added succinate (Fig. 6C), consistent with the experimental results. The model suggests that if this high level of succinate accumulation were achieved, more αKG would be needed to counter the effect. Thus, the balance of αKG to succinate is a key determinant of HIF1α levels.

Lastly, we used our model to predict that the ability of αKG to overcome the effects of succinate accumulation in HEK293 cells. The model predicts that small changes in αKG concentration (relative to succinate), should have dramatic effects on HIF1α concentration (colored curves in Fig. 6D). Addition of 10 μM αKG restores HIF1α concentration to the level of control cells at both 2% and 10% oxygen (blue dotted line, Fig. 6D). The addition of 5 μM αKG is likewise predicted to have large effects (red dot-dash line, Fig. 6D), and even addition of 1 μM αKG may have therapeutic effect when combined with hyperoxia. Future studies will subject these predictions to experimental validation.

Motivated by our in vitro experimental studies, we here utilize first principles of enzyme kinetics to generate a mathematical model (with or without feedback regulation) that can rationalize the observed synergistic inhibition of PHD by succinate and hypoxia in an SDH-loss cell culture model of PGL. Using experimental values for intracellular succinate and αKG concentrations (Her et al., 2015), we showed that this model confirms the observed synergistic inhibition of PHD by succinate and hypoxia and also showed how this inhibition can be overcome by addition of αKG. Of particular interest is the predicted range of PHD inhibition by succinate and hypoxia. As shown with our model, SDH-loss cells cultured in 21% oxygen do not stabilize HIF1α because PHD quickly hydroxylates HIF1α and the protein is rapidly degraded (Fig. 4). However, in the physiological oxygen range, our model predicts that HIF1α is stabilized by the combined inhibition of succinate and hypoxia (Fig. 4). The model results are consistent with experimental data and demonstrate how atmospheric oxygen level masks important biology that can be observed under physiologically-relevant conditions.

The next step in the development of this mathematical model will be to predict other experimental results that can then be validated by related experimental systems in our laboratory. It is likely that the ratio of PHD2 to HIF in different tissues and cell types will likely modulate the extent of HIF accumulation, and this avenue will be explored in the future. Additionally, simulations of acute hypoxia versus chronic will be used to refine the model moving forward. Relevant to the SDH-loss model of PGL, predictions of αKG and hyperoxia levels required for therapeutic treatment of PGL are an important avenue of future study. Our model can also predict effects of succinate accumulation and hypoxia on HIF2α levels. Additionally, other PHD isoforms can be readily incorporated into the model. More experiments will need to be performed to determine if the model confirms western blot data for HIF2α and HIF3α as well as PHD1 and PHD3 obtained in such cell culture models. While the model currently uses values from shRNA SDHB knockdown cells and Sdhc knockout cells (Her et al., 2015), future simulations with parameters more relevant to neuroendocrine cells will be performed. Additionally, iron will be explicitly incorporated into the model, allowing for simulations of competitive inhibitors such as cobalt or iron chelators such as deferoxamine. In summary, this mathematical model framework for the PHD dioxygenases can be applied to gain insight into effects of iron, succinate, hypoxia, and αKG on PHD and HIFα, guiding future experimental designs.

The system of nonlinear differential equations was solved with R software (Version 2.13.0) using the lsoda ODE solver in the package deSolve (Soetaert et al., 2010). Initial concentrations of all species are given in Table 1. From the rate equations, there are 15 unknown rate parameters: k_on,αKG, k_off,αKG, k_on,O2, k_off,O2, k_on,HIF1α, k_off,HIF1α, k_cat,HIF1α, k_on,HIF1α-OH, k_{off,HIF1α-OH}, k_on,CO2, k_off,CO2, k_on,Su, k_off,Su, k_{production,HIF1α}, and k_{turnover,HIF1α}. Estimates of these parameters were determined by nonlinear least squares fitting of the model to an estimation of the experimental data in Fig. 2A (see below) using the simplex and inductive search hybrid (SIH) algorithm (Offord and Bajzer, 2001). Qualitatively for the scr. control: at 21% and 10% oxygen levels HIF1α is essentially unchanged from its starting value but sharply rises at 2% oxygen approaching the theoretical maximum value q_max. For the shRNA case: HIF1α is essentially unchanged from its starting value only at 21% oxygen and is substantially larger at 10% and nearly maximal at 2% oxygen. In practice, this was accomplished by least squares fitting to the vector of values (0,0,0.5,0,0.25,0.85)×q_max. The fitted values of the rate constants are reported in Table 2.

Next, we reasoned from the qualitative data in Fig. 2 that the maximal level of HIF1α attained is decreased in the knockout cells as compared to the knockdown cells and that the basal level of HIF1α at 21% oxygen is slightly increased. Not having quantitative information about these differences, we arbitrarily adjusted the three model parameters most directly responsible for these effects by a scale factor: increased k_{production,HIF1α} by factor of 2.5, increased k_{turnover,HIF1α} by a factor of 6.25, and decreased k_on,HIF1α by a factor of 100 (Table 2).

The last step in model validation was to confirm that the parameter estimates are both plausible and reasonable and to determine if the kinetics can effectively be simplified in some instances. For example, the rate of CO₂ re-entering the complex after catalysis is negligible (2.78×10⁻⁵ µM⁻¹ min⁻¹) as expected. This indeed does lead to a high K_d,CO2 value (3.5 M), which is consistent with crude estimates of the K_i for CO₂ on the order of at least mM (Myllyla et al., 1977). Additionally, while the value of K_d,HIF1α is low due to the small off rate, the value of K_m,HIF1α (3 µM) is consistent with reported values (Hirsila et al., 2003). Finally, it has been reported that the upper limit of a diffusion-limited rate could reach 10¹⁰ M⁻¹ s⁻¹ (or 6×10⁵ µM⁻¹ min⁻¹), so that k_on,HIF1α in the knockdown cells while large (3310 min⁻¹ μM⁻¹) has not reached this limit.

We thank Z. Bazjer for valuable discussions.