# Robust simplifications of multiscale biochemical networks

- Ovidiu Radulescu
^{1}Email author, - Alexander N Gorban
^{2, 6}, - Andrei Zinovyev
^{3, 4, 5, 6}and - Alain Lilienbaum
^{7, 8}

**2**:86

https://doi.org/10.1186/1752-0509-2-86

© Radulescu et al; licensee BioMed Central Ltd. 2008

**Received: **12 April 2008

**Accepted: **14 October 2008

**Published: **14 October 2008

## Abstract

### Background

Cellular processes such as metabolism, decision making in development and differentiation, signalling, etc., can be modeled as large networks of biochemical reactions. In order to understand the functioning of these systems, there is a strong need for general model reduction techniques allowing to simplify models without loosing their main properties. In systems biology we also need to compare models or to couple them as parts of larger models. In these situations reduction to a common level of complexity is needed.

### Results

We propose a systematic treatment of model reduction of multiscale biochemical networks. First, we consider linear kinetic models, which appear as "pseudo-monomolecular" subsystems of multiscale nonlinear reaction networks. For such linear models, we propose a reduction algorithm which is based on a generalized theory of the limiting step that we have developed in [1]. Second, for non-linear systems we develop an algorithm based on dominant solutions of quasi-stationarity equations. For oscillating systems, quasi-stationarity and averaging are combined to eliminate time scales much faster and much slower than the period of the oscillations. In all cases, we obtain robust simplifications and also identify the critical parameters of the model. The methods are demonstrated for simple examples and for a more complex model of NF-*κ* B pathway.

### Conclusion

Our approach allows critical parameter identification and produces hierarchies of models. Hierarchical modeling is important in "middle-out" approaches when there is need to zoom in and out several levels of complexity. Critical parameter identification is an important issue in systems biology with potential applications to biological control and therapeutics. Our approach also deals naturally with the presence of multiple time scales, which is a general property of systems biology models.

## Keywords

## Background

Model reduction techniques are used to reduce the dimensionality of complex dynamics. Applications of model reduction techniques in chemical engineering (coarse graining in phase transitions, reactors, combustion [2–8]), in ecology [9] or climatology, are well developed. A collection of reviews in model reduction for kinetic problems can be found in [10]. In systems biology, ad hoc reduction methods have been applied to signal transduction [11] and to clocks [12, 13]. Combinatorial complexity of receptors and scaffolds can be reduced by exact lumping [14, 15].

We may distinguish among three classes of model reduction techniques. *Trajectory based techniques* use the integration of the dynamical equations and look for a small number of reduced variables [16]. The empirical orthogonal eigenfunctions (EOF), also called Proper Orthogonal Decomposition (POD), or Karhunen-Loève expansion (KL) method, consists in finding a low dimension linear (flat) manifold, containing (or sufficiently close to) the trajectories [17, 18]. *Singular perturbations techniques* eliminate fast variables whose dynamics is slaved by the slower variables. The Computational Singular Perturbation (CSP) method provides approximations of a low dimensional invariant manifold, containing the dynamics [2, 3]. Invariant manifolds can be calculated by various other methods [4–8]. Slow-fast or more general master-slave splittings (splittings with no feed-back) were discussed by [19, 20]. Chemical enzymatic kinetics beyond quasi-stationarity and quasi-equilibrium has been studied in [21]. Averaging has been used to eliminate rapid oscillations of microscopic degrees of freedom and to obtain smaller models [22–24]. *Aggregation or lumping techniques* have been proposed by many authors [9, 14, 15, 25]. Reaction graph contraction methods such as Clarke's [26] replace the reactions mechanism by simpler mechanisms in which some intermediate species are absent.

Normally, identification of two well separated time scales is enough to reduce the system by using slow/fast decompositions [20]. However, the biochemical networks used to model cell physiology are multiscale, i.e. they have many, well separated time scales. For example, changing gene expression programs can take hours and even days while protein complex formation goes on the second scale and post-translational protein modifications take minutes to happen. Protein life half-times can vary from minutes to days. This important observation applies not only to time scales but also to concentration values of various species in these networks. mRNA copy numbers can change from some units to tens of thousands, and the dynamic concentration range of biological proteins can reach up to five orders of magnitude.

The aim of our paper is to propose model reduction methods well adapted to this situation. The mathematical techniques that we use (limitation, averaging, quasi-stationarity) have a long history. However, their combination into practical recipes that we propose is original and well adapted for the study of multiscale biochemical networks. Our most important development is the concept of dominant subsystem (that we also call limit simplification).

The idea of dominant subsystems in asymptotic analysis of dynamical systems is due to Newton and developed by Kruskal [27]. There are several ways to obtain dominant subsystems. These can be leading terms in power expansions of small parameters. Thus, multiscale expansions are standard techniques in perturbation theory [28]. Asymptotic theories using powers of small parameters were applied to study spectral properties of multiscale matrices [27, 29–31]. In [1] we have proposed a different approach to dominant subsystems. This approach exploits the reaction network structure to select dominant pathways and to obtain simplified reaction mechanisms. The simplifications are robust because are valid for a large range of parameters.

Understanding the functioning of large networks of biochemical reactions could rely on having a hierarchy of such simplifications, ie a set of models that can be obtained one from another by model reduction. Molecular networks are designed to fulfill many simple tasks. For each one of this tasks, the system scans only a small part of its high dimensional phase space. Geometrically speaking, it evolves on a stable low dimensional invariant manifold with branching in the fast directions [5]. Changing tasks, the system can jump from one stable branch of the manifold to another one. These represent jumps from one simplification (dominant subsystem) to another one. Finding the set of simplifications of a molecular network means providing the set of functioning modes for the network.

Thus, dominant subsystems provide an answer to a very practical question: how to describe the dynamics of a multiscale network? During almost all time this could be simplified and the system behaves as a small one. Our methods show how to obtain the small dominant subsystem from the topology of the network and from the orders of magnitude of kinetic constants and species concentrations. In multiscale systems, concentration orders can change dynamically and the small system may change at discrete times. The whole system walks along small subsystems. The discrete dynamics of this walk supplements the dynamics of individual small subsystems.

Dominant subsystem can be used to answer another important question: given a network model, which are its critical parameters? Many of the parameters of the initial model are no longer present in the dominant subsystem: these parameters are non-critical. Parameters of dominant subsystems indicate putative targets to change the behavior of the large network.

Finally, dominant subsystems can be used to compare models. Systems biology model repositories contain models of various degree of complexity. To be compared, or to be integrated into larger ones, models must be simplified to a common level of complexity.

Our methods perform well when we have total or partial separation of time and/or concentration scales. Total separation of time scales means the following: picking two timescales at random *τ*_{
i
}, *τ*_{
j
}one has either *τ*_{
i
}<<*τ*_{
j
}or *τ*_{
i
}>> *τ*_{
j
}with probability close to one. It is easy to construct a totally separated linear system. Choose constants of biochemical reactions independently and distributed uniformly over a large interval in logarithmic scale: picking two timescales at random *τ*_{
i
}, *τ*_{
j
}one has either *τ*_{
i
}<<*τ*_{
j
}or *τ*_{
i
}>> *τ*_{
j
}with probability close to one. This situation has been studied in detail in [1]. Though, it is difficult to have total separation in non-linear systems. For these, even if reactions constants are independent, timescales are not. Our methods for robust simplifications of nonlinear systems functions also when scales are partially separated: in this case we gather terms of the same order in the quasi-stationarity and averaged steady state equations.

The models that we study here are deterministic. Reduction methods for stochastic multiscale biochemical kinetics can be found in [32, 33].

The structure of this paper is the following. In the first section we present how to compute dominant subsystems for totally separated linear networks of (pseudo)monomolecular reactions. These appear as subsystems in analysis of multiscale networks of nonlinear biochemical reactions. This method uses the theory of limitation developed in [1]. In the second section, we show how to obtain dominant subsystems of non-linear systems. The technique is based on a method for identification of quasi-stationary and non-oscillating species and on dominant approximations of the quasi-stationarity and averaged steady-state equations for these species. In the third section, we introduce and analyze a new high dimensional model for the NF-*κ* B signalling.

## Methods

### Reduction of linear hierarchical models

#### Introductory notes

In this section we present a general algorithm for finding dominant subsystems and critical parameters for linear systems with completely separated time scales. Linear systems represent a special situation when all the interactions in the reaction network are monomolecular, i.e., have the form *A* → *B*.

Although systems biology models are nonlinear and contain also multimolecular reactions, it is nevertheless useful to have an efficient algorithm for solving linear problems. First, as we shall see in the next section, non-linear systems can include linear subsystems, containing reactions that are pseudo(monomolecular) with respect to species internal to the subsystem (at most one internal species is reactant and at most one is product). Second, for reactions *A* + *B* → ..., if concentrations *c*_{
A
}and *c*_{
B
}are well separated, say *c*_{
A
}>> *c*_{
B
}, then we can consider this reaction as *B* → ... with rate constant proportional to *c*_{
A
}which is practically constant, or changes only slowly. We can assume that this condition is satisfied for all but a small fraction of genuinely non-linear reactions (the set of non-linear reactions changes in time but remains small). Thus, linear models can serve as very effective approximations of behavior of non-linear models in certain windows of time: in this way, non-linear behavior can be approximated as a sequence of linear dynamics, followed one each other in a sequence of "phase transitions". Third, linear networks represent the case when very large reaction networks models can be approached analytically, and some intuition and design principles can be learned and partially generalized to the non-linear case. As an example, see the robustness study made in [34]. The linear case offers nice simple illustrations of the concepts of dominant subsystem, critical monomials and critical parameters.

The algorithm presented here in its "recipe" form ready for computational implementations, is developed in detail elsewhere [34], with many examples and rigorous justifications.

The structure of linear (monomolecular) reaction networks can be completely defined by a simple digraph, in which vertices correspond to chemical species *A*_{
i
}, edges correspond to reactions *A*_{
i
}→ *A*_{
j
}with kinetic constants *k*_{
ji
}> 0. For each vertex, *A*_{
i
}, a positive real variable *c*_{
i
}(concentration) is defined.

"Pseudo-species" (labeled ∅) can be defined to collect all degraded products, and degradation reactions can be written as *A*_{
i
}→ ∅ with constants *k*_{0i}. Production reactions can be represented as ∅ → *A*_{
i
}with rates *k*_{i 0}.

or in vector form: *ċ* = *K*_{0} + *K* *c*.

The advantage of linear dynamics is that it is completely specified by the eigenvectors and the eigenvalues of the kinetic matrix *K*.

The system has an unique bounded steady state *c*^{
s
}= *K*^{-1} *K*_{0} if and only if the matrix *K* is non-singular.

where *λ*_{
k
}, *l*^{
k
}, *r*^{
k
}, *k* = 1,..., *n* are the eigenvalues, the left eigenvectors (vector-rows) and the right eigenvectors (vector-columns) of the matrix *K*, respectively, i.e.

*K r*^{
k
}= *λ*_{
k
}*r*^{
k
}, *l*^{
k
}*K* = *λ*_{
k
}*l*^{
k
}.

with the normalization (*l*^{
i
}, *r*^{
j
}) = *δ*_{
ij
}, where *δ*_{
ij
}is Kronecker's delta.

*K*

_{0}= 0 (no production reactions, although degradation is permitted). Close systems are conservative if the matrix

*K*is singular (a particular case is when there is no degradation at all). Then, the left kernel of

*K*provides a set of conservation laws (if

*l K*= 0, then quantities (

*l, c*) are conserved). Solution of the homogeneous linear equations are simply:

If all reaction constants *k*_{
ij
}would be known with precision then the eigenvalues and the eigenvectors of the kinetic matrix can be easily calculated by standard numerical techniques. Furthermore, singular value decomposition can be used for model reduction. But in systems biology models often one has only approximate or relative values of the constants (information on which constant is bigger or smaller than another one). In the further we will consider the simplest case: when all kinetic constants are very different (separated), i.e. for any two different pairs of indices *I* = (*i*, *j*), *J* = (*i'*, *j'*) we have either *k*_{
I
}>> *k*_{
J
}or *k*_{
J
}>> *k*_{
I
}. In this case we say that the system is hierarchical with timescales (inverses of constants *k*_{
ij
}, *j* ≠ 0) *totally separated*.

*K*(originated from a reaction graph) allows us to exploit the strong relation between the dynamics (1) and the topological properties of the digraph. Big advantage of the fully separated network is that the possible values of ${l}_{i}^{k}$ are 0, 1 and the possible values of ${r}_{i}^{k}$ are -1, 0, 1 with high precision [34]. Thus, if we can provide an algorithm for finding non-zero components of ${l}_{i}^{k}$, ${r}_{i}^{k}$, based on the network topology and the constants ordering, then this will give us a good approximation to the problem solution (2).

#### Some basic notions

Two vertices of a graph are called adjacent if they share a common edge. A path is a sequence of adjacent vertices. A graph is connected if any two of its vertices are linked by a path. A maximal connected subgraph of graph *G* is called a connected component of *G*. Every graph can be decomposed into connected components.

A directed path is a sequence of adjacent edges where each step goes in direction of an edge. A vertex *A* is *reachable* from a vertex *B*, if there exists an oriented path from *B* to *A*.

A nonempty set *V* of graph vertexes forms a *sink*, if there are no oriented edges from *A*_{
i
}∈ *V* to any *A*_{
j
}∉ *V*. For example, in the reaction graph *A*_{1} ← *A*_{2} → *A*_{3} the one-vertex sets {*A*_{1}} and {*A*_{3}} are sinks. A sink is minimal if it does not contain a strictly smaller sink. In the previous example, {*A*_{1}}, {*A*_{3}} are minimal sinks. Minimal sinks are also called ergodic components.

A digraph is strongly connected, if every vertex *A* is reachable from any other vertex *B*. Ergodic components are maximal strongly connected subgraphs of the graph, but the reverse is not true: there may exist maximal strongly connected subgraphs that have outgoing edges and, therefore, are not sinks. If the digraph has no branching (each vertex has only one successor), then we can define a deterministic flow (discrete dynamical system) on the set of its vertices. Every vertex is the origin of an unique directed path.

#### Basic procedure for approximating eigenvectors

*κ*

_{ i }=

*k*

_{ ij }. Also it is useful to introduce a map

*ϕ*(see Fig. 1):

#### Acyclic non-branching network

In this case, for any vertex *A*_{
i
}there exists an eigenvector. If *A*_{
i
}is a sink vertex (i.e. *ϕ*(*i*) = *i*) then this eigenvalue is zero. If *A*_{
i
}is not a sink (i.e. *ϕ*(*i*) ≠ *i* and reaction *A*_{
i
}→ *A*_{ϕ(i)}has nonzero rate constant *κ*_{
i
}) then this eigenvector corresponds to eigenvalue -*κ*_{
i
}. For left and right eigenvectors of *K* that correspond to *A*_{
i
}we use notations *l*^{
i
}(vector-row) and *r*^{
i
}(vector-column), correspondingly.

*A*

_{ f }is a sink vertex of the network. Its associated right and left eigenvectors corresponding to the zero eigenvalue are given by:

Generally, right eigenvectors can be constructed by recurrence starting from the vertex *A*_{
i
}and moving in the direction of the flow. The construction is in opposite direction for left eigenvectors.

*r*

^{ i }only coordinates ${r}_{{\varphi}^{k}(i)}^{i}$ (

*k*= 0, 1, ..) can have nonzero values, and

*l*

^{ i }coordinate ${l}_{j}^{i}$ can have nonzero value only if there exists such

*q*≥ 0 that

*ϕ*

^{ q }(

*j*) =

*i*(this

*q*is unique because the system of reactions has no cycles), and

*κ*

_{ i }. Thus, for the left eigenvectors ${l}_{i}^{i}$ = 1 and, for

*i*≠

*j*,

*κ*

_{ f }= 0 for a sink vertex

*A*

_{ f }. Then ${r}_{i}^{i}$ = 1 and

*r*

^{ i }has at most two non-zero coordinates. The formula (10) means that to find the -1 component in

*r*

^{ i }, one should find the first vertex

*j*downstream of

*i*with

*κ*

_{ j }<

*κ*

_{ i }("bottleneck" vertex): there ${r}_{j}^{i}$ = -1. Following (10,9) we find that for the example at Fig. 1a

#### Non-branching network with a unique simple cyclic sink

In this case we have a reaction network with components *A*_{1}, ... *A*_{
n
}and last *τ* vertices (after some change of enumeration) form a reaction cycle: *A*_{n-τ+1}→ *A*_{n-τ+2}→ ... *A*_{
n
}→ *A*_{n-τ+1}. We assume that the limiting step in this cycle (reaction with minimal constant) is *A*_{
n
}→ *A*_{n-τ+1}.

*b*= ∑

_{ i }

*c*

_{ i }is the total (conserved) mass, then the steady state is:

for *j* ∈ [*n* - *τ* + 1, *n*] and zero elsewhere.

*κ*

_{ n }≪

*κ*

_{ i },

*i*≠

*n*) then this expression in the first order is simplified to

which means that most of substance is concentrated just before the "bottleneck" *A*_{
n
}→ *A*_{n-τ+1}(*c*_{
n
}≫ *c*_{
i
}, *i* ≠ *n*).

To approximate the dynamics of the reaction network for *κ*_{
n
}≪ *κ*_{
i
}, *i* ≠ *n*, it is sufficient to remove the slowest step of the cycle *A*_{
n
}→ *A*_{n-τ+1}. After removing, we will have acyclic non-branching system of reactions with eigenvalues and eigenvectors that can be computed from the formulas in the previous section. These formulas give *n* - 1 eigenvector sets corresponding to *n* - 1 non-zero eigenvalues *λ*_{
i
}= -*κ*_{
i
}, *i* = 1..*n* - 1. For example, removing *A*_{8} → *A*_{6} step at Fig. 1b converts the reaction network to the Fig. 1 a whose dynamics approximates the dynamics of the simple cyclic network.

#### Auxiliary reaction network and auxiliary dynamical system

Now let us consider an arbitrary linear reaction network with well-separated constants. For each *A*_{
i
}, let us define *κ*_{
i
}as the maximal kinetic constant for reactions *A*_{
i
}→ *A*_{
j
}: *κ*_{
i
}= max_{
j
}{*k*_{
ji
}}. For correspondent *j* we use notation *ϕ*(*i*): *ϕ*(*i*) = arg max_{
j
}{*k*_{
ji
}}. The function *ϕ*(*i*) is defined under condition that for *A*_{
i
}outgoing reactions *A*_{
i
}→ *A*_{
j
}exist. If there exist no such outgoing reactions then let us define *ϕ*(*i*) = *i*.

*A*

_{ i }→

*A*

_{ϕ(i)}with kinetic constants

*κ*

_{ i }. The correspondent kinetic equation is

The auxiliary network also defines a auxiliary discrete dynamical system *i* → Φ (*i*) that is used to compute the eigenvectors of the kinetic matrix. The auxiliary network can have several connected components. In each connected component the minimal sink is an attractor of the auxiliary dynamical system, hence it is either a node, or a cycle.

### General algorithm for calculating the dominant behavior of the linear dynamics

#### Preprocessing reaction network

- 1)
Let us consider a reaction network $\mathcal{W}$ with a given structure and fixed ordering of constants that are well separated. Using this ordering let us construct the auxiliary reaction network $\mathcal{V}={\mathcal{V}}_{\mathcal{W}}$.

- 2)
If the auxiliary network does not contain cycles then the auxiliary network kinetics (16) approximates relaxation of the initial network $\mathcal{W}$. To obtain the solution, we use directly formulas (7,8) to calculate the eigenvectors (if all

*κ*_{ i }are known) or (9,10) to obtain the 0–1 asymptotics (if only the ordering of*κ*_{ i }is known). - 3)
In general case, let the system $\mathcal{V}$ have several cycles

*C*_{1},*C*_{2}, ... with periods*τ*_{1},*τ*_{2}, ... > 1.

By "gluing" cycles into points, we transform the reaction network $\mathcal{W}$ into ${\mathcal{W}}^{1}$ as follows. For each cycle *C*_{
i
}, we introduce a new vertex *A*^{
i
}. The new set of vertices is ${\mathcal{A}}^{1}=\mathcal{A}\cup \{{A}^{1},{A}^{2},\mathrm{...}\}\backslash ({\cup}_{i}{C}_{i})$ (we delete cycles *C*_{
i
}and add vertices *A*^{
i
}).

*A*→

*B*(

*A*,

*B*∈ $\mathcal{A}$). They can be separated into 5 groups:

- 1.
both

*A*,*B*∉ ∪_{ i }*C*_{ i }; - 2.
*A*∉ ∪_{ i }*C*_{ i }, but*B*∈*C*_{ i }; - 3.
*A*∈*C*_{ i }, but*B*∉ ∪_{ i }*C*_{ i }; - 4.
*A*∈*C*_{ i },*B*∈*C*_{ j },*i*≠*j*; - 5.
*A*,*B*∈*C*_{ i }. - 1.
Reactions from the first group ("transitive" reactions) do not change.

- 2.
Reactions from the second group ("entering to cycles") transform into

*A*→*A*^{ i }(to the whole glued cycle) with the same constant. - 3.
Reactions of the third type ("exiting from cycles") change into

*A*^{ i }→*B*with the rate constant renormalization: let the cycle*C*^{ i }be the following sequence of reactions*A*_{1}→*A*_{2}→ ... ${A}_{{\tau}_{i}}$ →*A*_{1}, and the reaction rate constant for*A*_{ i }→*A*_{i+1}is*k*_{ i }(${k}_{{\tau}_{i}}$ for ${A}_{{\tau}_{i}}$ →*A*_{1}). For the limiting reaction of the cycle*C*_{ i }we use notation*k*_{lim i}. If*A*=*A*_{ j }and*k*is the rate reaction for*A*→*B*, then the new reaction*A*^{ i }→*B*has the rate constant*kk*_{lim i}/*k*_{ j }. This corresponds to a quasistationary distribution on the cycle (15). It is obvious that the new rate constant is smaller than the initial one:*kk*_{lim i}/*k*_{ j }<*k*, because*k*_{lim i}<*k*_{ j }due to definition of limiting constant. If after gluing, several reactions*A*^{ i }→*B*appear, then only the one with the maximal constant should be kept. - 4.
The same constant renormalization is necessary for reactions of the fourth type ("between cycles"). These reactions transform into

*A*^{ i }→*A*^{ j }. - 5.
Reactions of the fifth type ("inside cycles") are discarded.

- 4)
After the new network ${\mathcal{W}}^{1}$ is constructed, we assign $\mathcal{W}:={\mathcal{W}}^{1}$, $\mathcal{A}:={\mathcal{A}}^{1}$ and iterate the algorithm from the step 1) until we obtain an acyclic network and exit at step 2).

The algorithm produces an hierarchy of cycles. Notice that the algorithm is based on an asymmetry between entering reactions and outgoing reactions from cycles in the hierarchy. Indeed, some fluxes of $\mathcal{W}$ entering cycles *C*_{
i
}can be neglected when they are dominated by a stronger flux of $\mathcal{V}$ bifurcating from the same node (this occurs at the first step of the algorithm when constructing $\mathcal{V}$). The cycles *C*_{
i
}are minimal sinks in $\mathcal{V}$ (they are attractors of the auxiliary dynamical system). There are no reactions *A* → *B* in $\mathcal{V}$ such that *A* ∈ *C*_{
i
}, *B* ∉ *C*_{
i
}. Nevertheless, there may be such reactions in the initial network $\mathcal{W}$. These fluxes can not be neglected because there are no exiting fluxes of $\mathcal{V}$ to dominate them. The rule of thumb is: neglect any dominated flux except for the fluxes exiting some cycle in the hierarchy. This explains our algorithm and was rigorously justified in [34].

#### Constructing the dominant kinetic system

Now we show how to find an approximation of the dynamics of the reaction network $\mathcal{W}$. To construct this approximation, we produce a new acyclic reaction network with the initial set of vertices *A*_{
i
}∈ $\mathcal{W}$, *i* = 1..*n* which is called *dominant kinetic system*. Dynamics of this acyclic system can be computed from (7,8,9,10). To construct the dominant kinetic system, the following algorithm is applied:

- 1.
For ${\mathcal{V}}^{m}$ let us select the vertices ${A}_{1}^{m}$, ${A}_{2}^{m}$, ... that are glued cycles from ${\mathcal{V}}^{m-1}$.

- 2.
For each glued cycle node ${A}_{i}^{m}$:

- a)
Recall its nodes ${A}_{i1}^{m-1}\to {A}_{i2}^{m-1}\to \dots {A}_{i{\tau}_{i}}^{m-1}\to {A}_{i1}^{m-1}$; they form a cycle of length

*τ*_{ i }, - b)
Let us assume that the limiting step in ${A}_{i}^{m}$ is ${A}_{i{\tau}_{i}}^{m-1}\to {A}_{i1}^{m-1}$,

- c)
Remove ${A}_{i}^{m}$ from ${\mathcal{V}}^{m}$,

- d)
Add

*τ*_{ i }vertices ${A}_{i1}^{m-1},{A}_{i2}^{m-1},\dots {A}_{i{\tau}_{i}}^{m-1}$ to ${\mathcal{V}}^{m}$, - e)
Add to ${\mathcal{V}}^{m}$ reactions ${A}_{i1}^{m-1}\to {A}_{i2}^{m-1}\to \dots {A}_{i{\tau}_{i}}^{m-1}$ (that are the cycle reactions without the limiting step) with correspondent constants from ${\mathcal{V}}^{m-1}$,

- f)
If there exists an outgoing reaction ${A}_{i}^{m}$ →

*B*in ${\mathcal{V}}^{m}$ then we substitute it by the reaction ${A}_{i{\tau}_{i}}^{m-1}$ →*B*with the same constant, i.e. outgoing reactions ${A}_{i}^{m}$ → ... are reattached to the heads of the limiting steps, - g)
If there exists an incoming reaction in the form

*B*→ ${A}_{i}^{m}$, find its prototype in ${\mathcal{V}}^{m-1}$ restore it in ${\mathcal{V}}^{m}$. - 3.
If in the initial ${\mathcal{V}}^{m}$ there existed a "between-cycles" reaction ${A}_{i}^{m}\to {A}_{j}^{m}$ then we find the prototype in ${\mathcal{V}}^{m-1}$,

*A*→*B*, and substitute the reaction by ${A}_{i{\tau}_{i}}^{m-1}$ →*B*with the same constant, as for ${A}_{i}^{m}\to {A}_{j}^{m}$ (again, the beginning of the arrow is reattached to the head of the limiting step in ${A}_{i}^{m}$). - 4.
Let

*m*←*m*- 1, and repeat steps 1–4 until no glued cycles left.

One has to notice that in the process of network preprocessing some reaction rates are substituted by monomials of the initial reaction constants, i.e. expressions in the form ${k}_{new}=\frac{{k}_{{i}_{1},{i}_{2}}\mathrm{...}{k}_{{j}_{1},{j}_{2}}}{{k}_{{l}_{1},{l}_{2}}\mathrm{...}{k}_{{m}_{1}}{,}_{{m}_{2}}}$. In totally separated case the values of these monomials are also well separated from the other constants with probability close to 1, however, the initial order of constants does not prescribe position of these monomials in the rates order. In this case the algorithm produces several dominant systems defined for all possible position of new monomial rate constants in the order. An example of this will be given later in this section. Such a situation can happen during the network preprocessing when maximum reaction constant should be chosen, or in the process of dominant system creation, when determining the limiting step in a cycle.

#### Finding stationary distributions

The dominant kinetic system fully describes the relaxation modes of the network. The construction of this system depends only on the matrix *K* and does not depend on the production reactions *K*_{0}. In particular, relaxation times do not depend on the system being closed or not. However, the stationary distribution *c*^{
s
}and the sequence of relaxation events depends on production reactions (see Eq.(2)).

For closed systems, steady states are solutions of the linear homogeneous equation *K**c* = 0, therefore they are determined up to multiplication by positive constants. They form a *k*-dimensional cone where *k* is the multiplicity of the zero eigenvalue of the matrix *K*, also the number of minimal sinks of the network.

*k*sink vertices of the auxiliary network ${\mathcal{V}}^{m}$. Let

*A*

_{ i },

*i*= 1..

*n*be vertices in the initial network $\mathcal{W}$. Below we describe a procedure for finding the basis of all stationary distributions of a closed network:

- 1.
Let us take the

*i*th sink vertex ${A}_{fi}^{m}$. - 2.
Define

*x*= ${A}_{fi}^{m}$,*b*= 1, and a null vector*b*^{ i }∈*R*^{ n }. - 3.
If

*x*is not a glued cycle then it corresponds to a vertex*A*_{ j }∈ $\mathcal{W}$, and the basis vector*b*^{ i }has components ${b}_{j}^{i}$ =*δ*_{ ij }; stop. - 4.
If

*x*is a glued cycle then recall all its vertices*x*^{1}, ...,*x*^{ τ }. - 5.
Determine the limiting (minimal) rate constant

*κ*_{ lim }= min_{s=1..τ}{${\kappa}_{{x}^{s}}$} and*s*_{ min }= arg min_{s=1..τ}{${\kappa}_{{x}^{s}}$}. - 6.
For each vertex

*x*^{ j }of the cycle repeat: - a)
Let ${c}_{j}=b\frac{{\kappa}_{lim}}{{\kappa}_{{x}^{j}}}$ if

*j*≠*s*_{ min }and ${c}_{j}=b\{1-{\displaystyle {\sum}_{s\ne {s}_{min}}\frac{{\kappa}_{lim}}{{\kappa}_{{x}^{s}}}}\}$ otherwise, - b)
if

*x*^{ j }corresponds to a simple vertex*A*_{ k }then ${b}_{k}^{i}$ =*c*_{ j }, - c)
if

*x*^{ j }corresponds to a glued cycle then do recursively steps 4–6 with*x*=*x*^{ j }and*b*=*c*_{ j }.

Any possible stationary distribution has form ${c}^{s}={\displaystyle {\sum}_{i=\mathrm{1..}k}{c}_{fi}^{m}{b}^{i}},{c}_{fi}^{m}\ge 0$, ${c}_{fi}^{m}$. The coefficients ${\tilde{c}}^{s}={\displaystyle {\sum}_{i=\mathrm{1..}k}{c}_{fi}^{m}{\tilde{b}}^{i}}$ are computed from initial data: they are equal to the total initial mass carried by vertices of ${\mathcal{V}}^{m}$ (when these are glued cycles we consider the total initial mass of the cycle) that are attracted by ${A}_{fi}^{m}$.

In brief, the distribution of the concentrations on any cycle is approximated by the first order expression (15), and this procedure is applied recursively for the vertices that represent glued cycles. The state thus obtained is equally a good approximation of the steady state of the dominant kinetic system.

Open systems can be reduced to closed ones by considering that all production reactions originate from the node ∅ that has concentration *c*_{0} = 1. The corresponding reaction are ∅ → *A*_{
i
}and the constants are the production rates *k*_{i 0}. The normalization *c*_{0} = 1 is possible for all bounded steady states, because these are determined up to multiplication by a constant. Furthermore, all steady states are bounded, provided that the following topological condition holds: if there exists an oriented path from ∅ to *A*_{
i
}, then there exists an oriented path from *A*_{
i
}to ∅. We suppose this condition to be always fulfilled. Applying the algorithm to the closed system we obtain ${\tilde{c}}^{s}={\tilde{c}}^{s}/({\displaystyle {\sum}_{i=\mathrm{1..}k}{c}_{fi}^{m}{\tilde{b}}_{0}^{i}})$ that is normalized to ${\tilde{c}}^{s}={\tilde{c}}^{s}/({\displaystyle {\sum}_{i=\mathrm{1..}k}{c}_{fi}^{m}{\tilde{b}}_{0}^{i}})$ in order to have ${c}_{0}^{s}$ = 1.

**Example**. Let us consider an example of the network $\mathcal{W}$ shown on Fig. 2 (1). Below we briefly detail every step of the algorithm.

- 2)
An auxiliary reaction network ${\mathcal{V}}_{\mathcal{W}}$ is constructed (this gives a non-branching network);

- 3)
The cycle

*A*_{3}→*A*_{4}→*A*_{5}→*A*_{3}in the auxiliary reaction network is glued in one vertex ${A}_{3}^{1}$ (shown by the circle node); In the initial network we find an "exiting from cycle" reaction*A*_{4}→*A*_{2}, renormalize its rate to ${k}_{32}^{1}=\frac{{k}_{24}{k}_{35}}{{k}_{54}}$ and insert in the new network ${\mathcal{W}}^{1}$; - 4)
The cycle

*A*_{3}→*A*_{2}→*A*_{3}in the auxiliary reaction network ${V}_{{\mathcal{W}}^{1}}$ (which is now coincides with ${\mathcal{W}}^{1}$ itself) is glued in one vertex ${A}_{2}^{2}$; now the network ${\mathcal{W}}^{2}$ is acyclic and we stop the network preprocessing.

Now if we restore the cycle *A*_{3} → *A*_{2} → *A*_{3} and try to determine the limiting step in it, we have two possibilies: $\frac{{k}_{24}{k}_{35}}{{k}_{54}}$ <*k*_{32} and $\frac{{k}_{24}{k}_{35}}{{k}_{54}}$ > *k*_{32}. Let us consider them separately:

*Case* $\frac{{k}_{24}{k}_{35}}{{k}_{54}}$ <*k*_{32}

3.1.1) Since ${A}_{3}^{1}$ <*k*_{32}, then we remove the limiting step ${A}_{3}^{1}$ → *A*_{2} and obtain 3.1.2).

3.1.3) We restore the glued cycle corresponding to ${A}_{3}^{1}$ and we recall that the reaction *A*_{2} → ${A}_{3}^{1}$ in ${\mathcal{W}}^{1}$ corresponds to *A*_{2} → *A*_{3} in $\mathcal{W}$.

3.1.4) We remove the limiting step reaction in the cycle *A*_{3} → *A*_{4} → *A*_{5} → *A*_{3} (it is *A*_{5} → *A*_{3}) and as a result we obtain acyclic dominant kinetic system shown at Fig. 2 (3.1.4).

*Case* $\frac{{k}_{24}{k}_{35}}{{k}_{54}}$ > *k*_{32}

3.2.1) Since $\frac{{k}_{24}{k}_{35}}{{k}_{54}}$ > *k*_{32}, then we remove the limiting step *A*_{2} → ${A}_{3}^{1}$ and obtain 3.1.2).

3.2.3) We restore the glued cycle corresponding to ${A}_{3}^{1}$ this time we should re-attach the reaction ${A}_{3}^{1}$ → *A*_{2} to the head of the limiting step in the cycle (it is *A*_{5} vertex); the rate of *A*_{5} → *A*_{2} is $\frac{{k}_{24}{k}_{35}}{{k}_{54}}$.

3.2.4) We remove the limiting step reaction in the cycle *A*_{3} → *A*_{4} → *A*_{5} → *A*_{3} (it is *A*_{5} → *A*_{3}) and as a result we obtain acyclic dominant kinetic system shown at Fig. 2 (3.2.4).

### Discussion and perspectives

Dominant approximations of hierarchical linear reaction network allow us to introduce some new concepts important for the dynamics of multiscale systems.

#### Hybrid and qualitative dynamics

Piecewise affine dynamics has been widely used to approximate dynamics of gene regulatory networks [35–37] as a sequence of discrete transitions between attractors of affine systems. This picture is based on threshold response of genes in models with steep regulation functions (Hill functions and other representations of sigmoidal response) and is not directly related to time scales. Here, we emphasize another possible way to obtain hybrid, or qualitative representations of dynamics, based on time separation.

Indeed, zero-one approximation of eigenvectors in hierarchical linear systems justifies a discrete coding of dynamics. Suppose that initial state is concentrated in *j*_{0}, *c*_{
i
}(0) ~ ${\delta}_{i,{j}_{0}}$. At times just larger than 1/*λ*_{
k
}an exponential vanishes in Eq. (2) and the state has a "jump" -*r*^{
k
}(*l*^{
k
}, *c*(0) -*c*_{
s
}). Let us consider that eigenvalues are ordered *λ*_{1} >> *λ*_{2} >> ... >> *λ*_{n-1}. Then, the sequence of right eigenvectors *r*^{
k
}such that (*l*^{
k
}, *c*(0) -*c*_{
s
}) ≠ 0 codes the dynamics starting in *c*(0). In other words, there is a sequence of well separated times *τ*_{1} = 1/*λ*_{1} <<*λ*_{2} = 1/*λ*_{2} << ... <<*τ*_{n-1}= 1/*λ*_{n-1}such that something happens (a state transition) between each one of these times. Left eigenvectors provides the lumping (several species cumulated to form pseudospecies) and right eigenvectors provide the sequence of state transitions. On timescales *τ*_{
k
}<*t* <*τ*_{k+1}one can observe a jump -*r*^{
k
}in state space provided that (*l*^{
k
}, *c*(0) - *c*^{
s
}) ≠ 0. On this timescale the dynamics is equivalent to the degradation of pseudospecies (*l*^{
k
}, *c*), $\frac{d({l}^{k},c)}{dt}$ = -*λ*_{
k
}(*l*^{
k
}, *c*).

#### Critical parameters and design principles

Our approach to dominant subsystems emphasizes some simple but important principles. First of all, dynamics of a hierarchical linear network can be specified if a) the topology of the network is given and if b) to each reaction we associate a positive integer representing order (1 for the most rapid reaction, 2 for the second most rapid reaction, and so on ...); c) for cyclic topologies, some monomials grouping constants of several reactions have also to be ordered in the same manner (which reactions depend both on topology and on initial ordering).

In the process of simplification some reaction pathways are dominated and do not appear in the dominant subsystem. Therefore, the corresponding constants are not critical for the system: although their ordering matters for establishing the simplification, their precise value have little importance. Because parameters of the dominant subsystem are generally monomials of parameters of the whole system, critical parameters are those parameters that occur in critical monomials. Our findings show rather counter-intuitive properties of critical and non-critical parameters, that can be useful as design principles. Thus, in cycles, the limiting step (slowest reaction) has little influence on dynamics (though is important for the steady state). Dynamically, a cycle with separated constants behaves like the chain obtained from the cycle by eliminating the limiting step. In particular, the slowest relaxation time of a cycle is the inverse of its *second* slowest constant [1, 34].

- 1)
*When some submechanisms of a complex and non-linear network are linear*, given fixed (or slowly changing) values of external inputs (boundaries); - 2)
*For approximating non-linear dynamics*. For multiscale nonlinear reaction networks the expected dynamical behaviour is to be approximated by the system of dominant networks. These networks may change in time but remain small enough. To give an example, we provided the Fig. 3S1–3S2 in Additional File 1 demonstrating that in a model of complex reaction network of NF-*κ*B pathway, containing 17 multimolecular reactions, only two reactions show genuinely non-linear behavior in some windows of time, with two more showing border-line behavior, and all others have well-separated reactant concentrations in any moment of time. The rigorous justification of these hybrid approximations for mass action reaction networks will be discussed elsewhere.

### Reduction of non-linear multiscale systems

Complex formation is a source of nonlinearity in biochemical networks. For instance, in signalling, ligand molecules form complexes with receptors. Transcription factors are often dimers or multimers or are sequestered by forming complexes with their inhibitors. In these examples, the reaction rates are non-linear functions of the concentrations of two or more molecules.

*A*

_{1}, ...,

*A*

_{ n }} and the list of reactions (the reaction mechanism):

where *j* ∈ [1, *r*] is the reaction number. Unless reactants and products belong to compartments of different volumes *α*_{
ji
}, *β*_{
jk
}are nonnegative integers (stoichiometric coefficients). Reactions involving components from different compartments have non-integer stoichiometry. For instance, a reaction translocating a molecule from nucleus to the cytosol has stoichiometry (..., 1, *k*_{
v
}, ...) where *k*_{
v
}is the volume ratio of cytosol to nucleus.

*ν*_{
j
}= *β*_{
j
}- *α*_{
j
}is the global stoichiometric vector. *S* is the stoichiometric matrix whose columns are the vectors *ν*_{
j
}. The reaction rates ${R}_{j}^{+/-}$ (*c*) are non-linear functions of the concentrations. For instance the mass action law reads ${R}_{j}^{+}(c)={k}_{j}^{+}{\displaystyle {\prod}_{i}{c}_{i}^{{\alpha}_{ji}}},{R}_{j}^{-}(c)={k}_{j}^{-}{\displaystyle {\prod}_{i}{c}_{i}^{{\beta}_{ji}}}$.

There are no simple rules to relate timescales to reaction constants of nonlinear models. The units of the inverse constants of bimolecular reactions are concentration multiplied by time and one needs at least one concentration value in order to construct a timescale. Generally, timescales are functions of many reaction constants and concentration variables. These functions are not necessarily smooth. Near bifurcations (for instance, near Hopf or saddle-node bifurcations), at least one timescale of the system diverges for finite changes of the reaction constants. However, nonlinear biochemical networks have wide distributions of time-scales, as can be shown by simple (Jacobian based) analysis of models.

Various reduction methods of nonlinear models are based on projection of the dynamics on a lower dimensional invariant manifold [4–8]. The reduced models are systems of differential equations, but no longer networks of chemical reactions. Quasi-equilibrium and quasi-stationarity methods keep the network structure of the model and propose lumped reaction mechanisms as dominant subsystems. This approach has some advantages. Indeed, it leads to more transparent analysis of the results and of design principles, produces hierarchies of models and facilitates model comparison. Graphical reduction methods using elementary modes, were proposed by Clarke [26] for chemical systems and more recently in systems biology by Klamt [38]. Similar methods can be found in [39], from which we have borrowed the terminology. The choice of the species to be eliminated and of the reactions to be aggregated, as well as the calculation of rates of elementary modes have no theoretical justification in these methods and their inappropriate use can alter dynamics (for instance, as Clarke noticed, the stability of limit cycles is not guaranteed). Thus, in order to have a complete practical recipe that applies to multiscale biochemical networks we need to solve three more problems: detection of rapid species, resolution of quasi-stationarity equations and calculation of reaction rates of the dominant mechanisms.

A major improvement in calculating dominant subsystems can be obtained by combining quasi-stationarity and averaging. Averaging techniques are widely used in physics and chemistry to simplify models by eliminating fast, oscillating (microscopic) variables [22–24]. Our use of averaging is different, because we employ it to obtain averaged stationarity equations for slow, non-oscillating variables and to eliminate these species. After choosing a "middle" time scale (corresponding to the time resolution of the experiment), we want to reduce all scales that are faster but also all scales that are slower than this middle scale. In order to do that we provide a unified framework for species elimination and reaction aggregation, either by quasi-stationarity (fast species) or by averaged stationarity (slow species).

Let *I* be the set of indices of intermediate components, that will be eliminated. ${\mathcal{R}}^{(I)}$ is the set of reactions that either produce or consume species from *I*. Rates of ${\mathcal{R}}^{(I)}$ depend on the concentrations of intermediate species and also on the concentrations of other species, which in the terminology of Temkin [39] are called terminal. Let *T* be the set of indices of terminal species. Terminal species represent the *frontier* between the rest of the system and the subsystems made of intermediate species. Although instead of terminal the name boundary species could be more appropriate, the latter term has already been employed in systems biology with a different meaning, which is species whose concentrations are fixed in a simulation.

Extracting from the matrix *S* the columns corresponding to the reactions ${\mathcal{R}}_{I}$ and the lines corresponding to the species $\mathcal{I}$ and $\mathcal{T}$ we obtain the intermediate stoichiometric matrix *S*_{
I
}and the terminal stoichiometric matrix *S*_{
T
}, respectively.

#### Eliminating fast species: quasi-stationarity

In multiscale biochemical systems, some components react much more rapidly to changes in the environment than others. The reasons for the existence of such fast species can be multiple. Thus, rapidly transformed or rapidly consumed molecules (for instance those taking part in metabolic chains or rapid chemical transformations such as phosphorylation), or promoter sites submitted to rapid binding/unbinding processes are examples of fast species. Fast species are good candidates for intermediate species. Indeed, it is easy to prove that they can be eliminated by quasistationarity. When production rates are not weak, fast species are those whose concentrations are small and well separated from the concentrations of other species. Though straightforward, the precise condition connecting quasi-stationarity and smallness of concentrations can not be easily found in literature, hence we briefly discuss it below.

Let *ϵ* be a small parameter, representing concentrations. Suppose for simplicity that reactions ${\mathcal{R}}^{(I)}$ are pseudo-monomolecular. This means that *S*_{
I
}*R*_{
I
}(*c*_{
I
}, *c*_{
T
}) = *K*_{
I
}(*c*_{
T
}) *c*_{
I
}+ ${K}_{I}^{0}$(*c*_{
T
}), where *R*_{
I
}is the restriction of the vector *R* to the intermediate species. An important assumption is ${K}_{I}^{0}$(*c*_{
T
}) = $\mathcal{O}$(1) meaning that the production of intermediate species is not weak.

Suppose that among the reactions ${\mathcal{R}}^{(I)}$ consuming intermediates, at least some have rates of order $\mathcal{O}$(1). This is current, because these reactions produce terminal species which have larger concentrations.

*c*

_{ I }= $\mathcal{O}$(

*ϵ*), it means that

*K*

_{ I }(

*c*

_{ T }) = $\mathcal{O}$(1/

*ϵ*). This leads to the following asymptotic:

*c*

_{ I }/

*ϵ*, ${\tilde{K}}_{I}$ =

*ϵ*

*K*

_{ I }= $\mathcal{O}$(1),

*I*

^{ c }is the complement of

*I*designating species other that

*I*. Intermediate species are fast and the system (19) can be reduced using Tikhonov's [40, 41] and Fenichel's [42] results. According to these results, after a short laps of time, the system evolves on an

*invariant manifold*(an invariant manifold is defined by the property that any trajectory starting in the manifold stays inside the manifold) which is at distance $\mathcal{O}$(

*ϵ*) from the quasi-steady state (QSS) manifold defined by the quasi-stationarity equations:

Quasi-stationarity equations can be used to express concentrations of the intermediate species as functions of the concentrations of terminal species. If matrix *K*_{
I
}has not full rank, conservation laws should be added to the quasi-stationarity equations in order to obtain a full rank system. Let *μ*_{1}, ..., *μ*_{
k
}be a basis of the left kernel of *S*_{
I
}(a complete set of conservation laws). We say that species of indices *I* are quasi-stationary if they approximately fulfill the equations:

*S*_{
I
}*R*_{
I
}(*c*_{
I
}, *c*_{
T
}) = 0

and exactly fulfill the conservation laws:

*μ*_{
i
}*i*_{
I
}= *C*_{
i
}, *i* = 1,..., *k*,

where *C*_{
i
}are real constants.

*c*

_{ T }) = $\mathcal{O}$(1), although informative for understanding of the dynamics, can not be used in practice. Furthermore, small concentration is not a necessary condition for quasi-stationarity. Therefore, our practical method for detection of fast, quasi-stationary species is based on the direct checking of Eqs.(22), (23) (see Fig 3a and the Results section for an example).

- 1.
Eliminate the intermediate concentrations by solving the equations (22), (23). Express

*c*_{ I }as function of*c*_{ T }. - 2.
Replace the mechanism ${\mathcal{R}}^{(I)}$ by "simple sub-mechanisms".

- 3.
Compute the rates of the simple sub-mechanisms as functions of

*c*_{ T }.

The simplicity criterion employed by Clarke does not follow from a physical principle. Nevertheless, in systems biology, biochemical reactions are simplified representations of complex physico-chemical processes. In the absence of detailed information, simplicity arguments are often employed. Elementary modes analysis widely used in metabolic control and gene network analysis [43–45] is based on exactly the same argument.

The same recipe applies also to model comparison, when we want to compare two models which differ in complexity (some species in one model are not present in the second). In this situation we declare the extra species intermediate and apply the three steps of the algorithm.

#### Simple sub-mechanisms and rates

Let us introduce some more definitions. A reaction route is a combination of reactions in ${\mathcal{R}}_{I}$ transforming terminal species into other terminal species and conserving the intermediate species. It is defined by a integer coefficient vector *γ* ∈ ℤ^{
s
} (the dimension *s* is the number of reactions in ${\mathcal{R}}_{I}$) satisfying the following three conditions:

*S*_{
I
}*γ* = 0

*γ*_{
i
}≥ 0, if the reaction *i* is irreversible

||*S*_{
T
}*γ*|| > 0

Reaction routes are usually defined [39] without the condition (26). By imposing this condition, we exclude internal cycles with zero terminal stoichiometry.

A sub-mechanism *M*(*γ*) is the set of all the reactions in the reaction route *γ*, *M*(*γ*) = {*i*|*γ*_{
i
}≠ 0}. A sub-mechanism is simple if it is minimal with respect to inclusion, i.e. if *M* (*γ'*) ⊂ *M* (*γ*) ⇒ *γ* = *γ'*. Simple sub-mechanisms are pathways with a minimal number of reactions, connecting terminal species without producing accumulation or depletion of the intermediate species. Thus, they are candidates for reduced reaction mechanisms. Simple sub-mechanisms are minimal dependent sets in oriented matroids [46], similar to elementary modes in flux balance analysis [43]. Algorithms for finding elementary modes can be applied for the search of simple sub-mechanisms [43–45].

In the reduced model, the reactions of the intermediate mechanism ${\mathcal{R}}_{I}$ are replaced by the sub-mechanisms *γ*_{1}, ..., *γ*_{
s
}.

*c*

_{ T }and

*c*

_{ I }satisfying (22),(23)

^{:}

where ${{R}^{\prime}}_{i}(c)$ are the rates of the simple sub-mechanisms.

*i*there is a terminal species

*j*such that

*S*

_{ T }

*γ*

_{ i }is the unique vector (among the

*s*different ones) having nonzero coordinate

*j*, (

*S*

_{ T }

*γ*

_{ i })

_{ j }≠ 0. Then, there is a straightforward solution for (27):

The above uniqueness condition is not fulfilled if there are two sub-mechanisms for which the terminal stoichiometries are proportional. This situation can be avoided by quotienting with respect to the following equivalence relation: *γ*_{
i
}and *γ*_{
j
}are equivalent iff *S*_{
T
}*γ*_{
i
}= *α S*_{
T
}*γ*_{
j
}, for some *α* = ≠ 0. After discarding some sub-mechanisms and keeping only one representative per class, we have a reduced set of simple sub-mechanisms for which rates can be calculated from (28).

#### Dominant solutions to the quasi-stationarity equations, multiscale ensembles

The most difficult part of the above algorithm is to solve the quasi-stationarity equations (22),(23). Even in the monomolecular case, symbolic solutions of the linear system (21) can involve long expressions. Furthermore, mass action law leads to polynomial equations in the binary or multi-molecular case. Symbolic methods for solutions of systems of polynomial equations are limited to a small number of variables.

In this subsection we show how the multi-scale nature of the system can be used to obtain approximate, dominant solutions of the quasi-stationarity equations.

*k*∈ [

*α*,

*β*], but most of the properties of this probability distribution will not be used here. The only property that we will use is the following: if

*k*

_{ i }>

*k*

_{ j }, then

*k*

_{ i }/

*k*

_{ j }≫ 1 (with probability close to one). It means that we can assume that

*k*

_{ i }/

*k*

_{ j }≫

*a*for any preassigned positive value of

*a*that does not depend on

*k*values. One can interpret this property as an asymptotic one for

*α*→ -∞,

*β*→ ∞. This property allows us to simplify algebraic formulas. For example,

*k*

_{ i }+

*k*

_{ j }can be substituted by max{

*k*

_{ i },

*k*

_{ j }} (with small relative error), or

for nonzero *a*, *b*, *c*, *d*.

Of course, some ambiguity can be introduced, for example, what is (*k*_{1} + *k*_{2}) - *k*_{1}, if *k*_{1} ≫ *k*_{2}? If we first simplify the expression in brackets, it is zero, but if we open brackets without simplification, it is *k*_{2}. This is a standard difficulty in use of relative errors for round-off. If we estimate the error in the final answer, and then simplify, we shall avoid this difficulty. Use of *o* and $\mathcal{O}$ symbols also helps to control the error qualitatively: if *k*_{1} ≫ *k*_{2}, then we can write (*k*_{1} + *k*_{2}) = *k*_{1}(1 + *o*(1)), and *k*_{1}(1 + *o*(1) - *k*_{1} = *k*_{1}*o*(1). The last expression is neither zero, nor absolutely small – it is just relatively small with respect to *k*_{1}.

It is slightly more difficult to solve equations. Some recipes were proposed such as Newton polyhedra for approximate solutions of polynomial systems of equations [47] but this type of methods suffers from combinatorial complexity. Here we use a simpler, but not so rigorous approach. In the case of pseudo-molecular subsystems, our algorithms for linear hierarchical systems are enough for this purpose. In general, we choose the dominant terms in the solutions as monomials of the parameters. This can be done either by educated guess, or by testing numerically the orders of various terms in the equations. The most frequent, truly non-linear simplification that occurs in biochemical models is the "min-funnel", which we present below.

Let us consider the production of a complex C from two proteins A and $\text{B}:\varnothing \underset{}{\overset{{k}_{A}}{\to}}\text{A}$ (production of A), $\varnothing \underset{}{\overset{{k}_{B}}{\to}}\text{B}$ (production of B), $\text{A}\underset{}{\overset{{k}_{deg,A}}{\to}}\varnothing $ (degradation of A), $\text{B}\underset{}{\overset{{k}_{deg,B}}{\to}}\varnothing $ (degradation of B), $\text{A}+\text{B}\underset{}{\overset{{k}_{c}}{\rightleftharpoons}}\text{C}$ (complex formation).

Supposing A, B quasi-stationary we have to find the positive solutions of the equations ${k}_{A}=\tilde{A}+{\tilde{k}}_{c}\tilde{A}\tilde{B}$, ${k}_{B}=\tilde{B}+{\tilde{k}}_{c}\tilde{A}\tilde{B}$, where $\tilde{A}={k}_{deg,A}A$, $\tilde{B}={k}_{deg,B}B$, ${\tilde{k}}_{c}={k}_{c}/({k}_{deg,A}{k}_{deg,B})$. We will consider two cases a) 1/${\tilde{k}}_{c}$ <<*k*_{
A
}<<*k*_{
B
}and b) 1/${\tilde{k}}_{c}$ <<*k*_{
B
}<<*k*_{
A
}. Both cases mean that degradation of *A*, *B* is weak and/or the propensity of complex formation is high. Case a) means also that *B* is in excess, the opposite being true in case b).

Let us consider the case a). We consider that the order of $\tilde{A}$ in the dominant solution is larger than the order of $\tilde{B}$, $\tilde{A}<<\tilde{B}$. From the linear equation *k*_{
A
}- *k*_{
B
}= $\tilde{A}-\tilde{B}$ we obtain $\tilde{B}$ = *k*_{
B
}and from the second nonlinear equation we obtain $\tilde{A}=\frac{{k}_{A}}{1+{\tilde{k}}_{c}{k}_{B}}\approx \frac{{k}_{A}}{{\tilde{k}}_{c}{k}_{B}}$. Finally, we have $\tilde{A}<<\tilde{B}$ consistently with the starting guess. The dominant solution in case b) is obtained by symmetry from the one in case a). The quantity of interest in this example, for which we want a reduced expression is the production rate of the complex *R*_{
c
}= *k*_{
c
}*AB*. Actually, the two solutions can be summarized by:

*R*_{
c
}= *min*(*k*_{
A
}, *k*_{
B
})

Using the exact solution of the system (after eliminating *A* from the linear equation we remain with a quadratic equation for *B*) we can show that the min-funnel approximation (29) is valid under less restrictive conditions. The only separation condition that we need is *min*(*k*_{
A
}, *k*_{
B
}) >> *k*_{
deg, A
}*k*_{deg,B}/*k*_{
c
}. We can easily identify the critical parameters *k*_{
A
}, *k*_{
B
}and the non-critical ones *k*_{deg,A}, *k*_{deg,B}, *k*_{
c
}. The validity of the expression (29) depends on order relations involving monomials of critical and non-critical parameters.

#### Eliminating slow species: averaging

Averaging is an useful model reduction technique for high-dimensional clocks or for other types of oscillating molecular systems (the activity of some transcription factors, among which NF-*κ* B, present oscillations under some conditions).

*x*, the non-oscillating species are

*y*, and ∈ is a small parameter, then we have:

*y*, the fast dynamics (30) has an attractive hyperbolic limit cycle

*x*=

*ψ*(

*τ*,

*y*), of period

*T*(

*y*):

*ψ*(

*τ*+

*T*(

*y*),

*y*) =

*ψ*(

*τ*,

*y*) (

*τ*=

*t*/

*ϵ*). Then, after a short transient, the slow variable satisfies the averaged equation:

The result can be extended to the case when *x* = *ψ*(*τ*, *y*) describes damped oscillations, with damping time much larger than *ϵ*, i.e. when the fast dynamics (30) has a stable focus and the eigenvalues of the Jacobian ∂*f*/∂*x* calculated at the focus are of the form -*λ* ± *iμ*, 0 <*λ* <<*μ* = $\mathcal{O}$(1/*ϵ*).

*y*:

If (32) has a stable steady state, we always reach this situation. In this case, the slow non-oscillating variables *y* are constant in time and can be considered to be conserved, which has two significant consequences.

First, Eq.(33) restores conservation. Slow variables are often the result of broken conservation laws. In fact, in biological open systems, nothing is conserved. Conservation laws result from balancing production and degradation either passively (slow processes) or actively (feed-back). Thus, we can ignore production and degradation of molecules whose level is rigorously controlled. Eq.(33) describes such a case.

Second, (33) are averaged steady state equations for the slow variables. If slow variables *y* reach stationarity, the only variables that change in time are the oscillating variables *x*. Eq.(30) describes the dynamics of *x*, considering that *y* satisfies (33). For oscillators, averaging provides a way to eliminate slow non-oscillating variables, which is formally equivalent to quasi-stationarity and represents a new case of applicability of Clarke's method. The difference between the two cases is that we eliminate fast variables by solving quasi-stationarity equations and we eliminate slow variables by solving averaged stationarity equations. Thus, intermediate non-oscillating variables are expressed in terms of only non-oscillating terminal variables. If there are no non-oscillating terminal variables, then non-oscillating intermediates become conserved quantities.

## Results and discussion

### Methodology

In this section, we demonstrate hierarchical model reduction, model comparison and critical parameter identification. Critical parameters are identified during the reduction procedure.

Model reduction starts with a complex model, from which we obtain a hierarchy of reduced models by eliminating various intermediate species. The intermediate species are either quasi-stationary species (in general), or non-oscillating species (for oscillators). The complexity of a model is quantified by three integers. A model with *n* species, *r* reactions and *p* parameters is designated by *M*(*n, r, p*). Our conception about systems biology models is summarized by the following idea. Instead of providing a single model, it is better to provide a hierarchy of models, and the relations between them. Depending on the application, we can choose the most appropriate model in the hierarchy or couple several simple models into a larger model.

The number of parameters in a model are obtained as follows. If the elementary reactions follow mass action law kinetics, there are *n*_{
k
}= 2*n*_{
r
}+ *n*_{
i
}kinetic constants, where *n*_{
r
}, *n*_{
i
}are the numbers of reversible and irreversible reactions. Reactions with kinetic constants zero are not considered in the counting. Each one of the *n*_{
c
}conservation laws adds an extra parameter, the value of the conserved quantity. These values follow from initial data and are important parameters for the dynamics. For multi-compartment models, the ratios of the compartments volumes (in the example below there is only one ratio *kv*, the cytoplasm to nucleus volume ratio) are extra parameters. Thus *p* = *n*_{
k
}+ *n*_{
c
}+ 1.

Model comparison has a similar flowchart. By model comparison we understand a) mapping one model to another one by model reduction or mapping each model to a third one, closest in some sense to both; b) compare predictions of the models (for instance, about how the system responds to perturbations) for sets of parameters related one to the other by the mapping at a). In this case, the choice of intermediate species is dictated by the differences between the models to be compared.

### Hierarchy of models for NF-*κ* B signalling

The transcription factor NF-*κ* B is involved in a wide diversity of domains such as the immune and inflammatory responses, cell survival and apoptosis, cellular stress and neuro-degenerative diseases, cancer and development. NF-*κ* B is sequestered in the cytoplasm by inactivating proteins named I*κ* B. Upon signalling, I*κ* B molecules are phosphorylated by a kinase complex, then ubiquitinylated, and finally degraded by the proteasomal complex. NF-*κ* B bound to I*κ* B molecules is then transported to the nucleus to activate its target genes. There are known five members of the NF-*κ* B family in mammals, Rel (c-rel), RelA (p65), RelB, NF-*κ* B1 (p50 and its precursor p105) and NF-*κ* B2 (p52 and its precursor p100). This generates a large combinatorial complexity of dimers, affinities and transcriptional capabilities. I*κ* B family comprises seven members in mammals (I*κ* B*α*, I*κ* B*β*, I*κ* B*ϵ*, I*κ* B*γ*, Bcl-3) [50]. All these inhibitors display different affinities for NF-*κ* B dimers, multiplying the combinatorial complexity. Moreover, the gene coding for I*κ* B*α*, is transcriptionally activated by NF-*κ* B. This negative feed-back loop can give rise to oscillations of the activity of NF-*κ* B [51, 52]. Phosphorylation of I*κ* B*α* upon signalling is provided by a kinases complex that includes IKK*α* and IKK*β* (I*κ* B Kinase, also named IKK1 and IKK2), associated to a regulating protein NEMO (NF-*κ* B Essential Modulator, also called IKK*γ*). Therefore, it is clear that understanding such a complex biological system requires modeling. Several mathematical models of NF-*κ* B have been published. The first model described a single NF-*κ* B molecule, which binds to I*κ* B*α*, I*κ* B*β* and I*κ* B*ϵ*. This work demonstrated oscillations in NF-*κ* B activity, confirmed by experimental data [51]. The model set by [53] included in addition an A20 molecule whose production is enhanced upon NF-*κ* B stimulation, and which negatively regulates IKK activity. A third model analyzed the critical parameters necessary for maintaining oscillations, with given amplitude and frequency [54]. In addition, a minimal simplified model was also set to study the oscillations of the NF-*κ* B module [55]. We propose here a fourth, new model, with more complex descriptions that takes into account transcription, translation and degradation of different NF-*κ* B units.

In our model, NF-*κ* B is considered to be made of two subunits, p50 and p65. All combinations of these subunits are allowed, including two homodimers p50:p50, p65:p65 and one heterodimer p50:p65. The three dimers of NF-*κ* B are characterized by different affinities for DNA sites, and associate differentially to I*κ* B*α*, *β* and *ϵ*, generating thus 9 species with different abundances and characteristics upon signalling and degradation. The production of the dimer p65 is considered under control of a transcription factor FTAy, which represents a simplification of many transcription factors supposed to activate this promoter. p50 is produced from a precursor molecule p105. The transcription factor FTAx binds to the promoter of p105 to activate its transcription at a basal level. Similarly to FTAy, this factor represents the sum of individual activities due to several transcription factors contributing to the basal activity of this promoter. As the p105 gene is activated by NF-*κ* B, this factor can also bind to the p105 promoter and activate the transcription above the basal level. Promoter of I*κ* B is controlled by NF-*κ* B and FTAz in a similar way as it is p105. In addition, it was supposed that nuclear I*κ* B can come and bind to NF-*κ* B when this is on the promoter of I*κ* B or of p105. Once the complex formed this can unbind from the DNA, taking NF-*κ* B away. The kinase activation/inactivation module including interactions with A20 was borrowed from [53]. Let us notice that transcription regulation modules are very simplified and do not take into account specificities of eukaryotic regulation (existence of several binding sites, enhancers, etc).

Initiation [56, 57] and elongation [58–60] for transcription and translation rates come from previous studies which were more recently re-examined [61]. Binding and unbinding constants for NF-*κ* B subunits come either from literature [62–64] or from previous models [51, 53].

We should signal large uncertainty concerning values of constants. For example, the rate of degradation of I*κ* B was assumed to be independent of the state of the molecule, either free or bound to NF-*κ* B. This led to a poor fit of computational simulations of the NF-*κ* B signalling module. The rate was newly measured in vivo and led to better fits of I*κ* B levels and basal NF-*κ* B activity [65]. This motivated us to determine which parameters of the model are critical and should be known with precision and which ones are not critical.

*κ*B

*α*inhibitor) is given in Table 1.

Model ℳ(39, 65, 90).

reaction | kinetic constants |
---|---|

IKK → IKK|active | k1 = 0.0025 |

IKK → null | k2 = 0.000125 |

null → IKK | k3 = 1e-005 |

IKK|active + A20 → A20 + IKK|inactive | k4 = 0.1 |

IKK|active → IKK|inactive | k5 = 0.0015 |

IKK|active → null | k6 = 0.000125 |

IKK|active + IkBa@csl → IKK|active:IkBa | k7 = 0.24 |

IKK|active:IkBa → IKK|active | k8 = 0.1 |

IKK|active + IkBa:p50:p65@csl → IKK|active:IkBa:p50:p65 | k9 = 1.2 |

IKK|active:IkBa:p50:p65 → IKK|active + p50:p65@csl | k10 = 0.1 |

IKK|inactive → null | k11 = 0.000125 |

IkBa:p50:p65@csl → p50:p65@csl | k12 = 2e-005 |

p50:p65@csl + IkBa@csl ⇌ IkBa:p50:p65@csl | kf13 = 0.5 kr13 = 0 |

p50:p65@ncl + IkBn ⇌ IkBa:p50:p65@ncl | kf14 = 0.5 kr14 = 0 |

p50:p65@csl ⇌ | kf15 = 0.0025 kr15 = 8e-005 |

mRNAA20 → mRNAA20 + A20 | k16 = 0.5 |

mRNAA20 → null | k17 = 0.0004 |

A20 → null | k18 = 0.0003 |

null → mRNAA20 | k19 = 0 |

p50:p65@ncl → p50:p65@ncl + mRNAA20 | k20 = 5e-007 |

IkBa@csl → null | k21 = 0.0001 |

mRNAIkBa → mRNAIkBa + IkBa@csl | k22 = 0.5 |

IkBa@csl ⇌ | kf23 = 0.001 kr23 = 0.0005 |

null → mRNAIkBa | k25 = 0 |

mRNAIkBa → null | k27 = 0.0004 |

| kf28 = 0.01 kr28 = 0 |

Prop105:RNAP + FTAx ⇌ Prop105:RNAP:FTAx | kf32 = 10 kr32 = 0.0001 |

Prop105:RNAP → Prop105:RNAP + RNAP1|active | k33 = 0.0005 |

Prop105:RNAP:FTAx → Prop105:RNAP:FTAx + RNAP1|active | k34 = 0.1 |

RNAP1|active → mRNAp105 | k35 = 0.01 |

mRNAp105 → mRNAp105 + p105 | k36 = 0.0041 |

mRNAp105 → null | k37 = 5e-005 |

p105 → null | k38 = 6e-005 |

p105 → p50 | k39 = 0.00013 |

p50 → null | k40 = 6.4e-005 |

FTAy + Prop65:RNAP ⇌ Prop65:RNAP:FTAy | kf41 = 10 kr41 = 0.0001 |

Prop65:RNAP → Prop65:RNAP + RNAP2|active | k42 = 0.0005 |

Prop65:RNAP:FTAy → Prop65:RNAP:FTAy + RNAP2|active | k43 = 0.1 |

RNAP2|active → mRNAp65 | k44 = 0.016 |

mRNAp65 → mRNAp65 + p65 | k45 = 0.0053 |

mRNAp65 → null | k46 = 5e-005 |

p65 → null | k47 = 6.4e-005 |

FTAz + ProIkBa:RNAP ⇌ ProIkBa:RNAP:FTAz | kf48 = 10 kr48 = 0.0001 |

ProIkBa:RNAP → ProIkBa:RNAP + RNAP3|active | k49 = 0.0005 |

ProIkBa:RNAP:FTAz → ProIkBa:RNAP:FTAz + RNAP3|active | k50 = 0.02 |

RNAP3|active → mRNAIkBa | k51 = 0.025 |

p50 + p65 ⇌ p50:p65@csl | kf52 = 0.003 kr52 = 0.001 |

p50:p65@csl → null | k53 = 0.0002 |

p50:p65@ncl → null | k54 = 0.0002 |

p50:p65@ncl + ProIkBa:RNAP ⇌ ProIkBa:RNAP:p50:p65 | kf55 = 0.62 kr55 = 0.00048 |

p50:p65@ncl + ProIkBa:RNAP:FTAz ⇌ ProIkBa:RNAP:FTAz:p50:p65 | kf56 = 0.62 kr56 = 0.00048 |

IkBn + p50:p65@nclProIkBa:RNAP ⇌ ProIkBa:RNAP:p50:p65:IkBa | kf57 = 18.4 kr57 = 0.055 |

IkBn + ProIkBa:RNAP:FTAz:p50:p65 ⇌ IkBnProIkBa:RNAP:FTAz:p50:p65 | kf58 = 18.4 kr58 = 0.055 |

ProIkBa:RNAP:p50:p65:IkBa ⇌ IkBa:p50:p65@ncl + ProIkBa:RNAP | kf59 = 0.0038 kr59 = 8e-013 |

IkBnProIkBa:RNAP:FTAz:p50:p65 ⇌ IkBa:p50:p65@ncl + ProIkBa:RNAP:FTAz | kf60 = 0.0038 kr60 = 8e-013 |

p50:p65@nclProIkBa:RNAP → p50:p65@nclProIkBa:RNAP + RNAP3|active | k61 = 0.06 |

ProIkBa:RNAP:FTAz:p50:p65 → ProIkBa:RNAP:FTAz:p50:p65 + RNAP3|active | k62 = 0.6 |

p50:p65@ncl + Prop105:RNAP ⇌ Prop105:RNAP:p50:p65 | kf63 = 0.62 kr63 = 0.00048 |

p50:p65@ncl + Prop105:RNAP:FTAx ⇌ Prop105:RNAP:FTAx:p50:p65 | kf64 = 0.62 kr64 = 0.00048 |

IkBn + Prop105:RNAP:p50:p65 ⇌ Prop105:RNAP:p50:p65:IkBa | kf65 = 18.4 kr65 = 0.055 |

IkBn + Prop105:RNAP:FTAx:p50:p65 ⇌ Prop105:RNAP:FTAx:p50:p65:IkBa | kf66 = 18.4 kr66 = 0.055 |

Prop105:RNAP:p50:p65:IkBa ⇌ IkBa:p50:p65@ncl + Prop105:RNAP | kf67 = 0.0038 kr67 = 8e-013 |

Prop105:RNAP:FTAx:p50:p65:IkBa ⇌ IkBa:p50:p65@ncl + Prop105:RNAP:FTAx | kf68 = 0.0038 kr68 = 8e-013 |

Prop105:RNAP:p50:p65 → Prop105:RNAP:p50:p65 + RNAP1|active | k69 = 0.006 |

Prop105:RNAP:FTAx:p50:p65 → Prop105:RNAP:FTAx:p50:p65 + RNAP1|active | k70 = 0.06 |

IkBa:p50:p65@csl → null | k71 = 0.0002 |

IkBa:p50:p65@ncl → null | k72 = 0.0002 |

### Model reduction

As an illustration of the model reduction flowchart, we obtain from the model proposed by Lipniacki [53] a series of simpler models. This model is ℳ(14, 25, 28) in our hierarchy: it contains only one reversible reaction and the total NF-*κ* B quantity is conserved *n*_{
c
}= *n*_{
r
}= 1. The description of the reactions can be found in Table 1 (Lipniacki's model is a submodel of our model).

The model was forced to function in a strongly oscillating regime. This situation is the most unfavorable for the simple version of Clarke's method which is doomed to shorten delays and to destabilize oscillations when intermediates are not appropriately chosen. Thus, it represents a good test for our method. First, we identify quasi-stationary and non-oscillating species. We define log-average concentration *c*_{log} = log <*c* > and the log-amplitude *a*_{log} = min(log max(*c*) - *c*_{log}, *c*_{log} - log min(*c*)) (the minimum is to avoid divergence when min(*c*) = 0). Species whose log-amplitudes are low and well separated from other values are declared non-oscillating. In order to detect quasi-stationary species, for each species *A*_{
i
}we compare two trajectories (concentrations as functions of time): a) the trajectory in the unreduced model b) the trajectory of *A*_{
i
}calculated from the trajectories of the species influencing *A*_{
i
}by using the quasi-stationarity equation (22) for *I* = {*i*}. The two trajectories must be close one to another for quasi-stationary species (except for a short transition region), see Fig. 3a). Hausdorff distance between the two trajectories can be used to detect quasi-stationary species for automatic computation. Non-oscillating species could also satisfy this criterion, but after a larger transition region, because they are slow (see the behavior of IKK|inactive in Fig. 3a)).

These procedures allow to identify 7 quasi-stationary species (IKK|active, IKK, IKK|active:IkBa, IKK|active:IkBa:p50:p65, IkBa@ncl, IkBa:p50:p65@ncl, p50:p65@csl) and one non-oscillating species (IKK|inactive). Two species with small concentration (mRNAA20, mRNAIkBa) are not quasi-stationary, as their relaxation time can be compared to the period of the oscillations. The smallness of their concentration is not a consequence of rapid consumption, but of small production (transcription) rate. Two species with large concentration are quasi-stationary (IkBa@ncl, p50:p65@csl).

*R*

_{1–11}, the complex release reaction

*R*

_{12}, the complex formation reaction

*R*

_{13}and the NF-

*κ*B translocation reaction

*R*

_{15}are replaced by two simple sub-mechanisms representing the modulated inhibitor degradation (IkBa → ∅), and summarizing the NF-

*κ*B release and translocation (IkBa:p50:p65@csl → p50:p65@ncl), respectively. The corresponding dominant rates are:

where *x*_{10} = [*IkBa*@*csl*], *x*_{8} = [*A* 20], *x*_{13} = [*IkBa* : *p* 50 : *p* 65@*csl*], *k*_{21p 1}= *k*_{3}*k*_{9}/*k*_{4}, *k*_{21}*p*_{2} = *k*_{15}*p*_{2} = *k*_{15}/*k*_{13}, *k*_{15}*p*_{2} = (*k*_{15}*k*_{3}*k*_{9})/(*k*_{4}*k*_{13}), *k*_{21}*p*_{3} = *k*_{15}*p*_{3} = *k*_{5}/*k*_{4}.

After reduction of the first module we obtain the model ℳ(8, 12, 19).

where *x*_{7} = [*p* 50 : *p* 65@*ncl*], *k*_{14p 1}= *k*_{23}, ${k}_{14p2}={k}_{v}{{k}^{\prime}}_{23}/{k}_{14}$.

This reduction step leads to the model ℳ(6, 10, 17). The dynamics (illustrated by trajectories in Fig. 3b) of the two new models is practically the same as the dynamics of the non-reduced model. One should not expect a perfect match because the method is based on asymptotic order relations between parameters. In establishing the expression of dominant rates we have considered that one parameter is much bigger than another one if the absolute value of their ratio is larger than ten. Of course, a more drastic criterion would produce more complex expressions, because less monomials could be simplified (separation of these monomials would not be large enough).

We have tested reduction of two more species that have small concentration but are not quasi-stationary. Reducing the species mRNAA20 leads to the model ℳ(5, 8, 15) Intermediate reactions (representing the transcription/translation module) are replaced by a single one (production of protein), of parameter *k*_{20p}= *k*_{16}*k*_{20}/*k*_{17}. This model has stable oscillations, but with slightly smaller period, and with different phase relations between oscillating species (A20 is almost in phase with nuclear NF*κ* B). Both period and phase changes result from the reduced delay on the negative feed-back loop containing A20. Reducing the species mRNAIkB has destabilizing effect on the oscillations. It is no longer possible to obtain self-sustained oscillations and damping times are generally smaller than for the non-reduced model. It is well known that delayed negative feed-back favors stable oscillations and that reducing the delay destabilizes oscillations. Our findings suggest that the delay along the IkBa negative feed-back loop is more important for the stability of the oscillations than the delay along the A20 loop.

Model reduction allows to identify critical and non-critical parameters. Parameters of reduced models are monomials of parameters of the non-reduced models (see Eqs.(34),(35),(36)). Some parameters of the non-reduced model may not occur in these monomials; these are non-critical parameters. Among monomials, only some are critical. Critical monomials are detected by sensitivity studies [66] performed on the reduced model. Critical parameters of the non-reduced model are contained in the critical monomials of the reduced model. The relation between critical parameters and critical monomials is hierarchical: monomials may be combined to form new monomials (in our example only two hierarchical levels are present). The degrees of critical monomials provide qualitative information on the influence of various critical parameters on the properties of the system. For instance, if two parameters have degrees of opposite signs in a critical monomial, their effects will be opposite.

*τ*be the studied quantity and

*k*the parameter (monomial). We say that

*k*is critical if ${\mathrm{sup}}_{|\mathrm{log}(s)|<A}\left|\frac{d\mathrm{log}\tau (s{k}_{0})}{d\mathrm{log}s}\right|>1$, where

*A*> 0 is some fixed constant and

*k*

_{0}some central value of the parameter. The sensitivity study is presented in Fig. 4. The relation between parameters of the initial and the reduced models is represented in Fig. 5. Damping time of the oscillations is most sensitive to parameters

*k*

_{14p 1},

*k*

_{18},

*k*

_{20p},

*k*

_{21p 1},

*k*

_{22},

*k*

_{26},

*C*

_{0}. By changing these parameters, the oscillations can be modified from damped to self-sustained. The above parameters are the critical monomials from which we get the critical parameters (with respect to damping time) of the unreduced model:

*k*

_{23},

*k*

_{18},

*k*

_{16},

*k*

_{20},

*k*

_{17},

*k*

_{3},

*k*

_{9},

*k*

_{4},

*k*

_{22},

*k*

_{26},

*C*

_{0}. The degrees of the critical monomials represent logarithmic sensitivities, therefore they provide both sign an strength of the influence of the critical parameters on the studied property. For instance, from

*k*

_{21p 1}=

*k*

_{3}

*k*

_{9}(

*k*

_{4})

^{-1}we can say that damping time can be increased (produce sustained oscillations) by reducing

*k*

_{3}, or by reducing

*k*

_{9}, or by increasing

*k*

_{4}), see also Fig. 3.

Critical parameters correspond to reactions affecting three targets: the kinase, A20, and the inhibitor, see Fig. 5. Four groups of critical monomials are easy to interpret. Increase of the monomial *k*_{20p}stands for increasing the NF-*κ* B dependent A20 production (changing *k*_{17}, *k*_{18} have the opposite effect, increase degradation). Increasing *k*_{26}, *k*_{22} stands for increasing the NF-*κ* B dependent I*κ* B production. The latter effect has been exploited in [52] to stabilize oscillations by transfecting HeLa cells with *κ* I*κ* B-EGFP vector. Decreasing *k*_{14p 1}stands for decreasing the nuclear concentration of the inhibitor, by reducing its translocation rate to the nucleus. It is possible that the experiment in [52] affected also this constant (in the right direction, ie towards decreased translocation rate) by attaching EGFP to the inhibitor. The critical monomial *k*_{21p 1}is more difficult to interpret in terms of putative targets. It gathers recovery (via *k*_{3}) and dynamical properties (via *k*_{9}), as well as the A20 dependent inactivation (via *k*_{4}) of the kinase IKK. Finally, increasing *C*_{0} means increasing the total concentration of NF-*κ* B (free or trapped).

The value of the period is remarkably robust. There are no critical monomials for the period.

Although the strongest effect on the oscillations has already been tested experimentally by increasing the NF-*κ* B dependent I*κ* B production [52], there are two remaining targets (the kinase and A20) that could be tested experimentally.

The sequence of reduction steps described above is illustrated on Fig. 1S in Additional File 1. A series of simplified models provided in SBML 2.1 [67] format and annotated by CellDesigner 3.5 [68] software are submitted to BioModels database http://www.ebi.ac.uk/biomodels/ with the following ids: MODEL7743386835, MODEL7743358405, MODEL7743315447, MODEL7743212613.

### Model comparison

To illustrate model comparison, we compare a version of our complex model (that employs only the most important member of the I*κ* B family, namely I*κ* B*α*) to the model ℳ(14, 25, 28), proposed by Lipniacki [53]. Our model is ℳ(39, 65, 90 (there are 39 species, 65 reactions, among which 18 are reversible, 6 conservation laws, though the total NF-*κ* B quantity is not conserved). The model ℳ(14, 25, 28) is a submodel of ℳ(39,65, 90) in the sense that all its species are included in our larger model. The description of our model has been sketched at the beginning of the section 3.2. A complete description is given in Table 1. To perform model comparison we define the set of intermediates *I* as the difference of the sets of species of the two models. There are 25 intermediate species and a small frontier: only 5 terminal species.

In order to verify that intermediates can be eliminated with no consequence on the dynamics we have used the method described in the previous section.

*κ*B, and min funnel production of the complex p50:p65@csl, see Fig 6. We found three categories of intermediates. There are 10 quasi-stationary species, 3 non-oscillating species and 7 buffered species (species in large excess whose concentrations are practically constant). The elimination of these is entirely justified and has no consequence on the oscillations. There are 5 non-quasistationary, oscillating species. Among these, 4 are low concentration species, representing the states of two promoters (Prop105:RNAP, PropIkBa:RNAP) free and singly occupied by transcription factors FTAx, FTAz, respectively. However, we can safely eliminate them because transcription initiation starts dominantly when both p50:p65@csl and FTAx (or FTAz) are on the promoter, therefore the non-quasistationary promoter states are not important. The last non-quasistationary, oscillating species is p50 who binds to p65 (another slow, but non-oscillating species) to produce p50:p65@csl via the min funnel. Concentrations of all quasi-stationary intermediates are small (see Fig. 7a)), (< 10

^{-4}

*μM*corresponding to less than 30 molecules per cell). The reduction that we propose is fully justified for a deterministic model, but one may ask if deterministic differential equations apply in this case. We have shown elsewhere [33] that deterministic approximation can be applied in two different situations. The first, well known situation is when the numbers of molecules are large; the law of large numbers applies. The second, less known situation, is when some species are in small numbers, but when the reactions involving these species are frequent. An example is the quick binding-unbinding of a transcription factor on a promoter site. In this case, we can consider that various states of the promoter are at stochastic equilibrium (meaning they have reached a time invariant probability distribution). Under some conditions (the intermediate reactions should be pseudo-monomolecular), stochastic averaging [69] of the remaining equations (describing the promoter activity) with respect to the invariant distribution is equivalent to applying quasi-stationarity to the fast concentrations in the deterministic approach.

*p*

_{50},

*p*

_{65}, and the mRNAIkBa. Thus, the reactions

*R*

_{41–46},

*R*

_{32–39},

*R*

_{63–70},

*R*

_{48–51},

*R*

_{55–62}are replaced by the simple submechanisms ∅ → p

_{65}, ∅ → p

_{50}, ∅ → mRNAIkBa, of rates ${{R}^{\prime}}_{45},{{R}^{\prime}}_{39},{{R}^{\prime}}_{26}$, respectively. The quasi-stationarity equations become linear after applying the strong binding, large concentration approximation for the transcription factors FTAx-y-z. The corresponding linear mechanisms ${\mathcal{R}}_{I}$ are represented in Fig. 6. The dominant solutions of the quasistationarity equations are obtained with techniques presented for linear subsystems. Using also Eq.(28) we find the following simple submechanism rates:

Where *x*_{7} = [p50 : p65@*ncl*], *x*_{11} = [*IkBa*@*ncl*], ${k}_{39p1}=\frac{{k}_{39}}{{k}_{38}+{k}_{39}}\frac{{k}_{36}{k}_{34}{k}_{68}{P}_{0}^{105}}{{k}_{37}{k}_{64}}$, ${k}_{39p2}=\frac{{k}_{39}{k}_{36}}{{k}_{38}+{k}_{39}}\frac{{{k}^{\prime}}_{66}{k}_{70}{P}_{0}^{105}}{{k}_{37}{k}_{66}}$, ${k}_{39p3}=\frac{{k}_{68}}{{k}_{64}}$, ${k}_{39p4}=\frac{{{k}^{\prime}}_{66}}{{k}_{66}}$, ${k}_{26p1}=\frac{{k}_{50}{k}_{60}{P}_{0}^{IkBa}}{{k}_{56}}$, ${k}_{26p2}=\frac{{{k}^{\prime}}_{58}{k}_{62}{P}_{0}^{IkBa}}{{k}_{58}}$, ${k}_{26p3}=\frac{{k}_{60}}{{k}_{56}}$, ${k}_{26p4}=\frac{{{k}^{\prime}}_{58}}{{k}_{58}}$. ${P}_{0}^{105}$, ${P}_{0}^{IkBa}$ are the concentrations of promoter sites of p105, IkBa, respectively.

The fourth step is a min funnel simplification of the production of the complex *p* 50 : *p* 65@*csl*.

*R*

_{40}, ${{R}^{\prime}}_{45}$,

*R*

_{47},

*R*

_{52}are replaced by ∅ → p50:p65@csl, of rate:

This leads to the model ℳ(14, 30, 41).

The fifth step, justified by averaging, introduces a new conservation law (the model ℳ(14, 30, 41) has no conservation law). Without the reactions ${{R}^{\prime}}_{52}$, *R*_{53–54}, *R*_{71–72}, that produce and consume *p* 50 : *p* 65, the total amount of *p* 50 : *p* 65 (free or in complexes with other species) would be conserved. Considering that the degradation reactions *R*_{53–54}, *R*_{71–72}, *R*_{72} have the same constant *k*, the total amount of NF-*κ* B (free, or in complexes) represented by the variable

*C*= [

*p*50 :

*p*65@

*csl*] + [

*p*50 :

*p*65@

*ncl*]/

*k*

_{ v }+ [

*p*50 :

*p*65 :

*IkBa*@

*csl*] + ... satisfies the equation :

*k*and a rapid timescale, the period of the oscillations (${{R}^{\prime}}_{52}$ is oscillating).

*C*is a slow, non-oscillating variable, it averages oscillations. Thus, the asymptotic, total amount of NF-

*κ*B is:

where the average is over a period of the oscillations.

In the fifth reduction step, reactions *R*_{52}, *R*_{53}, *R*_{54}, *R*_{71}, *R*_{72} are eliminated. Initial conditions of the system are chosen such that initial total NF-*κ* B is *C*_{0} (this is a conserved value and a new parameter of the model).

We obtain the model ℳ(14, 25, 33) that has the same species and reactions as Lipniacki's model ℳ(14, 25, 28), but slightly more parameters. The difference in the number of parameters comes from the more complex expressions of the mRNA I*κ* B transcription rate ${{R}^{\prime}}_{26}$ given by (39). In our model this rate is modulated by the nuclear I*κ* B concentration *x*_{11} (indeed, the inhibitor can unbind NF-*κ* B from DNA). This phenomenon is not taken into account in [53] where the corresponding rate is simply *R*_{26} = *k*_{26}*x*_{7}.

One important objective of model comparison is to obtain the parameter mapping. This allows to calculate the parameters of one model if the parameters of the other model are known. Then, dynamical properties of the models can be compared. In our example, all the parameters of ℳ(14, 25, 28) should be equal to the corresponding parameters of ℳ(39, 65, 90) except for *C*_{0} and *k*_{26} which are obtained by averaging (*C*_{0} is given by Eq.(42) and *k*_{26} is calculated in order to have equal average production rates of mRNAI*κ* B in the two models, see Fig. 7b). Dynamical comparison has been done in Fig. 7c. The model ℳ(14, 25, 33) is a reduction of ℳ(39, 65, 90), therefore the dynamics of these two models is very similar. The model ℳ(14, 25, 28) preserves the main features of the dynamics, except for the behavior I*κ* Ba. Without signal, the quantity of inhibitor in ℳ(14, 25, 28) is small (it is largely in excess in the other two models). With signal, the amplitude of the oscillations is higher in ℳ(14, 25, 28). These differences follow from the different kinetic laws for the transcription of I*κ* Ba. Basically, in *M*(14, 25, 33) and ℳ(39, 65, 90), I*κ* Ba has some negative influence on its own production (see Eq.(39)).

Our most complex model can account for phenomena that can not be studied by any of the conservative models ℳ(14, 25, 28), ℳ(14, 25, 33), namely it can take into account variations of the NF-*κ* B total quantity. Although this is not important in normal situations (when *C* is conserved), it could become important if one wants to cope with strong perturbations of NF-*κ* B activity.

Thus, using a more complex model depends on the experimental situation (number of variables that can be observed, or controlled, type of perturbation). The role of mathematics and modeling in quantitative biology is to predict the behavior of a system. Depending on which behavior, the simplest theory can change, and we want a hierarchy of models and model mapping methods. The process can go in both directions: reducing, or increasing the number of details.

Model mapping also allows to identify non-critical and critical parameters. Let us give only two examples. The constants or reactions 13,14 (formation of the complex) are not critical and one does not need to know them with precision. Actually, variations by a factor 100 in these constants do not change the dynamics.

The values that we use come from [51], who cites two other references [70, 71] that seem to propose very different values. Our analysis shows that this is not important. On the contrary, we show that the constant *C*_{0} (total concentration of NF-*κ* B) is a critical parameter. Reference [72] proposes 60000 molecules in a volume (of a fibroblast cell) of 2000 *μm*^{3}. This means *C*_{0} = 0.06 *μM*. Nevertheless, cell volume estimate is not really precise and errors can easily shift the model from a damped oscillatory to a sustained oscillatory dynamics. In this case, it is the comparison between model prediction and theoretical observation that can fix the value of the critical parameter.

The sequence of reduction steps described in this section is illustrated on Fig. 2S1–2S7 in Additional File 1. A series of models of decreasing complexity starting from ℳ(39, 65, 90) and upto ℳ(14, 25, 33) provided in SBML 2.1 [67] format and annotated by CellDesigner 3.5 software [68] are submitted to BioModels database http://www.ebi.ac.uk/biomodels/ with the following ids: MODEL7743656488, MODEL7743631122, MODEL7743608569, MODEL7743576806, MODEL7743528808, MODEL7743444866.

## Conclusion

We have presented a methodology for reducing and comparing systems biology models. We show how to produce a hierarchy of coarse grained models that can be used for understanding functioning of the biological systems. We show how models in the hierarchy can be mapped one onto another, thus allowing to decrease or to increase the number of details that are needed for the description of the system. Our method identifies the set of critical parameters of the system. This can be particularly useful for robustness studies (when robustness is understood as stability against parameter variability [34]) or for practical multi-target approaches in pharmacology.

We did not approach aspects of multi-scale modeling that occur in multi-organ physiology, or spatial aspects. Relation with stochastic modelling has been only briefly discussed and will be presented in detail elsewhere (Crudu et al., in preparation). The reduction methods presented here can be applied to systems of biochemical reactions modeling cell physiology [73] and can be usefully applied to various problems in signalling, metabolism, genetic regulation.

A central idea in our treatment is the hypothesis that biological systems are hierarchial, involving many separated time scales. The reduction methods were adapted to exploit this situation (we look for dominant subsystems, which lead to tremendous simplification). The hierarchical nature of the systems is not sufficiently exploited by more traditional approaches. For instance, singular perturbation copes with two time scales and eliminates the fastest. In biology, we are often interested in a "middle" time-scale, corresponding to a particular process that we study. We have shown how to eliminate both faster and slower variables. Another specificity of systems biology is the quest for critical parameters. Our approach offers naturally a solution: critical parameters are detected in the reduction process. It also extends the theory of limiting step to complex networks [1]. Showing how to find critical parameters and dominant simplifications is a first step towards a dynamical systems approach to physiology. Indeed, complex networks fulfill various tasks in simple ways by activating a few degrees of freedom. Dominant subsystems gather dynamical variables that are activated and can change when the system needs to perform a given task. Changing task could be represented as zooming in and out (change the number of degrees of freedom), or jumping laterally (change the set of degrees of freedom) in the hierarchy of models. As pointed out by Denis Noble [74], physiology should not be understood from the bottom upwards. Our approach suggests that not only the subjective understanding, but also the objective functioning of biological systems can be based on middle-out (meaning variable level of detail) pictures.

As future work we will improve our algorithms in order to propose fully automated reduction tools. At present, the automated sections of our methods are the calculation of dominant subsystems of pseudo-monomolecular subsystems and the calculation of simple sub-mechanisms stoichiometries and rates. The detection of quasi-stationary and non-oscillating species is semi-automated. The solutions of quasi-stationarity and averaged stationarity equations are not yet fully automated (except for the pseudo-monomolecular case).

We also plan to consider other applications such as high dimensional switches [75].

Concerning our model comparison methods, we would like to study hierarchies of kinetic models coming from various organisms, for which the conserved and the specific parts are the result of evolution.

## Declarations

### Acknowledgements

We acknowledge support from the French Ministry of Research program ACI IMPBio, from the British Council/French Foreign Affairs Ministry cooperation program Alliance (Partenariat Hubert Curien), from the French Complex Systems Institute ISC and EC-FP-7 (APO-SYS). AZ is member of the team "Systems Biology of Cancer "équipe labellisée par la Ligue Nationale Contre le Cancer. We thank Upinder Bhalla, Dennis Bray and John Reinitz for inspiring discussions. We also thank the students that contributed to some of the programs used in this work: Karine Yviquel, IFSIC intern and Debasis Panda, INRIA intern from IBAB.

## Authors’ Affiliations

## References

- Gorban AN, Radulescu O: Dynamic and static limitation in reaction networks, revisited. Advances in Chemical Engineering. 2008, 34: 103-173. http://arxiv.org/abs/physics/0703278View ArticleGoogle Scholar
- Lam SH, Goussis DA: The CSP Method for Simplifying Kinetics. International Journal of Chemical Kinetics. 1994, 26: 461-486.View ArticleGoogle Scholar
- Chiavazzo E, Gorban AN, Karlin IV: Comparisons of Invariant Manifolds for Model Reduction in Chemical Kinetics. Comm Comp Phys. 2007, 2: 964-992.Google Scholar
- Gorban AN, Karlin IV: Method of invariant manifold for chemical kinetics. Chem Eng Sci. 2003, 58: 4751-4768. 10.1016/j.ces.2002.12.001.View ArticleGoogle Scholar
- Gorban AN, Karlin IV: Invariant manifolds for physical and chemical kinetics, Lect. Notes. Phys. 660. 2005, Berlin, Heidelberg: SpringerGoogle Scholar
- Gorban AN, Karlin IV, Zinovyev AY: Invariant grids for reaction kinetics. Physica A. 2004, 333: 106-154. 10.1016/j.physa.2003.10.043.View ArticleGoogle Scholar
- Roussel MR, Fraser SJ: On the geometry of transient relaxation. J Chem Phys. 1991, 94: 7106-7113.View ArticleGoogle Scholar
- Krauskopf B, Osinga HM, Doedel EJ, Henderson ME, Guckenheimer J, Vladimirsky A, Dellnitz M, Junge O: A survey of method's for computing (un)stable manifold of vector fields. International Journal of Bifurcation and Chaos. 2005, 15: 763-791.View ArticleGoogle Scholar
- Auger P, de la Para RB, Poggiale JC, Sanchez E, Huu TN: Aggregation of variables and applications to population dynamics. Structured Population Models in Biology and Epidemiology, LNM 1936, Mathematical Biosciences Subseries. Edited by: Magal P, Ruan S. 2008, 209-263. Berlin: SpringerGoogle Scholar
- Gorban AN, Kazantzis N, Kevrekidis IG, Ottinger HC, : CT: Model Reduction and Coarse-Graining Approaches for Multiscale Phenomena. 2006, Berlin-Heidelberg-New York: SpringerGoogle Scholar
- Conzelmann H, Saez-Rodriguez J, Sauter T, Bullinger E, Allgower F, Gilles ED: Reduction of mathematical models of signal transduction networks: simulation-based approach applied to EGF receptor signalling. Syst Biol (Stevenage). 2004, 1 (1): 159-169.View ArticleGoogle Scholar
- Wang R, Zhou T, Jing Z, Chen L: Modelling periodic oscillations of biological systems with multiple timescales network. Syst Biol. 2004, 1: 71-84.View ArticleGoogle Scholar
- Indic P, Gurdziel K, Kronauer RE, Klerman EB: Development of a two-dimension manifold to represent high dimension mathematical models of the intracellular mammalian clock. J Biol Rhythms. 2006, 21: 222-232.View ArticlePubMedGoogle Scholar
- Borisov NM, Markevich NI, Hoek JB, Kholodenko BN: Signaling through Receptors and Scaffolds: Independent Interactions Reduce Combinatorial Complexity. Biophys J. 2005, 89: 951-966.PubMed CentralView ArticlePubMedGoogle Scholar
- Conzelmann H, Saez-Rodriguez J, Sauter T, Kholodenko BN, Gilles ED: A domain-oriented approach to the reduction of combinatorial complexity in signal transduction networks. BMC Bioinformatics. 2006, 7: 34-PubMed CentralView ArticlePubMedGoogle Scholar
- Reinhardt V, Winckler M, Lebiedz D: Approximation of slow attracting manifolds in chemical kinetics by trajectory-based optimization approaches. J Phys Chem A. 2008, 112: 1712-1718.View ArticlePubMedGoogle Scholar
- Jolliffe IT: Principal Component Analysis, Series: Springer Series in Statistics. 2002, XXIX: New York: Springer, 2Google Scholar
- Berkooz G, Holmes P, Lumley JL: The proper orthogonal decomposition in the analysis of turbulent flows. Annu Rev Fluid Mech. 1993, 25: 539-575.View ArticleGoogle Scholar
- Tresser C, Worfolk A, Baas H: Master-slave synchronization from the point of view of global dynamics. Chaos. 1995, 5: 693-699.View ArticlePubMedGoogle Scholar
- Pécou E: Splitting the dynamics of large interaction networks. J Theor Biol. 2005, 232: 375-384.View ArticlePubMedGoogle Scholar
- Schnell S, Maini PK: Enzyme Kinetics Far From the Standard Quasi-Steady-State and Equilibrium Approximations. Mathematical and Computer Modelling. 2002, 35: 137-144.View ArticleGoogle Scholar
- Bogoliubov NN, Mitropolski YA: Asymptotic Methods in the Theory of Nonlinear Oscillations. 1961, New York: Gordon and BreachGoogle Scholar
- Givon D, Kupferman R, Stuart A: Extracting macroscopic dynamics: model problems and algorithms. Nonlinearity. 2004, 17: R55-R127.View ArticleGoogle Scholar
- Acharya A, Sawant A: On a computational approach for the approximate dynamics of averaged variables in nonlinear ODE systems: Toward the derivation of constitutive laws of the rate type. J Mech Phys Sol. 2006, 54: 2183-2213.View ArticleGoogle Scholar
- Toth J, Li G, Rabitz H, Tomlin AS: The Effect of Lumping and Expanding on Kinetic Differential Equations. SIAM J Appl Math. 1997, 57: 1531-1556.View ArticleGoogle Scholar
- Clarke BL: General Method for Simplifying Chemical Networks while Preserving Overall Stoichiometry in Reduced Mechanisms. J Phys Chem. 1992, 97: 4066-4071.View ArticleGoogle Scholar
- Kruskal MD: Asymptotology. Mathematical Models in Physical Sciences. Edited by: Dobrot S. 1963, 17-48. New Jersey: Prentice-HallGoogle Scholar
- Holmes MH: Introduction to Perturbation Methods. 1995, New York: SpringerView ArticleGoogle Scholar
- Vishik MI, Ljusternik LA: Solution of some perturbation problems in the case of matrices and self-adjoint or non-selfadjoint differential equations. Russian Math Surveys. 1960, 15: 1-73.View ArticleGoogle Scholar
- Akian M, Bapat R, Gaubert S: Min-plus methods in eigenvalue perturbation theory and generalised Lidskii-Vishik-Ljusternik theorem. arXiv e-print math.SP/0402090. 2004Google Scholar
- White RB: Asymptotic Analysis of Differential Equations. 2006, London: Imperial College Press and World ScientificGoogle Scholar
- Ball K, Kurtz TG, Popovic L, Rempala G: Asymptotic analysis of multiscale approximations to reaction networks. Ann Appl Probab. 2006, 16: 1925-1961.View ArticleGoogle Scholar
- Radulescu O, Muller A, Crudu A: Théorémes limites pour des processus de Markov à sauts. Synthèse des resultats et applications en biologie moleculaire. Technique et Science Informatique. 2007, 26: 443-469. http://cat.inist.fr/?aModele=afficheN&cpsidt=18842024View ArticleGoogle Scholar
- Gorban AN, Radulescu O: Dynamical robustness of biological networks with hierarchical distribution of time scales. IET Systems Biology. 2007, 1: 238-246.View ArticlePubMedGoogle Scholar
- Glass L: Classification of biological networks by their qualitative dynamics. J Theor Biol. 1975, 54: 85-107.View ArticlePubMedGoogle Scholar
- Snoussi EH: Qualitative dynamics of piecewise-linear differential equations: a discrete mapping approach. Dyn Stab Syst. 1989, 4: 189-207.View ArticleGoogle Scholar
- de Jong H, Gouzé JL, Hernandez C, Page M, Sari T, Geiselmann J: Qualitative simulation of genetic regulatory networks using piecewise-linear models. Bull Math Biol. 2004, 66: 301-340.View ArticlePubMedGoogle Scholar
- Klamt S, Saez-Rodriguez J, Lindquist JA, Simeoni L, Gilles ED: A methodology for the structural and functional analysis of signaling and regulatory networks. BMC Bioinformatics. 2006, 7: 56-PubMed CentralView ArticlePubMedGoogle Scholar
- Temkin ON, Zeigarnik AV, Bonchev D: Chemical Reactions Networks. 1996, Boca Raton: CRC PressGoogle Scholar
- Tikhonov AN, Vasileva AB, Sveshnikov AG: Differential equations. 1985, Berlin: SpringerView ArticleGoogle Scholar
- Wasow W: Asymptotic Expansions for Ordinary Differential Equations. 1965, New York: WileyGoogle Scholar
- Fenichel N: Geometric Singular Perturbation Theory for Ordinary Differential Equations. J Diff Eq. 1979,31: 53-98.View ArticleGoogle Scholar
- Gagneur J, Klamt S: Computation of elementary modes: a unifying framework and the new binary approach. BMC Bioinformatics. 2004, 5: 175-PubMed CentralView ArticlePubMedGoogle Scholar
- Urbanczik R, Wagner C: An improved algorithm for stoichiometric network analysis: theory and applications. Bioinformatics. 2005, 21: 1203-1210.View ArticlePubMedGoogle Scholar
- Klamt S, Gagneur J, von Kamp A: Algorithmic approaches for computing elementary modes in large biochemical reaction networks. IEE Proc Syst Biol. 2005, 152: 249-55.View ArticleGoogle Scholar
- Bjorner A, Las Vergnas M, Sturmfels B, White N, Ziegler G: Oriented Matroids. 1999, Cambridge: Cambridge University Press, 2Google Scholar
- Bruno AD: Power Geometry in Algebraic and Differential Equations. 2000, Amsterdam: North-HollandGoogle Scholar
- Pontryagin LS, Rodygin LV: Approximate solution of a system of ordinary differential equations involving a small parameter in the derivatives. Soviet Math Dokl. 1960, 1: 237-240.Google Scholar
- Sari T, Yadi K: On Pontryagin-Rodygin's theorem for convergence of solutions of slow and fast systems. Electr J Diff Eq. 2004, 2004: 1-17.Google Scholar
- Ghosh S, Karin M: Missing pieces in the NF-
*κ*B puzzle. Cell. 2002, 109: S81-96.View ArticlePubMedGoogle Scholar - Hoffmann A, et al.: The I
*κ*B-NF-*κ*B signaling module: temporal control and selective gene activation. Science. 2002, 298: 1241-1245.View ArticlePubMedGoogle Scholar - Nelson DE, et al.: Oscillations in NF-
*κ*B Signaling Control the Dynamics of Gene Expression. Science. 2004, 306: 704-708.View ArticlePubMedGoogle Scholar - Lipniacki T, Paszek P, Brasier AR, Luxon B, Kimmel M: Mathematical model of NF-
*κ*B regulatory module. J Theor Biol. 2004, 228: 195-215.View ArticlePubMedGoogle Scholar - Ihekwaba AEC, et al.: Sensitivity analysis of parameters controlling oscillatory signalling in the NF-
*κ*B pathway: the roles of IKK and I*κ*Bα. Syst Biol. 2004, 1: 93-102.View ArticleGoogle Scholar - Krishna S, Jensen M, Sneppen K: Minimal model of spiky oscillations in NF-
*κ*B signaling. Proc Natl Acad Sci USA. 2006, 103: 10840-45.PubMed CentralView ArticlePubMedGoogle Scholar - Yean D, Gralla J: Transcription reinitiation rate: a special role for the TATA box. Mol Cell Biol. 1997, 17: 3809-16.PubMed CentralView ArticlePubMedGoogle Scholar
- Yie J, Senger K, Thanos D: Mechanism by which the IFN-
*β*enhanceosome activates transcription. Proc Natl Acad Sci USA. 1999, 96: 13108-13.PubMed CentralView ArticlePubMedGoogle Scholar - Dintzis HM: Assembly of the peptide chains of hemoglobin. Proc Natl Acad Sci USA. 1961, 47: 247-61.PubMed CentralView ArticlePubMedGoogle Scholar
- Jansen GMC, Moller W: Kinetic studies on the role of elongation factors 1
*β*and 1γ in protein synthesis. J Biol Chem. 1988, 263: 1773-8.Google Scholar - Narayan S, Widen SG, Beard WA, Wilson SH: RNA polymerase II transcription. J Biol Chem. 1994, 269: 12755-63.PubMedGoogle Scholar
- Wen JD, Lancaster L, Hodges C, Zeri AC, Yoshimura SG, Noller HF, Bustamante C, Tinoco I: Following translation by single ribosomes one codon at a time. Nature. 2008, 452: 598-603.PubMed CentralView ArticlePubMedGoogle Scholar
- Hart DJ, Speight RE, Cooper MA, Sutherland JD, Blackburn JM: The salt dependence of DNA recognition by NF-
*κ*B p50: a detailed kinetic analysis of the effects on affinity and specificity. Nucl Acids Res. 1999, 27: 1063-9.PubMed CentralView ArticlePubMedGoogle Scholar - de Lumley M, Hart DJ, Cooper MA, Symeonides S, Blackburn JM: A biophysical characterisation of factors controlling dimerisation and selectivity in the NF-
*κ*B and NFAT families. J Mol Biol. 2004, 339: 1059-75.View ArticlePubMedGoogle Scholar - Phelps CB, Sengchanthalangsy LL, Huxford T, Ghosh G: Mechanism of I
*κ*B*α*binding to NF-*κ*B dimers. J Biol Chem. 2000, 275: 29840-6.View ArticlePubMedGoogle Scholar - O'Dea E, Barken D, Peralta R, Tran T, Werner S, Kearns J, Levchenko A, Hoffmann A: A homeostatic model of I
*κ*B metabolism to control constitutive NF-*κ*B activity. Mol Syst Biol. 2007, 3: 1-7.View ArticleGoogle Scholar - Rabitz H, Kramer M, Dacol D: Sensitivity analysis in chemical kinetics. Annual Review of Physical Chemistry. 1983, 34: 419-461.View ArticleGoogle Scholar
- Hucka M, Finney A, Sauro HM, Bolouri H, Doyle J, Kitano H, Arkin AP, Bornstein BJ, Bray D, et al.: The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models. Bioinformatics. 2003, 19: 524-531.View ArticlePubMedGoogle Scholar
- Funahashi A, Tanimura N, Morohashi M, Kitano H: CellDesigner: a process diagram editor for gene-regulatory and biochemical networks. BIOSILICO. 2003, 1: 159-162.View ArticleGoogle Scholar
- Freidlin M, Wentzell A: Random perturbations of dynamical systems. 1984, New York: SpingerView ArticleGoogle Scholar
- Malek S, Huxford T, Ghosh G: I
*κ*B Functions through Diect Contacts with the Nuclear Localization Signal and the DNA Binding Sequences of NF-*κ*B. J Biol Chem. 1998, 273: 25427-25435.View ArticlePubMedGoogle Scholar - Carlotti F, Dower SK, Qwarnstrom EE: Dynamic Shuttling of Nuclear Factor
*κ*B between the Nucleus and Cytoplasm as a Consequence of Inhibitor Dissociation. J Biol Chem. 2000, 273: 41028-41034.View ArticleGoogle Scholar - Carlotti F, Chapman R, Dower SK, Qwarnstrom EE: Activation of nuclear factor
*κ*B in single living cells. J Biol Chem. 1999, 274: 37941-37949.View ArticlePubMedGoogle Scholar - Bray D: Protein molecules as computational elements in living cells. Nature. 1995, 376: 307-312.View ArticlePubMedGoogle Scholar
- Noble D: The Rise of Computational Biology. Nature Reviews Molecular Cell Bioloy. 2002, 3: 460-463.Google Scholar
- Bhalla US, Ram PT, Iyengar R: MAP Kinase Phosphatase as a Locus of Flexibility in a Mitogen-Activated Protein Kinase Signaling Network. Science. 2002, 297: 1018-1023.View ArticlePubMedGoogle Scholar
- Calzone L, Gelay A, Zinovyev A, Radvanyi F, Barillot E: A comprehensive modular map of molecular interactions in RB/E2F pathway. Molecular Systems Biology. 2008, 4: 174-View ArticleGoogle Scholar

## Copyright

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.