INTRODUCTION

TOP ABSTRACT INTRODUCTION MODEL DESCRIPTION AND SOLUTION... IDENTIFYING RELATIONSHIPS WITH... ANALYTICALLY IDENTIFYING THE... FORMULATION OF PHYSIOLOGICALLY... DYNAMIC RESPONSE OF RBC... DISCUSSION REFERENCES

Biology has progressed with the detailed study of various "model" experimental systems, such as Drosophilia melanogaster, Saccharomyces cerevisiae, and Escherichia coli. It is likely that the same approach will be used for developing methodologies for the computational analysis of complex biological processes. Metabolism represents a logical starting point because it has been studied through mathematical modeling and computer simulation for a number of years (Bailey, 1998; Fell, 1996; Heinrich and Schuster, 1996; Reich and Sel'kov, 1981). The red cell has played a special role in the development of mathematical models of metabolism given its relative simplicity and the detailed knowledge about its molecular components.

Mathematical models of red blood cell (RBC) metabolism have been studied since the early 1970s (Brumen and Heinrich, 1984; Heinrich et al., 1977; Heinrich and Rapoport, 1973; Holzhutter et al., 1985; Rapoport et al., 1981; Schauer et al., 1981; Schuster et al., 1988; Werner and Heinrich, 1985). The models have steadily grown in scope leading to comprehensive RBC metabolic models in the late 1980s and 1990s (Holzhutter et al., 1990; Jacobasch and Rapoport, 1996; Joshi and Palsson, 1989a; Mulquiney and Kuchel, 1999; Rae et al., 1990; Schuster et al., 1988, 1989, 1990).

Herein, we use simulations of a RBC metabolic model (Jamshidi et al., 2001) to develop a framework for extracting the dynamic characteristics of a biological network. This RBC metabolic model is comprised of 44 enzymatic reactions and membrane transport systems and 34 metabolites and ions. Using this RBC model, we show how dynamic phase planes, statistical time-lagged correlation analysis, and temporal decomposition can be used to decipher the relationship between the biochemical details and the overall metabolic function of the hRBC. The resulting analyses are then compared with existing biochemical knowledge.

	MODEL DESCRIPTION AND SOLUTION METHODS

TOP ABSTRACT INTRODUCTION MODEL DESCRIPTION AND SOLUTION... IDENTIFYING RELATIONSHIPS WITH... ANALYTICALLY IDENTIFYING THE... FORMULATION OF PHYSIOLOGICALLY... DYNAMIC RESPONSE OF RBC... DISCUSSION REFERENCES

Red blood cell metabolic model

The RBC metabolic model includes glycolysis, the Rapoport-Leubering shunt, the pentose phosphate pathway, nucleotide metabolism reactions, the sodium/potassium pump, and other membrane transport processes (see network representation in Fig. 1). Dynamic mass balance equations for the 34 metabolites constituting this metabolic system are represented as:

(1)

in which the vector v consists of the specific enzymatic reaction rates, the vector x is the specific concentrations of the various metabolites, and the matrix S is the stoichiometric matrix for the enzymatic reactions. The specific metabolites comprising the vector x and the enzymes comprising the vector v are summarized in Table 1.

View larger version (50K): [in this window] [in a new window]   FIGURE 1   The details of the hRBC metabolic map used for this study. It includes all the major metabolic pathways and events in the red cell. Nomenclature is defined in Table 1.

Kinetic expressions used for the enzymatic reaction rates (v = f(x)) in the RBC were based on rate equations previously assembled from the literature (Joshi and Palsson, 1990a; Lee and Palsson, 1991). Based on the rate laws, detailed dynamic mass balance equations were formulated.

Numerical solutions

The solutions to the model equations were obtained using Mathematica 4.0 (Wolfram Research, Champaign, IL) and MATLAB 6.0 (Math Works Inc., Natick, MA). The code can be downloaded from http://epicurus.che.udel.edu.

DYNAMIC PHASE PLANES TO IDENTIFY THE METABOLIC POOLS

The metabolic network in the RBC, although relatively simple compared with other cells, is sufficiently complex to make it difficult to conceptualize the interactions of the 34 concentration variables. Therefore, a reduced set of metabolites that can capture the behavior of the RBC metabolic network while still maintaining biochemical meaning is desired. Dynamic phase planes provide one avenue for determining aggregates of metabolites that help accomplish this.

General features of dynamic phase plane trajectories

Dynamic phase planes are traces of two metabolite concentrations, parameterized with respect to time, in an x, y plane. Thus, the trajectories in the diagram move chronologically from an initial state to a final state (e.g., Hubbard and West, 1991). The dynamic phase plane analysis can point to key relationships between metabolites of the metabolic network. For example, if a system is dynamically stable, the dynamic trajectory in a phase plane will converge to a single point in the plane, known as an attracting fixed point. A stable steady-state point biologically represents homeostasis. Conversely, if the system is unstable then the trajectories will not approach a fixed point but diverge away from it. The red cell model contains both stable and unstable steady states (Edwards and Palsson, 2000; Heinrich et al., 1977; Joshi and Palsson, 1990b). This study is based on perturbations to the biochemically relevant, stable steady state whose initial conditions are specified in the program code available at the web site.

Specific features of phase plane trajectories for step responses

A trajectory in the phase plane may indicate the presence of one or more general dynamic features. Namely, the shapes of the trajectories contain significant information about the dynamics of the system and its regulation. Some important features of phase plane trajectories for step responses are as follows (Fig. 2).

View larger version (14K): [in this window] [in a new window]   FIGURE 2   General features of phase planes. Dynamic phase planes are formed by graphing the time dependent concentrations of two metabolites (A and B) against one another. Phase planes for a step response can have certain characteristic features. (1) Conservation relationship; (2) pair of concentrations that could be in equilibrium with one another; (3) motion of the two metabolites dynamically independent of one another; (4) closed loop traces representing either a periodic motion or a return to the original steady state.

1) When the trajectory has a negative slope, it indicates that one metabolite is increasing while the other is decreasing. The metabolite concentrations are moving on the same time scales but in opposite directions, that is, one is consumed while the other is produced. This feature might represent the substrate concentration versus the product concentration of a given reaction.

2) When a trajectory in the phase plane between two metabolites is a straight line with a positive slope, it means that the two metabolites are moving in tandem; i.e., as one increases so does the other and vice-versa. This feature is observed when two or more metabolites move on the same time scales and are in equilibrium with one another. Such behavior helps define the aggregate concentration variables.

3) When a trajectory is vertical or horizontal, it indicates that one of the metabolite concentrations is changing while the other remains constant. This feature implies either that the motions of the metabolite concentrations during the trajectory are independent of one another or that the dynamic motions of the metabolite concentrations progress on different characteristic time scales.

4) When a trajectory forms a closed loop, it implies one of two possibilities. If the system never converges to a steady state over time, then it implies oscillatory motion. On the other hand, if the orbit begins at one point, moves away from it, then returns to the same point after a sufficiently long time interval, then it implies that a change in another variable in the system forced it away from its steady state temporarily, but it returned to the original steady state.

The qualitative characteristics of dynamic phase planes can provide insight into the inner workings of the RBC metabolic model. A trajectory may have more than one of these basic features. For instance, there can be a fast independent motion (i.e., a horizontal phase plane trajectory) followed by a line with a positive slope after an equilibrium state has been reached. Such features demonstrate the importance of timescale separation and the need for analysis procedures that allow one to study behavior on multiple time scales.

Analysis of the red blood cell metabolic network

The time scales in the red cell metabolic model vary over many orders of magnitude (see time constants in Table 2). There are several rapid reversible reactions in red cell metabolism, which quickly equilibrate between the reactants and products. The rapid equilibration leads to straight lines with positive slope in the phase plane (illustrated in Fig. 2) where the slope is the equilibrium constant of the reaction. These equilibria can take place on many different time scales and thus one would expect that the system could be decomposed in time to describe the dynamic events that take place on the various timescales. To visually inspect the relationship between the various metabolites, we have plotted all pair-wise phase planes in a matrix format (Figs. 3-5). The matrix of phase planes is symmetric; therefore, we have entered the correlation coefficients corresponding to a particular phase plane in the transpose position in the matrix; thus, providing a quantitative description of the phase plane linearity between the two variables. When the absolute value of the correlation coefficient is greater than 0.85, the slope of the resulting line is also included. For metabolites in equilibrium, the slope indicates the equilibrium constant between the two metabolites.

View larger version (97K): [in this window] [in a new window]   FIGURE 3   Matrix of pair-wise phase planes on a network wide scale. Each plot in the matrix corresponds to the phase planes formed by the metabolite in the corresponding row and column. The matrix is symmetric so that in the transpose position we show the correlation coefficient (r) for the corresponding trace (when  r  > 0.85, we also show the slope of the line, bottom number). The response is a composite of the entire dynamic response to a nonspecific, 0th order ATP load (VATP = 0.3 mM/h). Points on the phase planes are color coded to identify the location at specific times: initial point (red), 0.1 (green), 1 (blue), and 10 (yellow) h. Images of the various time scales studied can be found at http://epicurus.che.udel.edu/.

View larger version (101K): [in this window] [in a new window]   FIGURE 4   Matrix of pair-wise phase planes on a network wide scale for the response to a nonspecific, 0th order redox load (VNADPH = 4 mM/h). Description as in Fig. 3.

View larger version (41K): [in this window] [in a new window]   FIGURE 5   Matrix of pair-wise phase planes for the pooled variables. (A) Dynamically defined pools and (B) capacity, occupancy, and charge pools for step changes in ATP load and redox load. Description as in Fig. 3. Regardless of the disturbance, adenosine occupancy (AO) moves independently of capacity (AC) at shorter time scales. However, over longer time scales they are correlated. Glycolytic occupancy (GO) is highly coupled with capacity (GC) over both time scales. Redox occupancy (RO) is coupled with capacity (RC) for the ATP load, but for the redox load, the two move independently.

hRBC metabolism needs to respond to nonspecific loads corresponding to environmental changes during its normal functions. Therefore, the multitime scale response of the red cell metabolic network to an ATP load of 0.3 mM/h (any of a variety of environmental changes causing the cell to become less efficient with respect to ATP, implemented as a 0^th order additional ATPase activity term) is shown in Fig. 3. It should be noted that additional loads can be placed on the RBC, such as a redox load (Fig. 4). Points on the phase planes are color coded to visualize specific time points in the step response. In each case, the phase plane starts at the red point (0 h), proceeds through the green point (0.1 h), the blue point (1 h), and the yellow point (10 h), and ends at the black point (300 h). Additional images corresponding to the shorter time regions are available for download at http://epicurus.che.udel.edu. The pair-wise correlation coefficients shown represent the entire dynamic range. This data representation allows rapid analysis and interpretation of the dynamic behavior, as is immediately evident from the clusters of similar traces in the matrix of phase planes.

Dynamics faster than 0.1 h

Inspection of the phase planes visually reveals the formation of numerous metabolic pools. We have defined pools based on the value of the correlation coefficient, and metabolites with a correlation coefficient greater than 0.95 were pooled together. Additionally, we required that the metabolites satisfy this criterion under an ATP and redox load with a predicted equilibrium constant (the slopes of the linear relationships) within 20%.

The trace of G6P and F6P has a correlation coefficient of 1.0 and a slope that corresponds to the equilibrium constant of the PGI reaction (0.41). Similarly FDP, DHAP, and GAP all move together, as do PG3, PG2, and PEP (the PGs), ADP, AMP, and ADO from adenosine metabolism, RU5P, R5P, and X5P from pentose metabolism, and NADPH and GSH from redox potential. Thus, pools of metabolite concentrations are forming based on the rapid equilibria between them. Studying the dynamics in response to an ATP and redox load on this time scale leads us to a simplified representation of the reaction network as shown in Fig. 6 A.

View larger version (52K): [in this window] [in a new window]   FIGURE 6   hRBC metabolic map showing the simplification of the dynamic interpretation of the metabolic map resulting from metabolite pooling that occur over each dynamic region: (A) 0 to 0.1 h, (B) 0.1 to 1.0 h, (C) 1.0 to 10 h, and (D) 10 to 300 h. Each identified pool is shown with a different color coding.

Observing the dynamic phase plane behavior on additional time scales reveals several interesting phenomena. For example, with an ATP load, the G6P/F6P pool moves opposite the FDP/DHAP/GAP pool on a fast time scale (e.g., within the first 0.1 h); however, on slower time scales these metabolites form a larger pool. The time scale separation of the G6P/F6P pool and the FDP/DHAP/GAP pool can be described under the ATP load. When the ATP load is imposed, ATP levels will decrease, and less G6P will be produced, resulting in a decrease in the G6P/F6P pool. At the same time, there is a corresponding increase in AMP levels, leading to an increased activation of PFK. Thus, the metabolic network will consume F6P resulting in an increase in the FDP/DHAP/GAP pool. Therefore, the response to an ATP load is for the energy-rich glycolytic intermediates to increase transiently in concentration and then fall as the system adapts and is able to regenerate the ATP lost due to the additional load. This demonstrates the separation of these pools on a short time scale. However, on a longer time scale, these pools move together. All of this behavior is observed in the phase planes for the ATP load, but some differences are noticed using the redox load. Thus, to achieve a more accurate approximation of dynamically independent metabolite pools, one should observe the consensus dynamics observed upon application of different kinds of loads.

Dynamics between 0.1 and 1.0 h

The consensus of the dynamic phase planes on this time scale show that some of the pools change substantially, as shown in Fig. 6 B. For example, over this intermediate time scale, FDP no longer moves with DHAP and GAP (which still have an equilibrium constant of 0.05). There also is a coupling between RU5P, R5P, and X5P with R1P that did not exist on the shorter time scale. Further, there appears to be substantial coupling between the PGs and ADO, AMP, and ADP (although not quite sufficient to collapse all six metabolites into one pool). ATP also clearly moves opposite ADP, AMP, and ADO in a very predictable manner.

Dynamics on the remaining time scales

The dynamic phase planes indicate that, for each individual load on the 1.0- to 10-h time scale, G6P, F6P, FDP, DHAP, GAP, and DPG13 all are highly correlated. For a given load, these six metabolites all collapse into one pool, called the glycolytic pool. However, the behavior of the glycolytic pool is dependent on which load is imposed. From a biochemical standpoint, the high degree of correlation here makes sense, as does the dependence on the two separate inputs. Changing the redox load causes a change in the flux through the pentose phosphate pathway. This change in flux will affect the resulting steady-state relationships between G6P/F6P and FDP/DHAP/GAP/DPG13. Again, the differences between the two loads can provide insight into the underlying biochemical system. These results are consistent with previously recognized physiologically significant pools such as such as those in the glycolytic pathway (Heinrich et al., 1977; Palsson et al., 1987; Rapoport et al., 1976).

As seen in Fig. 6 C, both disturbances show that the PGs and PYR collapse together, and the similarities between the PGs and the adenosine metabolites (now consisting of ADO, INO, HX, and AMP) have disappeared. On the 10- to 300-h time scale the adenosine pool returns, containing ADO, AMP, ADP, ATP, INO, and HX (shown in Fig. 6 D). Finally, if the entire dynamic window is considered at once (as shown in Figs. 3 and 4), the following pools result: (G6P, F6P, GO6P), (DHAP, GAP, FDP, DPG13), (PG3, PG2, PEP), (RU5P, R5P, X5P, S7P, R1P), and (AMP, ADO). It is interesting to note that some of the pools that appear over the entire dynamic time range do not appear for any subset of the dynamic range explored (e.g., S7P grouped with the RU5P pool).

The key point of this discussion is to realize that the time scale of relevance for pooling depends on two considerations: 1) on what time scale does interesting behavior occur and 2) which of those predictions are required to accurately represent the overall system behavior under the conditions of interest. Phase plane analysis is capable of identifying equilibrium and opposing reactions when no delay exists between two metabolites. However, if two metabolites move in tandem with a delay, the phase plane analysis will miss this feature. One method for overcoming this would be to calculate new phase planes for all the metabolites with all possible delay combinations for the other metabolites. For obvious reasons, carrying out such an analysis is impractical. Statistical correlation metrics provide an alternative, more tractable approach.

	IDENTIFYING RELATIONSHIPS WITH CORRELATION METRICS

TOP ABSTRACT INTRODUCTION MODEL DESCRIPTION AND SOLUTION... IDENTIFYING RELATIONSHIPS WITH... ANALYTICALLY IDENTIFYING THE... FORMULATION OF PHYSIOLOGICALLY... DYNAMIC RESPONSE OF RBC... DISCUSSION REFERENCES

Like the visual interpretation of the data using the dynamic phase plane analysis, rigorous statistical methods can be applied to identify the relations between metabolites in a metabolic network (Arkin and Ross, 1995; Arkin et al., 1997). To perform a statistical analysis, we constructed the time-lagged correlation coefficient matrix, R:

(2)

in which S_ii is the variance of the dataset x_i(k) and S_ij() is the time-lagged covariance between the discrete, uniformly sampled datasets x_i(k) and x_j(k + ), determined as:

(3)

in which n is the number of data points in the series. Values of r_ij() range from 1 to 1, in which a value of 1 indicates the two datasets are perfectly anticorrelated with a delay of  time steps, a value of 1 indicates the two datasets are perfectly correlated with a delay of  time steps, and a value of 0 indicates that the two datasets are statistically uncorrelated at a delay of  time steps. The time-lagged correlation coefficient matrix can be used to generate a new variable, c_ij, which defines the strength of maximal correlation observed between two datasets:

(4)

Thus, when two datasets are perfectly correlated or anticorrelated for any sampling lag, , c_ij will approach a value of one. If the two datasets are perfectly random (i.e., absolutely no correlation), the maximal correlation metric will approach a value of zero.

For demonstrative purposes, Fig. 7 shows the maximal correlations calculated between each of five key intermediates (G6P, GAP, PEP, R5P, and ATP) and the other 34 metabolites found in the RBC. Calculated maximal correlations are shown for each of five distinct perturbations: glucose flux (µ = 5, ² = 2, ), NADH load (µ = 0.1, ² = 0.033, +), NADPH load (µ = 0.1, ² = 0.033, ), ATP load (µ = 0.1, ² = 0.033, ), and all four varied simultaneously (). In each case, a truncated random Gaussian perturbation (negative values were truncated to zero) was introduced every 10 h, and concentrations were sampled every 5 h for 10,000 h. Correlations were explored over a range of up to 250 sampling times. The results show a strong relation between the applied input(s) and the observed maximal correlations. Since a given load does not necessarily perturb every pathway of the metabolic network, the resulting dynamics would (incorrectly) indicate equilibrium. This incorrect identification of equilibrium pools results from leaving portions of the network unperturbed and is seen in the case of the NADH load (+). The NADH load plot shows that G6P, F6P, FDP, DHAP, GAP, DPG13, PG3, PG2, PEP, PYR, LAC, NADH, GO6P, NADPH, GSH, RU5P, R5P, X5P, S7P, E4P, and R1P all appear to be highly correlated. This indicates that all effects caused by the NADH load equally affect glycolysis and the pentose phosphate pathway. However, this is not the case for a general perturbation to the system.

View larger version (47K): [in this window] [in a new window]   FIGURE 7   Maximal time-lagged correlations. Calculated maximal correlations for each of five distinct perturbations: the glucose flux into the system (); an NADH load (+); an NADPH load (); an ATP load (); and all four inputs of ATP load, NADH load, redox load, and glucose flux were varied simultaneously (). For each disturbance, a random perturbation was introduced every 10 h, and samples were collected every 5 h.

Several differences are observed between the NADH load and the NADPH (redox) load. For example, while the PGs are all very close for every input sequence, their separation from the rest of glycolysis and from the pentose pathway is highly dependent on the input. Unlike the NADH load, the NADPH load demonstrates a difference in the dynamics of the PGs and the rest of the glycolytic intermediates. Further, an NADPH load seems to indicate that ADP is very closely related to the PGs, which is misleading. This clearly demonstrates the key role of the perturbations for determining the types and quality of interactions observed and indicates that further work needs to be done to determine combinations of inputs that sufficiently excite metabolic pathways.

Recall that this figure was generated from a discrete data set obtained by sampling every 5 h and perturbing every 10 h. Thus, the data set only contains information about frequencies of response that occur within the sampling time period. Any dynamics with a frequency of greater than ~0.2 h¹ will be missed by this sampling frequency. Based on the sampling times used for this analysis, the dynamics under consideration correspond to the 1- to 10-h time scale described under the phase plane analysis.

It is interesting to note certain pools define themselves by appearing together on this time scale regardless of the input used. These pools agree well with the phase plane analysis. From the pentose phosphate pathway, Fig. 7 shows that R5P is closely correlated to RU5P, X5P, and R1P for all inputs. As with the phase planes, these species always move in a concerted manner on the time scale sampled. The results show that measuring samples every 5 h indicates that ATP is uncorrelated with any other input, as predicted by the phase plane analysis. Similar to the phase plane analysis, the maximal correlation metrics need to be calculated for each of several different time scales to be totally effective.

As mentioned, there are several similarities between the statistical analysis and the phase plane analysis, but several differences also exist. Like the phase plane analysis, the statistical analysis identified G6P and F6P and FDP, DHAP, and GAP as pools. However, unlike the phase plane analysis, the statistical analysis determined that these pools remain separated over this time scale for the combined input. Further, the statistical analysis identified separation between PYR and the PG pool over this time scale. Finally, the statistical analysis identified PRPP and ADE as highly correlated (anticorrelated in this case), whereas the phase plane analysis did not.

The key points of this discussion are to realize that the insights from a statistical correlation analysis depend on: 1) what time scale the data was sampled and 2) whether the inputs used adequately perturb the system to excite the dynamics appropriately. The statistical analysis is capable of identifying time-lagged correlations and thus, for systems with delays, is able to provide additional insight beyond that gained by the phase plane analysis. Because the statistical analysis uses a much richer input sequence to perturb the system than the phase plane analysis, it is capable of better breaking down the dynamics and looking for places where the predictions of the phase plane may fall apart. One major drawback of the approach is that, as a result of the richer inputs used, the time evolutions of the various profiles can be more difficult to interpret from a physical standpoint.

Both of the statistical and phase plane techniques use perturbations to the metabolic network to identify metabolites that move together. The two methods compliment each other rather well with the phase plane analysis providing a first pass that is easily interpreted from a physical standpoint and the statistical analysis filling in additional details. However, it is also possible to use an analytical technique to look for metabolites that have similar dynamic behavior.

	ANALYTICALLY IDENTIFYING THE METABOLITE RELATIONS

TOP ABSTRACT INTRODUCTION MODEL DESCRIPTION AND SOLUTION... IDENTIFYING RELATIONSHIPS WITH... ANALYTICALLY IDENTIFYING THE... FORMULATION OF PHYSIOLOGICALLY... DYNAMIC RESPONSE OF RBC... DISCUSSION REFERENCES

Above we have discussed two approaches for identifying the relation between any two metabolites from data sets. However, if a kinetic model for the network exists, an analytical procedure can be used to decompose the metabolic network based on intrinsic time scales. The temporal decomposition of metabolic systems using modal analysis has been previously described (Heinrich and Schuster, 1996; Palsson, 1987; Palsson and Lightfoot, 1984). In brief, modal analysis involves linearizing the dynamic mass balance equations around a reference point (x_ref),

(5)

The linear system of equations described by Eq. 5 can be rewritten as a dynamically independent set of equations by diagonalizing (J = MM¹) the Jacobian matrix (J) (see Strang, 1988). A set of dynamically independent variables or modes (m), which move on distinct time scales, can be defined as:

(6)

in which M¹ is the modal matrix and the vector m contains the dynamically independent modes. This results in the transformation of Eq. 5 to:

(7)

The dynamic solution for each mode is m_i(t) = m_i,0exp(t/_i) in which _i, the time constants for each mode, are defined as the negative inverse of the eigenvalues of the Jacobian matrix (the elements of the diagonal  matrix). The time constants and the modal matrix for the RBC model used here are shown in Table 2.

Modal analysis can be used to transform the original concentration variables into linearly independent combination variables called modes. These independent modes can be used to decipher the dynamically independent motions in a complex system. Analysis of the modes provides meaningful poolings of metabolic concentrations (Palsson et al., 1987). Modal analysis results can be directly compared with the results from the phase plane analysis and the statistical analysis.

The modal matrix was used to identify metabolite pools at various timescales by the relative magnitude of contribution of each concentration metabolite to each mode. A pool was defined by lumping together those species whose ratio of contribution changed by less than 11% from one time step to the next. The 11% cutoff was identified by selecting a natural break in a histogram of the percent changes relative to each of several key variables. Pools were identified by normalizing all columns in the modal matrix by one column of interest. For example, we observed the G6P/F6P pool from our analysis. The modal matrix indicates that G6P and F6P contribute equally to each mode with a time-constant slower than ~10⁴ h.

On timescales slower than 10² h, FDP pools with GAP and DHAP on a 2:1:1 ratio, reflecting the stoichiometry of the aldolase reaction. On timescales longer than 10 h, the modal matrix indicates that the FDP/GAP/DHAP pool moves with DPG13, PG2, PG3, PEP, G6P, and F6P. Further analysis of the modal matrix also indicates that the PGs, DPG13, DHAP, GAP, FDP, GO6P, GL6P, RU5P, R5P, X5P, R1P, and INO pool together with G6P and F6P on timescales slower than 10¹ h.

The modal matrix identifies a key feature of DPG23 metabolism. DPG23 does not contribute significantly to modes faster than 10 h, suggesting that the DPG23 concentration only responds to perturbations on a 10-h timescale. The timescale response of DPG23 is indicative of its physiological function in the hRBC. On time scales greater than 10 h, the adenosine metabolites ADO, AMP, ADP, and ATP also form a pool in the modal matrix.

In conclusion, the modal matrix can be analyzed to identify the relevant pools in the RBC metabolic network based on an analytical analysis of the dynamic equations. Many of the predictions from the modal matrix are consistent with the phase plane and statistical approaches. However, there are some differences, such as the prediction that GO6P and GL6P move with RU5P, R5P, X5P, and R1P on slower time scales. The differences are due to nonlinear interactions in the model, which are not accounted for in the linearized modal analysis. However, under many circumstances, nonlinear interactions are important, as this analysis clearly demonstrates is the case for the RBC.

	FORMULATION OF PHYSIOLOGICALLY RELEVANT POOLS

TOP ABSTRACT INTRODUCTION MODEL DESCRIPTION AND SOLUTION... IDENTIFYING RELATIONSHIPS WITH... ANALYTICALLY IDENTIFYING THE... FORMULATION OF PHYSIOLOGICALLY... DYNAMIC RESPONSE OF RBC... DISCUSSION REFERENCES

Above, we have defined metabolite pools based on dynamic data or from dynamic analysis of a linearized model; however, one can also define physiologically meaningful metabolite pools a priori. The first choice of a metabolite pool of biological importance can be based on the pioneering ideas of Atkinson (1977), the energy charge:

(8)

As thoroughly discussed by Reich and Sel'kov (1981) the numerator is the "occupancy" of the adenosines moieties with high-energy phosphate bonds, whereas the denominator is a measure of the "capacity" to carry such high-energy bonds. The ratio is thus a "charge" measure. A ratio of 3 indicates that each adenosine has three high-energy phosphate bonds available for reaction, whereas a ratio of 0 indicates that no high-energy phosphate bonds are available.

It has been recognized that more general definitions of capacity and occupancy pools are consistent with the dynamic structure of RBC metabolism (Palsson et al., 1987). The entire RBC metabolic network can be viewed as an interaction of three separate charge pools: glycolytic, adenosine, and redox. Thus, rather than grouping the metabolites solely based on dynamic considerations (as discussed above), one can pool them according to the metabolic "value." This pooling is physiologically apt because the formation and hydrolysis of high-energy phosphate bonds and redox potential drives the metabolic "economy" of the cell.

The construction of the occupancy and capacity pools was based on the stoichiometry of the metabolic network (Table 3). The glycolytic high-energy phosphate occupancy pool was defined by the number of phosphates per six-carbons; for example, the occupancy of G6P is 1, and the occupancy of FDP is 2. The composition of all the pools is shown in Table 3. Note that we have not included the phosphates on the PPP intermediates. This is because the concentrations of these intermediates are much lower and will not significantly contribute to the high-energy phosphate occupancy. The redox occupancy was defined by the sum of the effects of the reduced electron carrier concentrations NADPH and NADH (Table 3).

Conceptually considering the red cell metabolic network purely in terms of energy and redox pools: one arrives at the reduced schematic of red cell metabolism shown in Fig. 8, which depicts the major functions of red cell metabolism. The glycolytic phosphate potential group interacts with all of the other groups. The NAD/NADH pool and the DPG23 pool will affect the Hb/met-Hb ratio and the affinity for Hb to bind oxygen. The redox pool combats oxidative stresses and the high-energy phosphate potential groups indicate how many phosphate equivalents the metabolites can transfer. It should be noted that the glycolytic capacity pool and the glycolytic phosphate occupancy pool serve as the backbone for the RBC metabolic network.

View larger version (18K): [in this window] [in a new window]   FIGURE 8   Simplified representation of the hRBC metabolic network using energy and redox pools.

	DYNAMIC RESPONSE OF RBC METABOLIC CURRENCY POOLS

TOP ABSTRACT INTRODUCTION MODEL DESCRIPTION AND SOLUTION... IDENTIFYING RELATIONSHIPS WITH... ANALYTICALLY IDENTIFYING THE... FORMULATION OF PHYSIOLOGICALLY... DYNAMIC RESPONSE OF RBC... DISCUSSION REFERENCES

The dynamic response of the metabolic network to perturbations can be analyzed in terms of the capacity, occupancy, and charge pools defined above. The response of the capacity, occupancy, and charge of each pool to two separate loads (ATP and redox) was studied. Due to the integrated nature of the metabolic network and the ability to convert one occupancy to another, the response of red cell metabolism to the two loads is similar. To demonstrate the types of behavior observed, Figs. 5 and 9 show the response of the red cell to ATP and redox loads. Fig. 5 shows the response to a low ATP or redox load as a phase plane for intermediate and long time intervals for each of the dynamic, data derived metabolite pools (Fig. 5 A) and the physiologically defined metabolite pools (Fig. 5 B). Fig. 9 shows the time domain responses of the red cell for the adenosine and redox pools to loads of varying magnitude. In all cases the initial adenosine and glycolytic transients are followed by a longer-term response that serves to stabilize the charge on each type of metabolic currency. Fig. 9 also shows that while capacity and occupancy vary by 10% to 15%, the metabolic charge varies less, changing by no more than 3% for reasonable loads. For the worst-case disturbances where charge does dramatically change, it should be noted that eventually, these simulations predict physiologically unobtainable states.

View larger version (22K): [in this window] [in a new window]   FIGURE 9   Time domain response to a step in ATP or redox loads. Time courses of the step responses of the capacity, occupancy, and charge (A) to high (0.4 mM/h), medium (0.3 mM/h), and low (0.2 mM/h) ATP loads and (B) to high (5 mM/h), medium (4 mM/h), and low (3 mM/h) redox loads.

Viewing the dynamic response of the red cell with these pooled variables several features become clear. 1) The occupancy and the capacity of the physiological pools move in a coordinated fashion, especially at the longer time intervals. The trace in the occupancy versus capacity phase plane forms a relatively straight line for all but the redox pools. 2) hRBC metabolism strives to keep the charge (occupancy/capacity) of its adenosine and glycolytic currencies a constant. Thus, the fluxes that move the occupancy and capacity pools respond to the imposed loads to maintain the physiological state of the cell. 3) The adenosine charge has robust, adaptive behavior. From Fig. 9, it is apparent that when a load (either of the two classes of loads shown) is imposed on the system, the adenosine charge is transiently affected, but on a longer timescale, the charge returns to near the original "set-point." Although not shown, similar behavior is observed for the glycolytic charge. 4) The occupancy and capacity pools can be interpreted from the modal matrix. If the physiological pools all moved independently, each pool would be equivalent to a row in the modal matrix. It is apparent that the pools are not strictly equivalent to the modes, but the motion of the pools can be interpreted from the modal matrix. From Fig. 9, it can be seen that the occupancy and capacity pools of the adenosine high-energy phosphate bonds moves on a different time scale than the redox occupancy/capacity and glycolytic occupancy/capacity pools, the timescale separation is indicated by the phase planes (as described in Fig. 2). From the modal matrix, it can be seen that the adenosine capacity pool is part of the long time constant modes (Table 2). The glycolytic high-energy phosphate occupancy and capacity pools are most similar to modes of intermediate time scale and the glycolytic pools should respond faster than the adenosine pools based on the modal matrix, and this is observed in Figs. 5 B and 9. 5) Although the hRBC metabolic model is robust to ATP and redox loads, there is a limit to that robustness. The metabolic network responds to nonspecific ATP and redox loads and is able to reject most of the disturbances with minor offset. 6) The phase plane shows that all metabolites respond to the imposed loads. However, the physiologically defined charges have a much more limited response to the loads than their individual components. Despite this, when the charge is perturbed beyond a critical value, the metabolic system "collapses" as demonstrated for the high ATP load and high redox loads in Fig. 9.

The relative motion of the occupancy and capacity of the physiologically relevant pools can be studied using the dynamic model of the hRBC. The time-parameterized phase planes clearly indicate the essential dynamic characteristics of the metabolic network and provide insight to the operation of the metabolic network.

	DISCUSSION

TOP ABSTRACT INTRODUCTION MODEL DESCRIPTION AND SOLUTION... IDENTIFYING RELATIONSHIPS WITH... ANALYTICALLY IDENTIFYING THE... FORMULATION OF PHYSIOLOGICALLY... DYNAMIC RESPONSE OF RBC... DISCUSSION REFERENCES

High-throughput experimental technologies and bioinformatics are placing an increased emphasis on mathematical modeling and computer simulation of cellular processes. The metabolic network in the hRBC will provide a useful "model" system for computational studies. Herein, we used a mathematical model of hRBC metabolism to study the use of computer simulation, modal analysis, phase plane analysis, and statistical analysis with the goal to unravel the complexities of this relatively simple metabolic network. The study showed 1) that these methods are well suited to study the relationship between all the detailed individual biochemical events and kinetic constants and the overall physiological function of red cell metabolism, 2) that network-wide definitions of metabolic pools and their charges with various metabolic currencies result in a few integrated measures that can be used to interpret overall behavior, 3) the metabolic charges are kept within a narrow range by the regulatory structure of the cell, indicating that the network strives to maintain homeostasis as well as it can given the constraints placed on it, and 4) the integrated function of a metabolic network can be described by a reduced set of variables.

Integrated models of red cell metabolism have been constructed for some 30 years (Ataullakhanov et al., 1981; Brumen and Heinrich, 1984; Heinrich and Rapoport, 1975; Jacobasch et al., 1983a,b; Lew and Bookchin, 1986; McIntyre et al., 1989; Rapoport and Heinrich, 1975; Schuster et al., 1988; Thorburn and Kuchel, 1987; Yoshida and Dembo, 1990). The result is a comprehensive, but not quite fully complete, model of all the major metabolic events in this simple cell. A full model contains not only the kinetic equations that describe enzyme kinetics, but also various constraints such as osmotic pressure and electroneutrality (Joshi and Palsson, 1989b; Werner and Heinrich, 1985). Given the numerical challenges associated with the physicochemical constraints, many studies, such as the present one, are focused on the kinetics alone.

Using phase planes, statistical time series, and modal analysis we can readily interpret fast motions in the network and pool metabolite concentrations into aggregate dynamic variables. This process is aided by representing the data as multiple matrices of pair-wise phase plane traces of metabolites representing dynamics on different time scales. Then using generalizations of Atkinson's energy charge (Atkinson, 1977) as a guide, further pooling leads to the definition of aggregates of concentrations that represent the network wide inventory of key metabolic currencies, such as high-energy phosphate bonds and two different types of redox potentials. Using these systemic variables, we can think of red cell metabolism as a battery that has a charge and capacity for these different metabolic currencies (e.g., see Reich and Sel'kov, 1981).

Dynamic simulation of the responses of the metabolic charges to various external loads is illuminating. The simulations show that, due to the interconnectivities of these pools, the cell responds in an analogous fashion to redox and energy loads by shifting its resources from one form to the other. The simulations also show that the network wide charges change numerically within narrow ranges in response to these loads. Thus, the overall regulatory scheme is coordinated to achieve homeostasis in terms of these pooled variables, but not necessarily in terms of the individual concentrations. If the loads placed on the cell are too severe, it cannot achieve this homeostasis, and its entire metabolism collapses. Therefore, the in silico analysis leads to the definition of what may be thought of as "emergent properties" of the network.

These in silico results call for further work on the hRBC. First, an experimental program should be designed to address the interpretations and predictions made herein. Further, a similar in silico study should be performed with the full model that includes all the physicochemical constraints on the cell. Second, an in silico analysis of the regulation that achieves the global properties described herein is needed. Deciphering the relative roles of network structure, mass action kinetics, and deliberate metabolic regulation of individual enzymes is needed. Finally, the lessons learned from the hRBC may be generalized to metabolic networks in cells that have less kinetic detail available by using the type of aggregate variables defined herein and the determination of the kinetic constants needed to construct such a description.