Tropical Abstraction of Biochemical Reaction Networks with Guarantees

TL;DR

This paper introduces a tropical analysis-based approximation method for biochemical reaction network models that guarantees bounds on species concentrations, simplifying complex ODE systems while maintaining accuracy.

Contribution

It presents a novel tropical abstraction technique that provides guaranteed bounds and hybrid models for high-dimensional biochemical ODE systems, improving analysis efficiency.

Findings

Method provides sound interval bounds for species concentrations.

Hybrid models adapt to changing dominance relations during system dynamics.

Tested successfully on multiple biochemical case studies.

Abstract

Biochemical molecules interact through modification and binding reactions, giving raise to a combinatorial number of possible biochemical species. The time-dependent evolution of concentrations of the species is commonly described by a system of coupled ordinary differential equations (ODEs). However, the analysis of such high-dimensional, non-linear system of equations is often computationally expensive and even prohibitive in practice. The major challenge towards reducing such models is providing the guarantees as to how the solution of the reduced model relates to that of the original model, while avoiding to solve the original model. In this paper, we have designed and tested an approximation method for ODE models of biochemical reaction systems, in which the guarantees are our major requirement. Borrowing from tropical analysis techniques, dominance relations among terms of each species' ODE are exploited to simplify the original model, by neglecting the dominated terms. As the dominant subsystems can change during the system's dynamics, depending on which species dominate the others, several possible modes exist. Thus, simpler models consisting of only the dominant subsystems can be assembled into hybrid, piecewise smooth models, which approximate the behavior of the initial system. By combining the detection of dominated terms with symbolic bounds propagation, we show how to approximate the original model by an assembly of simpler models, consisting in ordinary differential equations that provide time-dependent lower and upper bounds for the concentrations of the initial models species. Our method provides sound interval bounds for the concentrations of the chemical species, and hence can serve to evaluate the faithfulness of tropicalization-based reduction heuristics for ODE models of biochemical reduction systems. The method is tested on several case studies.

Equations

$$r_j:\sum _i{\alpha }_{ji}{x}_i\stackrel{k_j^{+},k_j^{-}}{\leftrightarrow }\sum _i{\beta }_{ji}{x}_i,$$

$$\frac{d\text{x}}{dt}=\sum _{j=1}^r({\mathbf{\beta }}_j-{\mathbf{\alpha }}_j)(R_j^{+}(\text{x})-R_j^{-}(\text{x})).$$

$$\begin{array}{rl}& R_j^{+}(\text{x})=k_j^{+}\prod _i{x}_i^{{\alpha }_{ji}}\\ & R_j^{-}(\text{x})=k_j^{-}\prod _i{x}_i^{{\beta }_{ji}},\end{array}$$

$$\frac{d{x}_i}{dt}=P_i^{+}(\text{x})-P_i^{-}(\text{x}),$$

$$\begin{array}{rl}& P_i^{+}(\text{x})=\sum _{{\beta }_{ji}-{\alpha }_{ji}>0}({\beta }_{ji}-{\alpha }_{ji})R_j^{+}(\text{x})+\sum _{{\beta }_{ji}-{\alpha }_{ji}<0}({\alpha }_{ji}-{\beta }_{ji})R_j^{-}(\text{x})\\ & P_i^{-}(\text{x})=\sum _{{\beta }_{ji}-{\alpha }_{ji}<0}({\alpha }_{ji}-{\beta }_{ji})R_j^{+}(\text{x})+\sum _{{\beta }_{ji}-{\alpha }_{ji}>0}({\beta }_{ji}-{\alpha }_{ji})R_j^{-}(\text{x})\end{array}$$

$$\frac{d{x}_i}{dt}=Dom(P_i^{+}(\text{x}))-Dom(P_i^{-}(\text{x}))$$

$$\mathrm{E}+\mathrm{S}\underset{k_{-1}}{\overset{k_1}{\rightleftharpoons }}\mathrm{E}:\mathrm{S}\stackrel{k_2}{\to }\mathrm{E}+\mathrm{P}$$

$$\{\begin{array}{l}\frac{d[S]}{dt}=k_{-1}[E:S]-k_1[E][S]\\ \frac{d[E]}{dt}=k_{-1}[E:S]+k_2[E:S]-k_1[E][S]\\ \frac{d[E:S]}{dt}=k_1[E][S]-k_{-1}[E:S]-k_2[E:S]\\ \frac{d[P]}{dt}=k_2[E:S]\end{array}$$

$$v=\frac{d[P]}{dt}=k_2\frac{[E{]}_T[S]}{K_M+[S]},$$

$$\mathrm{S}\stackrel{k_2\frac{[E{]}_T}{K_M+[S]}}{\to }\mathrm{P}$$

$$\mathrm{S}\stackrel{k_2\frac{[E{]}_T}{K_d+[S]}}{\to }P,$$

$$\{\begin{array}{ll}{k}_{-1}[E:S]-k_1[E][S]& \le \frac{d[S]}{dt}\le k_{-1}[E:S]-k_1[E][S]\\ k_2[E:S]-k_1[E][S]& \le \frac{d[E]}{dt}\le (1+ϵ)k_2[E:S]-k_1[E][S]\\ k_1[E][S]-(1+ϵ)k_2[E:S]& \le \frac{d[E:S]}{dt}\le k_1[E][S]-k_2[E:S]\\ k_2[E:S]& \le \frac{d[P]}{dt}\le k_2[E:S]\end{array}$$

$$\{\begin{array}{ll}{F}_{\underline{x_1}}^{\#}(\underline{\overline{x}})=& k_{-1}\underline{x_3}-k_1\overline{x_2}\underline{x_1}\\ F_{\overline{x_1}}^{\#}(\underline{\overline{x}})=& k_{-1}\overline{x_3}-k_1\underline{x_2}\overline{x_1}\\ F_{\underline{x_2}}^{\#}(\underline{\overline{x}})=& k_2\underline{x_3}-k_1\underline{x_2}\overline{x_1}\\ F_{\overline{x_2}}^{\#}(\underline{\overline{x}})=& (1+ϵ)k_2\overline{x_3}-k_1\overline{x_2}\underline{x_1}\\ F_{\underline{x_3}}^{\#}(\underline{\overline{x}})=& k_1\underline{x_1}\underline{x_2}-(1+ϵ)k_2\underline{x_3}\\ F_{\overline{x_3}}^{\#}(\underline{\overline{x}})=& k_1\overline{x_1}\overline{x_2}-k_2\overline{x_3}\\ F_{\underline{x_4}}^{\#}(\underline{\overline{x}})=& k_2\underline{x_3}\\ F_{\overline{x_4}}^{\#}(\underline{\overline{x}})=& k_2\overline{x_3}\end{array}$$

$$\underline{x_i}(t)\le x_i(t)\le \overline{x_i}(t)$$

$$\{\begin{array}{l}\frac{d{x}_2}{dt}\le (1+ϵ)k_2{x}_3-k_1{x}_2{x}_1\\ \frac{d{x}_2}{dt}\le (1+ϵ)k_2(e_0-x_2)-k_1{x}_2{x}_1\\ \frac{d{x}_2}{dt}\le (1+ϵ)k_2(s_0-x_1-x_4)-k_1{x}_2{x}_1\le (1+ϵ)k_2(s_0-x_1)-k_1{x}_2{x}_1\end{array}$$

$$\{\begin{array}{l}\frac{d\overline{x_{2,1}}}{dt}=(1+ϵ)k_2\overline{x_3}-k_1\overline{x_2}\underline{x_1}\\ \frac{d\overline{x_{2,2}}}{dt}=(1+ϵ)k_2(e_0-\overline{x_2})-k_1\overline{x_2}\underline{x_1}\\ \frac{d\overline{x_{2,3}}}{dt}=(1+ϵ)k_2(s_0-\underline{x_1})-k_1\overline{x_2}\underline{x_1}\end{array}$$

$$\min (\frac{d\overline{x_{2,1}}}{dt},\frac{d\overline{x_{2,2}}}{dt},\frac{d\overline{x_{2,3}}}{dt})=(1+ϵ)\cdot k_2\cdot \min (\overline{x_3},e_0-\overline{x_2},s_0-\underline{x_1})-k_1\overline{x_2}\underline{x_1}$$

$$\frac{d\overline{x_2}}{dt}=(1+ϵ)\cdot k_2\cdot \min (\overline{x_3},e_0-\overline{x_2},s_0-\underline{x_1})-k_1\cdot \max (\underline{x_1},0)\cdot \max (\overline{x_2},e_0-\overline{x_3})$$

$$\{\begin{array}{l}M+M\underset{k_{-1}}{\stackrel{k_1}{\rightleftharpoons }}{M}_2\\ M_2+DNA\underset{k_{-2}}{\stackrel{k_2}{\rightleftharpoons }}{M}_2.DNA\\ M_2.DNA\stackrel{k_3}{\to }DNA+P\end{array}$$

$$\{\begin{array}{l}\frac{d{x}_1}{dt}=-2k_1{x}_1^2+2k_{-1}{x}_2\\ \frac{d{x}_2}{dt}=-k_{-1}{x}_2-k_2{x}_2{x}_3+k_{-2}{x}_4+k_1{x}_1^2\\ \frac{d{x}_3}{dt}=k_{-2}{x}_4+k_3{x}_4-k_2{x}_2{x}_3\\ \frac{d{x}_4}{dt}=-k_{-2}{x}_4-k_3{x}_4+k_2{x}_2{x}_3\\ \frac{d{x}_5}{dt}=k_3{x}_4\end{array}$$

$$\{\begin{array}{l}{x}_1+2x_2+2x_4+2x_5=M_0\\ x_3+x_4=DN{A}_0\end{array}$$

$$\forall x_i\in \mathcal{S},\mathbb{F}(x_i)=P_i^{+}(\text{x})-P_i^{-}(\text{x}),$$

$$\{\begin{array}{ll}f(x_i)(0)=\mathbb{I}(0),& \forall x_i\in \mathcal{S}\\ \frac{df(x_i)}{dt}(t)=[[\mathbb{F}(x_i){]]}_{\mathcal{S}[x_i\mapsto f(x_i)(t)]},& \forall x_i\in \mathcal{S}\text{ and }t\in {\mathbb{R}}^{+}\end{array}$$

$$f(x_i)(t)=[[{\mathbb{E}}_b(x_i){]]}_{\mathcal{S}[x_i\mapsto f(x_i)(t)]}$$

$$\{\begin{array}{ll}{\mathbb{E}}_1(x_1)=M_0-2x_2-2x_4-2x_5;& {\mathbb{E}}_2(x_1)=M_0-2DN{A}_0-2x_2+2x_4-2x_5\\ {\mathbb{E}}_1(x_2)=\frac{M_0-x_1}{2}-x_4-x_5;& {\mathbb{E}}_2(x_2)=\frac{M_0-x_1}{2}-DN{A}_0+x_3-x_5\\ {\mathbb{E}}_1(x_3)=DN{A}_0-x_4;& {\mathbb{E}}_2(x_3)=DN{A}_0-\frac{M_0-x_1}{2}+x_2+x_5;\\ {\mathbb{E}}_1(x_4)=DN{A}_0-x_3;& {\mathbb{E}}_2(x_4)=\frac{M_0-x_1}{2}-x_2-x_5\\ {\mathbb{E}}_1(x_5)=\frac{M_0-x_1}{2}-x_2-x_4;& {\mathbb{E}}_2(x_5)=\frac{M_0-x_1}{2}-DN{A}_0+x_3-x_2\end{array}$$

$$r_i^{k,l}:=\{\begin{array}{ll}(M_k^{+},M_l^{-}),& \text{ if }k\le p,l\le n,Dom(P_i^{+})=M_k^{+},Dom(P_i^{-})=M_l^{-},\\ (P_i^{+},M_l^{-}),& \text{ if }k=p+1,l\le n,Dom(P_i^{-})=M_l^{-}\\ (M_k^{+},P_i^{-}),& \text{ if }k\le p,l=n+1,Dom(P_i^{+})=M_k^{+}\\ (P_i^{+},P_i^{-}),& \text{ if }k=p+1,l=n+1\end{array}$$

$${\mathbb{F}}_{\downarrow }^{k,l}(x_i):=\{\begin{array}{ll}{M}_k^{+}\dot{-}(1\dot{+}(n\dot{-}1)\dot{\cdot }ϵ)\dot{\cdot }{M}_l^{-},& \text{ if }k\le p,l\le n,Dom(P_i^{+})=M_k^{+},Dom(P_i^{-})=M_l^{-}\\ P_i^{+}\dot{-}(1\dot{+}(n\dot{-}1)\dot{\cdot }ϵ)\dot{\cdot }{M}_l^{-},& \text{ if }k=p+1,l\le n,Dom(P_i^{-})=M_l^{-}\\ M_k^{+}\dot{-}{P}_i^{-},& \text{ if }k\le p,l=n+1,Dom(P_i^{+})=M_k^{+}\\ P_i^{+}\dot{-}{P}_i-,& \text{ if }k=p+1,l=n+1\end{array}$$

$${\mathbb{F}}_{\uparrow }^{k,l}(x_i):=\{\begin{array}{ll}(1\dot{+}(p\dot{-}1)\dot{\cdot }ϵ)\dot{\cdot }{M}_k^{+}\dot{-}{M}_l^{-},& \text{ if }k\le p,l\le n,Dom(P_i^{+})=M_k^{+},Dom(P_i^{-})=M_l^{-}\\ P_i^{+}\dot{-}{M}_l^{-},& \text{ if }k=p+1,l\le n,Dom(P_i^{-})=M_l^{-}\\ (1\dot{+}(p\dot{-}1)\dot{\cdot }ϵ)\dot{\cdot }{M}_k^{+}\dot{-}{P}_i^{-},& \text{ if }k\le p,l=n+1,Dom(P_i^{+})=M_k^{+}\\ P_i^{+}\dot{-}{P}_i-,& \text{ if }k=p+1,l=n+1\end{array}$$

$$\{\begin{array}{l}{\mathbb{F}}_{\downarrow }^{2,1}(x_2)=k_{-2}{x}_4-(1+ϵ)k_{-1}{x}_2\\ {\mathbb{F}}_{\uparrow }^{2,1}(x_2)=(1+ϵ)k_{-2}{x}_4-k_{-1}{x}_2\end{array},$$

$$\forall r_i^{k,l},\forall x_j\in \mathcal{V},\{\begin{array}{l}{\mathbb{L}}_{i,b}^{k,l}(x_j):=\{\begin{array}{ll}{\mathbb{E}}_b(x_j),& \text{ if }\mathcal{V}={\mathcal{V}}_b\\ 0,& \text{ otherwise}\end{array}\\ {\mathbb{U}}_{i,b}^{k,l}(x_j):=\{\begin{array}{ll}{\mathbb{E}}_b(x_j),& \text{ if }\mathcal{V}={\mathcal{V}}_b\\ [[E_b(x_j){]]}_{{\mathcal{V}}_b\setminus \mathcal{V}[x_j\mapsto b_i^{k,l}(x_j)]},& \text{ otherwise }\end{array}\end{array}$$

Full text

In this paper, we have designed and tested an approximation method for ODE models of biochemical reaction systems, in which the guarantees are our major requirement. Borrowing from tropical analysis techniques, we look at the dominance relations among terms of each species’ ODE. These dominance relations can be exploited to simplify the original model, by neglecting the dominated terms. As the dominant subsystems can change during the system’s dynamics, depending on which species dominate the others, several possible modes exist. Thus, simpler models consisting of only the dominant subsystems can be assembled into hybrid, piecewise smooth models, which approximate the behavior of the initial system. By combining the detection of dominated terms with symbolic bounds propagation, we show how to approximate the original model by an assembly of simpler models, consisting in ordinary differential equations that provide time-dependent lower and upper bounds for the concentrations of the initial model’s species.

The utility of our method is twofold. On the one hand, it provides a reduction heuristics that performs without any prior knowledge of the initial system’s behavior (i.e., no simulation of the initial system is needed in order to reduce it). On the other hand, our method provides sound interval bounds for each species, and hence can serve to evaluate the faithfulness of tropicalization reduction heuristics for ODE models of biochemical reduction systems. The method is tested on several case studies.

††volume: NN††journal: Electronic Notes in Theoretical Computer Science††thanks: Email: \normalshapebeica@di.ens.fr††thanks: Email: \normalshapeferet@di.ens.fr††thanks: Email: \normalshapetatjana.petrov@uni-konstanz.de, Tatjana Petrov’s research is supported by the Ministry of Science, Research and the Arts of the state of Baden-Württemberg

1 Introduction

As biology becomes a data intensive science, due to advancements in high throughput molecular biology, the importance of in silico dynamical models of complex biological systems that are able to reproduce intricate behaviors observed in experimental settings increases. As such, modeling becomes a part of biological reasoning, but turns out to be a particularly challenging task. In particular, in models of biochemical networks, the number of possible chemical species is often subject to combinatorial explosion, due to the large number of species that may arise as a result of protein bindings and post-translational modifications [17]. As a consequence, mechanistic models of signaling pathways easily become very combinatorial. A common modeling approach describes the time-dependent evolution of concentrations of each of the modeled species through a system of coupled ordinary differential equations (ODEs). The combinatorial explosion of species and rich interaction scheme renders solving such a system of ODEs often prohibitory in practice, let alone the fact that is it already an approximation of its stochastic counter-part [20], as well as that the equations themselves do not transparently reflect the underlying mechanisms. Addressing the latter, formalisms allowing to write the mechanistic hypothesis in form of discrete transition steps have been proposed: Boolean networks[30], logical networks[32], Petri Nets[8], cellular automata[16], rule-based languages[10], to name the most common. Languages such as Kappa[10, 11] and BNGL[3] provide compact ways of describing models prone to combinatorial explosion, of simulating them [9], and even transforming them into ODEs [3]. However, the curse of dimensionality once again rises when trying to compute the system behavior.

A strategy to cope with such complexity is model reduction, in which certain properties of biochemical models are exploited in order to obtain simpler versions of the original complex model; these simpler models should preserve the important behavioral aspects of the initial system. An example of such a property is the multiscaleness of biochemical networks, with respect to both time-scales and species’ abundance. In the case of the former, it is known that biochemical processes governing network dynamics span over many well separated timescales: while protein complex formation occurs on the seconds scale, post-translational protein modification takes minutes, and changing gene expression can take hours, or even days. As for the latter, multiscaleness also applies to the abundance of various species in biochemical networks: the DNA molecule has one to a few copies, while mRNA copy numbers can vary from a few to tens of thousands. On the one hand, these widely different time- and concentration scales represent challenges for the estimation of rate constants, for the measurement of low-concentration species, and even for the numerical integration. On the other hand, they represent a feature that can be exploited for model reduction purposes, allowing to approximate the complete mechanistic description with simpler rate expressions, retaining the essential features of the full problem on the time scale or in the concentration range of interest. The dynamics of multiscale, large biochemical systems can be reduced to those of simpler models, called dominant subsystems [26], which contain less parameters and are easier to analyze. Dominant subsystems are chosen by comparing the time-scales of the large system. For example, the classical quasi steady-state (QSS) [5] and quasi-equilibrium (QE) approximations [23, 18] are conditions that lead to dominance, and represent popular methods for the computation of “first approximations” to the slow invariant manifold. Classical QSS is based on the small concentrations of highly reactive intermediate species (i.e., atoms, ions, enzymes and substrate-enzyme complexes)[4], while in the QE approximation the reduction of the full mechanism is done based on the existence of fast and slow reactions.

The multiscaleness property of biochemical network is by definition closely linked to the mathematical notion of dominance, captured in the framework of tropical analysis[1, 21]. Recently, a class of semi-formal methods for reducing and hybridizing models of biochemical networks has been developed, based on ideas from tropical analysis [26, 27, 28, 29]. These methods exploit the multiscaleness of biochemical networks, in order to deduce dominance relations among parameters and/or reaction rates, which can then be used to obtain a system of truncated ODEs (by eliminating the dominated terms). One of the advantages of using dominance relations in multi-scale networks is that it helps cope with parameter uncertainty: parameter values are replaced with their orders of magnitude, which are easier to determine. However, providing guarantees as to how the solution of the reduced model relates to the original one remains a challenge. Such is also the case in [2], where we proposed a reduction framework for rule-based models, based on time-scale separation. While this time-scale separation technique is justified by asymptotic convergence results, for any concrete parameter values, there is no information on the accuracy of the trajectories obtained by executing the reduced model.

In this paper, we design an approximation method for ODE models of biochemical networks, in which the guarantees are our major requirement. Our method combines abstraction and numerical approximation, and aims at providing better understanding/evaluation of tropical reduction methods. We abstract the solution of the original system of ODEs by a box, that over-approximates the state of the original system and provides lower and upper bounds for the value of each variable of the system in its current state. The simpler equations (which we call tropicalized) defining the hyperfaces of the box are obtained by combining the dominance concept borrowed from tropical analysis with symbolic bounds propagation. Mass invariants of the initial system of ODEs are used to refine the computed bounds, thus improving the accuracy of the method. The resulting (simplified) system provides a posteriori time-dependent lower and upper bounds for the concentrations of the initial model’s species, and thus bounds on numerical errors stemming from tropicalization. This means that no information on the original system’s trajectory is needed - the most important advantage of our approach. By contrast, the main difficulty of applying the classical QSS and QE reductions to biochemical models is that QE reactions and QSS species need to be specified a priori, which implies that some knowledge about the initial system’s behavior is necessary. This, in turn, means that significantly high-dimensional, non-linear systems cannot benefit from these reductions, as their analysis can be prohibitive in practice. An approach similar with respect to providing a posteriori time-dependent lower and upper bounds has been proposed in [7], where the differential semantics of rule-based models with non-contracting dynamics and unbounded sets of variables are treated. Rather than using dominance relations between ODE terms, a finite set of patterns is used in order to bound the number of occurrences of each pattern. Further related works, similar in the sense that they provide automatizible reduction methods with strong reduction guarantees are described in [13, 14]. However, both of these works are designed specifically for rule-based models, where they exploit the site-graph encoding of species’ structure, rather than the dominance regions.

Depending on the chosen granularity of mass-invariant-derived bounds, the method presented in this paper can either be used to reduce models of biochemical networks, or to quantify the approximation error of tropicalization reduction methods that do not involve guarantees. The guarantees of our method are obtained by formalizing the soundness relation between the original system of equations and the abstract system of ordinary differential equations operating on the coordinates of the hyper-faces of the box. The solution of a sound abstraction of an original system of differential equations, starting from a box that contains the initial state of the original system, defines a sound abstraction of the solution(s) of the original system. We apply our method to several case studies.

Outline The rest of this article is organized as follows. In Sect.2 we define the setting and concepts used in our approach, as well as introduce motivating examples. We then formally present and justify the method for deriving the system of reduced ODEs over the lower and upper bounds of species’ concentrations in Sect. 3. Also in Sect. 3, we present the two possible uses of our approach. We then discuss and conclude in Sect. 4.

2 Definitions and Motivating Example(s)

2.1 General Setting and Definitions

We define a dynamic reaction network over a set of species $\mathcal{S}=\{x_{1},\ldots,x_{s}\}$ as a reaction system of the form:

[TABLE]

where $1\leq j\leq r$ is the reaction number, and for each reaction $r_{j}$, $k_{j}^{+/-}$ are the non-negative reaction rate constants of the forward, respectively backward reaction.

In this paper, we focus on the ordinary differential equations (ODE) semantics of models of biochemical networks. The underlying assumptions are that the various species of the chemical network are highly abundant, that stochastic fluctuations are negligible, and that the reaction system is well-mixed. In these conditions, the state of the dynamic system (1) can be represented as a multiset of the species’ concentrations $\textbf{x}\in\mathbb{R}^{n}$, and its dynamics is described by a system of ODEs:

[TABLE]

For each reaction $r_{j}$, $R^{+/-}_{j}(\textbf{x})$ denotes the reaction (forward and backward) rate, and is a non-linear function of the concentrations. For example, the mass-action law reads:

[TABLE]

which in turn means that $\frac{d\textbf{x}}{dt}$ is a multivariate polynomial of the species’ concentrations. In other words, the ODE of the $i$-th species, $x_{i}$, under mass action kinetics reads as a sum of monomials, which can be split into production and consumption terms, according to the sign that precedes their occurence in the equation:

[TABLE]

where $P_{i}^{+/-}(\textbf{x})$ are Laurent polynomials with positive coefficients:

[TABLE]

For convenience purposes, we will denote $P_{i}^{+/-}(\textbf{x})=\sum\limits_{j}M_{i,j}^{+/-}(\textbf{x})$, where $M_{i,j}^{+/-}(\textbf{x})$ represent the production, respectively the consumption, monomials.

The reduction heuristics that use ideas from tropical analysis exploit the concept of dominance, which we borrow for our method. Let $M_{1}(\textbf{x})=c_{1}\textbf{x}^{\alpha_{1}}$ and $M_{2}(\textbf{x})=c_{2}\textbf{x}^{\alpha_{2}}$ be two (positive) monomials. We define $\epsilon$-dominance as the following partial order relation on the set of multivariate monomials defined on subsets of $\mathbb{R}_{+}^{n}$:

Definition 2.1**.**

($\epsilon$-dominance) For an $\epsilon\in[0,1]$, we say that $M_{1}$ dominates $M_{2}$ at a time point $t$, denoted by $M_{1}\succ_{\epsilon}M_{2}$, if $\epsilon\cdot M_{1}(\textbf{x}(t))\geq M_{2}(\textbf{x}(t))$.

In multiscale biochemical systems, the various monomials that compose the polynomials $P_{i}^{+/-}$ have different magnitude orders, such that at any given time there is only one or a few dominating monomials.

Definition 2.2**.**

(Dominant monomial of a polynomial) For a given $\epsilon\in[0,1]$, the dominant monomial of a polynomial $P_{i}(\textbf{x})=\sum\limits_{j=1}^{n}M_{i,j}(\textbf{x})$ is defined as $Dom(P_{i})=\{M_{i,j}\mid\forall 1\leq k\leq n,j\neq k,M_{i,j}\succ_{\epsilon}M_{i,k}\}$.

By using the max-plus algebra idea that the sum of positive, well separated terms, can be replaced by the maximum term, each of the two polynomials of (4) can be replaced by their dominant monomials. The result is a reduced model, consisting of a piecewise smooth function. As the dominant monomials of the $P_{i}^{+/-}$ can change from one concentration domain to another, the reduced model is a piecewise-smooth hybrid model.

Definition 2.3**.**

(Two-term tropicalization of the smooth ODE system) We call two-term tropicalization of the smooth ODE system (4) the following piecewise-smooth system:

[TABLE]

We note that a one-term tropicalization of the smooth ODE system, $Dom(\frac{dx_{i}}{dt})$, is also possible, but choosing only one dominant momomial instead of the production-consumption pair of dominant monomials leads to a less precise model reduction (as more information is discarded in the one-term tropicalization). Thus, in this paper, we choose to deal with the two-term method.

2.2 Motivating example: Michaelis-Menten

In Sect.1, we mentioned that the classical QSS[5] and QE[18, 23] approximations represent popular methods for the simplification of biochemical networks.

As such, our motivating example is the Michaelis-Menten mechanism. This enzymatic model illustrates how these two simple methods that use the idea of dominance can be useful for model reduction of nonlinear models with multiple timescales. The Michaelis-Menten mechanism consists of an enzyme that reversibly binds a substrate to form a complex, which in turn releases a product, while preserving the enzyme:

[TABLE]

The ODE system describing the evolution of the species’ concentration writes as:

[TABLE]

The Michaelis-Menten equation relates the rate of product formation to the substrate concentration:

[TABLE]

where $[E]_{T}=[E]+[E:S]$ is the total enzyme concentration, and $K_{M}=\frac{k_{2}+k_{-1}}{k_{1}}$ is the Michaelis-Menten constant. Eq.(9)can be interpreted as the reaction rate of a reduced reactive system, equivalent to system (7), in which the intermediary complex $[E:S]$ has been eliminated:

[TABLE]

The approximation of (8) by (9) is generally considered to be sufficiently good if the QSS assumption holds, that is if the total initial enzyme concentration is much lower than the total initial concentration of substrate: $[S]_{0}\succ_{\epsilon}[E]_{0}+[E:S]_{0}$. In this case, the complex $[E:S]$ is a low concentration fast species, whose concentration is dominated by that of the substrate; the value of $[E:S]$ almost instantly relaxes to a value determined by $[S]$. Thus, one can set $\frac{d[E:S]}{dt}=0$, and exploit this relation to pool the two reactions of the initial system (7) into an unique irreversible reaction (10). The QSS condition can also be stated as $k_{2}\gg k_{-1}$ [26].

The original MM analysis used the complementary QE approximation, which considers the complex formation reaction to be fast and reversible: $k_{-1}\gg k_{2}$. Thus, the term $k_{-1}[E:S]$ dominates the term $k_{2}[E:S]$ in Eq.(8), meaning the latter can be discarded from the ODE system, allowing for pooling of species, and resulting once again in a single step approximation that reads:

[TABLE]

with $K_{d}=\frac{k_{-1}}{k_{1}}$, if indeed $[S]\gg[E]+[E:S]$. We note that if the QE assumption is indeed valid, $K_{M}\approx K_{d}$.

One of the main difficulties of applying QSS and QE reductions to biochemical models is that the QE reactions and QSS species need to be specified a priori. Thus, simulation of the original model is sometimes 111In [27], the authors propose a formal method for the identification of QSS species and QE conditions, which follows from the calculation of the tropicalized system, and which does not require simulation of trajectories. needed in order to detect dominated species, which are either QSS species, or participate in QE reactions [26]. For high-dimensional non-linear systems, this requirement can represent an obstacle towards model reduction.

The issue regarding simulation of the initial system also arises when trying to quantify the efficiency of model reduction methods: ideally, the approximation errors resulting from the reduction should be computed without executing the original system.

Thus, in this paper we propose an approximation method for biochemical networks, in which no prior knowledge about the original system’s behavior is required. Our method combines the dominance concept with mass invariants of the original ODE system in order to compute inequality constraints on the species’ concentrations. These constraints are then combined with the original system of equations, in order to obtain a reduced system of ODEs that provides time-dependent lower and upper bounds on the species’ concentrations. Depending on the coarseness of detail we choose to incorporate in the mass invariant-generated inequalities, our approach can serve either as a reduction method, or to quantify the approximation errors of tropicalization reduction heuristics.

To achieve this, we abstract the original system by a box, the hyper-faces of which provide lower and upper bounds for the concentrations of the species. The two equations of the hyper-faces of a species represent simplified versions of the original differential equation of the species, in which only the dominant positive and negative monomials are considered. We refer to these equations as being tropicalized. Then, instead of interpreting the differential equations over the state of the original system, we will lift this interpretation conservatively over each hyper-face of the box. To do this, we will bound, for every hyper-face, the derivative of the corresponding coordinate in the solution of the original differential equation over the whole hyper-face. Our method should allow for formal evaluation of tropicalization approaches, and as such the bounds are derived using the dominance relations between monomials of the original ODE. Mass invariants of the original system will then be used to refine the bounds, and thus increase the accuracy of our method. By construction, the maximal solutions of the original, respectively tropicalized (i.e., abstracted) equations are related by the following soundness criterion: when both defined at time t, the state of the original system is within the hyper-box of the abstract system.

Example 2.4**.**

Let us consider the equations (2) of the Michaelis-Menten mechanism, under the QSS assumption: $k_{2}\gg k_{-1}$, i.e. $\epsilon\cdot k_{-2}\geq k_{-1}\geq 0$, for an $\epsilon\in[0,1]$. From (2.1), it follows that one can write (by extension): $k_{2}\succ_{\epsilon}k_{-1}$. Then, we can deduce the following lower and upper bounds (that we call tropicalized) on the concentration of $x_{2}$:

[TABLE]

For convenience purposes, we will use the notation $x_{1},x_{2},x_{3},x_{4}$ for the species’ concentrations, $[S],[E],[E:S],[P]$. We propose to approximate the state of the system by a box of $\mathbb{R}^{4}$. A box of $\mathbb{R}^{4}$ is a set of the form $\{(x_{1},x_{2},x_{3},x_{4})\mid\underline{x_{i}}\leq x_{i}\leq\overline{x_{i}},\forall 1\leq i\leq 4\}$, where $(\underline{x_{i}},\overline{x_{i}})$ are pairs of numbers satisfying $\underline{x_{i}}\leq\overline{x_{i}}$, $\forall 1\leq i\leq 4$. Intuitively, the real number $\underline{x_{i}}$ provides a lower bound to the value of the variable $x_{i}$, and denotes the face $\{(x_{1},x_{2},x_{3},x_{4})\in\mathbb{R}^{4}\mid x_{i}=\underline{x_{i}},\underline{x_{j}}\leq x_{j}\leq\overline{x_{j}},\forall 1\leq j\leq 4,i\neq j\}$. We will denote this face as $\mathcal{F}_{\underline{x_{i}}}(\underline{x_{1}},\overline{x_{1}},\underline{x_{2}},\overline{x_{2}},\underline{x_{3}},\overline{x_{3}},\underline{x_{4}},\overline{x_{4}})$. The other faces are defined in the same way, and the same reasoning applies to $\overline{x_{i}}$, which provides an upper bound to the same variable. For ease of notation, we shall use $\underline{\overline{x}}$ to denote the vector $(\underline{x_{1}},\overline{x_{1}},\underline{x_{2}},\overline{x_{2}},\underline{x_{3}},\overline{x_{3}},\underline{x_{4}},\overline{x_{4}})$ .

Next, let us consider the following functions:

[TABLE]

The abstraction of the concrete system of equations is then defined as $\begin{cases}\frac{d\underline{x_{i}}}{dt}=F^{\#}_{\underline{x_{i}}}(\underline{\overline{x}})\\ \frac{d\overline{x_{i}}}{dt}=F^{\#}_{\overline{x_{i}}}(\underline{\overline{x}})\end{cases}$,

$\forall 1\leq i\leq 4$.

If we fix the same initial conditions for both the concrete and the abstracted system, $\underline{x_{i}}(0)=\overline{x_{i}}(0)=x_{i}(0),\forall 1\leq i\leq 4$, we can relate the solution of the abstract system to that of the original one. For every $1\leq i\leq 4$, the real number $F^{\#}_{\underline{x_{i}}}(\underline{x_{1}},\overline{x_{1}},\underline{x_{2}},\overline{x_{2}},\underline{x_{3}},\overline{x_{3}},\underline{x_{4}},\overline{x_{4}})$ provides a lower bound to the value of the function $\frac{dx_{i}}{dt}$ over the face $\mathcal{F}_{\underline{x_{i}}}$, whereas the real number $F^{\#}_{\overline{x_{i}}}(\underline{x_{1}},\overline{x_{1}},\underline{x_{2}},\overline{x_{2}},\underline{x_{3}},\overline{x_{3}},\underline{x_{4}},\overline{x_{4}})$ provides an upper bound to the value of the function $\frac{dx_{i}}{dt}$ over the face $\mathcal{F}_{\overline{x_{i}}}$. That is to say, we have $\frac{d\underline{x_{i}}}{dt}\leq\frac{dx_{i}}{dt}$, for every pair $(x_{1},x_{2},x_{3},x_{4})\in\mathcal{F}_{\underline{x_{i}}}$, and $\frac{dx_{i}}{dt}\leq\frac{d\overline{x_{i}}}{dt}$, for every pair $(x_{1},x_{2},x_{3},x_{4})\in\mathcal{F}_{\overline{x_{i}}}$. Then, using the results of [19], we can conclude that, for every time point $t$, and $\forall 1\leq i\leq 4$, the bounds:

[TABLE]

are satisfied. Thus, the solution of the abstract system of equations provides lower and upper bounds for the value of the variables of the original system of equations.

Remark 2.5**.**

In the above example, in order to obtain safe lower/upper bounds on $x_{i}$’s concentration, we make the variables range over the hyper-faces. One notices that the variable $x_{i}$ is treated specifically in the derivatives of the variables $\underline{x_{i}},\overline{x_{i}}$ - any of its occurences is replaced by the variable corresponding to the hyper-face we want to bound. By contrast, the other variables, $x_{j}$, are replaced according to the sign of their occurence:

• •

in $\frac{d\underline{x_{i}}}{dt}$, $x_{j}$ is replaced with $\begin{cases}\overline{x_{j}},\text{ if }x_{j}\text{ occurs negatively, }\\ \underline{x_{j}},\text{ if }x_{j}\text{ occurs positively. }\end{cases}$

• •

in $\frac{d\overline{x_{i}}}{dt}$, $x_{j}$ is replaced with $\begin{cases}\overline{x_{j}},\text{ if }x_{j}\text{ occurs positively, }\\ \underline{x_{j}},\text{ if }x_{j}\text{ occurs negatively. }\end{cases}$

This comes from the fact that the derivative on $x_{i}$ is evaluated on the corresponding hyper-face, which allows for greatly reducing the loss of precision. For a formal proof of the soundness of this approach, the reader is referred to [19]. Intuitively, it is justified by the intermediate value theorem: given a family of functions $\{f_{i}\}$ over the real field, if one function $f_{i}$ does not take the highest value at time $t$, whereas it is the case at time $t^{\prime\prime}>t$, then necessarily, there exists a time $t^{\prime}$ such that $t<t^{\prime}\leq t^{\prime\prime}$ in which $f_{i}$ takes the highest value while crossing another function of the family.

In Example 2.4, the inequality constraints on the concentrations of species were determined based on reaction rates constants that verify the QSS condition. In Fig.1, we show the time-evolution of the bounds on the concentration of the 4 species in the Michaelis-Menten system, for an arbitrarily chosen set of reaction rate constants and initial concentrations that satisfy the QSS condition (i.e., $k_{2}\succ_{\epsilon}k_{-1}$, and $[S]\succ_{\epsilon}[E]+[E:S]$ at time $t=0$). Nonetheless, our model reduction is sound no matter the value of initial concentrations and reaction rate constants. The equations have been integrated using the solver $ode15s$ of Matlab[22]. Strictly speaking, numerical errors stemming from numerical integration may accumulate throughout the simulation, but herein we choose to ignore them.

In Fig.1, we notice that the bounds diverge at a fast rate from the original trajectory, despite the restriction of the derivative’s evaluation on the hyper-face of the box (as explained in Note 1). A way to improve accuracy is to take into account the original system’s mass invariants, when computing the bounds. In general, a biochemical system can have several conservation laws/mass invariants, which are linear functions $b_{1}(\textbf{x}),\ldots,b_{m}(\textbf{x})$ of the concentrations, that are constant in time. These equality constraints can be used to refine the bounds on the initial system’s species’ concentrations. We can safely further restrict the evaluation of the derivative of each coordinate to the intersection of the corresponding hyperface with the subspace delimited by the conservation laws containing the variable itself. Because a variable can appear in more than one mass invariant, we choose to keep the most optimistic bound that can be computed by intersecting the hyper-face with the mass invariant subspace: the greatest lower bound, respectively the smallest upper bound.

Example 2.6**.**

In the Michaelis-Menten system, the total number of enzymes is constant, and so is the overall number of substrates and product. The two conservation laws can be written as: $\begin{cases}x_{2}(t)+x_{3}(t)=e_{0}\\ x_{1}(t)+x_{3}(t)+x_{4}(t)=s_{0}\end{cases}$, with $e_{0}=x_{2}(0)+x_{3}(0)$, and $s_{0}=x_{1}(0)+x_{3}(0)+x_{4}(0)$.

Assuming once more that $k_{2}\gg k_{-1}$, by substituting $x_{3}$ by $e_{0}-x_{2}$ or $s_{0}-x_{1}-x_{4}$ into 12, three equivalent tropicalized upper bounds on the concentration of $x_{2}$ are obtained:

[TABLE]

Lifting the interpretation of the differential equations over the hyper-face corresponding to $\overline{x_{2}}$ results in three different expressions for the upper bound on $\frac{dx_{2}}{dt}$, of possibly different accuracies:

[TABLE]

The most accurate sound upper bound on $\frac{dx_{2}}{dt}$ then writes as:

[TABLE]

Note 2.7**.**

The choice to introduce $\min$ and $\max$ operations in the expressions of the computed bounds is accounted for by our initial motivation: because existing tropicalization reduction heuristics are not justified by rigourous estimates, we aim to provide a method for quantifying errors stemming from such tropicalization reduction approaches, at the same time creating a tropicalization approach with guarantees222We nonetheless stress that our goal is not to correct the faults of existing tropicalization-inspired reduction methods, but rather quantify them by proposing a more rigorous tropicalization approach, in which the dominated monomials are bounded, rather than discarded from the ODEs. As such, we aim at computing error bounds that are as precise as possible, hence the choice of using $\min$ and $\max$ operations for bound refinement, albeit with the disadvantage of using functions that are not $\mathbb{C}^{1}$, thus introducing non-smooth vector fields. The trade-off between smoothness and precision can be tuned according to the desired goal: less precise bounds can be obtained by choosing to use smooth functions. Moreover, smoothness of vector fields is generally not guaranteed during the numerical simulation of biochemical models: as the model variables represent biochemical species’ concentrations, a good practice is to call the numerical solvers used to approximate the system’s behavior using with the ’Non-Negative’ option, which amounts to introducing a $\max$ operation into the equations (i.e., $\max(0,x_{i})$), in order to prevent negative values of variables.

The same reasoning can be applied to all variables appearing in the expression of $\frac{dx_{2}}{dt}$, in order to obtain the most accurate upper bound:

[TABLE]

Note 2.8**.**

In (15), when computing the third bound, instead of substituting $x_{3}$ by its conservation law expression, $s_{0}-x_{1}-x_{4}$, we choose to bound its value by an expression not containing $x_{4}$. We do so in order to avoid introducing supplementary variables w.r.t. those present in the tropicalized original bound (i.e., $x_{1}$, $x_{2}$ and $x_{3}$). This method, in which mass invariant partial refinement is introduced after the tropicalized bounds have been computed, can be considered as a per se model reduction method, as no supplementary information/complexity is introduced by incorporating the conservation laws. By contrast, the approach in which all the information contained by the conservation laws is exploited in order to derive the most accurate bounds can constitute a method of error-estimation for tropicalization based reduction heuristics. We present the two different methods formally in Section 3.

The issue of specifying QSS species and QE reactions a priori, when performing model reductions, is circumvented by our method. Instead, the notion of region is used in order to eliminate monomials from the species’ ODEs. Our method uses static inspection of each ODE, in order to partition the state space into different regions according to which production, respectively consumption terms dominate the others. Using this partitioning, simplified expressions bounding the species concentrations are derived for each region of the state space, allowing symbolic simplification and limiting numerical approximations.

Example 2.9**.**

In the case of the Michaelis-Menten mechanism, there are three possible dominance regions for $\frac{dx_{2}}{dt}$ leading to three possible pairs of lower and upper bounds:

1. \normalshape(1)

Region 1, if $k_{-1}$ dominates $k_{2}$:

$\epsilon k_{-1}\geq k_{2}\Rightarrow k_{-1}x_{3}-k_{1}x_{2}x_{1}\leq\frac{dx_{2}}{dt}\leq(1+\epsilon)k_{-1}x_{3}-k_{1}x_{2}x_{1}$** 2. \normalshape(2)

Region 2, if $k_{2}$ dominates $k_{-1}$:

$\epsilon k_{2}\geq k_{-1}\Rightarrow k_{2}x_{3}-k_{1}x_{2}x_{1}\leq\frac{dx_{2}}{dt}\leq(1+\epsilon)k_{2}x_{3}-k_{1}x_{2}x_{1}$** 3. \normalshape(3)

Region 3, if there is no dominant rate (i.e., $k_{-1}$ and $k_{2}$ are of comparable magnitude):

$\begin{cases}\epsilon k_{-1}\leq k_{2}\\ \epsilon k_{2}\leq k_{-1}\end{cases}\Rightarrow(k_{-1}+k_{2})x_{3}-k_{1}x_{2}x_{1}\leq\frac{dx_{2}}{dt}\leq(k_{-1}+k_{2})x_{3}-k_{1}x_{2}x_{1}$**

The complete system of equations obtained using mass invariants refinement of bounds, for all the possible dominance regions, can be found in Example 2.10. The improvement of bound accuracy via mass invariants can be observed in Fig. 2. As expected, one can also observe in Fig.2 that results become more precise as the value of $\epsilon$ increases, i.e. as $k_{-1}$ and $k_{2}$ become more separated.

Example 2.10**.**

For convenience purposes, denote the species concentrations, $[S],[E],[E:S],[P]$, using $x_{1},x_{2},x_{3},x_{4}$. Then, the derivatives of the lower and upper bounds of the original system’s species’ concentrations write as:

• •

$\frac{d\underline{x_{1}}}{dt}=k_{-1}\max(\underline{x_{3}},e_{0}-\overline{x_{2}})-k_{1}\min(\underline{x_{1}},s_{0}-\underline{x_{3}})\min(\overline{x_{2}},e_{0}-\underline{x_{3}})$**

• •

$\frac{d\overline{x_{1}}}{dt}=k_{-1}\min(\overline{x_{3}},e_{0}-\underline{x_{2}},s_{0}-\overline{x_{1}})-k_{1}\min(\overline{x_{1}},s_{0}-\underline{x_{3}})\max(\underline{x_{2}},e_{0}-\overline{x_{3}})$**

• •

$\frac{d\underline{x_{2}}}{dt}=\underline{c_{+}}\max(\underline{x_{3}},e_{0}-\underline{x_{2}})-k_{1}\min(\overline{x_{1}},s_{0}-\underline{x_{3}})\min(\underline{x_{2}},e_{0}-\underline{x_{3}})$, with $\underline{c_{+}}=\begin{cases}k_{-1},&\text{ if }\epsilon k_{-1}\geq k_{2}\\ k_{2},&\text{ if }\epsilon k_{2}\geq k_{-1}\\ (k_{-1}+k_{2}),&\text{ otherwise }\end{cases}$

• •

$\frac{d\overline{x_{2}}}{dt}=\overline{c_{+}}\min(\overline{x_{3}},e_{0}-\overline{x_{2}},s_{0}-\underline{x_{1}})-k_{1}\max(\underline{x_{1}},0)\max(\overline{x_{2}},e_{0}-\overline{x_{3}})$, with $\overline{c_{+}}=\begin{cases}(1+\epsilon)k_{-1},&\text{ if }\epsilon k_{-1}\geq k_{2}\\ (1+\epsilon)k_{2},&\text{ if }\epsilon k_{2}\geq k_{-1}\\ (k_{-1}+k_{2}),&\text{ otherwise }\end{cases}$

• •

$\frac{d\underline{x_{3}}}{dt}=k_{1}\max(\underline{x_{1}},0)\max(\underline{x_{2}},e_{0}-\underline{x_{3}})-\underline{c_{-}}\min(\underline{x_{3}},e_{0}-\underline{x_{2}},s_{0}-\underline{x_{1}})$, with $\underline{c_{-}}=\begin{cases}(1+\epsilon)k_{-1},&\text{ if }\epsilon k_{-1}\geq k_{2}\\ (1+\epsilon)k_{2},&\text{ if }\epsilon k_{2}\geq k_{-1}\\ (k_{-1}+k_{2}),&\text{ otherwise }\end{cases}$

• •

$\frac{d\overline{x_{3}}}{dt}=k_{1}\min(\overline{x_{1}},s_{0}-\overline{x_{3}})\min(\overline{x_{2}},e_{0}-\overline{x_{3}})-\overline{c_{-}}\max(\overline{x_{3}},e_{0}-\overline{x_{2}})$,with $\overline{c_{-}}=\begin{cases}k_{-1},&\text{ if }\epsilon k_{-1}\geq k_{2}\\ k_{2},&\text{ if }\epsilon k_{2}\geq k_{-1}\\ (k_{-1}+k_{2}),&\text{ otherwise }\end{cases}$

• •

$\frac{d\underline{x_{4}}}{dt}=k_{2}\max(\underline{x_{3}},0)$**

• •

$\frac{d\overline{x_{4}}}{dt}=k_{2}\min(\overline{x_{3}},s_{0}-\overline{x_{4}})$**

The Michaelis-Menten system represents a particular, simple case study: the choice of reaction rate constants fixes the dominance region in which the system evolves. In general, the state of a biochemical network can traverse multiple such regions, as the dominant monomials can change from one concentration domain to another. Thus, we next introduce a case study in which the dominant monomials are concentration-dependent, which in turn means that the dominance region is no longer fixed. Our method is designed with this more general situation in mind: having computed the most accurate bounds for each region of the state space partitioning, and having no information regarding the region in which the original system evolves at a given time $t$, our approach chooses the least accurate local bound, in order to ensure global soundness.

2.3 Motivating example: A DNA model

We construct a simple extension of the Michaelis-Menten system, in which the product formation reaction is catalazyed by a dimer of an enzyme $M$. The reaction system and its ODE system333once again, we denote the species $M,M_{2},DNA,M_{2}.DNA,P$ by $x_{1},x_{2},x_{3},x_{4},x_{5}$ read:

[TABLE]

[TABLE]

The mass invariants write:

[TABLE]

Dominance regions become concentration dependent: for example, the dominant positive monomial in $\frac{dx_{2}}{dt}$ is determined by the dominance relations between both the concentrations of $x_{1}$ and $x_{4}$, and between reaction rate constants $k_{1}$ and $k_{-2}$.

This DNA example will serve as a case study for the remainder of our paper.

3 Combining ODEs and mass invariants

3.1 Model reduction using conservative numerical approximations

The guarantees of our method are a consequence of a carefully designed symbolic propagation of inequality constraints on the species’ concentrations. Thus, symbolic transformations have to be applied on numerical expressions, of which we introduce a syntax and semantics. We also introduce an alternative definition of a biochemical model to that presented in Sect. 2, which is then used to define and justify our approximation method.

Definition 3.1**.**

(Syntax of expressions) Let $\mathcal{S}$ be a set of variables. We define an $\mathcal{S}$-expression inductively, as follows444the syntactic operators are written using a superscript dot, in order to distinguish them from their associated mathematical functions:

1. \normalshape(1)

each positive real number $k\in\mathbb{R}_{+}$ is an $\mathcal{S}$-expression; 2. \normalshape(2)

each variable $x\in\mathcal{S}$ is an $\mathcal{S}$-expression; 3. \normalshape(3)

if $e$ is an $\mathcal{S}$-expression, then $(\dot{-}e)$ is an $\mathcal{S}$-expression; 4. \normalshape(4)

if $e_{1}$ and $e_{2}$ are $\mathcal{S}$-expressions, then $(e_{1}\dot{+}e_{2})$, $(e_{1}\dot{\cdot}e_{2})$, $\dot{\min}(e_{1},e_{2})$, $\dot{\max}(e_{1},e_{2})$ are all $\mathcal{S}$-expressions;

The set of $\mathcal{S}$-expressions is denoted as $Expr_{\mathcal{S}}$. Given an $\mathcal{S}$-expression $e$, we define its support, denoted $supp(e)$, as the set of variables it contains.

Definition 3.2**.**

(Semantics of expressions) Let $\mathcal{S}$ be a set of variables and $e$ be an $\mathcal{S}$-expression. The semantics of the expression $e$ is the function $\llbracket e\rrbracket_{\mathcal{S}}:\mathbb{R}^{\mathcal{S}}\rightarrow\mathbb{R}$, defined inductively as follows:

1. \normalshape(1)

$\forall c\in\mathbb{R},\llbracket c\rrbracket_{\mathcal{S}(\rho)}=c$** 2. \normalshape(2)

$\forall x\in\mathcal{S},\llbracket x\rrbracket_{\mathcal{S}(\rho)}=\rho(x)$** 3. \normalshape(3)

$\forall e\in Expr_{\mathcal{S}},\llbracket\dot{-}e\rrbracket_{\mathcal{S}(\rho)}=-\llbracket e\rrbracket_{\mathcal{S}(\rho)}$** 4. \normalshape(4)

$\forall e_{1},e_{2}\in Expr_{\mathcal{S}},\llbracket e_{1}\dot{+}e_{2}\rrbracket_{\mathcal{S}(\rho)}=\llbracket e_{1}\rrbracket_{\mathcal{S}(\rho)}+\llbracket e_{2}\rrbracket_{\mathcal{S}(\rho)}$** 5. \normalshape(5)

$\forall e_{1},e_{2}\in Expr_{\mathcal{S}},\llbracket e_{1}\dot{\cdot}e_{2}\rrbracket_{\mathcal{S}(\rho)}=\llbracket e_{1}\rrbracket_{\mathcal{S}(\rho)}\llbracket e_{2}\rrbracket_{\mathcal{S}(\rho)}$** 6. \normalshape(6)

$\forall e_{1},e_{2}\in Expr_{\mathcal{S}},\llbracket\dot{\min}(e_{1},e_{2})\rrbracket_{\mathcal{S}(\rho)}=\min(\llbracket e_{1}\rrbracket_{\mathcal{S}(\rho)}\llbracket e_{2}\rrbracket_{\mathcal{S}(\rho)})$** 7. \normalshape(7)

$\forall e_{1},e_{2}\in Expr_{\mathcal{S}},\llbracket\dot{\max}(e_{1},e_{2})\rrbracket_{\mathcal{S}(\rho)}=max(\llbracket e_{1}\rrbracket_{\mathcal{S}(\rho)}\llbracket e_{2}\rrbracket_{\mathcal{S}(\rho)})$**

for every environment $\rho\in\mathbb{R}^{\mathcal{S}}$.

We use Defs. 3.1 and 3.2 to define the notion of system of symbolic differential equations and symbolic equality constraints derived from conservation laws.

Definition 3.3**.**

(Symbolic ODE system)

A system of symbolic ordinary differential equations and equality constraints modeling a biochemical network is a tuple $(\mathcal{S},\mathbb{I},\mathbb{F},(\mathbb{E}_{b}))$, where:

• •

$\mathcal{S}=\{x_{1},\ldots,x_{s}\}$* is a set of variables, denoting species’ concentrations,*

• •

$\mathbb{I}:\mathcal{S}\rightarrow\mathbb{R}^{+}$* is a non-negative function, mapping each species to its initial concentration,*

• •

$\mathbb{F}:\mathcal{S}\rightarrow Expr_{\mathcal{S}}$* is a function describing the evolution of species’ concentrations, as described in Eq.(4):*

[TABLE]

with $P_{i}^{+/-}\in Expr_{\mathcal{S}}$, Laurent polynomials with positive coefficients,

• •

$\{\mathbb{E}_{b}\}$555the number $b$ indexes the different ways of expressing a species $x_{i}$, by using the mass invariants in which it appears* is a family of functions from the set $\mathcal{S}$ into the set $Expr_{\mathcal{S}}$, denoting equality constraints derived from conservation laws, such that $\forall f:\mathcal{S}\rightarrow\mathbb{R}^{\mathbb{R}^{+}}$ satisfying*

[TABLE]

the constraint

[TABLE]

is satisfied for every function $\mathbb{E}_{b}$ of the family $\{\mathbb{E}_{b}\}$, $\forall x_{i}\in\mathcal{S}$, and for every time $t\in\mathbb{R}^{+}$.

Example 3.4**.**

(A DNA example) In our running example, $\mathcal{S}=\{x_{1},x_{2},x_{3},x_{4}\}$, $\mathbb{F}$ is defined by the equations of (19), and the equality constraints derived from the conservation laws of (20) write:

[TABLE]

We partition the state space of each ODE into regions, each one defined by the corresponding pair of dominant monomials, ($Dom(P_{i}^{+}(\textbf{x})),Dom(P_{i}^{-}(\textbf{x}))$). At any given time $t$, several monomials can be dominant, which can lead to an exponential number of possible regions. To circumvent this issue and obtain a linear number of regions, we choose to replace each region that has more than one dominant term with the unique region in which no term is dominant: if $|Dom(P_{i}^{\pm}(\textbf{x}))|>1$, we choose to keep $P_{i}^{\pm}(\textbf{x})$ in the reduced ODE, instead of replacing it with $Dom(P_{i}^{\pm}(\textbf{x}))$. The following definition formalizes these concepts.

Definition 3.5**.**

(State partitioning of a symbolic ODE) Let $(\mathcal{S},\mathbb{I},\mathbb{F},\{\mathbb{E}_{b})\}$ be a symbolic ODE system, and $\epsilon\in[0,1]$ a scale separation constant. Then, for every variable $x_{i}\in\mathcal{S}$, if $P_{i}^{+}=\sum\limits_{j=1}^{p}M_{j}^{+}$ and $P_{i}^{-}=\sum\limits_{j=1}^{n}M_{j}^{-}$, its state space can be partitioned into $(p+1)\times(n+1)$ regions, each one determined by the corresponding pair of dominant monomials

[TABLE]

Example 3.6**.**

(A DNA example) In Eq.(19), the state space of $x_{2}$ can be partitioned in 9 regions, as its ODE contains 2 positive terms and 2 negative terms:

\begin{array}[]{ccc}r_{2}^{1,1}=(k_{1}x_{1}^{2},k_{-1}x_{2});\quad r_{2}^{2,1}=(k_{-2}x_{4},k_{-1}x_{2});\quad r_{2}^{3,1}=(k_{1}x_{1}^{2}+k_{-2}x_{4},k_{-1}x_{2})\\ r_{2}^{1,2}=(k_{1}x_{1}^{2},k_{2}x_{2}x_{3});\quad r_{2}^{2,2}=(k_{-2}x_{4},k_{2}x_{2}x_{3});\quad r_{2}^{3,2}=(k_{1}x_{1}^{2}+k_{-2}x_{4},k_{2}x_{2}x_{3})\\ r_{2}^{1,3}=(k_{1}x_{1}^{2},k_{-1}x_{2}+k_{2}x_{2}x_{3});\quad r_{2}^{2,3}=(k_{-2}x_{4},k_{-1}x_{2}+k_{2}x_{2}x_{3});\quad r_{2}^{3,3}=(k_{1}x_{1}^{2}+k_{-2}x_{4},k_{-1}x_{2}+k_{2}x_{2}x_{3})\end{array}**

We next use the dominance relations that define each region, in order to obtain region-specific lower and upper bounds on the ODE being considered. The next definition formalizes this procedure:

Definition 3.7**.**

(Region-specific tropicalized bounds) Given a symbolic ODE system $(\mathcal{S},\mathbb{I},\mathbb{F},(\mathbb{E}_{b}))$, and the set of regions $r_{i}^{k,l}$ for each species $x_{i}$, the dominance definition 2.1 can be used to define the following functions, for every region $r_{i}^{k,l}$:

[TABLE]

[TABLE]

Functions $\mathbb{F}^{k,l}_{\downarrow}$ and $\mathbb{F}^{k,l}_{\uparrow}$ provide symbolic tropicalized lower, resp. upper bounds for $\mathbb{F}(x_{i})$ on region $r_{i}^{k,l}$.

Example 3.8**.**

(A DNA example) In our running example, in region $r_{2}^{2,1}=(k_{-2}x_{4},k_{-1}x_{2})$, the dominant positive (production) monomial is $k_{-2}x_{4}$, and the dominant negative (consumption) monomial is $k_{-1}x_{2}$. Formally, this writes as $\epsilon\cdot k_{-2}x_{4}\geq k_{1}x_{1}^{2}\geq 0$, and $\epsilon\cdot k_{-1}x_{2}\geq k_{2}x_{2}x_{3}\geq 0$.

Thus, the $r_{2}^{2,1}$ specific tropicalized bounds write as:

[TABLE]

which by construction satisfy $\mathbb{F}_{\downarrow}^{2,1}(x_{2})\leq\frac{dx_{2}}{dt}\leq\mathbb{F}_{\uparrow}^{2,1}(x_{2})$.

The bounds of Def. 3.7 can further be refined by using the mass invariants given by the family of functions $\{\mathbb{E}_{b}\}$, as follows:

Definition 3.9**.**

(Region-specific refined tropicalized bounds) Given a symbolic ODE system $(\mathcal{S},\mathbb{I},\mathbb{F},(\mathbb{E}_{b}))$, the set of regions $r_{i}^{k,l}$ and the symbolic tropicalized bounds $\mathbb{F}^{k,l}_{\downarrow}(x_{i})$, $\mathbb{F}^{k,l}_{\uparrow}(x_{i})$ for each species $x_{i}$, we define the following bounds:

[TABLE]

with $\mathcal{V}=supp(\mathbb{F}^{k,l}_{\downarrow}(x_{i}))=supp(\mathbb{F}^{k,l}_{\uparrow}(x_{i}))$, $\mathcal{V}_{b}=supp({E}_{b}(x_{j}))$, for each function $\mathbb{E}_{b}$ of the family ($\mathbb{E}_{b}$) that applies to the variable $x_{j}$, and $b_{i}^{k,l}(x_{j})\in Expr_{\mathcal{V}}$ is either 0, or a bound generated by the dominating monomial inequality constraints.

Example 3.10**.**

(A DNA example)

When dealing with the tropicalized bounds of Ex.(3.8), one needs to refine the bounds of the variables in their support: $\mathcal{V}=\{x_{2},x_{4}\}$. We do so by using their respective equality constraints from (21): $\mathbb{E}_{1}(x_{2}),\mathbb{E}_{2}(x_{2}),\mathbb{E}_{1}(x_{4})$, and $\mathbb{E}_{2}(x_{4})$.

What’s more, the second dominance inequality of region $r_{2}^{2,1}$ in Ex.(3.8)can be rewritten as $x_{3}\leq\epsilon\frac{k_{-1}}{k_{2}}$. This allows for a new upper bound on variable $x_{3}$: $b_{2}^{2,1}(x_{3})=\epsilon\frac{k_{-1}}{k_{2}}\in Expr_{\mathcal{V}}$.

Using Def.3.9, the $r_{2}^{2,1}$-specific bounds on $x_{2}$ and $x_{4}$ write as:

\begin{array}[]{cccc}\mathbb{L}_{2,1}^{2,1}(x_{2})=0;\quad\mathbb{L}_{2,2}^{2,1}(x_{2})=0;\quad\mathbb{L}_{2,1}^{2,1}(x_{4})=DNA_{0}-\epsilon\frac{k_{-1}}{k_{2}};\quad\mathbb{L}_{2,2}^{2,1}(x_{4})=0\\ \mathbb{U}_{2,1}^{2,1}(x_{2})=\frac{M_{0}}{2}-x_{4};\quad\mathbb{U}_{2,2}^{2,1}(x_{2})=\frac{M_{0}}{2}-DNA_{0}+\epsilon\frac{k_{-1}}{k_{2}};\quad\mathbb{U}_{2,1}^{2,1}(x_{4})=DNA_{0};\quad\mathbb{U}_{2,2}^{2,1}(x_{4})=\frac{M_{0}}{2}-x_{2}\end{array}**

Using mass invariants to compute the most optimistic bound is done inductively over the $\mathcal{S}$ expressions of the candidate bounds, by applying usual formulae of interval arithmeticsto propagate the $\dot{\min}$ and $\dot{max}$ operators. The resulting evaluation functions, which we call $f_{\dot{\min}}$ and $f_{\dot{\max}}$ respectively, are detailed in Appendix A.

With all this in place, we can proceed to the definition of the reduced system.

Definition 3.11**.**

(Reduced system) Let $\mathcal{A}=(\mathcal{S},\mathbb{I},\mathbb{F},(\mathbb{E}_{b}))$ be a system of ordinary equations with equality constraints. The reduction of the system $\mathcal{A}$ is defined as the triple $(\mathcal{S}^{\#},\mathbb{I}^{\#},\mathbb{F}^{\#})$, with:

1. \normalshape(1)

$\mathcal{S}^{\#}=\{\underline{x_{i}}\mid x_{i}\in\mathcal{S}\}\cup\{\overline{x_{i}}\mid x_{i}\in\mathcal{S}\}$** 2. \normalshape(2)

$\mathbb{I}^{\#}:\mathcal{S}^{\#}\rightarrow\mathbb{R}^{+}$* is defined by $\mathbb{I}^{\#}(\underline{x_{i}})=\mathbb{I}^{\#}(\overline{x_{i}})=\mathbb{I}(x_{i}),\forall x_{i}\in\mathcal{S}$ h* 3. \normalshape(3)

$\mathbb{F}^{\#}:\mathcal{S}^{\#}\rightarrow Expr_{\mathcal{S}^{\#}}$, defined as:

[TABLE]

for every variable $x_{i}\in\mathcal{S}^{\prime}$, where:

• •

$\rho_{1}^{\downarrow}=\begin{cases}x_{j}\mapsto\dot{\max}(x_{j},\dot{\max\limits_{b}}(\mathbb{L}_{i,b}^{k,l}(x_{j}))),&\text{ if }x_{j}\in t_{\downarrow,+}^{k,l,i}\\ x_{j}\mapsto\dot{\min}(x_{j},\dot{\min\limits_{b}}(\mathbb{U}_{i,b}^{k,l}(x_{j}))),&\text{ if }x_{j}\in t_{\downarrow,-}^{k,l,i}\end{cases}$**

• •

$\rho_{2}^{\downarrow}=\begin{cases}x_{j}\mapsto\underline{x_{j}},&\text{ if }x_{i}=x_{j}\\ x_{j}\mapsto\underline{x_{j}},&\text{ if }x_{i}\neq x_{j},\text{ for positive polarity/sign occurences of }x_{j}\\ x_{j}\mapsto\overline{x_{j}},&\text{ if }x_{i}\neq x_{j},\text{ for negative polarity/sign occurences of }x_{j}\end{cases}$**

• •

$\rho_{1}^{\uparrow}=\begin{cases}x_{j}\mapsto\dot{\min}(x_{j},\dot{\min\limits_{b}}(\mathbb{U}_{i,b}^{k,l}(x_{j}))),&\text{ if }x_{j}\in t_{\uparrow,+}^{k,l,i}\\ x_{j}\mapsto\dot{\max}(x_{j},\dot{\max\limits_{b}}(\mathbb{L}_{i,b}^{k,l}(x_{j}))),&\text{ if }x_{j}\in t_{\uparrow,-}^{k,l,i}\end{cases}$**

• •

$\rho_{2}^{\uparrow}=\begin{cases}x_{j}\mapsto\overline{x_{j}},&\text{ if }x_{i}=x_{j}\\ x_{j}\mapsto\overline{x_{j}},&\text{ if }x_{i}\neq x_{j},\text{ for positive polarity/sign occurences of }x_{j}\\ x_{j}\mapsto\underline{x_{j}},&\text{ if }x_{i}\neq x_{j},\text{ for negative polarity/sign occurences of }x_{j}\end{cases}$**

Intuitively, for each region $(k,l)$ of species $x_{i}$, the reduction method first replaces $\mathbb{F}(x_{i})$ by the pair of tropicalized lower and upper bounds, $\mathbb{F}^{k,l}_{\downarrow}(x_{i})$ and $\mathbb{F}^{k,l}_{\uparrow}(x_{i})$, that result directly from the dominance inequalities that characterize the region. Then, $\mathbb{F}^{k,l}_{\downarrow}(x_{i})$ and $\mathbb{F}^{k,l}_{\uparrow}(x_{i})$ are refined, using the bounds on variables that can be deduced from the conservation laws of the original system. For example, replacing any occurence of a variable $x_{j}$ in $\mathbb{F}^{k,l}_{\downarrow}(x_{i})$ with one of its expressions $\mathbb{E}_{b}(x_{j})$ (or with its appropriate bound derived from $\mathbb{E}_{b}(x_{j})$666$\mathbb{L}_{i,b}^{k,l}(x_{j})$ for the positive occurences of $x_{j}$, and $\mathbb{U}_{i,b}^{k,l}(x_{j})$ for its negative occurences) results in another safe upper bound for $\mathbb{F}(x_{i})$. By choosing the minimum such candidate bound, one obtains the most accurate, locally safe upper bound. The same reasoning applies to the computation of lower bounds, but the $\min$ operation is replaced with $max$.

In order to obtain safe (i.e., correct) global bounds, the least precise local bounds are chosen: the miminal lower, resp. the maximal upper bounds.

Finally, the interpretation of the variables is lifted over the hyper-faces. Any occurence of $x_{i}$ is replaced with its analogue corresponding to the hyperface we want to bound, while the others are replaced to their analogue given by the polarity of their occurence, as explained in Note 2.5.

Theorem 3.12**.**

Let $\mathcal{A}=(\mathcal{S},\mathbb{I},\mathbb{F},(\mathbb{E}_{b}))$ be a system of ordinary equations with equality constraints. Let $(\mathcal{S}^{\#},\mathbb{I}^{\#},\mathbb{F}^{\#})$ be a reduction of the system $\mathcal{A}$.

Let $f$ be a function from the set $\mathcal{S}$ into the set $\mathbb{R}^{\mathbb{R}^{+}}$ s.t. for every variable $x_{i}\in\mathcal{S}$, we have:

[TABLE]

and $f\#$ be a function from the set $\mathcal{S}^{\#}$ into the set $\mathbb{R}^{\mathbb{R}^{+}}$ s.t. for very abstract variable $x_{i}\#\in\mathcal{S}^{\#}$, we have:

[TABLE]

Under these assumptions, we have that for every variable $x_{i}\in\mathcal{S}$ and every time $t\in\mathbb{R}^{+}$:

[TABLE]

i.e.**, the reduced system provides sound lower and upper bounds for the concentration of the original system’s species.

Example 3.13**.**

We apply our method on the DNA example constructed in Sect.2.3, for different values of the scale separating constant $\epsilon$, and for arbitrarily chosen reaction rate constants $k_{1}=k_{2}=k_{3}=0.1,k_{-1}=0.01,k_{-2}=0.00001$ and initial concentrations $[M]_{0}=1,[DNA]_{0}=0.05$. We show in Fig.3 the time evolutions on the bounds on the concentration of the product $[P]$, i.e. the variable $x_{5}$. We notice once again that the results become more precise as $\epsilon$ decreases, i.e. as the monomials become more separated. As an example, the ODE of the lower bound of species $x_{2}$ can be found in Appendix B.

4 Error estimates of tropicalized systems using conservative numerical approximations

Our approach also serves as a heuristics for quantifying errors of tropicalization approaches for biochemical model reduction, provided a slight modification is applied to Def.3.9. Instead of computing bounds using only variables from the support of the tropicalized bounds, one can use the equality constraints ${\mathbb{\{}E_{b}\}}$, to refine the accuracy of bounds. The resulting model presents a trade-off: it introduces new variables w.r.t. the support of the tropicalized bounds, albeit exclusively in the form of conservation laws which are always linear functions, but gains in bound accuracy. As such, the approximation error/accuracy of a given reduction method can be assessed by checking if the reduced trajectory lies between the lower and upper bounds computed by our method.

Example 4.1**.**

It is well known that the Michaelis-Menten reduction is valid only under the QSS and QE assumptions. In Fig.4, we simulate the reduced Michaelis-Menten system (10), as well as our modified reduced system, as presented above, for a set of initial conditions that no longer satisfy the QSS assumption, i.e. the total enzyme concentration is comparable to the total substrate concentration. As expected, the reduced Michaelis-Menten system no longer represents a good approximation of the initial enzymatic system (7); this is reflected by the fact that the trajectory of the reduced model does not lie between the lower and upper bounds computed by our approach.

4.1 Tyson’s Cell Cycle Model

The tropicalization heuristics can be difficult to justify by rigourous estimates, although this is possible in some cases[26]. For example, the existence of tropical varieties - the set of points $x\in\mathbb{R}^{n}$ where at least two monomials of $P^{-/+}$ are equal- can lead to sliding modes, which in turn represent challenges in providing accuracy justifications for hybrid models obtained using tropical ideas. Sliding modes are well known phenomena in ODEs with discontinuous vector fields, in which the dynamics can follow discontinuity hyper-surfaces where the vector field is not defined; what’s more, the conditions for the existence of sliding modes are usually intricate. As noted in [27], sliding modes can have a nefarious effect on the behavior of the tropicalized system: tropical varieties (i.e. tropical curves) decompose the state space into sectors corresponding to the smooth modes of the hybrid tropicalized system, which passes from one type of smooth dynamics to another intrinsically, when the trajectory attains the tropical curve. However, if certain conditions w.r.t. the sliding modes are fulfilled, the trajectory can continue along some tropical curve instead of changing sector, which further deviates the reduced system’s trajectory from the original one (see Figure 1 in [27], for an example).

In [27], such phenomena become apparent when tropicalization is applied to the minimal cell cycle model proposed by Tyson[31], in order to obtain a reduced hybrid model. The Tyson model describes the interplay between cyclin and cyclin dependent kinase cdc2 during the progression of the cell cycle, and demonstrates the existence of three possible regimes, that can be associated to different phases in the cell life: the biochemical system can either function as an oscillator, converge to a steady state, or behave as an excitable switch. The three possible behaviours can be associated to early embryos rapid division, arrest of unfertilised eggs and growth controlled division of somatic cells, respectively.

The dynamics of this non-linear model with rational reaction rates contains 6 species and 9 reactions, and is described by the following system of polynomial differential equations:

[TABLE]

and has the conservation law $y_{1}(t)+y_{2}(t)+y_{3}(t)+y_{4}(t)=1$, where the value 1 denoting the total initial concentration of kinase cdc2 (i.e. $y_{1}(0)+y_{2}(0)+y_{3}(0)+y_{4}(0)$) was chosen by convenience. The values of the reaction rates constants are fixed as to have the model display the oscillatory behavior: $k_{1}=0.015,k_{3}=200,k_{4}=180,k^{\prime}_{4}=0.018,k_{6}=1,k_{8}=10^{3},k_{9}=10^{6}$.

In [26, 27], a hybrid model of the Tyson cell cycle is obtained by detecting and eliminating QSS species of the original model, pruning dominated monomials, and then ultimately tropicalizing the reduced-size model. Besides having the inconvenient of analyzing trajectories of the original model in order to detect QSS species, the reduced model suffers from the the sliding mode-related issues mentioned above: although both the smooth (original) and the reduced system exhibit oscillating behavior and have stable periodic trajectory (i.e. limit cycle), the period of the tropicalized limit cycle is different with respect to that of the smooth cycle, due to the tropicalized trajectory sliding along the tropical manifolds instead of changing sectors. Having different oscillation periods means in turn that assessing the accuracy of the tropicalized reduced model is not a trivial question, as the distance between original and tropicalized trajectories is variable from cycle to cycle (as can be seen in Figure5). What’s more, it can also provide an indication of the poor performance of tropicalization based reduction methods when dealing with more complex systems, such oscillating systems.

Indeed, by applying our method to the original Tyson model, we are able to effectively provide guarantees on the reduced model, albeit not very strong ones: this could be interpreted as an indication of the poor accuracy of the tropicalized Tyson model. In Figure 5, we plot the bounds for the concentration of species $y_{4}$ obtained using our method, in order to compare the trajectory of the original model to the one of the hybrid one that can be found in [26]. The lack of oscillating behavior in the computed bounds could intuitively be explained by the difference in period of the original and reduced systems, that causes a shift at every cycle in the tropicalized trajectory w.r.t. the original behavior. Nonetheless, the obtained bounds accurately capture the amplitude of the tropicalized model. One also notes that the time points where the upper bound, respectively the tropicalized system, begin to diverge w.r.t. the original trajectory coincide.

From a more practical point of view, we note that while simulation of the tropicalized Tyson model proposed in [27] was performed in 354.876706 seconds, on a 2.2 GHz Intel Core i7 processor and 16 GB 1600 MHz DDR3 memory, simulation of the model obtained via our method was performed in only 9.775511 seconds using the same numerical integration method (i.e. Matlab’s ode15s), thus providing a significant improvement in computation time.

We note that an alternative reduced model is obtained in [26], using tropical equilibration, that circumvents the need to simulate the original system. We plan to include tropical equilibration techniques in future work.

5 Comparison with existing methods

We mentioned previously that numerical errors stemming from numerical integration are ignored herein. Indeed, numerical integration methods, while heavily used, only provide approximations of the solution of the initial value problem (IVP) of ODE systems. Even when using variable-step size methods, there are no guarantees that the approximate solution computed by the chosen method is close to the actual solution. In order to solve the drawbacks associated to traditional ODE solvers/numerical solutions of IVP, interval numerical methods for IVP are used for computing validated enclosures of the solution of an IVP for an ODE. For example, the VNODE-LP[24] C++ solver proves that a unique solution to a problem exists, and then computes rigorous bounds that are guaranteed to contain it. Such bounds can then be used to help prove theoretical results, check if a solution satisfied a condition in a safety-critical calculation, or simply to verify the results produced by a traditional ODE solver. Another example of such software is the CAPD library [6]. Both represent well-established software for computing enclosures of generic ODE systems, and are integrated in various SMT solvers (e.g., iSAT[12], dReal[15]). For a more comprehensive state of the art on such methods, the reader is refered to [25].

Interval methods for IVPs for ordinary differential equations are typically based on Taylor series expansions, which require the computation of Taylor coefficients up to some order $k$. Given a final time point, the aim is to compute interval vectors that are guaranteed to contain the solution to a given IVP, at all intermediary points. In order to compute such interval vectors, interval propagation methods are used to enclose roundoff and truncation errors in the computed bounds, and thus obtain rigorous bounds on the true solution of the ODE.

In our approach, instead of interpreting the differential equations over the state of the system, the interpretation is lifted conservatively over each hyper-face of the hyper-box abstracting the system state (i.e., we over-approximate the derivatives only on the hyper-faces). When compared to our method, interval propagation methods over-approximate the partial derivatives of the function over the whole enclosing hyper-box, instead of doing so only on the hyper-faces. This in turn means that our approach computes tighter bounds than those computed by interval methods for IVPs.

We demonstrate our claim with the following example:

Example 5.1**.**

Let us consider the following initial value problem :

$\begin{cases}\frac{dx}{dt}=y\cdot(2-cos(y))-x\cdot(2-sin(y))\\ \frac{dy}{dt}=x\cdot(2−cos(y))−y\cdot(2−sin(y))\\ x(0)=y(0)=1\end{cases}$**

As presented in Section 3, our framework can be decomposed in two independents parts: the first part consists in synthesizing bounds on the derivatives of the original system, and the second part deals with the propagation of said bounds, in order to obtain the enclosing system.

As our goal is to better understand and evaluate tropicalization approaches for biochemical model reduction, so far we chose to focus on bounds obtained by using dominance relations between monomials. The second part of our method simply represents an improved alternative to existing ODE enclosure methods, as explained above, and as such can be used in such methods in order to get better enclosure results.

For example, in order to compare the performance of our method to that of VNODE-LP and CAPD, instead of using dominance relations to derive inequality constraints on species’ concentrations, we now use the Taylor Series expansion with $k$ terms ($k$ will serve as a parameter) for the functions $sin$ and $cos$, in order to derive bounds on $\frac{dx}{dt}$ and $\frac{dy}{dt}$:

$\begin{cases}sin(x)\approx x-\frac{x^{3}}{3!}+\frac{x^{5}}{5!}-\frac{x^{7}}{7!}+\frac{x^{9}}{9!}-\ldots=\sum\limits_{n=0}^{k-1}(-1)^{n}\frac{x^{2n+1}}{(2n+1)!}\\ cos(x)\approx 1-\frac{x^{2}}{2!}+\frac{x^{4}}{4!}-\frac{x^{6}}{6!}+\frac{x^{8}}{8!}-\ldots=\sum\limits_{n=0}^{k-1}(-1)^{n}\frac{x^{2n}}{(2n)!}\end{cases}$**

*Then, for a fixed order $k$ and an $\epsilon$, instead of using dominance-related inequalities with our method, one can use the following inequalities:

[TABLE]

,

where $\epsilon=(-1)^{2k+1}\frac{x^{2k+1}}{(2k+1)!}cos(c_{x})$ for $c_{x}\in[0,x]$ is the residual for the Taylor expansion, and can be bound by $\epsilon\leq\frac{|x|^{2k+1}}{(2k+1)!}$, which in turn is $\leq(2k+1)!^{-1}$ for $x\in[-1,1]$ .

Our method then proceeds as usual to the computation of ODEs for bounds on the concentrations of $x$ and $y$.

In Fig.6, we compare the accuracy of our method to that of VNODE-LP and CAPD, for different values of the order $k$. The accuracy is given by the tightness of bounds, which can be evaluated by computing the difference between the upper and lower bounds, during a simulation. The results indicate that, when compared to existing enclosure interval methods, our approach represents a consistent improvement of several orders of magnitude, across different values of $k$.

6 Conclusion and outlook

In this paper, we present an approximation method for biochemical networks, which can also serve as a technique for evaluating the faithfulness of existing tropicalization reduction methods that do not involve guarantees. Our approach relies on the multiscaleness property of biochemical systems. Tropical geometry offers a natural framework for studying such networks. Tropical approaches [28, 29] can guide model reduction of ODE systems, by using time- and concentration scales separation to identify and neglect equation terms whose values are significantly smaller than those of other terms of the same equation. This leads to partitioning the state space into different regions, according to which term dominates the others. A similar approach is employed in our method, but instead of neglecting the dominated terms, we propose to conservatively bound their value using an amortizing scale separation constant and the value of the dominant terms. These bounds can be further refined by incorporating the conservation laws of the initial system. The resulting approximated model is composed of two-term ODEs (which we call tropicalized), which by construction provide time-dependent lower and upper bounds for the concentration of the initial system’s species. As such, our approach can also serve to test the accuracy of other given reduction methods, while circumventing the execution of the original system: the suitability of a reduction will be confirmed if the reduced model’s trajectory lies between the bounds provided by our abstraction.

We have tested our approach on the classical Michaelis-Menten system, a simple extension of it, and Tyson’s cell cycle model. Our method can be easily automatized, either using static analysis, or existing symbolic math tools777for example, Matlab’s Symbolic Math Toolbox ; as such, Definitions 3.1-3.11 are written in an operational-semantics style, as to describe the different procedures composing the algorithm that implements our method. A tool that automatizes our approach is currently being developed. Further work also includes expanding the case studies to larger networks, possibly with no conservation laws.

Appendix A Symbolic propagation of min and max

The propagation of the $\min$ and $max$ operations over $\mathcal{S}$-expression is done using the evaluation functions $f_{\dot{\min}}$ and $f_{\dot{\max}}$, which are defined by mutual induction over the syntax of $\mathcal{S}$-expressions denoting monomials:

1. \normalshape(1)

$\forall e_{1},e_{2},\ldots,e_{k}\in Expr_{\mathcal{S}}$,

• •

$f_{\dot{\min}}(e_{1},e_{2},\ldots,e_{k})\triangleq\dot{\min}(e_{1},e_{2},\ldots,e_{k})$

• •

$f_{\dot{\max}}(e_{1},e_{2},\ldots,e_{k})\triangleq\dot{\max}(e_{1},e_{2},\ldots,e_{k})$ 2. \normalshape(2)

$\forall e_{1},e_{2}\in Expr_{\mathcal{S}}$,

• •

$f_{\dot{\min}}(\dot{-}e_{1},\dot{-}e_{2})\triangleq\dot{-}(f_{\dot{\max}}(e_{1},e_{2}))$

• •

$f_{\dot{\max}}(\dot{-}e_{1},\dot{-}e_{2})\triangleq\dot{-}(f_{\dot{\min}}(e_{1},e_{2}))$ 3. \normalshape(3)

$\forall e_{1},e_{2}\in Expr_{\mathcal{S}},\forall c\in\mathbb{R}$,

• •

$f_{\dot{\min}}(c\dot{\cdot}e_{1},c\dot{\cdot}e_{2})\triangleq\begin{cases}c\dot{\cdot}f_{\dot{\min}}(e_{1},e_{2}),\text{ if }c\geq 0\\ c\dot{\cdot}f_{\dot{\max}}(e_{1},e_{2}),\text{ if }c<0\end{cases}$

• •

$f_{\dot{\max}}(c\dot{\cdot}e_{1},c\dot{\cdot}e_{2})\triangleq\begin{cases}c\dot{\cdot}f_{\dot{\max}}(e_{1},e_{2}),\text{ if }c\geq 0\\ c\dot{\cdot}f_{\dot{\min}}(e_{1},e_{2}),\text{ if }c<0\end{cases}$ 4. \normalshape(4)

$\forall e,e_{1},e_{2}\in Expr_{\mathcal{S}}$,

• •

$f_{\dot{\min}}(e_{1}\dot{\pm}e,e_{2}\dot{\pm}e)\triangleq f_{\dot{\min}}(e_{1},e_{2})\dot{\pm}e$

• •

$f_{\dot{\max}}(e_{1}\dot{\pm}e,e_{2}\dot{\pm}e)\triangleq f_{\dot{\max}}(e_{1},e_{2})\dot{\pm}e$ 5. \normalshape(5)

$\forall e,e_{1},e_{2}\in Expr_{\mathcal{S}}$,

• •

$f_{\dot{\min}}(e\dot{-}e_{1},e\dot{-}e_{2})=e\dot{-}f_{\dot{\max}}(e_{1},e_{2})$

• •

$f_{\dot{\max}}(e\dot{-}e_{1},e\dot{-}e_{2})=e\dot{-}f_{\dot{\min}}(e_{1},e_{2})$

Appendix B A DNA model: bound equations

Below, we give the equation of the lower bound on species $x_{2}$ from our running DNA model example. According to Def.3.11, the derivative of lower bound on the concentration of $x_{2}$ is computed by selecting the minimum region-dependent (i.e., local) lower bound, out of the 9 possible cases:

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]