Genome-Scale Methods Converge on Key Mitochondrial Genes for the Survival of Human Cardiomyocytes in HypoxiaClinical Perspective
Background—Any reduction in myocardial oxygen delivery relative to its demands can impair cardiac contractile performance. Understanding the mitochondrial metabolic response to hypoxia is key to understanding ischemia tolerance in the myocardium. We used a novel combination of 2 genome-scale methods to study key processes underlying human myocardial hypoxia tolerance. In particular, we hypothesized that computational modeling and evolution would identify similar genes as critical to human myocardial hypoxia tolerance.
Methods and Results—We analyzed a reconstruction of the cardiac mitochondrial metabolic network using constraint-based methods, under conditions of simulated hypoxia. We used flux balance analysis, random sampling, and principal component analysis to explore feasible steady-state solutions. Hypoxia blunted maximal ATP (−17%) and heme (−75%) synthesis and shrank the feasible solution space. Tricarboxylic acid and urea cycle fluxes were also reduced in hypoxia, but phospholipid synthesis was increased. Using mathematical optimization methods, we identified reactions that would be critical to hypoxia tolerance in the human heart. We used data regarding single-nucleotide polymorphism frequency and distribution in the genomes of Tibetans (whose ancestors have resided in persistent high-altitude hypoxia for several millennia). Six reactions were identified by both methods as being critical to mitochondrial ATP production in hypoxia: phosphofructokinase, phosphoglucokinase, complex II, complex IV, aconitase, and fumarase.
Conclusions—Mathematical optimization and evolution converged on similar genes as critical to human myocardial hypoxia tolerance. Our approach is unique and completely novel and demonstrates that genome-scale modeling and genomics can be used in tandem to provide new insights into cardiovascular genetics.
Systems biology uses mathematical and computational methods to describe and explore complex biological networks. An important recent trend in systems biology has been the development and application of constraint-based modeling.1 Constraint-based modeling provides 3 significant advantages compared with traditional mathematical approaches for the study of large and complex biochemical systems. First, large models (up to thousands of reactions) can be accommodated. Thus, the entire metabolic network of a mitochondrion, for example, can be modeled. Second, precise descriptions of the behavior of each enzyme in the system (ie, rate laws) are not required. Finally, detailed information regarding the activity of a single protein (eg, whether an enzyme is allosterically modified or not) is not necessary. Thus, unlike traditional kinetic models, constraint-based models do not rely on, nor do they require, detailed knowledge of an enzyme’s phosphorylation status, for example, nor the abundance of substrates and products.
Clinical Perspective on p 415
Constraint-based modeling is able to confer these advantages because the underlying models and assumptions are simple. The basic unit for constraint-based modeling is a network model, similar to the London Underground map. In the case of a biochemical network, this is constructed using (1) the known presence or absence of reactions based on genomic, proteomic, or biochemical data and (2) the known (species-specific) stoichiometry of all the chemical reactions included in the network. To this basic model are added a series of constraints (from which the method derives its name), including reaction directionality, mass and charge balancing, and absolute limits to metabolite uptake and excretion. Unlike traditional enzyme kinetic parameters (eg, Michaelis–Menten midpoints), many of the underlying assumptions in constraint-based models are robust to variations in physical environment (such as temperature). However, constraint-based models are not able to simulate the exact behavior of a biochemical system. Instead, by keeping the underlying assumptions as simple and robust as possible, they attempt to mirror the constraints which the true network faces in vivo. Nevertheless, one can predict the most likely behavior of the system (using Monte Carlo methods) or predict the behavior of the network at an optimum value of some assumed physiological objective. A more detailed description of the approach can be found in the Methods section and in the Data Supplement and in many excellent reviews.1–4 Regarding its use: constraint-based modeling, using genome-scale metabolic networks, has been used to successfully predict the metabolic signatures of human inherited diseases5–8 and to permit the in silico design of tumor-specific toxins9 and aid in the design of microbial strains for the purposes of metabolic engineering.10
Myocardial ischemia and hypoxia, whether cause or consequence, are common features of the failing heart; understanding the mitochondrial response to hypoxia is key to understanding ischemia tolerance. Myocardial hypoxia can be due to any number of factors, but is most commonly caused by coronary artery or microvascular heart disease, exacerbated by increased oxygen demand from ventricular remodeling. Ischemic heart disease remains the leading cause of death in the developed world; therefore, gaining new insights into the mechanisms whereby heart cells can survive hypoxia of any duration is a matter of considerable importance.
Hypoxia, consequent upon a reduction in barometric pressure, is also a consistent environmental challenge for human populations at high altitude, where it has led to a robustly detectable degree of genetic and phenotypic divergence over evolutionary timescales.11,12 Thus, human populations at high altitude offer a unique opportunity to study the genetic response to hypoxia. We used a novel combination of genome-scale modeling, mathematical optimization, and genome-wide analysis of single-nucleotide polymorphisms (SNPs) in humans to study the response of cardiac mitochondria to hypoxia. In particular, we sought to test our hypothesis that if evolution is an optimization process, then mathematical optimization methods, when applied to a metabolic model, would converge on the same set of reactions, critical to environmental (in this case hypoxic) performance. By comparing information from natural (evolution) and mathematical optimization methods, we sought to identify key genes and reactions that underlie cardiac tolerance to hypoxia.
The reconstruction of the human cardiac mitochondrial metabolic network from proteomic and biochemical data was described previously.13 Briefly, proteomic and transcriptomic data were used to derive an organelle metabolic parts list (ie, a list of all metabolic proteins known to be associated with a cardiac mitochondrion). These parts were connected by their species-specific stoichiometric chemical equations. The draft reconstruction was extensively tested and manually curated. The final model used herein comprised 195 reactions, 235 metabolites, and 25 exchange reactions (for a full description see the articles by Vo et al13 and Thiele et al14 and the Data Supplement). The exchange reactions did not represent genuine biochemical reactions but instead described the exchanges that were necessary between the network and its environment so that a steady-state could be achieved. Having reconstructed the network, a series of limits or constraints were added, all of which constrained the upper and lower limits of metabolite exchange of the model with its environment (eg, oxygen and glucose). These constraints are in Table I in the Data Supplement and represent maximum and minimum flux rates in human heart mitochondria in vivo. To simulate hypoxia, we reduced the upper constraint on oxygen uptake in the model to 25% of baseline values (normoxia), from 39.1 μM min−1 (g mitochondrial protein)−1 (henceforth shortened to U) to 9.775 U. Our choice of simulating severe hypoxia was motivated by an intention to highlight any effects; however, it is worth noting that complete anoxia can occur in regions of ischemic myocardium (eg, during acute myocardial infarction).
Computational analysis of network models rarely leads to a single set of predicted fluxes. Instead, methods are used to analyze the possible combinations of fluxes that allow a steady state, given the applied constraints. The solutions together are termed the feasible steady-state solution space. Alternatively, one can use linear optimization to compute a set of fluxes that optimize the value of a given objective function, an approach typically referred to as flux balance analysis (FBA). For example, this method would return a set of fluxes that correspond to the highest possible rate of ATP production by the network, if ATP production was the objective function. When conducting FBA, we optimized the mitochondrial network for 3 objective functions (phospholipid, heme, and ATP synthesis).13 Alternate optimal solutions (ie, other sets of fluxes that also gave an optimal objective) were accounted for via flux variability analysis (see below). We also studied the optimization of all 3 objectives simultaneously (see Data Supplement for details). This method comprises placing the objective functions under study into hierarchical order (eg, heme biosynthesis then phospholipid biosynthesis then ATP synthesis). The network is optimized for the first objective and then optimized for the second with the first held at optimal value and so forth.
We used 2 computational methods to identify reactions that are critical to hypoxia tolerance in the mitochondria metabolic network: shadow prices15 and flux spans.16 Shadow prices have been used in metabolic network analysis before.15,17 Shadow prices (also known as Lagrange multipliers) are measures of the degree to which the value of the objective function is affected by the availability of a particular resource. For example, if ATP synthesis were the objective function, a shadow price of 1.0 for glucose, for example, would indicate that a 1-unit increase in glucose availability would lead to an equivalent increase in ATP synthesis. A shadow price of 2.0 would indicate that a unit increase in the availability of glucose would result in a 2-unit increase in ATP synthesis, and so forth. We reasoned, therefore, that reactions for which metabolites with large, positive shadow prices were either substrates or products would be crucial to hypoxic performance (at least, for the objective function under investigation). To assess the likelihood that our method had outperformed chance, we used simple permutation testing.
Our second approach was to use flux spans. Using flux variability analysis,18 we computed the range of values that flux through each reaction could take at an optimum (computed using FBA). Taken together, these ranges delineate the set of alternate optimal solutions (ie, different sets of fluxes that result in the same optimal value of the objective).18 By calculating the difference between the upper and lower feasible fluxes, we derived a flux span for each reaction. Here, we express these as a relative ratio. Hence, a reaction with flux=10 U and with the lower and upper feasible fluxes being 8 and 12 U, respectively, would have a relative flux span of 0.4 (or 40%). We reasoned that reactions with the smallest relative flux spans would be critical to hypoxia tolerance and hypoxic performance. Again, we used permutations to estimate the probability that our method had outperformed chance alone.
We used data from a genome-wide allelic differentiation scan comparing SNP frequencies of Tibetans (n=35) residing at 3200 to 5000 m, with 84 individuals from the founder population. Subjects were recruited from 3 distinct regions of China: the North Western region of Yunnan province and Mag Xiang and Zhaxizong Xiang (both in the Tibet Autonomous Region). Genotypic data from the HapMap Phase III Han population were also included. These data have been analyzed previously,11 and full details can be found in this earlier publication. Each gene was assigned a genome-wide P value that serves as an estimate of the degree of selective pressure applied to that gene (through differences in SNP frequency). Details regarding the calculation of these P values can also be found elsewhere.11 We extracted the P values corresponding to the genes in our model and ranked genes by smallest genome-wide allelic differentiation scan P value first, producing a list of nuclear-encoded mitochondrial genes with an accompanying measure of selective pressure in humans living in persistent hypoxia.
Where appropriate, means and SDs are given. However, modeling results are often a single datum point (eg, differences in optimal ATP synthesis rate, determined using FBA, under hypoxia and normoxia) and are therefore given as such.
We optimized the network for 3 physiological targets (objective functions)—ATP, heme, and phospholipid biosynthesis13—using FBA.1 Hypoxia reduced the optimal ATP synthesis rate by 13%, from 45.8 to 36.6 U. Figure 1 shows a quantized heatmap of the accompanying differences in flux. There were reductions in flux through many reactions comprising the tricarboxylic acid (TCA) cycle and oxidative phosphorylation. There were also reductions in flux through most reactions comprising fatty acid uptake, transport, activation, and oxidation although some were maintained due to the imposition of a minimum uptake rate (this is a physiological constraint imposed by the ability of fatty acids to diffuse freely across membranes). Glycolytic rates were similar, which was expected as maximal ATP synthesis was the objective. The flux through multiple reactions required for phospholipid biosynthesis were increased, and the demand reaction was activated in hypoxia. To ensure that the degree of simulated hypoxia affected our results quantitatively but not qualitatively, we performed additional FBA experiments at various intermediate oxygen uptake rates. The results are in Figure I in the Data Supplement. Briefly, maximal ATP synthesis was progressively reduced by increasing hypoxia. Consistent with our interpretation, phospholipid biosynthesis was not activated until O2 uptake dropped below a critical level, at which point a sink for fatty acid carbons was required. There was no evidence of qualitative shifts in carbon flux as maximal O2 uptake rate was progressively reduced.
Heme synthesis was blunted by 75% in hypoxia (hypoxia: 0.650 versus normoxia: 2.44 U). The pattern of flux differences between the optimized network in normoxia and hypoxia was similar to that with ATP synthesis as the objective. Flux through reactions comprising the TCA cycle, oxidative phosphorylation, the urea cycle, and heme synthesis itself were suppressed. There were increases in long-chain (C20:4 and C22:6) activation and phospholipid biosynthesis. A heatmap of differences in flux across the network under these conditions (heme biosynthesis as the objective function in hypoxia versus normoxia) is given in Figure II in the Data Supplement. However, when phospholipid biosynthesis was the objective, it was unchanged by oxygen restriction, at 22.8 U.
We then performed multiple objective analyses with 3 different hierarchies of objective functions. 1. ATP≥heme≥phospholipid: In normoxia, and with ATP synthesis fixed at its optimal value of 45.8 U, heme and phospholipid synthesis were eliminated. In hypoxia, with ATP synthesis fixed at its optimal value of 36.6 U, heme synthesis was still eliminated; however, optimized phospholipid biosynthesis was now nonzero, although reduced ≈100-fold at 0.265 U. 2. Heme≥phospholipid≥ATP: In normoxia, with heme biosynthesis at its optimal rate of 2.44 U, both phospholipid and ATP synthesis were abolished. In hypoxia, with heme biosynthesis at 0.650, phospholipid biosynthesis was possible and optimized to 0.867 U; ATP synthesis was abolished. 3. Phospholipid≥heme≥ATP: As maximal phospholipid biosynthesis was unaffected by hypoxia, it was fixed at 0.867 U for both conditions (normoxia/hypoxia). In both normoxia and hypoxia, heme biosynthesis gained optimal values the same as those where it was the only objective function considered (normoxia: 2.44 U versus hypoxia: 0.650 U). In normoxia, ATP biosynthesis was subsequently limited to 8.85 U; in hypoxia, it was reduced far less, to 19.3 U. A summary of all the optima is given in Table II in the Data Supplement.
We next used uniform random sampling,19 a method that characterizes the steady-state solution space without requiring an objective function. Figure 1B shows a quantized heatmap of differences in median flux. Consistent with the FBA results, fluxes through reactions comprising oxidative phosphorylation, the TCA cycle, and fatty acid metabolism were reduced in hypoxia. Without the requirement to maximize ATP synthesis in normoxia that was elsewhere imposed by FBA, glycolytic flux increased in hypoxia. The heatmap shows a reduction in flux through lactate dehydrogenase and the lactate transporter; however, this represents a reversal in flux, from uptake to efflux. Also noteworthy is a reduction in urea cycle flux in hypoxia.
We analyzed the sampled data using principal component analysis, allowing us to visualize patterns of change. We modeled the sampled data together and found that 5 components captured 65% of the total variance. When plotted, the scores on these components revealed that hypoxia substantially reduced the dimensions of the solution space, reducing the flexibility of the metabolic network even though the dimensions of the space were the same. This was especially apparent in principal components 1 and 2 (Figure 2), with principal component 1 being dominated by reactions related to gas exchange, the TCA cycle, and oxidative phosphorylation and principal component 2 being dominated by reactions related to iron transport and heme biosynthesis.
Given that optimal phospholipid biosynthesis was unaffected by oxygen restriction, we continued by studying those metabolites and reactions that limited optimal ATP and heme biosynthesis in the mitochondrial metabolic network under hypoxia. We first computed shadow prices for all metabolites in the model when optimizing the network for either ATP or heme synthesis. Table 1 gives metabolites with the largest positive shadow prices and the corresponding 20 discrete reactions for each objective function. Three classes of metabolite (and reaction) dominated: long-chain (>20C) fatty acid transport, glycolysis, and heme biosynthesis. When optimizing for ATP synthesis, the largest shadow prices were several fold larger than when optimizing for heme synthesis [eg, 4.0 for cytosolic fructose diphosphate, fdp(c), versus 1.2 for cytosolic arachidonic acid, c206(c)]. When optimizing the network with ATP synthesis as the objective function, the metabolites with large positive shadow prices were mainly related to glycolysis, oxidative phosphorylation, and the TCA cycle.
We next computed flux spans using flux variability analysis. The reactions with the smallest relative flux spans when optimizing the network for either heme or ATP synthesis are given in Table 2. As with shadow prices, the magnitude of the parameter (in this case, relative flux span) was several fold different when optimizing ATP versus heme synthesis; in both cases, the difference (larger shadow prices and smaller flux spans) was consistent with ATP synthesis being more tightly restricted by hypoxia. Interestingly, reactions related to oxidative phosphorylation had narrow flux spans regardless of whether ATP or heme synthesis was the objective. In particular, complex IV of the respiratory chain was ranked in the top two (Table 2). Perhaps unsurprisingly, reactions related to protoheme synthesis and iron transport also had small relative flux spans when optimizing for heme synthesis. When optimizing the network for ATP synthesis, reactions from the TCA cycle and glycolysis were highly represented.
We generated a complete list of nuclear-encoded mitochondrial genes with a corresponding measure of selective pressure at high altitude. The 20 most selected genes (ie, smallest P value) are given in Table 3. When conducting shadow price analysis with ATP synthesis as the objective, we selected metabolites with the largest positive shadow prices resulting in the 20 discrete reactions shown in Table 1. Of these 20 reactions, two (phosphofructokinase and phosphoglycerate kinase) carried flux and corresponded to genes in Table 3. However, permutation testing (as detailed in the Data Supplement) suggested that our method performed little better than chance. We observed that 2 of the only 3 metabolites with shadow prices of 3 or greater when optimizing the network for either objective (cytosolic fructose 6-phosphate and fructose diphosphate, f6p(c) and fdp(c), respectively, in Table 1) are both substrate and product for phosphofructokinase, the most heavily selected gene in Table 3. We next compared the reactions highlighted by shadow price analysis while optimizing the network for heme synthesis. Three reactions were common with those in Table 3: hydroxymethylbilane synthase, porphobilinogen synthase, and phosphoglycerate kinase. Once again, permutations suggested that shadow pricing had not outperformed chance when identifying genes under pressure.
We used a similar approach to assess the performance of flux span analysis. Table 2 shows that flux span analysis identified 3 (heme synthesis) and 6 (ATP synthesis) reactions that were common with those that were the most heavily selected genes in Table 3. With ATP synthesis as the objective, we used permutation testing to assess whether modeling had outperformed chance when predicting genes under pressure. Random selections only matched flux span analysis 1675 out of 100 000 times, offering evidence that this approach had outperformed chance at P<0.05.
We repeated this process with heme synthesis as the objective function. Using 100 000 permutations, random selection matched the model performance approximately half the time (P=≈0.525). Hence, using flux span analysis with heme as the objective had not outperformed chance.
Myocardial hypoxia can be either acute or chronic and occurs whenever oxygen delivery is insufficient to meet the needs of the contracting myocardium. This can be due to any combination of reduced O2-carrying capacity due to anemia, reduced hemoglobin saturation (whether environmental or pathological), poor cardiac output or compromised blood flow due to coronary artery or microvascular heart disease, and increased oxygen demand associated with stress or structural remodeling (eg, ventricular hypertrophy). Ischemic heart disease remains the leading cause of death in the developed world. A notable feature of heart failure is that, once left ventricular dysfunction has been established, patients suffer a relentless apoptotic loss of viable cardiomyocytes that some investigators believe to be due to repeated, transient ischemic and hypoxic events.20 Therefore, understanding the mechanisms whereby heart cells can survive either transient or sustained hypoxia and ischemia is a matter of considerable importance. Here, we present an entirely new approach to this question using systems biology methods that encompass genomics, metabolic modeling, and mathematical optimization.
We first studied the effect that hypoxia had on the solution space (the set of all feasible fluxes) of the reconstructed cardiac mitochondrial metabolic network using 2 complementary methods: optimization (FBA) and Monte Carlo sampling. Optimization requires an objective function; in keeping with previous work, we studied 3 objectives that are central to mitochondrial function: the synthesis of ATP, heme, and mixed phospholipids.13 Although maximal phospholipid synthesis was unaffected by hypoxia, both heme and ATP synthesis were reduced (by 75% and 13%, respectively). The reduction in maximal ATP synthesis was accompanied by reductions in TCA cycle flux, oxidative phosphorylation, and fatty acid uptake and processing but a seemingly paradoxical increase in phospholipid biosynthesis.
The degree to which maximal heme synthesis was blunted in hypoxia was striking. This 75% loss of protoheme synthesis capacity was accompanied by many of the metabolic features observed when optimizing ATP production in hypoxia. Heme is a major component of hemoglobin, itself substantially increased in response to hypoxia to enhance systemic oxygen transport.21 Thus, network stoichiometry forms a constraint to heme biosynthesis that may partly define the speed with which heme can be synthesized in hypoxia. Non–iron-deficient anemia is a common feature in patients with heart failure, yet its cause is unknown.22 Furthermore, many studies have shown that reduced hemoglobin is an independent predictor of risk in patients with heart failure23,24 although again the mechanism remains poorly understood. Our simulations suggest that hypoxia itself can cause significant reductions in protoheme synthesis, both in the heart and elsewhere, and that hypoxia of any kind could lead to a vicious cycle of blunted heme synthesis, reduced O2-carrying capacity in blood, and subsequently worsened hypoxemia. Furthermore, in cultured human neurons, heme deficiency causes a decrease in (heme-containing) complex IV (cytochorome c oxidase) and activation of nitric oxide synthase.25 Given that complex IV release is a key component of the p53 apoptotic cascade,26 this suggests an intriguing new avenue for investigations into hypoxia-induced myocyte apoptosis.
We studied the effect on mitochondria of forcing them to balance competing objectives in hypoxia. When ATP production was given hierarchical superiority (as is almost certainly the case in vivo), heme synthesis was completely abolished in normoxia and hypoxia. Thus, heme and ATP synthesis compete for the same resources; moderate reductions in O2 supply, coupled with increases in ATP demand, might lead to profound reductions in heme synthesis capacity due to stoichiometric constraints in the metabolic network. Interestingly, phospholipid biosynthesis was abolished in normoxia under this hierarchy but was active in hypoxia. This was likely due to competition with ATP synthesis for lipids; in hypoxia, ATP synthesis was diminished freeing up lipids for phospholipid synthesis.
We next studied the set of feasible solutions using random sampling. Multiple random samples of the solution space allow the generation of probability density functions for flux through every reaction; the most probable value can often predict the measured in vivo rate.14 As with linear optimization (FBA), we observed reductions in TCA cycle and oxidative phosphorylation reactions, in addition to a reduction in urea cycle flux. This last is intriguing because flux through arginosuccinate synthetase (ASS) was decreased in simulated hypoxia. ASS been elsewhere reported as a target of the von Hippel-Lindau tumor suppressor gene (VHL),27 an inverse regulator of hypoxia-inducible factor-1α. Manipulation of VHL expression led to corresponding changes in ASS levels in RCC10 renal cell carcinoma cells.
Once again, we observed that many of the reactions required for phospholipid biosynthesis were increased and the biosynthesis reaction itself was activated in hypoxia. Without this redirection of fatty acid flux, the imposition of hypoxia, combined with minimum uptake rates for fatty acids, would have led to an accumulation of unoxidized fatty acids in the model and a loss of homeostasis. Cardiac mitochondria face an identical challenge in vivo and redirect fatty acids to storage (away from oxidation) when ischemic.28 Similarly, Langendorff-perfused hearts exposed to acute hypoxia increase phospholipid biosynthesis to maintain lipid homeostasis.29 We also observed an increase in glycolytic flux (Figure 1B).
Overall, the pattern of change in metabolic flux in our simulations was strikingly consistent with experimental observations of cellular responses to hypoxia, including the reduction in flux through mitochondrial pyruvate dehydrogenase in vivo that is brought about by modulation of pyruvate dehydrogenase kinase.30 It is interesting to note that the reduction in flux through mitochondrial pyruvate dehydrogenase in our simulations directly resulted from network stoichiometry, without any additional explanation or control. Although the notion that glycolytic flux is increased in hypoxia is certainly not new, altered (particularly increased) lipid biosynthesis in response to hypoxia is a less often considered component of hypoxia tolerance. Previous investigators have reported both increased29,31 and decreased32 lipid synthesis in hypoxia in model systems. These discrepancies may be due to differences in isotope labeling strategy (eg, acetate versus glycerol versus palmitate) or outcome measure. However, there is no question that lipids accumulate in the heart in response to hypoxia and ischemia.28 It should be noted that the details of whole-heart lipid handling in hypoxia and ischemia may be different to the mitochondrial response considered in isolation. It is also interesting that lipotoxicity—defined as a chronic mismatch between oversupply of acetyl-CoA from lipid breakdown and its subsequent mitochondrial oxidation—is a stoichiometric disorder and can be as readily caused by impaired oxidative phosphorylation (eg, by hypoxia) as lipid oversupply. The consistency of our simulations with experimental observations reinforced to us the notion that our methods were both robust and relevant.
We sought to test whether mathematical optimization had converged on the same reactions that human evolution had identified as being critical to optimal hypoxic function. We used data from a genome-wide allelic differentiation scan comparing SNP frequencies of Tibetans residing at 3500 m (and whose ancestors have lived high for >10 000 years)33 with individuals from the HapMap Phase III Chinese Han sample, who are closely related but have resided at sea level throughout.11 Tibetans were ideal for this study because, despite systemic adjustments (eg, increased breathing rates), they continue to have lower arterial oxygen content than sea-level dwellers.34
It is interesting that the largest shadow prices were recorded when optimizing the mitochondrial metabolic network for ATP synthesis in hypoxia. This suggests that, even in the case of competing objectives, increasing the supply of these metabolites would be especially advantageous when oxygen supply is limited (either by environment or pathology). The metabolite with the largest shadow price in any analysis we conducted was fructose diphosphate, a product of phosphofructokinase. However, permutation testing failed to support the notion that shadow prices and evolution had converged on similar reactions.
The results gained by examining flux spans were more compelling. Flux spans are the range of values within which a reaction rate can lie at a computed optimum. We reasoned that reactions with narrow flux spans would be under greater selective pressure. We generated a list of reactions with the smallest flux spans (yet which carried flux) and compared these with the SNP data. When we optimized for ATP synthesis, the results supported the notion that mathematical optimization and evolution had converged on similar reactions (where 6/20 reactions were common between the 2 selections). The common reactions selected by flux span analysis and evolution were related to heme synthesis (although only when optimizing for heme synthesis), glycolysis (phosphofructokinase and phosphoglycerate kinase), the TCA cycle (aconitase and fumarase), and oxidative phosphorylation (complexes II and IV). We propose that our combined method has identified reactions that are especially important in maintaining or increasing mitochondrial ATP synthesis in the hypoxic heart. This view is supported by the existing literature. For example, it was recently reported that mice exposed to 3 weeks of normobaric hypoxia had reduced complex II, complex IV, and aconitase activity in cardiac mitochondria35 while fumarate accumulation leads to pseudohypoxic activation of hypoxia-inducible factor-1α,36 suggesting that many of the same reactions highlighted here indeed have important roles in hypoxic adaptation and, hence, survival. Our combined approach also yielded an unexpected benefit: Computational analysis was able to provide suggestions as to whether genes were under positive or negative selective pressure (an important distinction to which traditional genome-wide analytic techniques are blind).
A final note regarding phosphofructokinase: basic biochemistry textbooks all highlight the importance of phosphofructokinase as a key regulatory step in glycolysis (eg, page 444 in the book by Berg et al37). Yet, there is a tautology here: phosphofructokinase is heavily regulated biologically (eg, by ATP/AMP and fructose 2,6-bisphosphate37). However, its heavy regulation is evidence for, not an explanation of, its importance. We note that in our simulations, using multiple objectives and alternative analytic strategies, phosphofructokinase was repeatedly highlighted as being an important determinant of the objective. Our model contained no information whatsoever regarding biological regulation (eg, allosteric modulation by other small molecules). In other words, our simulations suggest that phosphofructokinase is important because it occupies a critical point in the metabolic network due to network topology and nothing more. By extension, this protein is likely to be under strong evolutionary selective pressure in many environments, leading to complex phenotypic properties. Once again, this was supported in the genetic data, at least in hypoxia.
Our main hypothesis that evolution and mathematical optimization would converge on similar targets was supported. In doing so, we generated a list of genes that the 2 methods independently highlighted as potentially important for hypoxic survival. Although we think that the nature of our combined approach adds additional support to the significance of these genes, we wish to stress that genes identified by any genome-wide method should be treated as candidates only. Direct experimental evidence will always be required to clarify the function of each. Of course, for some of the genes identified by our approach, overwhelming evidence already exists confirming their importance (eg, pyruvate dehydrogenase).30,38,39
A second limitation relates to possible differences in the Han versus Tibetan environment beyond simply altitude (eg, diet and temperature). Several points are pertinent:
Although temperatures may differ between the 2 locations, most high-altitude populations descend lower in winter.
Diet may differ; however, many essential elements (reliance on vegetables and use of rice) are similar.
Multiple studies have used the Han versus Tibetan genome comparison. All have found the same primary hit (EPAS1), which is a gene regulating expression of a hypoxia-responsive transcription factor.
The candidates in the present study were chosen because computational analysis of a separate network model suggested their role in hypoxia. This makes it more likely that this was indeed the cause and is, potentially, another benefit of our approach.
Sources of Funding
Dr Thiele was supported, in part, by an Attract program grant (FNR/A12/01) from the Luxembourg National Research Fund (Fonds National de la Recherche).
H. Montgomery was, from 2011 to 2013, contracted as a consultant to GlaxoSmithKline relating to development of a drug in the field of hypoxia. However, no involvement was needed and he received no payment.
The Data Supplement is available at http://circgenetics.ahajournals.org/lookup/suppl/doi:10.1161/CIRCGENETICS.113.000269/-/DC1.
- Received July 4, 2013.
- Accepted April 29, 2014.
- © 2014 American Heart Association, Inc.
- Beall CM,
- Cavalleri GL,
- Deng L,
- Elston RC,
- Gao Y,
- Knight J,
- et al
- Yi X,
- Liang Y,
- Huerta-Sanchez E,
- Jin X,
- Cuo ZX,
- Pool JE,
- et al
- Vo TD,
- Greenberg HJ,
- Palsson BO
- Thiele I,
- Price ND,
- Vo TD,
- Palsson BØ
- Varma A,
- Boesch BW,
- Palsson BO
- Martin D,
- Windsor J
- Coats AJ
- Sharma R,
- Francis DP,
- Pitt B,
- Poole-Wilson PA,
- Coats AJ,
- Anker SD
- Go AS,
- Yang J,
- Ackerson LM,
- Lepper K,
- Robbins S,
- Massie BM,
- et al
- Atamna H,
- Killilea DW,
- Killilea AN,
- Ames BN
- Schuler M,
- Bossy-Wetzel E,
- Goldstein JC,
- Fitzgerald P,
- Green DR
- Beall CM
- Ashrafian H,
- O’Flaherty L,
- Adam J,
- Steeples V,
- Chung YL,
- East P,
- et al
- Berg JM,
- Tymoczko JL,
- Stryer L
- Mora A,
- Davies AM,
- Bertrand L,
- Sharif I,
- Budas GR,
- Jovanović S,
- et al
Systems biology offers a host of new tools to study complex biological systems. Instead of studying individual components (eg, a single protein), the approach promotes studying the effects of environmental or pathological perturbations on the system as a whole. An integral factor is the reconstruction of the biological system into mathematical form and the application of computational and statistical methods to estimate the effects of both internal (eg, genetic) and external (environmental) perturbations on the function of the model. Currently, metabolic systems biology is at the forefront of the development of systems biology. This article describes how several systems biology methods were used together to find the critical steps of mitochondrial metabolism involved in cellular tolerance to hypoxia. Both genomic screening and computational analysis of a metabolic network converged on similar reactions as being critical to hypoxia tolerance. Identifying these key reactions can aid the clinician in understanding interindividual variability in response to hypoxia and highlight potential new targets for treatment of patients experiencing acute and more chronic forms of cellular hypoxia. This can hopefully aid in the management of critically ill patients experiencing acute cellular hypoxia such as after cardiac arrest or in various forms of shock where cellular supply of oxygen is decreased. It might also have implications for patients experiencing severe hypoxia from acute lung disease (such as patients with acute respiratory distress syndrome) or chronic hypoxia (such as patients with chronic obstructive or restrictive lung disease).