New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Analysis of optimality in natural and perturbed metabolic networks

Edited by Philip P. Green, University of Washington School of Medicine, Seattle, WA, and approved September 18, 2002 (received for review June 10, 2002)
Abstract
An important goal of wholecell computational modeling is to integrate detailed biochemical information with biological intuition to produce testable predictions. Based on the premise that prokaryotes such as Escherichia coli have maximized their growth performance along evolution, flux balance analysis (FBA) predicts metabolic flux distributions at steady state by using linear programming. Corroborating earlier results, we show that recent intracellular flux data for wildtype E. coli JM101 display excellent agreement with FBA predictions. Although the assumption of optimality for a wildtype bacterium is justifiable, the same argument may not be valid for genetically engineered knockouts or other bacterial strains that were not exposed to longterm evolutionary pressure. We address this point by introducing the method of minimization of metabolic adjustment (MOMA), whereby we test the hypothesis that knockout metabolic fluxes undergo a minimal redistribution with respect to the flux configuration of the wild type. MOMA employs quadratic programming to identify a point in flux space, which is closest to the wildtype point, compatibly with the gene deletion constraint. Comparing MOMA and FBA predictions to experimental flux data for E. coli pyruvate kinase mutant PB25, we find that MOMA displays a significantly higher correlation than FBA. Our method is further supported by experimental data for E. coli knockout growth rates. It can therefore be used for predicting the behavior of perturbed metabolic networks, whose growth performance is in general suboptimal. MOMA and its possible future extensions may be useful in understanding the evolutionary optimization of metabolism.
The enormous number of components and interactions in a cell, together with the uncertainty about many parameters describing cellular dynamics, greatly hinder the task of performing accurate whole cell simulations. Consequently, computational efforts based on conceptual shortcuts are essential. One area in which such simplifications have proved extremely useful is metabolic flux analysis (1–7). Notably, flux balance analysis (FBA) (8–11), a method for studying the capabilities of metabolic networks at steady state, constitutes an example of how the knowledge of a restricted set of parameters in a system, combined with the application of fundamental thermodynamic and evolutionary principles, can generate quantitative predictions and testable hypotheses. In FBA, the constraints imposed by stoichiometry in a chemical network at steady state are treated analogously to Kirchoff's law for the balance of currents in electric circuits (2, 12). Thus, for each of M metabolites in a network, the net sum of all production and consumption fluxes, weighted by their stoichiometric coefficients, is zero: Here, S_{ij} is the element of the stoichiometric matrix S corresponding to the stoichiometric coefficient of metabolite i in reaction j. The flux v_{j} is the rate of reaction j at steady state, and is the jth component of an Ndimensional flux vector v, where N is the total number of fluxes. In addition to internal fluxes, which are associated with chemical reactions, v includes exchange fluxes that account for metabolite transport through the membrane. The steadystate approximation is generally valid because of the fast equilibration of metabolite concentrations (seconds) with respect to the time scale of genetic regulation (minutes) (1, 6).
Additional constraints, including those that relate to the availability of nutrients or to the maximal fluxes that can be supported by enzymatic pathways, can be introduced as inequalities For example, for a substrate uptake flux v_{j}, one can set α_{j} and β_{j} equal to the corresponding measured or imposed value. Eq. 2 can also be used to distinguish reversible and irreversible reactions, where α_{j} = 0 for the latter. Additional constraints are invoked to represent the requirement for metabolic homeostasis, and can be expressed in terms of linear relationships similar to Eq. 1 (8, 13).
All flux vectors that satisfy the constraints mentioned above define a feasible space, Φ. For an underdetermined system, as is typically the case in FBA models of cellular metabolic networks (11), Φ is a convex set in the Ndimensional space of fluxes (14). Because of the linear nature of Φ, it is possible to use linear programming (15) to characterize the points in Φ that maximize or minimize a given linear objective function. A natural choice for an objective function in metabolic models of prokaryotes is biomass production (8, 9, 13), as it is reasonable to hypothesize that unicellular organisms have evolved toward maximal growth performance. This process is formalized by introducing a growth flux that transforms a linear combination of fundamental metabolic precursors into biomass (see Methods).
The availability of annotated genomic sequences has led to FBA metabolic models of various microorganisms, including Haemophilus influenzae (16), Escherichia coli (13, 17, 18), and Helicobactyer pylori (19). The existence of ≈100 fully sequenced and annotated genomes (and many more in pipeline) paves the way for widescale application of flux analysis of the corresponding metabolic networks (5).
The theoretical basis of FBA is supported by several experiments. These include empirical validation of growth yield and flux predictions (8, 9), measurements of uptake rates around the optimum under various conditions (18), as well as results from largescale gene deletion experiments (20). Additional strong support based on intracellular flux comparisons is presented here in Fig. 4 A, D, and G.
An important application of FBA is the prediction of phenotypic effects arising from complete or partial metabolic gene deletions (13, 17, 21). A complete gene deletion is implemented by constraining the corresponding flux to zero. Linear programming provides then the flux distribution and maximal growth yield for the new genotype. Crucially, this approach assumes that the mutant bacteria display an optimal metabolic state; yet, mutants generated artificially in the laboratory are generally not subjected to the same evolutionary pressure that shaped the wild type. Therefore knockouts probably do not possess a mechanism for immediate regulation of fluxes toward the optimal growth configuration. To better understand the flux states of mutants, we introduce here the method of minimization of metabolic adjustment (MOMA) (Fig. 1A), which is based on the same stoichiometric constraints as FBA, but relaxes the assumption of optimal growth flux for gene deletions. A mutant is likely to initially display a suboptimal flux distribution that is somehow intermediate between the wildtype optimum and the mutant optimum. MOMA provides a mathematically tractable approximation for this intermediate suboptimal state, based on the conjecture that the mutant remains initially as close as possible to the wildtype optimum in terms of flux values. In other words, through MOMA, we test the hypothesis that the real knockout steady state is better approximated by the flux minimal response to the perturbation than by the optimal one. Predicting a metabolic phenotype by MOMA involves a different optimization problem than FBA, namely distance minimization in flux space. To this end, because distance in flux space is a quadratic function, we employ quadratic programming (QP) (see also ref. 22). In the present work we compare FBA and MOMA predictions for E. coli with experimental data from different sources, and find that MOMA can serve as a more precise method for predicting the metabolic phenotype of gene knockout bacteria. In addition to proposing MOMA as a computational tool for the analysis of artificially engineered bacteria, we suggest that this approach can help understand how a cell adapts to the loss of a gene by regulation and evolutionary optimization.
Methods
For the stoichiometric analysis of the metabolic network of E. coli, we have used the reconstruction by Edwards and Palsson (13). The list of metabolic reactions, and the 436 (metabolites) by 720 (fluxes) stoichiometric matrix (available at http://gcrg.ucsd.edu) were compiled using data from public databases and literature. This in silico model of E. coli MG1655 metabolism encompasses central carbon metabolism, transmembrane transport reactions, carbon source utilization pathways, as well as the metabolic pathways responsible for the synthesis and degradation of amino acids, nucleic acids, vitamins, cofactors, and lipids (13). As in previous FBA formulations, we use inequalities (Eq. 2) to limit nutrient uptake and to implement reaction irreversibility, and we include all fluxes related to maintenance requirements (8, 18). Reconstruction of other E. coli strains analyzed in this work, JM101 and PB25 (4), required minimal changes with respect to MG1655 (see supporting information, which is published on the PNAS web site, www.pnas.org).
Flux Balance Analysis.
In FBA the maximization of biomass production is implemented by defining an additional flux v_{gro} associated with cell growth. For this flux, the stoichiometric factors of the reactants are the experimentally known proportions c_{i} of metabolite precursors X_{i} contributing to biomass production (13, 23): The search for the flux vector maximizing v_{gro} under the constraints of Eqs. 1 and 2 is solved by the simplex algorithm. We denote here by v^{WT} the solution of FBA for the wildtype organism in its feasible space Φ^{WT}. A knockout of the enzyme catalyzing reaction j is obtained by imposing the extra constraint v_{j} = 0. Such condition defines the feasible space Φ^{j} for the mutant strain. The solution for knockout j is denoted by v^{j}.
Our perl implementation of FBA, whose results we verified matching previous FBA calculations, uses the open source GNU Linear Programming kit (24).
Minimization of Metabolic Adjustment.
In MOMA we search for a point in Φ^{j} that has minimal distance from a given flux vector w (Fig. 1A). The goal is to find the vector x ∈ Φ^{j} such that the Euclidean distance is minimized (22). Stated as a standard QP problem (15, 25), the aim is to minimize under a set of linear constraints. In Eq. 5, the vector L of length N and the N × N matrix Q define the linear and quadratic part of the objective function, respectively, and x^{T} represents the transpose of x. By observing that minimizing the function D of Eq. 4 is equivalent to minimizing its square, and that constant terms can be omitted from the objective function, one can choose Q to be an N × N unit matrix and set L = −w, and hence reduce the minimization of D to the minimization of f(x). Here we are interested in the case w = v^{WT}, i.e., in finding the point u^{j} in Φ^{j}, which is closest to the wildtype point (as calculated with FBA). Although a formal proof is not offered in this context, it can be intuitively inferred that, if a solution for the wildtype FBA problem v^{WT} exists, and if the space Φ^{j} is not empty (i.e., the constraint v_{j} = 0 is compatible with the other constraints), then a solution to this problem always exists. Its uniqueness is guaranteed by the convexity of f(x), which in turn derives from Q being semipositive definite (15). We emphasize that in MOMA, in contrast to FBA, the objective function does not explicitly depend on biomass production. Its linear part reflects the vector of fluxes of the wild type, whereas its quadratic part is the square of the Euclidean norm of x.
For QP, the IBM QP Solutions library (26) was used. The growth yield for the mutant can be evaluated after minimization by inspection of the corresponding flux (u). Time required for calculating one FBA or MOMA solution is of the order of 2 sec on a Pentium III computer running Linux.
Robustness Analysis.
Two methods were used to assess the stability of MOMApredicted flux points. The first (Fig. 1B) is analogous to the robustness analysis previously described for FBA (21), whereby one follows the gradual change in growth flux as v_{j} is varied from v to zero. In MOMA one analogously can study the gradual variation of the suboptimal flux point u^{j}, and in particular of the growth flux u. Comparative FBA and MOMA examples of robustness analysis can be found in the supporting information. Interestingly, there are cases in which FBA and MOMA predictions are identical for a complete knockout, but they differ from each other for partial deletions.
A second kind of robustness analysis (Fig. 1C) is unique to MOMA, and addresses the variability of the knockout prediction as the point w is displaced around the wildtype optimum v^{WT}. To obtain a map of this variability in the multidimensional flux space, it is necessary to perform a systematic analysis of the effects of displacements along all possible edges departing from v^{WT}. For this purpose, one can sample the flux space around v^{WT} by limiting one by one each of the fluxes to a given fraction (e.g., 90%) of its wildtype FBA value; on each such displacement one can recalculate the FBA optimum and use the result as the input w vector for MOMA. This sensitivity analysis can also be used for studying the dependence of MOMA predictions on the possible presence of alternative FBA optima (e.g., points a and a′ in Fig. 1C). The results of this analysis can be found in the supporting information.
Results
Prediction of Essential Genes in E. coli Central Carbon Metabolism.
We compared the predicted growth yield values of E. coli central carbon metabolism enzyme knockouts obtained by FBA and MOMA. For metabolic reactions carried out by multisubunit complexes, all members of the complex were deleted; similarly, for reactions catalyzed by multiple isoenzymes, all corresponding genes were deleted (13).
The normalized growth yields obtained with the two methods are plotted against each other in Fig. 2. All points lie in the lower triangle of the graph because MOMA necessarily gives growth yields smaller or equal to the optimal ones predicted by FBA. In 26 of 47 cases, the two predictions differ by <5%. Of the remaining 21 knockouts, 7 display a difference that exceeds 50%; among those whose predictions differ most between the two methods are fructose1,6bisphosphatate aldolase (fba), triosphosphate isomerase (tpiA), and phosphofructokinase (pfk). These were predicted to be nonessential by FBA, whereas the literature reported that each of them is essential for growth on glucose (see ref. 13 and references therein). Some flux change differences between FBA and MOMA predictions for the tpiA knockout are shown in Fig. 3B, using the pathway color code of Fig. 3A. On deletion of tpiA, the FBA optimal solution would tend to redistribute fluxes so as to increase significantly the use of the pentose phosphate pathway, whereas the MOMA point would tend to reflect the initially less prominent usage of this pathway. A number of related flux differences are observed, including a drastic divergence in the fluxes for the serine pathway genes serA, serB, and serC. For the wild type, these fluxes amount to ≈14% of the glucose uptake flux; in the FBA prediction for tpiA knockout the values are reduced to 3%, whereas in the MOMA prediction they are zero. This finding explains the corresponding zero final growth yield, because serine is one of the essential biomass components. A similar situation was observed for deletion of fba and pfkAB.
Comparison with Measured Growth Performance of Competing Insertional Mutants.
The FBA prediction of lethality for metabolic gene knockouts has been recently compared with the data obtained in an experiment of competitive growth (20). In this experiment, an E. coli transposon library was subjected to competitive growth in rich and minimal media, and 488 genes were classified according to whether their insertional inactivations negatively affected growth in the competition. The significance of the FBA prediction was assessed by a χ^{2} test, comparing the observed values of negatively and nonnegatively selected genes in the different predicted classes (essential, nonessential, and reduced growth), with the outcome expected if those genes were randomly distributed among the classes. The predicted trend of enrichment was observed (20). We repeated the same test with MOMA, and found the same trend, now with even higher statistical significance. The P value in the χ^{2} test for such enrichment trend was 10^{−5} in MOMA, as compared with the value of 4 × 10^{−4} in the FBA prediction (see Table 2, which is published as supporting information on the PNAS web site).
Comparison with Experimental Flux Measurements of an E. coli Pyruvate Kinase (pyk) Knockout.
Recent measurements of metabolic fluxes in a pyk knockout (4) allowed us to test the MOMA predictions directly. Emmerling et al. (4) empirically determined a collection of intracellular fluxes for E. coli central carbon metabolism by combining NMR spectroscopy in ^{13}C labeling experiments, physiological data measurement, and numerical data fitting. For both the wildtype strain (JM101) and the pyk mutant (PB25, pykA:kan pykF:cat), two glucoselimited (low and high concentration), and one nitrogenlimited experiments were performed, and six resulting sets of fluxes were measured (see Tables 2 and 3). For a discussion of the experimentally observed differences between the knockout and the wildtype phenotype, see ref. 4. A mathematical model of these strains had been discussed in ref. 35.
In the FBA and MOMA models, uptake fluxes for glucose and nitrogen were matched to the ones experimentally measured. We first compared the FBA prediction for the wild type with the corresponding experimental result. The remarkably high, and statistically significant correlation coefficients confirm the predictive value of FBA for wildtype E. coli, as seen in Fig. 4 A, D, and G and in Table 1. Correlation coefficients range between 0.78 and 0.97, with P values ranging from 7 × 10^{−5} to 8 × 10^{−12}. The lowest correlation is observed for the nitrogenlimited case. To our knowledge, this is the first time that several FBApredicted intracellular fluxes have been compared with experimental flux data for an array of different strains and conditions. Next, we compared the FBA and MOMA predictions for the pyk mutant with the corresponding experimental results. Overall, the correlation coefficients were lower than those of the wild type. Interestingly, however, the MOMA predictions in the two carbonlimited experiments were significantly more accurate than those of FBA (Fig. 4 B, C, E, and F). The lowest correlation (−0.064, P value = 6 × 10^{−1}) is found in the FBA prediction for the mutant under low glucose uptake conditions (Fig. 4B). MOMA in the same case (Fig. 4C) improves the correlation to 0.59, with P value = 7.4 × 10^{−3}. When repeating the correlation analysis for relative flux changes between knockout and wild type (u − v vs. v − v for flux i in knockout of gene j) we found that MOMA scores significantly better than FBA in all three conditions (see Table 1).
A peculiar distribution is observed for fluxes 16 (ppc) and 17 (pck) in Fig. 4 B, C, H, and I. These deviate from the main diagonal by comparable amounts, suggesting that the discrepancy may be related to the circular topology of their subnetwork (see Fig. 3A and ref. 4).
The overall higher predictive power of MOMA relative to FBA in the analysis of the pykdeleted strain is consistent with the notion that the assumption of optimality is less appropriate for knockout strains than for the wild type.
Discussion
Natural selection can be extremely efficient at generating systems that are optimally adapted to their environment (27–30). Optimal metabolic adaptation for growth under various environmental conditions is especially important for many bacteria (23). Optimality can be achieved by changing the topology of the metabolic network and by tuning enzymatic and regulatory parameters. FBA predicts growth yield and metabolic fluxes based on the assumption that growth efficiency has evolved to an optimal point. The high correlation that we find between FBA predictions for wildtype E. coli and some recent experimental flux measurements (4) confirms earlier evidence that the assumption of optimality is valid in this case. The slightly lower correlation observed in the nitrogenlimited environment relative to the carbonlimited cases could indicate that E. coli evolved toward optimality primarily under carbonlimited conditions.
In contrast to FBA, MOMA does not assume optimality of growth or of any other metabolic function. Instead, for perturbations such as gene deletions, MOMA approximates metabolic phenotype by performing distance minimization in flux space. We demonstrated that MOMA correctly predicts the lethality of some E. coli gene deletions that FBA failed to recognize, and we analyzed flux changes possibly crucial in this prediction. It will be interesting to analyze the flux distributions of other gene deletions whose points in Fig. 2 either show a discrepancy between FBA and MOMA predictions (e.g., eno and nuoA), or deviate from experimental results in all methods (e.g., aceE). Fig. 2 could serve in general as a map for suggesting maximally informative sets of experimental knockout flux data. It is likely that future experiments will produce more data for the intracellular fluxes of several E. coli gene knockouts, and hence more facts against which MOMA predictions could be tested. As part of this analysis, a systematic comparison of FBA and MOMA flux values would provide quantitative constraints on the complex metabolite and enzyme concentration changes expected on optimization of growth performance. To understand whether a transition from a suboptimal to an optimal metabolic state requires Darwinian evolution (33), or could be accomplished gradually through gene regulation, high timeresolution experiments are required, in which fluxes are measured before and after an enzyme is abruptly inactivated. Interestingly, in line with the basic hypothesis underlying MOMA, a previous study by Schuster and Holzhütter (31) revealed relative stability of erythrocyte metabolic fluxes against alterations in enzyme activity, such as the ones caused by gene mutation. A complete and detailed kinetic model, as explored for erythrocytes (31, 32), is probably beyond reach for E. coli at this stage (35). However, FBA and MOMA studies of erythrocytes, in parallel to kinetic models, could help relate the transition between different steady states to enzyme expression time courses. In addition, flux balance results and enzyme expression data can be in principle used to infer kinetic constants for a whole metabolic network, provided that enough independent nutrient conditions are considered.
The wildtype point that is used in MOMA does not necessarily need to be an FBA result. For example, a wildtype point determined experimentally could be used similarly. In an iterative process of gene deletions, several MOMA projections could be applied in a chain, and compared with the MOMA result for simultaneous deletion of multiple genes. The concept of suboptimal mutants could be also applied to heterologous gene additions (3), as well as to the prediction of metabolic flux distributions for bacteria grown under unnatural environmental conditions (33). Possible variants of the MOMA approach presented here could involve alternative flux space metrics. NonEuclidean distances, in alternative to Eq. 4, may be considered more suitable to representing some biological processes. For example, in alternative to the absolute flux changes currently used in the distance minimization algorithm, it is possible to use normalized flux changes. Such a method would remove the dependence of fluxes on constant factors, and hence may better capture the effects of gene and metabolite regulation on fluxes in mutant strains.
Among the various possible methods for identifying a single flux point given an underdetermined set of constraints (11), FBA employs the biological concept of optimal growth. Other biologically sound assumptions have been proposed for this purpose (5, 34). MOMA is an additional option, based on a simple hypothesis about the response to metabolic alterations. This approach seems especially relevant for analyzing gene deletions, but its possible future extensions could help understand metabolic networks for a wider range of perturbations.
Acknowledgments
We are grateful to J. Edwards for sharing with us his FBA experience, and to B. Palsson and U. Sauer for helpful comments. Many thanks to N. Reppas, Y. Pilpel, A. Dudley, K. Leptos, P. D'haeseleer, A. FalcovitzSegrè, and the whole Lipper center for invaluable help with manuscript preparation, and to D. Gamarnik and A. Egun for QP advice. This work was suported by grants from the U.S. Department of Energy and Defense Advanced Research Planning Agency.
Footnotes

↵* To whom correspondence should be addressed. Email: church{at}arep.med.harvard.edu.

This paper was submitted directly (Track II) to the PNAS office.
Abbreviations

MOMA, minimization of metabolic adjustment

FBA, flux balance analysis

QP, quadratic programming
 Received June 10, 2002.
 Copyright © 2002, The National Academy of Sciences
References
 ↵
 Heinrich R.
 ↵
 ↵
 Stephanopoulos G.
 ↵
 Emmerling M.
 ↵
 ↵
 Fell D.
 ↵
 ↵
 Varma A.
 ↵
 ↵
 ↵
 Edwards J. S.
 ↵
 ↵
 Vanderbei R. J.
 ↵
 ↵
 ↵
 ↵
 Schilling C. H.
 ↵
 ↵
 ↵
 ↵
 Neidhardt F. C.
 ↵
 Makhorin A.
 ↵
 Gill P. E.
 ↵
IBM, (1998) QP Solutions Library (IBM, White Plains, NY).
 ↵
 Freeland S. J.

 Akashi H.
 ↵
 ↵
 Jamshidi N.
 ↵
Ibarra, R. U., Edwards, J. S., Ramos, L. S. & Palsson, B. O. (2002) Nature, in press.
 ↵
 ↵