Hypoxia induces the expression of genes that alter metabolism through the hypoxia-inducible factor (HIF). A theoretical model based on differential equations of the hypoxia response network has been previously proposed in which a sharp response to changes in oxygen concentration was observed but not quantitatively explained. That model consisted of reactions involving 23 molecular species among which the concentrations of HIF and oxygen were linked through a complex set of reactions. In this paper, we analyze this previous model using a combination of mathematical tools to draw out the key components of the network and explain quantitatively how they contribute to the sharp oxygen response. We find that the switch-like behavior is due to pathway-switching wherein HIF degrades rapidly under normoxia in one pathway, while the other pathway accumulates HIF to trigger downstream genes under hypoxia. The analytic technique is potentially useful in studying larger biomedical networks.
Molecular oxygen is the terminal electron acceptor in the mitochondrial electron transport chain. Hypoxia, or oxygen deficiency, induces a number of metabolic changes with rapid and profound consequences on cell physiology. A hypoxia-induced shortage of energy alters gene expression, energy consumption, and cellular metabolism to allow for continued energy generation despite diminished oxygen availability. A molecular interaction map of the hypoxia response network has been proposed [1–3] on the basis of analyzing conserved components between nematodes and mammals. The key element in this network, hypoxia-inducible factor (HIF), is a master regulator of oxygen-sensitive gene expression [4–6]. HIF is a heterodimeric transcription factor which consists of one of the three different members (HIF-1α, HIF-2α, and HIF-3α) and a common constitutive ARNT subunit which is also known as HIFβ. The system also includes an enzyme family: prolyl hydroxylases (PHDs), which directly sense the level of oxygen and hydroxylate HIFα by covalently modifying the HIFα subunits. It is very likely that reactive oxidative species (ROS), which are a byproduct of mitochondrial respiration, are also involved in oxygen sensing by neutralizing a necessary cofactor, Fe2+, for the hydroxylation of HIFα by a PHD [7–10]. There are three members in this enzyme family: PHD1, PHD2, and PHD3. The hydroxylated HIFα is then targeted by the von Hippel-Lindau tumor-suppressor protein (VHL) for the ubiquitination-dependent degradation. Hypoxia response element (HRE) is the promoter of the hypoxia-regulated genes, and the occupancy of HRE controls the expression levels of these genes. The cascade in Figure 1 (reproduced from Figure 2 of [1]) consists of an input (the concentration of oxygen) and an output (the activation of promoters that are under control of HREs) as the core network. The network is characterized by a switch-like behavior, namely the sharp increase of HIFα when the oxygen decreases below a critical value, followed by a sharp increase of HRE occupancy. It was observed experimentally on many cell lines including Hela cells [11] and Hep 3B cells [12] that HIFα increases exponentially as the oxygen concentration decreases. The past two decades have seen a growing body of work on the use of mathematical modeling to help uncover both general principles behind molecular networks and to provide quantitative explanations of particular network phenomena [13] that may one day have sufficient predictive power to accurately model large subnetworks of the cell. In this sense, Kohn et al. [1] have successfully modeled the switch-like response characteristics of HRE occupancy, by numerically integrating a system of ordinary differential equations (ODEs) involving a score of molecular species related to hypoxia. The large model, however, does not identify the smaller components that are actually responsible for the switch-like response and that may occur in other such networks. Furthermore, a numerical solution does not provide the type of insight that mathematical formulas can. At the same time, it is virtually impossible to solve symbolically the type of nonlinear differential equations that model reactions. In this context, methods are desirable that are both tractable, that reduce a system to its key components, and that are not solely reliant on numerical solution. Extreme pathway analysis (EPA) is one such recently developed method [14–16]. In this method, the dynamics of interactions between species are formulated as a Boolean network in which the state of a gene is represented as either transcribed or not transcribed. Upregulation and downregulation of genes are captured through an appropriate sign (plus or minus) and a scaling constant. The Boolean network is then formulated as a matrix of interaction rules that is then analyzed to help reveal key components and their contributions to the dynamic behavior [16]. The theory of matrices then allows us to look for vectors that characterize the matrix in ways that are helpful for further analysis. The EPA technique, in particular, finds vectors (extreme pathways) that correspond to the boundaries of the space of steady-state solutions to the differential equations. We note that similar methods, such as flux balance analysis (FBA) and elementary modes analysis, have been developed in other contexts [17–19]. They essentially yield the same results [18], which have been verified by ExPa [20] and CellNetAnalyzer [21,22]. These methods provide a way out of the intractable complexity of sizable molecular networks [23–26]. Our contribution is to go beyond this type of matrix approach and provide a detailed quantitative analysis that explains the observed behavior in the models. This is achieved by combining elementary pathway identification via EPA, which depends solely on the network topology, and the detailed analytical as well as numerical analysis of the governing differential equations in the model, which allows studies of the phase space spanned by the mostly unknown rate constants in the differential equations. Specifically, EPA is first used in our approach to decompose the original network into several underlying pathways. Following this, we make some reasonable approximations to facilitate analytic solution. We show that this analytic solution, in the case of the hypoxia network, explains the switch-like behavior. This explanation is confirmed by comparing the numerical output of the simplified model with the numerical output of the complete (and complex) differential equation model. A second contribution of this paper is to highlight a particular mechanism of pathway-switching or pathway branching effect [27] that appears to cause the sharp response to oxygen concentration. In particular, we examine the flux redistribution among the elementary pathways as a function of oxygen concentration. We also identify the key molecular species involved in the subcomponent of the network and show quantitatively how the response of this subcomponent exactly matches the overall response and thus is responsible for it. For hypoxia, our analysis suggests that the cycle of abundant production and efficient degradation of HIFα plays the main role in the sharp response.
For consistency and ease of understanding, we use the notation and nomenclature in [1] and use their 23-species network and differential equation model as the starting point. With this background, the original network shown in Figure 2 of [1] can be further reduced in the following way. Kohn et al. [1] have shown that the feedback of mRNA is not necessary for the switch-like behavior. We therefore eliminate this feedback loop (reaction k32). Hence, Transcript intermediate 1,2, and 3 (Species 8,9, and 10) can also be dropped, as well as the associated reactions: k7, k8, k9, k10, and k11, because they do not affect the dynamics of the network. Species 23 is only the joint name of HIFα: ARNT: HRE (Species 7) and HIFαOH: ARNT: HRE (Species 22); therefore, it is dropped. HIFα precursor (Species 1) is a constant and is thus dropped, because its information can be simply encoded in the reaction k1. The degradation products (Species 2) are also eliminated because they are assumed to leave the network immediately after their production and do not affect the dynamics. Similarly, species 19,20, and 21 do not contribute to the dynamics and are therefore removed. The resulting network is summarized in Tables 1 and 2, where there are 13 molecular species and 19 reactions in total. The system can be described by the following set of ODEs where [Sn] stands for the concentration of species n as tabulated in Table 1 and [O2] indicates the input cellular oxygen concentration. Table 2 shows the specific reaction each rate constants kn represents. The real values of kn are from [1]. Note that the ODE system below is typical: the terms are based on mass-action principles and, taken together, result in complex behavior not readily discernible by examining the equations. We also dropped the precursor concentration [S1] since it is set to unity. This section assumes some familiarity with linear algebra. The 13 × 19 stoichiometric matrix Φ of the reduced hypoxia response network (Tables 1 and 2) is shown in Figure 2 (for details, see Materials and Methods). The rank of Φ is computed and shown to be 9, indicating that there are only nine independent molecular species to serve as constraints for the analysis. Therefore, the dimension of the corresponding null space is 10. The linearly independent basis B vectors are generated by Matlab 6. 5 and are shown in Figure 3. According to the constraint that no negative values are allowed in the basis vectors, we can uniquely transform basis B into basis P as shown in Figure 4. Both b8 and b9 have negative terms. b8 has to be transformed first, otherwise there will be no +1 to cancel out the −1 at the twelfth row of b8. Each −1 in b8 is canceled out through the operation b8 + b9. In the second step, one has to use b7 to cancel −1 at the ninth row of b9. In this way, we obtain the set of basis vectors P. The above analysis indicates that the dimension of this null space is the same as the number of edges for its corresponding convex cone [28], which is the algebraic basis for extreme pathways [14] and elementary modes [29]. The ten basis vectors of P represent ten underlying pathways of the hypoxia network. They are illustrated in Figures 5 and 6, from which one finds two distinct patterns: p1, p7, and p9 belong to the HIFα degradation pathways (Figure 5) and the others belong to the simple association–dissociation pathways (Figure 6). More specifically, through p1, HIFα is directly degraded by reaction k2, a presumably oxygen-independent degradation pathway; whereas in oxygen-dependent pathways p7 and p9, the hydroxylated HIFα is recognized by the VHL that channels it through a ubiquitin degradation component that is shown as the dotted box in Figure 5. Even though p1, p7, and p9 are all elementary modes [29], they can share certain reactions of the network. For example, the total influx for HIFα synthesis from a precursor can thus be decomposed into three parts with the overall rate constant k1 being given by γ1k1, γ2k1, and γ3k1, where γ1 + γ2 + γ3 = 1. The pathways p7, and p9 are almost the same. The only difference is that HIFα is associated with ARNT in the middle part of the pathway p9. Therefore, ARNT must be functionally very important, otherwise it would be hard to explain why the two underlying pathways, which should play significantly different roles, look so similar. Indeed, HIFα degrades differently through the two pathways. The k-sets were selected in [1] on the basis that they produced a switch-like behavior. For all the three k-sets, it was observed that HIFα has high affinity to PHD and low affinity to ARNT. The former is consistent with the usual case of high enzyme-substrate binding affinity. This implies that p9 is not the major degradation pathway because HIFα does not bind with ARNT very well. Moreover, p9 is immediately adjacent to HRE, which suggests its major role is to deliver signals to activate the promoters of hypoxia-regulated genes. As a signal transducer, the rate constant γ3k1 itself need not to be high; what the downstream genes are sensitive to is d (γ3k1) /dt. Therefore, we hypothesize that there is a negligible flux through p9 (or γ3 ≈ 0). To verify our hypothesis, we calculate the γ1, γ2, and γ3 values, as the indication of the relative importance of p1, p7, and p9 in HIFα degradation, at different oxygen levels. The results are given in Table 3. Note that [O2] = 0. 1 and [O2] = 1. 0 represent typical low and high oxygen levels according to [1]. One sees that the pathway p9 is always much less important than the other two as far as HIFα degradation is concerned. The majority of HIFα gets degraded via either p1 at low oxygen or p7 at high oxygen. The comparison of hypoxia response network and heat shock response network [30] as in Table 4 shows the similarity between these two networks with respect to the issue of affinity. The huge difference in the affinity can clearly separate the underlying pathways and assign different functions to them. This is also the basis for the Goldbeter-Koshland model [31]. We tested our method to all three parameter sets (k-sets 1,2, and 3) in [1], and find that the analytical results are almost identical with those of the direct simulations of the entire network, which strongly validates our approximation. For the rest of the paper thereafter, we only report numerical results for k-set 1. The EPA method gives us a starting point from which to analyze our reaction network in greater and more revealing detail. The verification of our hypothesis implies that the pathways associated with p9 could be neglected in the first place. The following equations describe the combination p1, p4, p6, and p7, which constitute the oxygen sensing mechanism: Note the differences between Equations 1,6, and 8, and Equations 14,15, and 17. The p9 related elements have been omitted due to their smallness. By setting the left-hand sides of the above equations to zero, one obtains the steady state equations: The total amount of PHD is conserved: PHD is either in the form of PHD (S12) or HIFα: PHD (S13). This implies where is the initial concentration of PHD. By some derivation, the following equation is obtained: where a = k2k12, b = b1 + b2[O2], c = c1 +c2[O2], where Since c < 0, Equation 25 has one and only one reasonable root Note that none of the species and reactions in the degradation box is present in Equation 26, which indicates that the components in the degradation box are not responsible for the sharp response curve. Once is determined, the analysis of the remaining network (p2, p3, p5, p6, p8, p9, and p10) becomes straightforward, and the results are given in the section “Additional results. ” In fact, these results can be further simplified (see Equations 30 and 43 for and). Figure 7A and 7B shows the steady-state values of [HIFα] () and [HIFα: ARNT: HRE] () at different oxygen values. The red lines depict the simulation results obtained by the numerical integration of the ODE system (1–13) until the steady state is reached. The black lines depict the analytical solutions that are obtained by the algebraic Equations 26 and 56. To better determine the critical point of pathway switching, we calculate ∂ /∂[O2] and ∂ /∂[O2]. The results are shown in Figure 7C and 7D. One sees that both ∂ /∂[O2] and ∂ /∂[O2] change abruptly in a very narrow region of [O2], with the rest of the values almost zero. One observes that the critical point is about [O2]c = 0. 65. We thus show that the sharpness of the response curve can be determined analytically, instead of exhaustively enumerating [O2] values combined with time-consuming numerical integration of a large number of ODEs at each [O2] value.
EPA is a powerful, yet simple tool that can significantly reduce the complexity of the original network and thus make further analytical effort feasible. In this paper, the additional analysis explains precisely the sharp reaction to oxygen of the network as a whole. The clear separation of p7 and p9 indicates their different functions: the pathway p7 and its other associated pathways constitute the sensing of ambient molecular oxygen; in contrast, the pathway p9 and its associated other pathways are responsible for the signal transduction to form the promoters of hypoxia-regulated genes. Most importantly, the simplification allows for a complete explanation of the switch behavior and a clear presentation of the relations between and [O2], and [O2], and. The first step below explains the sharp HIFα stabilization. The second step explains the sharp HRE occupancy that is induced by the HIFα stabilization. HIFα stabilization. This involves the dissociation of pathways p1, p4, p6, and p7 from the whole network, due to the fact that the flux through p9 is always small and can be neglected. It reveals a critical value that corresponds to the switching between pathways p1 and p7. Since an abrupt change often relates to the notion of singularity in mathematics, we proceed to see if a singularity can be found. Under nomoxia, ≈ 0, and can be neglected in Equation 25, which yields One immediately finds the singularity For k-set 1 in [1], one obtains [O2]c = 0. 64, which is exactly the critical value found in Figure 7. When the oxygen level decreases to a value close to [O2]c, the denominator in Equation 27 becomes very small and a can no longer be ignored. Moreover, the term c in Equation 25 can be ignored compared with the large value of. One thus has This explains the linear decrease of versus [O2] increasing in Figure 7A. In summary, one has One can check that k2 can be ignored in the upper branch of Equation 30 due to its smallness, which again demonstrates that the pathway p1 is not important under nomoxia. The very smallness of k2 reflects the importance of p1 under hypoxia, for k2 exists at the denominator of the lower branch of Equation 30. HRE occupancy. The remaining pathways reveal how HIFα stabilization triggers a sharp increase of HIFα: ARNT: HRE, namely the sigmoid curve of versus [O2]. We conclude that the magnitude of HIFα is crucial for the sharpness of the curve. To show this, we need Equations 47,48,53, and 54. By removing the terms that are negligible, these equations turn into Equation 31 holds because k15 and k16 are far less than k3 and k4. Equation 33 holds because is far less than and. Equation 34 holds because, , and are far less than, , and. The validity of the above approximations can be easily checked. For example, for [O2] = 0. 1, one finds k3 = k4 = 1. 66, k15 = k16 = 0. 005, = 1. 19, = 2. 69, = 0. 11, = 0. 89, = 0. 23, = 0. 0016, and = 0. 0005. From Equations 31–34, one obtains where α = α1 + α2/, α1 = + + k6/k5, and α2 = k4k6/ (k3k5). has one and only one reasonable solution where x = 2 /α. Taylor expanding, yields No matter what value is, x < 2 /α1 = 0. 74, so 1 − x2/2 > 0. 73 and x4/8 < 0. 0375. Therefore ≈1 − x2/2 and Under nomoxia is small, so α ≈ α2/ and 1/α ≈ 0. From Equation 37 one has which is also small. Under hypoxia, is large, and thus By substituting Equation 39 into Equation 37, one obtains the important relation between and: where and By substituting Equation 29 into Equation 40, one obtains where m = /β1 and λ = β2b2/a. It is well-known that Equation 41 represents a sigmoid curve with m controlling the saturation value and λ controlling the sharpness. By ignoring the small term k6/k5 in the expression of α1, one finds m is a function of and only. One also finds Here The association (disassociation) constants k3, k5 (k4, k6) exist in the numerator (denominator) of the term k3, k5/ (k4, k6), which implies that the higher the affinity, the sharper the response. The third term b2/a is proportional to the HIF level. In summary, our analysis yields One sees that HIFα is the key to triggering the HIFα: ARNT: HRE response. As long as the oxygen level is greater than [O2]c, HIFα is efficiently degraded by the pathway p7 and maintains a very low level, and the HIFα: ARNT: HRE level is also low (Equation 38). When the oxygen level drops below [O2]c, the system switches to the pathway p1, and HIFα stabilizes with a large concentration (b2/a large). This triggers the sharp increase of HIFα: ARNT: HRE. The smaller k2 is, the larger b2/a, and the sharper the HRE occupancy response (see Figure 8). Also, the validity of our analytical approximation is justified by the close resemblance of Figure 8B and 8C. The three major results of Kohn [1] involve HRE occupancy as a function of the oxygen concentration. The dependence of the curve on ARNT, VHL, and PHD are obtained by both simulation (Figure 9A) and analysis (Figure 9B) and are explained as follows. ARNT dependence. We need only to analyze the pathway p9. The amount of ARNT does not affect the shape of the response curve or the location of the sharp transition, because p9 is the pathway for HRE expression, while p7 is responsible for the sharpness. HIFα: ARNT would not be generated without ARNT and there would not be expressions of HRE for any level of oxygen. At high oxygen levels, the concentrations are similar because HRE occupancy is low anyway. At low oxygen levels, low levels of ARNT will give low HIFα: ARNT and then low HRE occupancy. VHL dependence. VHL is present in both p7 and p9. At high oxygen levels, one should analyze p7 because it is the major pathway for HIFα degradation. The VHL source will affect the upstream HIFα. If VHL concentration is low, it cannot degrade HIFα fast, and the system yields high HRE occupancy. At low oxygen levels, p1 is the major pathway for HIFα degradation, which does not depend on VHL. PHD dependence. One interesting property relates to the different locations of the transition. Using the criterion identified by the alternative model without reaction k2, we can calculate the transition locations as shown in Figure 9B3 for different PHD values. The results are the same as in Figure 9A3. As a matter of fact, considering the p9 pathway only yields a much simpler analytical solution that is also accurate. This further simplification is due to the fact that the expression of HIFαOH: ARNT: HRE is always negligible. The present model of the hypoxia response network is probably an oversimplified one. Nevertheless, it serves as an important starting point, from both theoretical and experimental perspectives, before a more detailed model can be understood. The present model will be gradually expanded and analyzed, with the input of more quantitative data from future experiments. One advantage of EPA is that the method can easily incorporate mechanistic details as soon as they become available [16]. There are various molecular interactions that can be added to the model. For example, it was demonstrated that HIF influences mitochondrial function by inducing pyruvate dehydrogenase kinase 1 (PDK1) to suppress the tricarboxylic acid (TCA) cycle and thus the aerobic respiration. Then the respiration shifts to be anaerobic, whereby the oxygen resource can be preserved to promote cell survival under hypoxic environment [32–34]. Another subject we are interested in is the inclusion of ROS in the network. It has been established that ROS affects HIFα degradation through Fe2+ [7–10]. The direct hydroxylation of HIFα by oxygen requires Fe2+. Under hypoxia conditions, however, ROS increases dramatically and consistently removes Fe2+ via oxidation to Fe3+. Together with the shortage of oxygen, this makes the HIFα degradation through the VHL pathway even more difficult. Consequently, the transition would be faster and sharper. Analysis should focus on the explanation of the coexistence of two oxygen sensing components, a matter that does not appear to be settled as yet. To obtain the dynamical response when the oxgen changes continuously in time, Equations 1–13 (the full model) are integrated. Figure 10 shows the temporal changes of [HIFα] and [HIFα: ARNT: HRE] as responses to the oxygen decreasing from 1. 0 to 0. 1 with different rates. One sees that [S3] and [S7] increase prominently only after [O2] decreases below [O2]c. The faster [O2] decreases, the more rapid the responses are. In particular, when [O2] abruptly jumps from 1. 0 to 0. 1, the responses ensue promptly. However, it is worth noting that one cannot tell practically how fast the responses are because the time scale is unknown. Indeed, the model is dimensionless and no units are given. Nevertheless, the sharp curves illustrate that the responses are very sensitive to the oxygen concentration and imply that the system can provide a timely response under hypoxia. Physiologists have long been puzzled by the ceaseless HIFα cycle, characterized by both abundant generation and efficient degradation, which seems to be a highly wasteful process. Our analysis provides a reasonable explanation. To deal with a sudden environmental change from nomoxia to hypoxia, an organism must respond in time to trigger the genes necessary for adapting to the new environment. To achieve such a sharp response, a high HIFα generation potential is necessary. Since the hypoxia conditions are rare, an efficient degradation pathway has to be designed to maintain a low HIFα under nomoxia. The HIFα cycle is indeed uneconomic, but it appears useful in helping the cell respond to sudden, unpredictable changes in its environment. In summary, we have obtained an accurate analytical solution to the hypoxia response network and have provided a complete explanation of the switch-like behavior first observed and modeled in [1]. The first step of our analysis applied the EPA technique to a reduced, yet complete system that resulted in exposing ten independent pathways, allowing us to focus on analyzing the pathways relevant to HRE occupancy. The analysis showed that the sharp response of HRE occupancy is due to the switch between the pathway p7 (p1) that degrades HIFα under nomoxia (hypoxia).
The remaining network can be solved analytically with HIFα () already determined. The ODE description is as follows: Together with the three constraints of, , and, the steady state equations are expressed as follows: where and have already been determined from the analysis of p7. By some reasoning, a quartic equation is derived, from which can be obtained. and can then be expressed as functions of: where
Protein accession numbers as listed in Table 1 are from http: //www. ncbi. nlm. nih. gov/entrez/query. fcgi? DB=protein.
A complex biomolecular network utilizes different pathways to perform different functions. However, the interactions within the network are typically so complicated that the pathway structure is usually hidden. By some mathematical techniques, the pathways can be identified and possibly decoupled, whereby the insightful details of the network can be exposed. As an example, we study in this paper the hypoxia response network that manifests a dramatic switch-like behavior for certain sets of rate constants: a slight change of the oxygen concentration close to a critical value will lead to distinct reaction patterns. By a technique called extreme pathway analysis, the network is decoupled into three major and some minor pathways. Flux distribution among these pathways can thus be measured by integrating the ordinary differential equations for any given set of rate constants. For the sets of rate constants where the switch-like behavior is observed, we found that such a behavior is due to the switching of flux between two of the three major pathways.
computational biology
biophysics
animals
Pathway Switching Explains the Sharp Response Characteristic of Hypoxia Response Network