A computational model of intracellular oxygen sensing by hypoxia-inducible factor HIF1α

The transcription factor HIF (hypoxia-inducible factor) plays a crucial role in mammalian response to oxygen (O₂) levels. HIF1, the first characterized member of the HIF family, transcriptionally activates hundreds of genes associated with angiogenesis in cancer, exercise and ischemia; energy metabolism; nutrient transport; and cell migration (Semenza, 2004; Wang et al., 1995).

Angiogenesis, the formation of blood vessels from preexisting vessels, is regulated in part by local tissue O₂ levels. This regulatory pathway links cell and tissue metabolic demand with vascular oxygen supply. The pathway is intricately governed by HIF1 regulation and HIF1 transcriptional activation of angiogenic factors.

HIF1 is a heterodimer, comprised of subunits HIF1α and HIF1β. The beta subunit is constitutively expressed in cells. Expression of the alpha subunit may be induced by a number of pathways, and its degradation is highly sensitive to O₂ levels. Called a `master switch for hypoxic gene expression' (Powell, 2003; Semenza, 2004), intracellular HIF1α in normoxia is experimentally undetectable; during hypoxia, it rapidly accumulates in the cell nucleus, and triggers gene expression. Molecular players involved in this process have come to light over the past six years. In normoxia, enzymes called prolyl hydroxylase domains (PHDs) react with HIF1α (Fig. 1). PHDs hydroxylate HIF1α at Pro402 and Pro564 in the oxygen-dependant degradation domain. The activity of PHDs depends on the amount of oxygen available. Three isoforms of HIF PHDs exist: PHD1, PHD2 and PHD3. Each isoform performs a separate function, with different kinetic properties and primary cellular locations (Appelhoff et al., 2004). PHD2 is the most abundant prolyl hydroxylase isoform in the cell cytoplasm during normoxia, and it has been credited as a controller of steady-state HIF1α concentrations under these conditions in a range of cell types (Berra et al., 2003). Conversely, in normoxia, intracellular PHD1 and PHD3 levels are low, if experimentally detectable (Appelhoff et al., 2004).

Following the reaction with PHD, the hydroxylated HIF1α is free to bind to a von Hippel Lindau (VHL) ubiquitin ligase complex, which tags HIF1α for proteasomal destruction. Like PHD2, VHL is found primarily in the cell cytosol in many cell types (Los et al., 1996), though VHL trafficking between nucleus and cytoplasm may be necessary in HIF1α degradation during reoxygenation from hypoxia (Groulx and Lee, 2002). VHL forms a stable complex with proteins including Elongin B, Elongin C, Cul2 and Rbx1 (Kamura et al., 1999). This VHL ligase complex binds to hydroxylated HIF1α and tags the protein with a polyubiquitin tail (Ivan et al., 2001; Kamura et al., 2000). A multi-protein complex, called the proteasome, recognizes this tail and destroys HIF1α.

In hypoxia, HIF1α escapes hydroxylation, accumulates and enters the cell nucleus, where it binds to HIF1β (known also as ARNT) (Fig. 1). The dimer transcriptionally activates a host of genes, including those encoding the angiogenic protein vascular endothelial growth factor (VEGF) and its receptor, Flk-1 or VEGFR2 (Milkiewicz et al., 2003); platelet-derived growth factor (Bos et al., 2005); and erythropoetin (Marti, 2004). By activating these genes, HIF1 contributes to angiogenesis, which provides nutrients to facilitate tumor growth or to extend muscle contraction, for example. As prominent players in the cell response to hypoxia and the onset of angiogenesis, HIF1 and its related pathways are attractive therapeutic targets in cancer and ischemia (Hewitson and Schofield, 2004).

In vitro studies have shown how the hypoxic response varies in tumors, based on their vascular microenvironment (Blouw et al., 2003). A balance of HIF1α levels and HIF1α activity seems necessary to achieve health (Josko and Mazurek, 2004; Koshiji and Huang, 2004). The underlying molecular mechanisms of how this balance is achieved and how the system responds to its microenvironment are not fully understood. One hypothesis is that HIF1α acts literally as a `switch' - an on/off mechanism for the onset of hypoxia-induced angiogenesis when a critical O₂ level is reached (Kohn et al., 2004). Fundamental issues include understanding how HIF1α acting as a generic switch would be correlated to the varied hypoxic responses found in tumor cells. Alternatively, if HIF1α and its pathways do not act as a switch, the observed sensitivity to oxygen and the rapid induction of hypoxic genes would need to be otherwise explained. We address these questions, by developing a detailed model of HIF1α degradation, which allows molecular mechanisms to be tested quantitatively. The one known existing computer simulation related to HIF1, a network representation of cell hypoxic response, has included a subset of core HIF1 pathways and tested the hypothetical dependency of gene expression on HIF1α synthesis and degradation rates (Kohn et al., 2004). This model led to the above HIF1α switch hypothesis, which to our knowledge, has not been tested further. No computational model has explored the biochemical kinetics of the HIF1 pathways in detail or in a quantitative relationship to experimental data. The current model is the first molecular level, mechanistic model of HIF1 hydroxylation and degradation. We used the model to study the effects of different intracellular molecular compositions on hypoxic response, where the cellular microenvironment is currently inaccessible in vivo, and is only measurable in isolated instances in vitro.

From the model, we predict several key characteristics of the mechanisms involved in the HIF1 pathway. We show that HIF1-activated cellular responses can be divided into two categories depending on the molecular environment: a switch-like response to O₂ levels, and a gradual one. We found conditions where iron and PHD2 are individual sensors of oxygen and determinants of the hypoxic response; and we showed the combined effects of three highly oxygen-sensitive compounds. From these studies, we compare two proposed therapeutic strategies targeting the HIF1 pathway, iron supplementation and PHD2 targeting, and predict under what microenvironments either would most efficiently increase HIF1α hydroxylation. These observations contribute to a better understanding of the hypoxic response at the molecular level and should stimulate further computational and experimental exploration, with particular applications to therapies that target cofactors in HIF1α hydroxylation.

Double reciprocal plots of the hydroxylation reaction before binding to HIF1α were consistent with experimental results from the collagen PHDs (Myllyla et al., 1977); this confirmed that the model represents the uncompetitive binding of Fe²⁺, O₂ and 2-oxoglutarate (2-OG) to PHD2 (Fig. 2A shows the example of iron binding). The reaction can proceed without ascorbate (Fig. 2B,C), however, high ascorbate concentrations (>100 μM) significantly increase the reaction rate (Fig. 2B).

The model was compared with independent experimental data from several sources to validate the oxygen dependency and the time course of HIF1α hydroxylation by PHD2 (Fig. 3). One form of validation was the amount of VHL captured at different O₂ levels. Tuckerman et al. measured the activity of endogenous HIF-PHDs (from MDA-MB-435 cell extracts) by a pVHL capture assay (Tuckerman et al., 2004). Oxygen concentrations in the simulations, all ⩽200 μM, are below the K_m values of PHD2 for O₂ of 250 μM (Table 2). This is reflected in Fig. 3A, where the model is compared with published data (Tuckerman et al., 2004). The 21% oxygen values show a near-linear increase in the fraction of the maximal pVHL capture with time; it is reasonable to conclude (Tuckerman et al., 2004), that all of the experiments were conducted below the threshold for oxygen saturation.

Model calculations for HIF1α half-life in normoxia fall within the range of experimental data (Fig. 3B). Experimental half-life of HIF1α was estimated as 5-8 minutes from reoxygenated (1-2% O₂ to 20% O₂) CCL39 and NHE-1 cell lysates (Berra et al., 2001). [Using oxygen solubility in water at 37°C = 1.30 μmol/l/mmHg (Tuckerman et al., 2004), exposure to 1% oxygen corresponds to an O₂ level of 9.9 μM; 21% to an O₂ level of 207 μM.] This complements independent finding in lysates from Hep3B cells exposed to 21% oxygen after hypoxia (1% oxygen) where the half-life of HIF1α was found to be <5 minutes (Huang et al., 1998); and from HeLaS3 cells reoxygenated at 20% from hypoxia (0.5% oxygen), where the half-life of HIF1α was 8 minutes (Jewell et al., 2001).

A third form of validation was the comparison of the predicted relative HIF1α accumulation of the model at different oxygen levels, with data from HeLa cell nuclear extracts (Jiang et al., 1996) (Fig. 3C). Studies have shown that HIF1 expression is maximal at low oxygen concentrations in vivo [e.g. at 0 or 1% O₂, in normal ferret lung ventilated for 4 hours (Yu et al., 1998)], and in vitro assays indicated the most pronounced changes in HIF1 expression occur at O₂ levels considered physiologically relevant (0-5%) (Jiang et al., 1996). To compare the computational model with these experiments, a constant [HIF1α]₀ was assumed, and [HIF1α] that was not hydroxylated was taken as a relative measure of HIF1α nuclear accumulation.

Fourth, the effects of ascorbate and iron in the model are qualitatively comparable to HIF1α expression observed in human prostate adenocarcinoma (PC3) cells (Knowles et al., 2003). Fig. 3E shows the model results of supplementing the system with 2000 μM ascorbate. Relative HIF1α values are a fraction of the maximum HIF1α expression during hypoxia, without supplementation. For different cell types, including PC3 cells, the PHD2:HIF1α concentration ratio has not been quantified. Fig. 3E provides an example of how the ratio affects the accumulation of HIF1α, during the first 3 hours of exposure to different O₂ levels. A PHD2:HIF1α ratio of 0.004 is an estimate from measurements in breast carcinoma cells (Tuckerman et al., 2004). In normoxia, ascorbate and iron supplementation have similar effects on suppressing HIF1α expression (Fig. 3F). Comparable experimental results in PC3 cells show decreased expression of HIF1α over time with ascorbate supplemented at 2000 μM, and no appreciable expression of HIF1α in normoxia after 4 hours of supplementation with 25 μM ascorbate or >40 μM FeCl₂ (40 μM FeCl₂ added to medium containing ∼26 μM Fe²⁺), see figure 2C,D in Knowles et al. (Knowles et al., 2003).

Sensitivity analysis was performed to confirm estimates for the unknown kinetic rate constants. The parameter of interest was varied over a minimum range of 1000-fold, while the remaining parameters were held constant. Calculated HIF1α half-life values were compared with experimental data, and this was used to narrow the range of reasonable parameter values. Details of the analysis are in provided in the Materials and Methods. Estimated values for five kinetic binding rates were determined: k_cat,Hα=0.098-0.164 minute^-1; k_off,Fe2=36 minute^-1; k_off,DG=10.8 minute^-1; k_off,O2=10.8 minute^-1; k_off,AS=3.6 minute^-1; k_off,Hα=0.7 minute^-1. Binding of Fe²⁺, 2-OG and O₂ to PHD2 are the more reversible steps in the hydroxylation reaction, with significant off-rates relative to on-rate binding. The final step in the hydroxylation, binding to HIF1α is largely irreversible, as indicated by a low k_off,Hα value and a significant k_cat,Hα.

A recent study reported the apparent K_m for Fe²⁺ of 0.03 μM (Hirsila et al., 2005). Using this estimate, as opposed to the K_m of 2 μM, estimated from binding of factor-inhibiting HIF (FIH) with Fe²⁺ in the hydroxylation of HIF1α (Koivunen et al., 2004), the model predicts a higher specific activity for PHD2 than found in vitro (Hirsila et al., 2005) (see supplementary material Fig. S1). This discrepancy may reflect different concentrations of PHD2 and HIF1α relative to other compounds in the hydroxylation reaction; the model used initial conditions where quantitative concentrations of PHD2 were reported (Tuckerman et al., 2004).

Sensitivity analysis for all kinetic parameters found from experiments, was also performed using the protocol described for the estimated parameters (see Materials and Methods). Over a wide range of feasible K_m,Fe2, K_m,DG, K_m,O2, K_m,Asc and K_m,Hα values (supplementary material Figs S2, S3), the changes in oxygen sensitivity were consistent, and the model features described below can be considered robust over these values.

In many in vitro cell extract experiments monitoring HIF1α reactions, there is an excess concentration of initial 2-OG, iron, ascorbate and PHD2. When any of these compounds was limiting, the sensitivity to oxygen in the model was uniform at all O₂ levels from 0-200 μM. Fig. 4 shows how initial reactant concentrations affect HIF1α hydroxylation at different O₂ concentrations. When [PHD2]₀ is in excess, the response to decreasing O₂ results in a steep change in hydroxylation upon reaching hypoxia - whereas 21 and 10% O₂ levels result in similar amounts of HIF1α_hydroxylated; at hypoxic levels of 1%, the amount hydroxylated at 20 minutes is half of that at normoxia (Fig. 4A). However, when [PHD2]₀ is low, the response to a 20% drop in O₂ levels is far less sensitive (Fig. 4D). Sensitivity to oxygen is measured by the slope of the [HIF1α_hydroxylated] vs [O₂] curve. A constant slope represents a uniform sensitivity across O₂ levels. Iron has a similar effect on HIF1α at different O₂ concentrations. When iron is available in excess, the response curve is steep. At initial iron concentrations of 0.05 μM, or one-thousandth of the default value, the amount of HIF1α hydroxylated is linearly related to the O₂ level (Fig. 4C). Changes in the concentration of ascorbate did not have such a significant effect on hydroxylation (Fig. 4F). The steepness of the response decreased when two or more compounds were limiting factors in the HIF1α hydroxylation (Fig. 5). Table 2 shows the theoretical range of initial conditions where steep, switch-like changes in hydroxylation occur.

The assumed maximal HIF1α half-life range of 5-8 minutes in normoxia is consistent with at least three experiments (Berra et al., 2001; Huang et al., 1998; Jewell et al., 2001). However, the half-life of HIF1α upon reoxygenation depends on conditions such as duration of hypoxic exposure. For example, exposure to low oxygen levels beyond 6 hours, appears to decrease HIF1α half-life (G. Semenza, Johns Hopkins University, Baltimore, MD, personal communication). For short-term hypoxic exposure, the model assumes a maximum half-life of 5-8 minutes. The cited HIF1α half-life studies exposed cells to hypoxia for 1 hour (Jewell et al., 2001), 4 hours (Berra et al., 2001), and 4-6 hours (Huang et al., 1998) before reoxygenation. Beyond 4-6 hours of hypoxia, synthesis of HIF1α and PHD2 proteins occur. Fig. 6 shows the variability of HIF1α hydroxylation under conditions of chronic hypoxia, where a synthesis production term was added to the mass balance equations for HIF1α and PHD2 in the model. The simulated curves represent both the synthesis of HIF1α and its hydroxylation by increasing amounts of PHD2. In vivo, systems can adapt to chronic conditions, decreasing HIF1α expression within days of hypoxic exposure. A balance of HIF1α and PHD2 synthesis is a possible contributing mechanism.

Therapeutic strategies to enhance HIF1α hydroxylation Increased HIF1α nuclear expression has been associated with poor prognosis in several cancers (Nomura et al., 2004; Zagzag et al., 2000), and decreased susceptibility to radiotherapy (Vordermark and Brown, 2003; Vordermark et al., 2004). By enhancing HIF1 hydroxylation, the proteasome degradation rate of HIF1 can be increased, leading to a drop in nuclear accumulation of the dimer and associated hypoxia-dependent transcriptional activation. Targeting cofactors in the PHD reactions is one viable approach to increasing the hydroxylation rate. As our computational model predicted different oxygen sensitivity based on intracellular concentrations of key cofactors in the HIF1 reactions, we hypothesized that the microenvironment would also dictate the effectiveness of therapeutic strategies, and computationally, we would be able to predict the relative efficacy of each strategy at the molecular level.

Using the model, we tested two proposed approaches to increasing HIF1α hydroxylation: (1) increasing intracellular iron concentration (McCarty, 2003; Siddiq et al., 2005; Wartenberg et al., 2003) while simultaneously increasing intracellular ascorbate (Jones et al., 2006; Knowles et al., 2006; Knowles et al., 2003) and (2) increasing the expression of the main cytosolic prolyl hydroxylase, PHD2, directly (Fig. 7). Results show that increasing ascorbate is a proportionately more effective way to increase hydroxylation (Fig. 7A,B), when all other compounds are not limiting factors. At low levels of iron, doubling the amount of ascorbate from the standard in vitro level of 1000 μM, increases hydroxylation by as much as 60% at 50 μM O₂; when iron is above 2 μM, the effect of additional ascorbate diminishes to a constant 3% increase in hydroxylated HIF1α (Fig. 7C,D).

A rationale for developing the HIF model was to test several hypotheses on the signaling pathway, as it relates to hypoxic response, and to then use the results to evaluate therapeutic approaches targeting the onset of angiogenesis. An existing hypothesis of a switch-like change in HIF1α expression in response to a decrease in O₂ levels to a critical level was tested (Kohn et al., 2004). The model demonstrated that based on the molecular environment and characteristics of the cell type, the response to hypoxia varies considerably.

There are several ways that cells are hypothesized to sense and respond to oxygen. The sensor of O₂ levels has been independently attributed to iron (Postovit et al., 2005), HIF1α and HIF1α PHDs experimentally, but how their responses differ and when, has yet to be understood. The model shows two classes of HIF1 oxygen responses: a steep drop and a gradual drop in HIF1α hydroxylation in response to decreasing O₂ levels. When all hydroxylation reactants are in excess, a steep drop in hydroxylation occurs during hypoxia. When 2-OG, Fe²⁺ or PHD2 are limiting, model results show the gradual response - a near-linear relationship between HIF1α hydroxylation and O₂ level, i.e. a constant sensitivity to O₂; this reflects the saturation kinetics used to represent the binding of HIF1α to these compounds (Fig. 4C-E). HIF1α hydroxylation is reduced when two or more required compounds are limiting (Fig. 5); what is notable is not only a significant reduction in the net amount of HIF1α hydroxylated, but a significant decrease in the relative sensitivity to O₂ levels, defined by the steepness of the slope of the [HIF1α_hydroxylated] vs [O₂] curve.

Characterizing the role of ascorbate in the hydroxylation is complicated by its reaction with oxidized iron. At low levels of ascorbate (i.e. 1 μM, one-thousandth of the default initial value), HIF1α hydroxylation follows the curve shown in Fig. 4F. The decline in HIF1α hydroxylation at low O₂ levels is less steep than the case where all compounds are in excess, however, it is not approaching linearity. This reflects the dual role of ascorbate: reactivating PHD2, by reducing the accumulating Fe³⁺; and, binding independently to the saturable enzyme complex formed by the binding of PHD2 with Fe²⁺, 2-OG and O₂ (Majamaa et al., 1986). The decarboxylation of 2-OG without subsequent hydroxylation (termed uncoupled decarboxylation), which is catalyzed by PHDs and requires ascorbate, is not yet considered in the model. From the above observations, HIF1-activated cellular responses can largely be divided into two categories depending on the molecular environment: a steep, switch-like response to O₂ levels, and a gradual one. The model shows that independently Fe²⁺ and PHD2 can act as the determinant of which response occurs. The high sensitivity of HIF1α hydroxylation to these compounds, notably at low (<5 μM) levels of oxygen, suggests that both iron and PHD2 play the role of hypoxic sensor.

Cellular in vivo concentrations of these compounds are difficult, if not yet possible, to measure. The concentration of HIF1α in vivo is thought to be in the sub-nanomolar range (Tuckerman et al., 2004). Based on the relative concentrations of other reactants in the HIF1 system, the response to O₂ could be anticipated by the model. Although the quantitative values are largely unknown in vivo, it is known that iron and PHD2 concentrations are variable and affected by conditions such as hypoxia and anemia (Berra et al., 2003; Wardrop and Richardson, 1999). From the model, it would be predicted that this variability confers an advantage. With the HIF system being highly sensitive to the effects of multiple O₂-dependent compounds, the response to O₂ levels could be finely tuned, as well as robust.

Furthermore, the computational prediction that increasing ascorbate has significantly more effect at iron levels below 2 μM compared with higher iron levels, might later be extrapolated to characterizing tumor responsiveness to therapy based on in vivo microenvironments.

Using oxygen solubility in water, 37°C STP, 1.30 μmol/l/mmHg, [O₂] of 21% corresponds to an O₂ concentration of ∼200 μM for in vitro experiments using cell lysates (Tuckerman et al., 2004). This is below the K_m of oxygen reacting with PHD2, which was reported as 250 μM (Hirsila et al., 2003). Typical in vivo tissue concentrations are significantly lower, in the range of 6-25 μM, corresponding to 5-20 mmHg tissue pO₂. The high K_m indicates that the HIF1α reaction response remains highly sensitive to tissue O₂ levels all the way from 0% oxygen to normoxia, under certain conditions. The model supports research on HIF1α nuclear expression changes being greatest at oxygen levels below 5% (Jiang et al., 1996), while showing how the most pronounced response by HIF1α to acute hypoxia can be achieved only with sufficient Fe²⁺, 2-OG and PHD2 (Fig. 4A,B). Varying the relative concentrations of the compounds involved in HIF1 hydroxylation, Fe²⁺ and PHD2 in particular, alters this sensitivity.

Other forms of variability and specificity could add to the flexibility of the HIF1 system in responding to oxygen. Although, as mentioned, PHD2 is the main PHD enzyme present in the cell cytoplasm during normoxia, the concentration of each PHD isoform varies by cellular microenvironment (Appelhoff et al., 2004; Berra et al., 2003). Hypoxia and estrogen change the relative concentrations of the PHD isoforms, and their relative contribution to the hydroxylation of HIF1α (Appelhoff et al., 2004).

While discussing other factors that influence HIF1α hydroxylation, it is appropriate to mention the limitations of the presented model. Five independent kinetic rates used in the model are unknown, and their values were estimated computationally. Their estimated values appear consistent with experimental results, and the key features of the hydroxylation reaction (steep and gradual responses to hypoxia) are robust over a range of all kinetic parameters. Additional experiments would be useful to validate the binding on- and off-rates. Kinetic rate constants taken from in vitro experiments may depend on slight differences in pH or other experimental variables (e.g. relative concentration of proteins) that have a limited degree of controllability within the cell. The enzyme used experimentally for the binding reactions was a minimal HIF1α peptide (residues 556-574) (Hirsila et al., 2003). Although increasing peptide length does not seem to affect the binding affinities in vitro to pVHL (Hon et al., 2002) or PHD2 (Hirsila et al., 2003) with HIF1α, the possibility that it may in vivo, cannot be ruled out. Additionally, estimates for binding kinetic rates are assumed sequential - this may or may not be what is actually measured by experiments, e.g. reactions may include PHD2-Fe²⁺ or PHD2-Fe²⁺-2-OG binding to HIF1α, etc., not just the completely modified PHD2 enzyme.

From the model, we predict several key characteristics about the mechanisms involved in the HIF1 pathway and apply the results to evaluate proposed therapies. Doing so, we provide a molecular representation of hypoxic response that merits further exploration experimentally and computationally. Future modeling studies include representing the independent hydroxylation of HIF1α on its Asn803 residue by factor-inhibiting HIF. To further approximate in vivo conditions, the model will need to represent the effect of pH, such as the acidic conditions found in tumors or in muscle during exercise, and its influence on VHL nucleolus sequestration. HIF2α and the different isoforms of PHD will become a part of the model, initially through modification of the binding rate constants and incorporation of spatial concentration heterogeneity. Incorporating the different sensitivity of the HIF1α degradation binding domains to iron and oxygen may yield a better mechanistic understanding of how the cell copes independently with a drop in iron or oxygen by altering PHD binding (Lee et al., 2005). Recently characterized compounds that influence the hydroxylation of HIF1α and its nuclear accumulation could be added, including OS-9, a protein that binds to HIF1α and the PHDs (Baek et al., 2005), and SUMO-1 protein, which covalently binds to HIF1α and affects its stability and transcriptional regulation (Bae et al., 2004). Although its role has yet to be fully elucidated, p53 binds to HIF1α in anoxia (∼0-0.2% O₂) and can promote its degradation (Fels and Koumenis, 2005); representing this binding in the model may account for the in vitro observation that a maximum in HIF1α nuclear expression occurs at 0.5% O₂ (Jiang et al., 1996). Research into this maximum would consider the relative concentrations of PHD isoforms, and their binding site specificity as a function of O₂ level (Chan et al., 2005), as well as the effect of reactive oxygen species (Schroedl et al., 2002).

It is the authors' anticipation that the computational predictions will stimulate new experiments. An immediate proposal for an in vitro assay would be to quantitatively compare the effects of iron, ascorbate and PHD2 enzyme levels on the hydroxylation rate across a spectrum of O₂ concentrations in a range of cell types. This would give an indication of whether indeed the response to O₂ is of two natures, and whether concentrations of HIF1α co-factors determine if there is a steep drop in hydroxylation or a gradual decrease. Fig. 7C,D provide valuable predictions for therapeutic studies. Experiments to assess the relationship of ascorbate supplementation (e.g. 25, 1000, 2000 μM) in relationship to iron availability at different O₂ levels (e.g. 0, 2, 5, 10, 20%) in cancer (e.g. MDA-MB-435) and endothelial cells (e.g. HUVEC, HBEC) would follow-up on the predictions of the model, and in vitro experiments (Jones et al., 2006; Knowles et al., 2006) to verify in which microenvironments the proposed therapies work effectively, and in what cells these microenvironments are found would be useful. Imaging techniques could be used to assess the variability present in tumor microenvironments - providing a basis for intra-tumor concentrations of compounds in the HIF1 pathway and allowing in vivo application of the computational model.

Different cell types may respond very differently to intracellular hypoxia than the model predicts, and the above experiments could help determine this. Deviations from predictions may reflect unaccounted mechanisms of the HIF1α pathway, and specifically, different reactions in the hydroxylation and degradation of HIF1α. There remains a question of whether the action of ascorbate is as described in the model (reducing Fe³⁺ to Fe²⁺ and binding the PHD2-2-OG-Fe²⁺-O₂ intermediate complex), and whether this varies by cell type and in vivo conditions. The possibility of ascorbate working as a pro-oxidant at low concentrations and an antioxidant at high concentrations intracellularly involves alternate reactions affecting Fe²⁺, H₂O₂, O₂ and HIF1α hydroxylation in the model, which could be assessed using techniques to determine ascorbate and H₂O₂ relationships (Kramarenko et al., 2006), in combination with VHL capture assays and HIF1α quantitative assessment.

The presented computational model adds a novel perspective to understanding the molecular details of how cells sense oxygen. It demonstrates how iron and PHD2 can determine whether the response to an acute hypoxic exposure is a pronounced or gradual change in the amount of HIF1α_hydroxylated. This knowledge is applied to predict the response to proposed therapies targeting the HIF1 pathway. The model provides evidence as to how a changing microenvironment can significantly alter cell susceptibility to drugs that target cofactors in the HIF1α hydroxylation reactions. Even before experimental techniques allow measurements sensitive enough to detect intracellular molecular concentrations changes, a HIF1 pathway computational model can be used to gain an understanding of the molecular mechanisms underlying pre-angiogenic hypoxic response.

From a comprehensive analysis of experimental data, we represent the hydroxylation of HIF1α by PHDs and the ubiquitylation of hydroxylated HIF1α by VHL (Fig. 1). In normoxia, PHD2 is the dominant PHD isoform that hydroxylates HIF1α and determines HIF1α concentrations in a range of cell types (Berra et al., 2003). We modeled the hydroxylation of HIF1α by PHD2 in the cell cytoplasm. The compounds involved in binding to PHD2 in preparation for the hydroxylation of HIF1α include iron, 2-oxoglutarate (2-OG), oxygen and ascorbate. The modified PHD2 then binds and hydroxylates HIF1α. Hydroxylated HIF1α is recognized and ubiquitylated by VBC, the complex that includes VHL bound to Elongins B and C, Cul2 and Rbx1. Equations 1, 2, 3, 4 describe the overall scheme of HIF1α degradation. This includes HIF1α hydroxylation (Eqn 1), independent reactions of iron and ascorbate (Asc) (Eqns 2 and 3), and the binding of HIF1α to VHL (Eqn 4). Table 3 lists the compounds included in the model.

The hydroxylation reactions follow enzyme-substrate binding kinetics. Governing equations are determined from mass balances surrounding the substrate and the intermediate enzyme-substrate complexes. Equation 5 shows an example of the kinetic reaction (Eqn 5A) and accompanying kinetic model (Eqn 5B) for initial steps in PHD2 hydroxylation of HIF1α (Eqn 5B). The full kinetic model can be found in Eqns 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26. A combination of enzyme-substrate saturation assumptions was used for the binding of iron, ascorbate, oxygen and 2-oxoglutarate to PHD2, PHD2 hydroxylation of HIF, and VHL-mediated ubiquitylation. Model inputs are initial compound concentrations, including cellular O₂ levels (Table 1). Output is HIF1α levels in the cell cytoplasm.

In the human type I collagen prolyl hydroxylase reactions with iron, 2-oxoglutarate, oxygen and the peptide substrate (Pro-Pro-Gly)₁₀, the reported K_m values are close to the dissociation constants, and the assumption of their equality (Eqn 7) appears valid (Hieta et al., 2003; Myllyla et al., 1977). The HIF-PHDs and the collagen PHDs share the same hydroxylation reaction sequence and binding cosubstrates. At first glance, it seems reasonable to make the same equilibrium binding assumptions for HIF-PHDs as for the collagen PHDs. However, the overall PHD2 hydroxylation of HIF1α is not at equilibrium and probably unidirectional (Chan et al., 2005; Willam et al., 2004). For the model presented here, the binding was assumed bi-directional between PHD2, iron, 2-oxoglutarate, oxygen (Epstein et al., 2001) and ascorbate, with negligible intermediates formed (Eqn 9) (Goda et al., 2003). For the intermediate complexes of PHD2 and its cofactors, this means the catalytic production terms involving k_cat,Fe2, k_cat,DG, k_cat,O2 and k_cat,AS were set to zero during model runs.

Effects of larger relative k_off rates were also tested, and the current model includes a final irreversible hydroxylation step that considers all three constants, k_on, k_off and k_cat, to be significant (Eqns 19, 20, 21).

Ascorbate is assumed to be a co-reactant with compounds in the hydroxylation reaction, as has been shown with collagen PHDs (Majamaa et al., 1986). Hydroxylation proceeds without ascorbate when there is sufficient iron in the form of Fe²⁺ present (Myllyla et al., 1984; Tuderman et al., 1977). To account for the possibility of the reaction proceeding without ascorbate, the ascorbate-independent binding of HIF1α to PHD2 · Fe²⁺ · 2-OG · O₂ and subsequent hydroxylation is part of the model (third term in Eqn 20). However, eventually without any ascorbate, Fe²⁺ can be oxidized to its Fe³⁺ (or Fe⁴⁺) form, leaving insufficient iron to bind to PHD2 and halting the hydroxylation reaction. The role of ascorbate in the model is to bind to O₂ and reduce Fe³⁺ back to Fe²⁺. We represent the Fe²⁺ reaction with H₂O₂ (Eqns 2, 10 and 22) and Fe³⁺ reduction by ascorbate (Eqns 3, 17 and 22) in the model.

Following hydroxylation, VHL ubiquitylation of HIF1α is likely also not at equilibrium. However VHL binding is notably reversible (k_cat,HαV=0, in Eqns 23, 24, 25, 26). The deubiquitylating enzyme, VDU2, interacts with pVHL, mediates the backward reaction and stabilizes HIF1α (Li et al., 2005). We make the kinetic assumption that proteosome degradation of HIF1α after ubiquitylation is approximately first order in ubiquitylated HIF1α. This simplification seems a valid approximation from reoxygenation experiments (Huang et al., 1998), although it leaves to future studies the exploration of how decay rate varies with the ratio of free to bound or modified HIF1α.

The model rate constants are given in Table 1. Ten of the constants are derived from experimental data (Buettner and Jurkiewicz, 1993; Hirsila et al., 2003; Hirsila et al., 2005; Hon et al., 2002; Kersteen et al., 2004; Lovstad, 2003; Tuckerman et al., 2004). K_m was available from experimental data while individual k_on and k_off were calculated from the model for all but the VHL binding step, where on and off rates were known experimentally. Values for the unknown parameters were estimated as described above and shown in supplementary figures. Default values for kinetic constants and initial conditions are shown in Table 1 with corresponding references. For experiments, O₂ levels are given in percentages, mmHg, or micromolar quantities. The last is reserved for cell culture experiments, where O₂ concentration is calculated in solution. Where appropriate for comparison, we converted model results to the measured experimental units. For cell culture experiments characterizing conditions in cell lysates, values given as percent oxygen or mmHg were equated to micromolar quantities based on an oxygen solubility in water (Tuckerman et al., 2004).

The system of nonlinear differential equations presented in Eqns 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26 was solved using Mathworks Matlab software. The ode23s solver, based on a modified Rosenbrock formula, was used to find a solution for the series of seventeen differential equations. For the time integration, the solver used adjustable time steps with default absolute error tolerance in the solution of 10^-6 μM.

Production terms q_o in Eqn 12 and q in Eqn 20 are nonzero only in chronic hypoxia (>6 hours, Fig. 6). Estimates of these functions are based on experiments (Appelhoff et al., 2004; Berra et al., 2003; D'Angelo et al., 2003).

Ten parameters (five of them independent), were unknown experimentally: k_on,Fe2, k_off,Fe2, k_on,DG, k_off,DG, k_on,O2, k_off,O2,, k_on,AS, k_off,AS, k_on,Hα and k_off,Hα. The following protocol was used to estimate these parameters. Initial rough estimates for k_on and k_off, for the reactions of iron, 2-OG and ascorbate with PHD2, were found from binding of these reactants with other substrates (collagen prolyl hydroxylase, lysyl hydroxylase). Fe²⁺ oxidization and reduction rate constants were estimated respectively from values for the Fenton reaction and from (Buettner and Jurkiewicz, 1996). Sensitivity analysis confirmed estimates for kinetic constants. With the assumption k_cat=0 for intermediate steps, k_on rates were determined by k_off/K_m. In the case of k_on,Hα and k_off,Hα, where k_cat,Hα was a significant value estimated from experiments, the relationship K_m=(k_cat+k_off)/k_on was used.

For all estimated parameters, the kinetic parameter of interest was varied over a minimum range of 1000-fold, while the remaining parameters were held constant. Calculated HIF1α half-lives were compared with experimental data (Berra et al., 2001; Huang et al., 1998; Jewell et al., 2001), to narrow the range of reasonable kinetic parameter values (supplementary material Fig. S4).

The estimated range for the catalytic constant k_cat,Hα (Eqn A12) was relatively small: 0.098-0.164 minute^-1; this corresponds well with a 0.1 minute^-1 approximation from experiments (Tuckerman et al., 2004). Testing ±100-fold differences in k_cat showed the model's high sensitivity to this parameter (supplementary material Fig. S4i). The remaining kinetic parameters estimated from the model were k_offs for the binding of iron (Eqns 10, 11, 12), 2-OG (Eqns 13,14), oxygen (Eqns 15,16) and ascorbate (Eqns 17,18) with PHD2. These were determined using the time of 5 minutes for half [HIF1α]₀ = 1 μM to be hydroxylated (this is based on experiments and the estimate that ubiquitylation of hydroxylated HIF1α takes ⩽3 minutes; see Model validation section). With this assumption, the estimated minimum k_Fe,off is 36 minute^-1; if the 5-8 minute HIF1α half-life range is considered, the range of k_Fe,off is 0.019-36 minute^-1 (supplementary material Fig. S4ii). Kinetic off-rate ranges were also determined for the PHD2 enzyme complex binding to 2-OG (supplementary material Fig. S4iii), 2-OG (supplementary material Fig. S4iv), and ascorbate (supplementary material Fig. S4v). To test the assumption of irreversibility in the binding of modified PHD2 with HIF1α, a range of k_off values were explored (supplementary material Fig. S4vi). The graph shows the feasibility of a nearly irreversible reaction, with k_off,Hα close to zero giving a half-life within the experimental range expected. However the sensitivity of the HIF1α half-life to k_off,Hα was low, and another comparison was needed to help limit its value.

A second means to confirm estimates for kinetic constants was a comparison of PHD2 specific activity with that from Tuckerman et al. (Tuckerman et al., 2004), which estimated a lower limit of 20 mol HIF1α hydroxylated/mol PHD2/minute, as a constant rate over the first 6 minutes in hypoxic MDA-MB-435 cell extracts. Using initial concentrations consistent with those reported in the experiment, [HIF1α]₀=1 μM; [Ascorbate]₀=1 mM; [α-ketoglutarate]₀=1 mM (model, [2-oxoglutarate]₀=1 mM); [FeCl₂]₀=50 μM (model, [Fe²⁺]₀=50 μM); [PHD2]₀=4 nM; [O₂]₀=209 μM (ambient, 21%), and varying the kinetic parameters of interest, the model was compared with experiments to determine minimum estimates for k_cat,Hα, k_off,Fe2, k_off,DG, k_off,O2, k_off,AS and k_off,Hα using the best fit determined by linear regression of the computational curves and least squares analysis comparisons to the data (supplementary material Fig. S5). The model's lower estimates for the kinetic rates corresponded well to the specific activity value from experiments. An exception was seen consistently at two minutes, where lower estimates for the model's kinetic constants underestimated the reported specific activity. Using the range of k_cat,Hα corresponding to the HIF1α half-life of 5-8 min, the higher k_cat,Hα values increased the predicted specific activity for the same k_offs.

After determining estimated kinetic constants, we investigated the sensitivity of the hydroxylation reaction to initial concentrations of reactants (supplementary material Fig. S6). Fig. S1 in supplementary material shows model estimates for PHD2 specific activity and HIF1α half-life using the newly reported K_m,Fe2 value of 0.03 μM, and the corresponding best fit k_on and k_off rates for the other binding reactions. Alternate choices for the set of k_on and k_off rates, using the maximum HIF1α half-life of 8 minutes, resulted in higher estimates for PHD2 specific activity.

Sensitivity analysis was also performed for all parameters found from experiments: K_m,Fe2, K_m,DG, K_m,O2, K_m,AS, K_m,Hα, k_Fe3, k_ASFe, k_on,VL and k_off,VL using the protocol described for estimated parameters (supplementary material Fig. S2 show examples of K_m,Fe2 and K_m,DG). Changes in k_Fe3 and k_ASFe, for the reactions of Fe²⁺ with H₂O₂ and ascorbate, had negligible effects on hydroxylation at default initial conditions; as determined by the model, constants for ubiquitylation (k_on,VL and k_off,VL) have no effect on hydroxylation. Sensitivity analysis shows the switch-like, steep drop in HIF1α hydroxylation at low O₂ levels is a consistent feature (supplementary material Fig. S3). The range of kinetic values used for analysis was 0.2K_m-4K_m, where K_m is the experimental value (Table 1). Comparison of PHD2 specific activity to experiments determined this was an appropriate range (supplementary material Fig. S2ii,iv). An exception was K_m,Fe2, where the range was 0.025K_m,Fe2-4K_m,Fe2 (supplementary material Fig. S2ii). Changes in K_m,Fe2 had minimal effects on the hydroxylation curve (supplementary material Fig. S3i). For K_m,DG, K_m,O2, K_m,AS, and K_m,Hα, increasing K_m values increase the time for hydroxylation, thereby decreasing the steepness of the HIF1α response (supplementary material Fig. S3ii). Increasing K_m values also make the system more sensitive to changes in oxygen at higher O₂ levels.

A script was written in Mathworks Matlab to run the model repeatedly and perform sensitivity analysis over the described parameter range.

This work was supported by NIH Grants HL79653 and HL18292. The authors thank G. Semenza, R. Pili, D. Qian, Z. Bhujwalla and J. Myllyharju for helpful discussions.