Theoretical Biology and Medical Modelling

Background: In order to improve understanding of metabolic systems there have been attempts to construct S-system models from time courses. Conventionally, non-linear curve-fitting algorithms have been used for modelling, because of the non-linear properties of parameter estimation from time series. However, the huge iterative calculations required have hindered the development of large-scale metabolic pathway models. To solve this problem we propose a novel method involving power-law modelling of metabolic pathways from the Jacobian of the targeted system and the steady-state flux profiles by linearization of S-systems. Results: The results of two case studies modelling a straight and a branched pathway, respectively, showed that our method reduced the number of unknown parameters needing to be estimated. The time-courses simulated by conventional kinetic models and those described by our method behaved similarly under a wide range of perturbations of metabolite concentrations. Conclusion: The proposed method reduces calculation complexity and facilitates the construction of large-scale S-system models of metabolic pathways, realizing a practical application of reverse engineering of dynamic simulation models from the Jacobian of the targeted system and steady-state flux profiles.


Background
Systematic modelling has emerged as a powerful tool for understanding the mathematical properties of metabolic systems. The rapid development of metabolic measurement techniques has driven advances in modelling, especially using data on the effects of perturbations of metabolite concentrations, which contain valuable information about metabolic pathway structure and regulation [1]. A power-law approximation for representing enzyme-catalyzed reactions, known as Biochemical Systems Theory, is an effective approach for understanding metabolic systems [2,3]. Generalized Mass Action (GMA) and S-systems [4,5], which are often used as power-law modelling approaches, have wide representational spaces that permit adequate expression of enzyme kinetics [6] in spite of their simple fixed forms. Moreover since S-system forms have a smaller number of parameters than GMA forms, the S-system is an appropriate modelling framework. The derivation of an S-system model from given experimental data is a powerful tool not only for understanding nonlinear properties but also for determining the regulatory structure of the system [7,8].
S-system modelling from time-course data is often difficult due to its non-linear properties. Non-linear-fitting algorithms, such as genetic algorithms or artificial neural networks, have been used to resolve this problem [9][10][11][12][13][14]. Although these methods can be applied to metabolic pathways, massive computing power is required in the case of targeted models involving a number of closelyconnected, underspecified parameters [13]. Moreover, the wider the range of targeted metabolic pathways, the more likely is the occurrence of local minima due to the expansion of the parameter search space. The network-structures-segmentation method can reduce the total parameter search range in genetic network modelling [15] but it is difficult to apply this to metabolic pathways because of the close relationships between reversible reactions and because of allosteric regulation. Diaz-Sierra and Fairén have proposed an approach, based on the steadystate assumption, that allows the construction of S-system models from a Jacobian matrix of the system [16]. Since the Jacobian constrains the search range of underspecified parameters at the optimization stage, these authors' method allows efficient parameter estimation. However, the problem of an excess number of parameters requiring estimation remains unsolved.
We present an approach to power-law modelling of metabolic pathways from the Jacobian of the targeted system and steady-state flux profiles with linearization of the Ssystem. This reduces the number of underspecified parameters. Two numerical experiments show that the S-system model generated by this method describes similar dynamic behaviour to that indicated by conventional kinetic models.

Retrieving the Jacobian from time-course data
As a first step, the Jacobian must be obtained from metabolic time-course data. In this section, we summarize the method of Sorribas et al [17], which we use in this work.
In biochemical systems, the Jacobian can be defined as: where J is the Jacobian matrix, and δX represents a small perturbation and contains the concentration X i as its elements. The elements of the Jacobian can be obtained from perturbed time-courses using linear least-squares fitting [17,18]. This method is based on the fact that transients yield linear responses to small perturbations under steady-state conditions [5]. The mathematical basis for this is that the linear representation constitutes the firstorder term of a Taylor series expansion, which is sufficiently accurate in this situation.

Determination of the kinetic orders of the S-system
The S-system is a power-law representation constructed of two terms: the production rate and the degradation rate: where α i and β i are rate constants, g ij and h ij are kinetic orders, and x i represents the concentration of a compound.
In steady state conditions, it can be expressed simply as: where Flux is the sum of the steady-state fluxes into x i , and V + and Vare production and consumption terms, respectively.
In steady state conditions, the production term and the consumption term in Eq.(2) can therefore be represented as: where x j,0 represents the steady-state concentration of x j .
Consequently α i and β i yield: Defining X i as: = y i (x 1 , ..., x n ) (i = 1, ..., n) (7) the Jacobian of Eq. (2) can be represented as: Once Jacobian J ij , x j,0 , and Flux are given, Eq. (10) is a constraint for the determination of the kinetic orders g ij and h ij . Savageau has described the linearization of S-system as an "F-factor" for stability analysis [19]. We use this representation to estimate parameter values.
In most cases g ij and/or h ij are available from the structure of the metabolic pathway; however, in the absence of known kinetic orders, parameter estimates are needed to determine them. In such cases Eq. (10) is adopted as limiting the parameter search range.

Case study 1: a linear biochemical pathway
In the first case study, we applied this approach to the published biochemical model of yeast galactose metabolism shown in Figure 1 [20]. This model consists of five metabolites, and four enzyme reactions mainly described by Michaelis-Menten equations. The kinetic equations, systems parameters, and initial conditions are listed in Appendix 1.
In this case, all parameters of the S-system model were determined, without the need for estimation. The biochemical model could be converted into the following Ssystem form.
Since X 1 is an independent variable, it was omitted from the listed rate equations.
In this modelling process, all the parameters were determined by simple calculations of this kind. The resulting Ssystem model was:

Case study 2: A branched biochemical pathway
The branched biochemical pathway shown in Figure 2 was tested in the second example. This pathway has the typical features of a biochemical pathway: branching, feedback regulation, and product inhibition (see Appendix 2). X 3 is the inhibitor of both V 3 and V 4 .
Most of the parameters were obtained by our method. However one parameter could not be determined by calculation and a parameter estimation method was used. The model could be converted into the following S-system form: The following constraint was provided by the pathway structure: h 11 = g 21 . Steady-state concentrations were: X 1,0 = 0.067, X 2,0 = 0.049, X 3,0 = 0.081, X 4,0 = 0.041, and steadystate fluxes were: V 1 = V 2 = 0.1, V 3 = V 5 = 0.043, and V 4 = V 6 = 0.057. The Jacobian was: More than half of the parameters were obtained by simple calculations using the steady-state concentrations, the steady-state fluxes, and the Jacobian. The following four parameters could not be determined: α 3 , g 33 , β 3 , and h 33 . For example, h 11 is a parameter which could be determined by: Because g 33 is a parameter that could not be determined by a calculation, an estimate provided by some constraint was needed. As α 3 , h 33 , and β 3 could be treated as dependent parameters, once g 33 was obtained all four parameters were determined: Moreover, the search range of g 33 could be limited by the definition that h 33 has a positive value since X 3 is the substrate of the reaction V 3 . The search range of g 33 was set from 0 to -0.95. When the number of underspecified parameters was one whose search range was limited, a lin- .
.  The branched pathway Figure 2 The branched pathway. X 3 is an inhibitor of both V 3 andV 4 .
ear-fitting algorithm could be used to fit the time-course data in the original kinetic model. In this modelling, g 33 was determined by a linear optimization method. Consequently, the following S-system model was generated: The list of parameters used in the above calculation is shown in Table 1.

Comparing transient dynamic responses after perturbation
To evaluate the versatility of this method, the transient dynamics of the S-system model in response to various perturbations were compared with those of the original Michaelis-Menten model, by calculating the mean relative errors (MRE): where n is the number of metabolites in the model, m, the number of sampling points in the time course, x' i (t) represents the time course calculated from the created S-sys-tem model, and x i (t), the time course calculated from the original Michaelis-Menten model. In this experiment there were 10 sampling points, the interval between sampling points was 0.5, and X 2 in case study 1 and X 1 in case study 2 were the targets of perturbations ranging from 0% to 200%. The time-courses of metabolites in response to perturbations of 100% are shown in Figure 3.
In both examples, similar dynamic behaviour was observed in the S-system model and the reference model as a response to the perturbation.
The changes of MREs in response to the perturbation range are shown in Figure 4. The initial MREs with no perturbation were within 1%. Although a slight increase in MREs was observed in both case studies, they remained within 4% at a perturbation range of 200%.

Discussion
Power-law modelling from time-course data of metabolite concentrations often requires parameter estimates that depend on the size of the target metabolic network. Especially when developing a large-scale metabolic model, this requires problem-specific simplifications. Our method can reduce the number of underspecified parameters by using steady-state flux profiles and the Jacobian of the targeted system derived from time courses of metabolites, and is thus suitable for large-scale power-law modelling.
To validate our methods we used two existing biochemical pathway models, the straight and branched pathway models, described by dynamic equations. In the S-system, modelling of a linear metabolic pathway with 12 parameters (case study 1), our method determines all parameters accurately, whereas the method developed by Diaz-Sierra and Fairén [16] leaves at least four parameters underspecified that include error correction parameters. In the case of the branched metabolic pathway with 16 parameters (case study 2), our method determines all but one parameter, whereas the method of Diaz-Sierra and Fairén has 12 underspecified parameters that include error correction parameters. The case studies demonstrate that our method of developing S-system models from time series can reduce the number of underspecified parameters more efficiently than the previously reported method [16]. Furthermore, the perturbation response experiments show that the models created by our method can reproduce dynamics similar to the reference models, since the MRE was around 3% when the perturbation range was 200% (Fig. 4). S-system models generated by our method can provide accurate simulations within a wide range of the steady-state point. This limitation does not prevent the modelling and analysis of metabolic pathways, as it does with many power-law metabolic models [4,5].
Robustness against experimental noise is an important requirement for the practical application of modelling dX dt X  methods. Since the Jacobian is rather sensitive to the noise in time series [18], a robust corrective response to the numerical errors contained in a Jacobian is important in developing a model involving its use. In the method of Diaz-Sierra and Fairén [16], error correction parameters are incorporated into the S-system model to reduce the effect of experimental noise. However, this approach may lead to an increase in the number of underspecified parameters. As an alternative approach to reducing the experimental noise in time-course data, the estimation of appropriate slopes using a non-linear neural network model is effective [10], and can provide an error-controlled time course that enables the Jacobian to be obtained with high accuracy. Methods for obtaining accurate esti-mates of the specific effects of general types of perturbation have been discussed [21] and might enable more precise analysis of data such as time-scale metabolite concentrations obtained under perturbations.
To assess the robustness of our method against experimental noise, we measured the calculation error when numerical errors were manually inserted into the Jacobian or steady state fluxes. Table 2 summarizes the MREs of the simulated time trajectories of the S-system and the original Michaelis-Menten model. It is evident that the Jacobian is quite sensitive to numerical error because of the direct effect on a kinetic order of a Jacobian. Furthermore, the steady-state flux profile data may include experimen-metabolite changes in response to perturbation of one metabolite by 100% Figure 3 metabolite changes in response to perturbation of one metabolite by 100%. (a) and (b) represent the changes in metabolites in case study 1 when X 2 is the perturbation target, derived using the reference and S-system model, respectively. (c) and (d) are the changes in case study 2 when X 1 is the perturbation target given by the reference and S-system model, respectively.
tal noise. However, the models created by our method are relatively robust to errors in the Jacobian and steady-state fluxes, indicating that these methods will be useful in practice.
Sorribas and Cascante proposed a method for identifying regulatory patterns using a given set of logarithmic gain measurements [7]. In their paper, they suggested that one strategy for selecting possible patterns is to perform per-turbation experiments and to measure the corresponding dynamic response. For practical application of this approach, it is crucial to develop appropriate ways of performing the required experiments; however, measuring the logarithmic gain resulting from various steady-state fluxes is not practical due to a lack of exhaustive measurement method. Our approach can be used to develop reasonably accurate models of metabolic pathways by using a single set of appropriate steady-state flux profiles. Although steady-state flux profile data remains difficult to be measured directly and comprehensively, several methods of measuring steady-state flux profiles by using isotopes have been developed [22,23]. Our method assumes that the Jacobian obtained is accurate. Therefore it is important to obtain time-course data reflecting transient dynamics after a suitably small perturbation in which the Jacobian behaves in a linear manner [17].
Comprehensive metabolome data will undoubtedly accumulate as a consequence of advances in metabolic measurement techniques [24][25][26]. In our laboratory we have developed a high-throughput technique using capillaryelectrophoresis mass spectrometry that provides effective time-course data involving a few hundred ionic metabolites [27,28]. Our method promises to provide highthroughput modelling of large-scale metabolic pathways by exploiting the accumulating metabolome and steadystate flux profile data along with the anticipated developments in metabolome measurement techniques.

Conclusion
Our method provides stable and high-throughput S-system modelling of metabolic pathways because it drastically reduces underspecified parameters by employing the  steady-state flux profile and Jacobian retrieved from timecourse data. S-system models generated by this method can provide accurate simulations within a wide range around the steady-state point. In combination with the metabolome measurement techniques it should permit high-throughput modelling of large-scale metabolic pathways.