A dynamic model of the hypoxia-inducible factor 1α (HIF-1α) network

Hypoxia induces a number of metabolic changes with rapid and profound consequences for cell physiology (Semenza, 2012). The regulation of oxygen homeostasis is a tightly controlled cellular process dependent on the master regulator HIF (hypoxia inducible factor) (Pouysségur et al., 2006; Semenza, 2012) which upregulates adaptive genes by binding to hypoxia response elements (HRE) in their promoters (Semenza, 2003). HIF is composed of an oxygen-regulated α-subunit (HIF-1α, HIF-2α and HIF-3α) and a constitutively expressed nuclear β-subunit (HIF-1β). In normoxia, HIF-1α protein levels are low, due to the action of oxygen-sensitive prolyl hydroxylases (PHDs), which hydroxylate HIF-1α at Pro-402 and/or Pro-564 (Schofield and Ratcliffe, 2004) and target it for ubiquitination-dependent degradation via the Von Hippel-Landau (VHL) protein (Bruick and McKnight, 2001; Epstein et al., 2001). Another level of control lies with the oxygen-sensitive asparaginyl hydroxylase FIH (factor inhibiting HIF), which hydroxylates HIF-1α at Asn-803 and inhibits the recruitment of the transcriptional co-activators p300 and CBP (Ebert and Bunn, 1998; Mahon et al., 2001; Lando et al., 2002; McNeill et al., 2002). In hypoxia, hydroxylase activity decreases, thus enabling HIF-1α to escape degradation and inactivation, translocate to the nucleus to form a transcriptional complex with the HIF1β/ARNT subunit and transcriptional co-activators and initiate gene expression to mount an effective, transcriptionally driven adaptive hypoxic response.

While many of the molecular components of the HIF pathway have been identified and characterized, the dynamics of their interaction within the network are less well understood. Knowing the components of the network is not, in itself, sufficient to understand the complexity of the system. A systems-level model of the HIF pathway will provide a dynamic and mechanistic understanding of how the physical and chemical processes interact to produce a complex cellular response to hypoxia, by suggesting explanations based on biologically plausible mechanisms and making experimentally testable predictions (Arkin and Schaffer, 2011). To this end, a number of mathematical models have been proposed for HIF signaling (Kohn et al., 2004; Kooner et al., 2005; Qutub and Popel, 2006; Yu et al., 2007; Dayan et al., 2009; Schmierer et al., 2010). Kohn and colleagues developed the first theoretical model of the HIF network, which led to a hypothesis that HIF activity behaves in a sharp switch-like manner in response to decreasing gradient of oxygen levels (Kohn et al., 2004). Subsequent models have verified this hypothesis and have provided quantitative explanations of the mechanism responsible for such switching behaviour (Qutub and Popel, 2006; Yu et al., 2007). These models were largely based on experimental data at steady state level, and have generally not considered the dynamics of the hypoxic response due to limitations in the availability of experimental data. Moreover, early models considered only PHD as the oxygen sensor in the HIF network and thus incompletely described regulation of HIF transcriptional activity. More recent models include FIH but were used to explore other features of the network such as the regulation of specific sets of HIF dependent genes (Dayan et al., 2009; Schmierer et al., 2010).

In this study, we have used a combination of mathematical and experimental analysis to understand the HIF-1α signalling network. This is achieved by constructing an iterative dynamic model which is validated by in-house experimental data and which has sufficient predictive power to accurately model the HIF transcriptional response to hypoxia both spatially and temporally. This model incorporates both PHD and FIH as major regulators of HIF-1α activity. It also considers cell compartmentalisation and feedback regulation which were lacking in previous models. In addition to manipulating the oxygen tension, we have used the pharmacological inhibitors dimethyloxalylglycine (DMOG) and JNJ-42041935 (abbreviated to JNJ1935 in this study), both of which are cell permeable 2-oxoglutarate analogues which inhibit the hydroxylation of HIF-1α by displacing the endogenous 2-oxoglutarate cosubstrate required by PHDs and FIH (Mole et al., 2003; Barrett et al., 2011). Importantly, JNJ1935 is a more potent and selective inhibitor of PHDs than FIH (Barrett et al., 2011), thus enabling the pharmacologic dissection of the activities of PHDs and FIH. Furthermore, the model captures the temporal dynamics of HIF transcriptional activity in response to hypoxia or to pharmacological inhibitors through the use of a secreted Gaussia luciferase reporter under the control of HIF (Bruning et al., 2012).

Thus, employing a systems approach of iterative experimentation and mathematical modelling, we have developed a spatiotemporal model of the HIF-1α signalling network which accurately predicts biological behaviour and have used this model to generate testable hypotheses. The model can distinguish between prolyl and asparaginyl hydroxylation and predicts that HIF-1α protein stabilisation does not always correlate with HIF transcriptional activity. Instead, the removal of the asparaginyl hydroxylation step is necessary for HIF-1α activity. Furthermore, through sensitivity-analysis of the model and qualitative analysis of the core module, we propose that asparaginyl hydroxylation confers upon HIF-1α resistance to proteosomal degradation. Using a combination of in vitro experimentation and in silico predictions, we confirm the network topology of the hypoxia response, establish the wirings controlling the dynamic regulation of HIF-1α transcriptional activity by hydroxylases and use the model to offer biologically plausible explanations to counter-intuitive experimental observations.

To provide a quantitative framework for understanding the HIF pathway, we have developed a dynamic, ordinary differential equations-based model from the validated and published core components of the HIF-1α network (Fig. 1). This model integrates our current understanding of the interactions between the known HIF-1α pathway components. The readout of the HIF-1α transcriptional activity used in our model is from the Gaussia luciferase reporter under the control of HIF (Bruning et al., 2012). Our model-based study follows an integrated cycle of model-driven hypothesis generation and experimental validation strategy which has not previously been reported for the HIF pathway.

The dynamic model incorporates key molecular interactions in the HIF-1α pathway and the molecular components and steps of the model are described in Fig. 1. Detailed discussion of the model including reactions, reaction rates and parameters are given in the Materials and Methods and in supplementary material Tables S1–S3. In normoxia, HIF-1α mRNA and protein are produced at a steady rate, but the protein is degraded by either non-specific protein turnover or by the oxygen-sensitive prolyl-hydroxylation, leading to Von Hippel-Lindau (VHL)-mediated proteosomal degradation (Kaelin, 2005). Since prolyl-hydroxylated HIF-1α protein is quickly ubiquitinated by VHL, we assume the prolyl-hydroxylation step is irreversible (Chan et al., 2005). Although the PHD isoforms have specific cytoplasmic or nuclear localisation (Metzen et al., 2003), they are considered as one entity for simplicity, and are assumed to be present in both compartments. HIF can also be inactivated by FIH through asparaginyl-hydroxylation (Mahon et al., 2001; Lando et al., 2002). Since FIH was also found in both compartments but mainly in the cytoplasm (Metzen et al., 2003), our model assumes a higher level of FIH in the cytoplasm than in the nucleus. We also assume that HIF-1α hydroxylation mediated by PHD and FIH as well as VHL-induced HIF degradation can occur in both cytoplasmic and nuclear compartments with comparable kinetics. We assume that asparaginyl-hydroxylated HIF-1α can be subsequently prolyl-hydroxylated, but prolyl-hydroxylated HIF is unlikely to be subsequently asparaginyl-hydroxylated due to the rapid VHL-mediated degradation of prolyl-hydroxylated HIF-1α (Chan et al., 2005). The hydroxylation steps are assumed to be mostly irreversible, although a small fraction could in theory be dehydroxylated (Lancaster et al., 2004).

Hydroxylated or unhydroxylated HIF-1α protein is also assumed to shuttle to the nucleus via the importin pathway (Depping et al., 2008). In the nucleus, unhydroxylated HIF-1α can dimerise with HIF-1β to form a transcriptional complex (Jiang et al., 1996). We did not include a separate step for the binding to the transcriptional co-activators p300/CBP and assume that the formation of the HIFα/β dimer results in a transcriptionally active complex. This complex can bind to hypoxia response elements (HRE) in the promoter region of hypoxia-responsive genes, such as those on the Gaussia reporter gene which is used to assay temporal HIF transcriptional activity (Bruning et al., 2012), and upregulate their expression. A negative feedback loop consisting of an upregulation of the HIF-regulated PHD (Marxsen et al., 2004; Stiehl et al., 2006; Minamishima et al., 2009) is also included. Furthermore, as PHDs are induced by hypoxia (Stiehl et al., 2006) and their levels can change dynamically, we assume PHDs can translocate in and out of the nucleus (Jokilehto et al., 2006; Steinhoff et al., 2009; Pientka et al., 2012).

HIF activation is well characterised in cells exposed to hypoxia (Semenza, 2006). Densitometric analysis of HIF-1α western blots show that the HIF protein was rapidly and transiently stabilised at 3% and more at 1% oxygen tensions (Fig. 2A,B), resulting in a corresponding HIF-1α transcriptional activity as measured in human embryonic kidney cells (HEK293) transfected with pGluc-HRE (Fig. 2C). No activity was detected in normoxia (21% oxygen tension). Thus we confirm previous studies demonstrating that HIF protein is stabilised in low oxygen tension, resulting in the corresponding transcriptional activity as reported by pGluc-HRE.

No obvious stabilisation of HIF-2α was observed in HEK293 cells in response to hypoxia over 12 hours (supplementary material Fig. S1). Therefore, in this model, the vast majority of the HIF-dependent transcriptional effects can be ascribed to the HIF-1 isoform. The HIF-1α stabilisation and transcriptional activity data was fitted to the HIF-1α model to estimate model parameters (supplementary material Table S3). Detailed discussion of the calibration of the model to oxygen level is given in the Materials and Methods. Fig. 2D shows simulations for the HIF-1α transcriptional activity in 21%, 3% and 1% oxygen which closely matched the experimental data (Fig. 2C cf. 2D). Thus our model can accurately predict HIF-1α transcriptional activity when the oxygen tension is changed.

A prediction of the model is that inhibition of both PHD and FIH are required for optimal HIF-1α-dependent transcriptional activity. The percentage of active PHD and FIH is decreased in the model to reflect an inhibition of either PHD alone or both PHD and FIH. Simulations for HIF stabilisation and transcriptional activity under these conditions are shown in Fig. 3A–D. The model predicts that HIF stabilisation is dependent on the abundance of active PHD (Fig. 3A,B) and that a reduction in active FIH is needed in order to observe a corresponding increase in HIF-1α transcriptional activity (Fig. 3C,D).

To test these predictions, we use two different hydroxylase inhibitors: DMOG, a pan-hydroxylase inhibitor; and JNJ-42041935, a prolyl-hydroxylase selective inhibitor (abbreviated to JNJ1935 in this paper). Low concentrations of JNJ1935 were previously shown to selectively inhibit PHDs, whereas higher concentrations inhibit all hydroxylases, including FIH (Barrett et al., 2011). For both inhibitors, we observed a dose-dependent and time-dependent stabilisation of HIF-1α protein (Fig. 3E). Densitometric analysis of the HIF-1α western blots show that 100 µM JNJ1935 stabilised HIF-1α at the same level as 1 mM DMOG in stabilising HIF after 12 hours (Fig. 3F,G). Consistent with the model prediction of PHD and FIH inhibition, we observed a HIF-1α transcriptional activity in response to increasing concentration of DMOG in HEK293 cells transfected with pGluc-HRE to (Fig. 3H). Furthermore, we observed that the transcriptional activity for cells exposed to 100 µM JNJ1935 is much smaller than those exposed to 1 mM DMOG (Fig. 3I cf. 3H), consistent with the model prediction (Fig. 3C,D). We also observed a stronger induction of the PHD2 protein [a HIF regulated gene (Metzen et al., 2005)] in hypoxia (1% O₂) and DMOG than in JNJ1935 (supplementary material Fig. S2). No effect from the DMSO solvent was observed on Gaussia activity (supplementary material Fig. S3). Thus we show that the predictions can be experimentally validated and our model supports a role for asparaginyl hydroxylation by FIH as an important step in regulating the transcriptional activity of stabilised HIF-1α.

Next we make predictions for the temporal effects of DMOG and JNJ1935 on the HIF-1α stabilisation and transcriptional activity. Both DMOG and JNJ1935 inhibitory efficiencies were chosen to generate a maximal level of HIF-1α stabilisation, and the model predicts similar patterns of stabilisation (Fig. 4A). However the transcriptional activities are predicted to be much less with JNJ1935 treatment than with DMOG treatment (Fig. 4B).

To test these predictions, we used DMOG and JNJ1935 at a concentration of 1 mM and 100 µM respectively which produced similar levels of HIF-1α protein stabilisation (Fig. 3F,G). The experimental data from HEK293 cells transfected with pGluc-HRE and treated with DMOG and JNJ1935 over a 12-hour period matched the predicted outcome (Fig. 4D cf. 4B). Thus we show that our mathematical model can predict the temporal dynamics of hydroxylase inhibition by DMOG and JNJ1935 on HIF-1α transcriptional activity.

FIH has previously been shown to have higher affinity for oxygen than PHD (Koivunen et al., 2004; Stolze et al., 2004). Our model predicts that a reduction in active FIH will result in an increase in HIF-1α activity (Fig. 5A) while an increase in FIH will cause a decrease in activity at 1% oxygen (Fig. 5B).

To experimentally test these predictions, we used either siRNA against FIH (FIH siRNA) or an overexpression plasmid for FIH (pcDNA3-FIH) in HEK293 cells transfected with pGluc-HRE. We show that we can manipulate the level of FIH protein expressed (Fig. 5C). Experimental data in cells exposed to 1% oxygen confirmed the predictions that that knocking down FIH increases the HIF-1α activity (Fig. 5D), while overexpressing FIH decreases the activity (Fig. 5E). No effect was observed when the experiments were performed in normoxia (supplementary material Fig. S4).

We next used the model to make predictions about the role of FIH during PHD inhibition using the pharmacological inhibitor JNJ1935. Our model simulations predict that, with JNJ1935 treatment, silencing FIH will result in a downregulation of total HIF-1α stabilisation while at the same time lead to an upregulation of HIF-1α transcriptional activity (Fig. 6A,B). This prediction was tested and validated experimentally (Fig. 6C,D). Thus, from mathematical predictions and experimental validation, we show that under JNJ1935, siRNA against FIH causes a reduction in HIF-1α stabilisation and an increase in HIF-1α transcriptional activity.

Next, we aimed to provide a possible explanation for these paradoxical observations described in Fig. 6 and identify a likely molecular mechanism. Given that JNJ1935 inhibits HIF-1α degradation mediated by PHD at 2 orders of magnitude higher than DMOG in vitro (Barrett et al., 2011), we assume that any upregulated prolyl hydroxylases will be inhibited and not contribute to the observed decrease in HIF-1α stabilisation. We isolated from the full model a core HIF-1α signalling module (Fig. 7A) which consists of unhydroxylated and asparaginyl-hydroxylated HIF-1α in the cytoplasm and in the nucleus (more details in the Materials and Methods). Prolyl-hydroxylated HIF-1α is assumed to be absent in the core module (the grey arrows in Fig. 7A) as the system is under JNJ1935 inhibition. The core module describes the dynamic transfer between these HIF-1α moieties and is thus sufficient in studying the qualitative behaviour of the system under JNJ1935-mediated inhibition.

In order to account for HIF-1α degradation by PHD-independent mechanisms (Kong et al., 2007; Liu et al., 2007; Koh et al., 2008), our full model currently assumes that only the cytoplasmic HIF-1α is subjected to direct protein turnover for simplicity (Reactions 1 and 2 in Fig. 1, Fig. 7A). To identify potential factors that may influence the behaviour of the total HIF-1α stabilisation and transcriptional activity in the core module, we carried out a parameter sensitivity analysis in which we systematically perturbed each reaction numbered in Fig. 7A by increasing or decreasing its rate by 10-fold and compare simulations under JNJ1935 with or without silencing FIH (supplementary material Figs S6, S7). We then compared the responses of total HIF-1α stabilisation and transcriptional activity to that of the unperturbed case. Although the perturbations appeared to alter the levels of total HIF-1α stabilisation and activity quantitatively, none of them resulted in a qualitative change. Specifically, we always observed a reduction of HIF-1α stabilisation and an increase in HIF-1α transcriptional activity under reduced FIH under all perturbations, including varying the HIF-1α abundance, consistent with the data (supplementary material Figs. S6, S7; Fig. 6). Similar observations were seen even when we varied the fold change in reaction rates up to 20 (data not shown), thus suggesting that the perturbed reactions are not essential in controlling the qualitative behaviour of HIF-1α under FIH silencing.

When we relaxed our assumption of oxygen-independent degradation and assumed that other forms of HIF-1α (HIF-1α_n, HIF-1α-aOH and HIF-1α_n-aOH) are also susceptible to such degradation (Reactions 2b–d in Fig. 7A,B), we found that, as long as the degradation rate of the asparaginyl-hydroxylated HIF-1α (rate [2c + 2d]) was less than that of the unhydroxylated HIF-1α (rate [2+ 2b]), silenced FIH would always lead to decreased HIF-1α stabilisation (Fig. 7C) but increased HIF-1α activity (Fig. 7D) as observed in our predictions and our experimental data (Fig. 6). On the other hand, if the asparaginyl-hydroxylated HIF-1α degraded faster than unhydroxylated HIF-1α by oxygen-independent mechanisms, reduced FIH would lead to increased HIF-1α stabilisation (Fig. 7E) and increased HIF-1α activity (Fig. 7F) instead. Analytical derivations further confirmed these findings mathematically, showing they were true regardless of the rate values of other reactions (see Materials and Methods for detailed derivation). Thus, analysis of the core signalling module reveals that the differential rates in PHD-independent degradation of the unhydroxylated HIF-1α and asparaginyl-hydroxylated HIF-1α directly control the qualitative response of HIF stabilisation and activity under JNJ1935 inhibition. Importantly, the finding suggests that asparaginyl-hydroxylated HIF-1α is less prone to degradation by oxygen-independent mechanisms, resulting in a model that better fits our experimental data.

We next tested the qualitative simulations from the core module by inhibiting proteosomal degradation in HEK293 cells transfected with siRNA against FIH and under PHD inhibition from JNJ1935. The use of MG132 increased the level of HIF-1α stabilisation (Fig. 7G), thus showing that HIF-1α is susceptible to proteosomal degradation. Furthermore, we also tested our simulations on HEK293 cells transfected with a full length HIF-1α that includes mutations to both proline residues required for hydroxylation by PHD [Dmut-Pro; (Hagen et al., 2003)] and siRNA against FIH. This Dmut-Pro HIF-1α is constitutively stable and can be detected by antibodies against its V5-tag or against HIF itself. Silencing FIH reduced its expression, and MG132 could reverse this decrease (Fig. 7H).

Thus, using a core module of the HIF-1α model and experimental data, we show that the decrease in HIF stabilisation due to FIH silencing in JNJ1935-treated cells arises from non-PHD dependent proteosomal degradation of unhydroxylated HIF. Our data suggests that FIH confers protection to HIF-1α from proteosomal degradation.

Here we provide a data-driven iterative mathematical model which describes the dynamics of HIF-1α signalling at a detailed mechanistic level. We validated the model using experimental data and generated predictions which were subsequently pursued and validated by experiments. Furthermore, we use the predictions made by the model, supported by experimental data, to propose that hydroxylation by FIH confers a degree of protection to the HIF-1α protein from proteosomal degradation.

Our predictive and adaptive model predicts testable hypothesis that are in agreement with experimental data, thus confirming the generally accepted HIF-1α network topology. It is supported by iterative experimental data on the HIF-1α stabilisation and transcriptional activity dynamics. The later is a key end point of the signalling pathway, when genes are switched on in response to hypoxia. This crucial step has been poorly represented in previous models due to the limited amount of experimental data. Our model uses data from the Gaussia luciferase reporter which provides temporal information on the HIF-1α transcriptional activity from a population of cells (Bruning et al., 2012). Furthermore, the use of the pharmacological hydroxylase inhibitors DMOG and JNJ1935 enabled use to selectively dissect out the role of PHDs and FIH. Importantly, we show that HIF-1α stabilisation does not directly correlate to HIF-1α transcriptional activity due to the inhibitory activity of FIH. Our model also includes a negative feedback loop through the upregulation of PHD by HIF-1α (Stiehl et al., 2006) to reflect the decrease in HIF-1α protein observed after 6–8 hours in hypoxia. The miR-155 negative feedback loop was not considered in this model as this is proposed to occur during prolonged hypoxia (Bruning et al., 2011). Our model assumes that only an unhydroxylated form of HIF-1α can be transcriptionally active as prolyl-hydroxylation will result in rapid VHL-mediated degradation (Chan et al., 2005) and asparaginyl-hydroxylation will prevent binding of the transcriptional co-activators p300 and the CREB binding protein (CBP) (Ebert and Bunn, 1998; Lando et al., 2002). Our data also supports theses assumptions, because we observe that the HIF-1α transcriptional activity is more reduced in the presence of JNJ1935 (which inhibits only PHDs) than in the presence of the pan-hydroxylase inhibitor DMOG.

Inhibition of FIH increases HIF-1α transcriptional activity in hypoxia while exogenous delivery of FIH reduces it, thus confirming the role of FIH in controlling the HIF-1α transcriptional response. Moreover, silencing FIH during inhibition of prolyl-hydroxylation reduces HIF-1α protein stabilisation but paradoxically increases HIF-1α transcriptional activity, for which we propose a novel role for FIH in conferring protection from non-PHD-mediated degradation. A possible explanation for the observed decrease in HIF-1α stabilisation could be the negative feedback upregulation of PHD enzymes. However this is unlikely as the cells were incubated in JNJ1935, which would inhibit any upregulated prolyl hydroxylases, given that its potency is of 2 orders of magnitude over DMOG (Barrett et al., 2011). Furthermore, we can also observe this FIH-mediated protection from proteosomal degradation using HEK293 cells transfected with a HIF-1α construct where both its proline residues for PHD hydroxylation were mutated (Hagen et al., 2003). Furthermore, the proteosomal inhibitor MG132 reversed the decrease in HIF stabilisation, indicating the decrease is at the protein level.

Thus we have a counter-intuitive observation: both the model and the experimental data show that a decrease in HIF-1α protein occurs parallel to an increase in HIF-1α activity when the PHD enzymes are inhibited and when FIH is silenced. A comparable decrease in HIF-1α protein stabilisation was observed in HeLa and U-2OS cells exposed to hypoxia and transfected with siRNA against FIH (Stolze et al., 2004), which was explained through upregulation of PHD2 via FIH activity. In our study, we go further and explore the effect of FIH activity in an environment where PHD enzymes are inhibited with JNJ1935 (Barrett et al., 2011). Through parameter sensitivity analysis, we have narrowed down the reactions to the degradation rate of unhydroxylated and asparaginyl-hydroxylated HIF-1α protein. Using a reduced core model, we show that our model simulation can match the experimental observation only if the rate of degradation of unhydroxylated HIF-1α protein is higher than the rate of degradation of asparaginyl-hydroxylated HIF-1α. We also provide mathematical proof that this is the most likely explanation. We propose that loss of FIH makes it more susceptible for HIF-1α to be degraded via non-PHD-mediated degradation. The most likely explanation is that asparaginyl-hydroxylation confers protection to this degradation. We speculate that the asparaginyl-hydroxyl residue on HIF-1α might be interfering a non-PHD-mediated degradation mechanism, of which several have been proposed, such as RACK1 (Liu et al., 2007), COMMD1 (van de Sluis et al., 2009) and DEC2/SHARP1 (Montagner et al., 2012) dependent degradation. It is probably not dependent on the interaction of FIH to VHL (Mahon et al., 2001), as VHL was not necessary for non-PHD degradation (Kong et al., 2007; Liu et al., 2007; Koh et al., 2008; van de Sluis et al., 2009).

In summary, our study illustrates how biological experiments coupled to mathematical modelling can synergize to provide a better understanding of a complex signalling pathway. The process is iterative and incremental; leading to a robust predictive and adaptive model whose predictions can be tested experimentally. Our in vitro data and in silico predictions provide new insights into the molecular mechanisms linking hypoxia to gene expression. Specifically, we propose that non-PHD-mediated degradation is an important step in controlling the HIF-1α response during hypoxia, and that FIH confers protection from this degradation mechanism.

Human embryonic kidney cells (HEK293) were cultured in DMEM high-glucose medium supplemented with 10% FCS and 100 U/ml penicillin-streptomycin. Human hepatocellular carcinoma cells (HepG2) were grown in minimum essential medium containing 10% FCS, 2 mM L-glutamine, non-essential amino acids, and 100 U/ml of penicillin-streptomycin. Cells were exposed to hypoxia using pre-equilibrated media and maintained in standard normobaric hypoxic conditions (1% or 3% O₂, 5% CO₂ and 94% or 92% N₂) in a hypoxia chamber (Coy Laboratories, Grass Lake, Michigan, USA). Normoxic controls were exposed to pre-equilibrated normoxic media and maintained at atmospheric O₂ levels (21% O₂, 5% CO₂) in a tissue culture incubator.

Whole-cell extracts were generated in either normoxia or hypoxia according to previously published protocol (Agbor et al., 2011). Protein concentration was quantified using a Lowry assay (Bio-Rad, Hertfordshire, UK), and samples were normalised accordingly. Samples were separated by SDS-PAGE and immunoblotted as described previously (Bruning et al., 2011) using the following primary antibodies and dilutions: HIF-1α (1∶250; BD Pharmingen, Oxford, UK), PHD1 (1∶1000, Novus Biological, Cambridge, UK), PHD2 (1∶1000, Novus Biological, Cambridge, UK), PHD3 (1∶1000, Novus Biological, Cambridge, UK), FIH (1∶2000, Abcam, Cambridge, UK) and β-actin (1∶10,000; Sigma, Wicklow, Ireland).

All transfection was performed using Lipofectamine 2000 (Invitrogen, Paisley, UK) according to the manufacturer's protocol. A custom FIH siRNA (Eurofins MWG Operon, Ebersberg, Germany) using published sequence (F1) (Cockman et al., 2006) was used to downregulate FIH expression. A pcDNA3-FIH-V5 tag construct (kind gift of Prof Eric Metzen, University of Duisburg-Essen) was used to overexpress FIH. Confirmation of knockdown or overexpression was obtained by immunoblotting for FIH protein. A full length HIF-1a construct with mutations at both proline residues required for hydroxylation by PHD (pcDNA3-HIF-DM-Pro-V5 tag, abbreviated to Dmut-Pro, kind gift of Prof Thilo Hagen, National University of Singapore). The pGluc-HRE Gaussia luciferase vector (containing four copies of the EPO HREs in the right orientation) was transfected into HEK293 cells as previously described (Bruning et al., 2012). The mammalian expression vector pGluc-TK (NEB) contains the coding sequence for Gaussia luciferase under the control of the Herpes Simplex Virus thymidine kinase (TK) promoter for constitutive activity. It does not contain the EPO HRE sequence. Gaussia luciferase activity was measured using the Biolux Gaussia luciferase Flex Assay kit (NEB, Hertfordshire, UK) in a plate reader (Synergy HT, Biotek, Bedfordshire, UK).

The cell permeable pan-hydroxylase inhibitor DMOG (Cayman Chemicals, Michigan, USA), the prolyl-hydroxylase inhibitor JNJ-42041935 [abbreviated to JNJ1935 in this paper; kind gift of Dr Mike Rabinowitz (Janssen Research & Development, LLC)] and the proteosomal inhibitor MG132 (Sigma, Wicklow, Ireland) were dissolved in DMSO (DMSO; Sigma, Wicklow, Ireland).

As many of the kinetic parameters in the HIF/PHD/FIH system are not measured at the present time, we needed to estimate unknown model parameters using a total of six experimental data sets. These include HIF stabilisation data under DMOG and JNJ1935 hydroxylase inhibitors, and transcriptional activity data under DMOG, JNJ1935, 1% and 3% oxygen tension, which constrained the model parameters well. Knowledge of the rate constants of reactions and initial concentrations of model species were required to describe the temporal (and steady-state) behaviours of the system. Where possible, estimates for model parameters were obtained from the literature (supplementary material Table S3). In the absence of such information, the kinetic rate constants and initial concentrations were set to intermediate values within physiologically plausible ranges before fitting so as to optimise model performance. Specifically, the association and dissociation rates were restricted to be within the typical ranges for protein-protein interactions. The association of protein molecules into dimers or larger complexes occurs with typical rate constants of the order of 10⁻⁴ to 10⁻¹ nM⁻¹ s⁻¹ (Eltis et al., 1991; Northrup and Erickson, 1992; Kholodenko et al., 1999). In addition, the reaction rates were always constrained not to be faster than the diffusion limit. The fitted parameter set is given in supplementary material Table S3. Although parameter fitting constitutes an important step in model development, it is worthwhile to note that the prime purpose of computational modelling is to provide a basis for guiding experimental analysis and testing explicit hypotheses; a model by itself is not an objective ‘truth’, but it can be used to falsify or confirm a specific hypothesis.

For realistic simulations of the hypoxia and inhibitors experiments, we employed the step function, which can be used to describe the sudden depletion of either oxygen level or PHD and/or FIH levels due to the inhibitors (DMOG and JNJ1935) while the system is in normoxia condition prior to the start of the experiments.

As a result, O₂ is replaced by Hypo3(t) or Hypo1(t) in the reaction rates v₃, v₉ and v₁₁ in supplementary material Table S1 to describe a model for 3% and 1% oxygen tension, respectively.

To model experiments using inhibitors, we assume that DMOG efficiently reduces both PHD and FIH levels to very low level (1% of the abundance of active enzymes prior to treatment) and that JNJ1935 reduces PHD to low level (10% of its pre-treated level) while FIH is slightly affected (reduced to 80% of its pre-treated level (Barrett et al., 2011). siRNA and overexpression experiments were modelled by adjusting the initial value of a species (FIH for example) to the measured extent of protein depletion or protein overexpression as determined by quantitative western blotting.

Here we consider the core module of the HIF signaling network extracted from the full model as explained in the main text. This reduced model retains the essential features of the HIF-1α dynamic signaling under JNJ1935 where PHD-mediated degradation of HIF-1α is prevented due to PHD inhibition. Thus, the reduced model can be used to analyse the dynamical properties of HIF-1α signalling under this particular inhibitor. In addition to simulations, we present the analytical derivation showing that the differential rates in oxygen-independent degradation of the unhydroxylated HIF-1α and asparaginyl–hydroxylated HIF-1α directly control the qualitative response of HIF stabilization and activity under JNJ1935 inhibition.

Now, consider three scenarios:

(i) k₂ = k₂*: degradation rate [HIFα]_total = degradation rate [HIFα-aOH]_total

Since the right-hand side of Eqn 1 depends only on HIF_tot, HIF_tot is not affected by change in any parameters other than k₁ and k₂. As a result, silencing FIH by reducing both FIH and FIH_n would not change the total HIF stabilisation but increase HIF activity due to more HIFα_n are available, as illustrated in supplementary material Fig. S8.

(ii) Finally k₂>k₂*: degradation rate [HIFα]_total >degradation rate [HIFα-aOH]_total

In this case, the last term of Eqn 1 is negative. When FIH is silenced, the sum (HIFα + HIFα_n) is increased, but in contrast to case (ii) this makes the right-hand side of Eqn 1 smaller due to the negative factor (k_2' − k₂) and so the total HIF stabilisation is decreased. This can be seen in the simulation given in Fig. 7C,D.

(iii) Finally k₂<k₂*: degradation rate [HIFα]_total<degradation rate [HIFα-aOH]_total

In this case, the last term of Eqn (1) is positive. When FIH is silenced, less HIFα and HIFα_n are converted into HIFα-aOH and HIFα_n-aOH, thereby increasing the sum (HIFα + HIFα_n), making the right-hand side of Eqn 1 bigger and so increasing the total HIF stabilisation. Moreover, since HIFα_n is larger, HIF activity increases. This can be seen in the simulation given in Fig. 7E,F.

All experiments were carried out a minimum of n = 3 independent times unless otherwise indicated and data are expressed as the mean±s.e.m. Statistical significance was calculated by one-way ANOVA followed by Bonferroni's multiple comparison test in Prism (GraphPad, California, USA).