Open Access

New time-scale criteria for model simplification of bio-reaction systems

BMC Bioinformatics20089:338

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 [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 quasi-steady state approximation (QSSA) and a partial equilibrium approximation, have been applied to relieve the stiffness characteristics [37]. 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) [1618], 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 time-scales of the system. Currently, ILDM based methods have been applied to reduce complexity of biochemical systems [2123]. 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

In homogeneous chemical kinetics, the dynamic model can be written in the following form of the ODE.
d y i d t = f i ( y ) = j = 1 m s i j v j i = 1 , ... , n , 0 t t f . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeWabiqaaaqaaKqbaoaalaaabaGaemizaqMaemyEaK3aaSbaaeaacqWGPbqAaeqaaaqaaiabdsgaKjabdsha0baakiabg2da9iabdAgaMnaaBaaaleaacqWGPbqAaeqaaOGaeiikaGIaeCyEaKNaeiykaKcabaqbaeqabeWaaaqaaiabg2da9maaqahabaGaem4Cam3aaSbaaSqaaiabdMgaPjabdQgaQbqabaGccqWG2bGDdaWgaaWcbaGaemOAaOgabeaaaeaacqWGQbGAcqGH9aqpcqaIXaqmaeaacqWGTbqBa0GaeyyeIuoaaOqaaiabdMgaPjabg2da9iabigdaXiabcYcaSiabc6caUiabc6caUiabc6caUiabcYcaSiabd6gaUjabcYcaSaqaaiabicdaWiabgsMiJkabdsha0jabgsMiJkabdsha0naaBaaaleaacqWGMbGzaeqaaOGaeiOla4caaaaaaaa@5E13@
Since a chemical reaction system generally consists of production and loss terms, the ODE can be rewritten as:
d y i d t = P i ( y ) L i ( y ) = s P , i j v j s L , i j v j s P , i j = max ( s i j , 0 ) s L , i j = s P , i j s i j , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeWabqqaaaaabaqcfa4aaSaaaeaacqWGKbazcqWG5bqEdaWgaaqaaiabdMgaPbqabaaabaGaemizaqMaemiDaqhaaOGaeyypa0Jaemiuaa1aaSbaaSqaaiabdMgaPbqabaGccqGGOaakcqWH5bqEcqGGPaqkcqGHsislcqWGmbatdaWgaaWcbaGaemyAaKgabeaakiabcIcaOiabhMha5jabcMcaPaqaaiabg2da9maaqaeabaGaem4Cam3aaSbaaSqaaiabdcfaqjabcYcaSiabdMgaPjabdQgaQbqabaGccqWG2bGDdaWgaaWcbaGaemOAaOgabeaaaeqabeqdcqGHris5aOGaeyOeI0YaaabqaeaacqWGZbWCdaWgaaWcbaGaemitaWKaeiilaWIaemyAaKMaemOAaOgabeaakiabdAha2naaBaaaleaacqWGQbGAaeqaaaqabeqaniabggHiLdaakeaacqWGZbWCdaWgaaWcbaGaemiuaaLaeiilaWIaemyAaKMaemOAaOgabeaakiabg2da9iGbc2gaTjabcggaHjabcIha4jabcIcaOiabdohaZnaaBaaaleaacqWGPbqAcqWGQbGAaeqaaOGaeiilaWIaeGimaaJaeiykaKcabaGaem4Cam3aaSbaaSqaaiabdYeamjabcYcaSiabdMgaPjabdQgaQbqabaGccqGH9aqpcqWGZbWCdaWgaaWcbaGaemiuaaLaeiilaWIaemyAaKMaemOAaOgabeaakiabgkHiTiabdohaZnaaBaaaleaacqWGPbqAcqWGQbGAaeqaaOGaeiilaWcaaaaa@81E5@
or with a matrix-vector notation,
y ˙ = S v ( y ) = S P v S L v , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeWabiqaaaqaaiqbhMha5zaacaGaeyypa0JaeC4uamLaeeiiaaIaeCODayNaeiikaGIaeCyEaKNaeiykaKcabaGaeyypa0JaeC4uam1aaSbaaSqaaiabdcfaqbqabaGccqqGGaaicqWH2bGDcqGHsislcqWHtbWudaWgaaWcbaGaemitaWeabeaakiabbccaGiabhAha2jabcYcaSaaaaaa@41D3@

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.

By chain rule,
d P i d t = ( P i y ) T f d L i d t = ( L i y ) T f , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeWabiqaaaqaaKqbaoaalaaabaGaemizaqMaemiuaa1aaSbaaeaacqWGPbqAaeqaaaqaaiabdsgaKjabdsha0baakiabg2da9maabmaajuaGbaWaaSaaaeaacqGHciITcqWGqbaudaWgaaqaaiabdMgaPbqabaaabaGaeyOaIyRaemyEaKhaaaGccaGLOaGaayzkaaWaaWbaaSqabeaacqWGubavaaGccqWHMbGzaeaajuaGdaWcaaqaaiabdsgaKjabdYeamnaaBaaabaGaemyAaKgabeaaaeaacqWGKbazcqWG0baDaaGccqGH9aqpdaqadaqcfayaamaalaaabaGaeyOaIyRaemitaW0aaSbaaeaacqWGPbqAaeqaaaqaaiabgkGi2kabdMha5baaaOGaayjkaiaawMcaamaaCaaaleqabaGaemivaqfaaOGaeCOzayMaeiilaWcaaaaa@558A@
where f = (f1, f2, . . .) T . Let δt i be a short time period after which P i and L i balance each other. Then, we have following relationship:
L i 0 + d L i d t | y 0 δ t i P i 0 + d P i d t | y 0 δ t i , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemitaW0aa0baaSqaaiabdMgaPbqaaiabicdaWaaakiabgUcaRmaaeiaajuaGbaWaaSaaaeaacqWGKbazcqWGmbatdaWgaaqaaiabdMgaPbqabaaabaGaemizaqMaemiDaqhaaaGccaGLiWoadaWgaaWcbaGaeCyEaK3aaWbaaWqabeaacqaIWaamaaaaleqaaOGaeqiTdqMaemiDaq3aaSbaaSqaaiabdMgaPbqabaGccqGHijYUcqWGqbaudaqhaaWcbaGaemyAaKgabaGaeGimaadaaOGaey4kaSYaaqGaaKqbagaadaWcaaqaaiabdsgaKjabdcfaqnaaBaaabaGaemyAaKgabeaaaeaacqWGKbazcqWG0baDaaaakiaawIa7amaaBaaaleaacqWH5bqEdaahaaadbeqaaiabicdaWaaaaSqabaGccqaH0oazcqWG0baDdaWgaaWcbaGaemyAaKgabeaakiabcYcaSaaa@58AC@
where superscript 0 indicates the reference value. Rearranging the equation gives:
δ t i L i 0 P i 0 d P i d t | y 0 d L i d t | y 0 . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqiTdqMaemiDaq3aaSbaaSqaaiabdMgaPbqabaGccqGHijYUjuaGdaWcaaqaaiabdYeamnaaDaaabaGaemyAaKgabaGaeGimaadaaiabgkHiTiabdcfaqnaaDaaabaGaemyAaKgabaGaeGimaadaaaqaamaaeiaabaWaaSaaaeaacqWGKbazcqWGqbaudaWgaaqaaiabdMgaPbqabaaabaGaemizaqMaemiDaqhaaaGaayjcSdWaaSbaaeaacqWH5bqEdaahaaqabeaacqaIWaamaaaabeaacqGHsisldaabcaqaamaalaaabaGaemizaqMaemitaW0aaSbaaeaacqWGPbqAaeqaaaqaaiabdsgaKjabdsha0baaaiaawIa7amaaBaaabaGaeCyEaK3aaWbaaeqabaGaeGimaadaaaqabaaaaiabc6caUaaa@530B@

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.

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.
δ τ j ( d ln v j 0 d t ) 1 . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqiTdqMaeqiXdq3aaSbaaSqaaiabdQgaQbqabaGccqGHijYUcqGHsisldaqadaqcfayaamaalaaabaGaemizaqMagiiBaWMaeiOBa4MaemODay3aa0baaeaacqWGQbGAaeaacqaIWaamaaaabaGaemizaqMaemiDaqhaaaGccaGLOaGaayzkaaWaaWbaaSqabeaacqGHsislcqaIXaqmaaGccqGGUaGlaaa@43A2@
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 |,
r i = | f i | max ( P i , L i ) . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemOCai3aaSbaaSqaaiabdMgaPbqabaGccqGH9aqpjuaGdaWcaaqaamaaemaabaGaemOzay2aaSbaaeaacqWGPbqAaeqaaaGaay5bSlaawIa7aaqaaiGbc2gaTjabcggaHjabcIha4jabcIcaOiabdcfaqnaaBaaabaGaemyAaKgabeaacqGGSaalcqWGmbatdaWgaaqaaiabdMgaPbqabaGaeiykaKcaaiabc6caUaaa@4396@

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. 1.

    Both production and loss terms must exist.

  2. 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.

Table 1

QSSA applicability.


0 ≤ δt i ϵ t

-ϵ t <δt i ≤ 0

|δt i | > | ϵ t

r i ϵ r




r i > ϵ r




Applicability of QSSA to y i

Iterative balancing

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.

Once δt i and r i are computed, identifying the fast variables, one can decompose the original stoichiometric matrix into two sub-matrices, S f and S s . Equation (2) is rewritten as:
y ˙ s = S s v y ˙ f = S f v y = ( y s y f ) , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeWabmqaaaqaaiqbhMha5zaacaWaaSbaaSqaaiabdohaZbqabaGccqGH9aqpcqWHtbWudaWgaaWcbaGaem4CamhabeaakiabhAha2bqaaiqbhMha5zaacaWaaSbaaSqaaiabdAgaMbqabaGccqGH9aqpcqWHtbWudaWgaaWcbaGaemOzaygabeaakiabhAha2bqaaiabhMha5jabg2da9maabmaabaqbaeqabiqaaaqaaiabhMha5naaBaaaleaacqWGZbWCaeqaaaGcbaGaeCyEaK3aaSbaaSqaaiabdAgaMbqabaaaaaGccaGLOaGaayzkaaGaeiilaWcaaaaa@4831@

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:
0 = S y k v ( y f k ) , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeGimaaJaeyypa0JaeC4uam1aa0baaSqaaiabdMha5bqaaiabdUgaRbaakiabhAha2jabcIcaOiabhMha5naaDaaaleaacqWGMbGzaeaacqWGRbWAaaGccqGGPaqkcqGGSaalaaa@3AD0@
where the superscript k is the number of iterations. From y f k MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeCyEaK3aa0baaSqaaiabdAgaMbqaaiabdUgaRbaaaaa@3035@ , δ t i k MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqiTdqMaemiDaq3aa0baaSqaaiabdMgaPbqaaiabdUgaRbaaaaa@31D2@ and r i k MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemOCai3aa0baaSqaaiabdMgaPbqaaiabdUgaRbaaaaa@3029@ are iteratively calculated and the updated equation (7) is solved until |yk+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:
y ˙ s = S s v 0 S f v , y 0 = y ( 0 + ) . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeWabiGaaaqaaiqbhMha5zaacaWaaSbaaSqaaiabdohaZbqabaGccqGH9aqpcqWHtbWudaWgaaWcbaGaem4CamhabeaakiabhAha2bqaaaqaaiabicdaWiabgIKi7kabhofatnaaBaaaleaacqWGMbGzaeqaaOGaeCODayNaeiilaWcabaGaeCyEaK3aaSbaaSqaaiabicdaWaqabaGccqGH9aqpcqWH5bqEcqGGOaakcqaIWaamdaahaaWcbeqaaiabgUcaRaaakiabcMcaPiabc6caUaaaaaa@4634@

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 [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 non-negative signs or only non-positive 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,
y ˙ f k = S f c k v . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafCyEaKNbaiaadaqhaaWcbaGaemOzaygabaGaem4AaSgaaOGaeyypa0JaeC4uam1aa0baaSqaaiabdAgaMjabdogaJbqaaiabdUgaRbaakiabhAha2jabc6caUaaa@3966@

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. 1.

    Calculate δt i , δτ j , and r i using equations (3) through (5).

  2. 2.

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

  3. 3.

    Perform iterative balancing using equation (7).

  4. 4.

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

  5. 5.

    Repeat steps 1 to 4 until the current t reaches the user-defined t f .


Results and Discussion

The Michaelis-Menten kinetics without inhibition,
e + s k 2 k 1 c 1 k 4 k 3 c 2 k 5 p + e , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemyzauMaey4kaSIaem4Cam3aa4GdaSqaaiabdUgaRnaaBaaameaacqaIXaqmaeqaaaWcbaGaem4AaS2aaSbaaWqaaiabikdaYaqabaaakiaawcCicaGL9gcacqWGJbWydaWgaaWcbaGaeGymaedabeaakmaao4aaleaacqWGRbWAdaWgaaadbaGaeG4mamdabeaaaSqaaiabdUgaRnaaBaaameaacqaI0aanaeqaaaGccaGLahIaayzVHaGaem4yam2aaSbaaSqaaiabikdaYaqabaGcdaWfWaqaaiabgkziUcWcbaaabaGaem4AaS2aaSbaaWqaaiabiwda1aqabaaaaOGaemiCaaNaey4kaSIaemyzauMaeiilaWcaaa@4E60@
and with inhibition,
e + s k 2 k 1 e s 1 k 4 k 3 e s 2 k 5 p + e e + i k 2 k 1 e i , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeaabiWaaaqaaiabdwgaLjabgUcaRiabdohaZbqaamaao4aaleaacqWGRbWAdaWgaaadbaGaeGymaedabeaaaSqaaiabdUgaRnaaBaaameaacqaIYaGmaeqaaaGccaGLahIaayzVHaaabaGaemyzauMaem4Cam3aaSbaaSqaaiabigdaXaqabaGcdaGdoaWcbaGaem4AaS2aaSbaaWqaaiabiodaZaqabaaaleaacqWGRbWAdaWgaaadbaGaeGinaqdabeaaaOGaayjWHiaaw2BiaiabdwgaLjabdohaZnaaBaaaleaacqaIYaGmaeqaaOWaaCbmaeaacqGHsgIRaSqaaaqaaiabdUgaRnaaBaaameaacqaI1aqnaeqaaaaakiabdchaWjabgUcaRiabdwgaLbqaaiabdwgaLjabgUcaRiabdMgaPbqaamaao4aaleaacqWGRbWAdaWgaaadbaGaeGymaedabeaaaSqaaiabdUgaRnaaBaaameaacqaIYaGmaeqaaaGccaGLahIaayzVHaaabaGaemyzauMaemyAaKMaeiilaWcaaaaa@5FC0@
are considered in this study. The parameters are (k1, k2, k3, k4, k5) = (500000, 5, 1000, 100, 0.16) and the initial values are (e0, s0, c10, c20, p0) = (1, 100, 0, 0, 0) and (e0, s0, c10, c20, p0, i0, ei0) = (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 c1 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 c1 and c2 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, c1, and c2 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 time-scales 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 semi-log plot of the full model simulation (Figure 2) supports this expectation.
Figure 1

Simplification results (I). Concentrations of e, c1, and c2 for the Michaelis-Menten system; (a) full ODE model solution and (b) reduced model solution.

Figure 2

Existence of multiple boundary layers. (a) Semi-log plot of s and (b) that of e, c1, and c2 for the Michaelis-Menten system. The existence of two boundary layers at the initial region are observed.

Table 2

δt values.






3.96 × 10-8

-3.09 × 10-13

3.57 × 10-13


3.96 × 10-8

-3.09 × 10-13

3.51 × 102

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

δt values during the initial iterative process of two complex Michaelis-Menten kinetics

Table 3

Comparison of time-scales.


Sums of δt

Mathematical scale


3.96 × 10-8

2.00 × 10-8


3.51 × 102

6.25 × 102

c 1

1.81 × 10-3

1.00 × 10-3

c 2

0.36 × 10-2

1.00 × 10-2

Mathematically derived time scales and sums of δt i of two complex Michaelis-Menten kinetics

The second example, the Michaelis-Menten 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 non-inhibition 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.
Figure 3

Simplification results (II). Concentrations of e, es1, es2, and ei for the Michaelis-Menten system with inhibition; (a) full ODE model solution and (b) reduced model solution.

The third example is the caspase activation model given by equation (12) in [26].
d [ c 8 ] d t = v 2 v 9 d [ c 8 ] d t = v 2 v 5 v 11 d [ c 3 ] d t = v 1 v 10 d [ c 3 ] d t = v 1 v 3 v 6 d [ IAP ] d t = v 3 v 4 v 8 d [ c 3 ~ IAP ] d t = v 3 v 7 d [ BAR ] d t = v 11 v 12 d [ c 8 ~ BAR ] d t = v 11 v 13 MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeWabGqaaaaaaeaajuaGdaWcaaqaaiabdsgaKjabcUfaBjabbogaJjabiIda4iabc2faDbqaaiabdsgaKjabdsha0baakiabg2da9iabgkHiTiabdAha2naaBaaaleaacqaIYaGmaeqaaOGaeyOeI0IaemODay3aaSbaaSqaaiabiMda5aqabaaakeaajuaGdaWcaaqaaiabdsgaKjabcUfaBjabbogaJjabiIda4iabgEHiQiabc2faDbqaaiabdsgaKjabdsha0baakiabg2da9iabdAha2naaBaaaleaacqaIYaGmaeqaaOGaeyOeI0IaemODay3aaSbaaSqaaiabiwda1aqabaGccqGHsislcqWG2bGDdaWgaaWcbaGaeGymaeJaeGymaedabeaaaOqaaKqbaoaalaaabaGaemizaqMaei4waSLaee4yamMaeG4mamJaeiyxa0fabaGaemizaqMaemiDaqhaaOGaeyypa0JaeyOeI0IaemODay3aaSbaaSqaaiabigdaXaqabaGccqGHsislcqWG2bGDdaWgaaWcbaGaeGymaeJaeGimaadabeaaaOqaaKqbaoaalaaabaGaemizaqMaei4waSLaee4yamMaeG4mamJaey4fIOIaeiyxa0fabaGaemizaqMaemiDaqhaaOGaeyypa0JaemODay3aaSbaaSqaaiabigdaXaqabaGccqGHsislcqWG2bGDdaWgaaWcbaGaeG4mamdabeaakiabgkHiTiabdAha2naaBaaaleaacqaI2aGnaeqaaaGcbaqcfa4aaSaaaeaacqWGKbazcqGGBbWwcqqGjbqscqqGbbqqcqqGqbaucqGGDbqxaeaacqWGKbazcqWG0baDaaGccqGH9aqpcqGHsislcqWG2bGDdaWgaaWcbaGaeG4mamdabeaakiabgkHiTiabdAha2naaBaaaleaacqaI0aanaeqaaOGaeyOeI0IaemODay3aaSbaaSqaaiabiIda4aqabaaakeaajuaGdaWcaaqaaiabdsgaKjabcUfaBjabbogaJjabiodaZiabgEHiQiabc6ha+jabbMeajjabbgeabjabbcfaqjabc2faDbqaaiabdsgaKjabdsha0baakiabg2da9iabdAha2naaBaaaleaacqaIZaWmaeqaaOGaeyOeI0IaemODay3aaSbaaSqaaiabiEda3aqabaaakeaajuaGdaWcaaqaaiabdsgaKjabcUfaBjabbkeacjabbgeabjabbkfasjabc2faDbqaaiabdsgaKjabdsha0baakiabg2da9iabgkHiTiabdAha2naaBaaaleaacqaIXaqmcqaIXaqmaeqaaOGaeyOeI0IaemODay3aaSbaaSqaaiabigdaXiabikdaYaqabaaakeaajuaGdaWcaaqaaiabdsgaKjabcUfaBjabbogaJjabiIda4iabgEHiQiabc6ha+jabbkeacjabbgeabjabbkfasjabc2faDbqaaiabdsgaKjabdsha0baakiabg2da9iabdAha2naaBaaaleaacqaIXaqmcqaIXaqmaeqaaOGaeyOeI0IaemODay3aaSbaaSqaaiabigdaXiabiodaZaqabaaaaaaa@D50E@

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 [26]. There are also two boundary layers at the initial and internal regions, but with a much larger δt relative to the former cases; δtO(10-1) at the initial area and δtO(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 4

Simplification results (III). Concentrations of c8 and c8* for the caspase system; (a) the solution profile of c8, from the full ODE model (c8) and from the reduced model (c8_red) and (b) the solution profile of c8*, from the full ODE model (c8*) and from the reduced model (c8*_red).

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.

Authors’ Affiliations

Department of Chemical and Biomolecular Engineering (BK21 Program), KAIST
Center for Systems and Synthetic Biotechnology, Institute for the BioCentury, KAIST


  1. 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
  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.View ArticleGoogle Scholar
  3. 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
  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.View ArticleGoogle Scholar
  5. 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
  6. 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
  7. Okino MS, Mavrovouniotis ML: Simplification of mathematical models of chemical reaction systems. Mathl Comput Modelling 1998, 98(2):391–408.Google Scholar
  8. Roussel MR, Fraser SJ: On the geometry of transient relaxation. J Chem Phys 1991, 94(11):7106–7113.View ArticleGoogle Scholar
  9. Roussel MR, Fraser SJ: Invariant manifold methods for metabolic model reduction. Chaos 2001, 11: 196–206.View ArticlePubMedGoogle Scholar
  10. Gorban AN, Karlin IV: Method of invariant manifold for chemical kinetics. Chem Eng Sci 2003, 58(21):4751–4768.View ArticleGoogle Scholar
  11. 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
  12. 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
  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.View ArticleGoogle Scholar
  14. Gorban AN, Karlin IV, Zinovyev AY: Constructive methods of invariant manifolds for kinetic problems. Phys Rep 2004, 396: 197–403.View ArticleGoogle Scholar
  15. 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
  16. Lam SH: Using CSP to understand complex chemical kinetics. Combust Sci Technol 1993, 89: 375–404.View ArticleGoogle Scholar
  17. Lam SH, Guossis DA: The CSP method for simplifying kinetics. Int J Chem Kinet 1994, 26: 461–486.View ArticleGoogle Scholar
  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.View ArticleGoogle Scholar
  19. Maas U, Pope SB: Simplifying chemical kinetics: intrinsic low dimensional manifolds in composition space. Combust Flame 1992, 88: 239–264.View ArticleGoogle Scholar
  20. Deuflhard P, Heroth J: Dynamic dimension reduction in ODE models. Technical report, Konrad-Zuse-Zentrum fur Informationstechnik Berlin, PreprintSC95–29 1995.Google Scholar
  21. 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
  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.Google Scholar
  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.Google Scholar
  24. Hayashi K, Sakamoto N: Dynamic analysis of enzyme systems. Japan Scientific Societies Press/Springer-Verlag; 1986.View ArticleGoogle Scholar
  25. Segel LA, Slemrod M: The quasi-steady state assumption: a case study in perturbation. SIAM Rev 1989, 31(3):446–477.View ArticleGoogle Scholar
  26. 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


© Choi et al; licensee BioMed Central Ltd. 2008

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 (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.