# A constraint solving approach to model reduction by tropical equilibration

- Sylvain Soliman
^{1}Email author, - François Fages
^{1}and - Ovidiu Radulescu
^{2}

**9**:24

https://doi.org/10.1186/s13015-014-0024-2

© Soliman et al.; licensee BioMed Central Ltd. 2014

**Received: **4 April 2014

**Accepted: **4 November 2014

**Published: **4 December 2014

## Abstract

Model reduction is a central topic in systems biology and dynamical systems theory, for reducing the complexity of detailed models, finding important parameters, and developing multi-scale models for instance. While singular perturbation theory is a standard mathematical tool to analyze the different time scales of a dynamical system and decompose the system accordingly, tropical methods provide a simple algebraic framework to perform these analyses systematically in polynomial systems. The crux of these methods is in the computation of tropical equilibrations. In this paper we show that constraint-based methods, using reified constraints for expressing the equilibration conditions, make it possible to numerically solve non-linear tropical equilibration problems, out of reach of standard computation methods. We illustrate this approach first with the detailed reduction of a simple biochemical mechanism, the Michaelis-Menten enzymatic reaction model, and second, with large-scale performance figures obtained on the http://biomodels.netrepository.

## Keywords

## Background

Model reduction is a central topic in systems biology and dynamical systems theory, for reducing the complexity of detailed models, finding important parameters, and developing multi-scale models for instance.

Indeed, for many of the problems in computation and analysis of complex systems, the upper limit on the size of the system that can be studied has been reached. This limit can be very low, namely tens of variables for system identification, symbolic calculation or bifurcation of attractors of large dynamical systems. For instance, the complexity of extant symbolic solvers of polynomial equations is exponential in the number of indeterminates and parameters, that sets a drastic limitation to the size of the models that can be analyzed [1],[2]. Some examples of computational difficulties that arise when trying to apply standard tools of algebraic geometry to systems biology models can be found in [3]. Model reduction is a way to bypass these limitations by replacing large scale models with models containing less parameters and variables, easier to analyse.

There are mathematical methods, based on singular perturbations or on the theory of invariant manifolds, allowing reduction of fully parametrized systems with separation of time scales. More precisely, in dissipative systems, fast variables relax rapidly to some low dimensional attractive manifold called invariant manifold [4] that carries the slow mode dynamics. A projection of dynamical equations onto this manifold provides the reduced dynamics. Numerical reduction methods, such as computational singular perturbation (CSP, [5]), intrinsic low dimensional manifold (ILDM, [6]) exploit the separation of timescales of various processes and compute approximations of the invariant manifold. Purely structural reduction methods can handle big models possibly with lack of kinetic information [7]. However, the case of biochemical models of intermediate size, with partially known parameters and that ask for symbolic analysis, is more open [8].

While singular perturbation theory is a standard mathematical tool to analyze the different time scales of a dynamical system and decompose the system accordingly, tropical methods provide a simple algebraic framework to perform these analyses systematically in polynomial systems, and in situations when model parameters are known only by their orders of magnitude. Differential equations describing kinetics of biochemical reactions are polynomial or become polynomial after decomposition of reaction mechanisms into elementary steps. For these models, quasi-equilibrium or quasi-steady state invariant manifolds allowing reductions are given by systems of algebraic equations [3]. A potentially crucial application of tropical mathematics is to enumerate and describe asymptotic solutions of algebraic systems of equations [9]. In particular, tropical solutions of polynomial equations provide the leading terms of their solutions via curves or in other terms via Newton-Puiseux series [10],[11]. At the basis of our method lies the idea that equilibration of fast variables on invariant manifolds implies, at lowest order, equilibration of at least two dominant monomials, one positive and the other negative in the right hand side of the differential equations. We have called such a condition, similar to Kapranov’s condition [11] for existence of Newton-Puiseux series with specified lowest order terms, tropical equilibration. The crux of our method lies in the computation of tropical equilibrations that define some reduced truncated systems with fewer parameters to identify, thus pointing to fewer experiments to (in)validate the model [12],[13]. Our method copes with uncertain parameters by replacing exact values by orders of magnitude and the reduction is performed symbolically in both variables and parameters. With respect to methods based on singular perturbations, this could be less precise at lowest order, but it is more general in implementation.

Solving the tropical equilibration problem boils down to solving a system of equations in the min-plus algebra (also known as the tropical semiring). For solving linear tropical systems there are pseudo-polynomial algorithms, i.e. whose complexity is polynomial in the size of the system and in the absolute values of its coefficients [14]. In the nonlinear case, the existence of tropical equilibrations, which is equivalent to the problem of the intersection of tropical varieties, was shown to be NP-complete [15]. In this paper we show that constraint-based methods, using reified constraints for expressing the equilibration conditions, make it possible to numerically solve non-linear tropical equilibration problems, out of reach of standard computation methods.

We first illustrate this approach with the detailed reduction of a simple biochemical mechanism, the Michaelis-Menten enzymatic reaction model. We detail the general procedure to obtain truncated systems by identification, through equilibration, of fast and slow species, and relate the obtained reduced systems to the usual notions of quasi-steady-state and quasi-equilibrium. Then, we demonstrate that the approach is computationally feasible, and scales up properly, by treating in an automatic way all the curated dynamical models of http://biomodels.netrepository [16].

## Model reduction by tropicalization

where *A*_{
i
}, *i*=1,…,*n* denote the chemical species and *α*_{
ji
}, *β*_{
jk
} are positive integers named stoichiometric coefficients defining which species are consumed and produced by the reaction *j*, 1≤*j*≤*r*, and in which quantities.

*x*

_{ i }, 1≤

*i*≤

*n*and read

where *k*_{
j
}>0 are kinetic parameters and we use the shorthand notation ${\mathit{x}}^{{\mathit{\alpha}}_{j}}={x}_{1}^{{\alpha}_{j1}}\dots {x}_{n}^{{\alpha}_{\mathit{\text{jn}}}}$.

where *S*_{
ij
}=*β*_{
ji
}−*α*_{
ji
} are entries of the stoichiometric matrix.

*ε*is a positive parameter much smaller than 1,

*γ*

_{ j }is an integer or, more generally, a rational number, and ${\stackrel{\u0304}{k}}_{j}$ has order unity. An approximate integer order can be obtained from any real positive parameter by

where round stands for the closest integer. Notice that in this representation, small quantities have large orders. Furthermore, the smaller *ε*, the better the separation between quantities of different orders, indeed ${lim}_{\epsilon \to 0}\frac{{k}_{i}}{{k}_{j}}=\infty $, if *γ*_{
i
}<*γ*_{
j
}.

We also define orders for species concentrations, using a vector of orders a =(*a*_{1},…,*a*_{
n
}), such that $\mathit{x}=\stackrel{\u0304}{\mathit{x}}{\epsilon}^{\mathit{a}}$. We suppose that various (*a*_{1},…,*a*_{
n
}) are integers or rational numbers with a common denominator. In our method we will calculate the concentration orders as solutions of the tropical equilibration problem (see below).

and <,> stands for the vector dot product.

*μ*

_{ j }indicate how large are these monomials, in absolute value. For sufficiently small

*ε*, monomials of different orders are well separated. For instance a monomial of smallest order ${\mu}_{j}<{\mu}_{{j}^{\prime}}$ dominates the other monomials ${k}_{j}\left|{S}_{\mathit{\text{ij}}}\right|{\mathit{x}}^{{\alpha}_{j}}\gg {k}_{{j}^{\prime}}\left|{S}_{i{j}^{\prime}}\right|{\mathit{x}}^{{\alpha}_{{j}^{\prime}}}$. One could see all these monomials as “forces” acting on the chemical species. At steady state, the resultant of all these forces should be naught. A consequence of this is that the orders of dominant positive and negative forces should be equal. This is exactly our notion of

*tropical equilibration*that we introduced in [13]. More precisely, we say the system (2) is

*tropically equilibrated*iff

The tropical equilibration problem consists in solving the system (4) for orders *a*_{
i
}, 1≤*i*≤*n*.

*α*

_{ j }are positive integers,

*γ*

_{ j }are rational powers, and

*b*

_{ j }are real coefficients. It is well known [10] that solutions of this equation can be expressed as Newton-Puiseux series, i.e. have the form

*c*

_{ i }are complex coefficients,

*a*

_{1}<

*a*

_{2}<… are integers,

*q*is a positive integer. By substituting $x\left(\epsilon \right)={c}_{1}{\epsilon}^{\frac{{a}_{1}}{q}}(1+{x}_{1}(\epsilon \left)\right)$ (where

*x*

_{1}(

*ε*) collects terms with positive orders in

*ε*) in (5) we get

*r*

_{1}(

*ε*) collects higher order terms. Necessary conditions for

*P*(

*x*,

*ε*)=0 read at lowest order

In order to satisfy (6), the minimum in (7) should be attained at least twice. Furthermore, if one looks for real solutions ${c}_{i}\in \mathbb{R}$, then from (6) it follows that at least two *b*_{
j
} corresponding to the minimum (7) should have opposite signs. This means that the lowest order *a*_{1}/*q* in the Newton-Puiseux series solution has to satisfy a tropical equilibration problem.

We must emphasize that the tropical equilibration condition is weaker than the steady state condition, and makes sense also away from a steady state. In systems with slow/fast variables, the fast variables are equilibrated by compensation of dominant forces whose orders result from the tropical condition (4). As a consequence, the fast variables can be expressed as functions of the slow variables. However, both fast and slow variables can slowly evolve under the influence of weaker, higher order forces. This picture is valid as long as the relative dominance relations between various monomial terms in Eqs.(2) are preserved. This is true if the rescaled concentrations ${\stackrel{\u0304}{x}}_{i}$ stay between bounds, whereas *ε* is allowed to tend to zero.

To summarize, the tropical equilibration is a necessary condition for the elimination of fast variables and model reduction. As we showed in [13], in order to become sufficient this condition should be combined with a boundedness condition, called permanency:

**Definition** **0.1**.

*C*

_{−}>0 and

*C*

_{+}>0, a set of renormalization exponents

*a*

_{ i }, and a function

*T*

_{0}of the initial conditions, such that after renormalization we have

## A simple example, the Michaelis-Menten reduction

where *S*,*E* *S*,*E*,*P* represent the substrate, the enzyme-substrate complex, the enzyme and the product, respectively.

where *x*_{1}= [ *S*], *x*_{2}= [ *E* *S*], *x*_{3}= [ *E*], *x*_{4}= [ *P*].

*x*

_{2}+

*x*

_{3})

^{′}=0, which implies

*e*

_{0}is a positive constant (the total amount of enzyme), and

*s*

_{0}is a positive constant (the total amount of substrate and product). These conservation laws can be used to reduce the model by elimination of

*x*

_{4}(by (10)) and

*x*

_{3}(by (9)). It follows:

The constraint *x*_{2}≤*e*_{0} resulting from the elimination is also to consider, but we will see that all equilibrations of the above equations already imply it.

### Tropical equilibration equations

By rescaling variables and parameters, we get ${x}_{i}={\stackrel{\u0304}{x}}_{i}{\epsilon}^{{a}_{i}}$, 1≤*i*≤2, ${k}_{1}={\stackrel{\u0304}{k}}_{1}{\epsilon}^{{\gamma}_{1}}$, ${k}_{-1}={\stackrel{\u0304}{k}}_{-1}{\epsilon}^{{\gamma}_{-1}}$, ${e}_{0}={\u0113}_{0}{\epsilon}^{{\gamma}_{e}}$.

*∞*plays the role of 0 and 0 plays the role of 1. The multiplicative inverse of

*a*is denoted −

*a*. Our tropical equilibration problem means solving a set of polynomial equations in this semi-field. Using these notations and properties of semi-field operation, the tropical equations become:

### Classical Michaelis-Menten reduction

*x*

_{2}for the complex concentration. Using (8) it follows that:

where *V*_{
m
}=*k*_{2}*e*_{0} is the maximum value of the production rate ${x}_{4}^{\prime}$, since *x*_{2}≤*e*_{0}.

The variable *x*_{2} satisfies equilibration relations and can be expressed as a function of a slow variable (either the substrate *x*_{1} when *x*_{2} is small, or the sum *x*_{1}+*x*_{2} in general) in two situations: quasi-stationarity and quasi-equilibrium.

*x*

_{2}that can be considered a fast species (radical). More precisely one has

*k*

_{1}

*x*

_{1}(

*e*

_{0}−

*x*

_{2})−(

*k*

_{−1}+

*k*

_{2})

*x*

_{2}=0, leading to

*x*

_{2}=

*x*

_{1}

*e*

_{0}/(

*K*

_{ m }+

*x*

_{1}), where

*K*

_{ m }=(

*k*

_{−1}+

*k*

_{2})/

*k*

_{1}, i.e.

*k*

_{1}

*x*

_{1}(

*e*

_{0}−

*x*

_{2})−

*k*

_{−1}

*x*

_{2}= 0, meaning zero net flux of the first reaction in the mechanism. This leads to

*x*

_{2}=

*x*

_{1}

*e*

_{0}/((

*k*

_{−1}/

*k*

_{1})+

*x*

_{1}), i.e.

This is justified by having a very fast transformations in the direct and reverse sense by the first reaction, much faster than the transformations by the second reaction. In this case both *x*_{1} and *x*_{2} are fast, but their sum *x*_{1}+*x*_{2} is slow.

We show next, in Section “Tropical equilibrations and model reductions”, that analysis of tropical equations provide the conditions for the asymptotic validity of quasi-stationarity and quasi-equilibrium approximations and also the exhaustive list of asymptotic regimes.

### Geometrical interpretation

*x*,

*y*) where at least two monomials of the considered polynomial are equal and larger than all the others. In logarithmic scale, this locus is made of lines, half-lines, or line segments [13],[18]. There is one tropical curve for each differential equation. For instance, the tropical curve defined by the polynomial −

*k*

_{1}

*e*

_{0}

*x*

_{1}+

*k*

_{1}

*x*

_{1}

*x*

_{2}+

*k*

_{−1}

*x*

_{2}is made of three half-lines with a common origin depicted in Figure 1, namely

*γ*

_{1}⊗

*a*

_{1})⊕

*γ*

_{−1}=(

*γ*

_{1}⊗

*a*

_{1}) and (

*γ*

_{1}⊗

*a*

_{1})⊕

*γ*

_{−1}=

*γ*

_{−1}, respectively. These are two half-lines in the plane of concentration orders:

The two branches of solutions can be also related to parts of the tropical curves corresponding to the equilibration of monomials of different signs. More precisely (21) corresponds to (18), and (22) corresponds to (20). The branch (19) of the tropical curve corresponds to the equality of two positive monomials and has no correspondence in the set of tropical equilibrations.

Similarly to computing steady states as intersections of null-clines, we are considering multiple tropical equilibrations as intersections of tropical curves.

*γ*

_{−1}⊕

*γ*

_{2}=

*γ*

_{−1}and

*γ*

_{−1}⊕

*γ*

_{2}=

*γ*

_{2}. In the first case the tropical equation (15) is equivalent to the tropical equation (14) (also, the tropical curves coincide). Therefore, the two solutions (21) and (22) equilibrate both equations. In the second case, the solutions of the tropical equation (15) form two branches, corresponding to (

*γ*

_{1}⊗

*a*

_{1})⊕

*γ*

_{2}=

*γ*

_{1}⊗

*a*

_{1}and (

*γ*

_{1}⊗

*a*

_{1})⊕

*γ*

_{2}=

*γ*

_{2}, respectively. They correspond to two half-lines in the plane of orders (

*a*

_{1},

*a*

_{2}), namely

*a*

_{2}=

*γ*

_{ e },

*a*

_{1}<

*γ*

_{2}−

*γ*

_{1}and

*a*

_{2}=

*a*

_{1}+

*γ*

_{1}+

*γ*

_{ e }−

*γ*

_{2},

*a*

_{1}>

*γ*

_{2}−

*γ*

_{1}. A simple graphical inspection of the relative positions of these half-lines with respect to the half-lines carrying solutions of the first tropical equation shows that there are four branches of tropical equilibrations:

The branch (23) equilibrates the two variables. The branch (25) equilibrates only the second variable, whereas the branches (24), (26) equilibrate only the first variable.

### Tropical equilibrations and conservation laws

This tropical problem is different from (14), (15) and leads to different solutions in general. Firstly, let us notice that elimination is not possible in semi-fields, because there is no additive inverse in general. Hence, (27), (28) (29) can not be reduced to an equivalent system of two tropical equations. Secondly, dominant monomial equilibration in the reduced model does not always correspond to monomial equilibrations in the unreduced model. A typical example is the monomial *x*_{1}*x*_{3} that becomes the difference *x*_{1}*e*_{0}−*x*_{1}*x*_{2} in the reduced model. The two monomials can equilibrate each other in the reduced model, but the same quantity is an unique, un-equilibrated monomial in the full model.

*γ*

_{−1}⊗

*a*

_{2}=

*γ*

_{−1}. In this case the two tropical equations (27), (28) are identical. Depending on

*a*

_{2}⊕

*a*

_{3}=

*a*

_{2}, or

*a*

_{2}⊕

*a*

_{3}=

*a*

_{3}we get:

These branches correspond to equilibrations of all the variables.

*γ*

_{−1}⊗

*a*

_{2}=

*γ*

_{2}the two tropical Eqs. 27, (28) are incompatible. Depending on

*a*

_{2}⊕

*a*

_{3}=

*a*

_{2}, or

*a*

_{2}⊕

*a*

_{3}=

*a*

_{3}and further choosing only one of the two tropical Eqs. 27, (28) we get the following branches:

In the branches (32), (34), the variables *x*_{2}, *x*_{3} are not equilibrated, whereas in the branches (33), (35), the variable *x*_{1} is not equilibrated.

Comparison of Eqs. (30)-(35) and Eqs. (21)-(26) proves that the tropical equations of the unreduced model have the same set of solutions as the reduced model. However, the branch of solutions (33) equilibrates all the variables in the reduced model and does not equilibrate the variable *x*_{1} in the reduced model. The reason is exactly the one given above: the monomial *x*_{1}*x*_{3} is dominant and un-equilibrated in the unreduced model, becomes *x*_{1}*e*_{0}−*x*_{1}*x*_{2} with equilibrated monomials in the reduced model.

### Tropical equilibrations and model reductions

Tropical equilibrations can be used for model reduction. The reduction starts by tropical truncation. We call tropically truncated model the model obtained by elimination of all dominated monomials from the r.h.s. of the ordinary differential equations. The next step is ordering the variables according to the values of the exponents *μ*_{
i
}−*a*_{
i
}. This allows to determine which variables are slow and fast.

An additional construction is needed in the case when the tropically truncated system of fast variables has conservation laws that are not conserved by the un-truncated system. The conservation laws define species pools that are supplementary slow variables. The pools follow differential equations involving previously dominated monomials.

For instance, in the two variables Michaelis-Menten model, we found essentially two types of reduced models, corresponding to quasi-equilibrium and quasi-stationarity approximations [19].

This truncated system has conserved quantity *z*=*x*_{1}+*x*_{2}. The variable *z* is not conserved by the full model described by (11). Indeed, addition of Eqs. (11) leads to *z*^{′}=−*k*_{−1}*x*_{2}. As the variable *x*_{2} can be eliminated from −*k*_{1}*x*_{1}*e*_{0}+*k*_{−1}*x*_{2}=0 and *x*_{1}+*x*_{2}=*z* we have the reduced dynamics *z*^{′}=−*k*_{
z
}*z*, where *k*_{
z
}=*k*_{−1}/(1+*k*_{−1}/(*k*_{1}*e*_{0})). For small *ε*, we can consider that ${k}_{z}\sim {\epsilon}^{{\gamma}_{z}}$, with *γ*_{
z
}=*γ*_{−1}− min(0,*γ*_{−1}−*γ*_{1}−*γ*_{
e
}). Because *μ*_{1}−*a*_{1}=*γ*_{1}+*γ*_{
e
}, *μ*_{2}−*a*_{2}=*γ*_{1}+*γ*_{
e
}+*a*_{1}−*a*_{2}=*γ*_{−1} the relation *k*_{
z
}>*μ*_{1}−*a*_{1}, *k*_{
z
}>*μ*_{2}−*a*_{2} are always satisfied guaranteeing that *z* is slower than *x*_{1}, *x*_{2}. The form (36) of the truncated equations and the conservation of *x*_{1}+*x*_{2} by the fast dynamics shows that this case corresponds to quasi-equilibrium of the first reaction in the Michaelis-Menten model, as described in Section “Classical Michaelis-Menten reduction”, equation 17.

These cases correspond to saturation of the enzyme (*x*_{2} has the same order as *e*_{0}). A slow variable *z*=*x*_{1}+*x*_{2} has to be introduced as before, but the reduced dynamics is *z*^{′}=−*k*_{−1}*x*_{2}=−*k*_{−1}*e*_{0}.

The variable *x*_{1} is not equilibrated, which still allows for model reduction because this variable is slow. The fast variable *x*_{2} is equilibrated, and the equilibration equation corresponds to the classical notion of quasi-stationary approximation, as described in Section “Classical Michaelis-Menten reduction”, equation 16. In this case, *μ*_{1}−*a*_{1}=*γ*_{1}+*γ*_{
e
}, *μ*_{2}−*a*_{2}=*γ*_{1}+*γ*_{
e
}+*a*_{1}−*a*_{2}=*γ*_{2}. The condition that *x*_{1} is slower than *x*_{2} reads *μ*_{1}−*a*_{1}>*μ*_{2}−*a*_{2} and we get the additional condition *γ*_{1}+*γ*_{
e
}>*γ*_{2}, which is a low enzyme concentration condition.

The variable *x*_{2} is not equilibrated, which is allowed only if this variable is slower than *x*_{1}. In this case, *μ*_{1}−*a*_{1}=*γ*_{1}+*γ*_{
e
}, *μ*_{2}−*a*_{2}=*γ*_{2}. The condition that *x*_{2} is slower than *x*_{1} reads *μ*_{1}−*a*_{1}<*μ*_{2}−*a*_{2} and leads to the additional condition *γ*_{1}+*γ*_{
e
}<*γ*_{2}, which is a high enzyme concentration condition.

The variable *x*_{1} is equilibrated, but it can not satisfy permanency. Indeed, at fixed *x*_{2} this variable will converge to zero. Therefore, this tropical equilibration does not lead to a reduced model.

## Tropical equilibration as a constraint satisfaction problem

As explained in Section “Model reduction by tropicalization”, given a biochemical reaction system with its Mass-Action kinetics, and a small *ε*, the problem of tropical equilibration is to look for a rescaling of the variables such that the dominating positive and negative term in each ODE *equilibrate* as per Eqs. (4), i.e., are of the same degree in *ε*.

Section “A simple example, the Michaelis-Menten reduction” showed that it is possible to iteratively reduce the equilibration problem to a linear system of equations for each possible pair of positive and negative dominating monomial. It is actually possible to consider fewer pairs by restricting that search to the pairs denoting edges of the Newton polygon [13]. Nevertheless, the number of linear systems to consider remains exponential in the number of species, and may lead to redhibitory computational costs, especially when handling biochemical systems with hub molecules, i.e., molecules involved in a high number of reactions (e.g., E2F, p53, cMyc in cell-cycle control or NF *κ* B in signalling), which corresponds to a Newton polygon with many vertices.

In order to tackle that complexity, we propose a numerical approach based on Constraint Programming, that encodes the equilibration problem as a Constraint Satisfaction Problem (CSP) [20]-[22] and uses reified constraints to prune the search space. Constraint Programming is a paradigm that has showed great success at solving combinatorial decision or optimization problems, in particular for real-world instances of NP-hard problems, e.g., in the field of production planning and scheduling. It is therefore an interesting way to approach the combinatorial explosion described above.

In presence of invariants (conservation laws) in the original system, Section “Conservation law constraints” has shown that some constraints related to rescaling need be added. We have shown in [23] that finding these conservation laws can be efficiently solved by constraint methods. Here we will thus assume that the conservation laws are given in input. In our prototype implementation, both the computation of conservation laws and the following equilibration are performed for a given system.

### Reified constraints

Key to the modeling of tropical equilibration problems as CSP are reified constraints. Reified constraints are special constraints that link in a bidirectional way the value of a boolean variable to the satisfaction of another constraint. They allow for powerful cuts in the search space by propagating the truth value of some constraints of the problem to the truth value of the Boolean variable, and vice versa.

states that the Boolean variable *B* is true (i.e. equal to 1) if and only if the constraints *X*=*Y* is satisfied. That constraint posts the constraint *X*=*Y* (resp. *X*≠*Y*) as soon as *B* gets value 1 (resp. 0), and vice versa, sets *B*=1 (resp. *B*=0) as soon as *X*=*Y* (resp. *X*≠*Y* i.e. when the domains of *X* and *Y* become disjoint).

For the tropical equilibration problem, these constraints are at the core of our representation of the minimum constraints as they enforce the propagation of existing knowledge before branching on the two possible values. Indeed, if *A* is the minimum of *B* and *C*, you can derive many things on the domains of *A*, *B* and *C* before eventually trying *A*=*B* or *A*=*C*. For instance it is safe to add that *A*≤*B* and *A*≤*C*, but also if you have, from other equations, that *B*≥*m* *i* *n*_{
B
} and *C*≥*m* *i* *n*_{
C
} then you can add the fact that *A*≥*m* *i* *n*(*m* *i* *n*_{
B
},*m* *i* *n*_{
C
}). If later you obtain that actually *A*=*B* then you can enforce *C*≥*B*, etc. Section “Minimum constraints” shows in more detail how reified constraint do precisely this kind of conditional addition of cuts and can therefore be used to handle minimum constraints while postponing enumerative search as much as possible.

### Variables and domains

For each original equation *d* *x*_{
i
}/*d* *t*, 1≤*i*≤*n* is introduced a variable ${a}_{i}\in \mathbb{Z}$ that is used to rescale the system by posing ${x}_{i}={\epsilon}^{{a}_{i}}\stackrel{\u0304}{{x}_{i}}$. These are the variables of our CSP. Note that they require a solver handling like for instance SWI-Prolog [24],[25] with the clpfd clpfd library by Markus Triska, which we used for our implementation.

### Tropical equilibration constraints

For each differential equation that should be equilibrated is a list of positive monomials ${M}_{i}^{+}$, and a list of negative monomials ${M}_{i}^{-}$. The degrees in *ε* of all these monomials are integer linear expressions in the *a*_{
i
}. Now, to obtain an equilibration one should enforce for each *i* that the minimum degree in ${M}_{i}^{+}$ is equal to the minimum degree in ${M}_{i}^{-}$. This corresponds to the Eqs. 4. We will see how they can be implemented with reified constraints, but for now, let us assume a constraint min(L, M)| that enforces that the variable M of is the minimum value of a list L of linear expressions over variables of . We have in our CSP, for each 1≤*i*≤*n*, min(PositiveMonomialDegrees, M) and min(NegativeMonomialDegrees, M).

### Conservation law constraints

The second kind of constraint comes from conservation laws. Each conservation law is an equality between a linear combination of the *x*_{
i
} and a constant *c*_{
i
}. By rescaling, we obtain a sum of rescaled monomials equal to ${\epsilon}^{log\left({c}_{i}\right)/log\left(\epsilon \right)}\stackrel{\u0304}{{c}_{i}}$. We want this equality to hold when *ε* goes to zero, which implies that the minimal degree in *ε* in the left hand side is equal to (the round of) the degree of the right hand side. Since once again the degrees on the left are linear combinations of our variables *a*_{
i
}, this is again a constraint of the form: min(ConservationLawDegrees, K) where K is equal to round(log(*c*_{
i
})/ log(*ε*)). This corresponds to the tropical equation (29).

### Minimum constraints

*ε*goes to zero, but at least twice. This is the case for the example treated in [12]. The constraint we need is therefore slightly more general than min/2: we need the constraint min(L, M, N) which is true if M is smaller than each element of L and equal to N elements of that list. Note that using CLP notation, we have:

*reified constraints*to propagate information between L, M and N in all directions, without enumeration. Using SWI-Prolog notations, the implementation of min/3 by reified constraints is as follows:

The translation of this predicate into words is roughly as follows, first ignoring the counts: M is smaller than a list with head H and tail T, if it is smaller than the tail T and it is smaller than the head, i.e., M ≤H. Now, we also impose that the value M is reached C times as follows: it is reached CC times in the tail and C = B + CC where B is a variable equal to 1 iff M is equal to the head and 0 otherwise. Note that this latest statement is enforced through a reified constraint, it will therefore not lead to immediate branching but to the propagation of as much information as possible (e.g., if some values in the list are already known to be strictly greater than others, the corresponding boolean for each of them will be set to 0, and thus the sum will by necessity enforce some other values to be equal to the required minimum).

This concise and portable implementation will probably improve when the minimum and min_{n} global constraints are available (see [27] for a reference). However it already proves very efficient as demonstrated in the next section.

When C is equal to one, we can fall back to using the built-in min construct in a constraint (e.g., M #= min(L1, min(L2, L3))). Some preliminary benchmarking showed that the reified version is more efficient if the length of the list is greater than 3 or 4.

### Enumeration strategy

Constraints over finite domains come with domain filtering algorithms which dynamically prune the domain of variables when the domain of other variables change in a constraint. However this strategy is not complete and must be combined with a search procedure for virtually enumerating all possible values of the variables. For this application we obtained good performances with dichotomic search by bissecting the domain of the variables (bisect option in SWI-Prolog) without any particular heuristics for choosing the variables.

Note that since this approach is numerical, contrary to solving *symbolically* an exponential but finite number of linear systems as done in Section “A simple example, the Michaelis-Menten reduction” and in [13], there can be an infinite number of solutions. This situation denotes an under-constrained linear system and remains to be interpreted biologically. In practice bounds are put on variables in order to force finiteness. This is not a restriction in practice since biochemical species’ concentrations usually do not vary by more than a hundred of magnitude orders.

Furthermore, in order to speed-up the computation of all solutions in such large domains, we used an iterative domain expansion strategy: the problem is first tried with a domain of [−2,2] for all variables, i.e., equilibrations are searched by rescaling in the 10^{−2},10^{2} interval. If that fails, the domain is doubled and the problem tried again until a limit of 10^{−128},10^{128}.

## Computation results on Biomodels.net

To benchmark our approach, we applied it systematically to all the dynamical models of the curated part of the http://biomodels.net repository [16] of biological systems, with *ε* set arbitrarily to 0.1.

*r24*from 2012-12-12 which includes 436 curated models. Among them, only 55 models have non-trivial purely polynomial kinetics (ignoring

*events*if any). Our computational results on those are summarized in Table 1, where the first column indicates whether a complete equilibration was found, and the times are in seconds.

**Number of models of the BioModels repository, with a polynomial kinetics, for which tropical equilibrations were found or not, with corresponding size of the model and computation time**

Found | # models | Variables (avg/min/max) | Time in seconds |
---|---|---|---|

(avg/min/max) | |||

yes | 23 | 17.348/3/ 86 | 0.486/0.004/2.803 |

no | 32 | 17.812/1/194 | 0.099/0.000/1.934 |

The domain expansion strategy coupled with dichotomic search by domain bisections allowed us to gain two orders of magnitude of computation time on the biggest models.

Only one of the models (number 002) used values far from 0 in the equilibration (up to *ε*^{40}) and has no complete equilibration if the domain is restricted to [−32,32]. This is because the model is written with units such that the initial concentrations are of the order 10^{−21}, translating the search accordingly. We thus do not believe that enlarging the domains even more would lead to more equilibrations. Nevertheless, choosing a smaller *ε* might increase the number of equilibrations.

18 of the 23 models for which there is a complete equilibration are actually under-constrained and appear to have an infinity of such solutions (typically linear relations between variables). For the 5 remaining ones, we computed all complete equilibrations as shown in Table 2.

**Number of equilibrations and computation time for the models of the BioModels database with finitely many numerical solutions**

Model | # variables | # equilibrations | Total time (s) |
---|---|---|---|

BIOMD0000000002 | 13 | 18 | 53 |

BIOMD0000000122 | 14 | 9 | 4.1 |

BIOMD0000000156 | 3 | 1 | <1ms |

BIOMD0000000229 | 7 | 1 | 0.07 |

BIOMD0000000413 | 5 | 5 | 0.4 |

## Conclusions

In this paper we have shown that constraint-based methods can be efficiently used to numerically solve tropical equilibration problems in biological models of real-size in the BioModels.net repository. These calulations are important for model reduction and for determining the unknown orders of the variables. Once the orders of the variables are known, the rapid variables can be identified and the system reduced to a simpler one. This truncation, described in Section “Tropical equilibrations and model reductions” coupled with the proposed constraint-based method for finding equilibrations therefore provides an automatic way to reduce models and to identify fast/slow variables. We have started the application of such technique on non-trivial models provided by biologists and modellers and hope to be able to improve both the understanding, through that identification of fast/slow variables, and the analysis, through the size reduction, of those models.

Even with the progress of high-throughput technologies, having more focused models, with fewer species and parameters to measure, will definitely permit an improvement in the quality and speed of development of the models. Furthermore, the structural methods for comparing models in model repositories, such as [7], can be refined by filtering the structural reduction relationships according to the kinetics of the reactions and the tropical reasoning on the magnitude orders.

In many cases, it makes sense biologically to only look for partial equilibrations. Strategies to decide when such decision has to be made remain unclear. Nevertheless the framework of partial constraint satisfaction and more specifically Max-CSP [28] would allow us to easily handle the maximization of the number of equilibrated variables.

One of the limits of this approach, is that it is not particularly well suited to equilibration problems with an infinite number of solutions. As discussed at the end of previous section, in such situations symbolic solutions would be more appropriate. Nevertheless, even the approximate detection of such a case by the very high number of (bounded) numerical solutions was shown to be not very costly in practice.

## Declarations

### Acknowledgements

This work has been supported by the French ANR BioTempo, CNRS Peps ModRedBio, EPIGENMED Excellence Laboratory and OSEO Biointelligence projects.

## Authors’ Affiliations

## References

- Grigoriev D, Vorobjov N:Solving systems of polynomial inequalities in subexponential time. J Symbolic Computat. 1988, 5: 37-64. 10.1016/S0747-7171(88)80005-1.View ArticleGoogle Scholar
- Grigoriev D:Complexity of quantifier elimination in the theory of ordinary differential equations. Lect Notes Comput Sci. 1989, 18: 11-25. 10.1007/3-540-51517-8_81.View ArticleGoogle Scholar
- Pantea C, Gupta A, Rawlings JB, Craciun G: The QSSA in chemical kinetics: as taught and as practiced. In Discrete and Topological Models in Molecular Biology. Berlin: Springer; 2014:419–442., Google Scholar
- Gorban A, Karlin I:Invariant manifolds for physical and chemical kinetics. Lect Notes Phys. 2005, 660: 1-491. 10.1007/978-3-540-31531-5_1.View ArticleGoogle Scholar
- Lam S, Goussis D:The CSP method for simplifying kinetics. Int J Chem Kinet. 1994, 26 (4): 461-486. 10.1002/kin.550260408.View ArticleGoogle Scholar
- Maas U, Pope SB:Simplifying chemical kinetics: intrinsic low-dimensional manifolds in composition space. Combustion Flame. 1992, 88 (3): 239-264. 10.1016/0010-2180(92)90034-M.View ArticleGoogle Scholar
- Gay S, Soliman S, Fages F:A graphical method for reducing and relating models in systems biology. Bioinformatics. 2010, 26 (18): i575-i581. [Special issue ECCB’10],View ArticlePubMedPubMed CentralGoogle Scholar
- Radulescu O, Gorban AN, Zinovyev A, Noel V:Reduction of dynamical biochemical reactions networks in computational biology. Front Genet. 2012, 3: 131-[http://www.frontiersin.org/bioinformatics_and_computational_biology/10.3389/fgene.2012.00131/abstract],View ArticlePubMedPubMed CentralGoogle Scholar
- Sturmfels B: Solving systems of polynomial equations, Volume 97, American Mathematical Soc: Providence; 2002., Google Scholar
- Walker RJ: Algebraic curves, New York: Springer; 1978., Google Scholar
- Einsiedler M, Kapranov M, Lind D:Non-archimedean amoebas and tropical varieties. J für die reine und angewandte Mathematik (Crelles J). 2006, 2006 (601): 139-157.Google Scholar
- Noel V, Grigoriev D, Vakulenko S, Radulescu O:Tropical geometries and dynamics of biochemical networks application to hybrid cell cycle models. Electron Notes Theor Comput Sci. 2012, 284: 75-91. 10.1016/j.entcs.2012.05.016.View ArticleGoogle Scholar
- Noel V, Grigoriev D, Vakulenko S, Radulescu O: Tropicalization and tropical equilibration of chemical reactions. In Tropical and Idempotent Mathematics and Applications, Volume 616 of Contemporary Mathematics. Edited by Litvinov G, Sergeev S: American Mathematical Society; 2014:261–277., Google Scholar
- Grigoriev D:Complexity of solving tropical linear systems. Comput Complexity. 2013, 22: 71-88. 10.1007/s00037-012-0053-5.View ArticleGoogle Scholar
- Theobald T:On the frontiers of polynomial computations in tropical geometry. J Symbolic Comput. 2006, 41 (12): 1360-1375. 10.1016/j.jsc.2005.11.006.View ArticleGoogle Scholar
- le Novère N, Bornstein B, Broicher A, Courtot M, Donizelli M, Dharuri H, Li L, Sauro H, Schilstra M, Shapiro B, Snoep JL, Hucka M:BioModels Database: a free, centralized database of curated, published, quantitative kinetic models of biochemical and cellular systems. Nucleic Acid Res. 2006, 1 (34): D689-D691. 10.1093/nar/gkj092.View ArticleGoogle Scholar
- Cohen G, Gaubert S, Quadrat J:Max-plus algebra and system theory: where we are and where to go now. Ann Rev Control. 1999, 23: 207-219. 10.1016/S1367-5788(99)90091-3.View ArticleGoogle Scholar
- Viro O:From the sixteenth Hilbert problem to tropical geometry. Jpn J Math. 2008, 3 (2): 185-214. 10.1007/s11537-008-0832-6.View ArticleGoogle Scholar
- Gorban AN, Radulescu O, Zinovyev AY: Asymptotology of chemical reaction networks Chem Eng Sci. 2010, 65 (7): 2310-2324. 10.1016/j.ces.2009.09.005. [International Symposium on Mathematics in Chemical Kinetics and Engineering]View ArticleGoogle Scholar
- Mackworth AK:Consistency in networks of relations. Artif Intell. 1977, 8: 99-118. 10.1016/0004-3702(77)90007-8.View ArticleGoogle Scholar
- Meseguer P:Constraint satisfaction problems: an overview. A.I. Commun. 1989, 2: 3-17.Google Scholar
- Kumar V:Algorithms for constraint- satisfaction problems: a survey. A.I. Mag. 1992, 13: 32-44.Google Scholar
- Soliman S:Invariants and other structural properties of biochemical models as a constraint satisfaction problem. Algorithms Mol Biol. 2012, 7 (15): 15-View ArticlePubMedPubMed CentralGoogle Scholar
- Wielemaker J, Schrijvers T, Triska M, Lager T: SWI-Prolog Theory Prac Logic Program. 2012, 12 (1-2): 67-96. 10.1017/S1471068411000494.View ArticleGoogle Scholar
- Wielemaker J:
*SWI-Prolog 6.3.15 Reference Manual*; 1990. [], http://www.swi-prolog.org/pldoc/refman/ - Radulescu O, Gorban A, Zinovyev A, Noel V:Reduction of dynamical biochemical reaction networks in computational biology. Front Bioinformatics Comput Biol. 2012, 3: 131-Google Scholar
- Beldiceanu N, Carlsson M, Demassey S, Petit T: Global constraints catalog. Tech. Rep. T2005-6, Swedish Institute of Computer Science 2005Google Scholar
- Freuder EC, Wallace RJ:Partial constraint satisfaction. Artif Intell. 1992, 58: 21-70. 10.1016/0004-3702(92)90004-H.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/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.