New timescale criteria for model simplification of bioreaction systems
 Junwon Choi^{1},
 Kyungwon Yang^{1},
 Taiyong Lee^{1}Email author and
 Sang Yup Lee^{1, 2}
https://doi.org/10.1186/147121059338
© Choi et al; licensee BioMed Central Ltd. 2008
Received: 06 March 2008
Accepted: 12 August 2008
Published: 12 August 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.
Keywords
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
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.
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.
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 }.
 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.
QSSA applicability.
0 ≤ δt_{ i }≤ ϵ_{ t }  ϵ_{ t }<δt_{ i }≤ 0  δt_{ i } >  ϵ _{ t }  

r_{ i }≤ ϵ_{ r }  Yes  Yes  No 
r_{ i }> ϵ_{ r }  Yes  No  No 
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.
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 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).
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
δt values.
Iteration  0  1  2 

e  3.96 × 10^{8}  3.09 × 10^{13}  3.57 × 10^{13} 
s  3.96 × 10^{8}  3.09 × 10^{13}  3.51 × 10^{2} 
c _{1}  3.96 × 10^{8}  1.81 × 10^{3}  4.57 × 10^{11} 
c _{2}  Not available  1.82 × 10^{3}  1.82 × 10^{3} 
Comparison of timescales.
Sums of δt  Mathematical scale  

e  3.96 × 10^{8}  2.00 × 10^{8} 
s  3.51 × 10^{2}  6.25 × 10^{2} 
c _{1}  1.81 × 10^{3}  1.00 × 10^{3} 
c _{2}  0.36 × 10^{2}  1.00 × 10^{2} 
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).
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.
Declarations
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.
Authors’ Affiliations
References
 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.View ArticlePubMedGoogle Scholar
 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.View ArticleGoogle Scholar
 Schauer M, Heinrich R: Quasisteadystate approximation in the mathematical modelling of biochemical reaction networks. Math Biosci 1983, 65: 155–170.View ArticleGoogle Scholar
 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.View ArticleGoogle Scholar
 Schuster R, Holzhutter HG: Rapidequilibrium approximation applied to mathematical models of tracer dynamics in biochemical systems. Mathl Comput Modelling 1994, 19: 241–253.View ArticleGoogle Scholar
 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.View ArticleGoogle Scholar
 Okino MS, Mavrovouniotis ML: Simplification of mathematical models of chemical reaction systems. Mathl Comput Modelling 1998, 98(2):391–408.Google Scholar
 Roussel MR, Fraser SJ: On the geometry of transient relaxation. J Chem Phys 1991, 94(11):7106–7113.View ArticleGoogle Scholar
 Roussel MR, Fraser SJ: Invariant manifold methods for metabolic model reduction. Chaos 2001, 11: 196–206.View ArticlePubMedGoogle Scholar
 Gorban AN, Karlin IV: Method of invariant manifold for chemical kinetics. Chem Eng Sci 2003, 58(21):4751–4768.View ArticleGoogle Scholar
 Lebiedz D: Computing minimal entropy production trajectories: An approach to model reduction in chemical kinetics. J Chem Phys 2004, 120(15):6890–6897.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 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.View ArticleGoogle Scholar
 Gorban AN, Karlin IV, Zinovyev AY: Constructive methods of invariant manifolds for kinetic problems. Phys Rep 2004, 396: 197–403.View ArticleGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 Lam SH: Using CSP to understand complex chemical kinetics. Combust Sci Technol 1993, 89: 375–404.View ArticleGoogle Scholar
 Lam SH, Guossis DA: The CSP method for simplifying kinetics. Int J Chem Kinet 1994, 26: 461–486.View ArticleGoogle Scholar
 Zagaris A, Kaper HG, Kaper TJ: Fast and slow dynamics for the computational singular perturbation method. Multiscale Model Simul 2004, 2(4):613–638.View ArticleGoogle Scholar
 Maas U, Pope SB: Simplifying chemical kinetics: intrinsic low dimensional manifolds in composition space. Combust Flame 1992, 88: 239–264.View ArticleGoogle Scholar
 Deuflhard P, Heroth J: Dynamic dimension reduction in ODE models. Technical report, KonradZuseZentrum fur Informationstechnik Berlin, PreprintSC95–29 1995.Google Scholar
 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.View ArticleGoogle Scholar
 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.Google Scholar
 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.Google Scholar
 Hayashi K, Sakamoto N: Dynamic analysis of enzyme systems. Japan Scientific Societies Press/SpringerVerlag; 1986.View ArticleGoogle Scholar
 Segel LA, Slemrod M: The quasisteady state assumption: a case study in perturbation. SIAM Rev 1989, 31(3):446–477.View ArticleGoogle Scholar
 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.View ArticlePubMedGoogle 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.