New time-scale criteria for model simplification of bio-reaction systems
© Choi et al; licensee BioMed Central Ltd. 2008
Received: 06 March 2008
Accepted: 12 August 2008
Published: 12 August 2008
Quasi-steady state approximation (QSSA) based on time-scale analysis is known to be an effective method for simplifying metabolic reaction system, but the conventional analysis becomes time-consuming 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.
This work derived new time-scale factor by focusing on the problem structure. By mathematically reasoning the balancing behavior of fast species, new time-scale 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 time-scaling 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.
By successive calculation of the newly derived time-scale 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.
The dynamic simulation of bio-reaction 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 bio-reaction system consists of multiple reactions with various enzymes that have different turn-over 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 bio-systems strongly coupled and have multiple time-scales. 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 .
For this purpose, apart from a numerical analysis approach, the kinetic field's specific solution methods have been required . Traditionally, simplification of an original complex model, such as a quasi-steady 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) , the minimal entropy production trajectory (MEPT) based methods [11, 12], and the nonlinear model reduction method . Gorban et al. collected and reviewed such kinds of methods .
The others suggested the appropriate simplification criteria: generalized sensitivity analyses  used singular values of the sensitivity matrix as the scale factor, while computational singular perturbation (CSP) [16–18], intrinsic low dimensional manifold (ILDM) , 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 time-scales 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(n3) 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.
Criteria for balancing
where y ∈ ℝ+nis 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 quasi-steady 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 non-reducing 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 .
Both production and loss terms must exist.
For certain ϵ r ≪ 1 and ϵ t ≪ t f , the applicability of QSSA to the i th variable can be summarized as in Table 1.
0 ≤ δt i ≤ ϵ t
-ϵ t <δt i ≤ 0
|δt i | > | ϵ t
r i ≤ ϵ r
r i > ϵ r
The above-mentioned 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 sub-ODE systems and successively solves them . 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).
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
Calculate δt i , δτ j , and r i using equations (3) through (5).
If no species can be approximated, then go to step 4, or else go to step 3.
Perform iterative balancing using equation (7).
Solve the ODE/DAE model with equation (8).
Repeat steps 1 to 4 until the current t reaches the user-defined t f .
Results and Discussion
3.96 × 10-8
-3.09 × 10-13
3.57 × 10-13
3.96 × 10-8
-3.09 × 10-13
3.51 × 102
3.96 × 10-8
1.81 × 10-3
4.57 × 10-11
1.82 × 10-3
1.82 × 10-3
Comparison of time-scales.
Sums of δt
3.96 × 10-8
2.00 × 10-8
3.51 × 102
6.25 × 102
1.81 × 10-3
1.00 × 10-3
0.36 × 10-2
1.00 × 10-2
The reaction rate equations for equation (12) are written as v1 = k1[c8*][c3], v2 = k2[c3*][c8], v3 = k3[c3*][IAP] - k-3 [c3*~IAP], v4 = k4[c3*][IAP], v5 = k5[c8*], v6 = k6[c3*], v7 = k7[c3*~IAP], v8 = k8[IAP] - k-8, v9 = k9[C8] - k-9, v10 = k10[c3] - k-10, v11 = k11[c8*][BAR] - k-11 [c8*~BAR], v12 = k12[BAR] - k-12 and v13 = k13[c8a~BAR], where the kinetic constants are listed in . 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 k1, k2, k3, and k4 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.
This work proposed new criteria for the time-scale 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.
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.
- Hoffmann A, Levchenko A, Scott ML, Baltimore D: The I κ B-NF- κ 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: Quasi-steady-state 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: Rapid-equilibrium 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 Ziggler-Natta 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 trajectory-based 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, Brandt-Pollmann 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, Konrad-Zuse-Zentrum fur Informationstechnik Berlin, PreprintSC95–29 1995.Google Scholar
- Zobeley J, Lebiedz D, Kammerer J, Ishmurzin A, Kummer U: A New time-dependent 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/Springer-Verlag; 1986.View ArticleGoogle Scholar
- Segel LA, Slemrod M: The quasi-steady 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 receptor-induced apoptosis. J Biol Chem 2004, 279(35):36892–36897.View ArticlePubMedGoogle Scholar
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.