 Methodology article
 Open Access
 Published:
New timescale criteria for model simplification of bioreaction systems
BMC Bioinformatics volume 9, Article number: 338 (2008)
Abstract
Background
Quasisteady state approximation (QSSA) based on timescale analysis is known to be an effective method for simplifying metabolic reaction system, but the conventional analysis becomes timeconsuming and tedious when the system is large. Although there are automatic methods, they are based on eigenvalue calculations of the Jacobian matrix and on linear transformations, which have a high computation cost. A more efficient estimation approach is necessary for complex systems.
Results
This work derived new timescale factor by focusing on the problem structure. By mathematically reasoning the balancing behavior of fast species, new timescale criteria were derived with a simple expression that uses the Jacobian matrix directly. The algorithm requires no linear transformation or decomposition of the Jacobian matrix, which has been an essential part for previous automatic timescaling methods. Furthermore, the proposed scale factor is estimated locally. Therefore, an iterative procedure was also developed to find the possible multiple boundary layers and to derive an appropriate reduced model.
Conclusion
By successive calculation of the newly derived timescale criteria, it was possible to detect multiple boundary layers of full ordinary differential equation (ODE) models. Besides, the iterative procedure could derive the appropriate reduced differential algebraic equation (DAE) model with consistent initial values, which was tested with simple examples and a practical example.
Background
The dynamic simulation of bioreaction pathways is becoming more important as the kinetic information of various pathways is revealed. Moreover, the necessary data for the specific pathways are easily obtained through various channels, including the internet. However, there are fundamental difficulties in the numerical solution of the differential equation system: model stiffness is one. a bioreaction system consists of multiple reactions with various enzymes that have different turnover numbers, meaning various magnitudes of reaction rates. Furthermore, the metabolites involved in one reaction can also participate in other reactions in the same system. These characteristics make typical biosystems strongly coupled and have multiple timescales. Therefore, ordinary differential equations (ODEs) based on the dynamic modeling of a metabolic system are usually stiff.
The stiffness problem requires unnecessary effort to track the boundary layer solutions, hence, the computational efficiency decreases. Furthermore, when the simulation is concurrent with the experiment, the calculation accuracy is closely related to the measurement interval. If the measurement interval is modest so that it is impossible to find some specific parameters, then the numerical result does not need to be precise. In these situations, by sacrificing the accuracy modestly, simplifying the model structure is necessary, which is often the case for parameter estimations and sensitivity analyses [1].
For this purpose, apart from a numerical analysis approach, the kinetic field's specific solution methods have been required [2]. Traditionally, simplification of an original complex model, such as a quasisteady state approximation (QSSA) and a partial equilibrium approximation, have been applied to relieve the stiffness characteristics [3–7]. However, since these approaches require the practitioner's intuition and experience, computational methods have been developed, especially vividly in combustion engineering fields. There are two important procedures in computational model simplification: the determination of the simplification criteria and the determination of the slow invariant manifolds.
Some computational methods concentrated on deriving the correct slow manifolds. They are iterative trajectory based methods [8, 9], the method of invariant manifold (MIM) [10], the minimal entropy production trajectory (MEPT) based methods [11, 12], and the nonlinear model reduction method [13]. Gorban et al. collected and reviewed such kinds of methods [14].
The others suggested the appropriate simplification criteria: generalized sensitivity analyses [15] used singular values of the sensitivity matrix as the scale factor, while computational singular perturbation (CSP) [16–18], intrinsic low dimensional manifold (ILDM) [19], and dynamic dimension reduction, which is a modified version of ILDM [20, 21], used the eigenvalue analysis of the system's Jacobian matrix with focusing on the timescales of the system. Currently, ILDM based methods have been applied to reduce complexity of biochemical systems [21–23]. Since CSP or ILDM based methods use the dynamic properties of a system, they can give the dynamically useful information of the system.
Even if there are several differences between the CSP and ILDM, the most important ideas are similar: the determination of the speed ranking is based on the eigenvalue calculation of the Jacobian matrix, which require at least O(n^{3}) flops of computation, and the derivation of solution is based on a linear transformation of the original system. Furthermore, they share common barriers in producing an explicitly reduced model, generally as a result of the system transformation. This feature is important when model simplification approaches are related to not only computational efficiency but also parameter estimation.
This work suggests an automatic method for speed ranking that directly uses the Jacobian matrix without a system transformation. Due to its simplicity, this approach requires O(n) flops of computation, being more physically intuitive. In addition to the scaling procedure, the decision of the differential and algebraic variables in a slow dynamic regime after relaxation of the fast dynamics is also introduced. The result of the proposed process is an explicitly phrased model, which can be a route to distinguish the meaningful parameters from the unobservable ones.
Methods
Criteria for balancing
In homogeneous chemical kinetics, the dynamic model can be written in the following form of the ODE.
Since a chemical reaction system generally consists of production and loss terms, the ODE can be rewritten as:
or with a matrixvector notation,
where y ∈ ℝ^{+n}is a concentration vector, P_{ i }: ℝ^{+n}→ ℝ^{+} is a production term, and L_{ i }: ℝ^{+n}→ ℝ^{+} is a loss terms. S is a stoichiometric matrix and v is a reaction rate vector. The subscripts P and L denote the production and loss, respectively.
Generally, it can be said that if y_{ i }exhibits a quasisteady state behavior, such behavior is observed after a short period of time for the corresponding P_{ i }and L_{ i }to balance each other. In a normal computational environment, the period can be readily estimated from equation (1), but the applicability of QSSA cannot be directly determined based on the estimation. The goal of this work is to determine the proper criteria to determine the applicability of the QSSA based on the estimation of balancing period. At the moment either P_{ i }or L_{ i }enlarge, the period during which P_{ i }balances with L_{ i }can be evaluated in a simple manner.
By chain rule,
where f = (f_{1}, f_{2}, . . .)^{T}. Let δt_{ i }be a short time period after which P_{ i }and L_{ i }balance each other. Then, we have following relationship:
where superscript 0 indicates the reference value. Rearranging the equation gives:
If the denominator on the righthand side of equation (3) is not zero, we can compute the time scale δt_{ i }from this equation. If the magnitude of the time scale is large, namely, δt_{ i } > ϵ_{ t }for some ϵ_{ t }> 0, y_{ i }is considered to exhibit slow dynamics and QSSA is not applied.
If δt_{ i } <ϵ_{ t }, y_{ i }exhibits fast dynamics, QSSA can be applied depending on the sign of δt_{ i }.

If δt_{ i }is positive, it will reach a balancing state quickly and QSSA can be applied to y_{ i }.

If δt_{ i }is negative, the applicability of QSSA cannot be determined by δt_{ i }alone.
If the denominator is zero, the numerator L_{ i } P_{ i }should be considered for following three cases:

P_{ i }≠ L_{ i }indicates a nonreducing state between production and loss. QSSA should not be applied.

P_{ i }= L_{ i }= 0 indicates that no dynamics has occurred yet for y_{ i }. QSSA should not be applied.

P_{ i }= L_{ i }≠ 0 indicates that δt_{ i }should be set to zero since complete balance has been obtained.
The time scale of the element reaction can be also estimated using the same method, which will be used to determine a closed subsystem later.
Another scale factor should be used to determine whether the i'th variable is balanced or not; the ratio of f_{ i } to the larger one of P_{ i }and L_{ i },
If ϵ_{ t }= δt_{ i }< 0, the information about r_{ i }becomes important to determine the dynamics of y_{ i }. If r_{ i }is large, namely r_{ i }> ϵ_{ r }for some ϵ_{ r }> 0, the production and loss are neither balanced nor can be balanced soon, hence QSSA is not applicable to such y_{ i }.
Based on the values calculated from (3) to (5), the criteria to separate QSSA variables are as follows.

1.
Both production and loss terms must exist.

2.
For certain ϵ_{ r }≪ 1 and ϵ_{ t }≪ t_{ f }, the applicability of QSSA to the i th variable can be summarized as in Table 1.
Iterative balancing
The abovementioned scale factors δt_{ i }and r_{ i }are locally determined at a certain moment, initially in the innermost boundary layer. Therefore, if there are multiple boundary layers, the separation process should be applied iteratively to exit the fastest dynamic regime and move to the next fastest one, jumping to the next boundary layers, and finally to the slow dynamic regime, outer region, where fast dynamics of the inner regions are fully relaxed. The iterative procedure locates the variable values of the fast dynamics, y(0), to the outer region, y(0^{+}), which become the consistent initial values of the reduced system.
Once δt_{ i }and r_{ i }are computed, identifying the fast variables, one can decompose the original stoichiometric matrix into two submatrices, S_{ f }and S_{ s }. Equation (2) is rewritten as:
where subscript s indicates slow variables and f indicates fast ones. To derive consistent initial values for the outer region, the following algebraic equation should be solved:0 = S_{ f }v.
The solution of equation (6) can give other scale factors that also satisfy the approximation criteria, which shows the possibility of multiple boundary layers with different time scales. The existence of multiple boundary layers with different time scales corresponds to the cascaded nested hierarchy concept of inertial manifolds [8]. For these events, equations for the iterative approach can be written as:
where the superscript k is the number of iterations. From ${y}_{f}^{k}$, $\delta {t}_{i}^{k}$ and ${r}_{i}^{k}$ are iteratively calculated and the updated equation (7) is solved until y^{k+1} y^{k} satisfies the convergence criteria. The convergence criterion of this iterative procedure can be considered as the partial equilibrium among the fast variables. This is conceptually similar to the equilibrium value convergence of Lebiedz's work [11]. After convergence is achieved, the following reduced models with the modified initial values are derived:
The solution of the algebraic equation (7) can be computed using the appropriate numerical method such as Newton's method or other similar methods. However, those methods sometimes give physically meaningless solutions. To overcome this, Mass and Pope applied an iterative technique that successively constructs subODE systems and successively solves them [19]. This work applied a similar approach to that of Mass and Pope to find the plausible initial guesses for the algebraic equation (7).
From equation (7), the columns of S_{ f }, say s_{ j }, give information about the closing and opening of the subsystem by the j 'th reaction. If s_{ j }is composed of only nonnegative signs or only nonpositive signs, this system is opened by the j'th reaction. If the time scale of the opening j'th reaction, δτ_{ j }, is large, say δτ_{ j }> ϵ_{ t }, the reaction is excluded from the subsystem, resulting in a closed subsystem,
Since the closed system always has steady state, the k'th solution of the transient subsystem reaching a steady state is used as the initial guess for equation (7).
Overall process
The criteria are used to detect the existence of a boundary layer and the iterative balancing computes the initial values of the reduced model at the border of the outer region. There are two categories of QSSA possibilities. The first category is that the large deviation between P_{ i }and L_{ i }decreases very fast, which gives a relatively larger r_{ i }but a small δt_{ i }. The second is that P_{ i }and L_{ i }are almost balanced by a coincidental initial value. In this case, the approximation is also dependent on how small δt_{ i }is. Then the iterative balancing relocates the initial values toward the outer region. If there are multiple boundary layers, a few more iterations of the iterative balancing are required to search for the proper initial values for the reduced model. Once the iterative balancing converges, indicating every possible balancing is completed with small values for both of r_{ i }and δt_{ i }, it is assumed that the fast dynamics are relaxed and slow dynamics of the outer region begin with the updated initial values.
This process can also identify the boundary layers that occur internally, not only at the initial point, because the detection checks the possibility of QSSA at every step of the calculation. These features are illustrated in the next section. In summary, the overall detecting and balancing process is as follows:
Iterative balancing process

1.
Calculate δt_{ i }, δτ_{ j }, and r_{ i }using equations (3) through (5).

2.
If no species can be approximated, then go to step 4, or else go to step 3.

3.
Perform iterative balancing using equation (7).

4.
Solve the ODE/DAE model with equation (8).

5.
Repeat steps 1 to 4 until the current t reaches the userdefined t_{ f }.
Results and Discussion
The MichaelisMenten kinetics without inhibition,
and with inhibition,
are considered in this study. The parameters are (k_{1}, k_{2}, k_{3}, k_{4}, k_{5}) = (500000, 5, 1000, 100, 0.16) and the initial values are (e_{0}, s_{0}, c_{10}, c_{20}, p_{0}) = (1, 100, 0, 0, 0) and (e_{0}, s_{0}, c_{10}, c_{20}, p_{0}, i_{0}, ei_{0}) = (1, 100, 0, 0, 0, 1000, 0) [24]. As can be seen, there are two boundary layers at the initial region and near t = 900 (see Figure 1). Since p is only produced, its dynamics are not considered when searching for the fast balancing species. In the initial region, the estimated value of δt ≈ 3.96 × 10^{8} and species e, s, and c_{1} were selected as fast variables, as expected. The subsystem composed of the chosen species was opened by the second reaction, hence the second reaction was removed from the subsystem. Based on the solution from the first iteration, species c_{1} and c_{2} were selected. Since the δt of e and s remains small, there indices were maintained as fast variables. The updated subsystem was opened by the third reaction, and consequently, the reaction was excluded. After the second iteration, the solution converged and was stored, then the open system was solved in a small interval period. Since a very large δt of s was identified in this refining step, s was excluded from the fast variable set. Finally, e, c_{1}, and c_{2} were selected as QSSA variables before the second boundary layer. The values of δt_{ i }at each iteration after three iterations of the iterative process are listed in Table 2. At the second boundary layer, the predefined criteria gave another iterative process and relocated the solution toward the outer area in the same manner as described above. For comparison with the conventional manual QSSA approach, the time scales of each species were also derived by mathematical balancing [25]. The meaning of the timescales of the fast variables from the conventional derivation is the time to exit the boundary layer. Therefore, the summation of δt_{ i }for every iteration until the species enters the outer region is the direct comparative value of the conventional scales (see Table 3). A similar tendency is observed between the sums of δt_{ i }and the mathematical scales. As in Table 3, the sums of δt_{ i }and the mathematically scaled values indicate that there are two boundary layers near t ≈ 10^{8} and t ≈ 10^{3} before exiting the initial regions. A semilog plot of the full model simulation (Figure 2) supports this expectation.
The second example, the MichaelisMenten kinetics with inhibition, shows a boundary layer at the initial area only (see equation (11) and Figure 3) with a similar scale of δt to that of the noninhibition case. The dynamic behavior of the second model in the inner region of the initial boundary layer is more complex because of the effect of the inhibition. These complex dynamics of the second example require a few more iterations than that of the first example to exit the initial boundary layer.
The third example is the caspase activation model given by equation (12) in [26].
The reaction rate equations for equation (12) are written as v_{1} = k_{1}[c8*][c3], v_{2} = k_{2}[c3*][c8], v_{3} = k_{3}[c3*][IAP]  k_{3} [c3*~IAP], v_{4} = k_{4}[c3*][IAP], v_{5} = k_{5}[c8*], v_{6} = k_{6}[c3*], v_{7} = k_{7}[c3*~IAP], v_{8} = k_{8}[IAP]  k_{8}, v_{9} = k_{9}[C8]  k_{9}, v_{10} = k_{10}[c3]  k_{10}, v_{11} = k_{11}[c8*][BAR]  k_{11} [c8*~BAR], v_{12} = k_{12}[BAR]  k_{12} and v_{13} = k_{13}[c8a~BAR], where the kinetic constants are listed in [26]. There are also two boundary layers at the initial and internal regions, but with a much larger δt relative to the former cases; δt ≈ O(10^{1}) at the initial area and δt ≈ O(1) at the internal boundary layer (Figure 4).
Although it is difficult to recognize the initial boundary layer in the figure because of the small change in the concentrations relative to that of the internal boundary layer, the initial boundary layer exists in the third model, which can be overlooked by the simple observation of the full model simulation. The existence of the initial boundary layer of the third model could be detected by the proposed algorithm.
Figure 1 and Figure 3 show that the exact value of the parameters k_{1}, k_{2}, k_{3}, and k_{4} cannot be properly identified using measurement intervals larger than O(10^{8}). Figure 4 also illustrates that it is impossible to obtained the specific parameters of the original model with a measurement interval larger than O(1). Only the ratios of some species's concentrations can be observed. Therefore, the ratios are meaningful for these situations, not the detailed dynamics occurred within the period smaller than the measurement interval. The suggested scheme improves computational efficiency in the stiff inner region and extracts the information of experimentally meaningful QSSA concentrations for the corresponding species. Besides, by determining a measurement scale, the level of simplification can be controlled easily.
Conclusion
This work proposed new criteria for the timescale analysis and iterative balancing approach to develop an automatic simplification. This approach has different consistent initial values after model reduction and successfully found consistent initial values of the simplified DAE model using iterative balancing. With some examples, from small systems to practical systems, this scheme gave a successful reduction and found consistent initial values. If a whole cell is the system of the dynamic simulation, the network of reaction pathways will be more complex, the number of variables will be increased, and the simulation will be more difficult. Henceforth, it may be also important to relate the derived model to the experimental view, and this approach can give the criteria to classify the meaningful values from the original model.
References
 1.
Hoffmann A, Levchenko A, Scott ML, Baltimore D: The I κ BNF κ B signaling module: temporal control and selective gene activation. Science 2002, 298(5596):1241–1245.
 2.
Dahlquist G, Edsberg L, Sköllermo G, Söderlind G: Are the numerical methods and software satisfactory for chemical kinetics? Lect Notes Math 1982, 968: 149–164.
 3.
Schauer M, Heinrich R: Quasisteadystate approximation in the mathematical modelling of biochemical reaction networks. Math Biosci 1983, 65: 155–170.
 4.
Schuster R, Schuster S, Holzhutter HG: Simplification of complex kinetic models used for the quantitative analysis of nuclear magnetic resonance or Radioactive Tracer Studies. J Chem Soc Faraday Trans 1992, 88(19):2837–2844.
 5.
Schuster R, Holzhutter HG: Rapidequilibrium approximation applied to mathematical models of tracer dynamics in biochemical systems. Mathl Comput Modelling 1994, 19: 241–253.
 6.
Lee TY, Nitirahardjo S, Lee S: An analytic approach in kinetic modelling of ZigglerNatta polymerization of butadiene. J Appl Polym Sci 1994, 53: 1605–1613.
 7.
Okino MS, Mavrovouniotis ML: Simplification of mathematical models of chemical reaction systems. Mathl Comput Modelling 1998, 98(2):391–408.
 8.
Roussel MR, Fraser SJ: On the geometry of transient relaxation. J Chem Phys 1991, 94(11):7106–7113.
 9.
Roussel MR, Fraser SJ: Invariant manifold methods for metabolic model reduction. Chaos 2001, 11: 196–206.
 10.
Gorban AN, Karlin IV: Method of invariant manifold for chemical kinetics. Chem Eng Sci 2003, 58(21):4751–4768.
 11.
Lebiedz D: Computing minimal entropy production trajectories: An approach to model reduction in chemical kinetics. J Chem Phys 2004, 120(15):6890–6897.
 12.
Reinhardt V, Winckler M, Lebiedz D: Approximation of slow attracting manifolds in chemical kinetics by trajectorybased optimization approaches. J Phys Chem A 2008, 112: 1712–1718.
 13.
Vora N, Daoutidis P: Nonlinear model reduction of reaction systems with multiple time scale dynamics. American Control Conference, 2001. Proceedings of the 2001 2001, 6: 4752–4757.
 14.
Gorban AN, Karlin IV, Zinovyev AY: Constructive methods of invariant manifolds for kinetic problems. Phys Rep 2004, 396: 197–403.
 15.
Lebiedz D, Kammerer J, BrandtPollmann U: Automatic network coupling analysis for dynamical systems based on detailed kinetic models. Phys Rev E Stat Nonlin Soft Matter Phys 2005, 72(4 Pt 1):041911.
 16.
Lam SH: Using CSP to understand complex chemical kinetics. Combust Sci Technol 1993, 89: 375–404.
 17.
Lam SH, Guossis DA: The CSP method for simplifying kinetics. Int J Chem Kinet 1994, 26: 461–486.
 18.
Zagaris A, Kaper HG, Kaper TJ: Fast and slow dynamics for the computational singular perturbation method. Multiscale Model Simul 2004, 2(4):613–638.
 19.
Maas U, Pope SB: Simplifying chemical kinetics: intrinsic low dimensional manifolds in composition space. Combust Flame 1992, 88: 239–264.
 20.
Deuflhard P, Heroth J: Dynamic dimension reduction in ODE models. Technical report, KonradZuseZentrum fur Informationstechnik Berlin, PreprintSC95–29 1995.
 21.
Zobeley J, Lebiedz D, Kammerer J, Ishmurzin A, Kummer U: A New timedependent complexity reduction method for biochemical Systems. Transactions on Computational Systems Biology 2005, 1: 90–110.
 22.
Surovtsova I, Sahle S, Pahle J, Kummer U: Approaches to complexity reduction in a systems biology research environment (SYCAMORE). In WSC '06: Proceedings of the 38th Conference on Winter Simulation. Winter Simulation Conference; 2006:1683–1689.
 23.
Vallabhajosyula RR, Sauro HM: Complexity reduction of biochemical networks. In WSC '06: Proceedings of the 38th Conference on Winter Simulation. Winter Simulation Conference; 2006:1690–1697.
 24.
Hayashi K, Sakamoto N: Dynamic analysis of enzyme systems. Japan Scientific Societies Press/SpringerVerlag; 1986.
 25.
Segel LA, Slemrod M: The quasisteady state assumption: a case study in perturbation. SIAM Rev 1989, 31(3):446–477.
 26.
Eissing T, Conzelmann H, Gilles ED, Allgöwer F, Bullinger E, Scheurich P: Bistability analysis of a caspase activation model for receptorinduced apoptosis. J Biol Chem 2004, 279(35):36892–36897.
Acknowledgements
This work was supported by the Korean Systems Biology Project from the Ministry of Education, Science and Technology through KOSEF. We thank Ms. Trisha Poole for correcting the linguistic and typing errors of the manuscript.
Author information
Affiliations
Corresponding author
Additional information
Authors' contributions
JC developed and implemented the algorithm and prepared the manuscript. K–WY implemented the subroutines of the proposed algorithm. T–YL corrected the mathematical expression and participated in the manuscript preparation. SYL participated in the design of the study. 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
Cite this article
Choi, J., Yang, Kw., Lee, Ty. et al. New timescale criteria for model simplification of bioreaction systems. BMC Bioinformatics 9, 338 (2008). https://doi.org/10.1186/147121059338
Received:
Accepted:
Published:
Keywords
 Boundary Layer
 Fast Variable
 Differential Algebraic Equation
 Fast Dynamic
 Internal Boundary Layer