Application of biochemical systems theory to metabolism in human red blood cells. Signal propagation and accuracy of representation.

Human erythrocytes are among the simplest of cells. Many of their enzymes have been characterized kinetically using steady-state methods in vitro, and several investigators have assembled this kinetic information into mathematical models of the integrated system. However, despite its relative simplicity, the integrated behavior of erythrocyte metabolism is still complex and not well understood. Errors will inevitably be encountered in any such model because of this complexity; thus, the construction of an integrative model must be considered an iterative process of assessment and refinement. In a previous study, we selected a recent model of erythrocyte metabolism as our starting point and took it through three stages of model assessment and refinement using systematic strategies provided by biochemical systems theory. At each stage deficiencies were diagnosed, putative remedies were identified, and modifications consistent with existing experimental evidence were incorporated into the working model. In this paper we address two issues: the propagation of biochemical signals within the metabolic network, and the accuracy of kinetic representation. The analysis of signal propagation reveals the importance of glutathione peroxidase, transaldolase, and the concentration of total glutathione in determining systemic behavior. It also reveals a highly amplified diversion of flux between the pathways of pentose phosphate and nucleotide metabolism. In determining the range of accurate representation based on alternative kinetic formalisms we discovered large discrepancies. These were identified with the behavior of the model represented within the Michaelis-Menten formalism. This model fails to exhibit a nominal steady state when the activity of glutathione peroxidase is decreased by as little as 9%. Our current understanding, as embodied in this working model, is in need of further refinement, and the results presented in this paper suggest areas of the model where such effort might profitably be concentrated.

Human erythrocytes are among the simplest of cells. Many of their enzymes have been characterized kinetically using steady-state methods in vitro, and several investigators have assembled this kinetic information into mathematical models of the integrated system. However, despite its relative simplicity, the integrated behavior of erythrocyte metabolism is still complex and not well understood. Errors will inevitably be encountered in any such model because of this complexity; thus, the construction of an integrative model must be considered an iterative process of assessment and refinement. In a previous study, we selected a recent model of erythrocyte metabolism as our starting point and took it through three stages of model assessment and refinement using systematic strategies provided by biochemical systems theory. At each stage deficiencies were diagnosed, putative remedies were identified, and modifications consistent with existing experimental evidence were incorporated into the working model. In this paper we address two issues: the propagation of biochemical signals within the metabolic network, and the accuracy of kinetic representation. The analysis of signal propagation reveals the importance of glutathione peroxidase, transaldolase, and the concentration of total glutathione in determining systemic behavior. It also reveals a highly amplified diversion of flux between the pathways of pentose phosphate and nucleotide metabolism. In determining the range of accurate representation based on alternative kinetic formalisms we discovered large discrepancies. These were identified with the behavior of the model represented within the Michaelis-Menten formalism. This model fails to exhibit a nominal steady state when the activity of glutathione peroxidase is decreased by as little as 9%. Our current understanding, as embodied in this working model, is in need of further refinement, and the results presented in this paper suggest areas of the model where such effort might profitably be concentrated.
The primary physiological function of the red blood cell is to bind oxygen in the lung and transport it to oxygen-consuming cells elsewhere in the body. In addition, there are several metabolic functions that a red blood cell must perform for its own survival. These include generation of metabolic energy (e.g. ATP), generation of reducing agents (e.g. NADH and NADPH), generation of 2,3-diphosphoglycerate, and maintenance of ionic and concentration gradients across the cell membrane (for detailed reviews, see e.g. Refs. 1 and 2). These functions involve a minimum set of metabolic pathways, including the glycolytic sequence, the pentose phosphate pathway, the 2,3-diphosphoglycerate shunt, and the pathways of nucleotide metabolism; the ability to synthesize RNA, protein, lipid, and purines has been lost. Over a period of 3 decades most of the molecular components that constitute the pathways of erythrocyte metabolism have been experimentally characterized. The kinetic data that have been accumulated are derived for the most part from measurements made in vitro, rather than in vivo, and they have been obtained under various experimental conditions (e.g. pH and temperature) from many independent laboratories. These data have formed the basis for a number of kinetic models dealing with various aspects of the human erythrocyte (e.g. [3][4][5][6][7][8][9][10][11][12][13][14][15][16]. The elucidation of a complex biochemical system like the metabolism of human erythrocytes is no trivial task (17). It is not sufficient to have detailed kinetic information for the individual reactions considered in isolation, even if this information were to reflect accurately the actual conditions in vivo. One needs in addition a systematic framework for integrating this type of information into a model of the intact system and rigorous methods for extracting the systemic implications latent within such a model. Moreover, if this approach is to succeed, one must demonstrate that the model is of sufficient quality to warrant our confidence.
In developing an integrated model of metabolism in human erythrocytes we have shown that the need for a systematic framework and rigorous methodology can be met by the powerlaw formalism 1 and biochemical systems theory (18,19). In the process, a comprehensive model of metabolism in human erythrocytes (12)(13)(14)(15) has been taken through three stages of assessment and refinement in an effort to develop a model that would justify a detailed systemic analysis. We found that the original model did not exhibit an appropriate steady state; this deficiency was identified with several inconsistencies, and these were resolved to produce a new reference model (18). Next we found that the steady state was unstable; a systematic search for interactions that could stabilize the steady state revealed one for which there is independent experimental support, and this was included in the model (18,19). Finally, we found that the model lacked robustness; the pattern of parameter sensi-tivities identified the pathways of nucleotide metabolism with this problem, and the addition of previously overlooked mechanisms to the model improved robustness (18,19). The result to emerge from this process (Model III) now will be subject to systemic analysis.
In this paper we shall first analyze the pattern of biochemical signal propagation that determines the steady-state behavior of metabolism in erythrocytes. This analysis reveals the importance of glutathione peroxidase, transaldolase, and the concentration of total glutathione in determining the systemic concentrations and fluxes, as well as the critical behavior of nucleotide metabolism. Second, we analyze the range of accurate representation by comparing the behavior of alternative representations. This analysis shows that the nominal steady state for this model exists only over a very narrow range of variation in some of the enzyme activities and suggests that the model is still in need of refinement. Taken together, the results suggest areas of the model where additional experimental effort might be focused.

MODEL OF METABOLISM IN HUMAN RED BLOOD CELLS
A schematic diagram of Model III representing metabolism in human red blood cells is given in Fig. 1. The fundamental equations that characterize any such metabolic network are Kirchhoff's node equations, which also are referred to as mass balance equations. These can be written in the following form 2 where i ϭ 1, 2, . . . n; an individual metabolite concentration is represented by X i , an individual flux from X i to X j is represented by ij , aggregate influx and aggregate efflux for the X i pool are given by V ϩi and V Ϫi , and the numbers of dependent and independent concentration variables are n and m (20). Kirchhoff's node equations for the specific system in Fig. 1 are given in Table I, using both the above notation, which correlates the flux with its substrate and product, and common mnemonic abbreviations. The 33 dependent concentration variables of this system are enumerated in the third column of Table I. The individual fluxes entering and leaving each node are given in the fourth column. The numbering in this case involves 58 independent concentration variables as well as the 33 dependent concentration variables. The grouping of individual fluxes indicates an aggregation strategy in which the net flux through a reaction is considered as a single flux. This is one of the most common forms of aggregation in the biochemical literature (21). Other aggregation strategies will be discussed in subsequent sections. The names of the dependent metabolites abbreviated in the first column of Table I are given in Table II. The subscripts of the s in the second column are abbreviations for enzymes or reactions (see Table III for full names and EC numbers). The rate laws for these reactions are represented in most cases by conventional expressions within the Michaelis-Menten formalism. The detailed rate laws are given elsewhere (19). not the case, one can preface the reaction velocities with the appropriate stoichiometry as in Table I.  Table II. The full names of enzymes or reactions are given in Table III. The concentration of reduced glutathione (GSH ϩ ) X 92 * is a constrained variable whose value must equal the difference between the total concentration of glutathione (G T ) X 80 , which is an aggregate variable that is also independent, and the concentration of oxidized glutathione (GSSG) X 4 , which is a dependent variable in this model. Other independent variables corresponding to fixed, extracellular, or aggregate metabolites, some of which are shown explicitly in this diagram, include the concentrations of glucose (GLC) X 77 , adenine (ADE) X 78 , CO 2 X 79 , total magnesium (MG T ) X 81 , total nicotinamide adenine dinucleotide (N T ) X 82 , total nicotinamide adenine dinucleotide phosphate (NP T ) X 83 , inorganic phosphate (P i ) X 84 , extracellular lactate (LAC e ) X 85 , extracellular pyruvate (PYR e ) X 86 , extracellular potassium (K e ) X 87 , extracellular sodium (Na e ) X 88 , and extracellular hypoxanthine (HX e ) X 89 . Arrows pointing in one direction represent essentially irreversible reactions; arrows pointing in both directions represent reversible reactions. The rate laws, which account for regulatory interactions not shown in this diagram, are given in Ref. 19. In addition to the 33 dependent variables of this model, there are 43 independent variables corresponding to the activities of the reactions numbered 34 -76, another 15 independent variables corresponding to the concentration variables numbered 77-91 (see Table IV), and one constrained variable corresponding to the concentration of reduced glutathione (numbered 92), whose value is determined by the difference between the total concentration of glutathione, which is an independent variable, and oxidized glutathione, which is a dependent variable.
The steady-state values for the concentrations of this model are calculated by conventional methods (22) and are shown in Table II along with the nominal steady state that is observed experimentally. The steady-state values calculated for aggregate fluxes through the reactions of this model are shown in Table III. METHODS We follow the customary procedure in systems analysis of distinguishing between dependent and independent variables and between steady state and transient behavior (23). The values for the dependent concentrations and fluxes of a biochemical system are determined by the internal dynamics of the system as dictated by Kirchhoff's node laws. The values for the independent concentrations and fluxes are determined independently by agents outside of the system; these variables can be fixed by experimental means, and they characterize the environment of the system. For our purposes the steady-state behavior is that which persists in a constant environment. The transient behavior is that which is exhibited by the system as it progresses from one steady state to another.
Characterization of the Systemic Steady State-When the aggregate rate laws in Equation 1 are represented by products of power-law functions one obtains the S-system representation within the powerlaw formalism (19, 20) Ϫhij and the subscript 0 indicates that the results are evaluated at the nominal steady state. The rate constants of the aggregate rate laws are ␣ i and ␤ i ; the kinetic orders of the aggregate rate laws are given by g ij and h ij . These variables and parameters are all nonnegative quantities, except for the kinetic orders, which have negative values when they represent inhibitory interactions. The explicit S-system representation for Model III is given in the Appendix. The steady-state solution of Equation 2 can be obtained explicitly (24) and expressed in matrix form as where (in matrix form)  10,11 Ϫ v 11,10 ) ϩ (v 8,9 Ϫ 14,15 Ϫ v 15,14 ) Ϫ (v 15,16 Ϫ v 16,15 ) Ϫ v 15 15,16 Ϫ v 16,15 ) ϩ v 18,19 ϩ (v 27,26 Ϫ v 26,27 ) Ϫ v 28, 10 Ϫ v 11,12 Ϫ v 25, 26 Ϫ v 28, 27 Ϫ v 24,23 Ϫ v 6 Within the S-system representation, steady-state fluxes can be expressed in the following matrix notation.
Henceforth, we shall refer to the flux through the X i pool in steady state simply as V i , since V ϩi ϭ V Ϫi ϭ V i under these conditions. Other fluxes can be defined by appropriate rate laws within the power-law formalism, and their steady-state values can be determined in exactly the same fashion from the steady-state values of the concentration variables. The explicit solution in Equations 3 and 4 completely determines the steady-state values for the dependent concentrations and fluxes through pools in terms of the values for the independent metabolite concentrations and the parameters of the S-system representation.
Logarithmic Gains-A change in the value of an independent variable (e.g. the concentration of an aggregate variable, the concentration of an environmental constituent, or the level of an enzyme) represents a signal that will be propagated throughout the system and cause a corresponding change in the dependent variables (metabolite concentrations and fluxes). If the resulting change in a dependent variable is greater than the original signal, then one speaks of gain (or amplifica-tion) in signal during propagation; if the resulting change is less, then one speaks of attenuation. In either case, the quantitative measure of signal propagation from one point in a biochemical network (the input) to another (the output) is a gain factor representing the ratio of output to input signal. These gain factors are important systemic properties that define the influence exerted by independent variables on the dependent variables of the system. In biochemical systems theory, the natural variables are the logarithms of metabolite concentrations and fluxes, and the gain factors are referred to as logarithmic gain factors.
The logarithmic gain in a dependent variable X i with respect to change in an independent variable X j is defined as follows where i ϭ 1, 2, . . . n; j ϭ nϩ1, nϩ2, . . . nϩm; and again the subscript 0 indicates values determined at the steady state. To obtain the theoretically expected value, one simply differentiates the explicit solution for the dependent variable (Equation 3) with respect to the independent variable. The full set of these logarithmic gain factors can be written in the form of an n by m matrix.
͓L͔ ϭ Ϫ͓M͔͓A͔ id (Eq. 6) The elements of the M and A matrices contain all of the kinetic orders and only the kinetic orders. Thus, these logarithmic gains are clearly systemic properties that in general depend upon all of the molecular interactions within the system (25). Similarly, the logarithmic gain in a dependent variable V i with respect to change in an independent variable X j is defined as follows  where i ϭ 1, 2, . . . n and j ϭ nϩ1, nϩ2, . . . nϩm.
To obtain the theoretically expected value, one simply differentiates the explicit solution for the dependent variable (Equation 4) with respect to the independent variable. The full set of these logarithmic gain factors also can be written in the form of an n by m matrix where, again, the subscripts d and id signify that the corresponding matrices contain kinetic orders for dependent and independent variables, respectively. As noted above, other fluxes can be defined by appropriate rate laws within the power-law formalism, and their logarithmic gains can be determined in exactly the same fashion from the logarithmic gains in the concentration variables. For example, given that the forward and reverse processes for the enzyme transketolase 1 have rate laws (19) The corresponding logarithmic gains in flux that result from a change in the concentration of the enzyme glutathione peroxidase (X 34 ) are given by where the logarithmic gains in concentration are obtained from the primary solution (Equation 3). The logarithmic gain for the net flux in the forward direction can then be obtained by the appropriate average of the logarithmic gains for the fluxes in the forward and reverse directions.
In general there is no explicit steady-state solution for the Michaelis-Menten model and, hence, no explicit expression for the logarithmic gains. Nevertheless, for purposes of comparison one can approximate the logarithmic gains of this model by an empirical procedure. First, make a small finite change in an independent variable. Second, determine the resulting change in the steady-state value of the dependent variable. We use the Newton-Raphson method for nonlinear systems (20,22) to calculate a numerical approximation to the steady-state concentrations and fluxes. Finally, take the appropriate ratios. The results can be expressed as follows where the subscript 0 indicates values determined at the steady state. The difference between X j and X j0 is a small finite value ⑀. As ⑀ approaches zero, the empirical values for these logarithmic gains converge to the theoretically expected values obtained for the S-system model (20). A logarithmic gain with a magnitude greater than 1 implies amplification of the original signal; a magnitude less than 1 indicates attenuation. A positive sign for the logarithmic gain indicates that the changes are in the same direction, both increase in value or both decrease. A negative sign indicates that the changes are in the opposite direction.
The influence of a given independent variable on a particular dependent variable is given by the magnitude of the corresponding logarithmic gain, e.g. ͉L(V i ,X j )͉ or ͉L(X i ,X j )͉. The total influence over a particular dependent variable is given by the sum of the individual influences by each independent variable. Accuracy Analysis-We have shown elsewhere (19) that within computational precision the S-system model and the Michaelis-Menten model exhibit the same nominal steady state. This is no surprise since the parameters of the S-system model were determined at the nominal steady state of the Michaelis-Menten model. A more critical test is whether a model can accurately predict changes in the steady-state values of the dependent variables. Unfortunately, experimental data of this type that might be used to evaluate the models are not yet available. In the interim, we shall compare the alternative models with each other. Where they differ we can attempt to determine which of the two models is the more biologically realistic.
We shall use two alternative measurements of fidelity for a kinetic model. First is the error of representation, which is defined as the difference between the output signals generated by the model and the reference in response to a small change in their input signals. In operational terms, one makes a 2% change in an independent variable, determines the response of a dependent variable, and measures the percentage difference between the predicted response of the model and the actual response of a reference system. Second is the range of accurate representation, which is defined as the range of variation in input signals within which the output signals of the model and the reference system being represented differ by less than a specified tolerance. In this case one makes changes of increasing size in an independent variable, determines the corresponding responses of a dependent variable, measures the percentage difference between the predicted response of the model and the actual response of a reference system, and determines the range of variation in the independent variable for which the percentage difference is within 10%.

RESULTS
A schematic representation of metabolism in human erythrocytes is given in Fig. 1. Although the regulatory interactions are not shown in this schematic, they are incorporated into the detailed rate laws, which are summarized in the Appendix. The names of the 33 dependent metabolite concentrations are given in Table II; those of the 58 independent variables are given in Tables III and IV. These variables are grouped into four categories in an effort to facilitate the recognition of significant patterns in the profiles of logarithmic gains and accuracy.
• Group I contains variables related to the pentose pathway.
Dependent variables are X 1 to X 9 and V 1 to V 9 ; independent are X 34 to X 43 .
• Group II contains variables related to glycolysis. Dependent variables are X 10 to X 24 and V 10 to V 24 ; independent are X 44 to X 61 . • Group III contains variables related to nucleotide metabolism. Dependent variables are X 25 to X 33 and V 25 to V 33 ; independent are X 62 to X 76 . • Group IV contains variables related to aggregate and extracellular metabolites. Independent variables are X 77 to X 91 .
The logarithmic gains characterize the propagation of biochemical signals throughout the system, and they reflect the influence exerted by independent variables on the dependent variables of the intact system. These gains are obtained from the explicit solution for the steady state that is available within the framework of the S-system representation of biochemical systems theory (see under "Methods"). The results are summarized in Figs. 2-4.
Logarithmic Gains in Concentration-The influence of the 58 independent variables (X 34 through X 91 ) on the 33 metabolite concentrations has the distribution shown in the three-dimensional bar chart of Fig. 2. The height of each bar represents the magnitude of the corresponding logarithmic gain, ͉L(X i ,X j )͉. Magnitudes less than 10 Ϫ3 (which means the input signal will be attenuated by a factor of 1/1000) are shown as black squares. There are five important patterns revealed in this figure.
First, there are five metabolites whose concentrations are unaffected by change in most independent variables (X 1 , X 4 , X 20 , X 22 , and X 19 ). This can be seen most clearly in the twodimensional bar chart projected on the right in Fig. 2, where the magnitudes for a given metabolite concentration X i with respect to change in the independent variables X j are summed over all j. The solid portion of each bar represents the sum of the positive logarithmic gains, whereas the hatched portion represents the sum of the negative logarithmic gains. The concentration of gluconolactone 6-phosphate X 1 has a pattern of influence that is indicative of a rate-limiting step in the pentose phosphate pathway. Gluconolactone 6-phosphate is the substrate of an essentially irreversible reaction, and, with the exception of changes in the activity of the enzyme that catalyzes this reaction, its concentration should vary directly with the flux through the reaction. Aside from the exception, X 1 is influenced only by X 34 (1.00) and X 80 (1.00). This indicates that the flux through the proximal portion of the pentose phosphate pathway is limited by the activity of glutathione peroxidase X 34 and by the concentration of total glutathione X 80 . This is to be expected because the pool of reduced glutathione (GSH) is at its maximum, essentially equal to the concentration of total glutathione (G T ) X 80 , and the rate law for glutathione peroxidase is represented as a simple first-order process (the product of substrate concentration and a rate constant proportional to enzyme activity). The concentration of gluconolactone 6-phosphate X 1 varies inversely with changes in the level of the enzyme 6-phosphogluconolactonase X 37 (Ϫ1.00), which is expected in this case because the flux is fixed. The concentration of oxidized glutathione X 4 is affected only by changes in X 34 FIG. 2. Influence of independent variables on metabolite concentrations as determined by the magnitudes of the logarithmic gains. The three-dimensional plot displays the magnitudes of the logarithmic gains as a function of the metabolite concentrations X i and the independent variables X j . Logarithmic gains with magnitudes less than 1/1000 are shown as black squares. The two-dimensional projection on the right gives the magnitudes for a particular metabolite concentration X i summed over all the independent variables X j . The twodimensional projection on the left gives the magnitudes for an independent variable X j summed over all the metabolite concentrations X i . The solid bars in each projection represent the sum of the positive gains, and the hatched bars represent the sum of the negative gains. The three-dimensional plot displays the magnitudes of the logarithmic gains as a function of the fluxes V i and the independent variables X j . Logarithmic gains with magnitudes less than 1/1000 are shown as black squares. The two-dimensional projection on the right gives the magnitudes for a particular flux V i summed over all the independent variables X j . The twodimensional projection on the left gives the magnitudes for an independent variable X j summed over all the fluxes V i . The solid bars in each projection represent the sum of the positive gains, and the hatched bars represent the sum of the negative gains.
(0.944), X 80 (0.997), and the level of glutathione reductase X 35 (Ϫ0.979). The explanation is the same as that given above; X 4 varies directly with the independent variables that determine the flux and inversely with the level of the enzyme for which it is a substrate. The concentration of lactate X 20 is influenced primarily by the level of the lactate transport process X 58 (Ϫ0.120) and by the concentration of extracellular lactate X 85 (0.879). The concentration of NADH X 22 is influenced by the concentrations of extracellular lactate X 85 (0.575) and extracellular pyruvate X 86 (Ϫ0.662), and by the concentration of total nicotinamide adenine dinucleotide X 82 (0.992). The concentration of pyruvate X 19 is unaffected by changes in almost all independent variables. This is indicated by the row of black squares for the dependent variable X 19 . This invariance can be understood as follows. In steady state the concentration of pyruvate within the cell must equal the concentration of pyruvate outside the cell, which is a constant in this model, otherwise there would be some net flux entering or leaving the intracellular pool X 19 . If there were such a flux, then the flux through glyceraldehyde phosphate dehydrogenase would not equal the flux through lactate dehydrogenase, and there could be no steady state for the pools of NAD and NADH. Thus, the only independent variable that has an influence on the intracellular concentration of pyruvate is the concentration of extracellular pyruvate X 86 (1.00).
Second, the seven metabolite concentrations most influenced by change in the independent variables are those for gluconate 6-phosphate X 2 , nicotinamide adenine dinucleotide phosphate X 3 , sedoheptulose 7-phosphate X 8 , glucose 6-phosphate X 10 , fructose 6-phosphate X 11 , 5-phosphoribosyl-1-pyrophosphate X 29 , and inosine X 31 . This can be seen most clearly in the two-dimensional bar chart projected on the right in Fig. 2. There are large negative logarithmic gains as well as positive ones. The distribution of total influence exerted over each of the metabolite concentrations varies from a low of 1.00 -3.48 for concentrations associated with the first pattern (X 1 , X 4 , X 20 , X 22 , and X 19 ) to a high of 22.7-35.2 for concentrations associated with the second pattern (X 2 , X 3 , X 8 , X 10 , X 11 , X 29 , and X 31 ). The large gradient in the changes between sedoheptulose 7-phosphate X 8 (35.2) and erythrose 4-phosphate X 9 (4.97), and the large increases in 5-phosphoribosyl-1-pyrophosphate X 29 (22.7), suggest that the flux through the distal portion of the pentose phosphate pathway is restricted, and that into nucleotide metabolism is augmented. Further support for this hypothesis is given below.
Third, the independent variables that have the greatest total influence on the metabolite concentrations are the activities of the enzymes glutathione peroxidase X 34 (42.1) and transaldolase X 43 (40.4) and the concentration of total glutathione X 80 (42.1). This can be seen in the two-dimensional projection on the left, where the magnitudes of the logarithmic gains for the dependent concentrations X i with respect to change in a given independent variable X j are summed over all i. Again, contributions from positive and negative logarithmic gains are shown separately by solid and hatched bars. These three variables determine the flux through the pentose phosphate pathway and its diversion into nucleotide metabolism. Those variables with the next greatest level of total influence are the concentrations of total phosphate X 84 (31.4) and total magnesium X 81 (23.9) and the activities of the enzymes hexokinase X 44 (27.9) and adenosine monophosphate phosphohydrolase X 62 (26.4), all of which are concerned with nucleotide metabolism.
Fourth, there are many independent variables that have almost no effect on any of the metabolite concentrations in the model. The distribution of total influence exerted by each of these independent variables ranges from 0 to 2.08. As can be seen in the two-dimensional projection on the left in Fig. 2, these are associated with the pentose phosphate pathway (X 35 to X 42 , X 79 , and X 83 ) glycolysis (X 45 , X 48 , X 49 , X 53 , X 54 , X 56 to X 58 , X 77 , X 82 , and X 87 ), and nucleotide metabolism (X 65 , X 67 to X 70 , X 73 , X 75 , X 76 , X 78 , X 89 , and X 90 ). This list includes the activities for 13 of the 16 fast reactions; it also includes the activities for four reactions that follow a rate-limiting step and three reactions involved in transport of intra-and extracellular metabolites that are near equilibrium.
Fifth, if one sums the signed logarithmic gains with respect to the enzyme concentrations (X 34 through X 76 ), one sees that the value is zero for each metabolite concentration (data not shown). Although this is not generally true (26), it is to be expected for a model, like the one analyzed in this paper, in which each of the rate laws is assumed to be independent and proportional to the concentration of the corresponding enzyme (26 -28).
Logarithmic Gains in Flux-The influence of the independent variables on the fluxes through the metabolite pools has the distribution shown in the three-dimensional bar chart of Fig. 3. Again, logarithmic gains that have magnitudes less than 10 Ϫ3 are shown as black squares, and there are five important patterns to be emphasized in this figure.
First, five of the fluxes are unaffected by change in most of the independent variables. This can be clearly seen from the two-dimensional projection on the right, where the magnitudes of the logarithmic gains for a given flux V i with respect to change in the independent variables X j are summed over all j. The solid portion of each bar represents the sum of the positive logarithmic gains, whereas the hatched portion represents the sum of the negative logarithmic gains. The fluxes through the pools of gluconolactone 6-phosphate V 1 and glutathione V 4 are affected only by two independent variables, the activity of glutathione peroxidase X 34 (0.992) and the intracellular con- Only key branch points and net fluxes are represented. The three aggregate nodes, which also must satisfy Kirchhoff's node law, include the proximal portion of the pentose phosphate pathway, the distal portion of the pentose phosphate pathway and most of glycolysis, and nucleotide metabolism. The net fluxes can be identified with the following enzymes: F 1 , hexokinase; F 2 , glucose-6-phosphate dehydrogenase; F 3 , xylulose-5-phosphate isomerase; F 4 , adenine phosphoribosyl transferase; F 5 , hypoxanthine transport process; F 6 , inosine transport process; F 7 , phosphoglucoisomerase; F 8 , ribose-5-phosphate isomerase; F 9 , transketolase I; F 10 , adenosine transport process; and F 11 , enolase. The numbers associated with the net fluxes are their logarithmic gains in response to a change in the concentration of glutathione peroxidase X 34 , except in the cases of F 6 and F 10 . In these two cases the net flux is zero at steady state, so the logarithmic gains are undefined. Instead, we have shown with an asterisk (*) the logarithmic gain for the efflux, which is indicative of the change in net flux because the influx is unchanged. See Fig. 1 for a more detailed metabolic map of this system. centration of total glutathione X 80 (0.992). This is because the pool of reduced glutathione (GSH) is at its maximum, essentially equal to the concentration of total glutathione (G T ), and thus the flux through the pentose phosphate pathway is limited by the activity of glutathione peroxidase X 34 and the intracellular concentration of total glutathione X 80 . The flux through the pool of sodium V 24 is a function only of the level of sodium leakage X 60 (1.00) and of the concentration of external sodium X 88 (0.419) and is therefore insensitive to change in the value of all other independent variables. The flux through the pool of potassium V 23 is only affected by changes in the activities of the pump X 61 (0.322) and the leakage mechanisms for sodium X 60 (0.170) and potassium X 59 (0.494); all other influences that would be mediated via changes in ATP have no effect because the ATP concentration is far in excess of its K m for the pump (18). The flux through 2,3-diphosphoglycerate V 21 is largely a function of diphosphoglycerate phosphatase levels X 52 (0.959) because this enzyme is operating at nearly its maximum rate with substrate concentrations far in excess of its K m .
Second, the fluxes that are most influenced by change in the independent variables are those through the pools of nucleotide metabolism (adenosine monophosphate V 26 , adenosine diphosphate V 27 , adenosine triphosphate V 28 , inosine V 31 , hypoxanthine V 32 , and ribose 1-phosphate V 33 ) followed by those through glycolytic intermediates (glucose 6-phosphate V 10 and fructose 6-phosphate V 11 ) and the pentose phosphate pathway (ribose 5-phosphate V 6 ). This can be seen from the two-dimensional projection on the right, which shows large negative logarithmic gains as well as positive ones. The distribution of total influence exerted over each of the fluxes varies from a low of 1.28 -2.40 for fluxes associated with the first pattern (V 1 , V 4 , V 21 , V 23 , V 24 ) to a high of 19.5-29.1 for the fluxes associated with the second pattern (V 26 , V 27 , V 28 , V 31 , V 32 , V 33 , V 10 , V 11 , V 6 ).
Third, the independent variables that have the greatest total influence on the dependent fluxes are the activities of the enzymes glutathione peroxidase X 34 (38.1) and transaldolase X 43 (33.4) and the concentration of total glutathione X 80 (38.1). This can be seen in the two-dimensional bar chart projected on the left. Increases in X 34 and X 80 and decreases in X 43 lead to large increases in flux through the lower portions of the pentose phosphate pathway and nucleotide metabolism. This will be examined below from another perspective (see Fig. 4). The independent variables with the next greatest total influence are total phosphate X 84 (27.1) and total magnesium X 81 (20.2) and the activities of the enzymes hexokinase X 44 (21.9) and adenosine monophosphate phosphohydrolase X 62 (25.5). This is the same qualitative pattern that was observed for the metabolite concentrations (see Fig. 2).
Fourth, many independent variables have almost no effect on any of the fluxes in the model (see the two-dimensional projection on the left in Fig. 3). Among these are metabolites and activities of enzymes associated with the pentose phosphate pathway (X 35 to X 42 , X 79 , and X 83 ). This is to be expected since the flux through this pathway is limited by the activity of glutathione peroxidase X 34 , as noted above. Also among these are metabolites and activities associated with glycolysis (X 45 , X 48 , X 49 , X 53 , X 54 , X 56 , X 58 , X 59 , X 77 , and X 87 ) and nucleotide metabolism (X 67 , X 68 , X 73 , X 75 , X 76 , X 78 , X 89 , and X 90 ). Many, though not all, of these reactions are fast, and some fast reactions are not included in this list. Again, the pattern is very similar to that for the metabolite concentrations.
Fifth, if one sums the signed logarithmic gains with respect to the enzyme concentrations (X 34 through X 76 ), one observes that the value is unity for each flux (data not shown). Again, this is to be expected for the special case of a model that is assumed to be homogeneous of degree one in the enzyme concentrations, although it is of no general biochemical significance. There are negative as well as positive values in each of these sums, and the negative values can be quite large, which is frequently the case (26, 28 -33) and contrary to what some investigators have claimed. This provides further evidence against the notion that a sum of unity is equivalent to a distribution function for influence, with some enzymes required to have less if others have more. In general, these summation relationships are neither necessary (26) nor sufficient (32) for a valid steady-state analysis.
Logarithmic Gains in Net Flux-The distribution of flux at the branch points of the pentose phosphate pathway is critical. This can be seen more clearly if one considers an alternative representation of Kirchhoff's node equations (21) that allows one to focus on the key fluxes and branch points. Following the procedure described in Shiraishi and Savageau (20), one can consider "aggregate nodes" for which Kirchhoff's node laws must apply in steady state. In Fig. 4 we show three such aggregate nodes: one enclosing the proximal portion of the pentose phosphate pathway, a second enclosing the distal portion of the pentose phosphate pathway and most of glycolysis, and a third enclosing nucleotide metabolism. The relevant net fluxes are numbered F 1 to F 11 (although the last five of these can be expressed in terms of the first six). The logarithmic gains in these fluxes in response to a change in the concentration of glutathione peroxidase X 34 can be readily calculated from the rate laws for the individual processes, the logarithmic gains in the associated metabolic pools, and ratios involving the forward and reverse fluxes. An example is given under "Methods." The results for the 11 net fluxes are shown in Fig. 4. The response to an increase in level of glutathione peroxidase shows that the increase in flux into the pentose phosphate pathway F 2 (1.00) is preferentially diverted from glycolysis F 3 (0.897) to nucleotide metabolism F 8 (1.21). Moreover, this results in a net increase in flux into nucleotide metabolism since the increase in F 8 (1.21) is greater than the increase in F 9 (0.897), which returns flux to glycolysis. The same pattern can be seen in the response to an increase in the concentration of total glutathione or a decrease in the level of transaldolase. Hence, in each case an increased proportion of the flux through ribose 5-phosphate X 6 is being diverted into nucleotide metabolism. This result shows that the distribution of flux at the branch points of the pentose phosphate pathway is a critical feature of this model. The diversion of flux from glycolysis to nucleotide metabolism is accompanied by a dissociation in the nominal rates of ATP production and utilization by the glycolytic kinases. The difference is rectified by appropriate changes in the rates of adenosine kinase and adenosine triphosphate phosphohydrolase and to a lesser extent by changes in the rates of Na/K-ATPase, adenylate kinase, and phosphoribosyl pyrophosphate synthetase.
For comparison, we also have determined the logarithmic gains empirically within the framework of the Michaelis-Menten representation. This involves making a small change in an independent variable, iteratively solving for the new steady state, and then calculating the logarithmic gains by taking finite differences. The empirically determined results are identical, within computational precision, to the theoretically expected results in Figs. 2, 3, and 4 when independent variables are changed by Ϯ 10 Ϫ8 %. However, this method is computationally inefficient; it takes 172 times longer to compute the logarithmic gains for the Michaelis-Menten model than it does for the S-system model. 3 Accuracy Analysis-We compare the power-law and Michaelis-Menten representations of Model III by solving for the steady-state behavior of their dependent variables in response to changes in their independent variables. The solution for the S-system representation within the power-law formalism is straightforward. Equations 3 and 4 completely characterize the relationship between input signals (changes in independent variables Y] id ) and output signals (changes in dependent variables Y] d and log V]). There is no such explicit solution within the Michaelis-Menten formalism. Instead, one needs to employ more complicated and time-consuming iterative algorithms, such as the Newton-Raphson method for nonlinear systems (22), to obtain a solution.
We determined the error of representation as described under "Methods" by changing each independent variable by Ϫ2% and then calculating the change produced in each dependent variable for both the S-system model and the Michaelis-Menten model considered as a reference. Only enzyme levels are considered in this analysis because the Michaelis-Menten model failed to realize a steady state when independent variables other than enzyme levels were changed by as little as 2% (data not shown).
The three-dimensional plot in Fig. 5 shows the magnitude of the error as a function of the independent and dependent concentration variables. Differences with magnitude less than 2.0 ϫ 10 Ϫ5 are shown as black squares. The changes produced by both models are in most cases identical within computational precision (black squares). The average magnitude of all the percentage differences is 0.0104%. The largest differences are produced by changes in the levels of glutathione peroxidase X 34 and transaldolase X 43 . The metabolic pools that exhibit the greatest differences are gluconate 6-phosphate X 2 , nicotina-mide adenine dinucleotide phosphate X 3 , sedoheptulose 7-phosphate X 8 , glucose 6-phosphate X 10 , fructose 6-phosphate X 11 , and ATP X 28 . A very similar pattern of results was obtained when independent variables were changed by ϩ2% (data not shown), and the average magnitude of all the percentage differences is 0.00995%.
To explore further the nature of the deviation between these alternative representations we have examined the range of accurate representation (see "Methods") in response to changes in the level of glutathione peroxidase X 34 . The four most significant differences are plotted in Fig. 6. The Michaelis-Menten model shows drastic changes when the level of glutathione peroxidase is reduced by as little as 9%. Beyond this point the Michaelis-Menten model does not exhibit a steady state. This behavior may be the result of a positive feedback effect: A decrease in glutathione peroxidase level lowers the influx to the pentose pathway, which in turn leads to a diversion of flux away from nucleotide metabolism. The resulting depletion of the ATP pool causes a reduction in the rate of synthesis of glucose 6-phosphate. The pool of glucose 6-phosphate becomes depleted and limits both the activity of phosphoglucoisomerase, which leads to the depletion of fructose 6-phosphate, and the activity of glucose-6-phosphate dehydrogenase, which completes the cycle by further lowering the influx to the pentose pathway and causing the accumulation of NADP. DISCUSSION We have previously taken a model of metabolism in human red blood cells through three stages of assessment and refinement (18,19) to produce the model used in this paper (Model III). With this model we have addressed two issues. First, we have examined the propagation of biochemical signals within this metabolic network. This requires an integrated systems approach, and the relevant concepts in this context are logarithmic gains in concentrations and fluxes. Second, we have examined the accuracy of representation. Here we have used a comparative approach involving alternative kinetic models to identify areas of metabolism that might not be well represented. The results of these studies point to specific areas of 3 The timing information refers to CPU time obtained from the standard function call time() in ANSI C. The C source code is converted to the corresponding object code with the Think C 5.02 compiler and then run on a Macintosh IIx computer. The standard Newton-Raphson method of finding the roots for a system of nonlinear algebraic equations was used to obtain the steady-state solution for the Michaelis-Menten model.

FIG. 5. Error of representation based on a comparison of power-law and Michaelis-Menten formalisms.
Independent variables are changed by Ϫ2%, and the Michaelis-Menten representation is considered as the reference system. The three-dimensional plot displays the magnitudes of the differences as a function of metabolite concentrations X i and of independent variables X j . Differences with magnitude less than 2.0E-5 are shown as black squares. The two-dimensional projection on the right gives the average of the magnitudes for a particular metabolite concentrations X i summed over all the independent variables X j . The two-dimensional projection on the left gives the average of the magnitudes for an independent variable X j summed over all the metabolite concentrations X i . erythrocyte metabolism that are in need of further experimental investigation.
Signal Propagation-The analysis of signal propagation identified metabolite concentrations whose response to a change of input signal was either highly attenuated or highly amplified (Fig. 2). Metabolite pools with a highly attenuated response to change in most independent variables include those for gluconolactone 6-phosphate X 1 , oxidized glutathione X 4 , lactate X 20 , NADH X 22 , and pyruvate X 19 . The first two are substrates of nearly irreversible reactions located within a pathway that is rate-limited by its initial reaction. The last must be in equilibrium with the extracellular pool of pyruvate to maintain proper coupling of the dehydrogenase reactions in glycolysis. The other two reflect primarily the influences of their immediate environment. Metabolite pools that exhibit highly amplified responses include those for gluconate 6-phosphate X 2 , nicotinamide adenine dinucleotide phosphate X 3 , sedoheptulose 7-phosphate X 8 , glucose 6-phosphate X 10 , fructose 6-phosphate X 11 , 5-phosphoribosyl-1-pyrophosphate X 29 , and inosine X 31 . These are metabolites at key points in the pentose phosphate pathway and nucleotide metabolism that are responding to and are part of the redistribution of flux between these two subsystems.
A similar pattern of systemic responses was identified among the fluxes of the system (Figs. 3 and 4). The flux through the proximal portion of the pentose phosphate pathway (e.g. V 1 ) and the fluxes through the pools of sodium V 24 , potassium V 23 , and 2,3-diphosphoglycerate V 21 have highly attenuated responses to nearly all of the independent variables. These fluxes are fixed by a rate-determining step or, in the case of potassium, are isolated from other influences by saturation of the key enzyme through which these influences would otherwise be communicated. On the other hand, fluxes that exhibit highly amplified responses include those through the pools of nucleotide metabolism (V 26 through V 28 , and V 31 through V 33 ) the early portion of glycolysis (V 10 and V 11 ), and the lower portion of the pentose phosphate pathway (V 6 ). Again, this manifests the redirection of flux into or out of nucleotide metabolism.
The analysis of signal propagation also identified independent variables (input signals) whose influence on the metabolite concentrations and fluxes of the system is either very little or very great. One half of the independent variables have almost no influence on the metabolite concentrations or fluxes of the system. These include enzyme levels for fast reactions operating near equilibrium and for reactions that follow a rate-limiting step in a pathway. On the other hand, three independent variables are highly influential in affecting the metabolite concentrations and fluxes of the system: the activities of the enzymes glutathione peroxidase X 34 and transaldolase X 43 , and the concentration of total glutathione X 80 (Figs. 2, 3, and 4). Each of these three variables has a major influence on the distribution of flux between the pentose phosphate pathway and nucleotide metabolism.
We have not ruled out a physiological role for this diversion of flux in the context of the whole cell. If this behavior were to be part of a normal response, then the maintenance of ATP levels in the cell would require a minimal flux through glutathione peroxidase. Attempts to maintain ATP levels during blood storage by eliminating the oxidative stress responsible for this flux (34) might then be counterproductive (18). On the other hand, in animals with a dietary deficiency for selenium, large changes in the activity of glutathione peroxidase (a selenium-dependent enzyme) occur without a major disruption of erythrocyte metabolism (35); whereas in the model, small changes (9%) in the activity of glutathione peroxidase cause diversions leading to the loss of the nominal steady state, which suggests that there are still unsolved issues here.
Accuracy of Representation-Some investigators have questioned the use of the power-law formalism because of perceived limitations in accuracy beyond the steady state (e.g. see Ref. 36). It is clear from these accounts that the authors are not questioning the canonical character of the power-law formalism when used as a fundamental representation or as a recast representation (37,38), contexts in which the question of accuracy is less relevant. The concerns are restricted to the use of the power-law formalism as a local representation, and they appear to be based on previous experience with local linear approximations, which are of more limited accuracy. However, the power-law representation is a nonlinear approximation that generally has a much wider range of accurate representation. Evidence that bears on this question can be summarized as follows. Worst case scenarios for single enzymes suggest that the range of accurate representation is greater than a 2-fold change in the concentration of reactants (39,40). Studies with simple systems of increasing complexity have shown that on average the range of accurate representation increases with FIG. 6. Range of accurate representation based on a comparison of power-law and Michaelis-Menten formalisms. The level of the independent variable glutathione peroxidase is considered 100% at the nominal steady state. As this independent variable is being reduced the Michaelis-Menten model undergoes drastic changes in the concentration of several metabolites. After a reduction of only about 9% there is an explosive accumulation of NADP ϩ , while the concentrations of ATP, glucose 6-phosphate, and fructose 6-phosphate plummet; beyond this point the Michaelis-Menten model fails to exhibit a steady state. See ''Results'' for further discussion. the size of the system being considered from 5-fold (41) to 20-fold (42). From another perspective, examination of the actual ranges of variation for a large number of biochemical variables in humans suggests that on average the range is about 4-fold (42), which is no larger than the range of accurate local representation observed with the simpler systems. Thus, the evidence suggests that the range of accurate local representation is sufficiently broad to be of use for many biochemical systems, and, in those situations where it might not be, one can increase the accuracy as necessary by using either the fundamental or recast representations within the power-law formalism.
We have used a Michaelis-Menten model of metabolism in human erythrocytes as a reference to test the predicted accuracy of the local power-law representation in the context of a larger system. Our previous analysis of parameter sensitivities led us to predict (19) that Model III would accurately represent the Michaelis-Menten model only over a limited range of variation in its concentration variables. In this paper, the question of accuracy was examined directly by comparing alternative representations involving the Michaelis-Menten and powerlaw formalisms. We first determined the error in response to a small change in the independent variables as a function of independent and dependent concentration variables (Fig. 5). On average the discrepancies between these alternative representations is small (0.0104%), which suggests that much of the system is well represented. However, the discrepancies that do exist are large. The greatest errors occur in the pentose phosphate pathway and in nucleotide metabolism, and these occur in response to changes in the activities of glutathione peroxidase X 34 and transaldolase X 43 . Furthermore, in determining the range of accurate representation we discovered that the large discrepancies are due to the behavior of the Michaelis-Menten model, which fails to exhibit a nominal steady state when the activity of glutathione peroxidase is decreased by as little as 9% (Fig. 6). Thus, the narrow range of accurate representation that was predicted on the basis of parameter sensitivities (19) has been confirmed here and associated with the lack of robustness of the reference model.
The comparative approach has demonstrated the ability of the theory to predict areas of problematic representation and has contributed to the diagnosis of problems in the Michaelis-Menten model. However, because of the inadequacies in the reference model, this approach could not yield useful information about the actual range of variation over which the biological system might be accurately represented by the power-law formalism.
Other Comparative Studies-The best test for the accuracy of local representation would be provided by a detailed comparison of values predicted by theory and values experimentally observed in situ for a large realistic network. Unfortunately, we still lack appropriate data of this type on which the range of accurate representation might be critically evaluated. Instead, mathematical models, which in turn are based on extensive data obtained in vitro, have been used as an indirect indicator of concentrations and fluxes in situ. Comparative studies using this indirect approach have been conducted recently with three different systems: the tricarboxylic acid cycle in Dictyostelium discoideum (20,24,32,37,43), glycolytic fermentation in Saccharomyces cerevisiae (33,44,45), and metabolism in human erythrocytes (18,19, this paper). In each of these cases no definitive answer regarding local accuracy has been forthcoming, not because of limitations in the power-law formalism, but because the conventional representations based on kinetic data obtained from experiments in vitro have exhibited anomalies. Another indication of concern with kinetic representations based on data obtained in vitro is the fact that ad hoc adjust-ments in the data often are needed to make these models consistent with data obtained in vivo (45)(46)(47)(48). Nevertheless, in each of these three cases analysis using the local representation led to the diagnoses of problems and to suggestions for areas in need of further study; in two cases it led to modifications that improved the model. These results, taken together with recent advances in our understanding of the intracellular milieu (49 -51), suggest that more definitive assessments of accuracy will require kinetic data obtained under conditions that more accurately reflect the natural conditions in vivo. Attempts to develop new theoretical approaches to this problem recently have been described (52,53). However, the inability to measure important biochemical signals dynamically in situ represents the most serious bottleneck at this time (17).
Critical Examination of the Theory-Given the absence of appropriate data with which to make direct tests, and the limitations of indirect tests using mathematical models based on data obtained in vitro to infer the values of concentrations and fluxes in situ, one should be appropriately critical of any theory. Wright and Field (54) have recently presented criticisms (which were repeated in Albe and Wright (55)) of the approach used in this paper, but the criticisms in this case are without foundation (56). Wright and Field specifically criticized the interpretation, methodology, and theory presented in the papers by Shiraishi and Savageau (20,24,32,37,43). They claim that the results of Albe and Wright (57) and Shiraishi and Savageau (20) should be interpreted as similar, whereas Shiraishi and Savageau (20,32) had shown them to be different. Examination of the original data shows that Wright and Field are mistaken because they have not considered comparable data from the two studies (56). Wright and Field also claim that the method used by Shiraishi and Savageau (20) for calculating logarithmic gains is less accurate than the method used by Albe and Wright (57). Shiraishi and Savageau (20) use an analytical method that gives the exact result based on infinitesimal changes about a given steady state; Albe and Wright (57) use an empirical method that gives an approximate result based on finite changes between two steady states. The empirical results can be made to closely approximate the analytical results, provided the finite differences are made sufficiently small and one is dealing with valid steady states (20; see also "Methods"). However, Albe and Wright (57) were not dealing with valid steady states, and this is the reason for the different results they obtained (32,56). Finally, Wright and Field claim that the theory underlying the work of Albe and Wright (57) is different from that underlying the work of Shiraishi and Savageau (20). It has been rigorously demonstrated (26,58) that the theories are identical within the domain considered by Albe and Wright. Albe and Wright (57) should have obtained exactly the same results as Shiraishi and Savageau (20); the fact that they did not is the result of failure to apply the theory correctly (32). Although the criticisms of Wright and Field are unwarranted, the approach we have used in this paper should continue to be criticized by testing its predictions both for their qualitative relevance and, eventually, for their quantitative accuracy in situ. Such testing will require further experimental work.
Areas Requiring Further Investigation-There are several areas in which experimental work is needed to improve our model of metabolism in human erythrocytes. When the results in this paper are combined with those of previous work (19) one finds a consistent pattern pointing to problems with the representation of nucleotide metabolism. (a) The predicted and observed steady-state values for many of the intermediates of nucleotide metabolism differ significantly (19). (b) The parameter sensitivities associated with nucleotide metabolism are the largest of the system (19). (c) The failure to maintain a steady state in the face of minor changes in independent variables is associated with the depletion of intermediates of nucleotide metabolism (this paper). These problems seem to result from a highly amplified diversion of flux at a key branch point between the pentose phosphate pathway and nucleotide metabolism. Earlier studies of the tricarboxylic acid cycle in D. discoideum have shown that the distribution of flux at a branch point can be critical (43) and that explosive increases in metabolite concentrations can occur when the distribution is not flexibly regulated (20). The physiological solution to this set of problems in the erythrocyte model is not obvious. The scanning algorithm described elsewhere (19) could be adapted to suggest putative processes or regulatory interactions that have the potential to improve the distribution of flux between the pentose phosphate pathway and nucleotide metabolism.
Another issue in need of further attention is the kinetic characterization of the fast reactions that operate near equilibrium. It is frequently assumed that fast reactions near equilibrium can be described by simple mass-action kinetics and the choice of large values for the rate constants. In principle, the results should be independent of the choice if the values are sufficiently large and their ratio is equal to the equilibrium constant for the reaction. However, when these values are not made sufficiently large, then an arbitrary choice of values can have significant effects on the state of the system (18). We tried to eliminate these arbitrary effects by assigning the same large number 10 5 for all forward rate constants. We could not use larger numbers because then the Michaelis-Menten model failed to maintain a steady state. After having made various refinements in the model we have gone back and reexamined this issue. We assumed the same value for the forward rate constants of all of these fast reactions, varied this value over the range 10 5 to 10 8 , and found the steady state to be largely unchanged. Nevertheless, this procedure remains questionable. The choice of very large values for the rate constants forces reactions closer to equilibrium than they may actually be in the organism, and this introduces error into the model. On the other hand, too low a value may not bring the reactions close enough to equilibrium, and this also introduces error. Choosing just the right set of values in an ad hoc fashion to make the model fit the data is unacceptable. The only satisfactory solution to this problem is better data and a more complete kinetic characterization for these fast reactions.
Conclusions-As emphasized elsewhere (19) and throughout this study, model formulation, assessment, and refinement are an iterative process, particularly in the case of large complex systems. Our purpose has been to provide further examination of the utility of biochemical systems theory by carrying out a few such iterations involving detailed analysis of a relatively large and realistic system: metabolism in human red blood cells. This system was chosen for two reasons. First, a substantial amount of experiment data concerning metabolism in human red blood cells already exists. Second, these data already have been synthesized into kinetic models of the integrated system. In the iterative process of applying the theory to this system (18) we encountered inconsistencies, omissions, and erroneous assumptions in the literature, which is to be expected for such a complex system. We also found that our earlier models were neither stable nor robust (19), systemic properties whose elucidation is perhaps more subtle. Our previous work on this system can be viewed as a case study illustrating how these types of difficulties might be overcome during the construction, assessment, and refinement of an integrated model. Although the refinements we made were modest and left substantial room for further improvement, this case study nevertheless identified several biochemically relevant features of the model that had been overlooked. It also pointed to nucleotide metabolism as the area most in need of further experimental study. This conclusion has been further reinforced here by the finding that Model III exhibits a highly amplified diversion of flux between the pathways of pentose phosphate and nucleotide metabolism in response to small changes in the rate of glutathione oxidation. Thus, biochemical system theory can be seen to provide a rigorous framework and systematic methodology that allows one not only to evaluate the quality of a model, to diagnose its deficiencies and to predict modifications that can improve the model, but also to explore the integrated behavior of a complex biological system. V Ϫ11 ϭ 1.52Dϩ4 X 7 Ϫ5.48DϪ5