Thermodynamic-based computational profiling of cellular regulatory control in hepatocyte metabolism

Daniel A. Beard1 and Hong Qian2

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
 TOP
 ABSTRACT
 BIOCHEMICAL THERMODYNAMICS
 METHODS AND RESULTS
 DISCUSSION
 GRANTS
 REFERENCES
 
Thermodynamic-based constraints on biochemical fluxes and concentrations are applied in concert with mass balance of fluxes in glycogenesis and glycogenolysis in a model of hepatic cell metabolism. Constraint-based modeling methods that facilitate predictions of reactant concentrations, reaction potentials, and enzyme activities are introduced to identify putative regulatory and control sites in biological networks by computing the minimal control scheme necessary to switch between metabolic modes. Computational predictions of control sites in glycogenic and glycogenolytic operational modes in the hepatocyte network compare favorably with known regulatory mechanisms. The developed hepatic metabolic model is used to computationally analyze the impairment of glucose production in von Gierke's and Hers' diseases, two metabolic diseases impacting glycogen metabolism. The computational methodology introduced here can be generalized to identify downstream targets of agonists, to systematically probe possible drug targets, and to predict the effects of specific inhibitors (or activators) on integrated network function.

network thermodynamics; nonequilibrium steady state; biochemical network


CONSTRAINT-BASED APPROACHES to modeling biochemical systems provide powerful computational tools for predicting the metabolic capabilities of whole genome cell systems (14, 15, 1719). Yet, the most widely used computational method for large-scale systems, flux-balance analysis (FBA), typically cannot provide precise predictions of biochemical fluxes. Rather, FBA uses the known network stoichiometry to impose a steady-state mass balance constraint that yields a feasible flux space within which an infinite number of equivalent solutions can exist (13, 29, 32). The recently introduced theory of energy-balance analysis (EBA) augments FBA by characterizing the thermodynamic constraints on reaction-free energies that arise from the network stoichiometry (8, 9, 28). This theory provides the framework for introducing thermodynamic constraints into biochemical network modeling, but still does not identify enough constraints to provide reliable estimates of free energies, concentrations, and fluxes for a whole cell system, a problem of intense current interest.

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
 TOP
 ABSTRACT
 BIOCHEMICAL THERMODYNAMICS
 METHODS AND RESULTS
 DISCUSSION
 GRANTS
 REFERENCES
 
In the treatment of biochemical thermodynamics in this paper, it is useful to formally distinguish reactants from species in chemical reactions (1–6). By use of this formalism, reactants are considered to be made up of sums of particular species. For example, the reactant ATP exists as several distinct species, ATP4–, HATP3–, and H2ATP2–, which are assumed to rapidly interconvert. Thus, using these definitions, the chemical reaction

(1)
is written in terms of biochemical reactants, each composed of sums of several species.

To characterize the thermodynamics of reactions such as Eq. 1, transformed thermodynamic quantities {Delta}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 {Delta}G' and K' are analogous to those for {Delta}G and K, where {Delta}G' and K' are expressed in terms of reactants rather than species. Therefore,

(2)
where R is the gas constant, and T is the absolute temperature. The activity of water is defined to be unity. In Eq. 2, the term [ATP], for example, denotes the summed concentration of all ATP species: [ATP] = [ATP4–] + [HATP3–] + [H2ATP2–]. Alberty (5) details the approach to biochemical thermodynamics presented in terms of sums of species in the recently published textbook.

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
 TOP
 ABSTRACT
 BIOCHEMICAL THERMODYNAMICS
 METHODS AND RESULTS
 DISCUSSION
 GRANTS
 REFERENCES
 
Flux Balance and Energy Balance Constraints on Fluxes

The central flux balance conservation statement is given by the equation

(3)
where {nfr_n} is the vector of n fluxes occurring in the network, S {nfr_mn} is the stoichiometric matrix, and m is the number of reactants in the system (notation is used throughout to distinguish a vector).

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)
where the matrix N is made up of vectors that form a basis of the null space of , and is the stoichiometric matrix of internal reactions (excluding boundary fluxes). The vector {Delta}' lists the free energies associated with each of the reaction fluxes. Vectors in the null space of can be thought of as internal loops or summed reactions for which there is no net change in the participating reactants. Thus Eq. 4 is analogous to Kirchoff's loop law: the potential drop summed over closed loop is zero.

The second law of thermodynamics takes the form of an inequality constraint for each flux:

(5)
which ensures that entropy production is positive for each reaction. Equations 4 and 5 represent thermodynamic constraints that must be considered in addition to the flux balance constraints for the analysis to be physically meaningful. The second law provides a link that couples mass balance and energy balance and constrains the feasible flux space more strictly than Eq. 3 alone. Predicted fluxes obtained with FBA alone can be thermodynamically infeasible, because no feasible free energy vector exists for these fluxes. EBA provides the equivalent of Kirchoff's voltage law for biochemical networks and facilitates the estimation of chemical potentials, as well as prediction of thermodynamically feasible fluxes, for biochemical systems.

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.



View larger version (27K):
[in this window]
[in a new window]
 
Fig. 1. Hepatic fluxes in glycogenic and glycogenolytic operational models. Top: reactions in glycogenesis and glycogenolysis are illustrated. For simplicity, the mitochondrial reactions of the tricarboxylic acid cycle and oxidative phosphorylation are not shown. As illustrated, reactions in the cytosol and in the mitochondrial space are coupled through the phosphate and redox potentials. Bottom: metabolic fluxes in glycogenic (mode 1) and glycogenolytic (mode 2) operating modes are plotted for all reactions tabulated in Table 1.

 

View this table:
[in this window]
[in a new window]
 
Table 1. Reactions of the hepatic cell model

 

View this table:
[in this window]
[in a new window]
 
Table 2. Values of transformed free energies of formation and predicted concentrations for metabolites in the hepatic cell model

 
In Table 1, certain reactions (ace and suc1) are expressed in a form slightly different from the standard stoichiometry found in most biochemistry texts. Water molecules appear on the lefthand side of the reactions because CO2 is treated as the total concentration of CO2 in the aqueous phase, including H2CO3, HCO3, and CO formed via the hydration of CO2 by carbonic anhydrase. Correct treatment of the transformed thermodynamic free energies and chemical equilibria for these reactions requires that H2O appear as a reactant (4, 5).

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·min–1·g wet wt–1 (12), and the rate of glucose uptake was set to 1.5 µmol·min–1·g wet wt–1 (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.


View this table:
[in this window]
[in a new window]
 
Table 3. Fluxes and predicted free energies in the hepatic cell model

 
Glycogenolytic mode (mode 2). As for the glycogenic mode, the rate of nonglycogenic ATP use (ATP hydrolysis for other cell functions) is set to 1.44 µmol·min–1·g wet wt–1 (12). The rate of glucose production is set to 3.80 µmol·min–1·g wet wt–1 on the basis of measurements made in isolated rat liver during glucagon infusion (16). The uniquely defined network fluxes, constrained by Eqs. 35, are plotted in Fig. 1 and tabulated in Table 3. In this case the glycogen phosphorylase flux is 3.834 µmol·min–1·g wet wt–1, meaning that ~99.1% of the glucose 1-phosphate generated from glycogen is converted to glucose. The boundary flux value of 3.80 µmol·min–1·g wet wt–1 for glucose represents the plateau value measured following transient (5-min) infusion with glucagon (16) and not necessarily a prolonged steady state. Thus this value provides an upper-limit estimate on the glycogenolytic flux.

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)
where {Delta}' and {Delta}'0 are vectors listing the steady-state and equilibrium standard reaction free energies for all of the reactions in the system; the vector lists the transformed standard free energies of formation for the reactants (5); R is the gas constant; and T is the absolute temperature. The vector ln contains the logarithms of concentration of the m reactants, and T is the transpose of the stoichiometric matrix of internal chemical reactions. Equation 6 can be rearranged as

(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 {Delta} 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)
where P is a diagonal matrix, with diagonal entries corresponding to the partition coefficients, or fractional intracellular spaces, associated with each metabolite in the system. The left null space of the matrix P–1S can be stored in a matrix L, such that

(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)
where mode 1 and mode 2 are the concentration vectors for the different steady metabolic states. Note that there are typically far more reactants than conserved metabolic pools. Thus Eq. 10 does not uniquely determine concentrations in one state, or metabolic operating mode, based on known values in another. Equation 10 does serve as an additional constraint that two or more feasible steady-state concentration vectors must satisfy. These constraints are applied below to analyzing the regulation of different operating modes in a hepatocyte.

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)
where k+1 and k–1 are the forward and backward reaction rate constants, and Ct is the total concentration, Ct = [A] + [B]. Similarly the steady-state flux through the reversible Michaelis-Menten mechanism,

(14)
can be expressed:

(15)
where E0 is the total enzyme concentration. The Michaelis-Menten constants for binding of A and B are Km,A = (k+2 + k–1)/k+1 and Km,B = (k+2 + k–1)/k–2, respectively. In the limits Ct << Km,A and Ct << Km,B, Eq. 15 can be expressed in a form similar to that of Eq. 13:

(16)

For reactions involving multiple species with the general form

(17)
we define Ct as

(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.



View larger version (17K):
[in this window]
[in a new window]
 
Fig. 2. Generalized flux-potential relationship. The flux-potential relationship of Eq. 19 is plotted as vs. for the case where , for various values of activity, , as indicated in the figure. See methods and results for explanation of terms.

 
Thermodynamic Parameters

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 {Delta}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.



View larger version (23K):
[in this window]
[in a new window]
 
Fig. 3. Predicted free energies and activities for glycogenesis (mode 1) and glycogenolysis (mode 2) obtained via minimization of Eq. 21. Top: predicted reaction potentials are plotted for the reactions in the system. Bottom: predicted differences in activities between the 2 modes are plotted as the ratio on a log scale. When , the activity jth enzyme is higher during glycogensis than during glycogenolysis. When , the activity jth enzyme is higher during glycogenolysis than during glycogensis.

 
The unrealistic result that nearly all enzymes in the system are actively controlled in switching between operational modes allows us to reject minimization of f1 as a useful objective. A more reasonable result is obtained when the differences between the predicted activities are minimized. This hypothesis is formalized by minimizing the objective:

(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.



View larger version (24K):
[in this window]
[in a new window]
 
Fig. 4. Predicted free energies and activities for glycogenesis (mode 1) and glycogenolysis (mode 2) obtained via minimization of Eq. 22. Top: predicted reaction potentials are plotted for the reactions in the system. Bottom: predicted differences in activities between the 2 modes are plotted as the ratio on a log scale. When , the activity jth enzyme is higher during glycogensis than during glycogenolysis. When , the activity jth enzyme is higher during glycogenolysis than during glycogensis. Five enzymes are identified as major regulatory sites for control of glycogensis and glycogenolysis.

 
Reconstruction of Cellular Regulatory Mechanisms

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.



View larger version (19K):
[in this window]
[in a new window]
 
Fig. 5. Predicted concentrations for the glycogenic and glycogenolytic modes are compared. {bullet}, Concentrations computed with the objective of Eq. 21, which minimized the differences in the concentrations; {circ}, concentrations of the minimal control scheme computed via the objective of Eq. 22. Two main outliers, PPi and UDPG, are indicated in the figure.

 
Disease Phenotyping and Target Identification

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.



View larger version (25K):
[in this window]
[in a new window]
 
Fig. 6. Computational profiling of metabolic diseases. Top: sites of action of 2 metabolic diseases impacting glycogenolysis are indicated. Von Gierke's disease arises from malfunctioning glucose-6-phosphatase enzyme (g6pt); Hers' disease is characterized by malfunctioning glycogen phosphorylase (glp). Bottom: predicted glucose production is plotted as a function of fractional activity of the impaired enzyme in both von Gierke's and Hers' diseases. Predicted glucose production is relatively insensitive to glp activity, but highly sensitive to g6pt activity.

 
These results, which match the clinical finding that patients with Hers' disease display milder symptoms than those with von Gierke's (30), are not intuitive from examination of the network structure illustrated in Fig. 1. Thus the FBA/EBA constraint-based model of hepatic metabolism, with the normally functioning enzymes parameterized on the basis of the data from isolated rat liver, predicts clinically significant behavior of two important metabolic diseases.


    DISCUSSION
 TOP
 ABSTRACT
 BIOCHEMICAL THERMODYNAMICS
 METHODS AND RESULTS
 DISCUSSION
 GRANTS
 REFERENCES
 
Thermodynamics and Constraint-Based Modeling

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 {Delta}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)
where Ci and Km,i are the concentration and Michaelis-Menten constant for the ith metabolite. If a given flux can be expressed in terms of metabolite concentrations in an arbitrary functional form, and the associated parameters are known with confidence, then the corresponding steady-state flux expression can be used for this flux instead of the generalized potential-flux relationship. In this case, the flux expression that is used should reproduce the correct thermodynamics. For example, flux should go to zero when the metabolite concentrations are in thermodynamic equilibrium (25–27).

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
 TOP
 ABSTRACT
 BIOCHEMICAL THERMODYNAMICS
 METHODS AND RESULTS
 DISCUSSION
 GRANTS
 REFERENCES
 
This work was supported by National Institutes of Health Grant GM-068610.


    FOOTNOTES
 

Address for reprint requests and other correspondence: D. A. Beard, Biotechnology and Bioengineering Center, Dept. of Physiology, Medical College of Wisconsin, 8701 Watertown Plank Rd., Milwaukee, WI 53226 (E-mail: dbeard{at}mcw.edu)

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
 TOP
 ABSTRACT
 BIOCHEMICAL THERMODYNAMICS
 METHODS AND RESULTS
 DISCUSSION
 GRANTS
 REFERENCES
 

  1. AlbertyRA. Equilibrium compositions of solutions of biochemical species and heats of biochemical reactions. Proc Natl Acad Sci USA 88: 3268–3271, 1991.
  2. AlbertyRA. Degrees of freedom in biochemical reaction systems at specific pH and pMg. J Phys Chem 96: 9614–9621, 1992.
  3. AlbertyRA. Levels of thermodynamic treatment of biochemical reaction systems. Biophys J 65: 1243–1254, 1993.
  4. AlbertyRA. The role of water in the thermodynamics of dilute aqueous solutions. Biophys Chem 100: 183–192, 2003.
  5. AlbertyRA. Thermodynamics of Biochemical Reactions. Hoboken, NJ: Wiley & Sons, 2003.
  6. AlbertyRA and Goldberg RN. Standard thermodynamic formation properties for the adenosine 5'-triphosphate series. Biochemistry 31: 10610–10615, 1992.
  7. BassingthwaighteJB and Goresky CA. Modeling in the analysis of solute and water exchange in the microvasculature. In: Handbook of Physiology. The Cardiovascular System. Microcirculation. Bethesda, MD: Am. Physiol. Soc., 1984, sect. 2, vol. IV, pt. 1, chapt. 13, p. 549–626.
  8. BeardDA, Babson E, Curtis E, and Qian H. Thermodynamic constraints for biochemical networks. J Theor Biol 228: 327–333, 2004.
  9. BeardDA, Liang SD, and Qian H. Energy balance for analysis of complex metabolic networks. Biophys J 83: 79–86, 2002.
  10. BeardDA, Qian H, and Bassingthwaighte JB. Stoichiometric foundation of large-scale biochemical system analysis. In: Modelling in Molecular Biology, edited by Ciabanu G and Rozenberg G. Berlin: Springer, 2004, p. 1–19.
  11. BeardDA, Schenkman KA, and Feigl EO. Myocardial oxygenation in isolated hearts predicted by an anatomically realistic microvascular transport model. Am J Physiol Heart Circ Physiol 285: H1826–H1836, 2003.
  12. BeauvieuxMC, Tissier P, Couzigou P, Gin H, Canioni P, and Gallis JL. Ethanol perfusion increases the yield of oxidative phosphorylation in isolated liver of fed rats. Biochim Biophys Acta 1570: 135–140, 2002.
  13. BonariusHP, Ozemre A, Timmerarends B, Skrabal P, Tramper J, Schmid G, and Heinzle E. Metabolic-flux analysis of continuously cultured hybridoma cells using (13)CO(2) mass spectrometry in combination with (13)C-lactate nuclear magnetic resonance spectroscopy and metabolite balancing. Biotechnol Bioeng 74: 528–538, 2001.
  14. CarlsonR, Fell D, and Srienc F. Metabolic pathway analysis of a recombinant yeast for rational strain development. Biotechnol Bioeng 79: 121–134, 2002.
  15. CovertMW, Schilling CH, Famili I, Edwards JS, Goryanin II, Selkov E, and Palsson BO. Metabolic modeling of microbial strains in silico. Trends Biochem Sci 26: 179–186, 2001.
  16. DoiY, Iwai M, Matsuura B, and Onji M. Glucagon attenuates the action of insulin on glucose output in the liver of the Goto-Kakizaki rat perfused in situ. Pflügers Arch 442: 537–541, 2001.
  17. EdwardsJS, Ibarra RU, and Palsson BO. In silico predictions of Escherichia coli metabolic capabilities are consistent with experimental data. Nat Biotechnol 19: 125–130, 2001.
  18. EdwardsJS and Palsson BO. The Escherichia coli MG1655 in silico metabolic genotype: its definition, characteristics, and capabilities. Proc Natl Acad Sci USA 97: 5528–5533, 2000.
  19. EdwardsJS and Palsson BO. Robustness analysis of the Escherichia coli metabolic network. Biotechnol Prog 16: 927–939, 2000.
  20. GomisRR, Ferrer JC, and Guinovart JJ. Shared control of hepatic glycogen synthesis by glycogen synthase and glucokinase. Biochem J 351: 811–816, 2000.
  21. GuynnRW, Veloso D, Lawson JW, and Veech RL. The concentration and control of cytoplasmic free inorganic pyrophosphate in rat liver in vivo. Biochem J 140: 369–375, 1974.
  22. JiangG and Zhang BB. Glucagon and regulation of glucose metabolism. Am J Physiol Endocrinol Metab 284: E671–E678, 2003.
  23. KorzeniewskiB. Theoretical studies on the regulation of oxidative phosphorylation in intact tissues. Biochim Biophys Acta 1504: 31–45, 2001.
  24. MooreMC, Connolly CC, and Cherrington AD. Autoregulation of hepatic glucose production. Eur J Endocrinol 138: 240–248, 1998.
  25. MulquineyPJ, Bubb WA, and Kuchel PW. Model of 2,3-bisphosphoglycerate metabolism in the human erythrocyte based on detailed enzyme kinetic equations: in vivo kinetic characterization of 2,3-bisphosphoglycerate synthase/phosphatase using 13C and 31P NMR. Biochem J 342: 567–580, 1999.
  26. MulquineyPJ and Kuchel PW. Model of 2,3-bisphosphoglycerate metabolism in the human erythrocyte based on detailed enzyme kinetic equations: computer simulation and metabolic control analysis. Biochem J 342: 597–604, 1999.
  27. MulquineyPJ and Kuchel PW. Model of 2,3-bisphosphoglycerate metabolism in the human erythrocyte based on detailed enzyme kinetic equations: equations and parameter refinement. Biochem J 342: 581–596, 1999.
  28. QianH, Beard DA, and Liang SD. Stoichiometric network theory for nonequilibrium biochemical systems. Eur J Biochem 270: 415–421, 2003.
  29. SchillingCH, Schuster S, Palsson BO, and Heinrich R. Metabolic pathway analysis: basic concepts and scientific applications in the post-genomic era. Biotechnol Prog 15: 296–303, 1999.
  30. StryerL. Biochemistry. New York: WH Freeman, 1995.
  31. StumpelF, Scholtka B, and Jungermann K. Impaired glucose sensing by intrahepatic, muscarinic nerves for an insulin-stimulated hepatic glucose uptake in streptozotocin-diabetic rats. FEBS Lett 436: 185–188, 1998.
  32. VarmaA and Palsson BO. Stoichiometric flux balance models quantitatively predict growth and metabolic by-product secretion in wild-type Escherichia coli W3110. Appl Environ Microbiol 60: 3724–3731, 1994.
  33. VendelinM, Kongas O, and Saks V. Regulation of mitochondrial respiration in heart cells analyzed by reaction-diffusion model of energy transfer. Am J Physiol Cell Physiol 278: C747–C764, 2000.
  34. XueC, Aspelund G, Sritharan KC, Wang JP, Slezak LA, and Andersen DK. Isolated hepatic cholinergic denervation impairs glucose and glycogen metabolism. J Surg Res 90: 19–25, 2000.




This Article
Abstract
Full Text (PDF)
All Versions of this Article:
288/3/E633    most recent
00239.2004v1
Alert me when this article is cited
Alert me if a correction is posted
Citation Map
Services
Email this article to a friend
Similar articles in this journal
Similar articles in ISI Web of Science
Similar articles in PubMed
Alert me to new issues of the journal
Download to citation manager
Search for citing articles in:
ISI Web of Science (4)
Google Scholar
Articles by Beard, D. A.
Articles by Qian, H.
Articles citing this Article
PubMed
PubMed Citation
Articles by Beard, D. A.
Articles by Qian, H.


HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
Visit Other APS Journals Online
Copyright © 2005 by the American Physiological Society.