Article Information
PubMed
Related Articles
- Using the Topology of Metabolic Networks to Predict Viabilit...Using the Topology of Metabolic Networks to Predict Viability of Mutant Strains
Biophysical Journal, Volume 91, Issue 6, 15 September 2006, Pages 2304-2311
Zeba Wunderlich and Leonid A. Mirny
AbstractUnderstanding the relationships between the structure (topology) and function of biological networks is a central question of systems biology. The idea that topology is a major determinant of systems function has become an attractive and highly disputed hypothesis. Although structural analysis of interaction networks demonstrates a correlation between the topological properties of a node (protein, gene) in the network and its functional essentiality, the analysis of metabolic networks fails to find such correlations. In contrast, approaches utilizing both the topology and biochemical parameters of metabolic networks, e.g., flux balance analysis, are more successful in predicting phenotypes of knockout strains. We reconcile these seemingly conflicting results by showing that the topology of the metabolic networks of both and are, in fact, sufficient to predict the viability of knockout strains with accuracy comparable to flux balance analysis on large, unbiased mutant data sets. This surprising result is obtained by introducing a novel topology-based measure of network transport: synthetic accessibility. We also show that other popular topology-based characteristics such as node degree, graph diameter, and node usage (betweenness) fail to predict the viability of mutant strains. The success of synthetic accessibility demonstrates its ability to capture the essential properties of the metabolic network, such as the branching of chemical reactions and the directed transport of material from inputs to outputs. Our results strongly support a link between the topology and function of biological networks and, in agreement with recent genetic studies, emphasize the minimal role of flux rerouting in providing robustness of mutant strains.
Abstract | Full Text | PDF (208 kb) - Thermodynamics-Based Metabolic Flux AnalysisThermodynamics-Based Metabolic Flux Analysis
Biophysical Journal, Volume 92, Issue 5, 1 March 2007, Pages 1792-1805
Christopher S. Henry, Linda J. Broadbelt and Vassily Hatzimanikatis
AbstractA new form of metabolic flux analysis (MFA) called thermodynamics-based metabolic flux analysis (TMFA) is introduced with the capability of generating thermodynamically feasible flux and metabolite activity profiles on a genome scale. TMFA involves the use of a set of linear thermodynamic constraints in addition to the mass balance constraints typically used in MFA. TMFA produces flux distributions that do not contain any thermodynamically infeasible reactions or pathways, and it provides information about the free energy change of reactions and the range of metabolite activities in addition to reaction fluxes. TMFA is applied to study the thermodynamically feasible ranges for the fluxes and the Gibbs free energy change, Δ′, of the reactions and the activities of the metabolites in the genome-scale metabolic model of developed by Palsson and co-workers. In the TMFA of the genome scale model, the metabolite activities and reaction Δ′ are able to achieve a wide range of values at optimal growth. The reaction dihydroorotase is identified as a possible thermodynamic bottleneck in metabolism with a Δ′ constrained close to zero while numerous reactions are identified throughout metabolism for which Δ′ is always highly negative regardless of metabolite concentrations. As it has been proposed previously, these reactions with exclusively negative Δ′ might be candidates for cell regulation, and we find that a significant number of these reactions appear to be the first steps in the linear portion of numerous biosynthesis pathways. The thermodynamically feasible ranges for the concentration ratios ATP/ADP, NAD(P)/NAD(P)H, and are also determined and found to encompass the values observed experimentally in every case. Further, we find that the NAD/NADH and NADP/NADPH ratios maintained in the cell are close to the minimum feasible ratio and maximum feasible ratio, respectively.
Abstract | Full Text | PDF (676 kb) - Levels of thermodynamic treatment of biochemical reaction sy...Levels of thermodynamic treatment of biochemical reaction systems
Biophysical Journal, Volume 65, Issue 3, 1 September 1993, Pages 1243-1254
R.A. Alberty
AbstractEquilibrium calculations on biochemical reaction systems can be made at three levels. Level 1 is the usual chemical calculation with species at specified temperature and pressure using standard Gibbs energies of formation of species or equilibrium constants K. Level 2 utilizes reactants such as ATP (a sum of species) at specified T, P, pH, and pMg with standard transformed Gibbs energies of formation of reactants or apparent equilibrium constants K'. Calculations at this level can also be made on the enzymatic mechanism for a biochemical reaction. Level 3 utilizes reactants at specified T, P, pH, and pMg, but the equilibrium concentrations of certain reactants are also specified. The fundamental equation of thermodynamics is derived here for Level 3. Equilibrium calculations at this level use standard transformed Gibbs energies of formation of reactants at specified concentrations of certain reactants or apparent equilibrium constants K". Level 3 is useful in calculating equilibrium concentrations of reactants that can be reached in a living cell when some of the reactants are available at steady-state concentrations. Calculations at all three levels are facilitated by the use of conservation matrices and stoichiometric number matrices for systems. Three cases involving glucokinase, glucose-6-phosphatase, and ATPase are discussed.
Abstract | PDF (1047 kb) - …more
Copyright © 2006 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 90, Issue 4, 1453-1461, 15 February 2006
doi:10.1529/biophysj.105.071720
Bioenergetics
Genome-Scale Thermodynamic Analysis of Escherichia coli Metabolism
Christopher S. Henry1, Matthew D. Jankowski1, Linda J. Broadbelt and Vassily Hatzimanikatis
, 
Address reprint requests to V. Hatzimanikatis, Tel.: 847-491-5357.1 C. S. Henry and M. D. Jankowski contributed equally to this work.
Abstract
Genome-scale metabolic models are an invaluable tool for analyzing metabolic systems as they provide a more complete picture of the processes of metabolism. We have constructed a genome-scale metabolic model of Escherichia coli based on the iJR904 model developed by the Palsson Laboratory at the University of California at San Diego. Group contribution methods were utilized to estimate the standard Gibbs free energy change of every reaction in the constructed model. Reactions in the model were classified based on the activity of the reactions during optimal growth on glucose in aerobic media. The most thermodynamically unfavorable reactions involved in the production of biomass in E. coli were identified as ATP phosphoribosyltransferase, ATP synthase, methylene-tetra-hydrofolate dehydrogenase, and tryptophanase. The effect of a knockout of these reactions on the production of biomass and the production of individual biomass precursors was analyzed. Changes in the distribution of fluxes in the cell after knockout of these unfavorable reactions were also studied. The methodologies and results discussed can be used to facilitate the refinement of the feasible ranges for cellular parameters such as species concentrations and reaction rate constants.Introduction
Thermodynamic analysis of reaction systems provides a means of characterizing and describing the equilibrium state of the reactions in the system. Metabolic pathways are open systems, and they cannot exist in a state of thermodynamic equilibrium. However, thermodynamic analysis is invaluable in establishing the limits of activity of metabolic systems, and these limits are important for constraints-based modeling 1,2 and for understanding the design and evolution of metabolism.
The most prevalent constraints-based modeling technique, flux balance analysis (FBA), is based on the quasi-steady-state assumption that the net accumulation of every metabolite in a cell is zero 3, and the mass balance equations of each metabolite are used to formulate a set of linear constraints. Metabolic systems typically involve far more reactions than metabolites making these systems underdetermined, and as a result, these mass balance constraints are insufficient to uniquely determine the flux through all of the reactions in the metabolic network. Despite this limitation, FBA can be used to test the feasibility of possible flux distributions, and it has been utilized extensively to interpret NMR data for estimating intracellular fluxes 4, to provide a reference state for metabolic control analysis 5, to analyze metabolite production and growth rates in cell cultures 6,7, and to predict the effect of gene knockouts 8,9,10.
One method for improving the quality and accuracy of flux quantification through FBA is to provide tighter constraints on the flux solution space 1,11. Thermodynamic analysis provides a means of accomplishing this goal. Currently thermodynamic analysis has found only limited application in the study of metabolic networks. Beard and Qian have conducted studies on the topic of eliminating internal flux cycles 2,12,13. These are sets of reactions for which the overall reaction is zero, such as A→B→C→A. According to the first law of thermodynamics, the overall thermodynamic driving force through this cycle must be zero, meaning no net-flux is possible through this cycle. These cycles are often referred to as type-3 extreme pathways 14. Through the introduction of the appropriate constraints, flux distributions from FBA will no longer involve any flux through type-3 extreme pathways. This analysis only requires that the stoichiometry of the system be known, but no quantitative information on the relative thermodynamic feasibility of the individual reactions and pathways in the metabolic chemistry is provided. Using the limited amount of experimental thermodynamic data currently available, Beard and colleagues also performed a study on the central carbon pathways of the hepatocyte cell, and they quantified the levels of metabolite concentrations and reaction fluxes using thermodynamic constraints 15.
Here we have applied thermodynamic analysis to study Escherichia coli metabolism described by the iJR904 genome-scale metabolic model of E. coli16. We employed the group contribution method of Mavrovouniotis 17,18 to estimate the thermodynamic feasibility of the reactions in E. coli metabolism. We utilized FBA to determine the thermodynamically unfavorable reactions that are essential for optimal growth yield, and we performed knockout studies of these reactions to determine the role these reactions play in cell growth and in the production of individual biomass precursors. We also studied the shift in the flux distribution when the activity of a thermodynamically unfavorable reaction was removed. The Methods and Results presented in this article are directly applicable to improving predictions of the effects of gene knockouts, refining the estimation of cellular parameters such as species concentrations or reaction rate constants, and analyzing a proposed pathway for thermodynamic infeasibilities.
Methods
Definition of ΔrG′m for the assessment of thermodynamic feasibility
The most common measure used for assessing the thermodynamic feasibility of reactions is the Gibbs free energy change of reaction, 
which can be calculated using Eq. (1),
 | (1) |
is the standard Gibbs free energy of formation of compound i, R is the universal gas constant, T is the temperature assumed to be 298K, m is the number of compounds involved in the reaction, xi is the activity of compound i, and ni is the stoichiometric coefficient of compound i in the reaction (ni is negative for reactants and positive for products). Although the activities of most compounds in biological systems are unknown, the mean activity in the cell is on the order of 1mM 19. Therefore, using 
for the assessment of the thermodynamic feasibility of metabolic reactions is not ideal, since this assumes the activity of every metabolite is 1M. We propose that a better measure of the thermodynamic feasibility of reactions in biological systems is the standard Gibbs free energy change of reaction based on a 1mM reference state, 
calculated by setting every xi value in Eq. (1) equal to 1mM. For a reaction with the same number of reactants and products (
), not including hydrogen or water, 
is equal to 
If 
is not equal to zero, 
and 
can be substantially different. For example, for a reaction with one product molecule and two reactant molecules, such as threonine aldolase, | (2) |
of −1.9 kcal/mol, the 
value is 2.2 kcal/mol or 4.1 kcal/mol greater than 
Based on 
this reaction is thermodynamically favorable, but based on 
the reaction is mildly unfavorable. A second example is methylthioadenosine nucleosidase, | (3) |
of 2.3 kcal/mol, although it is favorable at 1mM activities with a 
of −1.7 kcal/mol. Depending on the value of 
the difference between 
and 
can be generalized as shown in Fig. 1.
Figure 1
Effect of transformation of 
into 
The difference between 
and 
is shown for different reaction molecularities. The difference between 
and 
depends only on the difference between the number of reactant molecules and the number of product molecules.Group contribution theory and estimation of ΔrG′°
Although experimental measurements of 
are unavailable for most compounds in E. coli metabolism, the group contribution methodology of Mavrovouniotis 17,18 provides a means by which the 
of most metabolites can be estimated providing the estimated 
or 
Group contribution methods consider a single compound as being made up of smaller structural subgroups. The Gibbs free energy changes associated with the set of structural subgroups, 
commonly found in metabolites, are available in the literature along with special corrections for complex biochemical cofactors such as coA and NAD+/NADH 17,18. To estimate 
of the entire compound, the contributions of each of the subgroups to this property are summed along with an origin term
 | (4) |
is an origin term common to all compounds, Ngr is the number of subgroups, ni is the number of instances of subgroup i in the compound, and 
is the contribution of subgroup i to 
17. All 
values calculated using the group contribution methodology of Mavrovouniotis are based upon the standard condition of a solution with pH equal to 7 and with zero ionic strength.For any reaction taking place in aqueous media, reactants will dissociate into several ionic forms 15,20. For example, ATP will dissociate and interconvert between the ionic forms: ATP4−, HATP3−, and H2ATP2−. In the cellular environment, the total amount of ATP present is the sum of all of these dissociated forms. In the fitting of thermodynamic energies of formation in the group contribution method of Mavrovouniotis, the total amount of ATP is represented by the single most common charged form found in a pH 7 solution, ATP4−17,18. Thus, the reaction for the hydrolysis of ATP into ADP and phosphate will be written as
 | (5) |
Using the group contribution methodology, we were able to determine 
for 531 (85.9%) of the 618 compounds in the genome-scale iJR904 metabolic model of E. coli, which allowed the calculation of 
for 770 (82.6%) of the 932 reactions in the model. Values of 
could not be determined for 87 compounds because these molecules contain substructures for which no energy value has been provided by Mavrovouniotis.
Construction of the iHJ873 E. coli model
To obtain a metabolic model of E. coli for which 
of every reaction could be calculated, the reactions in the iJR904 model containing compounds for which 
could not be calculated were lumped into single reactions and these compounds were eliminated. For example, in the following series of reactions,
 | (6) |
of compound B is unknown, we add the reactions involving B such that B is eliminated creating the lumped reaction of | (7) |
Classification of iHJ873 model reactions
We performed flux variability analysis (FVA) 22 to determine the reactions involved in the maximum production of biomass from glucose in E. coli under aerobic conditions. Details of all flux analysis performed are listed in the Appendix . Under optimal growth conditions, reactions in E. coli may be classified as essential (requiring a nonzero flux for optimal growth to occur), substitutable (capable of carrying zero or nonzero flux at optimal growth), or blocked (do not carry any flux at optimal growth). In the iHJ873 model, 250 (28.6%) reactions are essential, 51 (5.8%) reactions are substitutable, and 572 (65.5%) reactions are blocked. The total number of essential and substitutable reactions (301), which represents the total set of all reactions that participate in every alternative solution that produces optimal growth, agrees well with the average number of essential and substitutable reactions (294) reported for optimal growth phenotypes of E. coli utilizing a variety of nutrient sources 23. FVA also provides the direction of flux through the essential and substitutable reactions, allowing the reactants and products of all of these reactions to be redefined according to the direction of flux required for optimal growth (every flux will be positive). If a reaction can be active in both directions at optimal growth, the reactants and products and, consequently the reference directionality of the reactions, are defined according to their conventional nomenclature 16,24,25. Calculating 
using this definition of reactants and products in the reaction means that a positive 
value is indicative of a reaction that is thermodynamically unfavorable in the direction of flux required for optimal growth to occur at 1mM activity conditions.
Results
Distribution of ΔrGest′m values for reactions in iHJ873
The distributions of 
values for the essential and substitutable reactions in iHJ873, shown in histogram form in Figure 2AB, indicate that 80.4% of the reactions have a 
that is less than or equal to zero. However, there is uncertainty in 
Uf,est, based on the group contribution methodology. The value Uf,est is given as ∼±4 kcal/mol 18, and the standard error is used for the uncertainty in 
Ur,est, which is calculated as the Euclidean norm of the uncertainty for 
of each compound involved in the reaction (blue error bars in Figure 2CD) 26:
 | (8) |
as well as the associated ranges of uncertainty depend on reaction molecularity (Figure 2CD).
Figure 2
Thermodynamic feasibility of the reactions in iHJ873. Histograms of 
values for the essential (A) and substitutable (B) reactions. The values of 
for the essential (C) and substitutable (D) reactions. The blue error bars indicate the uncertainty of 
calculated, given a 4 kcal/mol uncertainty in 
provided in the literature 17. The red error bars indicate the range of values that 
could take considering the uncertainty and concentrations in the cell ranging from 20mM to 10−2 mM. The solid arrows in C mark the reactions for which 
−Ur,est>0. These reactions must be unfavorable at the reference conditions.The reactions in iHJ873 can be categorized thermodynamically based on their 
value and the associated Ur,est. In category (i), 
; 321 (36.8%) of all of the reactions, and 90 (29.9%) of the essential and substitutable reactions in iHJ873 are in this category. In category (ii), 
and 
and this category contains 429 (49.1%) of all of the reactions and 152 (50.5%) of the required and substitutable reactions. In category (iii), 
and 
; this category consists of 114 (13.1%) of all reactions and 54 (17.9%) of the substitutable reactions. In category (iv), 
and five (0.6%) of all of the reactions and four (1.3%) of the essential and substitutable reactions are in this category. There are four different reactions that generate biomass in iHJ873, and these reactions are not part of any category, since the 
of these reactions cannot be calculated.
The 
values for the reactions in categories (ii) and (iii) are relatively close to zero, indicating that these reactions are close to equilibrium at reference conditions. Only the five reactions in category (iv) must be unfavorable at the standard conditions and millimolar metabolite activities. If we examine the distribution of the 
values instead, we find that 232, 496, 138, and 3 of all of the reactions are in categories (i), (ii), (iii), and (iv), respectively. The smaller portion of reactions in the extreme categories (i) and (iv) indicates that the distribution of 
values is narrower than the distribution of 
values.
The values of 
can deviate from 
depending on how different the metabolite activities are from the reference value of 1mM. Metabolite activities can range approximately between 10−2 mM and 20mM 19. Based on these considerations, the maximum and minimum values for 
were calculated using the equations
 | (9) |
 | (10) |
(Eq. (8)), xmin is the minimal metabolite activity assumed to be 10−2 mM, and xmax is the maximum metabolite activity assumed to be 20mM 19. Although metabolite activities can be lower than 10−2 mM, a decrease in the lower limit on metabolite activity will result in an increase in all 
and a decrease in all 
(Fig. 2). Of the five reactions in category (iv), which have the highest 
values, metabolite activity profiles exist that can reduce 
and thus make these reactions thermodynamically feasible. These cases are indicated in Figure 2C with arrows on the right side of the corresponding graphs.The large fraction of essential and substitutable reactions in categories (ii) and (iii) with 
values within the margin of error of the zero axis indicates that most reactions involved in growth are energetically balanced, and only small concentration gradients are required to make these reactions thermodynamically feasible. The large fraction of reactions with associated 
values that are near zero is advantageous to the cell, because this prevents reactant and product concentrations from rising to toxic levels or falling to levels that would limit reaction rates.
ATP synthase and transport reactions
The standard conditions of pH 7 solution and zero ionic strength upon which all 
values are based was applied to both the extracellular and intracellular environment when calculating 
for reactions involving the transport of metabolites across the cellular membrane. As a result, 
for these reactions is based on the assumption that the electrochemical potential, Δψ, and pH gradient, ΔpH (pHintracellular – pHextracellular), across the cell membrane is zero. For example, the ATP synthase reaction in E. coli is typically written in the form of
 | (11) |
for the portion of this reaction that takes place inside the cell, 
 | (12) |
is 12 kcal/mol, which agrees well with the experimentally measured 
of 10.4 kcal/mol 27.The energy contribution of the transmembrane transport portion of the ATP synthase reaction, 
 | (13) |
and the energy associated with the transport of an ion across the membrane, 
 | (14) |
for ATP synthase is 0.0 kcal/mol.The overall 
of a reaction energetically coupled to the transport of an ion across the cell membrane such as ATP synthase is
 | (15) |
for the ATP synthase reaction (Eq. (11)) at the standard conditions is 12 kcal/mol.However, under physiological conditions ΔpH, Δψ and 
are not zero. The value 
depends upon Δψ, which in turn depends on ΔpH according to the equations 28
 | (16) |
 | (17) |
depends only on ΔpH according to the equation 28 | (18) |
of ATP synthase is −15.6 kcal/mol, making the total 
of ATP synthase −3.6 kcal/mol. The value of 
for the ATP synthase reaction only becomes positive when the ΔpH is lower than −0.51, meaning the extracellular pH is higher than the intracellular pH and above the optimal pH for E. coli growth.Identification and characterization of unfavorable reactions
Only five of the 873 reactions in the iHJ873 model have a 
that is greater than Ur,est, which indicates that every possible value of 
given the uncertainty in the estimate, must be positive and these reactions are unfavorable at standard conditions and 1mM activities. These five reactions are listed in Table 1. Four of these five unfavorable reactions are classified as essential for optimal growth to occur. These four reactions are
1. ATP phosphoribosyltransferase:
2. ATP synthase (without accounting for membrane potential, as discussed earlier; see Eq. (11)).
3. Methylenetetrahydrofolate dehydrogenase:
4. Tryptophanase:
We simulated knockouts of each unfavorable reaction while maximizing the yield of each of the biomass precursors and the yield of biomass to study the effects of single-knockouts and simultaneous knockouts on the cell growth (see Methods). Although ATP synthase is typically an energetically favorable reaction due to the energy contribution of the pH gradient and electrochemical potential across the cell membrane, ATP synthase is also included in the knockout studies to investigate the response of the system in case ATP synthase becomes unfavorable due to the failure of the proton gradient coupling and transmembrane potential. | (19) |
2. ATP synthase (without accounting for membrane potential, as discussed earlier; see Eq. (11)).
3. Methylenetetrahydrofolate dehydrogenase:
 | (20) |
4. Tryptophanase:
 | (21) |
| Table 1 Unfavorable reactions |
| Name | Pathway |  | Classification | ||
|---|---|---|---|---|---|
| Tryptophanase | Tyrosine, tryptophan, and phenylalanine metabolism | 13 | Essential | ||
| ATP synthase | Oxidative phosphorylation | 12 | Essential | ||
| Methylene-tetra-hydrofolate dehydrogenase | Folate metabolism | 9.9 | Essential | ||
| ATP phosphoribosyltransferase | Histidine metabolism | 8.2 | Essential | ||
| 2-C-methyl-D-erythritol 2,4-cyclodiphosphate synthase | Cofactor and prosthetic group biosynthesis | 22 | Blocked | ||
ATP phosphoribosyltransferase knockout
Only the precursor histidine is affected by the knockout of ATP phosphoribosyltransferase, and the production of histidine is not possible without the activity of this reaction, making histidine the limiting component preventing any growth without ATP phosphoribosyltransferase (Figure 3A). Experimental evidence confirms that ATP phosphoribosyltransferase is essential for the production of histidine, and mutant strains lacking this enzyme cannot grow without a histidine supplement 29. Experimental evidence also confirms that this reaction is thermodynamically unfavorable 29. The value 
can range between 0.2 and 16.2 kcal/mol given the margin of uncertainty in the group contribution methodology. If the metabolite activities in the cell range between 20mM and 10−2 mM, then 
of this reaction can range between −0.81 kcal/mol and 17.2 kcal/mol. Therefore, the reactant to product activity gradients required to drive this unfavorable reaction are achievable within the range of the physiological intracellular activities.
Figure 3
Impact of unfavorable reactions on the yield of biomass precursors. The maximum yields of the individual biomass precursors and of biomass relative to the corresponding maximum yields in the wild-type cells for different reaction knockouts: (A) ATP phosphoribosyltransferase, (B) ATP synthase, (C) methylene-tetra-hydrofolate dehydrogenase, and (D) tryptophanase. The bold lines in the plot show the relative precursor yield required for optimal growth to occur.As the first step in the histidine metabolism pathway, ATP phosphoribosyltransferase is an important point of control for the production of histidine. A mechanism even exists in the cell for feedback inhibition of ATP phosphoribosyltransferase by histidine 30. The unfavorable thermodynamics of this reaction provides another mechanism for product-inhibition of this enzyme as a means of limiting the flux that enters the histidine metabolism pathway.
ATP synthase knockout
The knockout of ATP synthase affects the optimal production of 49 of the 53 biomass precursors in the iHJ873 model (Figure 3B). The production of the energy in the form of ATP during aerobic metabolism depends heavily upon the ATP synthase reaction, and without ATP synthase, the energy requirements for optimal growth are not satisfied. Although a lack of ATP synthase does not completely prevent cell growth, the cell can only grow at 39.2% of the optimal yield, and experimental evidence confirms this effect of ATP synthase on growth 31.
Methylene-tetra-hydrofolate dehydrogenase knockout
Although all biomass precursors can still be produced individually in sufficient quantity for optimal growth to occur with the knockout of methylenetetrahydrofolate dehydrogenase, this knockout does reduce the production of 14 biomass precursors by an average of 8.6% (Figure 3C). As a result, no precursors can be produced simultaneously in sufficient quantities for optimal growth to occur without this reaction, and growth yield is reduced to 97.3% of the optimum. Methylene-tetra-hydrofolate dehydrogenase is a key step in the folate-dependent one-carbon metabolism pathway. The one-carbon pool from folate cannot be synthesized without this reaction, and without this reaction, other sources of C1 in metabolism must be utilized. According to the literature, this reaction is thermodynamically unfavorable with a 
of 1.17 kcal/mol 32, which is within the margin of uncertainty of the group contribution 
estimate of 9.94 kcal/mol. Given the range of physiological intracellular metabolite activities, 
can deviate from 
of 1.17 kcal/mol and range between −7.8 kcal/mol and 10.2 kcal/mol. The typical NADP/NADPH ratio found in E. coli is 6, and this ratio alone is already sufficient to reduce 
of this reaction by 1.06 kcal/mol to a 
of 0.11 kcal/mol.
Tryptophanase knockout
Only the maximum yield of the precursor tryptophan is affected by the knockout of tryptophanase (Figure 3D), and it is only reduced by 3.8%. Knockout out of tryptophanase has a nearly negligible effect on growth, reducing the yield by 0.03%. According to experimental evidence found in the literature, the 
for this reaction is −4.98 kcal/mol 33, which transforms into a 
value of 3.21 kcal/mol for this three-reactant, one-product reaction, confirming that this reaction is thermodynamically unfavorable under mM activity conditions. Although the estimate of 
from group contribution theory, 13 kcal/mol, for this reaction is high relative to experimental values, the difference between the estimate and the experimental data, 9.8 kcal/mol, still falls near the standard uncertainty of this reaction, 8.9 kcal/mol.
Four-reaction knockout
To determine the cumulative effect on biomass production of knocking out multiple unfavorable reactions simultaneously, we performed knockout simulations in which the activities of all of the unfavorable reactions were removed in every possible combination (Figure 4I). The effect of the cumulative knockouts on energy production was also examined (Figure 4II). In the simultaneous knockout of ATP synthase and tryptophanase, the growth yield is the same as the lower growth yield from the single knockouts of the same reactions. In this case, the knockout is not additive and the reactions play independent roles in the production of biomass. However, in the simultaneous knockout of ATP synthase and methylene-tetra-hydrofolate dehydrogenase, the growth yield is lower than the yield achieved from either of the single knockouts of these reactions. The effect of the double knockout of these reactions is additive, demonstrating that the contribution of these reactions to growth is linked.
Figure 4
Impact of simultaneous knockouts of unfavorable reactions on growth and growth-limiting precursors. The maximum yields of biomass (I) and energy (II) for the combinations of simultaneous knockouts of the unfavorable reactions: A, ATP phosphoribosyltransferase, B, ATP synthase, C, methylene-tetra-hydrofolate dehydrogenase, and D, tryptophanase. The maximum yields of biomass (I) are scaled by the maximum yield of biomass in the wild-type. The maximum yields of energy (II) are scaled by the amount of energy required for optimal growth. Not all combinations including knockout of ATP phosphoribosyltransferase (A) are shown because zero growth is possible without ATP phosphoribosyltransferase and energy is not affected by knockout of ATP phosphoribosyltransferase.Effect of unfavorable reaction knockouts on reaction classification
Unlike ATP phosphoribosyltransferase, the activities of the unfavorable reactions ATP synthase, methylene-tetra-hydrofolate dehydrogenase, and tryptophanase are not essential for cell growth to occur. These reactions have a wide range of effects on the metabolism of the cell, indicated by the widespread effect on the production of the biomass precursors. To study the effect a knockout of these reactions has on the distribution of fluxes in E. coli, FVA was utilized to determine how the reactions that are essential for optimal growth change when the activity of these reactions is knocked out. Table 2 summarizes the results of this study.
| Table 2 Minimal reaction sets for optimal growth in knockout and wild-type metabolism |
| Phenotype | Growth yield (gm biomass/mmol glucose) | No. of essential reactions | No. of substitutable reactions | ||
|---|---|---|---|---|---|
| Wild-type | 0.0923 | 250 | 51 | ||
| ATP synthase KO | 0.0362 | 234 | 118 | ||
| Methylene-tetra-hydrofolate dehydrogenase KO | 0.0897 | 249 | 37 | ||
| Tryptophanase KO | 0.0923 | 248 | 54 | ||
The wild-type and ATP synthase knockout share 231 essential reactions in common. There are 19 reactions that are essential in the wild-type and nonessential in the ATP synthase knockout. Fourteen of these reactions become substitutable in the ATP synthase knockout. These reactions are primarily clustered in the citrate cycle, glycolysis, and oxidative phosphorylation pathways. The remaining five of the 19 nonessential reactions in ATP synthase knockout, including ATP synthase, are blocked in the ATP synthase knockout. These reactions are found in the folate metabolism and pentose phosphate pathways.
Three reactions that are essential in the ATP synthase knockout are blocked in the wild-type. Two of these reactions are involved in producing threonine while consuming one ATP, and the third reaction is an acetate transporter. Overall, the knockout of ATP synthase results in a deactivation of portions of the pentose phosphate pathways, citrate cycle, and glycolysis.
The wild-type and methylene-tetra-hydrofolate dehydrogenase knockout share 243 essential reactions in common. Six of the reactions that are essential in the methylene-tetra-hydrofolate dehydrogenase knockout are blocked in the wild-type. These reactions are involved in a variety of small-carbon metabolism pathways. Many of these reactions produce formate to compensate for the loss of the formate metabolism reactions with the knockout of methylene-tetra-hydrofolate dehydrogenase. Five reactions, in addition to methylene-tetra-hydrofolate dehydrogenase, that are essential in the wild-type are blocked in the methylene-tetra-hydrofolate dehydrogenase knockout. These reactions are involved in the alternate carbon metabolism, arginine and proline metabolism, and folate metabolism pathways. These reactions are associated with the decomposition of some small-carbon compounds and the production of formate and tetrahydrofolate. Overall, knockout of methylene-tetra-hydrofolate dehydrogenase results in the deactivation of the folate metabolism pathway and the activation of alternative pathways for the production of folate and other small carbon compounds.
Comparing the essential reactions in the wild-type at optimal growth to the tryptophanase knockout at optimal growth, every essential reaction in the tryptophanase knockout is also essential in the wild-type knockout. Other than tryptophanase, only the reaction tryptophan synthase from the aromatic amino acid metabolism pathways is essential in the wild-type and not essential in the tryptophanase knockout. This reaction becomes substitutable in the tryptophanase knockout.
Discussion and conclusions
The group contribution methodology of Mavrovouniotis is demonstrated to be an effective means of estimating the free energy change of biochemical reactions. This methodology was utilized to calculate 
for 82.6% of the reactions in the iJR904 genome-scale metabolic model of E. coli developed by Palsson and co-workers. The iJR904 model was modified to eliminate the 85 compounds for which no group contribution estimation of 
was possible to create the iHJ873 model, and 
was determined for all of the reactions in iHJ873.
The 
is an invaluable measure of the thermodynamic feasibility of the reactions in the metabolic pathways of the cell under physiological conditions. Four-hundred-and-twenty-nine (49.1%) of all of the reactions and 152 (50.5%) of the reactions that are essential or substitutable for optimal growth to occur have a negative 
such that 
+Ur,est>0. The majority of the reactions in the cell are thermodynamically favorable, with a 
that is relatively close to zero under standard conditions and 1mM metabolite activities. This result indicates that the cellular system is energetically buffered from large perturbations and a minimal thermodynamic driving force is utilized to drive reactions.
Only four reactions essential for optimal growth yield have a positive 
such that 
indicating that these reactions must be unfavorable at standard conditions and 1mM metabolite activity levels. These four reactions are ATP phosphoribosyltransferase in the histidine metabolism pathway, ATP synthase in the oxidative phosphorylation pathway, methylene-tetra-hydrofolate dehydrogenase in the folate metabolism pathway, and tryptophanase in the aromatic amino-acid metabolism pathway. Experimental data exists that confirms that ATP phosphoribosyltransferase, ATP synthase, methylene-tetra-hydrofolate dehydrogenase, and tryptophanase are unfavorable at standard conditions and mM activities. These unfavorable reactions represent crucial thermodynamic bottlenecks in the production of growth, constraining the activities of metabolites involved in these reactions so that sufficient activity gradient may be provided to drive the reactions.
The fact that out of the 250 reactions essential for growth, only four reactions have 
values that are sufficiently large that they must be unfavorable at the reference conditions, indicates that the reactions involved in metabolism in E. coli are thermodynamically optimized to a great extent. It is important to note, however, that 
values discussed in this article are estimates and not experimentally measured values, and in some cases, like tryptophanase, the estimates can differ from the experimental values. This emphasizes the importance of accounting for the uncertainty in the group contribution estimates before utilizing this data for any analysis. Although the group energy values upon which 
are based were obtained from a fitting of experimentally measured 
values, this experimental dataset consists of far fewer reactions than are involved in the genome-scale model of E. coli.
The thermodynamic data obtained from this methodology is essential for the determination of the thermodynamically feasible activity ranges for the metabolites involved in the active reactions in E. coli metabolism, as discussed in the literature 15,34. Such feasible ranges would be very useful for narrowing the constraints utilized in constraints-based models as well as the operating conditions explored in MCA and kinetic modeling 5,35,36. The 
may also be used to formulate additional thermodynamic constraints for metabolic flux analysis (MFA) to ensure that flux distributions generated are thermodynamically feasible. Addition of thermodynamic constraints would aid in improving the predictions by metabolic models of the effect of gene knockout or other perturbations to the cellular metabolism. The error analysis discussed here will form an integral part of such thermodynamic constraints.
Acknowledgments
We thank Professor Bernhard Palsson and colleagues at the University of California at San Diego for making the iJR904 model readily available.
The work is supported by the United States Department of Energy, Genomes to Life Program.
Appendix
Details of flux analysis
Metabolic flux analysis (MFA) defines the limits on the metabolic capabilities of a model organism under steady-state flux conditions 3. Steady-state flux conditions are described by constraining the net production of every metabolite in the system, given by the product of the stoichiometric matrix and flux vector, to 0, as shown in Eq. (22),
 | (22) |
MFA studies were performed under a specific set of constraints on the metabolites the cell could uptake from or excrete to the cell surroundings. The ability of E. coli to grow optimally under aerobic conditions was studied using glucose as a primary carbon source. The uptake of glucose and oxygen from the environment into the cell was restricted to 10 and 20 mmol/g per dw per h, respectively 7. The uptake and excretion of sulfate, phosphate, and ammonium, CO2, water, and hydrogen ion were left unrestricted and the ATP maintenance requirement was fixed at 7.6 mmol/g per dw per h 21,37,38. Under these conditions, the optimal growth on glucose was found to be 0.923g biomass/g per dw per h, with a yield of 0.0923gram biomass per mmol of glucose uptake (0.512g biomass/g glucose). This optimal growth yield agrees well with the optimal growth yields for E. coli under similar conditions reported in the literature from MFA and experiments 38.
Flux variability analysis was used to classify the behavior of the reactions in the model using the methods described in the literature 22. There are 17 internal flux loops, or type-3 extreme pathways 14, in the iJR904 model, and 13 of these internal flux loops also exist in the iHJ873. No flux should move through these internal flux loops in the FVA flux distributions, because no thermodynamic driving force can exist for such flux. To prevent any flux from moving through these internal flux loops, one reaction from each internal flux loop is blocked in the FVA (13 blocked reactions total). A list of the reactions that must be blocked is found in the iJR904 literature 23 and in the Supplementary Material .
Supplementary material
An online supplement to this article can be found by visiting BJ Online at http://www.biophysj.org.
Appendix B
- Supplemental (XLS 292 kb)
- Supplemental
References
1. (2003). Identifying constraints that govern cell behavior: a key to converting conceptual to computational models in biology. Biotechnol. Bioeng. 84, 763–772. CrossRef | PubMed
2. (2002). Energy balance for analysis of complex metabolic networks. Biophys. J. 83, 79–86. Abstract | Full Text | PDF (569 kb) | PubMed
3. (2004). Monte Carlo sampling can be used to determine the size and shape of the steady-state flux space. J. Theor. Biol. 228, 437–447. CrossRef | PubMed
4. (1999). Metabolic flux ratio analysis of genetic and environmental modulations of Escherichia coli central carbon metabolism. J. Bacteriol. 181, 6679–6688. PubMed
5. (2004). Metabolic control analysis under uncertainty: framework development and case studies. Biophys. J. 87, 3750–3763. Abstract | Full Text | PDF (427 kb) | CrossRef | PubMed
6. (1985). Equations and calculations of product yields and preferred pathways for butanediol and mixed-acid fermentations. Biotechnol. Bioeng. 27, 50–66. CrossRef | PubMed
7. (1994). Stoichiometric flux balance models quantitatively predict growth and metabolic by-product secretion in wild-type Escherichia coli W3110. Appl. Environ. Microbiol. 60, 3724–3731. PubMed
8. (2001). Probing the performance limits of the Escherichia coli metabolic network subject to gene additions or deletions. Biotechnol. Bioeng. 74, 364–375. CrossRef | PubMed
9. (2000). Robustness analysis of the Escherichia coli metabolic network. Biotechnol. Prog. 16, 927–939. CrossRef | PubMed
10. (2005). Construction of lycopene-overproducing E. coli strains by combining systematic and combinatorial gene knockout targets. Nat. Biotechnol. 23, 612–616. CrossRef | PubMed
11. (1997). Flux analysis of underdetermined metabolic networks: the quest for the missing constraints. Trends Biotechnol. 15, 308–314. Abstract | | CrossRef | PubMed
12. (2003). Stoichiometric network theory for nonequilibrium biochemical systems. Eur. J. Biochem. 270, 415–421. CrossRef | PubMed
13. (2004). Thermodynamic constraints for biochemical networks. J. Theor. Biol. 228, 327–333. CrossRef | PubMed
14. (2000). Theory for the systemic definition of metabolic pathways and their use in interpreting metabolic function from a pathway-oriented perspective. J. Theor. Biol. 203, 229–248. CrossRef | PubMed
15. (2005). Thermodynamic-based computational profiling of cellular regulatory control in hepatocyte metabolism. Am. J. Physiol. Endocrin. M 288, E633–E644. PubMed
16. (2003). An expanded genome-scale model of Escherichia coli K-12 (iJR904 GSM/GPR). Genome Biol. 4, , 54.51–54.12. PubMed
17. (1990). Group contributions for estimating standard Gibbs energies of formation of biochemical-compounds in aqueous-solution. Biotechnol. Bioeng. 36, 1070–1082. CrossRef | PubMed
18. (1991). Estimation of standard Gibbs energy changes of biotransformations. J. Biol. Chem. 266, 14440–14445. PubMed
19. (1990). Cellular concentrations of enzymes and their substrates. J. Theor. Biol. 143, 163–195. CrossRef | PubMed
20. (1998). Calculation of standard transformed Gibbs energies and standard transformed enthalpies of biochemical reactants. Arch. Biochem. Biophys. 353, 116–130. CrossRef | PubMed
21. (1988). Bacterial Metabolism. (New York: Springer-Verlag). PubMed
22. (2003). The effects of alternate optimal solutions in constraint-based genome-scale metabolic models. Metab. Eng. 5, 264–276. CrossRef | PubMed
23. (2004). Genome-scale in silico models of E. coli have multiple equivalent phenotypic states: assessment of correlated reaction subsets that comprise network states. Genome Res. 14, 1797–1805. CrossRef | PubMed
24. (2000). KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 28, 27–30. CrossRef | PubMed
25. (1999). KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 27, 29–34. CrossRef | PubMed
26. (1978). Statistics for Experimenters: An Introduction to Design, Data Analysis, and Model Building. (New York: Wiley).