## Abstract

The dynamic behavior of metabolic networks is governed by numerous regulatory mechanisms, such as reversible phosphorylation, binding of allosteric effectors or temporal gene expression, by which the activity of the participating enzymes can be adjusted to the functional requirements of the cell. For most of the cellular enzymes, such regulatory mechanisms are at best qualitatively known, whereas detailed enzyme‐kinetic models are lacking. To explore the possible dynamic behavior of metabolic networks in cases of lacking or incomplete enzyme‐kinetic information, we present a computational approach based on structural kinetic modeling. We derive statistical measures for the relative impact of enzyme‐kinetic parameters on dynamic properties (such as local stability) and apply our approach to the metabolism of human erythrocytes. Our findings show that allosteric enzyme regulation significantly enhances the stability of the network and extends its potential dynamic behavior. Moreover, our approach allows to differentiate quantitatively between metabolic states related to senescence and metabolic collapse of the human erythrocyte. We think that the proposed method represents an important intermediate step on the long way from topological network analysis to detailed kinetic modeling of complex metabolic networks.

## Visual Overview

### Synopsis

One of the most challenging goals of computational systems biology is the development of detailed kinetic models to simulate and predict the dynamic response of metabolic networks towards, for example, changes in kinetic parameters due to pharmacological interventions or variations of environmental conditions. However, for most metabolic networks and the participating cellular enzymes the governing regulatory mechanisms are at best qualitatively known, whereas explicit kinetic models are lacking. In this work we present a kinetic‐free computational approach based on structural kinetic modeling (Steuer *et al*, 2006) to explore the dynamic properties of metabolic networks. Our method is based on metabolic states that are characterized by a flux distribution and metabolic concentrations. The energy and redox metabolism of human erythrocytes, which have been subject to mathematical modeling for decades (Schuster and Holzhütter, 1995), are used to exemplify and validate our approach. The time‐dependent behavior of a metabolic network can be described by a set of differential equations, derived from the stoichiometry and enzyme‐kinetic reaction rates. At a steady state, the metabolite concentrations do not change over time and the dynamic properties and hence the stability of such a metabolic state are determined by the Jacobian matrix of the set of differential equations. Reformulation leads to a representation of the Jacobian as a product of two matrices, the first one fully specified by the metabolic state and the second one encoding the relative saturation of each reaction with respect to its substrates. The usually unknown entries of this second matrix are called saturation parameters and are drawn randomly from well defined intervals in an iterative process. Evaluating the Jacobian matrix repetitively allows to investigate the scope of dynamic behavior (Figure 1).

One aspect of our analysis is to elucidate the role of allosteric regulation with respect to stability without knowledge of the detailed reaction rates. A comparison of the metabolic network of human erythrocytes under suppressed and under allowed allosteric regulation reveals the stabilizing effects of the latter.

We propose three objective measures to quantify the impact of the individual saturation parameters on stability, based on correlation coefficients, mutual information and a KS test, respectively. All these measures are shown to be consistent with each other. Almost all high‐ranked saturation parameters are associated with the highly regulated and almost irreversible reactions hexokinase, PFK and pyruvate kinase. Furthermore, in case of suppressed allosteric regulation, mainly those parameters are ranked very high that affect reactions that are known to be allosterically regulated, even without incorporating such prior knowledge to our analysis.

We demonstrate the discriminatory power of our approach to detect changes in stability properties of metabolic states. Therefore, we compared the normal *in vivo* state with a second one, representing higher energy demand due to osmotic stress (Dariyerli *et al*, 2004) or mechanic stress exerted during passage of the cell through thin capillaries (Kodicek, 1986). These states cannot be discriminated based on stoichiometric considerations alone because both satisfy the flux balance equation.

We show that the subspace of the saturation parameter space, that guaranties stability, is larger under normal condition than under increased energy demand. This is in accordance with our previous finding that energy related reactions are highly influencing stability. Moreover, the normal state is more robust towards perturbations in the saturation parameters and hence in the enzyme‐kinetic parameters (Figure 7). The normal state is clearly superior to the state of increased energy demand in its ability to sustain and to compensate such perturbations without loss of stability, which is supported by the fact that high energy demand will lead to a metabolic collapse of the red blood cell.

Methods that seek to bridge between purely stoichiometry‐based approached and explicit kinetic models are becoming increasingly important to infer dynamic properties of cellular pathways. The weak preconditions and the semi‐automatic and straightforward manner of its implementation make our approach a suitable starting point to elucidate and detect changes in dynamic properties of metabolic networks for which the construction of detailed kinetic models is not yet feasible.

Independent of any explicit kinetic models we explore the dynamic behavior of human erythrocytes is a systematic way.

Using different objective measures, the influence of the individual network parameters on stability, as one of the most important dynamic property, is quantified.

Even without incorporating knowledge about allosteric regulation and irreversibility, our analysis of the network parameters reveals the highly regulated and almost irreversible reactions hexokinase, phosphofructokinase and pyruvate kinase to have the highest influence on stability.

Two different metabolic states (normal in vivo against increased energy demand), which cannot be discriminated on stoichiometric considerations alone, are clearly distinguishable by their changing stability properties.

## Introduction

One of the most challenging goals of computational systems biology is the development of detailed kinetic models to simulate and predict the dynamic response of metabolic networks toward, for example, changes in kinetic parameters due to pharmacological interventions or variations of environmental conditions. However, for complex metabolic networks comprised of several interwoven pathways, detailed kinetic modeling is usually not possible due to the inevitable lack of knowledge about the kinetic properties of the involved enzymes and membrane transporters. In this work, we extend a recently proposed method that bridges between topology‐based approaches and explicit kinetic models of metabolic networks (Steuer *et al*, 2006). In the face of lacking or incomplete enzyme‐kinetic information, we (i) derive and compare statistical measures for the relative impact of enzymatic reactions and parameters on the dynamic properties (such as local stability) of metabolic networks, (ii) evaluate the functional role of allosteric feedback regulation in the stabilization of metabolic networks and (iii) propose measures to quantitatively evaluate the stability and robustness properties of metabolic states.

Our approach is exemplified and validated using a representation of the metabolic network of the human erythrocyte. Due to the fundamental role of erythrocytes in the oxygen supply of cells, as well as the relative simplicity of its metabolism, erythrocytes have been subject to extensive experimental and theoretical research for decades. Numerous explicit mathematical models have been developed since the late 1970s (Rapoport *et al*, 1976; Ataullakhanov *et al*, 1981; Holzhütter *et al*, 1985; Joshi and Palsson, 1989; McIntyre *et al*, 1989; Ni and Savageau, 1996a; Mulquiney and Kuchel, 1999; Nakayama *et al*, 2005) providing a suitable benchmark to assess the reliability of our method.

Our approach is motivated by the increasing experimental accessibility of cellular characteristics, such as metabolic fluxes and concentrations of metabolic intermediates (Fernie *et al*, 2004; Goodacre *et al*, 2004; Sauer, 2004). Each metabolic state, characterized by a flux distribution and metabolite concentrations, is associated with a unique spectrum of dynamic properties, as defined by the ensemble of all possible kinetic models consistent with the respective state. Our main focus thus lies on a quantitative characterization and comparison of the stability properties of metabolic states.

In particular, transitions to instability, occurring via a loss of a stable steady state, were previously argued to play a crucial role in senescence and metabolic collapse of erythrocytes and may act as a primary signal for cell removal in patients with hemolytic anemia (Schuster and Holzhütter, 1995; de Atauri *et al*, 2006).While usually an investigation of such transitions necessitates the construction of explicit kinetic models, our approach allows to draw quantitative conclusions about the stability of metabolic states in response to an increased ATP demand, occurring, for example, under conditions of osmotic or mechanic stress (Kodicek, 1986; Dariyerli *et al*, 2004). It is demonstrated that different metabolic states, each satisfying the flux balance equation and thermodynamic constraints, can nonetheless show drastic differences in the ability to ensure stability and maintain metabolic homeostasis.

As our method requires no detailed information about enzyme‐kinetic rate equations and parameters and due to its computational simplicity it is applicable to large metabolic networks. In particular, as the construction of explicit kinetic model is usually not feasible, our method significantly extends previous approaches to metabolic robustness, often based on topological or stoichiometric considerations alone (Edwards and Palsson, 2000; Stelling *et al*, 2002; Deutscher *et al*, 2006; Tekir *et al*, 2006). We argue that dynamic aspects of metabolic networks are becoming more and more important in view of modern techniques like siRNA knockdowns or genetic modifications (Bailey, 1991; Becker *et al*, 2005) to modify the activity of individual enzymes *in vivo*. It has to be expected that such perturbations may give rise to fundamental changes in the dynamic behavior of the underlying network.

## Results

### The parameterization of metabolic states

A metabolic network is a set of coupled chemical reactions and transport processes. Neglecting spatial variations of the metabolite concentrations within the reactions compartments the time‐dependent changes of the metabolite concentrations can be described by a set of differential equations of the form =** Nν**(

**), where**

*S***denotes the**

*S**m*‐dimensional vector of metabolite concentrations,

**denotes the**

*N**m × r*‐dimensional stoichiometric matrix and

**ν**(

**) denotes an**

*S**r*‐dimensional vector of enzyme kinetic reaction rates. In the case of lacking or incomplete enzyme‐kinetic data and assuming the existence of a stationary state

*S*^{0}, the differential equation can be interpreted as linear equation for the stationary reaction rates

**ν**

^{0}=

**ν**(

*S*^{0}). The mass balance equation

*N*ν^{0}=0 provides the conceptual foundation for current stoichiometry‐based approaches to metabolic network analysis and brought forth a number of highly successful applications to determine the structure and function of metabolic networks (Varma and Palsson, 1994; Schuster

*et al*, 1999; Stelling

*et al*, 2002). Recently, the flux balance equation was supplemented with thermodynamic constraints, providing a link between feasible flux distributions and metabolite concentrations (Holzhütter, 2004; Henry

*et al*, 2006; Kümmel

*et al*, 2006; Hoppe

*et al*, 2007).

However, the mass balance equation itself, along with its associated thermodynamic constraints, does not allow to draw any conclusions about the possible dynamics or potential instabilities of a metabolic state. To obtain information about essential aspects of the dynamics, we thus augment the mass balance equation with the first‐order expansion of the differential equation. Given a metabolic system at a (possibly unknown and not necessarily unique) metabolic state, characterized by **ν**^{0} and *S*^{0}, the system of differential equations can be approximated by a Taylor series expansion.

The first term describes the steady‐state properties of the system, as exploited by flux–balance analysis to constrain the stoichiometrically feasible flux distributions. Along similar lines, taking the next term of the expansion into account, the structure of the Jacobian matrix ** J** constrains the possible dynamics of the system at each metabolic state.

Our method builds upon a statistical evaluation of the Jacobian matrix. Based on the formalism of structural kinetic modeling (Steuer *et al*, 2006), we construct a parametric representation of the Jacobian matrix, such that each element covers the comprehensive parameter space at a specific metabolic state. In particular, the Jacobian matrix can be written as the product of two matrices **Λ** and **θ**_{x}^{μ}. The elements of the matrix **Λ** are fully specified by the metabolic state of the system. In addition, the (usually unknown) elements of the matrix **θ**_{x}^{μ} specify the relative saturation of each enzyme with respect to its ligands and can be assigned to well defined intervals even when the explicit functional form of the rate equations is not known. In the following, these matrix elements are denoted by saturation parameters θ_{m}^{r}, where *r* stands for the reaction saturated by metabolite *m*. A brief mathematical synopsis is given in Materials and methods.

Evaluating the Jacobian matrix with respect to the (unknown) elements of the matrix **θ**_{x}^{μ} then defines the spectrum or scope of dynamic behavior at the respective metabolic state. The proposed workflow is summarized in Figure 1. First, the stoichiometry and a metabolic state **ν**^{0} and *S*^{0} are specified, based on available experimental data and existing mathematical models. Second, an ensemble of Jacobian models is generated by assigning random values to the elements of **θ**_{x}^{μ}, obeying the defined intervals. Evaluating the eigenvalues of the Jacobian matrix repetitively, allows to investigate and compare the scope of dynamic behavior under different preconditions, for example, suppressed or absent allosteric regulation. In this respect, especially the largest real part of the eigenvalues, denoted by λ_{Re}^{max}, is of interest, as it relates to the slowest timescale of the system and, if positive, implies (local) instability of the metabolic state. The metabolic state is stable only if all eigenvalues have a negative real part (see also Materials and methods).

### The role of regulation

Allosteric regulation is one of the main mechanisms to control enzyme activity. Since allosteric regulation occurs within metabolic networks as feedback or feed‐forward loops, it operates network‐wide and affects the dynamic properties at a systems level. The presented approach is used to analyze the effects of allosteric regulation on stability in a systematic way. Two sets of models under different preconditions are created, both corresponding to the normal *in vivo* conditions of the erythrocyte (see Figure 1 for a schematic representation of the energy and redox metabolism of the human erythrocyte and the abbreviations used). Within the first set of models *C*_{noreg} all saturation parameters associated with allosteric effectors are fixed to zero, corresponding to absence of regulatory interactions. In addition, a second set *C*_{reg} is constructed by assigning all saturation parameters, including those for allosteric regulation, to their respective intervals (see Materials and methods and the Supplementary information for details).

Each Jacobian model is evaluated according to its spectrum of eigenvalues, with λ_{Re}^{max}>0 implying instability of the metabolic state. In the case of absence of allosteric regulation, corresponding to the set *C*_{noreg}, the proportion of dynamically stable models is approximately 81% (see Figure 2 for the estimated probability density functions of the largest real part λ_{Re}^{max}). Thus, although vast majority of the models are stable, the proportion of unstable models cannot be neglected. As dynamic stability is mandatory for the existence of the metabolic state, it indicates a substantial risk for the unregulated network to be driven out of the observed steady state when changes of the binding constants for the substrates occur for genetic or pharmacological reasons.

Within the whole set, no (Jacobian) model is found exhibiting more than one eigenvalue larger than zero, suggesting that the occurrence of a Hopf bifurcation is, at least, rare under the precondition of suppressed allosteric regulation. Since a Hopf bifurcation indicates the transition to sustained oscillation, such dynamical behavior seems unlikely under these conditions. Looking at the set *C*_{reg}, thus including allosteric regulation, the proportion of stable models shifts to 91%, which is significantly higher than in case of suppressed regulation (see Materials and methods). Similar observations were made by Ni and Savageau (1996b), where additional regulation was introduced to a model to stabilize the steady state. Note that here regulation is only defined qualitatively, that is, the actual strength of each regulatory interaction is chosen randomly and varies between the individual samples. Nonetheless, even without specific fine‐tuning of parameters, allosteric regulation results in a higher proportion of stable networks. This is presumably evolutionary advantageous, since a larger parameter subspace corresponding to stable models increases the flexibility to optimize parameter towards additional requirements other than stability.

Within the set of unstable models, 592 out of 10^{6} samples in *C*_{reg} have two eigenvalues greater than zero, in each case exhibiting complex conjugate imaginary parts. A smaller fraction (41 samples) shows three eigenvalues larger than zero, pointing to bifurcations of higher co‐dimension. Although restricted to a very small region in parameter space, allosteric regulation thus expands the scope of dynamical behavior by increasing the region in parameter space where oscillatory or more complex dynamics can be expected.

Furthermore, we tested if the observed increase is specific for the actual set of regulation parameters or if it can be achieved by any randomly chosen set of allosteric regulation parameters. To this end, instead of the actual regulation parameters, 10 putative allosteric regulations were selected randomly and the increase in the percentage of stable models was recorded. Most random sets (∼85%) of regulation parameters lead to a decrease in the percentage of stable models, demonstrating that not every possible arbitrary allosteric regulation has a positive effect on stability.

### The ranking of parameters

Stability of a metabolic steady state is an emergent systemic property that is brought about by the kinetic properties of all enzymes. Nevertheless, changes in the kinetic parameters of individual enzymes may have quite differential impact on the stability of a given steady state. To evaluate and compare the degree of influence of individual parameters and reactions on the stability and response to perturbations, we rank the parameters according to several objective measures. Three distinct measures are used and compared with each other, namely the (Pearson) correlation coefficient, the mutual information and the Kolmogorov–Smirnov test (KS test) (see Materials and methods for explicit definitions).

All three measures were evaluated for all saturation parameters for both sets of models *C*_{noreg} and *C*_{reg}. Although the detailed ranking of the parameters is not identical, it is still consistent with respect to all different measures. Figure 3 depicts the significant parameters for both sets of models *C*_{noreg} and *C*_{reg}. Figure 4 exemplifies the influence of the most highly ranked parameters on the stability of the metabolic state. The percentage of stable models within the parameter space as a function of selected saturation parameters is shown. For a more detailed discussion and comparison of the different rankings see also the Supplementary information.

The ranking of parameters allows for several significant conclusions about the role of regulation within the metabolic network. First, almost all high‐ranked parameters are associated with reactions involved in ATP production or consumption. Especially the PFK, HK and PK play an important role in stabilizing the network. Interestingly, evaluating the set *C*_{noreg} reveals that, although no additional information about putative sites for allosteric regulation is included, mainly those parameters are ranked very high that affect reactions that are known to be allosterically regulated (see Materials and methods for a statistical verification of this assertion). We point out that for the construction of *C*_{noreg}, only the stoichiometry and the metabolic state under normal conditions were used. This emphasizes the usefulness of our approach to analyze the metabolic networks that are not as well studied as the one of erythrocytes, and where detailed information about allosteric regulation is not available.

Several more observations can be made from the ranking of the parameters. First, the high‐ranked parameters almost all belong to enzymes of the glycolytic pathway. Intriguingly, kinetic alterations of the allosterically regulated enzyme G6PD, which is known to control the flux through the pentose phosphate pathway, does not show significant impact on stability. This corresponds with results from de Atauri *et al* (2006), obtained from an explicit kinetic model. The authors show that the metabolic network breaks down for low concentrations of the glycolytic enzymes, whereas such a transition to instability does not occur if the enzymes of the pentose phosphate pathway are at very low concentrations. Second, all parameters associated with 23P2G as an allosteric effector are relatively low ranked. This indicates that the main function of these regulatory mechanisms is not to maintain or achieve stability.

### Comparison with the explicit model

The availability of a comprehensive and well‐established mathematical model of the erythrocyte metabolism (Schuster and Holzhütter, 1995) allows to validate our method by comparing the ranking of saturation parameters with results of metabolic control analysis (MCA; see Heinrich and Schuster, 1996). As metabolic instability of the erythrocyte may occur if the energy demand exceeds the glycolytic ATP production, we study the impact that changes in the kinetic parameters of the various enzymes of the network have on ATP utilization. The relative change of the rate of ATP utilization (ν_{ATPase}) elicited by a (small) change of the Michaelis constant characterizing affinity of metabolite M to enzyme E is given by the flux control coefficient.

Note that negative values of the flux coefficient indicate that decreasing value of the Michaelis constant (corresponding to increasing saturation) increases the rate of ATP production and thus stabilizes the steady state. Calculating the flux control coefficient for all 85 Michaelis constants occurring in the rate laws of the mathematical model and ranking them in ascending order reveals that only 10 affinity parameters each contribute significantly to the energetic stabilization and destabilization of the network (Figure 5). Changes in the Michaelis constants for binding of AMP and ATP to the phosphopfructokinase (PFK) have by far the highest impact on the ATP production. This underlines the well‐known central regulatory importance of this enzyme for red cell glycolysis. Remarkably, the set of 20 regulatory most relevant Michaelis constants determined by MCA of the basis of the full mathematical model comprises all saturation parameters identified by our random sampling method.

### Robustness of metabolic states

As yet our analysis has focused on the analysis of a single metabolic state corresponding to the normal *in vivo* conditions of the erythrocyte. However, the energy metabolism of this cell has to cope with large fluctuations of the ATP demand as the activity of the Na/K‐ATPase, accounting for about 70% of the total ATP utilization, is greatly enhanced under conditions of osmotic stress (Dariyerli *et al*, 2004) or mechanic stress exerted during passage of the cell through thin capillaries (Kodicek, 1986). Moreover, because of lacking *de novo* protein synthesis the erythrocyte is extremely susceptible to enzyme deficiencies which typically result in an impairment of glycolytic ATP production and subsequent break down (hemolysis) of the cell (Jacobasch and Rapoport, 1996).

To demonstrate the discriminatory power of our approach to detect changes in the stability properties of metabolic states, we thus consider a second metabolic state of the erythrocyte, characterized by an increased energy demand. To this end, we use the kinetic model to calculate fluxes and metabolite concentrations at a sixfold higher energetic load as compared with the normal reference state. Switching from the normal *in vivo* state (*k*_{ATPase}=1.6 mM/h) to the new steady state at increased energetic load (*k*_{ATPase}=10 mM/h), the glycolytic flux increases from 1.5 to 2.33 mM/h (=155%); whereas the ATP concentration decreases from 1.6 to 0.56 mm(=35%). These relative changes are in excellent agreement with experimental data (160% increase of glycolytic flux at 35% decrease of ATP) obtained by successively decoupling glycolysis from ATP consumption by means of arsenate titrations (Ataullakhanov *et al*, 1981).

The parameterization of the second state *S*_{2}^{0} is performed as described before, with all saturation parameters, including allosteric regulation, sampled randomly from their respective intervals. Figure 6 depicts the resulting distribution of the largest real eigenvalue λ_{Re}^{max} within the spectra of eigenvalues, as compared to the distribution under normal conditions (state *S*_{1}^{0}). While for normal *in vivo* conditions, the proportion of stable Jacobian models was approximately 91%; this value drops drastically to only about 13.2% for the second state *S*_{2}^{0}. This is in accordance with our earlier observation that energy related reactions are most crucial with respect to stability, as well as with the fact that in the detailed kinetic model of Schuster and Holzhütter (1995) the system is able to compensate an increased energy load only up to an upper critical value, but breaks down if the energy demand is increased any further. We emphasize that both states cannot be discriminated based on stoichiometric considerations alone: Both satisfy the flux balance equation and are consistent with thermodynamic constraints.

The proportion of unstable models alone, evaluated over the comprehensive parameter space of a metabolic state, does not necessarily imply actual instability of the respective flux distribution. However, the proportion of unstable models has significant consequences for the ability of the system to maintain the considered metabolic state at perturbations of enzyme‐kinetic parameters (Morohashi *et al*, 2002). To evaluate the robustness properties of both states quantitatively, we consider two distinct scenarios: First, for both metabolic states random instances of stable models are repeatedly selected from the parameter space and the set of parameters is subsequently perturbed within a given radius. The percentage of perturbations that remain stable, as a function of the magnitude or radius of the perturbations, then serves as a quantitative measure of robustness. See Figure 7A for a schematic representation. Second, to make the results for both metabolic states more comparable, the parameter space from which random instances of models are selected is restricted to a small interval with λ_{Re}^{max}∈[0, −0.01] for both states. Again each parameter set is perturbed with increasing radius and the frequency with which a given magnitude of perturbations leads to instability is recorded. The results are shown in Figure 7.

Clearly, for both metabolic states the probability that a perturbation results in a loss of stability increases with increasing magnitude of perturbations. However, starting (by construction) with initially 100% of stable models, the fraction of models that become unstable increases significantly faster for the state *S*_{2}^{0}. This effect is even more pronounced in the second scenario. Here, the stable models for both metabolic states are initially restricted to similar real parts within the spectra of the eigenvalues, and thus to similar distances to the bifurcation. Nonetheless, a perturbation of the (initially stable) metabolic state *S*_{2}^{0} is much more likely to result in a transition to instability than corresponding perturbations of the normal state *S*_{1}^{0}. In this sense, the *in vivo* metabolic state *S*_{1}^{0}, and concomitantly also its observed that flux distribution is more robust than the second state *S*_{2}^{0}. We emphasize again that our quantification of robustness does not involve any knowledge about the explicit functional form of the rate equations or kinetic parameters. Nonetheless, different metabolic states can be clearly differentiated, based only on information about metabolite concentrations and associated flux patterns. In this respect, our approach gives valuable insights on the qualitative and quantitative dynamic behavior of metabolic states that cannot be obtained by considering the stoichiometric balance equation alone and is also applicable to situations where detailed knowledge about the explicit rate equations is not yet available.

## Discussion and conclusions

A central goal of metabolic regulation is homeostasis, that is the maintenance of a stable quasi‐stationary state under largely varying external conditions. In this work, we present a kinetic‐free approach that enables the identification of those enzymes and putative allosteric regulators having the largest impact on the stability of experimentally observed metabolic steady states. Our analysis was focused on three different aspects: first, the role of allosteric regulation was elucidated by comparing the dynamic behavior of the network under suppressed and allowed allosteric regulation. The proportion of stable models is significantly increased by allosteric regulation, showing that feedback regulation has a stabilizing effect. Second, three statistical measures were introduced to quantify the influence of enzymatic reactions and saturation parameters on stability in a systematic way. The parameters were ranked according to these measures. Intriguingly, almost all high‐ranked parameters are involved in one of the three reactions HK, PFK or PK, corresponding to those reactions that are indeed highly regulated and almost irreversible. We note that these results also hold when knowledge about allosteric regulation and irreversibility is not presupposed in the initial analysis. Third, we provided a quantitative measure to analyze different metabolic states with respect to their robustness towards perturbations in parameters. We compared the *in vivo* state with a second metabolic state, corresponding to an increased energy demand of the cell. With respect to robustness, the *in vivo* state is clearly superior in accordance with the fact that high energy demand will lead to a metabolic collapse of the red blood cells.

Our approach is essentially based on the knowledge of stoichiometry and metabolic state of the system. While detailed kinetic models are available for only very few metabolic networks, knowledge of metabolite concentrations and flux distributions becomes increasingly experimentally accessible (Fernie *et al*, 2004; Goodacre *et al*, 2004; Sauer, 2004). In this respect, we have demonstrated that different metabolic states, only characterized by a flux distribution and metabolite concentrations, are indeed associated with a unique spectrum of dynamic capabilities and can be differentiated based on their stability properties. As our method specifically samples the parameter space associated with a given metabolic state, it can be directly related to experimental observations and thus improves methods based on a straightforward sampling of kinetic parameters within an explicit kinetic model (von Dassow *et al*, 2000).

In particular, we expect that recent efforts for biotechnological modifications of metabolic systems will concomitantly result in fundamental changes in the dynamic behavior of these networks. While a desired flux distribution might be stoichiometrically feasible, unanticipated changes in dynamic properties can lead to a failure of network function. The weak preconditions and the semi‐automatic and straightforward manner of its implementation thus make our approach a suitable starting point to elucidate and detect changes in dynamic properties of metabolic networks for which the construction of detailed kinetic models is not yet possible.

## Materials and methods

### Models of the human erythrocyte

To exemplify and validate our approach, we mainly draw upon a previously published model of Schuster and Holzhütter (1995), consisting of 30 metabolites and 31 reactions (see Figure 1 for a schematic representation). The model was slightly modified to account for free inorganic phosphate and additional transport reactions for the educt glucose, the intermediate phosphate and the end products phosphoribosyl pyrophosphate, pyruvate and lactate. Mg‐complexes were omitted. All reactions, except ATPase, GSHox and PRPPT, were treated as reversible. As ATPase and GSHox are merged overall reactions describing energy consumption and oxidative load, product inhibition for these reactions was not included, that is, ADP and GSSG have no influence on ATPase and GSHox, respectively.

The kinetic model was used to calculate the steady metabolic state of the human erythrocytes under normal *in vivo* conditions (state *S*_{1}^{0}). Metabolite concentrations and flux values are given in the Supplementary information. In addition to the normal *in vivo* state *S*_{1}^{0}, a second steady‐state *S*_{2}^{0} was calculated, corresponding to an increased energy demand of the cell. Analogous simulations were performed previously to explore senescence and metabolic collapse of erythrocytes (Schuster and Holzhütter, 1995; de Atauri *et al*, 2006); see also Tekir *et al* (2006) for an analysis of red blood cell enzymopathies based on stoichiometric analysis.

### Structural kinetic modeling

Our analysis is based on a decomposition of the Jacobian matrix of a metabolic system at a state *S*^{0} into a product of two matrices. Given a metabolic system consisting of *m* metabolites and *r* reactions, the set of differential equations ** S**=

**(**

*N*ν**), which describe the time‐dependent behavior of all metabolite concentrations**

*S**S*

_{i}(

*t*) can be rewritten as

where *S*_{i}^{0} and ν_{j}^{0}=ν_{j}(*S*^{0}) denote the metabolic state at which the system is to be evaluated. Using the definitions given in (3) and the variable transformation *x*_{i}(*t*)=*S*_{i}(*t*)/*S*_{i}^{0}, the Jacobian with respect to the normalized variables *x* is

The scaled Jacobian ** J_{x}** is related to the original Jacobian by a simple similarity transformation and it is fully specified by the parameter matrices

**Λ**and

**θ**

^{μ}

_{x}. The elements of

**Λ**describe the time scales of the system, as specified by the metabolic state

*S*^{0}and

**ν**(

*S*^{0}). The (usually unknown) elements of

**θ**

^{μ}

_{x}are defined as the normalized derivatives of the reaction rates and, analogous to the scaled elasticity coefficients of MCA, denote the effective kinetic order or normalized saturation of each reaction with respect to its substrates. Each element θ

_{S}

^{v}is constrained to the interval of [0,1] if the metabolite

*S*is a substrate and [0,−1] if

*S*is a product of the reaction ν(

*S*). Additional non‐zero terms arise from allosteric regulation. Allosteric regulation is included by assigning the corresponding parameter to intervals, such that θ

_{S}

^{ν}∈[0,−

*n*] for inhibition and θ

_{S}

^{ν}∈[0,

*n*] for activation of a reaction ν by

*S*, respectively. A detailed derivation is given elsewhere (Steuer

*et al*, 2006).

#### Statistical sampling of the parameter space.

The stability and possible dynamics of the metabolic network are evaluated at a given metabolic state, characterized by metabolite concentrations *S*^{0} and fluxes **ν**(*S*^{0}). The vector of reaction rates satisfies the steady‐state condition ** Nν**(

*S*^{0})=

**0**and is described by

*r*‐rank(

**)‐free parameters. The vector of metabolite concentrations**

*N*

*S*^{0}is restricted by thermodynamic constraints (Henry

*et al*, 2006; Kümmel

*et al*, 2006) and approximated by values adapted from Schuster and Holzhütter (1995). The metabolic state fully specifies the matrix

**Λ**. To evaluate the dynamic capabilities at a given metabolic state, the non‐zero elements of the matrix

**θ**

^{μ}

_{x}are sampled from their predefined intervals (Wang

*et al*, 2004; Steuer

*et al*, 2006), while the elements of the matrix

**Λ**are restricted to the respective metabolic state. The schematic workflow is shown in Figure 1.

Specifically, each reversible enzyme‐kinetic reaction ν_{i}(** S**) is split into a forward ν

_{i}

^{+}(

**) and backward rate ν**

*S*^{−}(

**) and described by the overall steady‐state flux ν**

*S*_{i}

^{0}, the flux ratio γ=ν

^{−}(

*S*^{0})/ν

^{+}(

*S*^{0}), the steady‐state

*S*^{0}, as well as by a set of saturation parameters

*θ*

_{si}

^{v}. Although the method does not presuppose a specific functional form of the rate equations, we illustrate the parameterization using a generic form of enzyme‐kinetic rate equations, such as

where *f* (*A*, *B*, *P*, *Q*) denotes a first order polynomial. The reaction is characterized by the net flux ν^{0}, the steady‐state concentrations *A*^{0}, *B*^{0}, *C*^{0}, *D*^{0}, as well as the flux ratio γ=ν^{−}/ν^{+}, relating to the (often accessible) equilibrium constant *K*_{eq}. The four unknown saturation coefficients apply to forward and backward rate separately, obeying the relationships θ_{A}^{ν+}∈[0, 1] and θ_{A}^{ν−}=θ_{A}^{ν+}−1∈[0,−1] for substrates and θ_{P}^{ν+}∈[0, −1] and θ_{P}^{ν−}=θ_{A}^{ν+}+1∈[0, 1] for products, respectively.

The model of the erythrocyte is parameterized by 87 saturation parameters for substrate and product dependencies of each reaction, as well as 10 additional parameters corresponding to allosteric regulation. See Supplementary information for a detailed listing. The parameters are denoted with reactions (superscript) and substrate (subscript) respectively, that is, θ_{Fru6P}^{PFK} denotes the dependence of the phosphofructokinase (PFK) on fructose‐6‐phosphate (Fru6P).

#### Stability and dynamics of metabolic states.

Our method is based upon a statistical evaluation of the Jacobian matrix. In particular, a metabolic state that satisfies the steady‐state condition ** Nν**(

*S*^{0})=0 must not necessarily be stable. Rather, its dynamic stability is determined by the eigenvalues of the Jacobian at the respective state. Each eigenvalue describes the behavior of the system after an (infinitesimal) perturbation of the concentrations (Heinrich and Schuster, 1996). The possible dynamics in the vicinity of a metabolic state are schematized in Figure 8: The metabolic state can either be (i) a stable (attracting) steady state, characterized by a largest real part of the eigenvalues λ

^{max}

_{Re}<0 (ii) an unstable (repelling) state, characterized by a positive largest real part λ

^{max}

_{Re}>0 within the spectrum of eigenvalues, or (iii) a stable (attracting) focus, characterized by non‐zero (complex conjugate) imaginary parts λ

^{max}

_{Re}±λ

^{max}

_{Im}with λ

^{max}

_{Re}<0 and λ

^{max}

_{Im}≠0, (iv) an unstable focus, characterized by a positive real part λ

^{max}

_{Re}>0 and complex conjugate eigenvalues λ

^{max}

_{Im}≠0.

Of particular interest are also transitions between the scenarios (bifurcations), most importantly the Hopf bifurcation, where a pair of complex conjugate eigenvalues cross the imaginary axis (stable → unstable focus) and bifurcations of the saddle‐node type, where the largest real part within the eigenvalues crosses the imaginary axis (stable node → unstable saddle). Further types of bifurcations are discussed in Steuer *et al* (2006). We emphasize that stability does not imply constancy of a metabolic state. Rather, local dynamic stability is mandatory for the existence of the state, but all actual states will fluctuate around their average values (Steuer *et al*, 2003). All reported results are robust against small deviations of metabolite concentrations and flux values, that is, an analysis with small alterations of the metabolic state under normal conditions yields identical results.

#### A simple example.

To illustrate our approach, we briefly consider the simple example pathway depicted in Figure 9.

Within our approach, not assuming any further knowledge of the explicit rate equations and parameters, the system is parameterized by the matrices **Λ** and **θ**_{A}^{μ}

where **Λ** defined the metabolic state of the system, constraint by ν_{1}^{0}=ν_{2}^{0}+ν_{3}^{0}. For simplicity, we assume linear dependence of ν_{2}(*A*) on its substrate *A*, and thus θ_{A}^{μ2}=1. The parameter θ_{A}^{μ3}∈[1−*n*, 1] includes possible non‐linear inhibition of ν_{3}(*A*) by its substrate *A*. As the units of times and concentrations are arbitrary, we set ν_{1}^{0}=1 and *A*^{0}=1 without loss of generality. The Jacobian at an observed metabolic state (specified by the matrix **Λ**) is thus given as

Figure 9 shows the region of stability of the observed state ν_{3}^{0} versus the (unknown) parameter θ_{A}^{μ3}. More importantly, the observed metabolic state, here only characterized by ν_{3}^{0}, restricts the stability properties of the state. For small flux the observed state is always stable, that is, there exists no set of parameters such that the state is unstable. However for high flux ν_{3}^{0} the system might lose stability and is stable only in a small region of the (unknown) parameter space. This behavior is again exemplified in Figure 10, using explicit differential equations. Starting in the vicinity of the metabolic state {*A*^{0}, ν_{3}^{0}}, all other parameters are chosen randomly from the comprehensive parameter space. For small ν_{3}^{0} (left plot), the system will always decay back to the state. Hence, the steady state remains stable, independent of the actual value of θ^{3}_{A}. However, for large ν_{3}^{0} (right plot) the probability of instability increases. For some perturbations of θ_{A}^{3} the steady state becomes unstable, i.e. the initial state cannot be restored. The systems transit into a new stable steady state with concentrations of **A** different from the initial value **A**=1. In this sense, the observed metabolic state puts constraints on the possible dynamics and allows to quantify the existent size of unstable regions in parameter space. The perturbation analysis shows that the metabolic state with small ν_{3}^{0} is more robust against changes in parameters than a metabolic state with large ν_{3}^{0}.

### The role of regulation

The evaluation of models (Jacobians) for *C*_{reg} and *C*_{noreg} was repeated 10^{3} times and the percentage of unstable instances recorded. In both cases, the variance of the values, due to finite sampling effects, was estimated numerically. A *t*‐test was used to test the average values for both cases *C*_{reg} and *C*_{noreg} against each other. The null hypothesis that the expectation value in case of allowed allosteric regulation is equal or lower than in case of suppressed allosteric regulation is rejected with a *P*‐value below 10^{−320}. So allosteric regulation significantly increases the frequency of stable models.

### Ranking of parameters

To assess the relative impact of enzyme‐kinetic parameters on the stability properties of a given metabolic state, we employed and compared several measures of dependency.

The most common choice to detect dependencies between variables is the (Pearson) correlation coefficient *r*, defined as

where *x*_{i} and *y*_{i} are *n* realizations of the random variables *X* and *Y*. Although the Pearson correlation only detects linear dependencies, it holds the advantage that its sign specifies whether a parameter must be increased or decreased to obtain a higher percentage of stable models. Nonetheless, the Pearson correlation suffers from several drawbacks, such as sensitivity to non‐Gaussian and skewed distributions, making more elaborate measures necessary (Kumar and Shoukri, 2007).

A more general measure of dependency is given by the mutual information (Shannon, 1948), defined as

where *H*(*X*)=∑_{k}*p*_{k}(*x*) log*p _{k}*(

*x*) denotes the entropy of the variable

*X*, measured from a binned histogram such that each bin occurs with probability

*p*. The entropy

_{k}*H*(

*X*,

*Y*) denotes the joint entropy of

*X*and

*Y*. Among its main advantages is that the mutual information is zero if and only if both variables are statistically independent. Evaluating the mutual information between network parameters and the resulting largest real part of the eigenvalue λ

_{Re}

^{max}thus accounts for arbitrary nonlinear dependencies and does not presuppose Gaussian or uniform distribution of the parameters. A detailed account of its numerical estimation is given elsewhere (Steuer

*et al*, 2002).

Based on a slightly different concept, the KS test is used to test for the equality of two distributions. The null hypothesis in the context of our analysis is as follows: if a saturation parameter has no impact on the stability of the metabolic system, then its distribution within the restricted subset of stable models equals (in a statistical sense) its initial distribution for the comprehensive set of models. On the other hand, if the distribution of a parameter within the restricted set of stable models shows a strong deviation from the initial distribution, a significant dependency can be expected. In this sense, the KS test determines whether two random variables *X* and *Y* have the same distribution. The cumulative frequencies *F*_{X} and *F*_{Y} and the maximal difference

are calculated. If the test statistic *D* is greater than the critical value for the sample size, the null hypothesis that both distributions are equal is rejected. Since *D* is always identically distributed, the KS test is independent of the distribution of *X* and *Y*. If the test rejects the null hypothesis, given a sufficiently small *P*‐value, the parameter under consideration has significant impact on stability.

All measures yielded consistent results, as depicted in Figure 11.

#### Significance of ranking.

To test for the significance of the correlation coefficient and the mutual information, we employed a permutation test: the values of λ_{Re}^{max} were randomly permuted in order to abolish any relationship between the saturation parameters and λ_{Re}^{max}. This yielded mean values of the correlation coefficient and the mutual information close to zero (<5 × 10^{−5}) and standard deviation of 3.1 × 10^{−3} and 2.5 × 10^{−4}, respectively. The correlation coefficient and the mutual information for the high‐ranked parameters are indeed significantly larger than those obtained in case of totally unrelated saturation parameters.

To verify the assertion that there is an enrichment of actual feedback parameters in the top‐ranking parameters, we conducted two statistical tests (i) for the ensemble of models with absent regulation (*C*_{noreg}), we verify that high‐ranking parameters are primarily associated with reactions that are actually allosterically regulated (PK, PFK, PK, G6PD and 6PDG). To this end, we record the fraction of parameters associated with regulated reactions within the top *k*‐ranked parameters. This number is then compared (statistically) to the number that must be expected if the high‐ranking parameters are indeed randomly distributed across all reactions. The results are depicted in Figure 12A (as a function of *k*) and show a clear significant enrichment of parameters associated with regulated reactions among the top‐ranking parameters; (ii) for the ensemble of models including allosteric regulation (*C*_{reg}), we evaluate if regulation parameters are statistically overrepresented in the set of top‐ranking parameters. Again, we record the expected number of regulation parameters within the *k* top‐ranking parameters, based on a purely random distribution of parameters. This value is compared with the actual number of regulation parameters within the set of top‐ranked parameters. The result is again significant and depicted in Figure 12B.

## Implementation

All calculations were made using MATLAB, Version 7.3.0.298. The code is available upon request.

## Acknowledgements

The work of SB was funded by the systems biology initiative of the Federal Ministry of Education and Research, Grant Nr 0313078. JS is supported by the *Go*FORSYS project funded by the Federal Ministry of Education and Research, Grant number 0313924.

## Supplementary Information

Supplementary Information [msb4100186-sup-0001.pdf]

## References

This is an open access article under the terms of the Creative Commons Attribution‐NonCommercial‐NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non‐commercial and no modifications or adaptations are made.

- Copyright © 2007 EMBO and Nature Publishing Group