Pathway Switching Explains the Sharp Response Characteristic of Hypoxia Response Network

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.


Introduction
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 hypoxiainduced 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][2][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][5][6]. HIF is a heterodimeric transcription factor which consists of one of the three different members (HIF-1a, HIF-2a, and HIF-3a) and a common constitutive ARNT subunit which is also known as HIFb. The system also includes an enzyme family: prolyl hydroxylases (PHDs), which directly sense the level of oxygen and hydroxylate HIFa by covalently modifying the HIFa 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, Fe 2þ , for the hydroxylation of HIFa by a PHD [7][8][9][10]. There are three members in this enzyme family: PHD1, PHD2, and PHD3. The hydroxylated HIFa 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 HIFa 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 HIFa 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][15][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][18][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][24][25][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  Author Summary 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.
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 HIFa plays the main role in the sharp response.

Model Description
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 k 32 ). Hence, Transcript intermediate 1, 2, and 3 (Species 8, 9, and 10) can also be dropped, as well as the associated reactions: k 7 , k 8 , k 9 , k 10 , and k 11 , because they do not affect the dynamics of the network. Species 23 is only the joint name of HIFa:ARNT:HRE (Species 7) and HIFaOH:ARNT:HRE (Species 22); therefore, it is dropped. HIFa precursor (Species 1) is a constant and is thus dropped, because its information can be simply encoded in the reaction k 1 . 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 [S n ] stands for the concentration of species n as tabulated in Table 1 and [O 2 ] indicates the input cellular oxygen concentration. Table 2 shows the specific reaction each rate constants k n represents.
The real values of k n 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 [S 1 ] since it is set to unity.

Network Decomposition by Extreme Pathway Analysis
This section assumes some familiarity with linear algebra. The 13 3 19 stoichiometric matrix U of the reduced hypoxia response network (Tables 1 and 2) is shown in Figure 2 (for details, see Materials and Methods). The rank of U 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 b 8 and b 9 have negative terms. b 8 has to be transformed first, otherwise there will be no þ1 to cancel out the À1 at the twelfth row of b 8 . Each À1 in b 8 is canceled out through the operation b 8 þ b 9 . In the second step, one has to use b 7 to cancel À1 at the ninth row of b 9 . 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: p 1 , p 7 , and p 9 belong to the HIFa degradation pathways ( Figure 5) and the others belong to the simple associationdissociation pathways ( Figure 6). More specifically, through p 1 , HIFa is directly degraded by reaction k 2 , a presumably oxygen-independent degradation pathway; whereas in oxygen-dependent pathways p 7 and p 9 , the hydroxylated HIFa 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 p 1 , p 7 , and p 9 are all elementary modes [29], they can share certain reactions of the network. For example, the total influx for HIFa synthesis from a precursor can thus be decomposed into three parts with the overall rate constant k 1 being given by c 1 k 1 , c 2 k 1 , and The pathways p 7 , and p 9 are almost the same. The only difference is that HIFa is associated with ARNT in the middle part of the pathway p 9 . 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, HIFa 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 HIFa has high affinity to PHD and low affinity to ARNT. The former is consistent with the usual case of high enzymesubstrate binding affinity. This implies that p 9 is not the major degradation pathway because HIFa does not bind with ARNT very well. Moreover, p 9 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 c 3 k 1 itself need not to be high; what the downstream genes are sensitive to is d(c 3 k 1 )/dt. Therefore, we hypothesize that there is a negligible flux through p 9 (or c 3 » 0). To verify our hypothesis, we calculate the c 1 , c 2 , and c 3 values, as the indication of the relative importance of p 1 , p 7 , and p 9 in HIFa degradation, at different oxygen levels. The results are given in Table 3. Note that [O 2 ] ¼ 0.1 and [O 2 ] ¼ 1.0 represent typical low and high oxygen levels according to [1]. One sees that the pathway p 9 is always much less important than the other two as far as HIFa degradation is concerned. The majority of HIFa gets degraded via either p 1 at low oxygen or p 7 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.

Analytical Solution to the Decomposed Hypoxia Response Network
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 p 9 could be neglected in the first place. The following equations describe the combination p 1 , p 4 , p 6 , and p 7 , which constitute the oxygen sensing mechanism: Note the differences between Equations 1, 6, and 8, and Equations 14, 15, and 17. The p 9 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 (S 12 ) or HIFa:PHD (S 13 ). This implies where ½S 0 12 is the initial concentration of PHD. By some derivation, the following equation is obtained: 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 ½S Ã 3 is determined, the analysis of the remaining network (p 2 , p 3 , p 5 , p 6 , p 8 , p 9 , and p 10 ) 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 ½S Ã 3 and ½S Ã 7 ). Figure 7A

Discussion
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 p 7 and p 9 indicates their different functions: the pathway p 7 and its other associated pathways constitute the sensing of ambient molecular oxygen; in contrast, the pathway p 9 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 ½S Ã 3 and [O 2 ], ½S Ã 7 and [O 2 ], ½S Ã 3 and ½S Ã 7 . The first step below explains the sharp HIFa stabilization. The second step explains the sharp HRE occupancy that is induced by the HIFa stabilization.
HIFa stabilization. This involves the dissociation of pathways p 1 , p 4 , p 6 , and p 7 from the whole network, due to the fact that the flux through p 9 is always small and can be neglected. It reveals a critical value that corresponds to the switching between pathways p 1 and p 7 . 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, ½S Ã 3 » 0, and a½S Ã 3 2 can be neglected in Equation 25, which yields One immediately finds the singularity For k-set 1 in [1], one obtains [O 2 ] c ¼ 0.64, which is exactly the critical value found in Figure 7. When the oxygen level Through p 1 , HIFa is degraded directly by k 2 . In p 7 , HIFa first binds with PHD after synthesis. HIFa:PHD is then dissociated into PHD and HIFaOH with the participation of oxygen. HIFaOH then binds with VHL to form HIFaOH:VHL. Finally, the dissociation of HIFaOH:VHL concludes HIFa degradation. The pathway p 9 differs from p 7 only in that HIFa first binds with ARNT after synthesization. doi: 10 This explains the linear decrease of ½S Ã 3 versus [O 2 ] increasing in Figure 7A. In summary, one has One can check that k 2 can be ignored in the upper branch of Equation 30 due to its smallness, which again demonstrates that the pathway p 1 is not important under nomoxia. The very smallness of k 2 reflects the importance of p 1 under hypoxia, for k 2 exists at the denominator of the lower branch of Equation 30.
HRE occupancy. The remaining pathways reveal how HIFa stabilization triggers a sharp increase of HIFa:ARNT:HRE, namely the sigmoid curve of ½S Ã 7 versus [O 2 ]. We conclude that the magnitude of HIFa 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 where a ¼ a 1 þ a 2 / ½S Ã 3 , a 1 ¼ ½S 0 4 þ ½S 0 6 þ k 6 /k 5 , and a 2 ¼ k 4 k 6 / (k 3 k 5 ). ½S Ã 7 has one and only one reasonable solution where x ¼ 2 ½S 0 4 ½S 0 6 /a. Taylor expanding No matter what ½S Ã 3 value is, x , 2 ½S 0 4 ½S 0 6 /a 1 ¼ 0.74, so 1 À x 2 /2 . 0.73 and x 4 /8 , 0.0375. Therefore ffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 À x 2 p »1 À x 2 /2 and Under nomoxia ½S Ã 3 is small, so a » a 2 / ½S Ã 3 and 1/a » 0. From Equation 37 one has which is also small. Under hypoxia, ½S Ã 3 is large, and thus By substituting Equation 39 into Equation 37, one obtains the important relation between ½S Ã 3 and ½S Ã 7 :   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 p 9 . The amount of ARNT does not affect the shape of the response curve or the location of the sharp transition, because p 9 is the pathway for HRE expression, while p 7 is responsible for the sharpness. HIFa: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 HIFa:ARNT and then low HRE occupancy.
VHL dependence. VHL is present in both p 7 and p 9 . At high oxygen levels, one should analyze p 7 because it is the major pathway for HIFa degradation. The VHL source will affect the upstream HIFa. If VHL concentration is low, it cannot degrade HIFa fast, and the system yields high HRE occupancy. At low oxygen levels, p 1 is the major pathway for HIFa 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 k 2 , 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 p 9 pathway only yields a much simpler analytical solution that is also accurate. This further simplification is due to the fact that the expression of HIFaOH: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][33][34]. Another subject we are interested in is the inclusion of ROS in the network. It has been established that ROS affects HIFa degradation through Fe 2þ [7][8][9][10]. The direct hydroxylation of HIFa by oxygen requires Fe 2þ . Under hypoxia conditions, however, ROS increases dramatically and consistently removes Fe 2þ via oxidation to Fe 3þ . Together with the shortage of oxygen, this makes the HIFa 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  ] 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 HIFa 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 HIFa generation potential is necessary. Since the hypoxia conditions are rare, an efficient degradation pathway has to be designed to maintain a low HIFa under nomoxia. The HIFa 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 p 7 (p 1 ) that degrades HIFa under nomoxia (hypoxia).

Materials and Methods
In following the law of mass conservation, a particular reactant through each reaction can be written in the form of homogeneous linear equations,