Skip to main content
  • Methodology article
  • Open access
  • Published:

A method for zooming of nonlinear models of biochemical systems

Abstract

Background

Models of biochemical systems are typically complex, which may complicate the discovery of cardinal biochemical principles. It is therefore important to single out the parts of a model that are essential for the function of the system, so that the remaining non-essential parts can be eliminated. However, each component of a mechanistic model has a clear biochemical interpretation, and it is desirable to conserve as much of this interpretability as possible in the reduction process. Furthermore, it is of great advantage if we can translate predictions from the reduced model to the original model.

Results

In this paper we present a novel method for model reduction that generates reduced models with a clear biochemical interpretation. Unlike conventional methods for model reduction our method enables the mapping of predictions by the reduced model to the corresponding detailed predictions by the original model. The method is based on proper lumping of state variables interacting on short time scales and on the computation of fraction parameters, which serve as the link between the reduced model and the original model. We illustrate the advantages of the proposed method by applying it to two biochemical models. The first model is of modest size and is commonly occurring as a part of larger models. The second model describes glucose transport across the cell membrane in baker's yeast. Both models can be significantly reduced with the proposed method, at the same time as the interpretability is conserved.

Conclusions

We introduce a novel method for reduction of biochemical models that is compatible with the concept of zooming. Zooming allows the modeler to work on different levels of model granularity, and enables a direct interpretation of how modifications to the model on one level affect the model on other levels in the hierarchy. The method extends the applicability of the method that was previously developed for zooming of linear biochemical models to nonlinear models.

Background

One of the main reasons for the rapid growth of the field of systems biology is that it makes extensive use of mathematical modeling [13]. This allows for a better handling of high complexity, which is an inherent property of all living systems. Using modeling, complex hypotheses can be formulated and tested in a more systematic manner than is possible using only biochemical reasoning [46]. However, even if one can obtain a detailed model of the system with a high predictive power, the model in itself does not automatically lead to a full understanding of the underlying biochemistry. One should for instance analyze the model to single out its essence, i.e., to identify those parts of the model that can be eliminated, while still preserving the model's crucial behavior. This latter task is referred to as model reduction, and it is the topic of this paper. There is an extensive literature available on the topic of model reduction. However, most of these studies have been done outside the field of systems biology, and since systems biology brings about new types of challenges, reduction of biochemical models is still in its early stages. Traditional engineering approaches like balanced truncation have focused on preserving the input-output profile in an optimal manner, both for linear [710], and for nonlinear [11] systems. However, these methods are not suitable for systems biology, because the reduced model has no natural interpretation in itself (nevertheless, some special cases where this problem can be circumvented have been identified [12, 13]). This lack of interpretation is a problem because systems biology models are usually developed to help characterizing the dominating parts and structure of the system, and not only to obtain a black-box predictor. Methods have therefore been developed with traditional chemical approaches that are more centered on reducing the internal dynamics of the system. These methods are typically based on a sensitivity analysis [1417], on time-scale separation [1821], or on the lumping of state variables [2226] (see [20] for a general review on model reduction). The perhaps most widely used method is lumping. Two of the main reasons for this are that an effective lumping scheme can be identified from basic properties of the model (e.g., the stoichiometry), and that lumped state variables are formed as easily interpretable pools of state variables in the original model. However, lumping does normally not come with the possibility of back-translation from the lumped state variables in the reduced model to the original state variables. In [27] we provided such relations. This means that we can take the result from a simulation of a reduced model, and without performing a new simulation, directly compute the corresponding trajectories of the desired original state variables. Because of this back-translation possibility, we refer to the resulting two models as two degrees of zooming of the same model. Nevertheless, like in other recent model reduction papers in systems biology [2832], the results in [27] were mainly developed with linear systems in mind. Linear systems virtually only appear in the cases of mono-molecular reaction networks and for models describing the probabilistic evolution of a single protein complex [27, 33]. However, already in [27] we proposed that zooming may in principle also be applicable to nonlinear models, but we did not derive formulae for back-translation. Note that a majority of the currently available systems biology models are in fact nonlinear.

With the method introduced in this paper, we provide the extension of the previously proposed method in [27] to nonlinear models. We show that new challenges arise due to the nonlinearities, but also how these challenges can be overcome, for instance with a wise choice of state variables in the reduced model. The method is demonstrated by application to two closed models of metabolic systems.

Methods

In this paper we present a more general version of the method that was introduced in [27], which is applicable to nonlinear models. We start with some basic definitions and key observations that are illustrated on a small example model, before we turn to the details of the method.

Basic Definitions and Assumptions

The method is developed for models of biochemical reaction systems on state space form that are based on nonlinear ordinary differential equations (ODEs)

= f ( x , p , u , t ) ,
(1)
y = h ( x , p , u , t ) ,
(2)

where t denotes time; the dot over x in Eq. (1) denotes derivative w.r.t. to time; the state vector x n ; the parameters p p ; the inputs u m ; the outputs y l ; and f and h are in general nonlinear functions. The state vector, whose individual elements are referred to as state variables, typically represents amounts or concentrations of chemical species, and the parameters commonly represent kinetic constants, initial conditions, or scaling factors. In this paper we are primarily interested in a comparison of the state variables between models (the original model and the reduced model), which means that the form of the nonlinear function in Eq. (2) is irrelevant for the application of our method. The right-hand side of Eq. (1) can be expressed as the stoichiometric matrix S n×qtimes a vector of reaction rates r = r(x, p, u, t); r q

= S r ( x , p , u , t ) .

The existence of separate time-scales are commonly utilized for reduction of biochemical models (e.g., by reduction of mass action kinetics to Michaelis-Menten kinetics). The typical approach is to investigate if subsets of the state variables are in steady state or in quasi-steady state (QSS). If state variable x i is in steady-state for t ≥ 0 it holds by definition that

i ( t ) = 0 ,
(3)

which implies that

x i ( t ) = x i ( 0 ) ,

which efficiently removes the state variable from the model, since it can be substituted for constant. If on the other hand the state variable x i is in QSS, there are terms on the right-hand side of the ODE that are much larger than the negligible term on the left-hand side. The approximation

f i ( x , p , u , t ) 0 .
(4)

is then commonly used to reduce the model. We refer to a state as fast in the time interval ≤ T0t < T1 if Eq. (4) is valid in this time interval, and holds for the class of all considered inputs to the system. Note that T1 = ∞ in the case that the systems remains in QSS, which may for example not be the case for models with switches (where, e.g., the values of a subset of the state variables may change when a certain condition is fulfilled) [34, 35].

Note that Eq. (3) (steady state) necessitates that Eq. (4) (quasi-steady state) is fulfilled, but not vice versa. Although QSS implies that (some of) the terms of the right-hand side of the ODE are large and leaves the left-hand side (derivative term) negligible, the derivative term may still be large enough for the state variable in QSS to change considerably during the time-span of a simulation; the key is that these changes mainly occur on a slow manifold.

Zooming of Linear Models

The concept of zooming was introduced in [27], and a method was presented that is applicable to linear time-invariant (LTI) models, which on state space form reads:

o : = A x + B u y = C x + D u ,

where A n × n, B n × m, C l × n, and D l × m. The method is based on the existence of at least one subset of state variables in o for which the internal dynamics is very fast with respect to the current time-scale of interest. An algorithm for automatic reduction of linear models that is based on the detection of such subsets, which are referred to as fast clusters, is presented in [27]. If the w state variables of a fast cluster are replaced by a single state variable x l L = i = 1 w x l i , we obtain a reduced version of the original model

r : r = A r x r + B r u y = C r x r + D u ,

where x r (n - w + 1), A r (n - w + 1) × (n - w + 1), B r (n - w + 1) × m, and C r l × (n - w + 1).

The fraction parameters, which are typically computed from QSS assumptions and mass conservation relations, take the form

η l i ( p ) = x l i x l L ,
(5)

The fraction parameters are used for back-translation of the lumped state variable to the original state variables. Note that the fraction parameters are functions of the model parameters only, and therefore time-invariant; as we will see, these fraction parameter properties do in general not hold for nonlinear models. By comparing the reactions of the original and reduced models, we see that

k j L = i = 1 w k j i η i ,
(6)

where k j L is the rate parameter in the reaction from state variable x l L to state variable x j in the reduced model, and k ji is the rate parameter in the reaction from state variable x i to state variable x j in the original model.

Finally, note that Eqs. (5) and (6) provide a link between o and r , which constitute two different levels of granularity. It is this link between the models that make us consider them as two different degrees of zooming, and the primary goal of this paper is to establish such a link also for nonlinear models.

Extension to Nonlinear Models - Initial Observations

We will now present some key observations that are used in the derivation of the method for zooming of nonlinear models.

First observe that the mass, which corresponds to a weighted (w.r.t. the molecular weight) sum of the state variables, of a closed (no exchange of matter with the surroundings) nonlinear model is conserved. However, the total number of molecules is in general not conserved in such a model as it is for linear models. This is for example due to the formation and dissociation of complexes, which alters the total number of molecules in the system. For instance, the binding of A to B reduces the number of molecules, as the product AB only counts as one molecule; binding reactions cannot occur in linear models. A second, and related, observation is that another type of conservations appears in nonlinear models; conserved moieties. A moiety is a specific functional part of a molecule, and the weighted sum of the number of molecules that contain this functional part is constant in a closed system. The presence of such a conserved moiety is equivalent to the existence of a row vector m n for which mS = 0, which also implies that

m = m S r ( x , p ) = 0 .
(7)

If we let the rank of S be denoted by n r , the number of linearly independent vectors for which Eq. (7) holds is equal to n - n r , which implies the existence of a matrix M

M S = 0 ,
(8)

where M n - n r × n .

Let us now make some remarks regarding fast state variables in a nonlinear model. Let x f n f be the vector of all fast state variables in T0t < T1. For simplification we will assume that there are no inputs to the system, although it would in principle be possible to incorporate inputs in the following discussion. The right-hand side of the ODEs for these fast state variables, if there are no inputs, can be separated into two parts. The first part contains reactions between fast state variables that are significant for the fast dynamics; r f (x f , p), and the second part contains all other reactions, r s (x, p), i.e.,

f = S s r s ( x , p ) + S f r f ( x f , p ) ,
(9)

where S f and S s are the corresponding stoichiometric matrices. Let us now consider the fast stoichiometric matrix, S f , and especially the conserved moieties that are implied by S f . Since these moieties are only (approximately) conserved on a fast enough time-scale, we refer to such moiety conservations as apparent conservations. Let M f be a matrix with a linearly independent rows such that

M f S f  = 0,
(10)

where M f a × n f . Each row of this matrix thus implies an apparent conserved moiety in the system. Let the sums of state variables that correspond to apparent moiety conservations (i.e., lumps of state variables) be denoted by l, so that

l  =  M f x f .
(11)

If we differentiate l with respect to time, we get

l ° = M f f = M f ( S s r s ( x , p ) + S f r f ( x f , p ) ) M f S s r s ( x , p ) ,
(12)

where Eqs. (9) and (10) were used.

It is interesting to note that the matrix M f is not unique, but that in fact any matrix M ^ =N M f can be used for lumping, where N a×ais non-singular. This observation allows us to choose a matrix M ^ for which a maximal number of rows in M ^ S s vanish, which results in the greatest possible reduction in the number of state variables. Finally note that Eq. (4), in the absence of inputs, gives

f f ( x , p , u t , t ) = f f ( x , p , 0 , t ) = (1) = S s r s ( x , p ) + S f r f ( x f , p ) (2) S f r f ( x f , p ) 0 (3)
(13)

since the term S f r f (x f , p) dominates the term S s r s (x, p). Note that Eq. (4) and consequently Eq. (13) only hold in T0t < T1, since the state is only known to be fast in this time span.

Eq. (12) defines the ODEs of the reduced model (the lumped state variables), and Eqs. (11) and (13) can in principle be used to calculate back-translation formulae, as is demonstrated with the small example model in the next section. However, as we shall see, this approach requires the explicit algebraic solution to a system of nonlinear equations, which is typically an infeasible task. Furthermore, there is not a clear one-to-one mapping between the state variables of the original and reduced models as in the case of proper lumping [27].

A Small Example Model

We will now present a small example model, with three fast state variables, which is reduced with the approach discussed above. An alternative approach is then demonstrated with the advantage that it scales better to larger models.

Consider the reversible formation of a complex C from a substrate A and an enzyme B

A + B k - 1 k 1 C

consisting of the fast state variables x f = (A B C) T , where the bullets (•) represent the slow state variables surrounding the three fast state variables in the model. The ODEs for the fast state variables take the form

f = - 1 - 1 1 ( k 1 A B - k - 1 C ) + S s r s ( x , p ) .
(14)

The three state variables in the model constitute a fast cluster with two apparent conserved moieties, which may be represented by the following relations

l = L 1 L 2 = 1 0 1 0 1 1 A B C = M f x f ,
(15)

where the lumped state variables L1 and L2 are introduced. Note that Eq. (12) defines the dynamics of the lumped state variables.

The distribution of mass among the fast state variables is given by Eq. (15) and by applying Eq. (13) to (14), which results in an equation system with the three fast state variables as unknowns

k 1 A B - k - 1 C 0 ,
(16)
A  +  C  =  L 1 ,
(17)
B  +  C  =  L 2 .
(18)

Analytic expressions for the fast state variables A, B, and C are given by the non-negative solution to Eqs. (16)-(18)

A 1 2 ( L 1 - L 2 - K 1 + ( L 1 + L 2 + K 1 ) 2 - 4 L 1 L 2 ) ,
(19)
B 1 2 ( - L 1 + L 2 - K 1 + ( L 1 + L 2 + K 1 ) 2 - 4 L 1 L 2 ) ,
(20)
C 1 2 ( L 1 + L 2 + K 1 + ( L 1 + L 2 + K 1 ) 2 - 4 L 1 L 2 ) ,
(21)

where K 1 = k - 1 k 1 .

We can now employ Eq. (12) to solve the ODEs for the lumped state variables L1 and L2, and use Eqs. (19)-(21) as back-translation formulae to compute the trajectories of the original state variables A, B, and C. However, note that for even slightly larger clusters of fast species than the one discussed here it would not be possible to calculate algebraic expressions of the original state variables with this approach, since it builds on the explicit solution of a system of nonlinear equations, which quickly becomes infeasible with growing problem size.

Alternatively, we can take an approach to the problem that is inspired by the method for linear systems in [27]. The first step is to express Eqs. (16) and (17) as a linear system w.r.t. the state variables A and C

k 1 B - k - 1 1 1 A C 0 L 1 .
(22)

The solution to Eq. (22) w.r.t. A and C is

A η A ( B , p ) L 1 K 1 B + K 1 L 1 ,
(23)
C η C ( B , p ) L 1 B B + K 1 L 1 ,
(24)

where K 1 = k - 1 k 1 , and the fraction parameters η A (B, p) and η C (B, p) are defined in Eqs. (23) and (24), respectively. The ODE for L1 is defined by Eq. (12), and the ODE for B can be derived by differentiation of L2 in Eq. (18), which gives that

d B d t ( 1 + K 1 L 1 ( B + K 1 ) 2 ) - 1 - B B + K 1 1 L 1 d t L 2 d t = (1) = ( 1 + K 1 L 1 ( B + K 1 ) 2 ) - 1 - B B + K 1 1 M f S s r s ( x , k ) ,
(25)

The reduced model consists of the two state variables L1 and B (note that L2 does not appear in the reduced model), and the dynamics is described by Eqs. (12) and (25), respectively. Note that the state variables A and C can be back-translated from the reduced model with Eqs. (23) and (24). This approach is a bit more intricate than the first, but comes with the advantage that we do not need to solve a system of nonlinear equations.

A Method for Zooming of Nonlinear Models

We will now step-by-step present a method that can be used to construct zoomable nonlinear biochemical models. This involves two sub-goals: i) to identify a reduced model that shares important characteristics with the original model, ii) to derive back-translation formulae that can be used to compute the original state variables and parameters from the reduced model.

In an initialization step of the method for a model o we first formulate mathematical equations for all conservation relations Eq. (8), state variables in steady-state Eq. (3), and quasi-steady state assumptions Eq. (4). If additional properties of the system are known, we also formulate the corresponding equations.

Step 1

The first step of the method is to identify the apparent conservation relations in the model.

Definition 1: Let S f be the stoichiometric matrix for the reactions r f (x f , p) as defined in Eq. (9). Each subset of state variables for which the corresponding rows of S f are linearly dependent constitutes an apparent conservation relation. Hence the apparent conservation relations lie in the left null space of S f and the dimension of this space is n - rank(S f ).

Note that the apparent conservation relations are defined in Eq. (11). It is trivial to identify the set of all linearly dependent rows of S f with a mathematical computing software (e.g., SBtoolbox for Matlab [36]).

Step 2

The second step of the method is to define the state variables of the reduced model, which we refer to as modified lumped state variables.

Definition 2: Let x be a lumped state variable corresponding to a subset of the state variables in an apparent conservation relation. Then x is a modified lumped state variable if the lumping scheme with respect to the state variables of the original model is proper.

Note that the original state variables have a clear interpretation in the reduced model (i.e., that the lumped variables form disjoint sets) if the lumping scheme is proper, i.e.,

l m = M m x f
(26)

where M m is a × n f matrix with elements equal to 0 or 1 and column sums equal to 1, and l m denotes the modified lumped state variables. We typically have a large freedom in the choice of M m . The number of state variables is maximally reduced if all exact conservation relations in the model are retained as modified lumped state variables (and replaced by constants).

Step 3

The third step of the method is to derive fraction parameters, which constitute the link between the reduced model and the original model. Let the original state variables that constitute the k:th modified lumped state variable l m k be denoted by x m k , so that

l m k = i = 1 w x m k i ,
(27)

A number of n m equations that are linear w.r.t. x m k , and linearly independent, are required to calculate fraction parameters. The existance of n m such equations results in an equation system

b k ( l m , p )  =  A ( l m , p ) x m k ,
(28)

where both A ( l m , p ) R n m × n m and b k ( l m , p ) R n m are known, although some of the equations may in general be approximate (e.g., QSS). The matrix A(l m , p) is invertible since the equations are linearly independent, and we have that

x m k = A - 1 ( l m , p ) b k ( l m , p ) .
(29)

The fraction parameters can then be calculated

η m k i ( l m , p ) = x m k i l m k = x m k i i = 1 n m x m k i ,
(30)

where we used Eq. (27) in the last step.

A modified lumped state variable for which an insufficient number of linear and linearly independent equations are available may still be used in the reduced model. However, the back-translation of the modified lumped state variable to the original state variables is then not possible, and step 3 of the method is ignored.

Step 4

The fourth step of the method is to derive the rate of change of the modified lumped state variables. Theorem: The dynamics of the modified lumped state variables is given by

l ° m = ( I + J ( l m , p ) ) - 1 l ° = (1) = ( I + J ( l m , p ) ) - 1 M f S s r s ( x , p ) , (2)
(31)

where S s and r s (x, p) were defined in Eq. (9), the matrix M f is defined in Eq. (11), and

J i j ( l m , p ) = k ( M f i k - M m i k ) g k ( l m , p ) l m j ,
(32)

where the matrix M m is defined in Eq. (26), and g i ( l m , p ) x f i is introduced to simplify the notation.

Proof:

First subtract Eq. (26) from Eq. (11)

l = l m + ( M f - M m ) x f = l m + ( M f - M m ) g ( l m , p ) ,

and differentiate l with respect to time, which gives

l °  =  l ° m  +  J ( l m ,p) l ° m  = ( I  +  J ( l m , p )) l ° m ,
(33)

where I is the identity matrix and J(l m , p) is the Jacobian of (M f - M m )g(l m , p) with respect to l m . The element J ij of J(l m , p) is given by

J i j ( l m , p ) = k ( M f i k - M m i k ) g k ( l m , p ) l m j .

From Eq. (33) it is straight-forward to derive l ° m , which takes the form

l ° m = ( I + J ( l m , p ) ) - 1 l ° = (1) = ( I + J ( l m , p ) ) - 1 M f S s r s ( x , p ) , (2)

where Eq. (12) was used in the last step. □

The matrix I + J(l m ; p) is symbolically invertible, but may in general contain singularities for particular combinations of parameters values and state variable values. However, the matrix is always invertible for the models discussed in this paper, since the corresponding determinants are strictly positive.

Step 5

The final step of the method is to back-translate the modified lumped state variables to the original state variables with the fraction parameters derived in step 3. This allows a comparison between the predictions by the reduced model to those of the original model.

The implementation of the method is straight-forward, and we have used Matlab (R2008b) together with the SBtoolbox [36] as computing software for the models in this paper.

Results

We will now demonstrate the method through application to two example models.

Enzyme Kinetics Model

The model below describes the process of conversion of a substrate, S, into a product, P, which is catalyzed by an enzyme, E.

S + E k - 1 k 1 C S k 2 C P k - 3 k 3 P + E

Note that the complexes C s and C p are formed by S bound to E, and P bound to E, respectively. This model is frequently occurring as part of larger models of biological systems, although the reaction from C S to C P is sometimes neglected, or reversible. The ODEs for the model are listed in Appendix A.1, where the three reactions are defined as: r1 = k1SE - k-1C S , r2 = k2C S , and r3 = k3C P - k-3PE.

The reaction terms r1(x, p) and r3(x, p) are assumed to be dominating, and the reaction term r2(x, p) to be insignificant in the ODEs. This results in that all state variables are in QSS, which gives

k 1 S E - k - 1 C S 0 ,
(34)
k 3 P E - k - 3 C P 0 .
(35)

We denote the sum of the state variables containing the enzyme by

L E = E + C S + C P ,
(36)

which is constant since the total amount of the enzyme E is conserved in the system.

Reduction of the Enzyme Kinetics Model

The first step of the method is to identify the apparent conservation relations from the matrix S f . Since r2(x, p) is dominated by r1(x, p) and r3(x, p) the model ODEs can be written on the form of Eq. (9)

Ė C S ° C P ° = S s r s ( x , p ) + S f r f ( x f , p ) = = 0 0 0 - 1 1 r 2 + - 1 0 - 1 1 0 1 1 0 0 - 1 r 1 r 3 ,

A basis of the left null space of S f is given by the row vectors of M f , which is defined by

l = L S L P L E = (1) = 1 0 0 1 0 0 0 1 0 1 0 1 0 1 1 S E P C S C P = (2) = M f x f , (3)
(37)

where L S and L P are apparent conservation relations and LE is an exact conservation relation.

The second step is to define the modified lumped state variables on the form of Eq. (26). The number of state variables is maximally reduced if L E is retained as a state variable in the reduced model (i.e., since L E , unlike L S and L P , can be replaced by a constant). The vector of modified lumped state variables then is defined

l m = S P L E = 1 0 0 0 0 0 0 1 0 0 0 1 0 1 1 S E P C S C P = = M f x f ,

In the third step of the method we calculate fraction parameters for the modified lumped state variable L E . There are five equations (Eqs. (34)-(35) and (37)) that are linear w.r.t. the state variables E, C S , and C P , which are lumped into the state variable L E . Note that only n m = 3 equations are required to derive fraction parameters, and we use Eqs. (34)-(36) to formulate an equation system as in Eq. (28) with the solution

E = η E L E 1 1 + M 1 S + M 3 P L E g 2 ( l m , p ) ,
(38)
C S = η C S L E M 1 S 1 + M 1 S + M 3 P L E g 4 ( l m , p ) ,
(39)
C P = η C P L E M 3 P 1 + M 1 S + M 3 P L E g 5 ( l m , p ) ,
(40)

where M 1 = k 1 k - 1 and M 3 = k 3 k - 3 . The two remaining modified lumped state variables correspond to S and P in the original model, so we define that l m 1 =S g 1 ( l m , p ) and l m 2 =P g 3 ( l m , p ) .

In the fourth step we derive the rate of change of the modified lumped state variables. Eq. (32) gives that L ° E =0, which is replaced by a constant, and

- k 2 M 1 S L E ( ( 1 + M 1 S + M 3 P ) 2 + M 3 L E ) ϕ ( S , P , L E , M 1 , M 3 ) ,
(41)
k 2 M 1 S L E ( ( 1 + M 1 S + M 3 P ) 2 + M 1 L E ) ϕ ( S , P , L E , M 1 , M 3 ) .
(42)

where

ϕ ( S , P , L E , M 1 , M 3 ) = ( ( 1 + M 1 S + M 3 P ) 3 + ( M 1 + M 3 + M 1 M 3 ( P + S ) ) ( 1 + M 1 S + M 3 P ) L E + M 1 M 3 L E 2 ) .

The two ODEs in Eqs (41)-(42) define the dynamics of the state variables in the reduced model. We finally note that the exact conservation relation for the substrate, L T = S + P + C S + C P , together with Eqs. (39)-(40) can be used to reduce the model further to a single state.

In the fifth step of the method we use the fraction parameters, defined in Eqs. (38)-(40), to back-translate the modified lumped state variables to the state variables of the original model. A comparison between predictions of the original state variables, from simulations of the original model and the reduced model, is presented in Figure 1. Implementations of the original model (Additional file 1), the reduced model (Additional file 2), and a script for simulation with SBtoolbox2 for MATLAB [36] (Additional file 3), are available in Additional files.

Figure 1
figure 1

Small example model. A comparison between the state variables of the original enzyme kinetics model and the backtranslated state variables of the reduced version of the same model.

The only assumption that was used in the derivation of the reduced model is that the reaction terms r1(x, p) and r3(x, p) dominate the reaction term r2(x, p), which results in that all state variables are in QSS. To assess the impact of these assumptions on the reduced model we compute the relative difference between the state variables in the original and in the reduced model

ε i ( t ) = | x i o ( t ) - x i r ( t ) | x i o ( t ) ,i=1,,n,

where x i o is state variable i in the original model, x i r is the corresponding back-translated state variable in the reduced model, and |x| denotes the elementwise absolute values of x. The maximal mean and infinity norm of ε i (t) in Eq. (43) over time is presented in Table 1 for parameter values over five orders of magnitude. In general, the reduced model appears to be robust to changes in the parameter values, although slightly more sensitive to some parameters (e.g., small values of k-1, large values of k1, or large values of k2, which violate the assumptions used in the reduction). However, note that the validity of the QSS assumption may also depend on the state variables, for example the total concentration of the enzyme. Interestingly, we observed that the reduced model can well approximate the original model over several order of magnitudes around the nominal enzyme concentration (L E = 1). It is well-known that the QSS approximation is only valid for sufficiently small enzyme concentrations, and as expected the performance of the reduced model starts to decrease for immense enzyme concentrations.

Table 1 Robustness of the reduced model for large deviations from the nominal parameter point are presented for the enzyme kinetics model, with a sampling frequency of 0.1 (starting from 0.1) time units.

Note that all the state variables of the original model have a direct biological interpretation also in the reduced model, and that Eqs. (38)-(40) can be used to back-translate the state variables. The reduced model may be depicted

S + η E L E k - 1 k 1 η C S L E k 2 η C P L E k - 3 k 3 P + η E L E

where the fraction parameters specify the distribution of the enzyme among the corresponding original state variables.

Glucose Transport in Budding Yeast

A model for the transport of glucose into a cell of baker's yeast (S. cerevisiae), which constitutes the first step of glycolysis, is presented in [37]. The inflow of glucose is modeled as a facilitated diffusion process, in which a carrier enzyme is responsible for the transport between the inner and outer regions of the cellular membrane. It is assumed that glucose 6-phosphate (G6P) has an inhibitory role in the glucose transport process by binding to the transporter. A graphical representation of the model is shown in Figure 2, and the ODEs for the state variables are listed in Appendix A.2.

Figure 2
figure 2

Glucose transport model. The original model for glucose transport in baker's yeast (S. cerevisiae). This figure was originally presented in [27].

In [27] we described how the calculation of fraction parameters, based on a set of assumptions, leads to the same reaction rates in the reduced model as were reported in [37]. The assumptions are that state variables participating in reactions for uptake and release of glucose and G6P across the cell membrane are in QSS, that the transporter is conserved, and that the concentrations of the transporter in the inner and outer regions of the cellular membrane are constant.

The assumption that the state variables x G l c e , x G l c i , x E - G 6 P i , and x E - G l c - G 6 P i , which participate in the uptake and release of G6P and glucose across the cell membrane, are in QSS gives that

k 1 x E e x G l c e - k - 1 x E - G l c e 0 ,
(43)
k 2 x E i x G l c i - k - 2 x E - G l c i 0 ,
(44)
k 4 x E i x G 6 P i - k - 4 x E - G 6 P i 0 ,
(45)
k 3 x E - G l c i x G 6 P i - k - 3 x E - G l c - G 6 P i 0 .
(46)

We have the following exact conservation relations in the model

L E = x E - G 6 P i + x E - G l c - G 6 P i + x E - G l c e + (1) x E - G l c i + x E e + x E i , (2) (3)
(47)
L G l c = x G l c e + x G l c i + x E - G l c - G 6 P i + x E - G l c e + x E - G l c i ,
(48)
L G 6 P = x G 6 P i + x E - G 6 P i + x E - G l c - G 6 P i ,
(49)

where L E , L Glc , and LG 6Pare constant over time. The assumption in [37] that the concentrations of the transporter in the inner and outer regions of the cell membrane are constant is formulated

α ( x E - G l c e - x E - G l c i ) + β ( x E e - x E i ) = 0 .
(50)

It is not clear how equations for back-translation of the state variables in the reduced model in [27, 37] can be derived. The reduced model has three state variables; external- and internal glucose and G6P, but only two differential equations for the in- and outflow of glucose, since the ODE for G6P is replaced by a representative function that is inferred from the G6P data. Our method does not rely on that such information is available, although it would in principle be possible to utilize data fitted functions for the state variables. Also note that the equations in the reduced model describe the total influx and efflux of glucose across the membrane [27], which cannot be interpreted w.r.t. the state variables of the original model. Other assumptions in [37] that complicates a comparison with our method are that the efflux of glucose is negligible, that the concentration of glucose in the cytosol is negligible, and that the concentrations of the transporter are constant in the inner and outer regions of the cell membrane (Eq. (50)).

It is not possible to generate a reduced model by direct substitution of the fraction parameters that were derived in [27] into the ODEs of the original model, since this would lead to the prediction that the state variables are constant (as discussed in [27]). We will now instead illustrate how our method can be used to derive a reduced and zoomable version of the glucose transport model.

Reduction of the Glucose Transport Model

Before applying our method to the glucose transport model we tried an alternative approach. Eqs. (43)-(46) were solved w.r.t. the state variables in QSS, and the resulting expressions were then substituted into the remaining ODEs. The details of the derivation of the reduced model are presented in Appendix A.3. The reduced model does not produce satisfactory predictions for any other state variable than x G 6 P i , which remains approximately constant during the simulation. Implementations of the original model (Additional file 4), the reduced model (Additional file 5), and a script for simulation with SBtoolbox2 for MATLAB [36] (Additional file 6), are available in Additional files.

Since the first approach turned out to be insufficient for reduction of the glucose transport model we applied our method to the same model. Following [37], we initially assumed constant transporter concentrations in the inner and outer regions of the cellular membrane, as defined by Eq. (50). For the details on the derivation of the reduced model we refer to Appendix A.4. The reduced model clearly performs better than the model resulting from the first approach, but it is still not satisfactory. However, the assumption of constant regional concentrations of the transporter may not be valid since the transport of glucose across the cell membrane is a rate limiting step in the model, and appears to be important for the state variable dynamics. We therefore decided to neglect Eq. (50) in the reduction process.

Implementations of the original model (Additional file 4), the reduced model (Additional file 7), and a script for simulation with SBtoolbox2 for MATLAB [36] (Additional file 8), are available in Additional files. In the first step of the method we identify the following apparent conservation relations

l = L G l c 1 L E 1 L G l c 2 L G 6 P L E 2 = M x = (1) = 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 1 0 0 1 0 1 0 0 1 0 0 0 0 1 1 1 0 0 0 0 0 0 1 1 0 0 1 0 1 x .
(51)

We note that there are two disjoint clusters of fast reactions in the model, corresponding to the outer- and inner parts of the cell membrane.

In the second step we define the modified lumped state variables. We decide to keep the lumped state variables L E 1 and LG 6Pas modified lumped state variables. The choice to keep LG 6Pleads to the largest possible reduction in the number of state variables, since the conservation of LG 6Pis exact (which is not true for any other state variable in l). The modified lumped state variables are defined true for any other state variable in l). The modified lumped state variables are defined

l m = x G l c e L E 1 x G l c i L G 6 P L E 3 = M m x = (1) = 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 1 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 1 0 1 x ,
(52)

where

x = x G l c e x G l c i x E - G 6 P i x E - G l c - G 6 P i x G 6 P i x E - G l c e x E - G l c i x E e x E i .

Note that two of the state variables in the original model, x G l c e and x G l c i , are also modified lumped state variables.

In the third step of the method we calculate fraction parameters for the modified lumped state variables that correspond to more than one of the original state variables (i.e., L E 1 , LG 6P, and L E 3 ). All of the modified lumped state variables satisfy the requirement that at least n m of Eqs. (43)-(49) and Eq. (52) are linear, and linearly independent, with respect to the corresponding original state variables. Eqs. (43) and (52) form a nonlinear equation system with the solution

x E e x E - G l c e = η E e η E - G l c e L E 1 (1) 1 ( K 1 + x G l c e ) K 1 x G l c e L E 1 .
(53)

Let us define that g 6 x E - G l c e and g 8 x E e . Similarly, the fractions of the two carrier state variables in the inner regions of the cell to the lumped state variable L E 3 can be computed from Eqs. (44) and (52) (Eq. (28))

x E - G l c i x E i = η E - G l c i η E i L E 3 (1) 1 ( x G l c i + K 2 ) x G l c i K 2 L E 3 .
(54)

We define that g 7 x E - G l c i and g 9 x E i . The fraction parameters for the G6P-state variables can be computed from Eqs. (45)-(46) and Eq. (28) with the solution

x E - G 6 P i x E - G l c - G 6 P i x G 6 P i = η E - G 6 P i η E - G l c - G 6 P i η G 6 P i L G 6 P (1) 1 ξ K 3 x E i K 4 x E - G l c i K 3 K 4 L G 6 P .
(55)

where ξ= K 3 K 4 + K 4 x E - G l c i + K 3 x E i . We define that g 3 x E - G 6 P i , g 4 x E - G l c - G 6 P i , and g 5 x G 6 P i . For the two original state variables that are kept as modified lumped state variables in the reduced model we define that g 1 x G l c e and g 2 x G l c i .

The fourth step of the method is to derive rate equations for the modified lumped state variables. Since the apparent conservations are separated into two disjoint clusters of fast state variables, we can treat the model for the inner and outer regions of the membrane separately. Let the modified lumped state variables corresponding to the outer region of the cell membrane be denoted by l m 1 = ( x G l c e L E 1 ) T , and the variables in the inner region of the cell membrane by l m 2 = ( x G l c i L E 3 ) T . Note that the state variable LG 6Pcan be replaced by a constant in the model, since L ° G 6 P =0. The ODEs of the modified lumped state variables are derived with Eq. (31). In the inner region the ODEs are

l ° m 1 = G l c e L ° E 1 = ( I + J ) - 1 l ° 1 = (1) = 1 + J 1 1 J 1 2 0 1 - 1 l ° 1 ,
(56)

where

l ° 1 = L ° G l c 1 L ° E 1 (1) - α ( x E G l c e - x E G l c i ) - α ( x E G l c e - x E G l c i ) - β ( x E e - x E i ) .
(57)

In the larger outer region of the cell membrane the ODEs take the form

l ° m 2 G l c i L ° E 3 = ( I + J ) - 1 l ° 2 = (1) = 1 + J 1 1 J 1 2 J 2 1 1 + J 2 2 - 1 l ° 2 ,
(58)

where

l ° 2 = L ° G l c 2 L ° E 2 (1) - α ( x E G l c i - x E G l c e ) - α ( x E G l c i - x E G l c e ) - β ( x E i - x E e ) .
(59)

In the fifth step of the method the reduced model, which is defined by the four ODEs in Eqs. (56)-(59), is simulated. The trajectories of the state variables of the reduced model can be back-translated to the original state variables with the fraction parameters defined in Eqs. (53), (54), and (55). The simulation results are shown in Figure 3 and Figure 4. All the state variables can be back-translated properly, which shows that the model properties that are important for recovery of the state variables are retained in the reduction. Implementations of the original model (Additional file 4), the reduced model (Additional file 9), and a script for simulation with SBtoolbox2 for MATLAB [36] (Additional file 10), are available in Additional files.

Figure 3
figure 3

Reduction with our method to four state variables. A comparison between the original glucose transport model and the model reduced to four state variables with our method, w.r.t. the state variables of the original model.

Figure 4
figure 4

Reduction with our method to four state variables. A comparison between the original glucose transport model and the model reduced to four state variables with our method, w.r.t. the state variables of the original model.

If we use L E 2 instead of LG 6Pas a modified lumped state variable, the reduced model will have the same state variables as the reduced model in [27, 37] (i.e., x G l c e , x G l c i , and x G 6 P i ) two additional state variables for the transporter. This gives a reduced model with five state variables, but equally many parameters as in the previous case. A comparison between the original model and the reduced model, w.r.t. the original state variables, is shown in Figure 5 and Figure 6. As can be seen the comparison is very good, in fact it is even slightly better than for the reduced model with four state variables. The details of the derivation of the reduced model are presented in Appendix A.5. Implementations of the original model (Additional file 4), the reduced model (Additional file 11), and a script for simulation with SBtoolbox2 for MATLAB [36] (Additional file 12), are available in Additional files.

Figure 5
figure 5

Reduction with our method to five state variables. A comparison between the original glucose transport model and the model reduced to five state variables with our method, w.r.t. the state variables of the original model.

Figure 6
figure 6

Reduction with our method to five state variables. A comparison between the original glucose transport model and the model reduced to five state variables with our method, w.r.t. the state variables of the original model.

Note that the only assumption used to derive the reduced model is that states that are involved in reactions at the membrane are in QSS. To investigate the parameter space region in which the QSS assumptions are valid we use the measure defined in Eq. (43). The maximal mean and infinity norm of the relative difference between the original and the reduced model in Eq. (43) over time is presented in Table 2. The reduced model appears to be relatively robust to changes in the parameters, although sensitive to small values of k-4 and to large values of k4. This is mainly due to that a large proportion of the transporter E is absorbed in x E - G 6 P i , which leads to that some of the QSS assumptions are invalid. We also observed that relative difference between the models is insensitive to the total concentration of the transporter for several orders of magnitude around the nominal value (L E = 0.01). However, note that this observation is specific to the studied model and may not be generalizable to other similar biochemical models.

Table 2 Robustness of the reduced model for large deviations from the nominal parameter point are presented for the glucose transport model, with a sampling frequency of 1 (starting from 1) time units.

Discussion

In this paper we have presented a novel method for reduction of biochemical models that is compatible with the concept of zooming. Several methods for reduction of biochemical models already exist in the literature. However, few of these methods result in biochemically interpretable models, and to our knowledge there are no nonlinear lumping methods for which the state variables and parameters of the reduced model can be back-translated (mapped) to the original model.

The application of the QSS assumption has been a commonly used tool in the modeling of biochemical networks since the late 1960s, and in chemical kinetics for more than 80 years [38]. The validity of the QSS approximation is well studied both for specific biochemical mechanisms [38, 39] and for more complex models [40, 41]. The resulting equations, together with conservation relations, are typically used to eliminate some of the state variables in the model (e.g., see [40]). However, with the examples in this paper we have showed that such an approach is not always sufficient, and we propose to use proper lumping of state variables in combination with back-translation.

Our method has several important advantages when applied to biochemical models. The most important advantage is that we end up with reduced models with a clear biological interpretation, meaning that each state variable of the original model corresponds to a fraction of exactly one of the state variables in the reduced model. A consequence is that neighboring species in the original model remain neighbors in the reduced model. Hence we can consider the original and reduced models as two different degrees of zooming; a concept that we discussed in some detail in [27] for linear models.

The work in this paper can be seen as an extension of the theory introduced for linear models in [27] to nonlinear models. The method is based on assumptions regarding the dynamics that result in a sufficient number of equations that are linear w.r.t. the state variables to be back-translated. Such equations are typically a natural result of QSS assumptions and conservations relations in models based on mass action kinetics [42], and in particular in models that involve transporters and enzymes (e.g., the models in this paper). However, note that our method may also be applicable to models with other types of reaction kinetics. We also note that if too few linear relations are available for calculation of fraction parameters for a part of a model, this part can still be reduced and the reduced model can be simulated, although we cannot back-translate the corresponding modified lumped state variables since no fraction parameters are available. However, depending on the purpose of the model it may be enough to calculate fraction parameters for a subset of the state variables in the reduced model. Linearization of the model around a steady state operating point may also be a feasible approach to calculate fraction parameters with the method in [27].

The proposed method enables mapping of the state variables and parameters of the reduced model to those of the original model. In [27] we referred to this mapping as back-translation. Back-translation is of great importance, since we can directly observe how modifications to the reduced model impact the original model. It also gives the modeler an opportunity to check whether the assumptions underlying the reduction are acceptable. To illustrate the power of back-translation we provide plots for comparison of simulations of the original and reduced models, w.r.t. the original state variables, for the models to which the method is applied in this paper.

Back-translation of state variables typically requires the solution of a system of nonlinear equations, which often results from the assumption of state variables in QSS and conservation relations. Unfortunately, analytic solutions to systems of nonlinear equations do in general not exist. An advantage with the proposed method is that such solutions are not required, since they are replaced by computation of the inverse of a matrix for each cluster of fast state variables, which is in general a more feasible task.

Our method was applied to a small model with five state variables that commonly appears as part of larger biochemical models, and to a previously published model for the transport of glucose in baker's yeast (S. cerevisiae) [37]. The first model was reduced from five to one state variable, and from five to three parameters. However, note that our focus has been on the reduction of the number of state variables and not the number of parameters, which are reduced as a side-effect of the QSS assumptions.

The model for glucose transport was first reduced with an approach in which the QSS equations and conservation relations were directly substituted into the remaining ODEs. The results of this approach are not satisfactory since the reduced model gives predictions that are different from the original model for most state variables. Our method was then applied to the same model both together with the assumption of equal concentration of the transporter in the inner and outer regions of the cell membrane used in [37], and without any additional assumptions. The application of our method together with the assumptions used in [37] results in a model with three state variables. The state dynamics is significantly better preserved than with the first approach, although still not satisfactory. We then decided to reduce the original model without the assumption regarding the localization of the transporter, with two different definitions of modified lumped state variables. While one of these definitions results in a reduced model with four state variables and gave rather accurate predictions, the other choice reduces the number of state variables to five and gives an excellent description of the state dynamics. It is therefore apparent that there is a tradeoff between accuracy and the number of state variables in the reduction process. The glucose transport model corresponds to the first part of glycolysis, in which glucose is transported into the cell. We therefore propose that it might be rewarding to carefully re-investigate the assumptions underlying the reaction rate equations in complete models of glycolysis (see [43] for one example).

We have also observed a few issues regarding the implementation of the method. The symbolic inversion of the matrix that is necessary to compute the dynamics of the modified lumped state variables may be expensive. However, this is typically only a practical limitation for large matrices, which result from large clusters of fast state variables. In our experience large clusters of fast state variables are relatively rare also in large biochemical models. Another option, if it is not practically feasible to invert the symbolic matrix, is to solve the system of linear equations in Eq. (31) numerically. We also observed that the symbolic right-hand side of the resulting differential equations may be long. However, these are usually not practical limitations for the applicability of the method, e.g., the simulations of all examples in this paper are very fast on a modern computer. Available methods to reduce the analytic reaction rate expressions include sensitivity analysis w.r.t. state variables and parameters, and the method proposed in [17].

There is still no consensus method for automatic identification of state variables in QSS, although criteria for the detection of state variables in QSS have been proposed, for example in [30]. A simple approach is to simulate the original model and investigate for which state variables the corresponding in- and outflow reaction rates are approximately equal. State variables for which this condition holds are then considered to be in QSS. Note that for the models in this paper it was already clear from the biochemical understanding of the corresponding systems which of the state variables that could be considered fast (see [37] for the glucose transport example). However, an appropriate general criterion for automatic identification of state variables in QSS is still lacking.

Although the theory presented in this paper constitutes a great leap forward for construction of zoomable models, more research is required to make the method fully automatic. An important challenge is to define a meaningful measure for the similarity between the hierarchical model layers (degrees of zooming). Another interesting, although trivial, observation that deserves further attention is that QSS assumptions typically do not hold in the whole parameter space. Although the reduced models in this paper appear to be robust to varying parameter values it may not be the case in general. It may therefore be revealing to compare the original model and the reduced model to characterize the parameter space regions in which the QSS-assumptions are valid.

Conclusions

We have presented a novel method for reduction of biochemical models that is compatible with the concept of zooming. Zooming allows the modeler to operate on different levels of model granularity, and enables a direct interpretation of how modifications to the model on one level affect the same model on other levels in the hierarchy. The proposed method is based on the application of proper lumping in combination with the identification of linear relations in nonlinear equations.

The method was applied to two example models. The first model is small and commonly occurring as a part of larger biochemical models. The second example is a model for glucose transport in baker's yeast, which constitutes the starting point for glycolysis. Both models could be significantly reduced with the proposed method, and the resulting state variables could be back-translated to the original state variables. The method that is presented in this paper constitutes an extension of the method that was previously developed for linear biochemical models to its nonlinear counterpart. Since most models in the systems biology community are in fact nonlinear, our method constitutes an important step towards zoomable biochemical models.

A Appendix

A.1 Appendix 1

The ordinary differential equations for the enzyme kinetics model take the form

= Ė C S ° C P ° = - k 1 S E + k - 1 C S - k 1 S E + k - 1 C S + k 3 C P - k - 3 P E k 3 C P - k - 3 P E k 1 S E - k - 1 C S - k 2 C S k 2 C S - k 3 C P + k - 3 P E = = - 1 0 0 - 1 0 1 0 0 1 1 - 1 0 0 1 - 1 r 1 r 2 r 3 = S r ( x , p ) ,

where r1 = k1SE - k-1C S , r2 = k2C S , and r3 = k3C P - k-3PE.

The parameters are set to values that satisfy the assumptions of dominating and insignificant reaction terms, with k1 = 1000, k-1 = 2000, k2 = 1, k3 = 3000, and k-3 = 1000, together with the initial conditions S(0) = E(0) = 1 and P(0) = C S (0) = C P (0) = 0. This gives the parameter values M1 = 0.5, M3 = 3 and L E = 1 in the reduced model. Initial conditions can in general be obtained from a short simulation of the original model, until the fast state variables reach QSS, but in this case Eq. (19) gives an analytic expression of S(0) (P(0) = 0)

S ( 0 ) = 1 2 ( L S - L E - M 1 - 1 + ( L S + L E + M 1 - 1 ) 2 - 4 L S L E ) = = 3 - 1 ,

where M 1 - 1 = K 1 , S = A, L S = L1, and L E = L2.

A.2 Appendix 2

The ordinary differential equations for the glucose transport model take the form

d x G l c e d t = - k 1 x E e x G l c e + k - 1 x E - G l c e , d x G l c i d t = - k 2 x E i x G l c i + k - 2 x E - G l c i , d x E - G 6 P i d t = k 4 x E i x G 6 P i - k - 4 x E - G 6 P i , d x E - G l c - G 6 P i d t = k 3 x E - G l c i x G 6 P i - k - 3 x E - G l c - G 6 P i , d x G 6 P i d t = - k 3 x E - G l c i x G 6 P i + k - 3 x E - G l c - G 6 P i - k 4 x E i x G 6 P i + k - 4 x E - G 6 P i , d x E - G l c e d t = α ( x E - G l c i - x E - G l c e ) + k 1 x E e x G l c e - k - 1 x E - G l c e , d x E - G l c i d t = α ( x E - G l c e - x E - G l c i ) - k 3 x E - G l c i x G 6 P i + k - 3 x E - G l c - G 6 P i + k 2 x E i x G l c i - k - 2 x E - G l c i , d x E e d t = β ( x E i - x E e ) - k 1 x E e x G l c e + k - 1 x E - G l c e , d x E i d t = β ( x E e - x E i ) - k 4 x E i x G 6 P i + k - 4 x E - G 6 P i - k 2 x E i x G l c i + k - 2 x E - G l c i ,

which was introduced in [27].

A.3 Appendix 3

In this section we investigate an alternative (naive) approach to reduce the glucose transport model. The first step is to identify state variables for which the QSS assumption holds, and the mass conservation relations in the model. In the second step of this approach we then substitute the corresponding system of equations into the ODEs corresponding to slow state variables.

Now consider the model for glucose transport in yeast. We assume that the state variables x G l c e , x G l c i , x E - G 6 P i and x E - G l c - G 6 P i are in QSS, which gives Eqs. (43)-(46). Note that Eqs. (45)-(46) indirectly imply that x G 6 P i is in steady state. The substitution of Eqs. (43) - (46) into the ODEs of the original model gives

d x E - G l c e d t α ( x E - G l c i - x E - G l c e ) ,
(61)
d x E - G l c i d t α ( x E - G l c e - x E - G l c i ) ,
(62)
d x E e d t β ( x E i - x E e ) ,
(63)
d x E i d t β ( x E e - x E i ) .
(64)

Note that the state variables x E - G l c e and x E - G l c i are decoupled from the state variables x E e and x E i in Eqs. (61)-(64).

There are three molecules (moieties) whose mass is conserved in the model as a whole, i.e., Glc, G6P, and E. However, we can not substitute any of the conservation relations into the remaining ODEs without re-introducing state variables that were already eliminated. So the final reduced model takes the form of Eqs. (61) - (64). However, the sum of the state variables x E - G l c e and x E - G l c i , and x E e and x E i is conserved in the reduced model, which makes it possible to reduce the model to two state variables.

Unfortunately, due to the form of the ODEs and the initial conditions of the state variables in the reduced model, the state variables x E - G l c e and x E - G l c i remain equal to zero at all times, and only the state variables x E e and x E i take non-zero values. We therefore decided to simulate the original model for a short time until the fast state variables reach QSS, and to use the final state variable values as initial conditions in the reduced model.

The solution to the equation system defined by Eqs. (43)-(46) and Eq. (49) is

x G l c e K 1 x E - G l c e x E e , x G l c i K 2 x E - G l c i x E i , x G 6 P i K 3 K 4 K 4 x E - G l c i + K 3 x E i + K 3 K 4 L G 6 P , x E - G 6 P i K 3 x E i K 4 x E - G l c i + K 3 x E i + K 3 K 4 L G 6 P , x E - G l c - G 6 P i K 4 x E - G l c i K 4 x E - G l c i + K 3 x E i + K 3 K 4 L G 6 P ,

which can be used for back-translation of the state variables of the reduced model to those of the original model.

The predictions of the state variables of the original model, resulting from simulations of the original model and the reduced model with the parameter values set as in [37], is not satisfactory for any other state variable than G 6P, which remains approximately constant over time. Implementations of the original model and the reduced model in SBtoolbox2 for MATLAB [36] are included in Additional files.

A.4 Appendix 4

In this section we apply our method to the glucose transport model, and following [37] we will assume that the concentrations of the transporter are constant in the inner and outer regions of the cellular membrane. With this assumption the distribution among the transporter state variables of the original model, which constitute the lumped state variables L E , is uniquely defined.

The first step of the method is to identify the apparent conservation relations in the model. We note that G6P and the transporter E are conserved, and apparent conserved glucose (see Definition 1) in the inner and outer regions of the membrane, respectively. The four apparent conservation relations take the form

l = L G l c 1 L G l c 2 L G 6 P L E = M x = = 1 0 0 0 0 1 0 0 0 0 1 0 1 0 0 1 0 0 0 0 1 1 1 0 0 0 0 0 0 1 1 0 1 1 1 1 x ,

where x= ( x G l c e x G l c i x E - G 6 P i x E - G l c - G 6 P i x G 6 P i x E - G l c e x E - G l c i x E e x E i ) T .

In the second step of the method we define the modified lumped state variables. We decide to keep L E in the reduced model since it corresponds to an exact conservation, and therefore results in the largest reduction possible (note that the exact conservation relations, L E and LG 6P, can not simultaneously be used since the lumping would then not be proper). The modified lumped state variables take the form

l m = x G l c e x G l c i x G 6 P i L E = M m x = = 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 1 0 1 1 1 1 x .

We note that Eqs. (43)-(47), and (50) are all linear w.r.t. the state variables that constitute state L E , so the requirement for the existence of at least 6 (n m ) linear relations is satisfied, which enables back-translation in step three of the method.

In the third step of the method we derive the fraction parameters for the lumped state variable L E . Eqs. (43)-(47), and (50) form an equation system, corresponding to Eq. (28)

A ( l m , p ) x m k = 0 0 1 0 - x G l c e 0 0 0 0 1 0 - x G l c i 1 0 0 0 0 - x G 6 P i 0 1 0 - x G 6 P i 0 0 0 0 α - α β - β 1 1 1 1 1 1 (1) x E - G 6 P i x E - G l c - G 6 P i x E - G l c e x E - G l c i x E e x E i 0 0 0 0 0 x L = b k ( l m , p ) , (2)
(65)

where x G l c e = k 1 k - 1 x G l c e , x G l c i = k 2 k - 2 x G l c i , x G 6 P i = k 2 k - 2 x G 6 P i and x G 6 P i = k 4 k - 4 x G 6 P i . The solution to Eq. (65) is given by Eq. (29)

x E - G 6 P i x E - G l c - G 6 P i x E - G l c e x E - G l c i x E e x E i = η E - G 6 P i η E - G l c - G 6 P i η E - G l c e η E - G l c i η E e η E i L E (1) x G 6 P i ( β + α x G l c e ) ζ x G 6 P i x G l c i p ( β + α x G l c e ) ζ x G l c e ( β + α x G l c i ) ( β + α x G l c e ) ζ x G l c i ( β + α x G l c e ) ζ ( β + α x G l c i ) ( β + α x G l c e ) ζ ( β + α x G l c e ) ζ L E (2)
(66)

where

ζ = ( β x G l c e + 2 x G l c e α x G l c i + 2 β + α x G l c i + α x G l c e + β x G l c i x G 6 P i + + α x G l c i x G 6 P i x G l c e + β x G l c i + β x G 6 P i + α x G l c e x G 6 P i )
(67)

and where the fraction parameters were calculated with Eq. (30). We note that the fraction parameters are functions of x G l c e , x G l c i , and x G 6 P i , which are state variables both in the original- and in the reduced model. In the fourth step of the method we derive differential equations for the modified lumped state variables. The ODE for the fourth state is l ° m 4 = L ° E =0, which is replace by a constant. The ODEs for the other states are

l ° m 1 : 3 = G l c e G l c i G 6 P i = 1 + J 1 1 J 1 2 J 1 3 J 2 1 1 + J 2 2 J 2 3 J 3 1 J 3 2 1 + J 3 3 - 1 l ° 1 : 3

where l m 1 : 3 denotes the first three states variables in l m . Note that there are three state variables in the reduced model, which is the same number as for the reduced model in [37].

In the fifth step of our method we compare predictions of the original state variables between the original and reduced models, where L E is back-translated with the fraction parameters defined in Eqs. (66)-(67).

The simulation results are clearly more accurate than with the approach in Appendix A.3, although still not satifying. We refer to Additional files for implementations of the original and reduced models in SBtoolbox2 for MATLAB [36].

A.5 Appendix 5

In this section we apply our method to the glucose transport model, but with an alternative definition of the modified lumped state variables. We do not use the assumption of constant regional concentrations of the transporter (Eq. (50)).

In the first step of the method we note that the apparent conservations are given by Eq. (51).

In the second step of our method we decide to keep state variable L E 2 , instead of LG 6P, in the reduced model. This leads to the following definition of the modified lumped state variables

l m = x G l c e L E 1 x G l c i x G 6 P i L E 2 = M m x = (1) = 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 1 0 0 1 0 1 x ,
(68)

Where x= ( x G l c e x G l c i x E - G 6 P i x E - G l c - G 6 P i x G 6 P i x E - G l c e x E - G l c i x E e x E i ) T , and we note that x G l c e , x G l c i , and x G 6 P i are state variables both in the original and reduced models. Also note that the requirement of at least n m equations, that are linear w.r.t. the original state variables and linearly independent, is satified for each of the modified lumped state variables by Eqs. (43)-(49).

In the third step of the method we calculate fraction parameters for the modified lumped state variables L E 1 and L E 2 , which correspond to more than one of the original state variables. The fraction parameters for state variable L E 1 are given by Eq. (53). We can now use Eqs. (44)-(46) and Eq. (68) to form an equation system corresponding to Eq. (28)

A ( l m , p ) x m k = - x G l c i K 2 0 0 - x G 6 P i 0 K 4 0 0 - x G 6 P i 0 K 3 1 1 1 1 x E i x E - G l c i x E - G 6 P i x E - G l c - G 6 P i 0 0 0 L E 2 = b k ( l m , p ) ,

with the solution given by Eq. (29)

x E i x E - G l c i x E - G 6 P i x E - G l c - G 6 P i = η E i η E - G l c i η E - G 6 P i η E - G l c - G 6 P i L E 2 (1) K 2 K 3 K 4 ξ K 3 K 4 x G l c i ξ K 2 K 3 x G 6 P i ξ K 4 x G 6 P i x G l c i ξ L E 2 ,
(69)

where

ξ= x G 6 P i ( K 4 x G l c i + K 2 K 3 ) + K 3 K 4 x G l c i + K 2 K 3 K 4

and the fraction parameters were calculated with Eq. (30).

The fourth step of the method is to derive ODEs for the modified lumped state variables. Since the apparent conservations are separated into two disjoint clusters of fast state variables, we can treat the model for the inner and outer regions of the membrane separately. The rate equations for the outer region are given by Eqs. (56)-(57). Eq. (31) gives us the ODEs of the modified lumped state variables in the inner region

l ° m 3 : 5 = G l c i G 6 P i L ° E 2 = ( I + J ) - 1 l ° = (1) = 1 + J 3 3 J 3 4 J 3 5 J 4 3 1 + J 4 4 J 4 5 0 0 1 - 1 l ° 3 : 5 ,
(70)

where l3:5 and l m 3 : 5 are the last three state variables of l and l m , respectively, and

l ° 3 : 5 = L ° G l c 2 L ° G 6 P L ° E 2 (1) - α ( x E - G l c i - x E - G l c e ) 0 - α ( x E - G l c i - x E - G l c e ) - β ( x E i - x E e ) ,
(71)

In the fifth step of the method we simulate the reduced model with Eqs. (56)-(57) and (70)-(71) and we then use Eqs. (53) and (69) for back-translation of the state variables. A comparison between the original model and the reduced model, w.r.t. the state variables of the original model, is presented in Figure 5 and Figure 6. The agreement between the models is very good. We refer to Additional files for implementations of the original and reduced models in SBtoolbox2 for MATLAB [36].

Additional files

The original and reduced versions of the models presented in this paper, and scripts for simulation and comparison between the original and reduced versions of the models. Note that the systems biology toolbox for MATLAB [36] and the symbolic math toolbox for MATLAB must be installed on the system for simulation of the attached models.

References

  1. Kitano H: Computational systems biology. Nature. 2002, 420 (6912): 206-10. 10.1038/nature01254

    Article  CAS  PubMed  Google Scholar 

  2. Janes K, Yaffe M: Data-driven modelling of signal-transduction networks. Nat Rev Mol Cell Biol. 2006, 7: 820-828. 10.1038/nrm2041

    Article  CAS  PubMed  Google Scholar 

  3. Aldridge B, Burke J, Lauffenburger D, Sorger P: Physicochemical modelling of cell signalling pathways. Nat Cell Biol. 2006, 8: 1195-1203. 10.1038/ncb1497

    Article  CAS  PubMed  Google Scholar 

  4. Cedersund G, Roll J: Systems biology: Methods for evaluation and comparison of competing explanations for given biological data. FEBS J. 2009, 276: 903-922. 10.1111/j.1742-4658.2008.06845.x

    Article  CAS  PubMed  Google Scholar 

  5. Becker V, Schilling M, Bachmann J, Baumann U, Raue A, Maiwald T, Timmer J, Klingmüller U: Covering a Broad Dynamic Range: Information Processing at the Erythropoietin Receptor. Science. 2010, 328: 1404-1408. 10.1126/science.1184913

    Article  CAS  PubMed  Google Scholar 

  6. Brännmark C, Palmer R, Glad T, Cedersund G, Stralfors P: Mass and information feedbacks through receptor endocytosis govern insulin signaling as revealed using a parameter-free modeling framework. J Biol Chem. 2010, 26: 20171-20179.

    Article  Google Scholar 

  7. Zhou K, Doyle J, Glover K: Robust and Optimal Control. 1996, Prentice Hall,

    Google Scholar 

  8. Ljung L, Glad T: Control theory - multivariable and nonlinear methods. 2000, Taylor and Francis,

    Google Scholar 

  9. Eriksson O: Simplicity within Complexity - Understanding Dynamics of Cellular Networks by Model Reduction. PhD thesis. 2008, Stockholm University,

    Google Scholar 

  10. Cedersund G: Core-box modelling. PhD thesis. 2006, Chalmers, Gothenburg, Sweden,

    Google Scholar 

  11. Hahn J, Edgar T: An improved method for nonlinear model reduction using balancing of empirical gramians. Computers Chem Eng. 2002, 26: 1379-1397. 10.1016/S0098-1354(02)00120-5.

    Article  CAS  Google Scholar 

  12. Liebermeister W, Baur U, Klipp E: Biochemical network models simplified by balanced truncation. FEBS J. 2005, 272: 4034-4043. 10.1111/j.1742-4658.2005.04780.x

    Article  CAS  PubMed  Google Scholar 

  13. Sandberg H, Murray R: Model reduction of interconnected linear systems. Optimal Control, Applications and Methods, Special Issue on Directions, Applications, and Methods in Robust Control. 2008,

    Google Scholar 

  14. Danø S, Madsen M, Schmidt H, Cedersund G: Reduction of a biochemical model with preservation of its basic dynamic properties. FEBS J. 2006, 273: 4862-77. 10.1111/j.1742-4658.2006.05485.x

    Article  PubMed  Google Scholar 

  15. Jacobsen E, Cedersund G: Structural robustness of biochemical network models-with application to the oscillatory metabolism of activated neutrophils. IET Syst Biol. 2008, 2: 39-47. 10.1049/iet-syb:20070008

    Article  CAS  PubMed  Google Scholar 

  16. AnalogInsydes. http://www.wolfram.com/products/applications/insydes/

  17. Schmidt H, Madsen M, Danø S, Cedersund G: Complexity reduction of biochemical rate expressions. Bioinformatics. 2008, 24: 848-54. 10.1093/bioinformatics/btn035

    Article  CAS  PubMed  Google Scholar 

  18. 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, Priami, Berlin: Springer-Verlag,

    Google Scholar 

  19. Scherpen J: Model reduction for nonlinear control systems. DISC model reduction coure notes. 2005,

    Google Scholar 

  20. Okino M, Mavrovouniotis M: Simplification of mathematical models of chemical reaction systems. Chemical Reviews. 1998, 98: 391-408. 10.1021/cr950223l

    Article  CAS  PubMed  Google Scholar 

  21. Gorban A, Karlin I: Method for invariant manifold for chemical kinetics. Chem Eng Sci. 2003, 58: 4751-4768. 10.1016/j.ces.2002.12.001.

    Article  CAS  Google Scholar 

  22. Toth J, Li G, Rabitz H, Tomlin A: The Effect of Lumping and Expanding on Kinetic Differential Equations. SIAM J Appl Math. 1997, 57: 1531-1556. 10.1137/S0036139995293294.

    Article  Google Scholar 

  23. Conzelmann H, Saez-Rodriguez J, Sauter T, Bullinger E, Allgöwer F, Gilles ED: Reduction of mathematical models of signal transduction networks: simulation-based approach applied to EGF receptor signalling. Syst Biol (Stevenage). 2004, 1: 159-169. 10.1049/sb:20045011.

    Article  CAS  Google Scholar 

  24. Conzelmann H, Saez-Rodriguez J, Sauter T, Kholodenko B, Gilles E: A domain-oriented approach to the reduction of combinatorial complexity in signal transduction networks. BMC Bioinformatics. 2006, 7: 34- 10.1186/1471-2105-7-34

    Article  PubMed Central  PubMed  Google Scholar 

  25. Koschorreck M, Conzelmann H, Ebert S, Ederer M, Gilles E: Reduced modeling of signal transduction -a modular approach. BMC Bioinformatics. 2007, 8: 336- 10.1186/1471-2105-8-336

    Article  PubMed Central  PubMed  Google Scholar 

  26. Conzelmann H, Fey D, Gilles E: Exact model reduction of combinatorial reaction networks. BMC Syst Biol. 2008, 2: 78-

    Article  PubMed Central  PubMed  Google Scholar 

  27. Sunnåker M, Schmidt H, Jirstrand M, Cedersund G: Zooming of states and parameters using a lumping approach including back-translation. BMC Syst Biol. 2010, 4: 28- 10.1186/1752-0509-4-28

    Article  PubMed Central  PubMed  Google Scholar 

  28. Wei J, Kuo J: A lumping analysis in monomolecular reaction systems: Analysis of the exactly lumpable system. Ind Eng Chem Fund. 1969, 8: 114-123. 10.1021/i160029a019.

    Article  CAS  Google Scholar 

  29. Gorban AN, Radulescu O: Dynamical robustness of biological networks with hierarchical distribution of time scales. IET Syst Biol. 2007, 1 (4): 238-246. 10.1049/iet-syb:20060083

    Article  CAS  PubMed  Google Scholar 

  30. Radulescu O, Gorban AN, Zinovyev A, Lilienbaum A: Robust simplifications of multiscale biochemical networks. BMC Syst Biol. 2008, 2: 86- 10.1186/1752-0509-2-86

    Article  PubMed Central  PubMed  Google Scholar 

  31. Gorban A, Radulescu O: Dynamic and Static Limitation in Reaction Networks, revisited. Adv Chem Eng. 2008, 34: 103-173.

    Article  CAS  Google Scholar 

  32. Dokoumetzidis A, Aarons L: Proper lumping in systems biology models. IET Syst Biol. 2009, 3: 40-51. 10.1049/iet-syb:20070055

    Article  CAS  PubMed  Google Scholar 

  33. Vredenberg W, Bulychev A: Photo-electrochemical control of photosystem II chlorophyll fluorescence in vivo. Bioelectrochemistry. 2002, 57: 123-128. 10.1016/S1567-5394(02)00062-2

    Article  CAS  PubMed  Google Scholar 

  34. Battogtokh D, Tyson JJ: Bifurcation analysis of a model of the budding yeast cell cycle. Chaos. 2004, 14 (3): 653-661. 10.1063/1.1780011

    Article  CAS  PubMed  Google Scholar 

  35. Chen KC, Calzone L, Csikasz-Nagy A, Cross FR, Novak B, Tyson JJ: Integrative analysis of cell cycle control in budding yeast. Mol Biol Cell. 2004, 15 (8): 3841-3862. 10.1091/mbc.E03-11-0794

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  36. Schmidt H, Jirstrand M: Systems Biology Toolbox for MATLAB: a computational platform for research in systems biology. Bioinformatics. 2006, 22 (4): 514-515. 10.1093/bioinformatics/bti799

    Article  CAS  PubMed  Google Scholar 

  37. Rizzi M, Theobald U, Querfurth E, Rohrhirsch T, Baltes M, Reuss M: In vivo investigations of glucose transport in Saccharomyces cerevisiae. Biotechnol Bioeng. 1996, 49 (3): 316-327.

    Article  CAS  PubMed  Google Scholar 

  38. Flach EH, Schnell S: Use and abuse of the quasi-steady-state approximation. Syst Biol (Stevenage). 2006, 153 (4): 187-191.

    Article  CAS  Google Scholar 

  39. Tzafriri AR, Edelman ER: The total quasi-steady-state approximation is valid for reversible enzyme kinetics. J Theor Biol. 2004, 226 (3): 303-313. 10.1016/j.jtbi.2003.09.006

    Article  CAS  PubMed  Google Scholar 

  40. Millat T, Bullinger E, Rohwer J, Wolkenhauer O: Approximations and their consequences for dynamic modelling of signal transduction pathways. Math Biosci. 2007, 207: 40-57. 10.1016/j.mbs.2006.08.012

    Article  PubMed  Google Scholar 

  41. Stoleriu I, Davidson FA, Liu JL: Effects of periodic input on the quasi-steady state assumptions for enzyme-catalysed reactions. J Math Biol. 2005, 50 (2): 115-132. 10.1007/s00285-004-0282-6

    Article  CAS  PubMed  Google Scholar 

  42. Schoeberl B, Eichler-Jonsson C, Gilles ED, Müller G: Computational modeling of the dynamics of the MAP kinase cascade activated by surface and internalized EGF receptors. Nat Biotechnol. 2002, 20 (4): 370-375. 10.1038/nbt0402-370

    Article  PubMed  Google Scholar 

  43. Teusink B, Passarge J, Reijenga CA, Esgalhado E, van der Weijden CC, Schepper M, Walsh MC, Bakker BM, van Dam K, Westerhoff HV, Snoep JL: Can yeast glycolysis be understood in terms of in vitro kinetics of the constituent enzymes? Testing biochemistry. Eur J Biochem. 2000, 267 (17): 5313-5329.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

This work was supported by grants from the European Commission 6th Framework Programme (BIOSIM, grant No 005137) and 7th Framework Programme (UNICELLSYS, grant No 201142), SystemsX.ch (RTD project YeastX), the Swedish Foundation for Strategic Research through the Gothenburg Mathematical Modelling Centre, the Swedish Research Council, and by Lions. We would also like to thank John Johansson for his contributions to the development of the method in his M.Sc. thesis at the Fraunhofer-Chalmers Centre.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Mikael Sunnåker.

Additional information

Authors' contributions

MS developed the theory, most of the aspects of the method, and did the calculations within this project, which was managed by MJ. MS wrote the major part of the paper, with contributions from GC and MJ. All authors read and approved the final manuscript.

Electronic supplementary material

Additional file 1:Model 1. The original enzyme kinetics model. (TXT 2 KB)

Additional file 2:Model 2. The reduced enzyme kinetics model. (TXT 3 KB)

12918_2011_751_MOESM3_ESM.TXT

Additional file 3:Script 1. Script for comparison between the original enzyme kinetics model and the reduced model. (TXT 1 KB)

Additional file 4:Model 3. The original glucose transport model. (TXT 6 KB)

Additional file 5:Model 4. The reduced glucose transport model with the alternative (naive) approach. (TXT 2 KB)

12918_2011_751_MOESM6_ESM.TXT

Additional file 6:Script 2. Script for comparison between the original glucose transport model and the reduced model with the alternative (naive) approach. (TXT 3 KB)

12918_2011_751_MOESM7_ESM.TXT

Additional file 7:Model 5. The reduced glucose transport model with our method and the assumption of constant concentrations of the transporter in the inner and outer regions of the cellular membrane. (TXT 2 MB)

12918_2011_751_MOESM8_ESM.TXT

Additional file 8:Script 3. Script for comparison between the original glucose transport model and the reduced model with our method and the assumption of constant concentrations of the transporter in the inner and outer regions of the cellular membrane. (TXT 3 KB)

Additional file 9:Model 6. The reduced glucose transport model with four state variables with our method. (TXT 151 KB)

12918_2011_751_MOESM10_ESM.TXT

Additional file 10:Script 4. Script for comparison between the original glucose transport model and the reduced model with four state variables with our method. (TXT 3 KB)

Additional file 11:Model 7. The reduced glucose transport model with five state variables with our method. (TXT 115 KB)

12918_2011_751_MOESM12_ESM.TXT

Additional file 12:Script 5. Script for comparison between the original glucose transport model and the reduced model with five state variables with our method. (TXT 3 KB)

Authors’ original submitted files for images

Rights and permissions

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.

Reprints and permissions

About this article

Cite this article

Sunnåker, M., Cedersund, G. & Jirstrand, M. A method for zooming of nonlinear models of biochemical systems. BMC Syst Biol 5, 140 (2011). https://doi.org/10.1186/1752-0509-5-140

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/1752-0509-5-140

Keywords