INTRODUCTION AND BACKGROUND

TOP ABSTRACT INTRODUCTION AND BACKGROUND THEORY OF ENERGY BALANCE... NONEQUILIBRIUM STEADY-STATE... APPLICATION OF ENERGY BALANCE... ANALYSIS OF ESCHERICHIA COLI... ESTIMATION OF REACTION... DISCUSSION AND CONCLUSIONS REFERENCES

Conservation principles impose constraints on the fluxes and chemical potentials associated with biochemical network reactions that are analogous to Kirchoff's current and voltage laws for electrical networks (Balabanian and Bickart, 1981). Flux balance analysis (FBA) Varma and Palsson, 1994; Bonarius et al., 1996, 1997; Pramanik and Keasling, 1997; Schilling et al., 1999, 2001; Schuster et al., 1999; Edwards and Palsson, 2000a, b, c; Edwards et al., 2001 invokes mass conservation, but does not consider the chemical potential associated with nonequilibrium steady-state chemical fluxes. The thermodynamic analysis presented here, energy balance analysis (EBA), provides additional constraints on the system that are analogous to the voltage loop law. In addition to predicting network fluxes that are thermodynamically feasible, EBA facilitates a detailed analysis of regulation of metabolic networks that is not available from FBA alone. These simple and fundamental theoretical results will have consequences in the analysis and manipulation of genomes, as we demonstrate for the genes involved in the Escherichia coli metabolism.

The central flux balance conservation statement is given by the equation

(1)

where J  ⁿ is the vector of n fluxes occurring in the network, S  ^m×n is the stoichiometric matrix, and m is the number of reactants in the system. The matrix S stores the stoichiometric coefficients associated with each flux in the network. In the above formulation, both internal fluxes and boundary fluxes, which transport material into or out of the system, are included in S. Typically, a number of inequalities are introduced to constrain the boundary (also called injection) fluxes depending upon the external media (Edwards and Palsson, 2000a, b, c; Edwards et al., 2001; Ramakrishna et al., 2001),

(2)

where _i and _i represent the upper and lower bounds on the flux J_i.

As recently implemented with significant success (Bonarius et al., 1997; Edwards et al., 2001; Ramakrishna et al., 2001; Schilling et al., 2001), biological networks are assumed to optimize a certain biologically meaningful objective function, which is a linear combination of the fluxes,

(3)

For example, in the analysis of E. coli intermediary metabolism, the objective is constructed from the production of biosynthetic precursors required to generate biomass (Edwards and Palsson, 2000a). Another example uses maximization of ATP production to simulate mitochondrial energy metabolism (Ramakrishna et al., 2001). Regardless of the application, optimization of a linear objective function (Eq. 3), together with the equality constraints (Eq. 1) and the inequality constraints (Eq. 2), represents a linear optimization problem, which can be solved with linear programming (Strang, 1986).

The power of FBA is illustrated by the tour de force assembly of the complete stoichiometry of the known reactions in E. coli metabolism provided by Edwards and Palsson (2000a, b, c). Edwards and Palsson show that the flux balance simulation of the organism-scale metabolic network predicts the metabolic capabilities of E. coli. Under various external constraints (e.g., aerobic versus anaerobic growth), FBA can distinguish between genes essential and not essential for growth in 68 of 79 cases studied (Edwards and Palsson, 2000a). This predictive capability of FBA is striking considering the tremendously complex problems that can be modeled using little or no free parameters.

	THEORY OF ENERGY BALANCE ANALYSIS

TOP ABSTRACT INTRODUCTION AND BACKGROUND THEORY OF ENERGY BALANCE... NONEQUILIBRIUM STEADY-STATE... APPLICATION OF ENERGY BALANCE... ANALYSIS OF ESCHERICHIA COLI... ESTIMATION OF REACTION... DISCUSSION AND CONCLUSIONS REFERENCES

Lacking from FBA is the explicit consideration of the energy balance and thermodynamics of the network reactions. Because biochemical networks are composed of multiple-species reactions, energy-balance loop equations cannot be obtained from topological loops in the network, as is done in electrical circuit analysis (Balabanian and Bickart, 1981). We will show that the energy-balance equations are obtained from an analysis of the network stoichiometry.

If we combine redundant fluxes and remove the columns from S that correspond to boundary fluxes, the remaining matrix, denoted S', represents the complete set of possible internal chemical transformations. Using the singular value decomposition (Strang, 1986), S' can be decomposed as

(4)

where  has the form,

(5)

The first m columns of  form a diagonal matrix where the diagonal entries are the singular values _i. The entries of the remaining columns are all equal to zero. If n_c is the number of columns of S' and r is the rank of S', then columns r + 1 through n_c of the matrix B provide the (n_c  r)-dimensional null space of S'. We introduce the matrix K, which stores the null space vectors of S',

(6)

and we denote the ith row of K as _i. Previous analyses based on the null space of stoichiometric matrices are presented by Alberty (1991) and Schilling and Palsson (1998).

Summing the n_c chemical reactions, each scaled by the corresponding entry of _i, results in a perfectly balanced reaction equation (with same reactants in equal proportion on either side of the equation). An example of this analysis for a simple five-metabolite network with four reactants is illustrated in Fig. 1.

View larger version (44K): [in this window] [in a new window]   FIGURE 1   Illustration of energy-balance equations for a network of five reactions. The first step is determination of the stoichiometric matrix from the reactions in the network. For this example, the stoichiometric matrix has rank r = 3, resulting in two independent null-space vectors. The null-space vectors provide mutually orthogonal solutions to S'J = 0. In addition, the null-space vectors are the rows of the matrix K, the energy-balance matrix for which Kµ = 0.

To study network thermodynamics, we consider the vector µ that lists the n_c chemical potential differences associated with the reaction fluxes. Because each _i provides weights for exactly balancing the chemical reaction equations (see Fig. 1), solutions to the equation

(7)

balance the global potential energy of the network. The second law of thermodynamics takes the form of an inequality constraint for each flux,

(8)

which ensures that entropy production is positive for each reaction (Qian, 2001a, b). Eqs. 7 and 8 represent thermodynamic constraints that should be considered in addition to the flux-balance constraints. Eq. 8 provides a link that couples mass balance and energy balance, and constraints the feasible flux space more strictly than Eqs. 1 and 2 alone.

The relationship between EBA and FBA is illustrated in Fig. 2 with an analogy to electric circuit analysis. Flux balance (or Kirchoff's current law) constrains the current through the two parallel passive resistors to sum to J₁ + J₂ = J, the total current into the circuit. Applying FBA alone, current can travel in either direction in either branch of the circuit. However, if we apply Kirchoff's loop law (energy balance), we find that the first two cases in Fig. 2 must be excluded, and the currents J₁ and J₂ must travel in the same direction.

View larger version (13K): [in this window] [in a new window]   FIGURE 2   Electric circuit analogy for energy balance analysis. EBA tells us that the currents through the two parallel passive resistors must be in the same direction. Therefore, the first two possibilities, with current flow in opposite directions, are excluded.

	NONEQUILIBRIUM STEADY-STATE THERMODYNAMICS

TOP ABSTRACT INTRODUCTION AND BACKGROUND THEORY OF ENERGY BALANCE... NONEQUILIBRIUM STEADY-STATE... APPLICATION OF ENERGY BALANCE... ANALYSIS OF ESCHERICHIA COLI... ESTIMATION OF REACTION... DISCUSSION AND CONCLUSIONS REFERENCES

In studying network thermodynamics, it is essential to differentiate between free energy of equilibrium and the chemical potential of a nonequilibrium steady state (NESS). Although the basic thermodynamic theory for the equilibrium case can be found in a textbook, it is only recently that the theory for nonequilibrium steady states has been cogently established. See Qian (2001a, b).

Consider the following examples of how chemical potential µ is balanced in a NESS system. For a given NESS reaction, AB, there is a nonzero flux J and nonzero µ = µ^o + k_BT ln([B]/[A]), where k_B represents Boltzmann's constant and T is absolute temperature. For each molecule of A that reacts to form B, µ equals the change in internal energy plus the change in entropy. The internal energy change is also exactly the amount of heat dissipated (or absorbed). Now consider a cycle, or loop,
(see Fig. 3). The summation of chemical potential µ₁ + µ₂ + µ₃ = 0 is obvious. Alternatively, if the reaction BC is coupled to the reaction
(i.e., the full reaction is B + DC + E) with clamped (externally fixed) concentrations for D and E, then the summation of chemical potential gives µ₁ + µ₂ + µ₃ = µ₄.

View larger version (15K): [in this window] [in a new window]   FIGURE 3   Energy balance for simple loops. The cycle ABCA is illustrated in the left panel. There can be no steady-state flux for this cycle because the chemical potentials µ1, µ2, and µ3 must have the same sign for the flux to be nonzero, while energy balance maintains µ1 + µ2 + µ3 = 0. When a driving force is added on the right panel, the energy-balance equation becomes µ1 + µ2 + µ3 = µ4, and a nonzero cycle flux is possible.

When a molecule undergoes the reaction AB, with fixed concentrations [A] and [B], then the amount of heat dissipated per turnover is precisely equal to the amount of work required to sustain constant concentrations [A] and [B]. For a given flux J, total heat dissipation is  J · µ, which is always positive. When the energy balance equation for the loop illustrated in the right panel of Fig. 3 is multiplied by the flux, we obtain the equation Jµ₁ + Jµ₂ + Jµ₃ = Jµ₄, which represents the balance between energy dissipated by the circuit and energy supplied to the system. The energy supplied by pumping D (and E) molecules into (and out of) the system is given by Jµ₄, which is dissipated by the heat associated with the cycle flux, Jµ₁ + Jµ₂ + Jµ₃. For a reaction cycle (or resistance loop), the amount of heat dissipated by all of the components is equal to the energy supplied by concentration clamping (or the battery in an electrical circuit).

To determine effective loops in a more complex reaction network, we turn to the null space decomposition of the stoichiometric matrix. The null space decomposes the internal reactions in the system into generalized loops with no net change in stoichiometry from the left side to the right side of the summed reaction (see Fig. 1). For the example in Fig. 1, summing the five reactions, weighted by the components of the first null space vector for this system, yields the exactly balanced reaction,

The µ for this reaction is obviously identically equal to zero, regardless of the concentrations of the reactants. Thus, Eq. 7 follows directly from the definition of chemical potential.

	APPLICATION OF ENERGY BALANCE ANALYSIS

TOP ABSTRACT INTRODUCTION AND BACKGROUND THEORY OF ENERGY BALANCE... NONEQUILIBRIUM STEADY-STATE... APPLICATION OF ENERGY BALANCE... ANALYSIS OF ESCHERICHIA COLI... ESTIMATION OF REACTION... DISCUSSION AND CONCLUSIONS REFERENCES

The EBA constraints are summarized as follows: For a flux vector J to be considered thermodynamically feasible, there must exist a vector µ, for which Kµ = 0; and µ_i · J_i < 0, for all µ_i and J_i, where µ_i and J_i represent the ith entries of µ and J, respectively, and the rows of the matrix K correspond to a complete basis of null-space vectors of the stoichiometric matrix of internal reactions S'. The EBA constraints must be applied in addition to the FBA constraints to ensure that the predicted flux vector is thermodynamically feasible.

As an example, we return to the simple reaction network illustrated in Fig. 1, and consider the following problem: Determine the maximum steady-state production of reactant D, for a given set of maximal input fluxes of reactants A, B, and C. (Three related problems are considered below. The first problem assumes that reactant A is the only available input substrate. For this case, the maximal input fluxes of A, B, and C are set to be 1, 0, and 0 mmol sec¹, respectively. The remaining two examples use B only and C only as input substrates.)

Optimal fluxes are obtained using the MATLAB (The Mathworks Inc., Natick, MA) software package. Production of reactant D is used as an objective function for optimization of flux vectors. For predictions based on FBA alone, linear programming is implemented as described by Edwards and Palsson (2000a). To implement both the linear FBA constraints and the nonlinear EBA constraints, we use the nonlinear constrained optimization package available in the MATLAB optimization package.

Predicted fluxes obtained using FBA, and FBA combined with EBA, are presented in Table 1. All fluxes are reported in units of mmol sec¹. The predicted optimal fluxes obtained with FBA alone are thermodynamically infeasible, because no feasible potential vector exists for these fluxes. When EBA and FBA are applied simultaneously, thermodynamically feasible mass-balanced fluxes are predicted. Because we consider the EBA constraints in addition to the FBA mass balance constraints, EBA reduces the size of the feasible flux space. However, we do find more than one optimal solution when we maximize the production of D. By varying the initial guess, the optimization routine arrives at different feasible solutions that maximize production. Therefore the fluxes reported in Table 1 do not represent unique solutions to the optimization problem, either for the FBA cases or for the FBA/EBA cases.

	ANALYSIS OF ESCHERICHIA COLI METABOLISM

TOP ABSTRACT INTRODUCTION AND BACKGROUND THEORY OF ENERGY BALANCE... NONEQUILIBRIUM STEADY-STATE... APPLICATION OF ENERGY BALANCE... ANALYSIS OF ESCHERICHIA COLI... ESTIMATION OF REACTION... DISCUSSION AND CONCLUSIONS REFERENCES

We have obtained the stoichiometric matrix provided by Edwards and Palsson (2000c), used to represent the flux balances in the E. coli metabolism. The web-posted supplementary information provides detailed descriptions of the reaction network, which contains 953 fluxes (735 internal; 218 boundary) acting on 536 metabolites. Again using the MATLAB (The Mathworks) optimization package, we reproduced the linear programming analysis presented in Edwards and Palsson (2000a), and optimized the production of biomass with glucose and oxygen uptake constrained to be less than or equal to 10 and 15 mmol g-DW¹h¹, respectively. The resulting flux produces a growth rate of 0.81 h¹ on glucose minimal media.

To compute the thermodynamic properties of this large-scale network, we first combine redundant fluxes (e.g., the phosphofructokinase A and B reactions are combined into a single column of S'), resulting in n_c = 617 internal reactions operating on 536 metabolites. The growth is then optimized under the flux balance constraints (Eqs. 1 and 2) and the constraint that the energy-balance equations (Eqs. 7 and 8) are satisfied. We impose the additional constraint that µ must be finite for all nonzero fluxes: the free energy can go to zero only if the flux is zero, implying that a given reaction is in chemical equilibrium. This analysis predicts the same optimal growth rate (0.81 h¹ on glucose minimal media) as that reported above for FBA. Yet, as for the simple example illustrated above, the FBA prediction does not represent a unique solution to the optimization problem because redundancies in the metabolic network allow the optimal growth rate to occur for an infinite number of possible internal flux distributions. Because the EBA-constrained solution is (at least in this case) able to achieve the same optimal growth rate as that obtained by considering only the FBA constraints, the fluxes obtained by EBA represent one particular optimal solution to the FBA linear programming problem. However, optimal flux distributions obtained by considering flux balance alone are not necessarily thermodynamically feasible. The EBA constraints further restrict the set of feasible fluxes, and provide a more physically realistic flux distribution. Selected fluxes from glycolysis, TCA cycle, and respiration are tabulated in Fig. 4. Fluxes from the wild type (WT) on glucose media under aerobic and anaerobic conditions are labeled WT (oxygen) and WT (anaerobic), respectively.

View larger version (89K): [in this window] [in a new window]   FIGURE 4   Regulation of reactions in glycolysis, TCA cycle, and respiration. Predicted biochemical fluxes (in µmol g-DW1h1kcal1), chemical potentials (in kcal mol1), and conductances (in 103 g-DW1h1kcal1) are reported for 64 central metabolic reactions, for WT aerobic and anaerobic growth, and for zwf, pgl, gnd, pyk, pgi, eno, and pfk knockouts. Predicted growth rate for each case is reported in units of hour one. Green and blue shading indicate up and down regulation relative to WT, respectively. Yellow indicates that a reversible flux has changed direction. Gray shading indicates that the flux cannot be identified as a site of regulation because, although the conductance has changed moderately (by less than a factor of 10), the magnitudes of the potential and the flux have both either increased or decreased. "N.P." denotes enzymes that are not present under given conditions. "Eq." refers to reactions that are at equilibrium, for which the conductance cannot be estimated. All simulations are performed for glucose-minimal media.

	ESTIMATION OF REACTION POTENTIALS

TOP ABSTRACT INTRODUCTION AND BACKGROUND THEORY OF ENERGY BALANCE... NONEQUILIBRIUM STEADY-STATE... APPLICATION OF ENERGY BALANCE... ANALYSIS OF ESCHERICHIA COLI... ESTIMATION OF REACTION... DISCUSSION AND CONCLUSIONS REFERENCES

Beyond providing the constraints necessary to produce thermodynamically feasible fluxes, EBA may make it possible to make quantitative estimates of the chemical potentials in the system. As a first attempt to approximately identify the chemical potentials, we introduce the "flux resistances," defined as r_i = µ_i/J_i. To avoid the unphysical situation in which µ_i = 0 for a finite J_i, which is equivalent to setting the flux resistance equal to zero, we assume that there exists a minimum flux resistance, r_min (which is equivalent to saying that there exists an upper limit to the effective reaction rate constants). Thus, the second law constraint is modified:

(9)

Realistically, each flux may have a different r_min. However, for our purposes, we find that a single value, r_min = 400 kcal mol²g-DW h, produces reasonable behavior and can be used to describe the entire network. The energy balance constraint can alternatively be written in terms of the chemical potentials (Eq. 7) or the flux resistances,

(10)

To estimate the chemical potentials, we assume that the biological system tends to operate in a steady-state regime where the magnitudes of the chemical potentials are minimized. With the fluxes fixed using the values obtained from the above analysis, we use quadratic programming to find an optimal µ that minimizes the norm of the chemical potential vector |µ|². In addition to the fluxes, Fig. 4 lists the predicted chemical potentials and the conductances, c_i = r, of each reaction. The choice of r_min = 400 kcal mol²g-DW h results in reasonable predictions of chemical potential. For example, µ_i = 9.55 kcal mol¹ for ATP hydrolysis. Small changes to the assumed value of r_min result in proportional changes in all of the predicted reaction potentials, while the relative values remain fixed.

The reaction conductance provides a measure of the activation level of the pathway. If c_i = 0, then the associated enzyme(s) is not present or is deactivated. Increases in flux or conductance, at a fixed chemical potential, indicate that a pathway is up regulated at either the expression level or the post-translational level. By changing the boundary constraints, we can study how the metabolic network responds to changes in substrate. In Fig. 4, we compare the predicted EBA fluxes and potentials for the WT cell grown under aerobic and anaerobic conditions. We identify a pathway as "down regulated" (colored blue) if the flux conductance decreases by a factor of 4 or more, and "up regulated" (colored green) when the conductance increases by more than a factor of 4. On the basis of this analysis, we identify three enzymes that require activation upon moving from aerobic to anaerobic conditions: fumarate reductace, pyruvate formate lysase, and pyridine nucleotide transhydrogenase. Other pathways show increases or decreases in flux compared to aerobic growth. However, these differences can be accounted for by changes in the potential and thus do not necessarily require regulation.

Following the work of Edwards and Palsson (2000a, b, c), we studied the effects of gene knockouts on the metabolic capabilities of E. coli. We found that zwf, pgl, and gnd knockouts result in the same predicted phenotype (Fig. 4), with only two up regulations compared to WT:succinate dehydrogenase and pyridine nucleotide transhydrogenase. Again, a number of predicted fluxes that differ from the WT do not require up regulation of the associated enzymes. The situation is similar for pyk and pgi knockouts. The pyk knockout requires an up regulation of phosphoenolpyruvate carboxylase to maintain growth of almost 99% that predicted for WT. No other significant regulations are predicted. The pgi knockout analysis predicts one significant up regulation and three down regulations and a similar growth rate. Thus, we find that much of the capacity for metabolic control is built in to the WT expression and activation levels. Few regulatory steps are required when nonessential genes are knocked out; the network is robust and tolerant to errors and deletions (Jeong et al., 2000).

However, the situation is different when essential genes are knocked out. The eno and pfk genes are examples of genes that FBA falsely identifies as nonessential to growth on glucose minimal media (Edwards and Palsson, 2000a). The fact that eno and pfk are essential for growth in vivo (Fraenkel, 1996) may be related to the demands that the knockouts of these genes place on metabolic regulatory mechanisms. These demands are much heavier than those imposed by nonessential knockouts. Our analysis predicts that maintaining growth with these knockouts requires a greater number of pathway regulations than for the nonessential knockouts. Of the 64 reactions listed in Fig. 4, 15 are predicted to be up regulated and 15 down regulated for the pfk knockout. Thus, the predicted phenotype is markedly different from the WT for nearly half of the reactions associated with glycolysis, TCA cycle, and respiration, whereas the pyk and pgi knockouts differ very little from WT.

	DISCUSSION AND CONCLUSIONS

TOP ABSTRACT INTRODUCTION AND BACKGROUND THEORY OF ENERGY BALANCE... NONEQUILIBRIUM STEADY-STATE... APPLICATION OF ENERGY BALANCE... ANALYSIS OF ESCHERICHIA COLI... ESTIMATION OF REACTION... DISCUSSION AND CONCLUSIONS REFERENCES

The predictions of reaction potentials and conductivities presented above are based on one major simplificationthat the entire network can be characterized by a single r_min, or equivalently a maximum pathway conductance. More realistically, constraints could be placed on each pathway based on a priori knowledge of the biochemistry. One promising aspect of EBA is that it is possible to incorporate as much, or as little, as is known about the individual reactions. For example, if we know that the ratio [ATP]/[ADP] in the cell is greater than the equilibrium ratio (Stryer, 1995), then we know that the free energy of ATP hydrolysis satisfies the inequality

(11)

In general, consider the reaction AB. If the ratio [A]/[B] is greater than 1, then the potential energy of that reaction is constrained µ_i < µ. Alternatively, if the concentration ratio is less than 1, then µ_i > µ. If the reactant concentrations are known, then the constraint becomes an equality: µ_i = k_BT ln(K_eq[A]/[B]), where K_eq is the equilibrium constant of the reaction.

The arbitrary parameter r_min was introduced to estimate reaction potentials. However, the EBA constraints used to determine the thermodynamic feasibility of a network flux do not depend on the choice of r_min. The feasible flux space allowed by EBA and FBA together is more restricted than the space predicted by FBA alone. Therefore, predictions made without considering energy balance in general do not satisfy the EBA constraints and are not physically realistic.

Together, flux balance and energy balance provide basic laws for the analysis of biochemical networks. These laws, akin to the basic principles of circuit theory for electrical networks, make the analysis and design of large-scale biochemical systems practical. The engineering approach to analysis and design of such complex systems is the identification of modular components that are separable within the system (Hartwell et al., 1999). Flux balance and energy balance analysis provide a basis for understanding how these modules function and interact.