New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Structural kinetic modeling of metabolic networks

Edited by Philip P. Green, University of Washington School of Medicine, Seattle, WA, and approved June 6, 2006 (received for review January 8, 2006)
Abstract
To develop and investigate detailed mathematical models of metabolic processes is one of the primary challenges in systems biology. However, despite considerable advance in the topological analysis of metabolic networks, kinetic modeling is still often severely hampered by inadequate knowledge of the enzyme–kinetic rate laws and their associated parameter values. Here we propose a method that aims to give a quantitative account of the dynamical capabilities of a metabolic system, without requiring any explicit information about the functional form of the rate equations. Our approach is based on constructing a local linear model at each point in parameter space, such that each element of the model is either directly experimentally accessible or amenable to a straightforward biochemical interpretation. This ensemble of local linear models, encompassing all possible explicit kinetic models, then allows for a statistical exploration of the comprehensive parameter space. The method is exemplified on two paradigmatic metabolic systems: the glycolytic pathway of yeast and a realisticscale representation of the photosynthetic Calvin cycle.
Cellular metabolism constitutes a complex dynamical system and gives rise to a wide variety of dynamical phenomena, including multiple steady states and temporal oscillations. The elucidation, understanding, and eventually prediction of the behavior of metabolic systems represent one of the primary challenges in the postgenomic era (1–5). To this end, substantial effort has been dedicated in recent years to develop and investigate detailed models of cellular metabolic processes (6, 7).
Once a mathematical model is established, it can serve a multitude of purposes: It can be regarded as a “virtual laboratory” that allows the building up of a characteristic description of the system and gives insights into fundamental design principles of cellular functions, such as adaptability, robustness, and optimality (8–10). Likewise, mathematical models of cellular metabolism serve as a basis to investigate questions of major biotechnological importance, such as the effects of directed modifications of enzymatic activities to improve a desired property of the system (11).
However, although there has been a formidable progress in the structural (or topological) analysis of metabolic systems (12, 13), and despite the long history of metabolic modeling, dynamic models of cellular metabolism incorporating a realistic complexity are still scarce.
This scarcity is owed to the fact that the construction of such models encompasses a number of profound difficulties. Most importantly, the construction of kinetic models relies on the precise knowledge of the functional form of all involved enzymatic rate equations and their associated parameter values. Furthermore, even if both are available from the literature, parameter values may (and usually do) depend on many factors such as tissue type or experimental and physiological conditions. Likewise, most enzyme–kinetic rate laws have been determined in vitro, and often there is only little guidance available whether a particular rate function is still appropriate in vivo.
In this work, we aim to overcome some of these difficulties and propose a bridge between structural modeling, which is based on the stoichiometry alone (12–14), and explicit kinetic models of cellular metabolism. In particular, we demonstrate that it is possible to acquire an exact, detailed, and quantitative picture of the bifurcation structure of a given metabolic system, without explicitly referring to any particular set of differential equations.
Our approach starts with the observation that in most circumstances an explicit kinetic model is not necessary. For example, to determine under which conditions a steady state loses its stability, only a local linear model of the system at this state is needed, i.e., we only need to know the eigenvalues of the associated Jacobian matrix. Note that by stating this assertion, and unlike related approaches to qualitative modeling (14, 15), we do not aim at an approximation of the system. The boundaries of an oscillatory region in parameter space that arise out of a Hopf (HO) bifurcation are actually and exactly determined by the eigenvalues of the Jacobian. Likewise, other bifurcations, including bifurcations of higher codimension, can be deduced from the spectrum of eigenvalues and give rise to specific dynamical behavior.
The basis of our approach thus consists of giving a parametric representation of the Jacobian matrix of an arbitrary metabolic system at each possible point in parameter space, such that each element is accessible even without explicit knowledge of the functional form of the rate equations. Once this representation of the Jacobian is obtained, it allows to give a detailed statistical account of the dynamical capabilities of a metabolic system, including the stability of steady states, the possibility of sustained oscillations, as well as the existence of quasiperiodic and chaotic regimes. Moreover, the analysis is quantitative, i.e., it allows the deduction of specific biochemical conditions under which a certain dynamical behavior occurs and allows the assessment of the plausibility or robustness of experimentally observed behavior by relating it to a quantifiable region in parameter space.
Structural Kinetic Modeling
The temporal behavior of a metabolic network, consisting of m metabolites and r reactions, can be described by a set of differential equations (6) where S denotes the mdimensional vector of biochemical reactants and N the m × r stoichiometric matrix. The rdimensional vector of reaction rates ν(S, k) consists of nonlinear (and often unknown) functions, which depend on the substrate concentrations S, as well as on a set of parameters k.
In the following, we will not assume explicit knowledge of the functional form of the rate equations but instead aim at a parametric representation of the Jacobian of the system. As the only mathematical assumption about the system, we require the existence of a positive state S ^{0} that fulfills the steadystate condition Nν(S ^{0}, k) = 0. Note that the state S ^{0} is not required to be unique or stable.
Using the definitions (16, 17) with i = 1, …, m and j = 1, …, r and applying the variable substitution S_{i} = x _{i} S _{i} ^{0}, the system can be rewritten in terms of new variables x(t)
The corresponding Jacobian of the normalized system at the steady state x ^{0} = 1 is Because the new variables x are related to S by a simple multiplicative constant, J_{x} can be straightforwardly transformed back into the original Jacobian.
Any further evaluation of the Jacobian now rests on the interpretation of the terms in Eq. 4 . We begin with an analysis of the matrix Λ: Its elements Λ_{ij} have the units of an inverse time and consist of the elements of the stoichiometric matrix N, the vector of steadystate concentrations S ^{0}, and the steadystate fluxes ν(S ^{0}). Provided a metabolic system is designated for mathematical modeling, we can assume that there exists some knowledge about the relevant concentrations, i.e., for each metabolite, we can specify an interval S _{i} ^{−} ≤ S _{i} ^{0} ≤ S _{i} ^{+}, which defines a physiologically feasible range of the respective concentration. Furthermore, the steadystate fluxes ν(S ^{0}) are subject to the massbalance constraint Nν(S ^{0}) = 0, leaving only r − rank(N) independent reaction rates (6). Again, an interval ν_{i} ^{−} ≤ ν_{i} ^{0} ≤ ν_{i} ^{+} can be specified for all independent reaction rates, defining a physiologically admissible flux space.
In the following, we denote S ^{0} and ν(S ^{0}), usually corresponding to an experimentally observed state of the system, as the “operating point” at which the Jacobian is to be evaluated. This information, together with the stoichiometric matrix N, fully specifies the matrix Λ.
The interpretation of the matrix θ _{x} ^{μ} in Eq. 4 is slightly more subtle because it involves the derivatives of the unknown functions μ(x) with respect to the new normalized variables at the point x ^{0} = 1. Nevertheless, an interpretation of these parameters is possible and does not rely on the explicit knowledge of the detailed functional form of the rate equations: Each element θ_{xi} ^{μj} of the matrix θ _{x} ^{μ} measures the normalized degree of saturation of the reaction ν_{j} with respect to a substrate S _{i} at the operating point S^{0} . In particular, the dependence of almost all biochemical rate laws ν_{j}(S) on a biochemical reactant S _{i} can be written in the form ν_{j}(S, k) = k _{v} S _{i} ^{n}/f _{m}(S, k), where n denotes an integer exponent and f _{m}(S, k) a polynomial of order m in S _{i} with positive coefficients k. All other reactants have been absorbed into k (6). After applying the transformation of Eq. 2 , we obtain
with α ∈ [0, 1] denoting a free variable in the unit interval. The limiting cases are always lim_{Si0→0α} = 1 and lim_{Si0→∞α} = 0. To evaluate the matrix θ _{x} ^{μ}, we thus restrict each saturation parameter to a well defined interval, specified in the following way: As for most biochemical rate laws n = m = 1, the partial derivative usually takes a value between zero and unity, determining the degree of saturation of the respective reaction. In the case of cooperative behavior with exponents n = m ≥ 1, the normalized partial derivative lies in the interval [0, n] and, analogously, in the interval [0, −m] for inhibitory interaction with n = 0 and m ≥ 1. For examples and proof of Eq. 5 , see, respectively, Materials and Methods and the Supporting Appendix and Figs. 8–11, which are published as supporting information on the PNAS web site.
The matrices θ _{x} ^{μ} and Λ, as defined above, fully specify the Jacobian of the system. In the following, both quantities are treated as free parameters, defining the physiologically admissible “parameter space” of the system. Importantly, our representation of the Jacobian fulfills the following three essential conditions. (i) The reconstructed Jacobian represents the exact Jacobian at this point in parameter space. There is no approximation involved. (ii) Each term in the Jacobian is either directly experimentally accessible, such as flux or concentration values, or has a well defined biochemical interpretation, such as a normalized degree of saturation of a given reaction. (iii) The Jacobian does not depend on a particular choice of specific rate functions. Rather, it encompasses all possible kinetic models of the system that are consistent with the considerations above. In this sense, the reconstructed Jacobian is exhaustive.
An Illustrative Example
Before an application to more detailed biochemical models, we exemplify our approach using a simple hypothetical pathway. Suppose the reaction scheme depicted in Fig. 1, consisting of two metabolites and three reactions, is designated for mathematical modeling. The starting point of our analysis is then usually an experimentally observed operating point, characterized by metabolite concentrations S^{0} = (G ^{0}, T ^{0}) and flux values ν^{0} = (ν_{1} ^{0}, ν_{2} ^{0}, ν_{3} ^{0}). Furthermore, an analysis of the stoichiometric matrix N reveals that there is only one independent steadystate reaction rate c, with ν_{1} ^{0} = ν_{2} ^{0} = c, and ν_{3} ^{0} = 2c. Thus, we only require knowledge of the average overall flux through the pathway, specifying the value c. This information already enables the construction of the matrix Λ, which defines the operating point at which the system is to be evaluated.
The only remaining parameters are now the elements of the matrix θ _{x} ^{μ}. Starting with the dependence of each reaction on its substrate and assuming conventional biochemical rate laws, we obtain θ_{G} ^{μ2}∈[0, 1], specifying the degree of saturation of ν_{2} with respect to its substrate glucose (G). Furthermore, θ_{T} ^{μ3}∈[0, 1] specifies the degree of saturation of ν_{3} with respect to ATP (T). Additionally, the known regulatory feedback of the metabolite T upon the reaction ν_{2} is incorporated by θ_{T} ^{μ2}∈[0, n], where n ≥ 1 denotes a positive integer. The matrix θ _{x} ^{μ} thus contains three nonzero values, each restricted to a well defined interval
We emphasize that the three elements of θ _{x} ^{μ} represent bona fide parameters of the system, specifying the Jacobian matrix no less unique and quantitative than a corresponding set of Michaelis constants. Given the elements of θ _{x} ^{μ} as free parameters, we thus have obtained a parametric representation of the Jacobian matrix, which encompasses all possible kinetic models consistent with the experimentally observed operating point. In the remainder of this work, we use our approach to evaluate the dynamical capabilities of two more complex examples of metabolic system.
Glycolytic Pathway
Among the most classical and probably best studied examples of a biochemical oscillator is the breakdown of sugar by means of glycolysis in yeast. Damped and sustained glycolytic oscillations have been observed for several decades and have triggered the development of a large variety of kinetic models (18–20). In the following, we will address some of the characteristic questions that led to the development of those earlier models and show that these questions can be readily answered by using the concept of structural kinetic modeling. Given a schematic representation of the pathway, as depicted in Fig. 2, the first and foremost question is to establish whether the proposed reaction mechanism indeed facilitates sustained oscillations at the observed operating point. And, if yes, what are the specific kinetic conditions under which sustained oscillations can be expected?
We start out by constructing the matrix Λ using the experimentally observed state S ^{0} and ν ^{0}, identified here with the average concentration and flux values reported in refs. 19 and 20. Additionally, the matrix of saturation coefficients θ _{x} ^{μ} has to be specified. For simplicity, we assume that all reactions are irreversible and depend on their respective substrates only, resulting in 13 free parameters. Based on our discussion of conventional biochemical rate laws above, the saturation coefficients are restricted to the unit interval θ_{S} ^{μ} ∈ [0, 1].
For the dependence of the PFK–HK reaction on ATP, we follow a previously proposed kinetic model (19) and assume linear activation because of its effect as a substrate and a saturable inhibition involving a positive exponent n ≥ 1. The corresponding parameter is thus θ_{ATP} ^{μ1} = 1 − ξ, with ξ ∈ [0, n]. No further assumptions about the detailed functional form of any of the rate equations are necessary. For an explicit representation of both matrices Λ and θ _{x} ^{μ}, see the Supporting Appendix. To investigate the possibility of sustained oscillation, we begin with the most simple scenario and set θ_{S} ^{μ} = 1 for all reactions, corresponding to bilinear massaction kinetics. Note, however, that the inhibition term is still assumed to be an unspecified nonlinear function. Fig. 3 shows the largest eigenvalue of the resulting Jacobian at the experimentally observed operating point as a function of the feedback strength ξ. For sufficient inhibition, the spectrum of eigenvalues passes through a HO bifurcation, and the system facilitates sustained oscillations. Importantly, for a HO bifurcation to occur at the observed operating point, an exponent n ≥ 2 is needed, irrespective of the detailed functional form of the rate equation.
We have to highlight one fundamental aspect of our analysis: Given our parametric representation of the Jacobian, the impact of the inhibition is decoupled from the steadystate concentrations and flux values the system adopts (the latter being solely determined by the matrix Λ). Thus, we specifically ask whether the assumed inhibition is indeed a necessary condition for the observation of oscillations at the experimentally observed operating point. In contrast to this fact, using a conventional kinetic model and reducing the influence of the regulation, i.e., by increasing the corresponding Michaelis constant, would concomitantly result in altered steadystate concentrations, thus not straightforwardly contributing to the answer to this question.
Furthermore, because glycolytic oscillations have no obvious physiological role and are only observed under rather specific experimental conditions, some questions concerning their possible functional significance have been raised. One assertion is that the observed oscillations might only be an unavoidable side effect of the regulatory interactions, optimized for other purposes (6). Indeed, as shown in Fig. 3, a varying feedback strength ξ allows for different dynamical regimes. In particular, an intermediate value speeds up the response time with respect to perturbations, as also frequently observed in explicit models of cellular regulation (21).
Statistical Analysis of the Parameter Space
Going beyond the case of bilinear kinetics, we now evaluate the properties of a Jacobian at the most general level. All saturation coefficients θ_{S} ^{μ} ∈ (0, 1] are allowed to take arbitrary values in the unit interval, encompassing all possible explicit kinetic models of the pathway shown in Fig. 2. The steadystate concentrations and flux values are again restricted to the experimentally observed operating point. To assess the dynamical properties of the system, the saturation coefficients θ_{S} ^{μ} ∈ (0, 1] are repeatedly sampled from a uniform distribution. For each random realization the Jacobian is evaluated, and the largest real part λ_{R} ^{max} of its eigenvalues is recorded. Fig. 4 shows the histogram of the largest real part within the spectrum of eigenvalues, with λ_{R} ^{max} > 0 implying instability of the operating point. In the absence of inhibitory feedback ξ = 0, the operating point is likely to be unstable, i.e., most realizations result in a spectrum of eigenvalues with at least one positive real part.
Two ways to circumvent this inherent instability are conceivable. First, we can ask about the dependence on particular reactions, that is, whether the saturation (or nonsaturation) of a specific reaction contributes to an increased stability of the system. To this end, the correlation coefficient between λ_{R} ^{max}, reflecting the stability of the system, and the saturation parameters θ_{S} ^{μ} was estimated. Indeed, several parameters θ_{S} ^{μ} show a strong correlation with λ_{R} ^{max}, indicating that their value essentially determines the stability of the system (for data see Supporting Appendix and Figs. 12–14, which are published as supporting information on the PNAS web site). Fig. 4 Left depicts the distribution of λ_{R} ^{max} under the assumption that these reactions are restricted to weak saturation. In this case, the resulting distribution is shifted toward negative values, corresponding to an increased probability of the system to operate at a stable steady state.
The second option to ensure stability of the system arises from the negative feedback of ATP upon the combined PFK–HK reaction. Fig. 4 Right shows the distribution of the largest real part λ_{R} ^{max} of the eigenvalues for a nonzero feedback strength ξ > 0. Again, the distribution is markedly shifted toward negative values, increasing the probability of a stable steady state.
To investigate the role of the feedback in more detail, Fig. 5 depicts the distribution of λ_{R} ^{max} as a function of the feedback strength ξ. As can be observed, in the absence of the regulatory feedback the system is prone to instability, i.e., it is not possible (or rather unlikely) for the observed operating point to exist as a stable steady state. Subsequently, as the feedback strength is increased, the probability of obtaining a stable steady state increases. For an intermediate value ξ = 1, the system is fully stable: Any realization of the Jacobian will result in a stable steady state, independent of the detailed functional form of the rate equations or their associated parameters. However, as the feedback is increased further, the operating point again loses its stability. This time the instability arises out of a HO bifurcation, indicating the presence of sustained oscillations.
Based on these findings, we can summarize some essential properties of the pathway depicted in Fig. 2: Given the experimentally observed metabolite concentrations and flux values, our results show that in the absence of the regulatory interaction it would not be possible (or highly unlikely) to observe either sustained oscillations or a stable steady state. However, for sufficiently large inhibitory feedback, the system will inevitably exhibit sustained oscillations. Furthermore, as the feedback strength ξ ∈ [0, n] is bounded by an exponent n of the (unspecified) rate equation, n ≥ 2 is required for the existence of sustained oscillations. As demonstrated, our method thus allows one to derive the likeliness or plausibility of the experimentally observed oscillations, as well as the specific kinetic requirements for oscillations to occur, without referring to the detailed functional form of the rate equations.
Photosynthetic Calvin Cycle
The CO_{2}assimilating Calvin cycle, taking place in the chloroplast stroma of plants, is a primary source of carbon for all organisms and is of central importance for many biotechnological applications. However, even when restricting an analysis to the core pathway, the construction of a detailed kinetic model already entails considerable challenges with respect to the required rate equations and kinetic parameters (22, 23).
In the following, we thus use a representation of the photosynthetic Calvin cycle, as adapted from earlier kinetic models (22, 23), to demonstrate the applicability of our approach to a system of a reasonable complexity. The structural kinetic model consists of 18 metabolites, subject to 2 conservation relations, and 20 reactions, including 3 export reactions, starch synthesis, and regeneration of ATP. For a schematic representation of the pathway, see Supporting Appendix and Figs. 15–17, which are published as supporting information on the PNAS web site.
We seek to describe a general strategy to extract information about the dynamical capabilities of the system, without referring to an explicit set of differential equations. Our agenda focuses on (i) the stability and robustness of the experimentally observed concentration and flux values, (ii) the relative impact or importance of each reaction on the dynamical properties of the system, (iii) the existence and quantification of different dynamical regimes such as oscillations and multistability, and (iv) the possibility of complex or chaotic temporal behavior.
The starting point is again an experimentally observed state, characterized by the vector of metabolite concentrations S^{0} and flux values ν^{0} , as reported by Petterson and RydePetterson (22). Although additional knowledge on the reactions is often available, for the moment we assume that all reactions depend only on their substrates and products, with parameters θ_{S} ^{μ} ∈ (0, 1] and θ_{P} ^{μ} ∈ [0, −1], respectively. This information, embedded within the matrices Λ and θ _{S} ^{μ}, constitutes the structural kinetic model of the Calvin cycle at the observed operating point.
As a first approximation, we commence with global saturation parameters, θ_{S} ^{μ} and θ_{P} ^{μ}, set equal for all reactions. Although clearly oversimplified, the resulting bifurcation diagram, depicted in Fig. 6, already reveals some fundamental dynamical properties of the system. First, the observed operating point is indeed a stable steady state for most parameters θ_{S} ^{μ} and θ_{P} ^{μ}. Interestingly, however, in the absence of product inhibition θ_{P} ^{μ} = 0, a steady state is no longer feasible. In particular, for pure irreversible massaction kinetics (θ_{S} ^{μ} = 1, θ_{P} ^{μ} = 0), corresponding to a nonenzymatic chemical system, the pathway could not operate at the observed steady state. Second, for low product saturation (θ_{P} ^{μ} close to zero), a HO bifurcation occurs. Although this result does not necessarily imply that this region within parameter space is actually accessible under normal conditions, it shows the dynamical capability of the system to generate sustained oscillations; i.e., there exists a region in parameter space that allows for oscillatory behavior. Additionally, for low values of the substrate saturation θ_{S} ^{μ}, a saddlenode (SN) bifurcation occurs. This result shows that the observed steady state will eventually lose its stability, i.e., there are conditions under which the observed steady state is no longer stable. And, indeed, both dynamical features have been observed for the Calvin cycle: Photosynthetic oscillations are known for many decades and have been subject to extensive experimental and numerical studies (24). Likewise, multistability was recently found in a detailed kinetic model of the Calvin cycle and verified in vivo (23).
To proceed with a systematic analysis, the next step is to drop the assumption of global saturation parameters. All individual parameters θ_{S} ^{μ} ∈ (0, 1] are now allowed to take arbitrary values in the unit interval, reflecting the full spectrum of possible dynamical capabilities of the metabolic system. For simplicity, all reactions are still restricted to weak saturation by their products θ_{P} ^{μ} = −1/3. Of foremost interest is again the stability of the experimentally observed operating point: Evaluating an ensemble of 5 × 10^{5} random realizations of the Jacobian at this operating point, the system gives rise to a stable steady state in ≈94.3% of all cases (see Supporting Appendix and Figs. 15–17 for convergence and dependence on ensemble size). Thus, the stability of the observed operating point is indeed generic and does not rely on a specific choice of the kinetic parameters.
As for the remaining ≈5.7% of models, corresponding to the case where the observed operating point is unstable, ≈5.1% give rise to a single positive eigenvalue. Only ≈0.6% correspond to a more complex situation, with two or more real parts >0. The latter case, although only restricted to a small region within parameter space, holds profound implications for the possible dynamics of the system. As a further step within our approach, the existence of certain bifurcations of higher codimension allows the prediction of specific dynamics (see Materials and Methods). Fig. 7 shows a bifurcation diagram of the Calvin cycle within a particular region of parameter space where such bifurcations occur. Here, the system gives rise to a Gavrilov–Guckenheimer (GG) bifurcation, implying the existence of quasiperiodic dynamics and making the existence of chaotic dynamics likely. In close vicinity of the GG bifurcation, we also find a double Hopf (DH) bifurcation, formed by the interaction of two codimension1 HO bifurcations. The generic existence of a chaotic parameter region close to the DH bifurcation can be proven (25, 26).
Thus, our results demonstrate the possibility of quasiperiodic and chaotic dynamics for the model of the photosynthetic Calvin cycle, without relying on any particular assumptions about the functional form of the kinetic rate equations. Furthermore, because it is a quantitative method, we can assert that complex dynamics at the operating point are confined to a rather small region in parameter space and that the experimentally observed steady state is generically stable.
Discussion and Conclusions
We have presented a systematic approach to explore and quantify the dynamic capabilities of a metabolic system. Based on a parametric representation of the Jacobian matrix, constructed in such a way that each element is either directly experimentally accessible or amenable to a clear biochemical interpretation, we look for characteristic bifurcations that give insight into the possible dynamics of the system. Our method then builds on the construction of a large ensemble of models, encompassing all possible explicit kinetic models, to statistically explore and quantify the parameter regions associated with a specific dynamical behavior.
One of the primary advantages of our approach is that available information, such as experimentally accessible concentration values, can be readily incorporated into the description of the system. Focusing on a particular observed operating point, our approach then allows for the identification of crucial reaction steps that predominantly contribute to the stability, and thus robustness, of the observed state and results in specific biochemical conditions for which certain dynamical behavior can be expected. Furthermore, by taking bifurcations of higher codimension into account, we go beyond the usually considered case and are able to predict the possibility of complex or chaotic dynamics, often a nontrivial task even if an explicit kinetic model is available.
Exemplified here with two paradigmatic examples of metabolic pathways, our approach holds a vast number of possible further applications. In particular, with respect to biotechnological applications, a desired flux distribution must not necessarily be stable. By using our approach, it is thus possible to incorporate dynamic aspects into the description of the system and explore the conditions that support the stability of directed modifications of the system. Likewise, structural kinetic modeling can serve as a prequel to explicit mathematical modeling, aiming to identify crucial reaction steps and parameters in best time.
Materials and Methods
Interpretation of the Saturation Matrix.
Our approach relies crucially on the interpretation of the elements of the matrix θ _{x} ^{μ}. As a simple example, consider a single bilinear reaction rate of the form ν(S _{1}, S _{2}) = v _{max} S _{1} S _{2}. Then, according to Eq. 2 , the normalized rate is μ(x _{1}, x _{2}) = x _{1} x _{2}; thus, In the case of Michaelis–Menten kinetics ν(S) = v _{max} S/(K _{M} + S), depending on a single substrate S, we obtain Clearly, the partial derivative θ_{x} ^{μ} ∈ [0, 1] measures the normalized degree of saturation of the reaction at the steady state S ^{0}. The limiting cases are lim_{s0→0}θ_{x} ^{μ} = 1 (linear regime) and lim_{s0→∞}θ_{x} ^{μ} = 0 (full saturation). This result implies that the saturation parameter indeed covers the full interval, which holds likewise for the general case of Eq. 5 . For additional instances of specific rate functions, see the Supporting Appendix and Figs. 8–10.
Note that, except for the change in variables, the saturation parameters θ _{x} ^{μ} are reminiscent of the scaled elasticity coefficients, as defined in the realm of metabolic control analysis (6). However, for our reasoning to hold, the analysis is restricted to unidirectional reactions, i.e., in the case of reversible reactions, forward and backward terms have to be treated separately. Because the denominator is usually preserved for both terms, no additional free saturation parameters arise.
Another close analogy to the saturation parameters is found within the powerlaw approximation, where each enzyme kinetic rate law is replaced by a function of the form ν_{j}(S)=α_{j}∏_{i} S _{i} ^{gij} with g_{ij} denoting the “effective kinetic order” of the reaction (6). In fact, the powerlaw formalism can be regarded as the simplest possible way to specify explicit nonlinear functions that are consistent with a given Jacobian. Applying the transformation of Eq. 2 , we obtain μ_{j}(x) = ∏_{i} x _{i} ^{gij}, thus θ_{xi} ^{μj} = g _{ij}. However, beyond the properties of the Jacobian itself, only little confidence can be placed in an actual numerical integration of these functions (6). Generally, it is possible to specify several classes of explicit functions that, by construction, result in a given Jacobian but have no, or only little, biochemical justification otherwise. Consequently, we opt for using the properties of the parametric representation of the Jacobian directly, instead of going the loop way by means of auxiliary ad hoc functions.
Dynamics and Bifurcations.
One of the foundations of our approach is the fact that knowledge of the Jacobian matrix alone is sufficient to deduce certain characteristic bifurcations of a metabolic system. In general, the stability of a steady state is lost either in a HO bifurcation or in a bifurcation of SN type, both of codimension1. Of particular interest for revealing insights about the dynamical behavior of systems are bifurcations of higher codimension, such as the Takens–Bogdanov (TB), the GG, and the DH bifurcation (16, 25). Each of these local bifurcations of codimension2 arises out of an interaction of two codimension1 bifurcations and has important implications for the possible dynamical behavior. For instance, the TB bifurcation indicates the presence of a homoclinic bifurcation and therefore the possibility of spiking or bursting behavior. The presence of a GG bifurcation shows that complex (quasiperiodic or chaotic) dynamics exist generically in a certain parameter space. In the same way, the DH bifurcation indicates the generic existence of a chaotic parameter region. For details, see refs. 16 and 25, the Supporting Appendix, and Fig. 10.
Acknowledgments
R.S. was supported by the State Brandenburg Hochschulund Wissenschaftsprogramm 2004–2006. T.G. and B.B. were supported by German VWStiftung and Deutsche Forschungsgemeinschaft Grant SFB 555. T.G. was supported by the Alexander von Humboldt Foundation.
Footnotes
 ^{‡}To whom correspondence should be addressed. Email: steuer{at}agnld.unipotsdam.de

Author contributions: R.S., T.G., J.S., and B.B. designed research; R.S. performed research; T.G. contributed new reagents/analytic tools; and R.S., T.G., and B.B. wrote the paper.

Conflict of interest statement: No conflicts declared.

This paper was submitted directly (Track II) to the PNAS office.
 Abbreviations:
 DH,
 double Hopf;
 GG,
 Gavrilov–Guckenheimer;
 HK,
 hexokinase;
 HO,
 Hopf;
 PFK,
 phosphofructokinase;
 SN,
 saddlenode.

Freely available online through the PNAS open access option.
 © 2006 by The National Academy of Sciences of the USA
References
 ↵
 ↵

↵
 Fernie A. R. ,
 Trethewey R. N. ,
 Krotzky A. J. ,
 Willmitzer L.
 ↵

↵
 Palsson B. O. ,
 Joshi A. ,
 Ozturk S. S.

↵
 Heinrich R. ,
 Schuster S.

↵
 Hashimoto K. ,
 Tomita M. ,
 Takahashi K. ,
 Shimizu T. S. ,
 Matsuzaki Y. ,
 Miyoshi F. ,
 Saito K. ,
 Tanida S. ,
 Yugi K. ,
 Venter J. C. ,
 Hutchison C. A., III
 ↵
 ↵

↵
 Angeli D. ,
 Ferrell J. E., Jr ,
 Sontag E. D.
 ↵

↵
 Famili I. ,
 Förster J. ,
 Nielsen J. ,
 Palsson B. O.
 ↵
 ↵
 ↵

↵
 Gross T.
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵

↵
 Kuznetsov Y. A.
 ↵
Citation Manager Formats
More Articles of This Classification
Physical Sciences
Applied Mathematics
Biological Sciences
Related Content
 No related articles found.
Cited by...
 Quantifying uncertainty in partially specified biological models: how can optimal control theory help us?
 Plant Metabolic Modeling: Achieving New Insight into Metabolism and Metabolic Engineering
 How to predict community responses to perturbations in the face of imperfect knowledge and network complexity
 The stubborn roots of metabolic cycles
 A probabilistic approach to identify putative drug targets in biochemical networks
 Simulating Plant Metabolic Pathways with EnzymeKinetic Models
 Deciphering Transcriptional and Metabolic Networks Associated with Lysine Metabolism during Arabidopsis Seed Development
 Use of Randomized Sampling for Analysis of Metabolic Networks
 Formulating genomescale kinetic models in the postgenome era
 Getting to grips with the plant metabolic network
 The stability and robustness of metabolic states: identifying stabilizing sites in metabolic networks