Algorithmic Reduction of Biological Networks With Multiple Time Scales
Niclas Kruff, Christoph Lüders, Ovidiu Radulescu, Thomas Sturm, Sebastian Walcher
AAlgorithmic Reduction of Biological Networks
With Multiple Time Scales
Niclas Kruff, RWTH Aachen University, Germany [email protected]
Christoph Lüders, University of Bonn, Germany [email protected]
Ovidiu Radulescu, University of Montpellier, and CNRS UMR5235 LPHI, France [email protected]
Thomas Sturm, CNRS, Inria, and the University of Lorraine, FranceMPI Informatics and Saarland University, Germany [email protected]
Sebastian Walcher, RWTH Aachen University, Germany [email protected]
October 2020
Abstract
We present a symbolic algorithmic approach that allows to compute invariant manifolds and cor-responding reduced systems for differential equations modeling biological networks which comprisechemical reaction networks for cellular biochemistry, and compartmental models for pharmacology,epidemiology and ecology. Multiple time scales of a given network are obtained by scaling, basedon tropical geometry. Our reduction is mathematically justified within a singular perturbationsetting using a recent result by Cardin and Teixeira. The existence of invariant manifolds is sub-ject to hyperbolicity conditions, which we test algorithmically using Hurwitz criteria. We finallyobtain a sequence of nested invariant manifolds and respective reduced systems on those manifolds.Our theoretical results are generally accompanied by rigorous algorithmic descriptions suitable fordirect implementation based on existing off-the-shelf software systems, specifically symbolic com-putation libraries and Satisfiability Modulo Theories solvers. We present computational examplestaken from the well-known BioModels database using our own prototypical implementations.
1. Introduction
Biological network models describing elements in interaction are used in many areas of biology andmedicine. Chemical reaction networks are used as models of cellular biochemistry, including gene regula-tory networks, metabolic networks and signaling networks. In epidemiology and ecology, compartmentalmodels can be described as networks of interactions between compartments. Both in chemical reactionnetworks and in compartmental models the probability that two elements interact is assumed propor-tional to their abundances. This property, called mass action law in biochemistry, leads to polynomialdifferential equations in the kinetics.For differential equations that describe the development of such networks over time a crucial questionis concerned with reduction of dimension. We illustrate such a reduction and the steps involved for1 a r X i v : . [ q - b i o . M N ] O c t he classical Michaelis–Menten system, an archetype of enzymatic reactions. The differential equationsfor the concentrations of relevant chemical species, which are substrate and enzyme-substrate complex,have the form ˙ y = − εk y + ( k y + k − ) y ˙ y = εk y − ( k y + k − + k ) y , involving a small parameter ε that represents the ratio of the total concentration of the enzyme to theconcentration of the substrate. The fact that this ratio is small is an assumption of the model that hasto be verified in applications. In a first step toward reduction, a scaling transformation y = x and y = εx yields ˙ x = ε ( − k x + ( k x + k − ) x )˙ x = k x − ( k x + k − + k ) x . In a second step, one uses singular perturbation theory to obtain the famous Michaelis–Mentenequation. It consists of two components: First, we obtain a one dimensional invariant manifold givenapproximately by the quasi-steady state condition k x − ( k x + k − + k ) x = 0. This considers thefast variable x to be at the steady state and lowers dimension from two to one. Second, we obtain areduced system for the slow variable: ˙ x = − ε k k x k x + k − + k . With our example, we paraphrased the approach in a seminal paper by Heineken et al. [28], whichwas the first one to rigorously discuss quasi-steady state from the perspective of singular perturbationtheory. Realistic network models may have many species and differential equations. Considerable efforthas been put into model order reduction, i.e., finding approximate models with a smaller number ofspecies and equations, where the reduced model can be more easily analyzed than the full model [46].The scaling of parameters and variables by a small parameter ε and the study of the limit ε → ε → ε ∗ , provided some technical conditions are satisfied. To determine scalingsof polynomial or rational vector fields that model biological networks, tropical equilibration methodswere introduced and developed in a series of papers by Noel et al. [44], Radulescu et al. [47], Samal etal. [51, 50], and others. These methods open a feasible path for biological networks of high dimension.For a given system they provide a list of possible slow-fast systems, which may or may not yield invariantmanifolds and reduced equations. Other methods due to Goeke et al. [26], and recently extended tomultiple time scales by Kruff and Walcher [32], determine critical parameter values and manifolds forsingular perturbation reductions.The principal purpose of the present paper is to complement scaling with an algorithmic test forthe existence of invariant manifolds and the computation of those manifolds along with correspondingreduced systems of differential equations. In the asymptotic limit, methods from singular perturbationtheory, principally developed by Tikhonov [58] and Fenichel [21], are available. A recent extensionto multiscale systems by Cardin and Teixeira [10] turns out to be a valuable tool for the systematiccomputation of reductions with nested invariant manifolds, and allows an algorithmic approach.In the language of systems biology the situation at a given time scale can be described as follows:Faster variables have relaxed and satisfy quasi-steady state conditions, a subset of variables evolvestoward quasi-steady state values, and all slower variables are constant. The sets of quasi-steady stateconditions for relaxed variables define invariant manifolds, more precisely, they provide the lowest orderapproximations to the invariant manifolds. As the set of relaxed variables and thus quasi-steady stateconditions increases, the respective invariant manifolds get nested so that later manifolds are contained2n earlier ones. Local linear approximations of these manifolds were proposed by Valorani and Paolucci[59] using numerical methods based on the local Jacobian. However, to the best of our knowledge,constructive approaches providing the nonlinear description of these manifolds and reduced models arestill missing.From a computer science point of view, we propose a novel symbolic computation-based algorithmicworkflow for the reduction process outlined above. This includes in particular the automatic verificationof certain hyperbolicity conditions required for the validity of the reductions. We restrict ourselves tothe case of polynomial differential equations that covers mass action chemical reaction networks andcompartmental models. We present a series of algorithms that takes as input a system of polynomialautonomous ordinary differential equations together with numerical information related to the desiredcoarse graining of the scaling. As output one finally obtains a collection of nested invariant manifoldsfor the input system, associated with smaller dimensional systems that govern the dynamics on thosemanifolds. This output establishes the reduced systems discussed above.The computationally hard parts of our methods are reduced to decision problems in interpreted first-order logic over various theories. It turns out that quantifier alternation can be entirely avoided, sothat the Satisfiability Modulo Theories (SMT) framework by Nieuwenhuis et al. [41] can be applied.Several corresponding SMT solvers are freely available and professionally supported [1, 12, 15, 17].It is remarkable that we arrive with our comprehensive algorithmic work here at SMT sub-problemsfor several different logics, viz. linear integer arithmetic, linear real arithmetic, and non-linear realarithmetic. The algorithms presented here are suitable for straightforward implementation providedthat a symbolic computation library, or computer algebra system, and an SMT solver are available. Toensure this, we have realized two independent prototypical realizations in software on our own, one inPython using freely available libraries, and one in Maple.The plan of the paper is as follows: In Sect. 2.1 we introduce an abstract scaling procedure, whichassumes, for given 0 < ε ∗ <
1, the existence of families of exponents c k,J and d k for scaling polynomialcoefficients and variables, respectively. From the scaled system, higher order terms are truncated,and the obtained system is partitioned into several time scales, ordered from fastest to slowest. Acorresponding generic algorithm uses black-box functions c and d . In Sect. 2.2 we make precise onepossible way to realize c and d , based on tropical geometry. So far, our transformations are mostlyof formal nature. On these grounds, we algorithmically determine in Sect. 3 invariant manifolds andcorresponding reduced systems, which makes the formal scaling meaningful in a mathematically preciseway. In general, this is possible only for a certain number ‘ of time scales, where ‘ is explicitlyfound and—in contrast to existing alternative approaches—often larger than 2. Technically, we applyrecent results by Cardin and Teixeira [10] based on Fenichel theory. In Sect. 4, we employ symboliccomputation techniques, specifically Gröbner basis theory, to equivalently simplify our reduced systems,which are still scaled in terms of ε ∗ , c , and d . In Sect. 5, we finally transform back to the principalscale of the original system while preserving the obtained multiple time scales and the structure of thecorresponding reduced systems. In particular, the various time scale factors remain explicit. Until here,the mathematical development of our framework has been accompanied by nine algorithms, and we givea tenth top-level algorithm and make precise how various modules are combined and interact with oneanother. In Sect. 6 we discuss various computational examples with software developed in the courseof the present work. We consider models from the BioModels database, a repository of mathematicalmodels of biological processes [34]. Our primary objective is to support the understanding of ouralgorithms, which naturally comes with a high ratio of negative examples, which do not have meaningfulreductions. This is counter-balanced by a collection of biologically interesting positive examples in theAppendix A. In Sect. 7, we wrap up and point at possible future research directions.
2. Scaling of Polynomial Vector Fields
In what follows, we adopt a rather general scaling formalism that has been used recently in [43, 44, 46,47, 51] and is generally rather recurrent in the literature on singular perturbations, see for instance [42,Sect. 3]. We use the convention that the natural numbers N include 0.3 .1. An Abstract Scaling Procedure Our starting point is a parameter dependent system S of polynomial differential equations˙ y k := d y k d t = X J γ k,J y J , ≤ k ≤ n, (1)where the summation ranges over multi-indices J = ( j , . . . , j n ) ∈ N n , γ k,J ∈ R , and only finitely many γ k,J are non-zero. We abbreviate y J = y j · · · y j n n , as usual. In terms of network models, y k representsthe concentration of either a chemical species or a type of individual in a compartment. Note thatwe use positive integers as indices, instead of concrete names for species and compartments. The realcoefficients γ k,J describe actions of other species or individuals on the species or individual k . If theseactions are activations one has γ k,J >
0, whereas for repressions one has γ k,J <
0. Several species mayinteract to produce an action on a given species k . This information is contained in the number of non-zero components of J . More precisely, the order of the action, defined as the number of species neededto produce that action, is the finite cardinality of the set { i ∈ { , . . . , n } | j i = 0 } . This terminology isinspired from chemical reactions, where the order represents essentially the number of reactant species.Throughout this paper, we require that positive y k remain positive as time progresses. In otherwords, the positive first orthant U = (0 , ∞ ) n ⊆ R n is positively invariant for system (1), which isthe case, e.g., in chemical reaction networks when γ k,J y J ≥ { ( y , . . . , y n ) ∈ R n | y k = 0 } with U .We fix some small ε ∗ ∈ (0 , γ k,J = ε c k,J ∗ ¯ γ k,J , (2)with rational numbers c k,J . The tacit understanding is that only nonzero γ k,J are being considered.The intuitive idea is that the ¯ γ k,J are close to one. Moreover, we introduce a positive parameter ε andconsider the system ˙ y k = X J ε c k,J ¯ γ k,J y J , ≤ k ≤ n (3)with ε -dependent coefficients. Notice that (3) matches (1) at ε = ε ∗ . By renormalizing y k = ε d k x k , d k ∈ Q , one obtains a system in scaled variables˙ x k = X J ε c k,J + h D,J i− d k ¯ γ k,J x J , ≤ k ≤ n, (4)with D = ( d , . . . , d n ) and the dot product in R n denoted by h· , ·i . This transformation preserves thepositive invariance of U . The scaling comes with the implicit assumption that for i , j ∈ { , . . . , n } , therelative order of y i with respect to y j is bounded by y i /y j = Θ( ε d i − d j ) for ε →
0, so that all x k get thesame order of magnitude. Continuing, we set ν k = min { c k,J + h D, J i − d k | ¯ γ k,J = 0 } to obtain˙ x k = ε ν k X J ε c k,J + h D,J i− d k − ν k ¯ γ k,J x J , ≤ k ≤ n, (5)where now all exponents of ε inside the sums are nonnegative. Finally one may perform a preliminarytime scaling τ = ε µ t , µ = min { ν , . . . , ν n } to arrive at x k := d x k d τ = ε ν k − µ X J ε c k,J + h D,J i− d k − ν k ¯ γ k,J x J , ≤ k ≤ n, (6)with all exponents nonnegative. We are interested in system (6) for variable ε >
0, in the asymptoticlimit ε → ν i − µ in vectors z , . . . , z m , where z k ∈ R n k for k ∈ { , . . . , m } , in ascending order of exponents and such that n + . . . + n m = n . We obtain a4ystem of the form z k = ε a k e f k ( z, ε ) = ε a k (cid:16) e f k ( z,
0) + ε a k, p k, + · · · + ε a k,wk p k,w k (cid:17) = ε a k (cid:16) e f k ( z,
0) + o (1) (cid:17) , ≤ k ≤ m, (7)where a k , a k,j ∈ Q , 0 = a < a < . . . < a m , 0 < a k,j , and p k,j are multivariate polynomials in z for1 ≤ k ≤ m and 2 ≤ j ≤ w k . Note that the case m = 1 is not excluded. By substituting δ := ε /q , witha sufficiently large positive integer q , one ensures that only nonnegative integer powers of δ appear: z k = δ b k b f k ( z, δ ) = δ b k (cid:16) b f k ( z,
0) + δ b k, p k, + · · · + δ b k,wk p k,w k (cid:17) = δ b k (cid:16) b f k ( z,
0) + o (1) (cid:17) , ≤ k ≤ m, (8)where b k , b k,j ∈ N , 0 = b < b < . . . < b m , 0 < b k,j for 1 ≤ k ≤ m and 2 ≤ j ≤ w k .Our idea is that the indices k correspond to different time scales δ b k τ . For m >
1, system (8), as δ →
0, may be thought of as separating fast variables from increasingly slow ones. It will turn out inSect. 3 that the exact number of time scales finally obtained by our overall approach can actually besmaller than m .Given certain conditions, which will be made explicit in Theorem 1 and with its application inSect. 3.1, we may formally truncate the right hand sides of (8) and keep only terms of lowest order in δ : z k = δ b k b f k ( z, , ≤ k ≤ m. (9)In the sequel, we refer to the transformation process from (1) to (8) as scaling . Strictly speaking, thiscomprises scaling in combination with partitioning . We refer to the step from (8) to (9) as truncating .Algorithm 1 reflects our discussions so far. It takes as input a list S of differential equations rep-resenting system (1) and a choice of 0 < ε ∗ < S as well as ε ∗ are taken from Q . Our algorithm is furthermore parameterized with afunction c mapping suitable indices to rational numbers and a constant function d yielding either a tuple D = ( d , . . . , d n ) of rational numbers or ⊥ . The black–box functions c and d reflect the mathematicalassumptions around (2) and (4) that suitable c k,J and d k exist, respectively. Suitable instantiationsfor the parameters c and d can be realized, e.g., using tropical geometry, which will be the topic ofSect. 2.2. It will turn out that instantiations of d can fail on the given combination of S and ε ∗ , whichis signaled by the return value ⊥ of d , and checked right away in l.1 of Algorithm 1. So far, the above transformations leading to (4) are a formal exercise. No particular strategy wasapplied for choosing ε ∗ ∈ (0 , ε ∗ as a power product in model parameters [28, 53].Here we discuss a different approach, based on tropical geometry [43, 44, 46, 47, 51, 50]. This approachstarts with a slightly different interpretation of the scaling problem. In this interpretation, the value ε ∗ is not dictated by physico-chemistry, but it is freely chosen to provide “power” parametric descriptionsof all the quantities occurring in the differential equations (parameters, monomials, time scales), ina similar way to describing curves by continuously varying real parameters. This interpretation isrooted in physics where it unravels scaling laws. It is also natural for any computation with orders ofmagnitude. In tropical geometry, it is encountered in several places: as Litvinov–Maslov dequantizationof real numbers leading to degeneration of complex algebraic varieties into tropical varieties [37, 61], orin the theory of Puiseux series in relation to tropical varieties and pre-varieties [4].The abstract scaling procedure leading to (6) is implemented with two additional requirements.Firstly, the orders c k,J are not freely chosen, but are computed from ε ∗ ∈ (0 ,
1) and γ k,J as c k,J = round( p log ε ∗ | γ k,J | ) p . (10)5 lgorithm 1 ScaleAndTruncate
Input:
1. A list S = ( d y d t = f , . . . , d y n d t = f n ) of autonomous first-order ordinary differentialequations where f , . . . , f n ∈ Q [ y , . . . , y n ];2. c : { , . . . , n } × { , . . . , n } n → Q ;3. d : () → Q n ∪ {⊥} ;4. ε ∗ ∈ (0 , ∩ Q Output:
1. A list ( T , . . . , T m ) where, abbreviating dd τ by a prime, T k = ( z k = δ b k f k ) with z k ⊆ ( x , . . . , x n ), S k z k = ( x , . . . , x n ), z , . . . , z m pairwise disjoint, b < · · · < b m ∈ N , and f k ⊆ Q [ x , . . . , x n ], or the empty list;2. A list ( P , . . . , P m ) of lists with P k ⊆ Q [ x , . . . , x n ][ δ ] and | P k | = | T k | for k ∈ { , . . . , m } ;3. A substitution σ for x , . . . , x n , τ , δ , and ε The first output ( T , . . . , T m ) contains differential equations z k = δ b k b f k ( z,
0) for k ∈ { , . . . , m } in terms of system (8). The second output ( P , . . . , P m ) contains the higher order terms in (8) aspolynomials p k = δ b k + b k, p k, + · · · + δ b k + b k,wk p k,w k . The last output is a substitution that undoesall substitutions applied for obtaining (8) from (1).This gives the following invariant: Denote e S = ( S mk =1 T k ⊕ P k ) σ , where ( x = g ) ⊕ p stands for x = g + p and is applied elementwise. Then e S is equal to S up to multiplication of the differentialequation ˙ y i = P J γ i,J y J in S with a positive scalar factor 1 /ε µ + d i ∗ .For q ∈ Q [ x , . . . , x n ]( δ ) we use deg δ ( q ) for the univariate degree of q in δ . Similarly, tmon δ ( q ) isthe trailing monomial in δ . if d () = ⊥ then return (), (), [ ] end if µ := ∞ q := 1 ( d , . . . , d n ) := d () ∈ Q n for k := 1 to n do h k := 0 for all monomials γy J in f k do ¯ γ := γ/ε c ( k,J ) ∗ ∈ Q η := c ( k, J ) + h ( d , . . . , d n ) , J i − d k ∈ Q µ := min( µ, η ) ∈ Q q := lcm( q, denom η ) ∈ N \ { } h k := h k + ε η ¯ γx J end for end for for k := 1 to n do h k := h k /ε µ h k := h k [ ε ← δ q ] ∈ Q [ x , . . . , x n ][ δ ] g k := tmon δ h k p k := h k − g k end for L := ( d x d τ = g , . . . , d x n d τ = g n ) ( b , . . . , b m ) := sort(deg δ g , . . . , deg δ g n ), ascending and removing duplicates for k := 1 to m do T k := ( d x d τ = g ∈ L | deg δ g = b k ) P k := ( p j ∈ { p , . . . , p n } | deg δ g j = b k ) end for σ := [ x ← y /ε d , . . . , x n ← y n /ε d n ] ◦ [ τ ← ε µ t ] ◦ [ δ ← ε /q ] ◦ [ ε ← ε ∗ ] return ( T , . . . , T m ), ( P , . . . , P m ), σ lgorithm 2 TropicalC
Input: k ∈ { , . . . , n } ;2. J ∈ { , . . . , n } n ;3. A list S = ( ˙ y = f , . . . , ˙ y n = f n ) of autonomous first-order ordinary differential equationswhere f , . . . , f n ∈ Q [ y , . . . , y n ];4. ε ∗ ∈ (0 , ∩ Q .5. p ∈ N \ { } Output: c ∈ Q γ := coeff( f k , y J ) ∈ Q c := round( p log ε ∗ | γ | ) /p ∈ Q return c Algorithm 3
TropicalD
Input:
1. A list S = ( ˙ y = f , . . . , ˙ y n = f n ) of autonomous first-order ordinary differential equationswhere f , . . . , f n ∈ Q [ y , . . . , y n ];2. ε ∗ ∈ (0 , ∩ Q .3. p ∈ N \ { } Output: D ∈ Q n ∪ {⊥} Π( a , . . . , a n ) := TropicalEquilibration( S, ε ∗ , p ) if not R | = ∃ a . . . ∃ a n Π then return ⊥ end if ( d , . . . , d n ) := one possible choice for a , . . . , a n return ( d , . . . , d n )Here, the rounding function rounds to nearest integer. The positive integer p controls the precision ofthe rounding step.Secondly, the orders D = ( d , d , . . . , d n ) satisfy certain constraints. These constraints result heuris-tically from the idea of compensation of dominant monomials [43]. Slow dynamics is possible if for eachdominant, i.e., much larger than the other, monomial on the right hand side of (6), there is at least oneother monomial of the same order, but with opposite sign. This condition, named tropical equilibrationcondition [43, 44, 46, 47, 51, 50], readsmin γ k,J > ( c k,J + h D, J i ) = min γ k,J < ( c k,J + h D, J i ) . (11)On these grounds, given system (1), the choice of ε ∗ boils down to defining orders of magnitude.Model parameters are coarse-grained and transformed to orders of magnitude in order to apply tropicalscaling. The result depends on which parameters are close and which are very different as dictated bythe coarse-graining procedure, i.e., by the choice of ε ∗ . Decreasing ε ∗ destroys details and parameterstend to have the same order of magnitude. Increasing ε ∗ refines details and parameters range overseveral orders of magnitude. For instance, using (10) and p = 1 parameters k = 0 . k = 0 . c = 1 and c = 2 for ε ∗ = 1 /
10, but c = c = 1 for ε ∗ = 1 /
50. This is the perspectivetaken in [43, 44, 51].It is noteworthy that in the context of singular perturbation methods (cf. Sect. 3), which provideasymptotic results as a small parameter approaches zero, there are independent arguments for choosing ε ∗ rather small.We are now ready to instantiate the black-box functions c and d in our generic Algorithm 1 withtropical versions as given in Algorithm 2 and Algorithm 3, respectively. To be precise, we use the IEEE 754 rounding rule round to nearest, ties to even. lgorithm 4 TropicalEquilibration
Input:
1. A list S = ( ˙ y = f , . . . , ˙ y n = f n ) of autonomous first-order ordinary differential equationswhere f , . . . , f n ∈ Q [ y , . . . , y n ];2. ε ∗ ∈ (0 , ∩ Q ;3. p ∈ N \ { } . Output:
A formula Π( a , . . . , a n ) describing a finite union of convex polyhedra in R n .We use h· , ·i to denote the standard scalar product in Q n +1 . A := (1 , a , . . . , a n ) ∈ Q [ a , . . . , a n ] n +1 for j := 1 to n do c := 0 for all monomials γy α · · · y α n n in f j do α := round( p log ε ∗ | γ | ) /p ∈ Q c := c + 1 Σ c := sgn γ ∈ {− , , } A c := ( α , α , . . . , α n ) ∈ Q × Z n ⊆ Q n +1 end for B j := ∅ for k := 1 to c do for ‘ := k + 1 to c do if Σ k Σ ‘ < then P := {h A k − A ‘ , A i = 0 } h A k − A ‘ , A i ∈ Q [ a , . . . , a n ] for m := 1 to c do P := P ∪ {h A m − A k , A i ≥ } h A m − A k , A i ∈ Q [ a , . . . , a n ] end for B j := B j ∪ { P } set of sets of constraints end if end for end for end for Π := DisjunctiveNormalForm( V nj =1 W P ∈ B j V P ) return ΠAlgorithm 2 explicitly uses, besides the parameters k and J specified for c in Algorithm 1, also theright hand sides of the input system (1) and the choice of ε ∗ . As yet another parameter it takes thedesired precision p for rounding in (10). Notice that the use of this extra information is compatiblewith the abstract scaling procedure in Section 2.1. Currying [16] allows to use Algorithm 2 in place of c in a formally clean manner.Similarly, Algorithm 3 takes parameters ε ∗ and p , while d is specified in Algorithm 1 to have noparameters at all. In l.1 we use Algorithm 4 as a subalgorithm for tropical equilibration. One obtainsa disjunctive normal form Π, which explicitly describes a set P = { p ∈ Q n | Π( p ) } as a finite unionof convex polyhedra, as known from tropical geometry. Every ( d , . . . , d n ) ∈ P satisfies (11). Thesatisfiability condition in l.2 tests whether P 6 = ∅ . We employ Satisfiability Modulo Theories (SMT) solving [41] using the logic
QF_LRA [2] for quantifier-free linear real arithmetic. The set P can get empty,e.g, when all monomials on the right hand side of some differential equation have the same sign. Suchan exceptional situation is signaled with a return value ⊥ in l.3. In the regular case P 6 = 0, the choice( d , . . . , d n ) in l.5 is provided by the SMT solver. From a practical point of view, the disjunctive normalform computation in Algorithm 4 is a possible bottleneck and requires good heuristic strategies [39].With applications in the natural sciences one often wants to make in l.5 an adequate choice for( d , . . . , d n ) lying in a specific convex polyhedron P ⊆ P , which technically corresponds to one conjunc-tion in Π. Such choices are subtle and typically require human interaction. For instance, when the chainof reduced dynamical systems ends with a steady state, it is interesting to consider the polyhedron P
3. Singular Perturbation Methods
The theory of singular perturbations is used to compute and justify theoretically the limit of system (8)when δ →
0. There are several types of results in this theory. The results of Tikhonov, further improvedby Hoppensteadt, show the convergence of the solution of system (8) to the solution of a differential-algebraic system in which the slowest variables z m follow differential equations and the remaining fastvariables follow algebraic equations [58, 29]. The results of Fenichel are known under the name of geometrical singular perturbations . He showed that the algebraic equations in Tikhonov’s theory definea slow invariant manifold that is persistent for δ > δ is needed in system (8).Samal et al. have noted that Tikhonov’s theorem is applicable to tropically scaled systems [51]. Forinstance, with δ = δ b , system (8) may be rewritten as z = b g ( z, δ ) , z = δ b g ( z, δ ) , . . . , z m = δ b g m ( z, δ ) . (12)However, this approach comes with certain limitations. To start with, it allows only two time scales.Furthermore, in case b >
1, there may be differentiability issues with respect to δ , and some care hasto be taken when one tries to apply to (12) also Fenichel’s results [21]. In this section, we are goingto generalize geometrical singular perturbations, and compute invariant manifolds and reduced modelsfor more than two time scales, introducing further δ , . . . , δ ‘ . Our generalization is based on a recentpaper by Cardin and Teixeira [10].Section 3.1 presents relevant results from [10] adapted to our purposes here and applied to our system(8). In contrast to the original article, which is based on a series of hyperbolicity conditions, we introducethe stronger notion of hyperbolic attractivity. In Sect. 3.2 we describe efficient algorithmic tests forhyperbolic attractivity. Section 3.3 gives sufficient algorithmic criteria addressing the above-mentioneddifferentiability issues. From now on we consider our system (8) over the positive first orthant U = (0 , ∞ ) n ⊆ R n . A recentpaper by Cardin and Teixeira [10] generalizes Fenichel’s theory to provide a solid foundation to obtainmore than one nontrivial invariant manifold. This allows, in particular, the reduction of multi-timescale systems such as system (8). Technically, the approach considers a multi-parameter system usingtime scale factors δ , δ δ , . . . instead of increasing powers of one single δ .We let ‘ ∈ { , . . . , m } and define β = b − b = b , . . . , β ‘ − = b ‘ − b ‘ − , (13)and furthermore δ = δ β , . . . , δ ‘ − = δ β ‘ − , and ¯ δ = ( δ , . . . , δ ‘ − ).These definitions allow us to express also all δ b k,j occurring in (8) as products of powers of δ , . . . , δ ‘ − ,with nonnegative but possibly non-integer rational exponents, via expressing each b k,j as a nonnegativerational linear combination of β , . . . , β ‘ − . This yields b g k ( z, δ , . . . , δ ‘ − ) = b f k ( z, δ ) , ≤ k ≤ m. (14)Moreover, we express δ b ‘ +1 = δ · · · δ ‘ − · η ‘ +1 , . . . , δ b m = δ · · · δ ‘ − · η m , via η k ( δ , . . . , δ ‘ − ) = δ b k − b ‘ , ‘ + 1 ≤ k ≤ m, (15)9hich is obtained by writing each b k − b ‘ as a nonnegative rational linear combination of β , . . . , β ‘ − .In these terms our system (8) translates to z = b g ( z, ¯ δ ) z = δ b g ( z, ¯ δ )... z ‘ = δ · · · δ ‘ − b g ‘ ( z, ¯ δ ) z ‘ +1 = δ · · · δ ‘ − η ‘ +1 (¯ δ ) b g ‘ +1 ( z, ¯ δ )... z m = δ · · · δ ‘ − η m (¯ δ ) b g m ( z, ¯ δ ) . (16)In terms of the right hand sides of (16) the application of relevant results in [10] requires that b g ,. . . , b g ‘ and η ‘ +1 b g ‘ +1 , . . . , η m b g m are smooth on an open neighborhood of U × [0 , ϑ ) × · · · × [0 , ϑ ‘ − )with ϑ >
0, . . . , ϑ ‘ − >
0. We are going to tacitly assume such smoothness here and address this issuefrom an algorithmic point of view in Sect. 3.3.We are now ready to transform our system into ‘ time scales as follows, where possibly ‘ > τ = τ, τ = δ τ, . . . , τ ‘ = δ · · · δ ‘ − τ. In time scale τ k , with 1 ≤ k ≤ ‘ , system (16) then becomes δ · · · δ k − d z d τ k = b g ( z, ¯ δ )... δ k − d z k − d τ k = b g k − ( z, ¯ δ )d z k d τ k = b g k ( z, ¯ δ )d z k +1 d τ k = δ k b g k +1 ( z, ¯ δ )...d z ‘ d τ k = δ k · · · δ ‘ − b g ‘ ( z, ¯ δ )d z ‘ +1 d τ k = δ k · · · δ ‘ − η ‘ +1 (¯ δ ) b g ‘ +1 ( z, ¯ δ )...d z m d τ k = δ k · · · δ ‘ − η m (¯ δ ) b g m ( z, ¯ δ ) . (17)For k = 1 and k = ‘ we obtain empty products, which yield the neutral element 1, as usual.Similarly to Sect. 2.1, we are interested in the asymptotic behavior for ¯ δ →
0, which is approximatedby the elimination of higher order terms. We are now going to introduce a construction required fora justification of this approximation, which also clarifies the greatest possible choice for ‘ ≤ m above.Define F = 0 and Z k = z ... z k , F k ( z, δ ) = b f ( z, δ )... b f k ( z, δ ) , ≤ k ≤ m. b f k ( z,
0) in favor of b g k ( z, , . . . , M k = (cid:0) F k ( z ∗ ,
0) = 0 (cid:1) , M k = { z ∗ ∈ U | F k ( z ∗ ,
0) = 0 } , ≤ k ≤ m. (18)The sets M k are obtained from varieties defined by the lists M k via intersection with the first orthant.Furthermore, U = M ⊇ M ⊇ · · · ⊇ M m establishes a chain of nested subvarieties, again intersectedwith the first orthant.We say that M is hyperbolically attractive on M , if M = ∅ and for all z ∈ M all eigenvaluesof the Jacobian D z b f ( z,
0) have negative real parts. Therefore M is a manifold. For k ∈ { , . . . , m } , M k is hyperbolically attractive on M k − , if M k = ∅ and the following holds. Recall that using thedefining polynomials F k − of M k − , the implicit function theorem yields a unique local resolution of Z k − as functions of z k , . . . , z m , and thus we obtain b f k ( z,
0) = b f ∗ k ( z k , . . . , z m ,
0) on M k − . We now require that for all z ∈ M k all eigenvalues of D z k b f ∗ k ( z k , . . . , z m ,
0) have negative real parts.Again, M k is a manifold. When M k is hyperbolically attractive on M k − we write M k − . M k , where Z k will be clear from the context.If we find for some ‘ ∈ { , . . . , m } that M . M , M . M , . . . , M ‘ − . M ‘ , then we simply write M . · · · . M ‘ , and call this a hyperbolically attractive ‘ -chain . Such a chain is called maximal if either ‘ = m or M ‘ . M ‘ +1 .Let M . · · · . M ‘ be a hyperbolically attractive ‘ -chain. Consider for each k ∈ { , . . . , ‘ } thefollowing differential-algebraic system:0 = F k − ( z, , d z k d τ k = b f k ( z, , d z k +1 d τ k = 0 , . . . , d z m d τ k = 0 . (19)In the limiting case ¯ δ = 0, this corresponds to system (17). Recall that τ k = δ · · · δ k − τ = δ b − b · · · δ b k − b k − τ = δ b k − b τ = δ b k τ, and equivalently rewrite (19) as a triplet ( M k − , T k , R k ) with entries as follows: F k − ( z,
0) = 0 , d z k d τ = δ b k b f k ( z, , d z k +1 d τ = · · · = d z m d τ = 0 . (20)For a given index k , we call ( M k − , T k , R k ) a reduced system on M k − , where the relevant hyperbolicattractivity relation is M k − . M k . In order to indicate the relevance of M . · · · . M ‘ we write( M , T , R ) . · · · . ( M ‘ − , T ‘ , R ‘ ) also for reduced systems, where M ‘ is not made explicit but relevantfor the last triplet. Slightly abusing language, we speak of a hyperbolically attractive ‘ -chain of reducedsystems, which is maximal if M . · · · . M ‘ is.The following theorem is a consequence of [10, Theorem A and Corollary A], specialized to thesituation at hand. Theorem 1.
Let ‘ ≥ . Assume that ( M , T , R ) . · · · . ( M ‘ − , T ‘ , R ‘ ) is a hyperbolically attractive ‘ -chain of reduced systems for system (16) . Let K ⊆ U be compact. Then for sufficiently small ¯ δ and all k ∈ { , . . . , ‘ } , system (16) admits invariant manifolds N k − that depend on ¯ δ and are ( δ + · · · + δ k − ) -close to M k − ∩ K with respect to the Hausdorff distance. Moreover, there exists T > such thatsolutions of system (16) on N k − in time scale τ k converge to solutions of ( M k − , T k , R k ) , uniformlyon any closed subinterval of (0 , T ) , as ¯ δ → . For k ∈ { , . . . , ‘ } , the M k − are critical manifolds, which contain only stationary points. Thesystems ( T k , R k ) of ordinary differential equations on M k − approximate invariant manifolds N k − inthe sense of the theorem. They furthermore approximate solutions in time scale τ k of system (16), whichis equivalent to our system (8). In other words, system (8) admits a succession of invariant manifolds,11 lgorithm 5 ComputeReducedSystems
Input:
Output of Algorithm 1:1. ( T , . . . , T m ), a list of lists z k = δ b k f k ;2. ( P , . . . , P m ), a list of lists of polynomials in Q [ x , . . . , x n ][ δ ];We denote ξ k := | T k | , Ξ k := P ki =1 ξ i , and X = ( x , . . . , x n ). Output:
A list (( M , T , R ) , . . . , ( M ‘ − , T ‘ , R ‘ )) of triplets where ‘ ∈ { , . . . , m } , or the empty list.For k ∈ { , . . . , ‘ } , M k − is a list of real constraints defining M k − ⊆ R n ; T k is a list of differentialequations; R k is a list of trivial differential equations x = 0 for all differential variables from T k +1 ,. . . , T m .The triplets ( M k − , T k , R k ) represent reduced systems according to (20). U := ( x > , . . . , x n > M , Z, F, A := () for k := 1 to m do z := ( x | x = δ b k g ∈ T k ) ⊆ X , | z | = ξ k f := ( g | x = δ b k g ∈ T k ) = ˆ f k ( z, ∈ Q [ X ] ξ k M k := M k − ◦ ( f = 0) = M ◦ ( F = 0) ◦ ( f = 0) ϕ, A := IsHyperbolicallyAttractive( U ◦ M k , Z, z, F, f, k, A ) if not ϕ then break end if R k := ( x = 0 | x = h ∈ T k +1 ∪ · · · ∪ T m ) Ξ k − + ξ k + | R k | = n Z := Z ◦ z ⊆ X , | Z | = Ξ k F := F ◦ f ∈ Q [ X ] Ξ k end for k , or we have k = m + 1 . ‘ := k − if ‘ < then return () end if if TestSmoothness(( T , . . . , T m ) , ( P , . . . , P m ) , ‘ ) = failed then print "Warning: differentiability requires further verification" end if return (( M , T , R ), . . . , ( M ‘ − , T ‘ , R ‘ ))on which the behavior in the appropriate time scale is approximated by the respective reduced equations(19) and, equivalently, (20). Note that only the δ b k b f k ( z,
0) without the higher order terms enter thereduced systems ( M k − , T k , R k ).Algorithm 5 now starts with the output ( T , . . . , T m ) of Algorithm 1, which represents the scaledsystem (9). Notice that each T k already meets the specification in (20). In l.1 we define U to containdefining inequalities of the first orthant U . Starting with k = 1, the for-loop in l.3–14 successivelyconstructs M k and R k such that in combination with T k from the input, ( M k − , T k , R k ) forms a reducedsystem as in (20). The loop stops when either k = m + 1 or a test for hyperbolic attractivity in l.7finds that M k − . M k . We are going to discuss this test in detail in the next Sect. 3.2. Notethat we maintain a matrix A for storing information between the subsequent calls of our test. Ineither case we arrive at a maximal hyperbolically attractive ( k − M , T , R ) , . . . , ( M k − , T k − , R k − )). Following the notational convention used throughoutthis section we set ‘ to k − ‘ ∈ { , . . . , m } at thebeginning of this section. Finally, l.20 uses the second input ( P , . . . , P m ) of the algorithm to addressthe smoothness requirements for system (16). We are going to discuss the corresponding procedure indetail in Sect. 3.3. It will turn out that this procedure provides only a sufficient test. Therefore we issue12n case of failure only a warning, allowing the user to verify smoothness a posteriori, using alternativealgorithms or human intelligence. One might mention that it is actually sufficient to consider weakerfinite differentiability conditions instead of smoothness, which can be seen by inspection of the proofsin [10].From an application point of view, attracting invariant manifolds are relevant in the context ofbiological networks, and our notion of hyperbolic attractivity holds for large classes of such networks[20]. This is our principal motivation for using hyperbolical attractivity here. From a computationalperspective, hyperbolic attractivity can be straightforwardly tested using the Hurwitz criterion, as weare going to make explicit in Sect. 3.2.The relevant results in [10], in contrast, are based on a series of hyperbolicity conditions, which aresomewhat weaker than hyperbolic attractivity. Hyperbolicity can be tested algorithmically as well,albeit with more effort. For approaches based on Routh’s work see, e.g., [23, Chapter V, §4], whichchecks the number of purely imaginary eigenvalues of a real polynomial via the Cauchy index of arelated rational function. Our definition of hyperbolic attractivity M k − . M k refers to the eigenvalues of the Jacobians of the b f ∗ k , which cannot be directly obtained from the Jacobians of the b f k [10, 11]. Generalizing work onsystems with three time scales [32], we take in this section a linear algebra approach to obtain therelevant eigenvalues without computing the b f ∗ k .To start with, recall the well-known Hurwitz criterion [30]:
Theorem 2 (Hurwitz, 1895) . Consider f = a x n + a x n − + · · · + a n ∈ R [ x ] , a > . For i ∈ { , . . . , n } define H i = a a a . . . a i − a a a . . . a i − a a . . . a i − . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . a i , ∆ i = | H i | . Then all complex zeros of f have negative real parts if and only if ∆ > , . . . , ∆ n > . Notice that ∆ n = a n ∆ n − , and therefore ∆ n > can be equivalently replaced with a n > . We call H n the Hurwitz matrix and ∆ i the i -th Hurwitz determinant of f . Furthermore, we refer toΓ = (∆ > ∧ · · · ∧ ∆ n − > ∧ a n >
0) as the
Hurwitz conditions for f .Our first result generalizes [32, Proposition 1 (ii)]. The proof is straightforward by induction. Lemma 3.
For k ∈ { , . . . , m } define J k = ... % ··· % k − ! · D Z k F k ( z,
0) = D z b f ( z, . . . D z k b f ( z, % D z b f ( z, . . . % D z k b f ( z, . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .% · · · % k − D z b f k ( z, . . . % · · · % k − D z k b f k ( z, . Let ‘ ∈ { , . . . , m } . Then M . · · · . M ‘ if and only if M ‘ = ∅ and for all k ∈ { , . . . , ‘ } , allsufficiently small % ∗ > , . . . , % ∗ k − > , and all z ∗ ∈ M k , all eigenvalues of J k ( % ∗ , . . . , % ∗ k − , z ∗ ) havenegative real parts.In particular, one can choose % ∗ = · · · = % ∗ k − = % ∗ with sufficiently small % ∗ and consider J k =diag(1 , . . . , % k − ) · D Z k F k ( z, . Let Γ k denote the Hurwitz conditions for the characteristic polynomial of J k . Then Lemma 3 allowsto state hyperbolic attractivity M . · · · . M ‘ as a first-order formula over the reals as follows: (cid:18) ∃ (0 < z ) : F ‘ ( z,
0) = 0 (cid:19) ∧ (cid:18) ‘ V k =1 ∃ (0 < σ ) ∀ (0 < % < σ ) ∀ (0 < z ) : F k ( z,
0) = 0 ⇒ Γ k ( %, z ) (cid:19) . (21)13n these grounds, any real decision procedure [57, 13, 63] provides an effective test for hyperbolicattractivity. However, our formulation (21) uses a quantifier alternation ∃ σ ∀ % in its second part. Froma theoretical point of view, such a bounded alternation does not affect the asymptotic worst-casecomplexity, which remains single exponential [27]. From a practical point of view, we would like tocontinue using SMT solving over a quantifier-free logic. Our next result allows a suitable first-orderformulation without quantifier alternation. Its proof combines [32, Lemma 3] with our Lemma 3. Proposition 4 (Effective Characterization of Hyperbolically Attractive ‘ -Chains) . Define A = D z b f ( z, . For k ∈ { , . . . , m } define (cid:18) A k − B k C k V k (cid:19) = (cid:18) D Z k − F k − ( z, D z k F k − ( z, D Z k − b f k ( z, D z k b f k ( z, (cid:19) , and note that (cid:16) A k − B k C k V k (cid:17) = A k . Let ‘ ∈ { , . . . , m } . Then M . · · · . M ‘ if and only if(i) M ‘ = ∅ ,(ii) for all z ∗ ∈ M all eigenvalues of W ( z ∗ ) , where W = A , have negative real parts,(iii) for all k ∈ { , . . . , ‘ } and all z ∗ ∈ M k , A k − ( z ∗ ) is regular and all eigenvalues of W k ( z ∗ ) , where W k = V k − C k A − k − B k , have negative real parts.Proof. Assume M . · · · . M ‘ . By Lemma 3 we have M ‘ = ∅ . For all z ∗ ∈ M , all eigenvaluesof the Jacobian W ( z ∗ ) have negative real parts by the definition of hyperbolic attractivity. Let now k ∈ { , . . . , ‘ } , z ∗ ∈ M k , and define P = diag(1 , . . . , % k − ). Using Lemma 3 we fix 0 < τ ∗ < < % ∗ < τ ∗ all eigenvalues of J k − ( % ∗ , z ∗ ) = P( % ∗ ) A k − ( z ∗ ) have negative real parts. It followsthat P( % ∗ ) A k − ( z ∗ ), P( % ∗ ), and A k − ( z ∗ ) are all regular. Next, consider J k = (cid:18) P A k − P B k % k − C k % k − V k (cid:19) . Using Lemma 3 once more, we find 0 < σ ∗ < τ ∗ such that for all 0 < % ∗ < σ ∗ also all eigenvalues of J k ( % ∗ , z ∗ ) have negative real parts. Now J k ( % ∗ , z ∗ ) satisfies condition (ii) of [32, Lemma 3] with δ = σ ∗ and ε = ( % ∗ ) k − , which allows us to conclude that all eigenvalues of ( V k − C k (P A k − ) − P B k )( % ∗ , z ∗ ) =( V k − C k A − k − P − P B k )( % ∗ , z ∗ ) = W k ( z ∗ ) have negative real parts as well.Assume, vice versa, that (i)–(iii) hold. We use induction on k to show M . · · · . M k for 1 ≤ k ≤ ‘ .For k = 1 we have M . M by definition of hyperbolic attractivity. Assume that 2 ≤ k ≤ ‘ and M . · · · . M k − . By Lemma 3 there exists 0 < τ ∗ such that for all 0 < σ ∗ < τ ∗ and all z ∗ ∈ M k − all eigenvalues of P( σ ∗ ) A k − ( z ∗ ) have negative real parts, where P = diag(1 , . . . , ( σ ∗ ) k − ).We rewrite W k = V k − C k (P A k − ) − P B k . Then W k ( z ∗ ) satisfies condition (i) of [32, Lemma 3] with A = P( σ ∗ ) A k − ( z ∗ ), B = P( σ ∗ ) B k ( z ∗ ), C = C k ( z ∗ ) and D = V k ( z ∗ ). Thus there exists 0 < δ such thatfor all 0 < ε < δ all eigenvalues of (cid:18) P( σ ∗ ) A k − ( z ∗ ) P( σ ∗ ) B k ( z ∗ ) εC k ( z ∗ ) εV k ( z ∗ ) (cid:19) have negative real parts. Choosing % ∗ = min { σ ∗ , k − √ ε } in Lemma 3 yields M . · · · . M k .From now on let Γ k denote the Hurwitz conditions for the characteristic polynomial of W k , which—incontrast to the ones used in (21)—do not depend on % anymore. Corollary 5 (Logic-Based Test for Hyperbolically Attractive ‘ -chains) . For k ∈ { , . . . , m } define ϕ k = (cid:0) ∃ (0 < z ) : F k ( z,
0) = 0) (cid:1) ,ψ k = (cid:0) ∀ (0 < z ) : F k ( z,
0) = 0 ⇒ Γ k ( z ) (cid:1) . Let ‘ ∈ { , . . . , m } . Then M . · · · . M ‘ if and only if R | = ϕ ‘ ∧ V ‘k =1 ψ k . lgorithm 6 IsHyperbolicallyAttractive
Input: M , 2. Z , 3. z , 4. F , 5. f , 6. k , 7. A , as in the calling Algorithm 5Knowing that M . · · · . M k − , we check here whether also M k − . M k . We denote ξ := | f | = | z | ,Ξ := | F | = | Z | , and X = ( x , . . . , x n ). In these terms, A ∈ Q [ X ] Ξ × Ξ . Output:
1. Boolean, 2. A ∈ Q [ X ] (Ξ+ ξ ) × (Ξ+ ξ ) if not R | = ∃ V M then return false, (cid:0) (cid:1) end if V := Jacobian( f, z ) ∈ Q [ X ] ξ × ξ if k = 1 then W := V A := V else B := Jacobian( F, z ) ∈ Q [ X ] Ξ × ξ C := Jacobian( f, Z ) ∈ Q [ X ] ξ × Ξ W := V − CA − B ∈ Q [ X ] ξ × ξ A := (cid:18) A BC V (cid:19) ∈ Q [ X ] (Ξ+ ξ ) × (Ξ+ ξ ) end if χ := λ ξ + · · · + a ξ := CharacteristicPolynomial( W ) ∈ Q [ X ][ λ ] H := HurwitzMatrix( χ ) ∈ Q [ X ] ξ × ξ for j := 1 to ξ − do ∆ j := det (cid:0) H r,s (cid:1) ≤ r,s ≤ j ∈ Q [ X ] end for Γ := { ∆ > , . . . , ∆ ξ − > , a ξ > } return R | = ∀ ( V M −→ V Γ), A Proof.
Assume M . · · · . M ‘ . Then Proposition 4 yields its conditions (i)–(iii). Now, ϕ ‘ holds asa formalization of (i). Furthermore, ψ holds as a formalization of (ii), and the validity of ψ , . . . , ψ ‘ follows directly from (iii). Hence R | = ϕ ‘ ∧ V ‘k =1 ψ k .Assume, vice versa, that R | = ϕ ‘ ∧ V ‘k =1 ψ k . We show M . · · · . M ‘ by induction on ‘ . If ‘ = 1,then ϕ formalizes (i) and ψ formalizes (ii) in Proposition 4, and we obtain M . M . Let now ‘ >
1. Then ϕ ‘ formalizes Proposition 4 (i). Our induction hypothesis yields M . · · · . M ‘ − . ByLemma 3 there exists 0 < τ ∗ such that for all 0 < σ ∗ < τ ∗ and all z ∗ ∈ M ‘ − ⊇ M ‘ all eigenvalues ofP( σ ∗ ) A ‘ − ( z ∗ ), where P = diag(1 , . . . , ( σ ∗ ) ‘ − ), have negative real parts. In particular, P( σ ∗ ) A ‘ − ( z ∗ )is regular and so is A ‘ − ( z ∗ ). Furthermore, the Hurwitz conditions in ψ ‘ guarantee for all z ∗ ∈ M ‘ that W ‘ ( z ∗ ) has only negative eigenvalues. Taking these observations together, Proposition 4 (iii) issatisfied, hence M . · · · . M ‘ .In contrast to (21), our first-order characterization (cid:18) ∃ (0 < z ) : F ‘ ( z,
0) = 0 (cid:19) ∧ (cid:18) ‘ V k =1 ∀ (0 < z ) : F k ( z,
0) = 0 ⇒ Γ k ( z ) (cid:19) (22)in Corollary 5 has no quantifier alternation. Note that the two top-level components of (22) establishtwo independent decision problems, addressing non-emptiness of the manifold and our requirement onthe eigenvalues, respectively.It is easy to see that for all ‘ ∈ { , . . . , m } and all k ∈ { , . . . , ‘ − } , ϕ ‘ entails ϕ k . Thus (22) can beequivalently rewritten as V ‘k =1 ( ϕ k ∧ ψ k ), explicitly: ‘ V k =1 (cid:0) ∃ (0 < z ) : F k ( z,
0) = 0 ∧ ∀ (0 < z ) : F k ( z,
0) = 0 ⇒ Γ k ( z ) (cid:1) . (23)15ur approach tests the conjunction in (23) using a for-loop over k in Algorithm 5. Technically, thisconstruction ensures with the test for M k − . M k in Algorithm 6 that M . · · · . M k − already holds,and exploits the fact that ψ k and ϕ k do not refer to smaller indices than k .In l.1–3 we test the validity of ϕ k . Using from the input the defining inequalities and equations M = U ◦ M k of M k along with Z = Z k − , z = z k , F = F k − , f = f k , and A = A k − , we constructin l.4–13 A = A k as noted in Proposition 4. In l.14–19 we construct the Hurwitz conditions Γ = Γ k according to Theorem 2. On the grounds of the validity of ϕ k tested in l.1, we finally test in l.20the validity of ψ k and return a corresponding Boolean value. We additionally return A = A k forreuse with the next iteration. The validity tests for ϕ k and ψ k in l.1 and l.20, respectively, amount toSMT solving, this time using the logic QF_NRA [2] for quantifier-free nonlinear real arithmetic. Recallthe positive integer parameter p used for the precision with both Algorithm 2 and Algorithm 3. For p > QF_NRA . When this happens, we catch the corresponding error from the SMT solver andrestart with floats.
Let us get back to the requirement in Sect. 3.1 that b g , . . . , b g ‘ and η ‘ +1 b g ‘ +1 , . . . , η m b g m occurring on theright hand sides of system (16) are all smooth on an open neighborhood of U × [0 , ϑ ) × · · · × [0 , ϑ ‘ − )with ϑ >
0, . . . , ϑ ‘ − >
0. An obvious criterion for smoothness is that all those expressions arepolynomials in z and ¯ δ .Recall the definitions of b g k for k ∈ { , . . . , m } in (14) and of η k for k ∈ { ‘ + 1 , . . . , m } in (15). For k ∈ { , . . . , m } and j ∈ { , . . . , w k } one finds nonnegative r , . . . , r ‘ − ∈ Q such that h ( β , . . . , β ‘ − ) , ( r , . . . , r ‘ − ) i = b k,j , and for k ∈ { ‘ + 1 , . . . , m } one finds nonnegative r , . . . , r ‘ − ∈ Q such that h ( β , . . . , β ‘ − ) , ( r , . . . , r ‘ − ) i = b k − b ‘ . Such representations always exist but are not unique in general. If one even finds suitable nonnegativeintegers r , . . . , r ‘ − ∈ N , which do not always exist, then one obtains b g , . . . , b g m as polynomials in z and ¯ δ , and η ‘ +1 , . . . , η m as polynomials in ¯ δ , which is sufficient for our criterion above.An improved but still only sufficient criterion uses similar constructions to directly verify the existenceof polynomial representations of the products η ‘ +1 b g ‘ +1 , . . . , η m b g m , in contrast to considering the factorsindependently. From an algorithmic point of view, we furthermore have to take into account that P ,. . . , P m obtained in Algorithm 1 do not contain b k,j but b k + b k,j . For k ∈ { , . . . , ‘ } and j ∈ { , . . . , w k } we try to find r , . . . , r ‘ − ∈ N such that h ( β , . . . , β ‘ − ) , ( r , . . . , r ‘ − ) i = b k,j = ( b k + b k,j ) − b k > , (24)and for k ∈ { ‘ + 1 , . . . , m } we try to find r , . . . , r ‘ − ∈ N such that h ( β , . . . , β ‘ − ) , ( r , . . . , r ‘ − ) i = ( b k − b ‘ ) + b k,j = ( b k + b k,j ) − b ‘ > . (25)Notably, such representations exist whenever 1 ∈ { β , . . . , β ‘ − } .On these grounds, we introduce Algorithm 7, which specifies the sufficient test applied in l.20 ofAlgorithm 5. The first two parameters ( T , . . . , T m ) and ( P , . . . , P m ) originate from Algorithm 1, whilethe last parameter ‘ originates from the calling Algorithm 5.In l.1–8 of Algorithm 7 we compute β , . . . , β ‘ − as defined in (13) and simultaneously obtain b ,. . . , b ‘ . In l.9–14 we compute the right hand sides of the conditions in (24) or (25), depending on thecurrent index k . For checking those conditions in l.16 we once more employ SMT solving, this timeusing the adequate logic QF_LIA [2] for quantifier-free linear integer arithmetic. Since we are aiming atnonnegative integer solutions, we introduce explicit non-negativity conditions r ≥
0, . . . , r ‘ − ≥
0. In16 lgorithm 7
TestSmoothness
Input: ( T , . . . , T m ), ( P , . . . , P m ), ‘ as in the calling Algorithm 5:1. ( T , . . . , T m ), a list of lists z k = δ b k f k ;2. ( P , . . . , P m ), a list of lists of polynomials in Q [ x , . . . , x n ][ δ ];3. ‘ ∈ N , ‘ ≥ Output: "true" or "failed" in terms of a 3-valued logic; b := 0 for k := 2 to ‘ do b k := the unique exponent of δ in T k β k − := b k − b k − if β k − = 1 then return true end if end for E := ∅ for k = 1 to m do for all p in P k do E := E ∪ { deg δ m − b min( k,‘ ) | m monomial of p } ⊆ N \ { } end for end for for all e ∈ E do if not Z | = ∃ r . . . ∃ r ‘ − ( r ≥ ∧ · · · ∧ r ‘ − ≥ ∧ h ( β , . . . , β ‘ − ) , ( r , . . . , r ‘ − ) i = e ) then return failed end if end for return truecase of unsatisfiability Algorithm 7 returns “failed” in l.17. When this happens, the calling Algorithm 5issues a warning but continues. This protocol is owed to the fact that our procedure provides onlya sufficient test, which could be supplemented with other software or human intuition. In case ofsatisfiability, in contrast, smoothness is guaranteed, we reach l.20, and return “true.” We remark thatthe computation time spent on E is negligible compared to the SMT solving later on. The constructionof the entire set E beforehand avoids duplicate SMT instances.
4. Algebraic Simplification of Reduced Systems
In the output ( M , T , R ), . . . , ( M ‘ − , T ‘ , R ‘ ) of Algorithm 5, the T k are taken literally from the input,and the M k − and R k are obtained via quite straightforward rewriting of the input. As a matter offact, the computationally hard part of Algorithm 5 consists in the computation of the upper index ‘ . We now want to rewrite the triplets ( M k − , T k , R k ) once more, aiming at less straightforward butsimpler and, hopefully, more intuitive representations. The principal idea is to heuristically eliminateon the right hand side of the differential equations in T k those variables whose derivatives have alreadyoccurred as left hand sides in one of the T , . . . , T k − . Of course, our simplifications will preserve allrelevant properties of ( M , T , R ), . . . , ( M ‘ − , T ‘ , R ‘ ), such as hyperbolic attractivity and sufficientdifferentiability. Technically, our next Algorithm 8 employs Gröbner basis techniques [9, 3].Recall that z k are the variables occurring on the left hand sides of differential equations in T k , and Z k − = ( z , . . . , z k − ). In l.1–5 we construct a block term order ω on all variables { x , . . . , x n } so thatvariables from Z k − are larger than variables from z k . This ensures that all multivariate polynomialreductions modulo ω throughout our algorithm will eliminate variables from Z k − in favor of variables17 lgorithm 8 SimplifyReducedSystems
Input:
A list (( M , T , R ) , . . . , ( M ‘ − , T ‘ , R ‘ )), the output of Algorithm 5, with entries correspondingto (20) Output:
A list (( M , T , R ) , . . . , ( M ‘ − , T ‘ , R ‘ )); M k − describes the same manifold as M k − in acanonical form; the system T k is equivalent to T k modulo M k − , its right hand sides are in acanonical normal form modulo M k − , possibly with fewer different differential variables than T k for k := 1 to ‘ do z k := { x | x = g ∈ T k } end for y := { x | x = 0 ∈ R ‘ } ω := a block term order with z (cid:29) · · · (cid:29) z ‘ (cid:29) y for k := 1 to ‘ do F := ( f | f = 0 ∈ M k − ) G := GroebnerBasis(Radical( F ) , ω ) M k − := ( g = 0 | g ∈ G ) T k := () for x = g in T k do T k := T k ◦ ( x = h ) where g −→ ∗ G h and h is irreducible mod G end for end for return (( M , T , R ) , . . . , ( M ‘ − , T ‘ , R ‘ ))from z k rather than vice versa. Prominent examples for such block orders are pure lexicographical orders,but ordering by total degree inside the z , . . . , z ‘ , y will heuristically give more efficient computations.Recall that the radical ideal p h F i of F is the infinite set of all polynomials with the same commoncomplex roots as F . In l.8, we compute a finite reduced Gröbner basis G modulo ω of that radical. Ifradical computation is not available on the software side, then the algorithm remains correct with aGröbner basis of the ideal h F i instead of the radical ideal, but might miss some simplifications.In l.9, the polynomials in G equivalently replace the left hand side polynomials of the equations in M k − . In l.12, reduction modulo ω , which comes with heuristic elimination of variables, applies oncemore to the reduction results h obtained from right hand sides g of differential equations in T k . Since G is a Gröbner basis, the reduction in l.11–13 furthermore produces unique normal forms with thefollowing property: if two polynomials g , g coincide on the manifold M k − defined by M k − , thenthey reduce to the same normal form h . In particular, if g vanishes on M k − , then it reduces to 0.We call the output of Algorithm 8 simplified reduced systems .
5. Back-Transformation of Reduced Systems
Let ‘ ∈ { , . . . , m } and k ∈ { , . . . , ‘ } . Recall that a triplet ( M k − , T k , R k ) obtained from Algorithm 5describes a reduced system according to (20). A corresponding simplified system ( M k − , T k , R k ) isobtained from Algorithm 8 via an equivalence transformation on the set of equations M k − and furtherequivalence transformations modulo M k − on the right hand sides of the differential equations in T k ,while the left hand sides of those differential equations remain untouched. It is not hard to see that forboth these outputs scaling can be reversed using the substitution σ = [ x ← y /ε d , . . . , x n ← y n /ε d n ] ◦ [ τ ← ε µ t ] ◦ [ δ ← ε /q ] ◦ [ ε ← ε ∗ ]obtained with Algorithm 1. For our discussion here, we use names M k − , T k , R k as in the unsimplifiedsystem.The application of σ to all components of ( M k − , T k , R k ) yields a raw back-transformation as follows: M k − σ = ( f j = 0 | x j ∈ Z k − ) σ = ( f j σ = 0 | x j ∈ Z k − ) , lgorithm 9 TransformBack
Input:
1. (( M , T , R ) , . . . , ( M ‘ − , T ‘ , R ‘ )), the output of either Algorithm 5 or Algorithm 8;2. σ , the output of Algorithm 1 Output:
A list (( M ∗ , T ∗ , R ∗ ) , . . . , ( M ∗ ‘ − , T ∗ ‘ , R ∗ ‘ )). for k := 1 to ‘ do M ∗ k − := M k − σ v := (cid:0) ( δ b k τ ) σ (cid:1) /t , extracting δ b k from T k = ε ( b k /q )+ µ ∗ T ∗ k := () for all x j = δ b k f j ∈ T k do h := ( y j f j /x j ) σ = ε d j ∗ ( f j σ ) T ∗ k := T ∗ k ◦ ( ˙ y j = vh ) = ε b k /q + µ + d j ∗ ( f j σ ) end for R ∗ k := ( ˙ y j = 0 | x j = 0 ∈ R k ) end for return (( M ∗ , T ∗ , R ∗ ) , . . . , ( M ∗ ‘ − , T ∗ ‘ , R ∗ ‘ )) T k σ = (cid:18) d x j d τ = δ b k f j (cid:12)(cid:12)(cid:12)(cid:12) x j ∈ z k (cid:19) σ = (cid:18) d y j d ε µ + d j ∗ t = ε b k /q ∗ ( f j σ ) (cid:12)(cid:12)(cid:12)(cid:12) x j ∈ z k (cid:19) ,R k σ = (cid:18) d x j d τ = 0 (cid:12)(cid:12)(cid:12)(cid:12) x j ∈ z k +1 ∪ · · · ∪ z m (cid:19) σ = (cid:18) d y j d ε µ + d j ∗ t = 0 (cid:12)(cid:12)(cid:12)(cid:12) x j ∈ z k +1 ∪ · · · ∪ z m (cid:19) . In T k σ , we multiply by ε µ + d j ∗ in order to arrive at differential equations in d y j d t . Furthermore, recallthat the explicit factor δ b k in the original T k corresponds to a time scale δ b k τ . The corresponding timescale in t is given by ( δ b k τ ) σ = ε b k /q + µ ∗ t , which we make explicit by equivalently rewriting T k σ as T ∗ k = (cid:16) ˙ y j = ε b k /q + µ ∗ ( ε d j ∗ f σ ) (cid:12)(cid:12)(cid:12) x j ∈ z k (cid:17) . Similarly, R k σ can be rewritten as R ∗ k = ( ˙ y j = 0 | x j ∈ z k +1 ∪ · · · ∪ z m ), and we set M ∗ k = M k σ .We call ( M ∗ , T ∗ , R ∗ ), . . . , ( M ∗ ‘ − , T ∗ ‘ , R ∗ ‘ ) back-transformed reduced systems . In terms of the defini-tions after (9) in Sect. 2.1 we have reverted the scaling but not the partitioning and not the truncating.Furthermore, we have preserved all information obtained with the computation of the reduced systemsin Sect. 3, where we keep the time scale factors explicit, and with their algebraic simplification in Sect. 4.Our back-transformation is realized in Algorithm 9. In l.3 we compute the time scale factor ε ( b k /q )+ µ ∗ for T ∗ k as described above, and in l.6 we compute its co-factor ε d j ∗ f σ as ( y j f /x j ) σ .Let us discuss what has been gained in (( M ∗ , T ∗ , R ∗ ) , . . . , ( M ∗ ‘ − , T ∗ ‘ , R ∗ ‘ )) for our original system S given in (1). To start with, notice that our decision at the beginning of Sect. 3.1 to limit ourselvesto the positive first orthant U using defining inequalities U = { x > , . . . , x n > } translates into U σ = { y /ε d ∗ > , . . . , y n /ε d n ∗ > } , which is again the positive first orthant U . The manifolds M k − described by M k − lead to manifolds M ∗ k − described by M ∗ k − , preserving the nestedness U = M ∗ ⊇ M ∗ ⊇ · · · ⊇ M ∗ ‘ − . Moreover, the system ( T ∗ k , R ∗ k ) defines differential equations on M ∗ k − .We are now faced with a discrepancy. On the one hand, we fix ε = ε ∗ . On the other hand, therequirement that ¯ δ be sufficiently small in Theorem 1 entails that ε be sufficiently small. It is ofcrucial importance whether invariant manifolds of (3), which do exist for sufficiently small ε , persist at ε = ε ∗ . We are not aware of any algorithmic results addressing this question. In particular, singularperturbation theory is typically concerned with asymptotic results, which are not helpful here.In case of persistence, there exist nested invariant manifolds N ∗ k − which are Hausdorff-close to M ∗ k − for system (1). Moreover, the differential equations T ∗ k associated with M ∗ m − correspond to the k th19 lgorithm 10 TropicalMultiReduce
Input:
1. A list S = ( ˙ y = f , . . . , ˙ y n = f n ) of autonomous first-order ordinary differential equationswhere f , . . . , f n ∈ Q [ y , . . . , y n ];2. ε ∗ ∈ (0 , ∩ Q ;3. p ∈ N \ { } Output:
A list (( M ∗ , T ∗ , R ∗ ) , . . . , ( M ∗ ‘ − , T ∗ ‘ , R ∗ ‘ )) of triplets where ‘ ∈ { , . . . , m } , or the empty list.For k ∈ { , . . . , ‘ } , M ∗ k − is a list of real constraints defining M ∗ k − ⊆ R n ; T ∗ k is a list of differentialequations; R ∗ k is a list of trivial differential equations ˙ y = 0 for all differential variables from T ∗ k +1 ,. . . , T ∗ m .The relevance of the output in terms of the input is discussed in Sect. 5. TropicalC
S,ε ∗ ,p := curry(TropicalC , S, ε ∗ , p ) TropicalC S,ε ∗ ,p is a binary function TropicalD
S,ε ∗ ,p := curry(TropicalD , S, ε ∗ , p ) TropicalD S,ε ∗ ,p is a constant function T, P, σ := ScaleAndTruncate( S, TropicalC
S,ε ∗ ,p , TropicalD
S,ε ∗ ,p , ε ∗ ) Σ := ComputeReducedSystems(
T, P ) = (( M , T , R ) , . . . , ( M ‘ − , T ‘ , R ‘ )) Σ := SimplifyReducedSystems(Σ) = (( M , T , R ) , . . . , ( M ‘ − , T ‘ , R ‘ )) Σ ∗ := TransformBack(Σ , σ ) = (( M ∗ , T ∗ , R ∗ ) , . . . , ( M ∗ ‘ − , T ∗ ‘ , R ∗ ‘ )) return Σ ∗ level in a hierarchy of time scales and approximate the flow on N ∗ k − . We have achieved a decompositionof (1) into ‘ systems of smaller dimension. At the very least, one obtains a well-educated guess aboutpossible candidates for invariant manifolds and reductions. For the investigation of those candidates onemay check the N ∗ k for approximate invariance using, e.g., numerical methods, or by applying criteriaproposed in [45].Algorithm 10 provides a wrapper combining all our algorithms to decompose input systems like (1)into several time scales. The underlying tropicalization is not made explicit, and the result is presentedon the original scale. Figure 5 explains the functional dependencies and principal data flow betweenour algorithms graphically.
6. Computational Examples
Based on our explicit algorithms in the present work, we have developed two independent softwareprototypes realizing all methods described here. The first one is in Python using SymPy [40] forsymbolic computation, pySMT [24] as an interface to the SMT solver MathSAT5 [12], and SMTcut forthe computation of tropical equilibrations [39]. The second one is a Maple package, which makes useof Maple’s built-in SMTLIB package [22] for using the SMT solver z3 [17]. For our computations herewe have used our Python code. Computation results are identical with both systems, and timings aresimilar. We have conducted our computations on a standard desktop computer with an 3.3 GHz 6-coreIntel 5820K CPU and 16 GB of main memory. Computation times listed are CPU times.In the next subsection, we will discuss in detail the computations for one specific biological input sys-tem from the BioModels database. The subsequent subsections showcase several further such examplesin a more concise style. The focus here is on biological results. For an illustration of our algorithms,we discuss in Appendix A examples where reduction stops at ‘ < m for various reasons.
We consider a model related to the transmission dynamics of subtype H5N6 of the influenza A virusin the Philippines in August 2017 [35]. That model is identified as
BIOMD0000000716 in the BioModelsdatabase, a repository of mathematical models of biological processes [34]. The model specifies fourspecies:
S_b (susceptible bird),
I_b (infected bird),
S_h (susceptible human), and
I_a (infected human),the concentrations of which over time we map to differential variables y , y , y , y , respectively. The20 lgorithm 10 TropicalMultiReduce
Algorithm 1
ScaleAndTruncate
Algorithm 5
ComputeReducedSystems
Algorithm 8
SimplifyReducedSystems
Algorithm 9
TransformBack
Algorithm 4
TropicalEquilibration
Algorithm 3
TropicalD
Algorithm 2
TropicalC
Algorithm 7
TestSmoothness
Algorithm 6
IsHyperbolicallyAttractive S , ε ∗ , pS , ε ∗ , pP , . . . , P m σ S = ( ˙ y = g , . . . , ˙ y n = g n )2. ε ∗ p S
2. TropicalC
S,ε ∗ ,p
3. TropicalD
S,ε ∗ ,p ε ∗ T , . . . , T m ,where T k is z k = δ b k b f k ( z, b < · · · < b m ( M , T , R ) . · · · . ( M ‘ − , T ‘ , R ‘ )( M , T , R ) . · · · . ( M ‘ − , T ‘ , R ‘ )( M ∗ , T ∗ , R ∗ ) . · · · . ( M ∗ ‘ − , T ∗ ‘ , R ∗ ‘ ) Figure 1: Functional dependencies (thin arrows) and principal data flow (thick arrows) between ouralgorithms 21nput system is given by S = (cid:2) dd t y = − y y − y + , dd t y = y y − y , dd t y = − y y − y + , dd t y = y y − y (cid:3) . We choose ε ∗ = , p = 1, and Algorithm 3 non-deterministically selects d = ( − , − , − , −
3) fromthe tropical equilibration. Algorithm 1 then yields the following scaled and truncated system with threetime scales: T = (cid:2) dd τ x = 1 · (cid:0) − x x + (cid:1)(cid:3) ,T = (cid:2) dd τ x = δ · (cid:0) x x − x (cid:1)(cid:3) ,T = (cid:2) dd τ x = δ · (cid:0) − x + (cid:1) , dd τ x = δ · (cid:0) x x − x (cid:1)(cid:3) . From this input, Algorithm 5 produces the following reduced systems: M = (cid:2) (cid:3) , T = (cid:2) dd τ x = 1 · (cid:0) − x x + (cid:1)(cid:3) ,R = (cid:2) dd τ x = 0 , dd τ x = 0 , dd τ x = 0 (cid:3) ,M = (cid:2) x x − (cid:3) , T = (cid:2) dd τ x = δ · (cid:0) x x − x (cid:1)(cid:3) ,R = (cid:2) dd τ x = 0 , dd τ x = 0 (cid:3) ,M = (cid:2) x x − , T = (cid:2) dd τ x = δ · (cid:0) − x + (cid:1) , x x − x = 0 (cid:3) , dd τ x = δ · (cid:0) x x − x (cid:1)(cid:3) ,R = (cid:2) (cid:3) . In that course, Algorithm 6 successfully tests all three scaled systems for hyperbolic attractivity. Fur-thermore, Algorithm 7 applies the sufficient smoothness test from Sect. 3.3 with ‘ = 3 , b = 3 , b = 3 , P = 1 · ( − δ · x ) , P = δ · ( − δ · x x ) . This yields E = { } , where 4 cannot be expressed as an integer multiple of 3. Thus the test fails,which causes a warning in Algorithm 5. Notice that in R , . . . , R ‘ the differential variables are orderedin the same way as in the scaled and truncated system T , . . . , T m . Incidentally, this coincides withlexicographic order in this example.Algebraic simplification through Algorithm 8 yields the simplified reduced systems M = (cid:2) (cid:3) , T = (cid:2) dd τ x = 1 · (cid:0) − x x + (cid:1)(cid:3) ,R = (cid:2) dd τ x = 0 , dd τ x = 0 , dd τ x = 0 (cid:3) ,M = (cid:2) x x = (cid:3) , T = (cid:2) dd τ x = δ · (cid:0) − x + (cid:1)(cid:3) ,R = (cid:2) dd τ x = 0 , d τ x = 0 (cid:3) ,M = (cid:2) x = , T = (cid:2) dd τ x = δ · (cid:0) − x + (cid:1) ,x = (cid:3) , dd τ x = δ · (cid:0) x − x (cid:1)(cid:3) ,R = (cid:2) (cid:3) . Notice that our implementations conveniently rewrite equational constraints as monomial equationswith numerical right hand sides when possible. This supports readability but is not essential for thesimplifications applied here, which are based on Gröbner basis theory. Comparing T with T , we seethat the equation for x x in M is plugged in. Similarly, M is simplified to M , which is in turn usedto reduce T to T .The back-transformed reduced systems as computed by Algorithm 9 read as follows: M ∗ = (cid:2) (cid:3) , T ∗ = (cid:2) dd t y = 1 · (cid:0) − y y + (cid:1)(cid:3) ,R ∗ = (cid:2) dd t y = 0 , dd t y = 0 , dd t y = 0 (cid:3) ,M ∗ = (cid:2) y y = (cid:3) , T ∗ = (cid:2) dd t y = · (cid:0) − y + (cid:1)(cid:3) ,R ∗ = (cid:2) dd t y = 0 , dd t y = 0 (cid:3) ,M ∗ = (cid:2) y = , T ∗ = (cid:2) dd t y = · (cid:0) − y + (cid:1) ,y = (cid:3) , dd t y = · (cid:0) y − y (cid:1)(cid:3) ,R ∗ = (cid:2) (cid:3) . We compare T ∗ , . . . , T ∗ to the input system S : In the equation for ˙ y , the monomial in y is identified asa higher order term with respect to δ and discarded by Algorithm 1. In the equation for ˙ y , the monomialin y y has been Gröbner-reduced to a constant modulo the defining equation in M . Similarly, theequation for ˙ y loses its monomial in y y by truncation of higher order terms, and in the equation for˙ y , the monomial in y y is Gröbner-reduced to a monomial in y .Notice the explicit constant factors on the right hand sides of the differential equations in T ∗ , . . . , T ∗ .They originate from factors δ b k in the respective scaled systems T , . . . , T , corresponding to (8). Theyare left explicit to make the time scale of the differential equations apparent. We see that the system T ∗ ◦ R ∗ is 125 times slower than T ∗ ◦ R ∗ , and T ∗ ◦ R ∗ is another 125 times slower.Figure 2 visualizes the direction fields of T ∗ ◦ R ∗ , . . . , T ∗ ◦ R ∗ on their respective manifolds M ∗ ,. . . , M ∗ along with their respective critical manifolds M ∗ , . . . , M ∗ , where M ∗ can be derived from M ∗ by additionally equating the vector field of T ∗ ◦ R ∗ to zero: M ∗ = (cid:2) y = , y = , y = , y = (cid:3) . This list M ∗ does not explicitly occur in the output. However, its preimage M is constructed inAlgorithm 5 and justifies the presence of ( M , T , R ) in the output there. The total computation timewas 0.906 s.This multiple time scale reduction of the bird flu model emphasizes a cascade of successive relaxationsof different model variables. First, the population of susceptible birds relaxes. As explained in theintroduction, by relaxation we mean that these variables reach quasi-steady state values. This relaxationis illustrated in Fig. 2(b). Then, the population of infected birds relaxes as shown in Fig. 2(c). Finally,the populations of susceptible and infected humans relax to a stable steady state as shown in Fig. 2(d),following a reduced dynamics described by T ∗ . 23igure 2: Critical manifolds and direction fields of our reductions of BioModel 716. (a) The sur-face is the critical manifold M ∗ ⊆ M ∗ = U projected from R into real ( y , y , y )-space.The line located at ( y , y ) ≈ (1 . , .
8) is the critical submanifold M ∗ ⊆ M ∗ . Thedot located at ( y , y , y ) ≈ (1 . , . , .
2) is the critical submanifold M ∗ ⊆ M ∗ .Both M ∗ and M ∗ extend to ±∞ in both y and y direction, and M ∗ is located near(1 . , . , . , . (b) The direction field of T ∗ ◦ R ∗ on M ∗ = U projected from R into real ( y , y )-space. The curve is the critical submanifold M ∗ ⊆ M ∗ . (c) The direc-tion field of T ∗ ◦ R ∗ on M ∗ projected from R into real ( y , y )-space. The line is the criticalsubmanifold M ∗ ⊆ M ∗ . The system here is slower than the one in (b) by a factor of 125. (d) The direction field of T ∗ ◦ R ∗ on M ∗ projected from R into real ( y , y )-space. The dot isthe critical submanifold M ∗ ⊆ M ∗ . The system here is slower than the one in (c) by anotherfactor of 125. 24 .2. Caspase Activation Pathway BIOMD0000000102 is a quantitative kinetic model that examines the intrinsic pathway of caspase acti-vation that is essential for apoptosis induction by various stimuli including cytotoxic stress [36]. Speciesconcentrations over time are mapped as follows:Species variable Differential variable Species A y APAF-1 C9 y Caspase 9
C9X y Caspase 9-XIAP complex X y XIAP
AC9X y APAF-1-Caspase 9-XIAP complex
AC9 y APAF-1-Caspase 9 complex C3 y Caspase 3
C3_star y Caspase 3 cleaved
C3_starX y Caspase 3 cleaved-XIAP complex
C9_starX y Caspase 9 cleaved-XIAP complex
C9_star y Caspase 9 cleaved
AC9_star y APAF-1-Caspase 9 cleaved complex
AC9_starX y APAF-1-Caspase 9 cleaved-XIAP complexThe input system is given by S = (cid:2) dd t y = − y y − y y − y y − y y − y + y + y + y + y + , dd t y = − y y − y y − y y − y + y + y + , dd t y = − y y + y y − y + y , dd t y = − y y + y − y y − y y − y y − y y − y + y + y + y + y + , dd t y = y y + y y − y , dd t y = y y − y y + y − y y − y , dd t y = − y y − y y − y y − y y − y + , dd t y = y y − y y + y y + y y + y y − y + y , dd t y = y y − y , dd t y = − y y + y y − y + y , dd t y = − y y + y y − y y + y − y + y , dd t y = y y − y y + y y − y + y , dd t y = y y + y y − y (cid:3) . We choose ε ∗ = , p = 1 and select d = ( − , , , , , , − , − , − , − , − , ,
0) from the tropicalequilibration. Our back-transformed reduced systems read as follows: M ∗ = (cid:2) (cid:3) , T ∗ = (cid:2) dd t y = 1 · (cid:0) − y y + (cid:1)(cid:3) ,R ∗ = (cid:2) dd t y = 0 , dd t y = 0 , dd t y = 0 , dd t y = 0 , dd t y = 0 , dd t y = 0 , dd t y = 0 , dd t y = 0 , dd t y = 0 , dd t y = 0 , dd t y = 0 , dd t y = 0 (cid:3) ,M ∗ = (cid:2) y y = (cid:3) , T ∗ = (cid:2) dd t y = · (cid:0) y y − y (cid:1) , d t y = · (cid:0) y y − y (cid:1) , dd t y = · (cid:0) y y − y (cid:1) , dd t y = · (cid:0) y y − y (cid:1)(cid:3) ,R ∗ = (cid:2) dd t y = 0 , dd t y = 0 , dd t y = 0 , dd t y = 0 , dd t y = 0 , dd t y = 0 , dd t y = 0 , dd t y = 0 (cid:3) ,M ∗ = (cid:2) y y = , T ∗ = (cid:2) dd t y = · (cid:0) − y y + (cid:1)(cid:3) ,y y − y = 0 , R ∗ = (cid:2) dd t y = 0 , dd t y = 0 , dd t y = 0 , dd t y = 0 , y y − y = 0 , dd t y = 0 , dd t y = 0 , dd t y = 0 (cid:3) . y y − y = 0 ,y y − y = 0 (cid:3) , The total computation time was 8.547 s, of which Algorithm 4 took 6.188 s.The multiple time scale reduction of the caspase activation model emphasizes a cascade of successiverelaxations. First, the inhibitor of apoptosis XIAP binds rapidly to the cleaved caspase. Then, the fourAPAF complexes are formed. Finally, the Caspase 9 is recruited to the apoptosome. β Pathway
BIOMD0000000101 , is a simple representation of the TGF- β signaling pathway that plays a central rolein tissue homeostasis and morphogenesis, as well as in numerous diseases such as fibrosis and cancer[60]. Concentrations over time of species RI (receptor 1), RII (receptor 2), lRIRII (ligand receptorcomplex-plasma membrane), lRIRII_endo (ligand receptor complex-endosome),
RI_endo (receptor 1endosome), and
RII_endo (receptor 2 endosome), are mapped to differential variables y , y , y , y , y , and y , respectively. The original BIOMD0000000101 has a change of ligand concentration at time t = 2500. For our computation here, we ignore this discrete event. The input system is given by S = (cid:2) dd t y = − y y − y + y + y + 8 , dd t y = − y y − y + y + y + 4 , dd t y = y y − y , dd t y = y − y , dd t y = y − y , dd t y = y − y (cid:3) . We choose ε ∗ = , p = 1, and select d = (0 , − , − , − , − , −
5) from the tropical equilibrium. Ourback-transformed reduced systems read as follows: M ∗ = (cid:2) (cid:3) , T ∗ = (cid:2) dd t y = 5 · (cid:0) − y y + (cid:1)(cid:3) ,R ∗ = (cid:2) dd t y = 0 , dd t y = 0 , dd t y = 0 , dd t y = 0 , dd t y = 0 (cid:3) ,M ∗ = (cid:2) y y = 800 (cid:3) , T ∗ = (cid:2) dd t y = 1 · (cid:0) − y + 8 (cid:1)(cid:3) ,R ∗ = (cid:2) dd t y = 0 , dd t y = 0 , dd t y = 0 , dd t y = 0 (cid:3) ,M ∗ = (cid:2) y y = 800 , T ∗ = (cid:2) dd t y = · (cid:0) − y + y (cid:1)(cid:3) ,y = (cid:3) , R ∗ = (cid:2) dd t y = 0 , dd t y = 0 , dd t y = 0 (cid:3) . The total computation time was 0.906 s. 26he multiple time scale reduction of the TGF- β model emphasizes a cascade of successive relaxationsof concentrations of different species. First, the concentration of receptor 1 relaxes rapidly. Then followsthe membrane complex, and, even slower, the relaxation of receptor 2. BIOMD0000000709 examines bird-to-human transmission of different strains of avian influenza A viruses,such as H5N1 and H7N9 [38]. Species concentrations over time of
S_a (susceptible avian),
I_a (infectedavian),
S_h (susceptible human),
I_h (infected human), and
R_h (recovered human) are mapped todifferential variables y , y , y , y , and y , respectively. The input system is given by S = (cid:2) dd t y = − y + y − y y − y , dd t y = y y − y , dd t y = − y y − y + 30 , dd t y = y y − y , dd t y = y − y (cid:3) . We choose ε ∗ = , p = 1, and select d = ( − , , − , , −
2) from the tropical equilibration. Ourback-transformed reduced systems read as follows: M ∗ = (cid:2) (cid:3) , T ∗ = (cid:2) dd t y = 1 · (cid:0) − y + y (cid:1)(cid:3) ,R ∗ = (cid:2) dd t y = 0 , dd t y = 0 , dd t y = 0 , dd t y = 0 (cid:3) ,M ∗ = (cid:2) y − y = 0 (cid:3) , T ∗ = (cid:2) dd t y = · (cid:0) y y − y (cid:1)(cid:3) ,R ∗ = (cid:2) dd t y = 0 , dd t y = 0 , dd t y = 0 (cid:3) . The total computation time was 0.578 s.The multiple time scale reduction of this avian influenza model emphasizes a cascade of successiverelaxations of different model variables. First, the susceptible bird population relaxes rapidly. Thereduced equation T and manifold M suggest that the bird population dynamics is of the Allee typeand evolves toward the stable extinct state. It follows the relaxation of infected human population thatalso evolves toward the extinct state, the end of the epidemics.
7. Concluding Remarks
We provided a symbolic method for automatic model reduction of biological networks described byordinary differential equations with multiple time scales. This method is applicable to systems withtwo time scales or more, superseding traditional slow-fast reduction methods that can cope with onlytwo time scales. We also proposed, for the first time, the automatic verification of the hyperbolicityconditions required for the validity of the reduction. Our theoretical framework is accompanied by rigor-ous algorithms and prototypical implementations, which we successfully applied to real-world problemsfrom the BioModels database [34].We would like to list some open points and possible extensions of our research here. Our reductionalgorithm is based on a fixed scaling (8) leading to a fixed ordering of the time scales of differentvariables. In our reduction scheme, different variables relax hierarchically, firstly the fastest ones, thenthe second fastest, and lastly the slowest ones, which justifies our geometric picture of nested invariantmanifolds. However, there are situations, e.g. in models of relaxation oscillations, when the orderingof time scales changes with time: variables that were fast can become slow at a later time, and viceversa. In order to cope with such situations, one would like to use different scalings for different timesegments. One attempt to implement such a procedure has been provided in [54].Although our proposed method identifies the full hierarchy of time scales, the subsequent reductionmay stop early in this hierarchy when hyperbolic attractivity is not satisfied at some stage. One27ossible reason is the presence of conservation laws, also known as first integrals, at the given reductionstage. Such conservation laws necessarily force an eigenvalue zero for the Jacobian. A theorem bySchneider and Wilhelm [52] can be employed to reduce such a setting to the hyperbolically attractivecase. As for the behavior of first integrals when proceeding to the reduced system, see the discussionof the non-standard case in [33] for two time scales; an extension to multiple time scales should bestraightforward. Work in progress is concerned with the introduction of novel slow variables, one forevery independent conservation law of the fast subsystem, applying this to networks with multiple timescales, and approximate linear and polynomial conservation laws.More generally, it is of interest to consider cases when hyperbolic attractivity fails but hyperbolicitystill holds: In such cases, Cardin and Teixeira show there still exist invariant manifolds [10]. Testingfor hyperbolicity is more involved than testing for hyperbolic attractivity, but in theory it is wellunderstood, and there exists an algorithmic approach due to Routh [23]. In the case of hyperbolicity,but not attractivity, the ensuing global dynamics may be quite interesting; for instance slow-fast cyclesmay appear.Concerning differentiability requirements, we checked for smoothness of the full system in Sect. 3.3.However, Fenichel’s results, and in principle also those by Cardin and Teixeira, require only sufficientfinite differentiability. Therefore, given a differential equation system and a scaling, invariant manifoldsand corresponding reduced systems exist for C p functions with fixed p < ∞ . Going through the detailswill involve intricate analysis that is left to future work.In the introduction we sketched a Michaelis–Menten system abstracting from the known numericalvalues for the reaction rate constants k , k − , k . It would be indeed interesting to work on suchparametric data. In the presence of parameters, one would consider effective quantifier eliminationover real closed fields [14, 63, 31, 55] as a generalization of SMT solving. Robust implementationsare freely available [8, 18] and well supported. They have been successfully applied to problems inchemical reaction network theory during the past decade [56, 62, 7, 19]. Such a generalization is notquite straightforward. With the tropical scaling in Sect. 2.2, Algorithm 2 would introduce logarithmsof polynomials in the parametric coefficients, which is not compatible with the logical framework usedhere. Similar tropicalization methods, which are unfortunately not compatible with our abstract viewon scaling in Sect. 2.1, require only logarithms of individual parametric coefficients [51]. Such a morespecial form would allow the use of abstraction in the logic engine.From a point of view of user-oriented software, it would be most desirable to develop automaticstrategies for determining good default values for ε ∗ and for choices of d from the tropical equilibrationin Algorithm 3. Acknowledgments
This work has been supported by the interdisciplinary bilateral project ANR-17-CE40-0036/DFG-391322026 SYMBIONT [5, 6].
References [1] Clark Barrett, Christopher L. Conway, Morgan Deters, Liana Hadarean, Dejan Jovanović, TimKing, Andrew Reynolds, and Cesare Tinelli. CVC4. In Ganesh Gopalakrishnan and Shaz Qadeer,editors,
Proc. CAV 2011 , volume 6806 of
LNCS , pages 171–177. Springer, July 2011. doi:10.1007/978-3-642-22110-1_14 .[2] Clark Barrett, Pascal Fontaine, and Cesare Tinelli. The SMT-LIB standard: Version 2.6. Technicalreport, Department of Computer Science, The University of Iowa, 2017.[3] Thomas Becker, Volker Weispfenning, and Heinz Kredel.
Gröbner Bases, a Computational Ap-proach to Commutative Algebra , volume 141 of
Graduate Texts in Mathematics . Springer, 1993. doi:10.1007/978-1-4612-0913-3 . 284] Tristram Bogart, Anders Nedergaard Jensen, David Speyer, Bernd Sturmfels, and Rekha R.Thomas. Computing tropical varieties.
J. Symb. Comput. , 42(1–2):54–73, January–February 2007. doi:10.1016/j.jsc.2006.02.004 .[5] François Boulier, François Fages, Ovidiu Radulescu, Satya S. Samal, Aandreas Schuppert,Werner M. Seiler, Thomas Sturm, Sebastian Walcher, and Anderas Weber. The SYMBIONTproject: Symbolic methods for biological networks.
F1000Research , 7(1341), August 2018. doi:10.7490/f1000research.1115995.1 .[6] François Boulier, François Fages, Ovidiu Radulescu, Satya S. Samal, Andreas Schuppert,Werner M. Seiler, Thomas Sturm, Sebastian Walcher, and Andreas Weber. The SYMBIONTproject: Symbolic methods for biological networks.
ACM Communications in Computer Algebra ,52(3):67–70, December 2018. doi:10.1145/3313880.3313885 .[7] Russell Bradford, James H. Davenport, Matthew England, Hassan Errami, Vladimir Gerdt, DimaGrigoriev, Charles Hoyt, Marek Košta, Ovidiu Radulescu, Thomas Sturm, et al. A case study onthe parametric occurrence of multiple steady states. In Michael Burr, editor,
Proc. ISSAC 2017 ,pages 45–52. ACM, 2017. doi:10.1145/3087604.3087622 .[8] Christopher W. Brown. QEPCAD B: A program for computing with semi-algebraic sets usingCADs.
ACM SIGSAM Bulletin , 37(4):97–108, December 2003. doi:10.1145/968708.968710 .[9] Bruno Buchberger.
Ein Algorithmus zum Auffinden der Basiselemente des Restklassenringes nacheinem nulldimensionalen Polynomideal . Doctoral dissertation, Mathematical Institute, Universityof Innsbruck, Innsbruck, Austria, 1965.[10] Pedro Toniol Cardin and Marco Antonio Teixeira. Fenichel theory for multiple time scale singularperturbation problems.
SIAM J. Appl. Dyn. Syst. , 16(3):1425–1452, January 2017. doi:10.1137/16m1067202 .[11] Pedro Toniol Cardin and Marco Antonio Teixeira. Corrigendum: Fenichel theory for multipletime scale singular perturbation problems.
SIAM J. Appl. Dyn. Syst. , 18(2):1223, 2019. doi:10.1137/19M1241660 .[12] Alessandro Cimatti, Alberto Griggio, Bastiaan Schaafsma, and Roberto Sebastiani. The MathSAT5SMT solver. In Nir Piterman and Scott A. Smolka, editors,
Proc. TACAS 2013 , volume 7795 of
LNCS , pages 93–107. Springer, 2013. doi:10.1007/978-3-642-36742-7_7 .[13] George E. Collins. Quantifier elimination for the elementary theory of real closed fields by cylindri-cal algebraic decomposition. In H. Brakhage, editor,
Automata Theory and Formal Languages. 2ndGI Conference , volume 33 of
LNCS , pages 134–183. Springer, 1975. doi:10.1007/3-540-07407-4_17 .[14] George E. Collins and Hoon Hong. Partial cylindrical algebraic decomposition for quantifierelimination.
J. Symb. Comput. , 12(3):299–328, September 1991. doi:10.1016/S0747-7171(08)80152-6 .[15] Florian Corzilius, Gereon Kremer, Sebastian Junges, Stefan Schupp, and Erika Ábrahám. SMT-RAT: An open source C++ toolbox for strategic and parallel SMT solving. In Marijn Heule andSean Weaver, editors,
Proc. SAT 2015 , volume 9340 of
LNCS , pages 360–368. Springer, 2015. doi:10.1007/978-3-319-24318-4_26 .[16] Haskell B. Curry and Robert Feys.
Combinatory Logic , volume I of
Studies in Logic and theFoundations of Mathematics . North Holland Publishing Company, Amsterdam, The Netherlands,1958. 2917] Leonardo de Moura and Nikolaj Bjørner. Z3: An efficient SMT solver. In C. R. Ramakrishnanand J. Rehof, editors,
Proc. TACAS 2008 , volume 4963 of
LNCS , pages 337–340. Springer, 2008. doi:10.1007/978-3-540-78800-3_24 .[18] Andreas Dolzmann and Thomas Sturm. Redlog: Computer algebra meets computer logic.
ACMSIGSAM Bulletin , 31(2):2–9, June 1997. doi:10.1145/261320.261324 .[19] Matthew England, Hassan Errami, Dima Grigoriev, Ovidiu Radulescu, Thomas Sturm, andAndreas Weber. Symbolic versus numerical computation and visualization of parameter re-gions for multistationarity of biological networks. In Gerdt V., Koepf W., Seiler W., andVorozhtsov E., editors,
Proc. CASC 2017 , volume 10490 of
LNCS , pages 93–108. Springer, 2017. doi:10.1007/978-3-319-66320-3_8 .[20] Martin Feinberg.
Foundations of Chemical Reaction Network Theory , volume 202 of
Applied Math-ematical Sciences . Springer, 2019. doi:10.1007/978-3-030-03858-8 .[21] Neil Fenichel. Geometric singular perturbation theory for ordinary differential equations.
J. Differ.Equations , 31(1):53–98, January 1979. doi:10.1016/0022-0396(79)90152-9 .[22] Stephen Forrest. Integration of SMT-LIB support into Maple. In
Proc. Satisfiability Checking andSymbolic Computation 2017 , volume 1974 of
CEUR Workshop Proceedings . CEUR-WS, Kaiser-slautern, Germany, July 2017. doi:10.1145/3055282.3055285 .[23] Feliks R. Gantmacher.
The Theory of Matrices. Vol. 2. Transl. from the Russian by K. A. Hirsch.
Providence, RI: AMS Chelsea Publishing, reprint of the 1959 translation edition, 1998.[24] Marco Gario and Andrea Micheli. PySMT: A solver-agnostic library for fast prototyping of SMT-based algorithms. In
SMT Workshop 2015. 13th International Workshop on Satisfiability ModuloTheories, Affiliated With the 27th International Conference on Computer Aided Verification , SanFrancisco, CA, July 2015.[25] Naama Geva-Zatorsky, Nitzan Rosenfeld, Shalev Itzkovitz, Ron Milo, Alex Sigal, Erez Dekel, TaliaYarnitzky, Yuvalal Liron, Paz Polak, Galit Lahav, and Uri Alon. Oscillations and variability inthe p53 system.
Mol. Syst. Biol. , 2(1):2006–0033, June 2006. doi:10.1038/msb4100068 .[26] Alexandra Goeke, Sebastian Walcher, and Eva Zerz. Determining “small parameters” for quasi-steady state.
J. Differ. Equations , 259(3):1149–1180, August 2015. doi:10.1016/j.jde.2015.02.038 .[27] Dima Grigoriev. Complexity of deciding Tarski algebra.
J. Symb. Comput. , 5(1–2):65–108,February–April 1988. doi:/10.1016/S0747-7171(88)80006-3 .[28] F. G. Heineken, H. M. Tsuchiya, and R. Aris. On the mathematical status of the pseudo-steadystate hypothesis of biochemical kinetics.
Math. Biosci. , 1(1):95–113, March 1967. doi:10.1016/0025-5564(67)90029-6 .[29] Frank Hoppensteadt. On systems of ordinary differential equations with several parametersmultiplying the derivatives.
J. Differ. Equations , 5(1):106–116, January 1969. doi:10.1016/0022-0396(69)90106-5 .[30] Adolf Hurwitz. Ueber die Bedingungen, unter welchen eine Gleichung nur Wurzeln mit negativenreellen Theilen besitzt.
Math. Ann. , 46:273–284, June 1895. doi:10.1007/BF01446812 .[31] Marek Košta.
New Concepts for Real Quantifier Elimination by Virtual Substitution . Doctoraldissertation, Saarland University, Germany, December 2016. doi:10.22028/D291-26679 .[32] Niclas Kruff and Sebastian Walcher. Coordinate-independent singular perturbation reduction forsystems with three time scales.
Math. Biosci. Eng. , 16(5):5062–5091, June 2019. doi:10.3934/mbe.2019255 . 3033] Christian Lax and Sebastian Walcher. Singular perturbations and scaling.
Discrete Cont. Dyn.-B ,25(1):1–29, January 2020. doi:10.3934/dcdsb.2019170 .[34] Nicolas Le Novere, Benjamin Bornstein, Alexander Broicher, Melanie Courtot, Marco Donizelli,Harish Dharuri, Lu Li, Herbert Sauro, Maria Schilstra, Bruce Shapiro, et al. BioModels database:A free, centralized database of curated, published, quantitative kinetic models of biochemical andcellular systems.
Nucleic acids res. , 34(suppl_1):D689–D691, January 2006. doi:10.1093/nar/gkj092 .[35] Hanl Lee and Angelyn Lao. Transmission dynamics and control strategies assessment of avianinfluenza A (H5N6) in the Philippines.
Infectious Disease Modelling , 3:35–59, 2018. doi:10.1016/j.idm.2018.03.004 .[36] Stefan Legewie, Nils Blüthgen, and Hanspeter Herzel. Mathematical modeling identifies inhibitorsof apoptosis as mediators of positive feedback and bistability.
PLoS Comput. Biol. , 2(9):e120,September 2006. doi:10.1371/journal.pcbi.0020120 .[37] Grigori L. Litvinov. Maslov dequantization, idempotent and tropical mathematics: A brief intro-duction.
J. Math. Sci. , 140(3):426–444, January 2007. doi:10.1007/s10958-007-0450-5 .[38] Sanhong Liu, Shigui Ruan, and Xinan Zhang. Nonlinear dynamics of avian influenza epidemicmodels.
Math. Biosci. , 283:118–135, January 2017. doi:10.1016/j.mbs.2016.11.014 .[39] Christoph Lüders. Computing tropical prevarieties with satisfiability modulo theories (SMT)solvers.
CoRR , abs/2004.07058, April 2020. arXiv:2004.07058 .[40] Aaron Meurer, Christopher P. Smith, Mateusz Paprocki, Ondřej Čertík, Sergey B. Kirpichev,Matthew Rocklin, Amit Kumar, Sergiu Ivanov, Jason K. Moore, Sartaj Singh, Thilina Rathnayake,Sean Vig, Brian E. Granger, Richard P. Muller, Francesco Bonazzi, Harsh Gupta, Shivam Vats,Fredrik Johansson, Fabian Pedregosa, Matthew J. Curry, Andy R. Terrel, Štěpán Roučka, AshutoshSaboo, Isuru Fernando, Sumith Kulal, Robert Cimrman, and Anthony Scopatz. SymPy: Symboliccomputing in Python.
PeerJ Computer Science , 3:e103, January 2017. doi:10.7717/peerj-cs.103 .[41] Robert Nieuwenhuis, Albert Oliveras, and Cesare Tinelli. Solving SAT and SAT modulo theories:From an abstract Davis–Putnam–Logemann–Loveland procedure to DPLL(T).
J. ACM , 53(6):937–977, November 2006. doi:10.1145/1217856.1217859 .[42] Kaspar Nipp. An algorithmic approach for solving singularly perturbed initial value problems. InUrs Kirchgraber and Hans-Otto Walther, editors,
Dynamics Reported , volume 1, pages 173–263.John Wiley & Sons and B. G. Teubner, 1988. doi:10.1007/978-3-322-96656-8_4 .[43] Vincent Noel, Dima Grigoriev, Sergei Vakulenko, and Ovidiu Radulescu. Tropical geometries anddynamics of biochemical networks application to hybrid cell cycle models. In Jérôme Feret andAndre Levchenko, editors,
Proc. SASB 2011 , volume 284 of
ENTCS , pages 75–91. Elsevier, 2012. doi:10.1016/j.entcs.2012.05.016 .[44] Vincent Noel, Dima Grigoriev, Sergei Vakulenko, and Ovidiu Radulescu. Tropicalization andtropical equilibration of chemical reactions. In G. L. Litvinov and S. N. Sergeev, editors,
Tropicaland Idempotent Mathematics and Applications , volume 616 of
Contemporary Mathematics , pages261–277. AMS, 2014. doi:10.1090/conm/616/12316 .[45] Lena Noethen and Sebastian Walcher. Quasi-steady state and nearly invariant sets.
SIAM J. Appl.Math. , 70(4):1341–1363, 2009. doi:10.1137/090758180 .[46] Ovidiu Radulescu, Alexander N. Gorban, Andrei Zinovyev, and Vincent Noel. Reduction of dy-namical biochemical reactions networks in computational biology.
Front. Genet. , 3:131, July 2012. doi:10.3389/fgene.2012.00131 . 3147] Ovidiu Radulescu, Sergei Vakulenko, and Dima Grigoriev. Model reduction of biochemical reactionsnetworks by tropical analysis methods.
Math. Model. Nat. Pheno. , 10(3):124–138, June 2015. doi:10.1051/mmnp/201510310 .[48] Dennis Reddyhoff, John Ward, Dominic Williams, Sophie Regan, and Steven Webb. Timescaleanalysis of a mathematical model of acetaminophen metabolism and toxicity.
J. Theor. Biol ,386:132–146, December 2015. doi:10.1016/j.jtbi.2015.08.021 .[49] Shigui Ruan. Modeling the transmission dynamics and control of rabies in China.
Math. Biosci. ,286:65–93, April 2017. doi:10.1016/j.mbs.2017.02.005 .[50] Satya S. Samal, Dima Grigoriev, Holger Fröhlich, and Ovidiu Radulescu. Analysis of reactionnetwork systems using tropical geometry. In V. Gerdt, W. Koepf, W. Seiler, and E. Vorozhtsov,editors,
Proc. CASC 2015 , volume 9301 of
LNCS , pages 424–439. Springer, 2015. doi:10.1007/978-3-319-24021-3_31 .[51] Satya S. Samal, Dima Grigoriev, Holger Fröhlich, Andreas Weber, and Ovidiu Radulescu. Ageometric method for model reduction of biochemical networks with polynomial rate functions.
B.Math. Biol. , 77(12):2180–2211, October 2015. doi:10.1007/s11538-015-0118-0 .[52] Klaus R. Schneider and Thomas Wilhelm. Model reduction by extended quasi-steady-state ap-proximation.
J. Math. Biol. , 40(5):443–450, May 2000. doi:10.1007/s002850000026 .[53] Lee A. Segel and Marshall Slemrod. The quasi-steady-state assumption: A case study in pertur-bation.
SIAM Rev. , 31(3):446–477, September 1989. doi:10.1137/1031091 .[54] Jasha Sommer-Simpson, John Reinitz, Leonid Fridlyand, Louis Philipson, and Ovidiu Radulescu.Hybrid reductions of computational models of ion channels coupled to cellular biochemistry. InEzio Bartocci, Pietro Lio, and Nicola Paoletti, editors,
Proc. CMSB 2016 , volume 9859 of
LNCS ,page 273. Springer, 2016. doi:10.1007/978-3-319-45177-0_17 .[55] Thomas Sturm. A survey of some methods for real quantifier elimination, decision, and satisfiabilityand their applications.
Math. Comput. Sci. , 11(3–4):483–502, December 2017. doi:10.1007/s11786-017-0319-z .[56] Thomas Sturm, Andreas Weber, Essam O. Abdel-Rahman, and M’hammed El Kahoui. Investigat-ing algebraic and logical algorithms to solve Hopf bifurcation problems in algebraic biology.
Math.Comput. Sci. , 2(3):493–515, March 2009. doi:10.1007/s11786-008-0067-1 .[57] Alfred Tarski. A decision method for elementary algebra and geometry. Prepared for publicationby J. C. C. McKinsey. RAND Report R109, August 1, 1948, Revised May 1951, Second Edition,RAND, Santa Monica, CA, 1957.[58] Andrei Nikolaevich Tikhonov. Systems of differential equations containing small parameters in thederivatives.
Mat. Sb. (N. S.) , 73(3):575–586, 1952.[59] Mauro Valorani and Samuel Paolucci. The G-scheme: A framework for multi-scale adaptive modelreduction.
J. Comput. Phys. , 228(13):4665–4701, July 2009. doi:10.1016/j.jcp.2009.03.011 .[60] Jose M. G. Vilar, Ronald Jansen, and Chris Sander. Signal processing in the TGF- β superfamilyligand-receptor network. PLoS Comput. Biol. , 2(1):e3, January 2006. doi:10.1371/journal.pcbi.0020003 .[61] Oleg Viro. Dequantization of real algebraic geometry on logarithmic paper. In Carles Casacuberta,Rosa Maria Miró-Roig, Joan Verdera, and Sebastià Xambó-Descamps, editors,
European Congressof Mathematics , volume 201 of
Progress in Mathematics , pages 135–146. Springer, 2001. doi:10.1007/978-3-0348-8268-2_8 . 3262] Andreas Weber, Thomas Sturm, and Essam O. Abdel-Rahman. Algorithmic global criteria for ex-cluding oscillations.
B. Math. Biol. , 73(4):899–916, April 2011. doi:10.1007/s11538-010-9618-0 .[63] Volker Weispfenning. Quantifier elimination for real algebra—the quadratic case and beyond.
Appl.Algebr. Eng. Comm. , 8(2):85–101, February 1997. doi:10.1007/s002000050055 .[64] Dominik Wodarz and Dean H. Hamer. Infection dynamics in HIV-specific CD4 T cells: Does aCD4 T cell boost benefit the host or the virus?
Math. Biosci. , 209(1):14–29, September 2007. doi:10.1016/j.mbs.2007.01.007 . A. Illustration of Some Border Cases With Our Algorithms
A.1. Failure of Tropicalization With an Unbalanced Monomial
BIOMD0000000609 describes the metabolism and the related hepatotoxicity of acetaminophen, a painkiller [48]. The species concentrations over time of
Sulphate__PAPS , GSH , NAPQI , Paracetamol_APAP ,and
Protein_adducts are mapped to differential variables y , y , y , y , and y , respectively. Theinput system is given by S = (cid:2) dd t y = − y y − y + , dd t y = − y y − y + , dd t y = − y y − y + y , dd t y = − y y + y − y , dd t y = 110 y (cid:3) . Since there is only one monomial on the right hand side of the equation for ˙ y , equilibration is impossible.This causes Algorithm 4 to return in l.23 a disjunctive normal form Π equivalent to “false”, whichdescribes the empty set. Hence Algorithm 3 returns ⊥ , and Algorithm 1 returns the empty list. Thetotal computation time was 0.006 s. A.2. Failure of Hyperbolic Attractivity due to an Empty Manifold in the FirstOrthant
BIOMD0000000726 examines the transmission dynamics of rabies between dogs and humans [49]. Themapping of the model variables over time and our differential variables is as follows:Species Differential variable Description
S_d y susceptible dogs E_d y exposed dogs I_d y infectious dogs R_d y recovered dogs S_h y susceptible humans E_h y exposed humans I_h y infectious humans R_h y recovered humansThe input system is given by S = (cid:2) dd t y = − y y − y + y + y + 3000000 , dd t y = y y − y , dd t y = y − y , d t y = y + y − y , dd t y = − y y − y + y + y + 15400000 , dd t y = y y − y , dd t y = y − y , dd t y = y − y (cid:3) . We choose parameters ε ∗ = , p = 1, and d = ( − , − , − , − , − , − , − , − M = [ − x x + 4608 x = 0 , x x − x = 0 , x x − x = 0] , i.e., the corresponding manifold M is empty. Consequently, Algorithms 5, 8, and 9 return empty lists.The total computation time was 0.921 s. A.3. Failure of Hyperbolic Attractivity in the First Step
BIOMD0000000156 examines the dynamics of a negative feedback loop between the tumor suppressorprotein p53 and the oncogene protein Mdm2 in human cells [25]. The species concentrations over timefor x (p53), y (Mdm2), and y0 (precursor Mdm2) are mapped to differential variables y , y , and y ,respectively. The input system is given by S = (cid:2) dd t y = − y y + 2 y , dd t y = − y + y , dd t y = y − y (cid:3) . We choose parameters ε ∗ = , p = 1, and point d = (2 , , M , defined by M = [ − x x + x = 0], is not empty in the first orthant, as has been ensured in l.1. Obviously, thesimplified and back-translated systems are empty lists as well. The total computation time was 0.453 s. A.4. Failure in SMT Solving for a Reduced System With Fractional Exponents
BIOMD0000000663 illustrates how CD4 T-cells can influence the spread of the HIV infection [64]. Speciesconcentrations over time for x (x_Tcell_infected), y (y_Tcell_uninfected), and v (v_free_virus) aremapped to variables y , y , and y , respectively. The input system is given by S = (cid:2) dd t y = − y y − y y y + y y − y , dd t y = − y y y + y y − y y + y y − y , dd t y = y − y (cid:3) . We choose ε ∗ = , p = 5, and d = (1 , , p = 5 causes fractional exponents in thescaled and truncated system, viz. T = (cid:2) dd τ x = δ · (cid:0) √ x x − √ x (cid:1)(cid:3) , T = (cid:2) dd τ x = δ · (cid:0) √ x x − √ x (cid:1)(cid:3) . However, the relevant SMT logics
QF_LRA and
QF_NRA do not accept fractional exponents. Recall fromSect. 3.2 that in such cases, we catch the corresponding error from the SMT solver, convert to floats,and restart.We then get into the special case that ‘‘