## Abstract

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, 17–19). 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

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, ATP^{4−}, HATP^{3−}, and H_{2}ATP^{2−}, 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 Δ*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) 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] = [ATP^{4−}] + [HATP^{3−}] + [H_{2}ATP^{2−}]. 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

### Flux Balance and Energy Balance Constraints on Fluxes

The central flux balance conservation statement is given by the equation (3) where *J̄* ∈ is the vector of *n* fluxes occurring in the network, **S** ∈ 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 **S̃**, and **S̃** is the stoichiometric matrix of internal reactions (excluding boundary fluxes). The vector Δ*Ḡ′* lists the free energies associated with each of the reaction fluxes. Vectors in the null space of **S̃** 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. 3*–*5*, 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.

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 CO_{2} is treated as the total concentration of CO_{2} in the aqueous phase, including H_{2}CO_{3}, HCO_{3}^{−}, and CO formed via the hydration of CO_{2} by carbonic anhydrase. Correct treatment of the transformed thermodynamic free energies and chemical equilibria for these reactions requires that H_{2}O 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. 3*–5, 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.

#### 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. 3*–*5*, 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 Δ*Ḡ′* and Δ*Ḡ*′^{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 *C̄* contains the logarithms of concentration of the *m* reactants, and **S̃**^{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, **S̃**^{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) 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**^{−1}**S** can be stored in a matrix **L**, such that (9)

(Note that the matrix **N** contains a basis for the *right* null space of **S̃**.) 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 *C̄*^{mode 1} and *C̄*^{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 *C _{t}* is the total concentration,

*C*= [A] + [B]. Similarly the steady-state flux through the reversible Michaelis-Menten mechanism, (14) can be expressed: (15) where

_{t}*E*is the total enzyme concentration. The Michaelis-Menten constants for binding of A and B are

_{0}*K*

_{m,A}= (

*k*

_{+2}+

*k*

_{−1})/

*k*

_{+1}and

*K*

_{m,B}= (

*k*

_{+2}+

*k*

_{−1})/

*k*

_{−2}, respectively. In the limits

*C*≪

_{t}*K*

_{m,A}and

*C*≪

_{t}*K*

_{m,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 *C _{t}* as (18)

In *Eqs. 17* and *18*, *x _{i}* represents the

*i*

^{th}metabolite, and μ

*and*

_{i}*v*are the stoichiometric coefficients associated with

_{i}*x*on the left- and righthand sides the reaction.

_{i}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.

### 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 CO_{2} (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 *C̄*^{mode 1} and *C̄*^{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)

*X _{j}* is the activity of the

*j*

^{th}enzyme,

*J*is the flux of the

_{j}*j*

^{th}reaction, and

*C*(

_{t,j}*C̄*) and Δ

*G′*(

_{j}*C̄*) 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 *f*_{1} 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 *f*_{1} 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.

The unrealistic result that nearly all enzymes in the system are actively controlled in switching between operational modes allows us to reject minimization of *f*_{1} 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 *f*_{2} so that it measures relative differences between reaction activities in the two metabolic modes.

The results obtained using *f*_{2} 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.

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

### 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, CO_{2}, 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.

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

### 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 H_{2}O = Glc + ADP + 2 P_{i}. 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. 3–5. 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) where *C _{i}* and

*K*are the concentration and Michaelis-Menten constant for the

_{m,i}*i*

^{th}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

This work was supported by National Institutes of Health Grant GM-068610.

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

- Copyright © 2005 by American Physiological Society