 Research
 Open Access
 Published:
Distinguishing enzymes using metabolome data for the hybrid dynamic/static method
Theoretical Biology and Medical Modellingvolume 4, Article number: 19 (2007)
Abstract
Background
In the process of constructing a dynamic model of a metabolic pathway, a large number of parameters such as kinetic constants and initial metabolite concentrations are required. However, in many cases, experimental determination of these parameters is timeconsuming. Therefore, for largescale modelling, it is essential to develop a method that requires few experimental parameters. The hybrid dynamic/static (HDS) method is a combination of the conventional kinetic representation and metabolic flux analysis (MFA). Since no kinetic information is required in the static module, which consists of MFA, the HDS method may dramatically reduce the number of required parameters. However, no adequate method for developing a hybrid model from experimental data has been proposed.
Results
In this study, we develop a method for constructing hybrid models based on metabolome data. The method discriminates enzymes into static modules and dynamic modules using metabolite concentration time series data. Enzyme reaction rate time series were estimated from the metabolite concentration time series data and used to distinguish enzymes optimally for the dynamic and static modules. The method was applied to build hybrid models of two microbial centralcarbon metabolism systems using simulation results from their dynamic models.
Conclusion
A protocol to build a hybrid model using metabolome data and a minimal number of kinetic parameters has been developed. The proposed method was successfully applied to the strictly regulated centralcarbon metabolism system, demonstrating the practical use of the HDS method, which is designed for computer modelling of metabolic systems.
Background
Since a biochemical network is essentially a nonlinear, nonequilibrium, nonsteadystate system, dynamic simulation is especially effective for analyzing or predicting its behaviour in a detailed and realistic manner. However, a large amount of experimental information, including reaction mechanisms of enzymes, kinetic constants, and initial concentrations of enzymes and metabolites, is required to construct a dynamic model of a metabolic pathway. Although a number of highthroughput technologies for obtaining comprehensive biochemical data have been developed [1–6], most experimental methods for determining enzyme kinetics are of the lowthroughput variety. Recently, several databases for enzyme kinetics have been published on the internet [7–9]. However, in many cases, the parameters in these databases are insufficient for building an accurate metabolic model. Moreover, although intracellular data can be collected from the published literature, experimental conditions and target strains are, in general, not uniform. Therefore, a huge amount of experimental work is currently needed to build an accurate dynamic model of a biochemical system. For this reason, a modelling method requiring less experimental effort needs to be developed.
Yugi et al. proposed a novel method for dynamic modelling of metabolism, the hybrid dynamic/static method (HDS method) [10]. The HDS method divides a dynamic system into a dynamic module and a static module. Enzyme reactions included in the dynamic module are represented by differential equations. Reaction rates of enzymes included in the static module are calculated by metabolic flux analysis (MFA) [11, 12]. Since MFA needs no kinetic information, the amount of experimental work required is dramatically reduced. According to Okino and Mavrovouniotis's classification [13], the HDS method can be regarded as a "linear transformation into standard twotimescale form," which is a timescale analysis method. The superior points of the HDS method are its simple architecture and the admissibility of multiple metabolites. Only relationships among enzyme reactions are employed in the HDS method; thus a model builder does not have to consider the problem of multiple timescale reactions of a given metabolite [14]. Since the MoorePenrose pseudoinverse [15, 16] of the stoichiometric coefficient matrix for the unknown variables (i.e. reaction rates of enzymes in the static module) is applied in performing the MFA, the stoichiometric coefficient matrix for the unknown variables does not have to be square and regular.
Although the HDS method has the aforementioned advantages, no method has been proposed for splitting a dynamic system into a dynamic module and a static module before completion of the initial model construction. Advanced measurement technologies have been developed that now enable researchers to obtain the metabolome, that is, comprehensive metabolite concentration data [17–19]. It is reasonable to expect that the indepth information of the metabolome contributes to the process for distinguishing dynamic and static enzymes in a metabolic system. In this study, we have developed a method of distinguishing dynamic and static enzymes based on metabolome data before construction of a complete model. The purpose of the proposed method is to provide the information (distinguishing dynamic from static enzymes) for initial HDS model construction required by the model builders without losing the advantage of the HDS method: reducing experimental efforts to obtain kinetic information of the modelled metabolic system. Identification of enzyme kinetic rate equations and the fitting of kinetic parameters using metabolite concentration data are outside the scope of this study. Moreover, biological meanings of the dynamic/static modules are not considered explicitly in the HDS method.
The proposed method consists of two parts. First, the enzyme reaction rate time series are estimated from metabolite concentration time series data. The dynamic and static enzymes are distinguished using the estimated enzyme reaction rate time series. The purpose of this study was to confirm that the proposed method can be used to construct accurate hybrid models, with accuracy comparable to that of a fully dynamic model. Therefore, we used pseudoexperimental data obtained from preliminarily constructed fully dynamic models. Two models of microorganisms, Escherichia coli [20] and Saccharomyces cerevisiae [21], were used for evaluation.
Methods
Hybrid dynamic/static method
The hybrid dynamic/static method (HDS method) is described in Yugi et al. [10]. Enzyme reaction rates in the static module are calculated by the following equation:v_{ static }(t) = S_{ static }^{#} · S_{ module boundary }· v_{ module boundary }(t)
where v_{static} is the static module enzyme reaction rate vector, v_{module boundary} is the module boundary enzyme reaction rate vector, S_{static}^{#} is the MoorePenrose pseudoinverse of the stoichiometric coefficient matrix for enzymes in the static module, and S_{module boundary} is the stoichiometric coefficient matrix for module boundary enzymes. The HDS method aims to describe a system in which a quasisteady state is attained in the static module at each instant, while the overall system (both the dynamic and the static modules) acts dynamically [10]. A transient value of the modelled system is calculated by an interaction between kineticbased dynamic models and MFAbased static models.
Estimation of internal enzyme reaction rates
To calculate reaction rates of enzymes from metabolite concentrations, we define a "system boundary enzyme" as an enzyme located on the border of the metabolic system and extending outside the system. The system boundary enzyme is not the same as the "module boundary enzyme" defined by Yugi et al. [10]. A non system boundary enzyme is defined as an "internal enzyme." The relationship among the dynamic module, static module, module boundary enzyme, system boundary enzyme, and internal enzyme is shown in Figure 1. Since all system boundary enzymes should be included in the dynamic module, we assumed that the kinetics of system boundary enzymes have already been determined and that the reaction rates of system boundary enzymes can be calculated from metabolite concentrations.
The reaction rates of internal enzymes were calculated from the slopes of metabolite concentrations and the reaction rates of the system boundary enzymes. With the definitions of the system boundary enzyme and the internal enzyme, a mass balance equation of a metabolic system under a dynamic transition state can be expressed as follows [22]:
where S is the stoichiometric coefficient matrix, diag(1) is a diagonal matrix (column number = row number = metabolite number), v_{system boundary}(t) is the system boundary enzyme reaction rate vector, v_{internal}(t) is the internal enzyme reaction rate vector, and C'(t) is the metabolite concentration slope vector.
If C(t) and v_{system boundary}(t) are known, the reaction rates of the internal enzymes can be estimated from Eq. (3), which is transformed from Eq. (2).
where S_{internal}^{#} is the MoorePenrose pseudoinverse of the stoichiometric coefficient matrix for internal enzymes, and S_{system boundary} is the stoichiometric coefficient matrix for system boundary enzymes [see Supplementary Text (see additional file 1) for an example of this procedure]. This procedure uses only the mass balance of the overall system and rate equations of the system boundary enzymes; thus, no information about regulation in the internal system is required beforehand. When Eq. (3) is applied to a determined system, the equation provides a true solution for v_{internal}, and when Eq. (3) is applied to an overdetermined system, the leastsquares estimation of v_{internal} is obtained [10]. In both cases, the solution is reasonable even if the modelled metabolic system has a complex network [10]. When Eq. (3) is applied to an underdetermined system, the equation provides the least norm solution. However, such a least norm solution is not always a physiologically optimal estimation of v_{internal}. This is a limitation of the current procedure.
Evaluation of estimated internal enzyme reaction rates
The accuracy of the estimated internal enzyme reaction rates was evaluated by means of the reproduced metabolite concentration time series, which were calculated from the estimated enzyme reaction rates. Since it is difficult to compare the true and estimated reaction rates, we compared the metabolic concentrations. If an enzyme catalyzes a reversible reaction, the sign of the sum of the forward and reverse reaction rates may change. Near such a sign change, the calculated relative error between the true reaction rate and the estimated reaction rate may at times be a very large value (see Eq. (4) below). When the value of a data point is close to zero, a large error will be obtained. However, in general, most metabolite concentrations have a sufficiently large positive value for the problem caused by a value close to zero to be avoided.The metabolite concentration time series slope was calculated from the reaction stoichiometric matrix and each estimated enzyme reaction rate time series. The metabolite concentration time series was calculated by numerical integration of the metabolite concentration slope time series obtained. The mean relative error (MRE) [23] between the true values (data) and the calculated values in the metabolite concentration time series was calculated by the following equation:
where C_{data,i,j}is the true concentration of the jth metabolite at the ith sampling point, C_{estimated,i,j}is the estimated (reproduced) concentration of the jth metabolite at the ith sampling point, n_{metabolite} is the number of metabolites, and n_{sampling point} is the number of sampling points.
In this study, the MRE between the true metabolite concentration data and the reproduced metabolite concentrations is called the "basal error".
Distinction of dynamic and static enzymes
The genetic algorithm (GA) [24] is employed to search for an optimal dynamic/static enzyme combination in a metabolic system. In this work, an individual code set for the GA was defined to represent the dynamic/static enzymes in a metabolic system. For example, DDSSDD represents a metabolic system consisting of six enzymes: the 1st, 2nd, 5th, and 6th enzymes for the dynamic module and the 3rd and 4th enzymes for the static module. In the GA calculation, the enzyme reaction rate time series in the static module were calculated from enzyme reaction rate time series in the dynamic module, which were derived from metabolite concentration time series data, by the same HDS method. Consequently, each metabolite concentration time series data point was calculated by the same method as that described in "Evaluation of estimated internal enzyme reaction rates". The fitness function defined in Eq. (5) was calculated for each code set; thereafter, propagation, crossover, and mutation followed. This procedure was repeated until the optimal solution, which minimizes Eq. (5), was found. A flowchart of the process for distinguishing dynamic/static enzymes is shown in Figure 2.
where C_{data,i,j}is true concentration of the jth metabolite at the ith sampling point, C_{estimated,i,j}is estimated concentration of the jth metabolite at the ith sampling point, n_{metabolite} is number of metabolites, n_{sampling point} is number of sampling points, n_{enzyme} is number of internal enzymes, n_{static enzyme} is number of enzymes included in static module, and w is weighting coefficient.
The first term in the fitness function represents the average error of the metabolite concentrations. For the fitness function, for the same reason as in the evaluation of estimated enzyme reaction rates, the metabolite concentrations rather than the enzyme reaction rates themselves were used. The second term in the fitness function evaluates the ratio of static enzymes included in the metabolic system; this term was added to adjust the number of enzymes in the static modules. The second term is multiplied by an adjusting parameter, a weighting coefficient, to control the balance between the model error and the static enzyme ratio.
9 different values of the weighting coefficient (w = 1.000, 0.750, 0.500, 0.250, 0.100, 0.075, 0.050, 0.025, and 0.010) were employed. The results of distinguishing dynamic and static enzymes were used to construct the hybrid models.
Error calculation
MRE of the metabolite concentration time series in a result of the process for distinguishing dynamic/static enzymes or in a hybrid model was calculated by Eq. (4). Finally, in the process for distinguishing dynamic/static enzymes, the "basal error", which originated from the incompleteness of the estimation of the enzyme reaction rates and from the error of the numerical integration of the enzyme reaction rates, rather than from the HDS calculation, was subtracted from the MRE.
Pseudo experiments
Two microbial centralcarbon metabolism models were chosen for testing: the E. coli model constructed by Chassagnole et al. [20] and the S. cerevisiae model constructed by Hynne et al. [21]. For the E. coli model, starting from a steady state for which the extracellular glucose concentration was 5.56 × 10^{2} mM, a glucose pulse was added. The concentration of the injected glucose pulse was 1.67 mM. In Chassagnole's original model, time series of nucleotides (ATP/ADP/AMP, NAD(H), NADP(H)) were expressed by timedependent functions [20]. However, in our study, the nucleotide concentrations were fixed as initial values. For the S. cerevisiae model, starting from a steady state for which the glucose concentration in the feed solution was 2.50 mM, the glucose concentration was shifted to 5.00 mM. The metabolite concentrations in both models at the steady state – that is, the initial concentrations for the dynamic simulations – are shown in Table S1 (see additional file 1). The running time after perturbation was set to 20 s for the E. coli model and 60 s for the S. cerevisiae model; these settings were chosen to allow time for the change from the original steady state to another steady state after the perturbation. The calculated metabolite concentration time series data sets were obtained at intervals of 1 s. These data sets were used as noisefree pseudometabolome data to calculate the slopes of the metabolite concentrations (C'(t)) and the reaction rates of the system boundary enzymes (v_{system boundary}(t)) in a performance test of the method. The slopes of the metabolite concentrations were obtained by firstorder differentiation of the interpolated metabolite concentration time series.
Noise addition to the pseudoexperimental data
To evaluate the practical use of the proposed method, artificial noise was added to each pseudoexperimental metabolite concentration data point. The coefficient of variance (CV) was assumed to be 15%, and the standard deviation (SD) of each pseudoexperimental data point was calculated by multiplying the CV by the noisefree value. A normally distributed random number around the noisefree value was generated for each data point using the SD obtained. Five noiseadded data points were generated for each noisefree data point as pseudoreplicated measurements. The average of the five noiseadded data points was used in the following smoothing procedure.
Smoothing of noisy pseudoexperimental data
Each noiseadded metabolite concentration time series pseudo data set was smoothed by fitting it to a polynomial or a rational function of time using the leastsquares method.
Calculation tools
MATLAB Release 2006a (MathWorks) was used for all calculations. Ordinary differential equations were solved by the ODE15s algorithm [25]. For interpolation, differentiation and smoothing of the metabolite concentration time series data, Curve Fitting Toolbox 1.1.5 (MathWorks) was used. Cubic spline interpolation was employed. For optimization, the Genetic Algorithm and Direct Search Toolbox 2.0.1 (MathWorks) was employed. In each GA calculation, the number of code set was set to 100. The other parameters were set to default values. Each optimal solution was taken after the fitness function converged to a constant value.
Results
Estimation of enzyme reaction rates using noisefree data
In the HDS method, reaction rates of enzymes in a dynamic module are used to estimate reaction rates of enzymes in a static module. If the true reaction rates of all enzymes in a metabolic system are known, they can be used directly for discriminating dynamic and static enzymes. However, the true reaction rates of enzymes in a cell cannot be determined in most cases. Therefore, we tried to estimate the reaction rates of enzymes from metabolite concentrations, which can be experimentally measured by highthroughput metabolome technologies.
We calculated the estimated reaction rates by using the metabolite concentration time series obtained from the E. coli and S. cerevisiae models to evaluate our method of estimating reaction rates. In this section, the noisefree pseudoexperimental data were used to obtain a clear assessment of the estimation method itself. In the true reaction rate time series of Tkb in E. coli, TA in E. coli, and AK in S. cerevisiae, some signchanging points were observed (Figure S1, see additional file 1). As predicted, around such points, huge relative errors between the true enzyme reaction rates and the estimated enzyme reaction rates were calculated (Figure S1). To avoid the undesired influence of such huge errors caused by using the reaction rates themselves, the reproduced metabolite concentrations were employed for the evaluation, as explained in the Methods. Therefore, the accuracy of the estimated reaction rates of the internal enzymes was assessed by the MRE between the original metabolite concentration time series and the reproduced metabolite concentration time series (Table 1). In the results for E. coli, the MRE was relatively large, mainly because of the large error in PGP. Errors in metabolites except for PGP were within approximately 10%; thus the estimation can be considered practically meaningful. For S. cerevisiae, errors of all metabolites were sufficiently small. On the whole, enzyme reaction rate time series data can be estimated from metabolite concentration time series data.
Distinction of dynamic and static enzymes using noisefree data
Using enzyme reaction rate time series data, we can apply the HDS method to calculate the reaction rates of static enzymes from the reaction rates of dynamic enzymes. These calculated static enzyme reaction rates can then be compared with the original reaction rate data. The errors between the estimated static enzyme reaction rates and the static enzyme reaction rate data can be used to find an optimal pattern for distinguishing dynamic from static enzymes. In this study, a fitness function (Eq. (5)) consisting of two terms was used for the optimization. In Eq. (5), the second term is multiplied by an adjusting parameter, a weighting coefficient (w). Even if the same data set is used, the result for distinguishing dynamic/static enzymes may vary for different w.
The E. coli and S. cerevisiae models and the estimated reaction rates obtained in the previous section (i.e., calculated from noisefree metabolite concentration data) were used to test this method for distinguishing enzymes, and the optimized patterns of dynamic and static enzymes shown in Table 2 were obtained as a result. As expected, the proportion of static enzymes decreased with decreasing w. The dynamic/static enzymes displayed on the metabolic map are shown in Supplementary Figure S2 (see additional file 1). The results obtained by using the noiseadded metabolite concentration data are shown in the following section.
In the next step, the estimated optimal results for distinguishing dynamic/static enzymes in Table 2 were used to convert the full dynamic models for E. coli and S. cerevisiae to hybrid models. In a process for distinguishing dynamic/static enzymes – that is, numerical integration of a given enzyme reaction rate timeseries curve – the calculated static enzyme reaction rates at one sampling point do not affect those calculated at the next sampling point. In contrast, in the HDS method – that is, the initial value problem of simultaneous differential equations – the calculated static enzyme reaction rates at one integration step affect the calculation in the next step. Accordingly, the error calculated in a process for distinguishing dynamic/static is not always equal to the error in the hybrid model. Thus, comparison of errors between these two types of calculations is required.
Figure 3 shows the relationship between the MRE of metabolite concentrations obtained by processes for distinguishing dynamic/static enzymes and the MRE of metabolite concentrations in the hybrid models for various weighting coefficients. The errors obtained by these two methods showed a high positive correlation (r = 0.948). This result indicates that the accuracy of the hybrid model constructed using the estimated distinguishing of dynamic/static enzymes exactly reflects the magnitude of the error estimated by processes for distinguishing dynamic/static enzymes. Therefore, the proposed method for distinguishing dynamic/static modules can be used to build a hybrid model.
The error in the hybrid models was higher than that obtained by processes for distinguishing dynamic/static enzymes. In particular, in the distinguishing of dynamic/static enzymes of S. cerevisiae with w = 0.250, a considerable degree of error enlargement was shown in the hybrid model. This result can be considered to have been caused by error propagation at each integration step, as expected.The relationship between w and the MRE of the metabolite concentration time series and that between w and the static enzyme ratio was examined (Figure 4). The two metabolic systems tested showed very similar results, perhaps because both models deal with centralcarbon metabolism. The dependency of the MRE and the static enzyme ratio on w showed a staircase pattern, rather than a pattern of simple linear increase (or decrease).
Evaluation of the total process using noiseadded data
In the previous sections, we used noisefree values to obtain a clear evaluation of the proposed method itself. However, real experimental data of metabolite concentrations are generally noisy. For practical use of the proposed method, the effect of noise on the process for distinguishing dynamic/static enzymes should be evaluated. Thus, we added noise to the noisefree data and then smoothed the noisy data for use in distinguishing the dynamic/static enzymes. In this study, simple smoothing by fitting to a polynomial or rational function of time was employed. The smoothing functions that were used and their parameters are shown in Supplementary Tables S2 and S3 (see additional file 1). Comparisons of noisefree values, noiseadded values, and smoothed curves of metabolites are shown in Supplementary Figure S3 (see additional file 1). The results of distinguishing dynamic/static enzymes from the noisy metabolite concentration data are shown in Table 2. In most cases, when noiseadded data were used, entirely or almost the same distinctions between dynamic/static enzymes were obtained as when noisefree data were used. However, in the results for S. cerevisiae obtained using smoothed noisy data, when w < 0.250, the number of static enzymes tended to be larger than in the results obtained using noisefree data. In the results for E. coli, the same tendency was observed when w = 0.010. Because the smoothing process of the metabolite concentration time series might result in loss of the highfrequency component of the time series data, the smoothed data might apparently change more slowly than is actually the case. Thus, when smoothed noisy data are used, the number of required dynamic enzymes in a HDS model tends to be smaller than the number needed when noisefree data are used. Because more precise metabolite concentrations need to be calculated when w is small, this tendency might be enhanced.
Discussion
Estimation of enzyme reaction rates
As shown in Table 1, the accuracy of the estimations of the enzyme reaction rates was confirmed by the reproduced metabolite concentrations, except for PGP in E. coli. Since the concentration of PGP was very low (average concentration, 3.60 × 10^{3} mM), even a slight error in the enzyme reaction rate had a large influence. In fact, the average errors between the true enzyme reaction rate time series and the estimated enzyme reaction rate time series for both GAPDH (PGPproducing enzyme) and PGK (PGPconsuming enzyme) in E. coli were adequately small, 2.44% and 1.46%, respectively. In the process for distinguishing dynamic/static enzymes, the average of the squared errors of all metabolite concentrations is used to calculate the fitness function (Eq. (5)); thus, an error in only one metabolite concentration has a limited effect. Actually, the results of distinguishing dynamic/static enzymes without the PGP time series (data not shown) were entirely the same as those shown in Table 2. However, if many metabolites with low concentrations are included in the modelled metabolic system, the processes for distinguishing dynamic/static enzymes may cause an erroneous conclusion to be drawn. This is a limitation of the current procedure. In comparison with the results for E. coli, errors for all metabolites for S. cerevisiae were adequately small, because the dynamics of the metabolic system in S. cerevisiae is relatively slow compared with the sampling frequency.
Another difficulty in applying the proposed method is that we assume that the concentrations of all metabolites are measurable. It is expected that highthroughput measurement techniques for detecting a huge number of metabolites, such as capillary electrophoresis combined with mass spectrometry (CEMS) [17–19], can be used for such comprehensive measurements. The 1s sampling interval employed in this study is feasible, because some rapidsampling instruments capable of drawing multiple samples within 1 s from a bioreactor have already been developed [26–28].
Distinction of dynamic and static enzymes
After a process for distinguishing dynamic/static enzymes is completed, the MRE in the corresponding hybrid model can be estimated using the linear relationship between the MRE in the process for distinguishing dynamic/static enzymes and the MRE in the hybrid model (Figure 3). This information helps to build a hybrid model that has the desired accuracy.
The staircase pattern of the relationships between the error and static enzyme ratio with decreasing w, observed in Figure 4, was probably caused by a property of metabolic systems. In a testing system, the number of enzymes that can potentially be allocated to the static module may be restricted. If w is greatly changed, the few potentially static enzymes would eventually start to be converted to static enzymes.
Weighting coefficient in the fitness function
The weighting coefficient in the fitness function (Eq. (5)) is a tuning parameter. Since a suitable value for the weighting coefficient (w) is not given a priori, we need to consider how to define the value.
As shown in Figure 4, with a w of 1.000, about half of the enzymes were discriminated to the static module. Thus, a large amount of experimental work can be saved because no kinetic information is required by the static module. The MRE at w = 1.000 was 15.2% for the E. coli hybrid model and 18.6% for the S. cerevisiae hybrid model (Figure 4). These errors are acceptable considering the accuracy of the experimentally measured metabolite concentrations. Thus, w = 1.000 may simply be chosen at the initial trial stage of model construction. When a more precise model is required, a smaller w can be used. Even if w is set to between 0.025 and 0.100, the proportion of static enzymes remains at about 30% for both the metabolic systems tested. Our recommendation for w for general modelling is 0.050. At around this w value, the sensitivity of the error to a change of w is low; thus, strict specification of w is not required. Moreover, even if the actual error in the constructed hybrid model becomes considerably higher than the expected value – as in the case of S. cerevisiae at w = 0.250 – the actual error remains low.
Noise in metabolome data
As shown in Table 2, almost the same results in distinguishing dynamic/static enzymes were obtained between the procedures using noisefree data and those using noiseadded data. This result could be predicted because most metabolite time series were successfully reproduced from the noisy data by the smoothing treatment, as shown in Figure S3. This result indicates that the proposed method for distinguishing dynamic/static enzymes can be applied to noisy measurements if a suitable noise reduction method is employed. To remove noise and obtain the slopes of metabolite concentration time series, a smoothing technique based on an artificial neural network, proposed by Voit et al. [29–31], is efficient. Many other noise cancellation techniques have been proposed for biochemical time series data [32–35]. For example, Rizzi et al. [36] obtained timecourse functions of metabolites from noisy metabolite concentration measurements and used those functions to tune the parameters in their dynamic model.
Toward construction of accurate hybrid models
In the HDS method, accurate kinetics should be known not only for system boundary enzymes but also for all enzymes assigned to the dynamic modules. For this reason, highthroughput techniques for determining accurate and detailed enzyme kinetics are needed for the efficient development of models of metabolic systems. A promising powerlaw approach, generalized mass action (GMA) [37, 38], may be used to solve this problem. This method has a large representational space that enables enzyme kinetics to be sufficiently expressed in spite of its simple fixed form. Although modelling that uses this kind of powerlaw approach from time series data is often difficult owing to their nonlinear properties, Polisetty et al. [39] have proposed a method employing branchandbound principles to find optimized parameters in GMA models. Using this method, the global optimal parameter set can be efficiently searched.
To ensure the validity of the predicting performance of an HDS model, careful perturbation experiments should be carried out to obtain the metabolome time series data to be used for distinguishing dynamic/static enzymes. The metabolite concentration variations used should be those considered to be of the maximum possible magnitude under the modelled conditions. To reproduce a rapidly changing metabolite concentration time series by an HDS model, a larger number of dynamic enzymes is required. Thus, if the number of dynamic enzymes included in the model is defined by using data showing the maximum possible variation in magnitude, that is, the model is constructed with the maximum possible number of dynamic enzymes, then the model can calculate all probable states of the system. For instance, consider building a metabolic model of cultured cells in a reactor, where the model has no mechanism for calculating gene expression levels or the consequent changes in protein concentrations (most proposed metabolic models are of this type). A substratepulse injection experiment giving the maximal substrate concentration that does not cause changes in gene expression levels in the cells (i.e., enzyme concentrations in the cells are kept constant) is useful for distinguishing dynamic/static enzymes. To determine the maximal permitted substrate concentration, many preliminary experiments may be required, and this seems to decrease the value of the HDS method, which aims to reduce experimental efforts. However, fundamentally speaking, such evaluation of the limits of a model's parameters is absolutely necessary for maintaining the accuracy of calculations in any kind of modelling, not only in HDS modelling. Therefore, this requirement for experiments to determine the maximal possible variation is not a specific disadvantage of the HDS method.
Conclusion
The proposed method of using metabolite concentration time series,i.e., experimentally measurable variables, enables us to discriminate dynamic/static enzymes to construct a hybrid model. In this method, the enzyme reaction rate time series are estimated from metabolite concentration time series data. Since this estimation relies on only the mass balance in the system, no kinetic information about internal enzymes is required. Therefore, the aim of employing the HDS method – to reduce the experimental effort required to obtain enzyme kinetics information – can be achieved. Two microbial centralcarbon metabolism models were used to evaluate our method. Centralcarbon metabolism has many feedback loops and is rigidly controlled to maintain homeostasis of a living cell. Since our method was successfully applied for such a strictly regulated system, we believe it will have wideranging applicability to many types of metabolic systems. Furthermore, the analysis using noisy metabolite concentration data demonstrated that, for the most part, the proposed method tolerates noise well.
Abbreviations
 M:

Metabolites
 2PG:

2phosphoglycerate
 3PG:

3phosphoglycerate
 6PG:

6phosphogluconate
 ACA acetaldehyde:

intracellular
 ACA_{x} acetaldehyde:

extracellular
 CN_{o} cyanide:

mixed flow
 CN_{x} cyanide:

extracellular
 DHAP:

dihydroxyacetone phosphate
 E4P:

erythrose 4phosphate
 EtOH ethanol:

intracellular
 EtOH_{x} ethanol:

extracellular
 F6P:

fructose 6phosphate
 FDP fructose 1:

6bisphosphate
 G1P:

glucose 1phosphate
 G6P:

glucose 6phosphate
 GAP:

glyceraldehyde 3phosphate
 Glc_{o} glucose:

mixed flow
 Glc_{x} glucose:

extracellular
 Glyc glycerol:

intracellular
 Glyc_{x} glycerol:

extracellular
 PEP:

phosphoenolpyruvate
 PGP 1:

3bisphosphoglycerate
 Pyr:

pyruvate
 Rib5P:

ribose 5phosphate
 Ribu5P:

ribulose 5phosphate
 Sed7P:

sedoheptulose 7phosphate
 Xyl5P:

xylulose 5phosphate
 E:

Enzymes/reactions
 ADH:

acetaldehyde dehydrogenase
 AK:

adenylate kinase
 ALDO:

aldolase
 consum:

ATP consumption
 ifACA:

ddiffusion of acetaldehyde
 difEtOH:

diffusion of EtOH
 difGlyc:

diffusion of glycerol
 ENO:

enolase
 G6PDH:

glucose6phosphate dehydrogenase
 GAPDH:

glyceraldehyde3phosphate dehydrogenase
 GlcTrans:

glucose transporter
 HK:

hexokinase
 lpGlyc:

lumped glycerol formation reaction
 lpPEP:

lumped PEP formation reaction
 PDC:

pyruvate decarboxylase
 PFK:

phosphofructokinase
 PGDH:

6phosphogluconate dehydrogenase
 PGI:

glucose6phosphate isomerase
 PGK:

phosphoglycerate kinase
 PGluMu:

phosphoglycerate mutase
 PGM:

phosphoglucomutase
 PK:

pyruvate kinase
 R5PI:

ribosephosphate isomerase
 Ru5P:

ribulosephosphate epimerase
 TA:

transaldolase
 TIS:

triosephosphate isomerase
 TKa transketolase:

reaction a
 TKb transketolase:

reaction b
References
 1.
van der Werf MJ, Jellema RH, Hankemeier T: Microbial metabolomics: replacing trialanderror by the unbiased selection and ranking of targets. J Ind Microbiol Biotechnol. 2005, 32 (6): 234252. 10.1007/s1029500502314.
 2.
Ishii N, Soga T, Nishioka T, Tomita M: Metabolome analysis and metabolic simulation. Metabolomics. 2005, 1 (1): 2937. 10.1007/s1130600511055.
 3.
Patton WF: Detection technologies in proteome analysis. J Chromatogr B Analyt Technol Biomed Life Sci. 2002, 771 (12): 331. 10.1016/S15700232(02)000430.
 4.
Bertone P, Snyder M: Advances in functional protein microarray technology. Febs J. 2005, 272 (21): 54005411. 10.1111/j.17424658.2005.04970.x.
 5.
Leung YF, Cavalieri D: Fundamentals of cDNA microarray data analysis. Trends Genet. 2003, 19 (11): 649659. 10.1016/j.tig.2003.09.015.
 6.
Johnson JM, Edwards S, Shoemaker D, Schadt EE: Dark matter in the genome: evidence of widespread transcription detected by microarray tiling experiments. Trends Genet. 2005, 21 (2): 93102. 10.1016/j.tig.2004.12.009.
 7.
Schomburg I, Chang A, Schomburg D: BRENDA, enzyme data and metabolic information. Nucleic Acids Res. 2002, 30 (1): 4749. 10.1093/nar/30.1.47.
 8.
Sundararaj S, Guo A, HabibiNazhad B, Rouani M, Stothard P, Ellison M, Wishart DS: The CyberCell Database (CCDB): a comprehensive, selfupdating, relational database to coordinate and facilitate in silico modeling of Escherichia coli. Nucleic Acids Res. 2004, 32 (Database issue): D2935. 10.1093/nar/gkh108.
 9.
Selkov E, Basmanova S, Gaasterland T, Goryanin I, Gretchkin Y, Maltsev N, Nenashev V, Overbeek R, Panyushkina E, Pronevitch L, Selkov E, Yunus I: The metabolic pathway collection from EMP: the enzymes and metabolic pathways database. Nucleic Acids Res. 1996, 24 (1): 2628. 10.1093/nar/24.1.26.
 10.
Yugi K, Nakayama Y, Kinoshita A, Tomita M: Hybrid dynamic/static method for largescale simulation of metabolism. Theor Biol Med Model. 2005, 2 (1): 4210.1186/17424682242.
 11.
Shimizu K: Metabolic flux analysis based on 13Clabeling experiments and integration of the information with gene and protein expression patterns. Adv Biochem Eng Biotechnol. 2004, 91: 149.
 12.
Aiba S, Matsuoka M: Identification of metabolic model: Citrate production from glucose by Candida lipolytica. Biotechnol Bioeng. 1979, 21 (8): 13731386. 10.1002/bit.260210806.
 13.
Okino MS, Mavrovouniotis ML: Simplification of mathematical models of chemical reaction systems. Chem Rev. 1998, 98 (2): 391408. 10.1021/cr950223l.
 14.
Gerdtzen ZP, Daoutidis P, Hu WS: Nonlinear reduction for kinetic models of metabolic reaction networks. Metab Eng. 2004, 6 (2): 140154. 10.1016/j.ymben.2003.11.003.
 15.
Moore EH: On the reciprocals of the general algebraic matrix. Bull American Math Soc. 1920, 26: 394395.
 16.
Penrose R: A generalized inverse for matrices. Proc Cambridge Philos Soc. 1955, 51 (406413):
 17.
Fiehn O, Kopka J, Dormann P, Altmann T, Trethewey RN, Willmitzer L: Metabolite profiling for plant functional genomics. Nat Biotechnol. 2000, 18 (11): 11571161. 10.1038/81137.
 18.
Soga T, Ueno Y, Naraoka H, Matsuda K, Tomita M, Nishioka T: Pressureassisted capillary electrophoresis electrospray ionization mass spectrometry for analysis of multivalent anions. Anal Chem. 2002, 74 (24): 62246229. 10.1021/ac0202684.
 19.
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): 488494. 10.1021/pr034020m.
 20.
Chassagnole C, NoisommitRizzi N, Schmid JW, Mauch K, Reuss M: Dynamic modeling of the central carbon metabolism of Escherichia coli. Biotechnol Bioeng. 2002, 79 (1): 5373. 10.1002/bit.10288.
 21.
Hynne F, Dano S, Sorensen PG: Fullscale model of glycolysis in Saccharomyces cerevisiae. Biophys Chem. 2001, 94 (12): 121163. 10.1016/S03014622(01)002290.
 22.
Delgado JP, Liao JC: Identifying ratecontrolling enzymes in metabolic pathways without kinetic parameters. Biotechnol Prog. 1991, 7 (1): 1520. 10.1021/bp00007a003.
 23.
Kitayama T, Kinoshita A, Sugimoto M, Nakayama Y, Tomita M: A simplified method for powerlaw modelling of metabolic pathways from timecourse data and steadystate flux profiles. Theor Biol Med Model. 2006, 3 (1): 2410.1186/17424682324.
 24.
Willett P: Genetic algorithms in molecular recognition and design. Trends Biotechnol. 1995, 13 (12): 516521. 10.1016/S01677799(00)890150.
 25.
Shampine LF, Reichelt MW: The MATLAB ODE Suite. SIAM J Scientific Comp. 1997, 18: 122. 10.1137/S1064827594276424.
 26.
Schaefer U, Boos W, Takors R, WeusterBotz D: Automated sampling device for monitoring intracellular metabolite dynamics. Anal Biochem. 1999, 270 (1): 8896. 10.1006/abio.1999.4048.
 27.
Buziol S, Bashir I, Baumeister A, Claassen W, NoisommitRizzi N, Mailinger W, Reuss M: New bioreactorcoupled rapid stoppedflow sampling technique for measurements of metabolite dynamics on a subsecond time scale. Biotechnol Bioeng. 2002, 80 (6): 632636. 10.1002/bit.10427.
 28.
Hoque MA, Ushiyama H, Tomita M, Shimizu K: Dynamic responses of the intracellular metabolite concentrations of the wild type and pykA mutant Escherichia coli against pulse addition of glucose or NH3 under those limiting continuous cultures. Biochem Eng J. 2005, 26 (1): 3849. 10.1016/j.bej.2005.05.012.
 29.
Voit EO, Almeida J: Decoupling dynamical systems for pathway identification from metabolic profiles. Bioinformatics. 2004, 20 (11): 16701681. 10.1093/bioinformatics/bth140.
 30.
Marino S, Voit EO: An automated procedure for the extraction of metabolic network information from time series data. J Bioinform Comput Biol. 2006, 4 (3): 665691. 10.1142/S0219720006002259.
 31.
Chou IC, Martens H, Voit EO: Parameter estimation in biochemical systems models with alternating regression. Theor Biol Med Model. 2006, 3: 2510.1186/17424682325.
 32.
Wood GA: Data smoothing and differentiation procedures in biomechanics. Exerc Sport Sci Rev. 1982, 10: 308362. 10.1249/0000367719820100000010.
 33.
Gazzani F: Comparative assessment of some algorithms for differentiating noisy biomechanical data. Int J Biomed Comput. 1994, 37 (1): 5776. 10.1016/00207101(94)90071X.
 34.
Eilers PH: A perfect smoother. Anal Chem. 2003, 75 (14): 36313636. 10.1021/ac034173t.
 35.
Tsai KY, Wang FS: Evolutionary optimization with data collocation for reverse engineering of biological networks. Bioinformatics. 2005, 21 (7): 11801188. 10.1093/bioinformatics/bti099.
 36.
Rizzi M, Baltes M, Theobald U, Reuss M: In vivo analysis of metabolic dynamics in Saccharomyces cerevisiae. II. Mathematical model. Biotechnol Bioeng. 1997, 55 (4): 592608. 10.1002/(SICI)10970290(19970820)55:4<592::AIDBIT2>3.0.CO;2C.
 37.
Voit EO, Ferreira AEN: Computational analysis of biochemical systems: a practical guide for biochemists and molecular biologists. 2000, Cambridge, NY , Cambridge University Press
 38.
Torres NV, Voit EO: Pathway analysis and optimization in metabolic engineering. 2002, Cambridge, NY , Cambridge University Press
 39.
Polisetty PK, Voit EO, Gatzke EP: Identification of metabolic system parameters using global optimization methods. Theor Biol Med Model. 2006, 3: 410.1186/1742468234.
Acknowledgements
The authors would like to thank Katsuyuki Yugi, Ayako Kinoshita and Yoshihiro Toya for insightful discussions. This work was supported in part by 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 grantinaid from the Ministry of Education, Culture, Sports, Science and Technology for the 21st Century Center of Excellence (COE) Program (Understanding and Control of Life's Function via Systems Biology).
Author information
Additional information
Competing interests
The author(s) declare that they have no competing interests.
Authors' contributions
NI contributed to the development of the proposed method and wrote this manuscript. YN provided the basic ideas and directed the project, and MT was the project leader. All authors read and approved the final manuscript.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
About this article
Received
Accepted
Published
DOI
Keywords
 Hybrid Model
 Metabolite Concentration
 Metabolic System
 Metabolic Flux Analysis
 Mean Relative Error