- Research
- Open Access

# Hybrid dynamic/static method for large-scale simulation of metabolism

- Katsuyuki Yugi†
^{1}, - Yoichi Nakayama†
^{1}Email author, - Ayako Kinoshita
^{1}and - Masaru Tomita
^{1}

**2**:42

https://doi.org/10.1186/1742-4682-2-42

© Yugi et al; licensee BioMed Central Ltd. 2005

**Received:**25 April 2005**Accepted:**04 October 2005**Published:**04 October 2005

## Abstract

### Background

Many computer studies have employed either dynamic simulation or metabolic flux analysis (MFA) to predict the behaviour of biochemical pathways. Dynamic simulation determines the time evolution of pathway properties in response to environmental changes, whereas MFA provides only a snapshot of pathway properties within a particular set of environmental conditions. However, owing to the large amount of kinetic data required for dynamic simulation, MFA, which requires less information, has been used to manipulate large-scale pathways to determine metabolic outcomes.

### Results

Here we describe a simulation method based on cooperation between kinetics-based dynamic models and MFA-based static models. This hybrid method enables quasi-dynamic simulations of large-scale metabolic pathways, while drastically reducing the number of kinetics assays needed for dynamic simulations. The dynamic behaviour of metabolic pathways predicted by our method is almost identical to that determined by dynamic kinetic simulation.

### Conclusion

The discrepancies between the dynamic and the hybrid models were sufficiently small to prove that an MFA-based static module is capable of performing dynamic simulations as accurately as kinetic models. Our hybrid method reduces the number of biochemical experiments required for dynamic models of large-scale metabolic pathways by replacing suitable enzyme reactions with a static module.

## Keywords

- Static Module
- Hybrid Model
- Metabolic Flux Analysis
- Boundary Reaction
- Metabolic Control Analysis

## Background

Recent progress in high-throughput biotechnology [1–3] has made advances in understanding of cell-wide molecular networks possible at the systems level [4, 5]. To reconstruct cellular systems using the high-throughput data that are becoming available on their components, computer simulations are being revisited as an integrative approach to systems biology. Mathematical modelling of biochemical networks has been attempted since the 1960s, and before genome-scale pathway information became available, they mostly employed numerical integration of ordinary differential equations for reaction rates [6–10]. This kind of dynamic simulation model provides the time evolution of pathway properties such as metabolite concentration and reaction rate. To create accurate simulations, dynamic models require kinetic parameters and detailed rate-laws such as the MWC model [11] and those derived using the King-Altman method [12]. However, with few exceptions such as human erythrocyte metabolism [13, 14], it is virtually impossible to collect a complete set of kinetic properties for large-scale metabolic pathways. Therefore, the applicability of the dynamic method has been limited to relatively small pathways.

Another approach, such as metabolic flux analysis (MFA) using stoichiometric matrices, has been employed for large-scale analyses of metabolism [4, 15, 16]. Assuming a steady-state condition, MFA provides a flux distribution as the solution of the mass balance equation without the need for rate equations and kinetic parameters [16, 17]. Since it is a "static" approach, the ability of MFA to predict the dynamic behaviour of metabolic pathways is limited. It provides a snapshot of a certain pathway in a single state, but is insufficient to predict the dynamic behaviour of metabolism [18]. Recently, this approach was extended to allow the prediction of dynamic behaviour. This extension, dynamic flux balance analysis (DFBA) [19], provides optimal time evolution based on pre-defined constraints, including kinetic rate equations. However, this extension was not intended to reduce the masses of information necessary for developing dynamic cell-scale simulation models. In addition, this DFBA study did not define the criteria for segmenting a whole metabolic pathway into parts defined by kinetic rate equations and a stoichiometric model. Therefore this effort does not suffice as a generic modelling approach.

Here we propose a method for dynamic kinetic simulation of cell-wide metabolic pathways by applying the kinetics-based dynamic method to parts of a metabolic pathway and the MFA-based static method to the rest. Because the static module does not require any kinetic properties except the stoichiometric coefficients, this method can drastically reduce the number of enzyme kinetics assays needed to obtain the dynamic properties of the pathway. We have evaluated the accuracy of the hybrid method in comparison to a classical dynamic kinetic simulation using small virtual pathways and an erythrocyte metabolism model.

## Results

### Evaluation of errors

^{-16}of the maximum for the reaction rates. All the metabolite concentrations in the hybrid model were identical to those in the dynamic model (Table 1). When the concentration of metabolite A was increased two-fold, the hybrid and the dynamic models displayed similar time evolutions (Figure 3a and 3b). The maximum one-step errors after this perturbation were 4.000 × 10

^{-11}and 8.889 × 10

^{-6}for metabolite concentrations and reaction rates, respectively (Table 1).

Errors between the dynamic model and the hybrid model of the pathway shown in Fig. 2a. The maximum errors were measured within one numerical integration step. "Perturbation" denotes whether the errors were measured under a steady-state condition (-) or after a two-fold increase of metabolite A (+)

Perturbation | Maximum error (concentration) | Maximum error (reaction rate) | |
---|---|---|---|

Boundary | - | 0 | 0 |

+ | 8.000 × 10 | 0 | |

Static part | - | 0 | 8.592 × 10 |

+ | 4.000 × 10 | 8.889 × 10 |

^{-12}for metabolite concentrations and 2.837 × 10

^{-6}for reaction rates (Table 2). The time courses after a two-fold increase in the concentration of metabolite A were very similar between the dynamic and the hybrid model (Figure 3c and 3d). The maximum errors at the first integration step after the perturbation were 3.575 × 10

^{-7}for the metabolite concentrations and 0.00120 for the reaction rates.

Errors between the dynamic model and the hybrid model of the pathway shown in Fig. 2b. The maximum errors were measured within one numerical integration step. "Perturbation" denotes whether the errors were measured under a steady-state condition (-) or after a two-fold increase of metabolite A (+)

Perturbation | Maximum error (concentration) | Maximum error (reaction rate) | |
---|---|---|---|

Boundary | - | 5.049 × 10 | 5.609 × 10 |

+ | 3.575 × 10 | 1.323 × 10 | |

Static part | - | 7.176 × 10 | 2.837 × 10 |

+ | 1.192 × 10 | 0.00120 (E_CD) |

_{f}= 0.01s

^{-1}, k

_{r}= 0.001s

^{-1}to k

_{f}= 0.1s

^{-1}, k

_{r}= 0.091s

^{-1}.

### Correlation between elasticity and errors

### Application to erythrocyte metabolism

^{-4}, was observed in the reaction rate of 6-phosphogluconate dehydrogenase (6PGODH) (Table 3). (Note that this was true only when the gluconolactone-6-phosphate (GL6P) concentration was excluded. Owing to its small initial concentration (7.572 nM), the error in GL6P was sensitive to small changes and was associated with a large error of 0.00780.) The error in the 6PGODH rate remained the maximum error when the FDP concentration was perturbed.

Comparisons of the dynamic model and the hybrid model of the erythrocyte pathway shown in Fig. 5. The maximum errors were measured within one numerical integration step. "Perturbation" denotes whether the errors were measured under a steady-state condition (-) or after a three-fold increase of FDP concentration (+).

Perturbation | Maximum error (concentration) | Maximum error (reaction rate) | |
---|---|---|---|

Boundary | - | 7.796 × 10 | 1.555 × 10 |

+ | 1.153 × 10 | 3.020 × 10 | |

Static part | - | 1.111 × 10 | 2.170 × 10 |

+ | 4.282 × 10 | 2.170 × 10 |

When the boundary reaction was relocated from G6PDH, which forms a bottleneck of dynamic response in a transient state and has low elasticity at steady state, to phosphoglucoisomerase (PGI), which has a larger elasticity, the time courses calculated by the hybrid model were different from those produced by the dynamic model.

## Discussion

In the simulation experiments using hypothetical pathways and an erythrocyte model, the discrepancies between the dynamic and the hybrid models were sufficiently small to prove that an MFA-based static module is capable of performing dynamic simulations as accurately as a kinetic model. The key idea behind our method is to distinguish between dependent and independent variables (reactions). Although independent reactions can be affected by other dependent/independent reactions through effectors such as ADP in the phosphofructokinase reaction, the time evolution of adjacent reaction rates are mainly determined by independent reactions which constitute bottlenecks of dynamic behaviour in the metabolic network. Therefore, static modules should consist of only such dependent reactions, whereas dynamic modules can include both independent and dependent reactions. Our hybrid method reduces the number of biochemical experiments required for dynamic models of large-scale metabolic pathways by replacing suitable enzyme reactions with a static module. The optimal conditions for this method are (a) a system with few bottleneck reactions in order to enlarge the static modules, (b) small fluctuations in the reaction rates in static modules, and (c) accurately identifiable bottleneck reactions. How can such enzymes be identified? One obvious criterion for the enzymes to be suitably modelled by a static module is not to incorporate a bottleneck reaction in a transient state. Thus, the enzymes should not reach the maximum velocity quickly or be restrained at lower activities by allosteric regulation. Although the model comprising dynamic and static modules as a whole can represent transient states, it is assumed that the reactions in the static modules achieve or nearly achieve steady states within one numerical integration step. The existence of one or more bottleneck reactions in the static module may cause inconsistencies, because the hybrid method solves algebraic equations for static modules under a steady state assumption, although metabolites will be accumulated or depleted in real cells. Therefore, bottleneck reactions must be excluded from static modules. Another situation that should be avoided involves reaction rates in static modules that are affected by changes in enzyme concentration, such as those caused by changing levels of transcriptional/ post-transcriptional control. Such reactions should be included in dynamic modules.

A similar cause of inconsistency is the reversibility of reactions. Since the hybrid method assumes that reactions in the static module are reversible, inclusion of an irreversible step may cause inconsistencies, particularly in the presence of a perturbation downstream of the irreversible step (data not shown).

The accuracy of the calculation can also be affected by a time lag. In the static module of the hybrid model, time lags between the upstream and downstream reactions are not represented because the boundary reactions affect all subsequent reactions in the static module within one integration step regardless of the number of enzyme reactions. Depending on the simulation time scale, the static module should be limited to minimize the influence of time lags. This influence can be estimated by the ratio of elasticities, which can be an important criterion for including a reaction in the static module.

The correlation between elasticity and one-step error (Figure 4) indicates that, to ensure the accuracy of the simulation, the static module of a pathway should include reactions with larger elasticities and should be surrounded by boundary reactions with small elasticities. A large elasticity indicates that the enzyme is capable of changing its reaction rate rapidly in response to changes in substrate concentrations [21]. The result shown in Figure 4 demonstrates that enzymes with large elasticity contribute to the accuracy of the static module. On the other hand, boundary reactions with small elasticities, large substrate concentrations and/or small reaction rates change their activities little in response to substrate concentrations over a short period of time; perturbations are thus dampened by boundary reactions before being transmitted to the static modules. As a result, the reaction rates in the static modules do not change much after perturbations. Such a moderate time evolution allows even reactions that are not very fast to realize a reaction-rate distribution, **v**, that can be calculated from **v** = **S**^{#}**b** in as little as one numerical integration step. This allows the hybrid model to produce results that are in agreement with the dynamic model when the boundary reactions weaken perturbation.

The results we obtained when we relocated the boundary of the static module in the erythrocyte model support the importance of elasticity ratios. When G6PDH was included inside the static module, PGI became the new boundary reaction instead of G6PDH. The elasticity of PGI is large (elasticity = -452.496) compared to its neighbour G6PDH (elasticity = 0.0955). The relocated boundary is therefore composed of a pair of reactions that might produce unacceptable calculation errors, and in fact led to inconsistencies between the hybrid and dynamic models. Thus, the analytical conclusion presented in Eqs. (2) and (3) also holds for complex pathways, and elasticity provides a criterion for identifying groups of enzymes that can be approximated with sufficient accuracy by static modules. However, a large amount of experimental data is still required to determine the elasticities of all enzymatic reactions. In addition, the demarcation of the static module using elasticities determined by conventional biochemical experiments is unrealistic with respect to their throughput. Hence, the comprehensive determination of bottleneck reactions is the key task in the construction of large-scale metabolic pathway models using the hybrid method. Recent advances in flux measurement, quantitative metabolomics and proteomics allow large-scale measurement of flux distributions [22], intracellular metabolite concentrations and amounts of enzymes [23].

Recently, a method for high-throughput metabolomic analyses using capillary electrophoresis assisted by advanced mass spectrometry (CE-MS) and LC-MS/MS has been developed by the metabolomics group at our institute [24–27]. This technology allows us to determine the concentrations of more than 500 different metabolites quantitatively in a few hours. Furthermore, we are developing a method to calculate whole reaction rates of metabolic systems. This method has already achieved preliminary successes in determining the reaction rates of glycolysis in *E. coli* and human red blood cells. Pulse-chase analyses using ^{13}C labeled molecules and the CE-MS/LC-MS high-throughput system have also been used successfully by the same metabolomics group to determine fluxes in the *E. coli* central carbon pathway.

Several approaches have been proposed to quantify elasticity and other coefficients of metabolic control analysis from experimental data such as flux rates, metabolite concentrations or enzyme concentrations [28–31]. Thus, the hybrid method, in combination with the 'omics' data of metabolism, enables a dynamic kinetic simulation of cell-wide metabolism.

## Conclusion

Using this hybrid method, the cost of developing large-scale computer models can be greatly reduced since precise modelling with dynamic rate equations and kinetic parameters is limited to bottleneck reactions. This drastically reduces the number of experiments needed to obtain the kinetic properties required for the dynamic simulation of metabolic pathways.

## Methods

### Calculation procedure

The hybrid method works within one numerical integration step as follows: (i) all the reaction rates in the dynamic module are calculated from dynamic rate equations (V_{1}, V_{2}, V_{9}, and V_{10} in Figure 1); (ii) the reaction rate distribution in the static module (V_{3}, V_{4}, V_{5}, V_{6}, V_{7}, and V_{8}) is derived from the balance equation **Sv** = **b**, where **S** denotes the stoichiometric matrix, **v** the flux distribution, and **b** the rates of the dynamic exchange reactions at the system boundary (V_{2}, V_{9}, and V_{10}) that are calculated in step (i); and (iii) the concentrations of the metabolites (X_{1}-X_{13}) are determined by numerical integration of the reaction rates calculated in steps (i) and (ii). All the reactions in the static module are assumed to be reversible.

The calculation of the reaction rate distribution in the static module is similar to that in the MFA method. The only difference is that the exchange reactions between the dynamic and static modules are represented by kinetic rate equations instead of constant fluxes. In this study, we term a dynamic exchange reaction of a static module a "boundary reaction". Dynamic boundary reactions provide quasi-dynamic changes in the reaction rate distribution in the static module. The reaction rate distribution in the static module is calculated at every integration step that refers to the boundary reaction rates, which are determined by concentrations of metabolites inside and outside the static module. The time evolution of the metabolite concentration in the static module is calculated at every integration step by numerical integration of the reaction rates as well as the metabolites in the dynamic module.

In step (ii), the Moore-Penrose pseudo-inverse is employed to calculate the reaction rate distribution of the static module at each numerical integration step. This should result in a smaller computational cost than linear programming, which is commonly used to determine the flux distribution of the underdetermined system. When the linear equation **Sv** = **b** is determined, **S**^{#}, the Moore-Penrose pseudo-inverse of **S**, is identical to **S**^{-1}, the inverse of **S**. Thus, the reaction rate distribution of the static module is solved uniquely as **v** = **S**^{-1}**b**. If the equation **Sv** = **b** is over-determined, **v** = **S**^{#}**b** provides the least squares estimate of the reaction rate distribution [32] which minimizes |**Sv-b**|^{2}. Through this procedure, the error is distributed equally among the reaction rates of the static module.

In the case of an underdetermined static module, the solution was chosen from the solution space of the balance equation **Sv** = **b** to minimize the error of the ideal reaction rate distribution specified by the user. The optimal solution **v**_{best} is represented in Eq. (1) below [see Supplementary Text 1 in Additional file 1 for the derivation]:-

**v**_{best} = **i** + **S**^{#} (**b** - **Si**) (1)

where **v**_{best} is the closest solution to the ideal reaction rate distribution i in the solution space [Figure 6 in Additional file 1].

### Evaluation of errors at steady state

To compare the accuracy of the hybrid method with the conventional dynamic kinetic method analytically, we first employed a pathway model comprising the three sequential reactions shown below. The whole pathway is assumed to be at a steady-state.

In the remainder of this report, a "dynamic model" refers to a metabolic pathway model that is represented by kinetic rate equations only. Let v_{1}, v_{2} and v_{3} be the reaction rates of the three sequential reactions. In the hybrid model, the reaction rate v_{2} was represented as a static module of this pathway. When the concentration of metabolite A, the substrate of v1, is perturbed, the discrepancy between v2 in the hybrid model and v2 in the dynamic model is as described below [see Supplementary Text 2 in Additional file 1 for the derivation]:

where v_{2d}, v_{2k}, [A], [B], *ε*^{v1}_{A} and *ε*^{v1}_{B} denote the reaction rate v_{2} in the dynamic model, v_{2} in the hybrid model, concentration of metabolite A, concentration of metabolite B, elasticity of v_{1} with respect to metabolite A, and elasticity of v2, respectively. The variables with Δ are increments after a small time step Δt. The parameter p represents a ratio of the reaction rate in the static module to the influx, as in Δ*v*_{2h} = *p* Δ*v*_{1}. The ratio p is determined by the stoichiometric matrix of the pathway.

In Eq. (2), the left bracket term on the right-hand side indicates the magnitude of the perturbation transmitted to the static module. This term indicates that the error between the hybrid and dynamic models is proportional to the increment of metabolites and the elasticity of the boundary reactions. The right bracket describes the susceptibility of the reaction rate v_{2} to v_{1}. When *ε*^{v2}_{B} satisfies the relationship below, v_{2} in the hybrid model exhibits identical time evolution to the dynamic model:

Since a small Δt (<<1.0s) is usually employed for accurate simulations of metabolic pathways, Eq. (3) implies that a reaction with large elasticity can be appropriately replaced by a static module.

For more complex pathways, such a theoretical analysis is not practical because large numbers of variables and parameters might impede clear discussions. Instead, simulation experiments were performed to compare the accuracy of hybrid models with dynamic models by numerical methods.

The accuracy of the hybrid model was evaluated numerically in comparison with a conventional kinetic model of the same metabolic pathway under two conditions: a steady-state condition and a time evolution after a two-fold increase of metabolites that are catalyzed by boundary reactions. The errors under steady-state conditions were employed as controls to evaluate discrepancies in dynamic behaviour. These computer simulations were performed using the E-Cell Simulation Environment version 1.1 or 3.1.102 for RedHat Linux 9.0/i386. The errors of reaction rates and metabolite concentrations were measured as below:-

where v_{d} and v_{h} denote either the reaction rates or the concentrations in the dynamic and hybrid models, respectively. The values of v_{d} and v_{h} were taken at the first numerical integration step, in which the concentration increase influences the initial steady-state values of the reaction rates and metabolite concentrations. In this article, this is termed "one-step error". We used one-step errors to represent the discrepancies between the two simulation methods in transient dynamics.

The one-step errors were evaluated using two simple pathways; the static module of one is determined, while the other is underdetermined (Figure 2a and 2b). All the reaction rates in these simple pathways were represented as v = k_{f}[S]-k_{r}[P] where v, k_{f}, k_{r}, [S], and [P] are a reaction rate, a forward rate constant, a reverse rate constant, a substrate concentration and a product concentration, respectively. In the pathway with the over-determined static module, the rate constants were k_{f} = 0.05s^{-1} and k_{r} = 0.091s^{-1} for E_CD and E_CE, k_{f} = 0.1s^{-1} and k_{r} = 0.091s^{-1} for E_DF and E_EG, and k_{f} = 0.01s^{-1} and k_{r} = 0.001^{-1} for the other reactions in the pathway of Figure 2a. The initial metabolite concentrations were 1.0 mM for A, B and C, and 0.5 mM for the other metabolites. Metabolite A was increased two-fold to evaluate the errors in transient dynamics. In the pathway with an underdetermined static module, the kinetic parameters were k_{f} = 0.01s^{-1} and k_{r} = 0.001s^{-1} for E_AB, E_BC, and E_FG; k_{f} = 0.1s^{-1} and k_{r} = 0.098s^{-1} for E_CD and E_DF; k_{f} = 0.1s^{-1} and k_{r} = 0.097s^{-1} for E_CE and E_EF; and k_{f} = 0.1s^{-1} and k_{r} = 0.96s^{-1} for E_CF. The steady-state flux distribution was employed for the ideal reaction rate distribution in the static module; the ideal reaction rates were 2 *μ* M/s for E_CD and E_DF, 3 *μ* M/s for E_CE and E_EF, and 4 *μ* M/s for E_CF. All the initial metabolite concentrations were 1.0 mM. The concentration of metabolite A was increased two-fold to evaluate the error.

### Correlation between elasticity and error

Elasticity is a coefficient used to quantify the sensitivity of the enzyme to its substrates and is defined as below in the context of metabolic control analysis [21]:

where [S] and v denote the substrate concentration and the reaction rate of the enzyme, respectively. Correlation between the one-step errors and elasticities of each enzyme at a steady state was examined using a linear pathway and a glycolysis model [13, 20] (Figure 2c and 2d, respectively). In the linear pathway model, the reaction rate v is represented by the same equation as in the two hypothetical models above. The kinetic parameters were k_{f} = 0.01s^{-1} and k_{r} = 0.009s^{-1} for E_AB, E_BC, and E_FG and k_{f} = 0.1s^{-1} and k_{r} = 0.099s^{-1} for E_CD, E_DE, and E_EF. All the initial metabolite concentrations were 1.0 mM. The two rate constants of reactions E_BC, E_CD, E_DE, and E_DF were altered within the range 0.01<k_{f}<1.0. The value of k_{r} was determined to satisfy k_{f}-k_{r} = 0.01 to sustain the initial steady-state concentrations. The concentration of metabolite A was increased two-fold to evaluate the errors. For error measurements in the glycolysis model, each enzymatic reaction was replaced, one by one, with a static module. The substrate concentrations of the boundary reactions were increased three-fold.

### Application to erythrocyte metabolism

A cell-wide model of erythrocyte metabolism [14] was employed to evaluate the applicability of the hybrid method in a more realistic and complex pathway. This erythrocyte model reproduces steady-state metabolite concentrations similar to experimental data. The static region was determined using a ratio of elasticities as below:

where *ε*^{b} and *ε*^{x} denote the elasticities of a boundary reaction and of reaction X, respectively. All the elasticities of the model were calculated by numerical differentiation of each rate equation. A group of enzymes with small r values were regarded as appropriate candidates for inclusion in a static module. The concentration of fructose-1,6-diphosphate (FDP) was increased three-fold to measure the errors in dynamic behaviours.

## Notes

## Declarations

### Acknowledgements

We thank Nobuyoshi Ishii for insightful discussions; Yoshihiro Toya for the preparation of one of the small virtual pathway models; Pawan Kumar Dhar, Yasuhiro Naito, Shinichi Kikuchi and Kazuharu Arakawa for critically reading the manuscript; and Kouichi Takahashi for providing technical advice. This work was supported in part by a grant from Leading Project for Biosimulation, Keio University, The Ministry of Education, Culture, Sports, Science and Technology (MEXT); a grant from CREST, JST; a grant from New Energy and Industrial Technology Development and Organization (NEDO) of the Ministry of Economy, Trade and Industry of Japan (Development of a Technological Infrastructure for Industrial Bioprocess Project); and a grant-in-aid from the Ministry of Education, Culture, Sports, Science and Technology for the 21 st Century Centre of Excellence (COE) Program (Understanding and Control of Life's Function via Systems Biology).

## Authors’ Affiliations

## References

- Blattner FR, Plunkett GIII, Bloch CA, Perna NT, Burland V, Riley M, Collado-Vides J, Glasner JD, Rode CK, Mayhew GF, Gregor J, Davis NW, Kirkpatrick HA, Goeden MA, Rose DJ, Mau B, Shao Y: The complete genome sequence of Escherichia coli K-12. Science. 1997, 277 (5331): 1453-1462. 10.1126/science.277.5331.1453.View ArticlePubMedGoogle Scholar
- Fiehn O, Kopka J, Dormann P, Altmann T, Trethewey RN, Willmitzer L: Metabolite profiling for plant functional genomics. Nature Biotechnology. 2000, 18 (11): 1157-1161. 10.1038/81137.View ArticlePubMedGoogle Scholar
- Wang Y, Liu CL, Storey JD, Tibshirani RJ, Herschlag D, Brown PO: Precision and functional specificity in mRNA decay. Proceedings of the National Academy of Sciences of the United States of America. 2002, 99 (9): 5860-5865. 10.1073/pnas.092538799.PubMed CentralView ArticlePubMedGoogle Scholar
- Edwards JS, Palsson BO: The Escherichia coli MG1655 in silico metabolic genotype: its definition, characteristics, and capabilities. Proceedings of the National Academy of Sciences of the United States of America. 2000, 97 (10): 5528-5533. 10.1073/pnas.97.10.5528.PubMed CentralView ArticlePubMedGoogle Scholar
- Shen-Orr SS, Milo R, Mangan S, Alon U: Network motifs in the transcriptional regulation network of Escherichia coli. Nature Genetics. 2002, 31: 64-68. 10.1038/ng881.View ArticlePubMedGoogle Scholar
- Bakker BM, Michels PA, Opperdoes FR, Westerhoff HV: What controls glycolysis in bloodstream form Trypanosoma brucei?. Journal of Biological Chemistry. 1999, 274 (21): 14551-14559. 10.1074/jbc.274.21.14551.View ArticlePubMedGoogle Scholar
- Barkai N, Leibler S: Robustness in simple biochemical networks. Nature. 1997, 387 (6636): 913-917. 10.1038/43199.View ArticlePubMedGoogle Scholar
- Bhalla US, Iyengar R: Emergent properties of networks of biological signaling pathways. Science. 1999, 283 (5400): 381-387. 10.1126/science.283.5400.381.View ArticlePubMedGoogle Scholar
- Cornish-Bowden A, Cardenas ML: Information transfer in metabolic pathways: effects of irreversible steps in computer models. European Journal of Biochemistry. 2001, 268 (24): 6616-6624. 10.1046/j.0014-2956.2001.02616.x.View ArticlePubMedGoogle Scholar
- Chance B, Garfinkel D, Higgins J, Hess B: Metabolic control mechanisms V: a solution for the equations representing interaction between glycolysis and respiration in ascites tumor cells. Journal of Biological Chemistry. 1960, 235 (8): 2426-2439.PubMedGoogle Scholar
- Monod J, Wyman J, Changeux JP: On the nature of allosteric transitions: a plausible model. Journal of Molecular Biology. 1965, 12: 88-118.View ArticlePubMedGoogle Scholar
- King EL, Altman C: A schematic method of deriving the rate laws for enzyme catalyzed reactions. Journal of Physical Chemistry. 1956, 60: 1375-1378. 10.1021/j150544a010.View ArticleGoogle Scholar
- Joshi A, Palsson BO: Metabolic dynamics in the human red cell: part I. a comprehensive kinetic model. Journal of Theoretical Biology. 1989, 141 (4): 515-528.View ArticlePubMedGoogle Scholar
- Ni TC, Savageau MA: Model assessment and refinement using strategies from biochemical systems theory: application to metabolism in human red blood cells. Journal of Theoretical Biology. 1996, 179 (4): 329-368. 10.1006/jtbi.1996.0072.View ArticlePubMedGoogle Scholar
- Henriksen CM, Christensen LH, Nielsen J, Villadsen J: Growth energetics and metabolic fluxes in continuous cultures of Penicillium chrysogenum. Journal of Biotechnology. 1996, 45: 149-164. 10.1016/0168-1656(95)00164-6.View ArticlePubMedGoogle Scholar
- Ibarra RU, Edwards JS, Palsson BO: Escherichia coli K-12 undergoes adaptive evolution to achieve in silico predicted optimal growth. Nature. 2002, 420 (6912): 186-189. 10.1038/nature01149.View ArticlePubMedGoogle Scholar
- Aiba S, Matsuoka M: Identification of metabolic model: citrate production from glucose by Candida lipolytica. Biotechnology and Bioengineering. 1979, 21 (8): 1373-1386. 10.1002/bit.260210806.View ArticleGoogle Scholar
- Varner J, Ramkrishna D: Mathematical models of metabolic pathways. Current Opinion in Biotechnology. 1999, 10 (2): 146-150. 10.1016/S0958-1669(99)80025-1.View ArticlePubMedGoogle Scholar
- Mahadevan R, Edwards JS, Doyle FJ: Dynamic flux balance analysis of diauxic growth in Escherichia coli. Biophysical Journal. 2002, 83 (3): 1331-1340.PubMed CentralView ArticlePubMedGoogle Scholar
- Mulquiney PJ, Kuchel PW: Model of 2,3-bisphosphoglycerate metabolism in the human erythrocyte based on detailed enzyme kinetic equations: equations and parameter refinement. Biochemical Journal. 1999, 342 (Pt 3): 581-596. 10.1042/0264-6021:3420581.PubMed CentralView ArticlePubMedGoogle Scholar
- Fell DA: Metabolic control analysis: a survey of its theoretical and experimental development. Biochemical Journal. 1992, 286: 313-330.PubMed CentralView ArticlePubMedGoogle Scholar
- Wittmann C, Heinzle E: Genealogy profiling through strain improvement by using metabolic network analysis: metabolic flux genealogy of several generations of lysine-producing corynebacteria. Applied and Environmental Microbiology. 2002, 68 (12): 5843-5859. 10.1128/AEM.68.12.5843-5859.2002.PubMed CentralView ArticlePubMedGoogle Scholar
- Ghaemmaghami S, Huh WK, Bower K, Howson RW, Belle A, Dephoure N, O'Shea EK, Weissman JS: Global analysis of protein expression in yeast. Nature. 2003, 425 (6959): 737-741. 10.1038/nature02046.View ArticlePubMedGoogle Scholar
- Soga T, Kakazu Y, Robert M, Tomita M, Nishioka T: Qualitative and quantitative analysis of amino acids by capillary electrophoresis-electrospray ionization-tandem mass spectrometry. Electrophoresis. 2004, 25 (13): 1964-1972. 10.1002/elps.200305791.View ArticlePubMedGoogle Scholar
- Soga T, Ohashi Y, Ueno Y, Naraoka H, Tomita M, Nishioka T: Quantitative metabolome analysis using capillary electrophoresis mass spectrometry. J Proteome Res. 2003, 2 (5): 488-494. 10.1021/pr034020m.View ArticlePubMedGoogle Scholar
- Soga T, Ueno Y, Naraoka H, Ohashi O, Tomita M, Nishioka T: Simultaneous determination of anionic intermediates for Bacillus subtilis metabolic pathways by capillary electrophoresis electrospray Ionization mass spectrometry. Analytical Chemistry. 2002, 74: 2233-2239. 10.1021/ac020064n.View ArticlePubMedGoogle Scholar
- Ishii N, Robert M, Nakayama Y, Kanai A, Tomita M: Toward large-scale modeling of the microbial cell for computer simulation. J Biotechnol. 2004, 113 (1-3): 281-294. 10.1016/j.jbiotec.2004.04.038.View ArticlePubMedGoogle Scholar
- Delgado J, Liao JC: Determination of flux control coefficients from transient metabolite concentrations. Biochemical Journal. 1992, 282 (Pt 3): 919-927.PubMed CentralView ArticlePubMedGoogle Scholar
- Acerenza L, Cornish-Bowden A: Generalization of the double-modulation method for in situ determination of elasticities. Biochemical Journal. 1997, 327 (Pt1): 217-223.PubMed CentralView ArticlePubMedGoogle Scholar
- de la Fuente A, Snoep JL, Westerhoff HV, Mendes P: Metabolic control in integrated biochemical systems. European Journal of Biochemistry. 2002, 269 (18): 4399-4408. 10.1046/j.1432-1033.2002.03088.x.View ArticlePubMedGoogle Scholar
- Wu L, Wang W, van Winden WA, van Gulik WM, Heijnen JJ: A new framework for the estimation of control parameters in metabolic pathways using lin-log kinetics. European Journal of Biochemistry. 2004, 271 (16): 3348-3359. 10.1111/j.0014-2956.2004.04269.x.View ArticlePubMedGoogle Scholar
- Stephanopoulos GN, Aristidou AA, Nielsen J: Metabolic Engineering: Principles and Methodologies. 1998, San Diego , Academic PressGoogle Scholar

## Copyright

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.