1Biotechnology and Bioengineering Center, Department of Physiology, Medical College of Wisconsin, Milwaukee, Wisconsin; and 2Department of Applied Mathematics, University of Washington, Seattle, Washington
Submitted 7 June 2004 ; accepted in final form 20 October 2004
![]() |
ABSTRACT |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
network thermodynamics; nonequilibrium steady state; biochemical network
In this work, physicochemical constraints on metabolite concentrations are incorporated into FBA and EBA, and these tools are applied to understanding metabolic regulation in the hepatocyte. Application of these constraints, which arise naturally from reaction thermodynamics, provide modeling results that 1) are more physically realistic and biologically meaningful than are provided by FBA alone, 2) facilitate direct comparison of computational predictions and experimental results, and 3) reveal insight into the regulatory and control mechanisms operating in cells. In this work, the developed tools are applied in silico to profiling of glycogenesis and glycogenolysis and the action of these pathways in certain metabolic diseases.
![]() |
BIOCHEMICAL THERMODYNAMICS |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
![]() | (1) |
To characterize the thermodynamics of reactions such as Eq. 1, transformed thermodynamic quantities G' and K' which refer to the apparent Gibbs free energy and equilibrium constant, respectively, are used. The apparent thermodynamic constants operate at constant pH; thus chemical reactions written in terms of sums of species in general do not stoichiometrically conserve hydrogen atoms. Thermodynamic relationships for apparent
G' and K' are analogous to those for
G and K, where
G' and K' are expressed in terms of reactants rather than species. Therefore,
![]() | (2) |
The stoichiometric thermodynamic constraints on fluxes that have been introduced in previous works (8, 9, 28) apply equally well to biochemical reactions expressed using species and those using reactants. For biochemical applications, it is often convenient to work in terms of reactants rather than species.
![]() |
METHODS AND RESULTS |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
The central flux balance conservation statement is given by the equation
![]() | (3) |
To enforce thermodynamic balance of free energies in the reaction network, we use the null space decomposition of the stoichiometric matrix, which 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 a reaction loop (8, 9, 28). The free energies are balanced by the equation
![]() | (4) |
The second law of thermodynamics takes the form of an inequality constraint for each flux:
![]() | (5) |
Identification of Hepatic Cell Fluxes
On the basis of the FBA/EBA constraints defined in Eqs. 35, we next consider the metabolic fluxes in the rat liver during glycogenic and glycogenolytic operational modes, as illustrated in Fig. 1. The complete system of reactions is listed in Table 1; abbreviations for metabolic species are listed in Table 2. In addition to glycogen synthase and glycogen phosphorylase reactions, the reactions considered include those of glycolysis, tricarboxylic acid (TCA) cycle, oxidative phosphorylation, and various transporters. The notation "_m" is used to denote a species present in the mitochondrial matrix. Thus "ATP" and "ATP_m" denote cytosolic and mitochondrial ATP, respectively.
|
|
|
From available measurements of rates of oxygen consumption, ATP utilization, and glucose uptake and output during glycogenic and glycogenolytic stimulation in the perfused rat liver (12, 16, 34), the reaction fluxes through the metabolic network can be quantified using the FBA and EBA constraints. For these cases, the fluxes through each reaction are uniquely identified on the basis of the experimental measures of the boundary fluxes and the FBA/EBA constraints. The flux measurements used to calculate the steady-state fluxes for the two modes are specified as follows.
Glycogenic mode (mode 1).
The rate of nonglycogenic ATP use (ATP hydrolysis for other cell functions) is set to 1.44 µmol·min1·g wet wt1 (12), and the rate of glucose uptake was set to 1.5 µmol·min1·g wet wt1 (31). Given these boundary fluxes, in addition to the mass-balance and thermodynamic constraints of Eqs. 35, the fluxes throughout the network are uniquely identified and plotted in Fig. 1 and tabulated in Table 3. The constraint-based model predicts that 93.7% of the glucose taken up is converted to glycogen. The remaining glucose is used to drive glycolysis and oxidative phosphorylation, providing the ATP source for glycogenesis and other cell functions.
|
The flux distributions for mode 1 and mode 2, as illustrated in Fig. 1, are uniquely determined by the constraints imposed by the input experimental data and the FBA and EBA laws; thus no optimality condition is necessary to determine the flux values. Note that, in addition to differences in the glycogen fluxes, the rates of glycolysis and oxidative phosphorylation are greater in the glycogenic mode than in the glycogenolytic mode due to the hydrolysis of two ATP molecules per glucose converted to glycogen.
Thermodynamic Constraints on Concentrations
To further improve the accuracy of constraint-based modeling, thermodynamic constraints on the metabolite concentrations are applied, which in turn constrain the allowable reaction free energies and fluxes and improve the accuracy of model predictions.
The free energies for a biochemical network can be expressed as a function of reactant concentrations in matrix-vector form:
![]() | (6) |
![]() | (7) |
For typical large systems such as Escherichia coli metabolism, where the number of reactions is larger than the number of reactants (18), Eq. 7 is overdetermined; that is, T has more rows than columns. The requirement that
satisfies the energy balance loop law (Eq. 4) ensures that concentrations that satisfy Eq. 7 exist. In addition, if certain concentrations are known a priori, or if constraints are enforced on the feasible concentration values, these constraints act through Eq. 7 on the feasible free energies.
Strictly speaking, Eqs. 6 and 7 should be formulated in terms of metabolite activities, rather than concentrations, as the intracellular environment is not an ideal dilute solution (10). Because the nonideal nature of intracellular compartments is not well characterized, the present analysis is formulated in terms of concentrations, as is common in the literature. In the limit that activity and concentration are proportionally related, concentration variables in the present analysis can be replaced by metabolite activities without changing any of the results.
In addition to Eqs. 6 and 7, a reaction network's stoichiometry reveals a set of constraints on certain conserved concentration pools, as introduced by Alberty (1, 2). These constraints follow from the equation for kinetic evolution of the metabolite concentration vector:
![]() | (8) |
![]() | (9) |
(Note that the matrix N contains a basis for the right null space of .) Thus, if we consider two different steady metabolic states, given a complete set of allowable reactions and transport fluxes included in S, the feasible concentrations will be constrained, so that
![]() | (10) |
Generalized Flux-Potential Relationship
To quantify enzyme activity in terms of reactant concentrations and reaction flux and free energy, consider a general uni-unimolecular reaction:
![]() | (11) |
The concentration ratio [B]/[A] can be expressed in terms of the transformed free energy:
![]() | (12) |
When the reaction kinetics is governed by simple mass action, the flux can be expressed as (28)
![]() | (13) |
![]() | (14) |
![]() | (15) |
![]() | (16) |
For reactions involving multiple species with the general form
![]() | (17) |
![]() | (18) |
In Eqs. 17 and 18, xi represents the ith metabolite, and µi and vi are the stoichiometric coefficients associated with xi on the left- and righthand sides the reaction.
In the sections below, we make use of the following definition for enzyme activity, X:
![]() | (19) |
For the Michaelis-Menten flux, activity is proportional to enzyme concentration. The behavior of as a function of
, predicted by the generalized potential-flux relationship of Eq. 19, is plotted in Fig. 2. Notice that, in the limit that
is large,
. Thus reaction fluxes saturate at a level that is proportional to enzyme activity.
|
Values for the transformed standard free energies of formation, , for the reactants included in the hepatic cell model are tabulated in Chapter 4 of Ref. 5, with the exception of values for
, and
Values for these five thermodynamic parameters were obtained as follows.
Assuming that the free energies of hydrolysis of UTP and GTP are similar to that of ATP, and
are set equal to
;
and
are set equal to
The free energy of formation of UDP-glucose, , is set so that the standard free energy of the udp reaction is
kcal/mol, based on a measure of
at
(21).
As it appears in the reaction stoichiometry (Table 1), GLYCOGEN represents the formation of a glycosidic bond via the reaction . The parameter
is set to 59.2 kcal/mol, which results in the standard free energy for the glycosidic bond formation reaction,
, of +5.5 kcal/mol.
The complete set of values for the reaction network is listed in Table 2. These values correspond to a parameterization of the standard transformed free energies of formation at fixed pH of 7. By definition, hydrogen ion in bulk solution has a transformed standard free energy of formation of
. The parameter
corresponds to the free energy for H+ transport across the mitochondrial membrane. Assuming that the mitochondrial membrane potential is 160 mV, the free energy is set to
kcal/mol. Flux through adenine nucleotide translocase involves the net movement of charge across the inner mitochondrial membrane. Thus the free energy for this flux includes a contribution from the electrostatic membrane potential and is computed
kcal/mol. According to this relationship, the net transport of an energetic phosphate from matrix to cytosol may be thermodynamically favored, whereas the ratio of ATP to ADP is smaller in the mitochondrial matrix than in the cytosol.
Additional Constraints on Concentrations
We assume that the variable metabolite concentrations do not attain values greater than 10 mM. In addition to the flux balance and thermodynamic constraints, the following constraints on the metabolite concentrations in the network are set.
Glucose concentration is fixed at 5 mM.
Total CO2 (aqueous carbon dioxide plus carbonic acid) is fixed at 25 mM.
Oxygen concentration is fixed at 27 µM.
The nucleoside diphosphokinase reactions (ndk1 and ndk2) are assumed to operate close to equilibrium, so that kcal/mol.
To ensure that cellular processes driven by the phosphate potential can operate normally, the following constraint is applied: kcal/mol.
These five constraints, combined with Eq. 5 (the Second Law of Thermodynamics) and Eq. 10 (the conservation of metabolic pools), represent the complete set of constraints used to define the feasible space for mode 1 and
mode 2, the metabolic concentrations in the glycogenic and glycogenolytic metabolic modes. Constraints 4 and 5 in the above list, which are expressed in terms of reaction free energies, are translated to constraints on reactant concentrations via Eq. 6.
Computational Profiling of Enzyme Activities
Although the metabolic fluxes in the glycogenic and glycogenolytic modes are uniquely identified based on the experimentally measured boundary fluxes, the sets of metabolite concentrations and enzyme activities that are feasible in these metabolic modes are not unique. In this section, criteria are introduced for determining metabolite concentrations and enzyme activities, given the constraints outlined above.
For a given feasible concentration vector, the enzyme activities are computed on the basis of the definition Eq. 19:
![]() | (20) |
Xj is the activity of the jth enzyme, Jj is the flux of the jth reaction, and Ct,j() and
G'j(
) are functions of the metabolite concentrations computed according to Eqs. 19 and 6.
Differences in activities between the two modes may be quantified by formulating a hypothetical objective function for reaction free energies (or, equivalently, reactant concentrations). As an example, we first explore the possibility that differences in reactant concentrations between the two operational modes are minimized. This hypothesis is formalized by minimizing the objective function
![]() | (21) |
The free energies and activities resulting from minimizing f1 are illustrated in Fig. 3. In Fig. 3, top, the predicted reaction free energies for each operating mode are plotted; the bottom plot illustrates the differences in predicted activities between the two modes. Minimization of the objective f1 of Eq. 21 results in free energies that are nearly identical between the two modes, because the differences between concentrations in the two modes are minimized. However, this analysis reveals that multiple up- and downregulations are required to achieve feasible operation with the concentration differences minimized. The required activations/deactivations manifested as differences in enzyme activities between the two modes, as illustrated in Fig. 3, bottom. For the thermodynamically feasible chemical driving forces plotted in Fig. 3, the majority of enzymes (36 of 41 internal reactions) are activated during glycogenesis, whereas two are upregulated during glycogenolysis.
|
![]() | (22) |
The denominator in Eq. 22 normalizes f2 so that it measures relative differences between reaction activities in the two metabolic modes.
The results obtained using f2 as an objective function are illustrated in Fig. 4. In this case, we obtain a different set of reaction free energies for each mode, with a relatively small number of enzymes regulated between the two modes. Although the fluxes in the two modes are different for each reaction in the system, changes in fluxes between modes 1 and 2 correspond to equivalent changes in driving forces for many of the reactions. Thus the activities remain constant for the majority of enzymes.
|
The enzymes identified in Fig. 4, bottom, represent a minimal set of enzymes that must be regulated in switching between the two metabolic modes. In other words, Fig. 4 represents the minimal scheme, in terms of number of enzymes regulated, in controlling switching between the two metabolic modes. We identify these enzymes, for which activity differs between the two modes (Fig. 4, bottom), as putative regulatory sites in the hepatic cell network.
These results are promising, because our predictions match the regulatory structure that has been experimentally determined for hepatic glucose/glycogen metabolism. Major sites activated by glucagon are glucose-6-phosphatase (g6pt) and glycogen phosphorylase (glp) (22, 24, 30), whereas the primary sites of insulin activation are glucokinase (glk) and glycogen synthase (gys) (20, 30). Validation of the identified upregulation of udp and pph during glycogenesis requires further investigation.
This successful identification of key regulatory sites supports the hypothesis that biological systems operate at or near optimal controllability (with a minimal number of regulatory sites in the network) within a given set of constraints. This hypothesis requires further testing in the context of both small-scale and large-scale systems.
In addition to the identification of the cellular regulatory scheme, the optimal solution (with controllability as the objective) provides estimates for reactant concentrations in hepatocytes in each metabolic mode. For the concentrations estimated by optimizing the objective defined in Eq. 22, computed free energies are tabulated in Table 3; computed metabolite concentrations are tabulated in Table 2.
The predicted concentrations for the glycogenic and glycogenolytic modes are compared in Fig. 5. Filled circles correspond to the concentrations computed with the objective of Eq. 21, which minimized the differences in the concentrations. Open circles correspond to the concentrations of the minimal control scheme computed via the objective of Eq. 22. Note that there is a positive correlation between the metabolite concentrations in each mode, even for concentrations associated with the minimal control scheme. Three major outliers, PPi, UDPG, and G1P are indicated in the figure. The computed control scheme predicts that the concentrations of these metabolites, which are important intermediates in glycogen synthesis, increase sharply in transition from glycogenic to glycogenolytic operation. In particular, the PPi concentration is predicted to increase three orders of magnitude.
|
The optimality condition used above allows us to compute a unique state (fluxes, concentrations, activities) in each mode. The corresponding activity values can be used to explore the consequences of impaired (or enhanced) activity at specific locations in the network. As an example of the type of analysis that is possible, we have studied a pair of metabolic diseases that impact glycogenolysis, von Gierke's disease and Hers' disease, which are characterized by malfunctioning glucose-6-phosphatase and glycogen phosphorylase, respectively (30). Here, a malfunctioning enzyme is modeled as a reduction in the enzyme's activity, relative to the normal case. All other enzymes are treated as normal, with activities based on the glycogenolytic mode (mode 2).
The behavior of the integrated network is simulated using Eq. 8, where the fluxes are computed on the basis of Eq. 19. Given a set of activities, Eq. 8 is integrated, holding the fixed concentrations constant (glucose, CO2, and oxygen), until a new steady state is reached. The impacts of the malfunctioning enzymes in von Gierke's and Hers' disease are quantified by comparing the predicted rate of glucose production to that of the normal state.
The plot in Fig. 6, bottom, illustrates that, as enzyme activity is diminished, the predicted glucose production in the glycogenolytic mode is reduced. The model predicts that glycogenolytic operation is significantly more sensitive to impairment of glucose-6-phosphatase than to impairment of glycogen phosphorylase. Figure 6 shows that, at 10% of normal function of glycogen phosphorylase (Hers' Disease), steady-state glucose production in the glycogenolytic mode is greater than 90% of that in normal function. The sensitivity of glycogenolytic operation to the impairment of g6pt is much greater than that of impairment of glp. The simulations predict that, at 10% of normal function of g6pt, glucose production is effectively shut down, with glucose production less than 5% of normal.
|
![]() |
DISCUSSION |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
In summary, we have introduced thermodynamic-based constraints on biochemical fluxes and concentration to augment the FBA modeling approach. The methodology allows us to make predictions not available from FBA alone, namely, predictions of reactant concentrations, reaction potentials, and enzyme activities. With these predictive capabilities, new applications of constraint-based modeling, such as reconstruction of cellular regulatory mechanisms and the profiling of metabolic diseases, are introduced.
The penalty that is paid for introducing greater functionality into constraint-based modeling is the requirement that thermodynamic data on biochemical reactions be obtained. These data are available in the form of physicochemical constants (albeit available in varying degrees of accuracy and consistency); they are not a set of model- and cell-dependent parameters (e.g., kinetic constants) that must be identified from physiological experiments. We anticipate that, in general, physicochemical constants may be identified with greater confidence than other model parameters.
Although these thermodynamic data are required for the estimation of metabolic concentrations and reaction free energies presented here, they need not be known for application of the EBA constraints of Eqs. 4 and 5, which constrain the feasible network fluxes. When Eqs. 4 and 5 are applied, an ad hoc set of irreversibility constraints, often used to constrain network behavior in flux balance analysis, is not needed. These EBA constraints may be imposed in constraint-based network analysis even when the equilibrium thermodynamic data are not known, because the EBA constraints on reaction directions require knowledge of only the stoichiometric matrix (8, 9, 28). However, to incorporate constraints on concentrations and to estimate concentrations and steady-state free energies, the thermodynamic data are required.
Model Assumptions
The biochemical network analyzed in this paper does not represent a complete set of all reactions coupled to glucose and glycogen metabolism. In fact, through network interactions no components of metabolism are uncoupled from any others, and no computational model is exhaustive in scope. In this case in particular, introducing knowledge of the reaction directions of the gluconeogenic pathway would provide further constraints on the feasible concentrations of glucose, lactate, ATP, ADP, and other key metabolites. It is possible that, in the experiment used to identify the glycogenolytic mode (16), gluconeogenesis from endogenous sources was responsible for a portion of the measured glucose production. In addition, the model assumes that there is no cycling between glucose and glycogen in either operating mode. In other words, the fluxes through g6p and glp are zero in mode 1, and the fluxes through udp1, pph, gys, and ndk1 are zero in mode 2. More realistically, is expected that there exists a small, finite flux through all six of these reactions, which combine for a net reaction stoichiometry of G6P + ATP + 2 H2O = Glc + ADP + 2 Pi. Thus the network and the associated boundary fluxes used above may be incomplete, resulting in inaccurate estimates of the network fluxes in this mode.
However, our results on identifying regulated enzymes are sensitive mainly to overall glucose flux (net uptake vs. net uptake) in the two operating modes, and not particularly to the magnitudes, of fluxes in the network. We find that the estimated flux magnitudes can be varied by as much as 50% without changing the results presented in Figs. 35. Of course, explicitly including additional reactions would increase the number of activities considered in the objective function for computing the minimal set of controlled enzymes and perhaps change the resulting predictions.
The major physicochemical assumption used in this work is the assumption that enzyme binding sites are not saturated with reactants in our development of a generalized potential-flux relationship; i.e., metabolite concentrations are assumed to be lower than the Michaelis-Menten constants for enzymes in the network. Although this assumption is not expected to hold for all metabolic reactions, it is a useful assumption to begin with when the Michaelis-Menten constants are not known. When a given enzyme is not saturated, the generalized potential-flux relationship, Eq. 19, involves only one unknown kinetic parameter, X. Only this one parameter per reaction is necessary, because reaction-specific thermodynamic information has been incorporated into G'0.
When in vivo values for Michaelis-Menten constants are known for certain enzymes, Eq. 19 can be modified to include saturation effects:
![]() | (23) |
In studying glycogenesis and glycogenolysis, in addition to the simplified potential-flux relationship, simplifying model assumptions were introduced into the reaction network (Table 1). In particular, the electron transport chain was collapsed to include a single electron acceptor Q; the individual steps (complexes I, II, III, and IV) are not explicitly considered. With this simplification, the overall stoichiometry is correct, although the thermodynamic constraints associated with oxidative phosphorylation may be quantitatively different when the intermediate steps are incorporated. In future work, we will introduce a more detailed representation of oxidative phosphorylation, including proton leak across the inner mitochondrial membrane (23) and a detailed model of the electron transport system. In addition, we will expand the model to include important components of energy metabolism, such as gluconeogenesis and fatty acid synthesis and oxidation.
A further assumption made in the present analysis is that each cellular compartment (mitochondrial, cytoplasmic) is treated as well mixed; that is, the potential for spatial gradients in concentration is ignored. It is well known that, for species that are rapidly consumed, such as oxygen, significant spatial gradients in concentration can be established (7). As the concentration of oxygen varies with intracellular location, mitochondrial metabolism may vary as well, although drastic changes in mitochondrial fluxes are not observed until the oxygen level drops below 1 mmHg in skeletal and cardiac muscle (11, 23, 33).
Possible Applications
Because the thermodynamic-based network analysis approach introduced here allows the prediction of reactant concentrations, it will facilitate the incorporation of experimental measures directly into mathematical modeling predictions. This application requires a modification of the traditional approach to data fitting and parameter estimation. Instead of varying parameter values used in a deterministic kinetic model to minimize the difference between measured and predicted values, the measured values may be used directly as additional constraints in the constraints-based model.
In typical applications, the concentration data may come from repeated measurements obtained over a number of experiments. One possible approach to imposing experimentally obtained concentration constraints will be to translate statistical distributions of measured values into representative ranges of allowable concentrations to apply as model constraints.
The advantage of this approach is that the resulting models will be automatically experimentally validated, since the available experimental data will be explicitly incorporated into the model design. The disadvantage will be that, in many cases, the available experimental data may not provide enough information to uniquely constrain the models. However, the introduction of explicit hypotheses, such as those introduced here, may provide unique predictions.
In addition to profiling metabolic diseases, the approach outlined here for computationally profiling the influence of individual malfunctioning enzymes can be applied to systematically probe possible drug targets and predict the effects of specific inhibitors (or activators) on the integrated network behavior. Applications to target identification and/or lead optimization will require 1) building specific disease models or sets of models for given diseases, 2) determining a number of effected targets in more than one tissue type, and 3) quantifying the disease state(s) in terms of enzyme activities. With a determined model in hand, profiling the disease state will be accomplished by simulating the integrated system behavior and comparing the model predictions to the healthy case, a generalization of the methodology introduced here for phenotyping metabolic diseases that act on single targets. Downstream targets affected by therapeutic agents may be identified using an approach analogous to that used to identify the regulatory structure in the hepatic cell network.
![]() |
GRANTS |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
![]() |
FOOTNOTES |
---|
The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked "advertisement" in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
![]() |
REFERENCES |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
|
HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
Visit Other APS Journals Online |