Algorithmic Reduction of Biological Networks With Multiple Time Scales

We present a symbolic algorithmic approach that allows to compute invariant manifolds and corresponding reduced systems for differential equations modeling biological networks which comprise chemical reaction networks for cellular biochemistry, and compartmental models for pharmacology, epidemiology and ecology. Multiple time scales of a given network are obtained by scaling, based on tropical geometry. Our reduction is mathematically justified within a singular perturbation setting using a recent result by Cardin and Teixeira. The existence of invariant manifolds is subject to hyperbolicity conditions, which we test algorithmically using Hurwitz criteria. We finally obtain a sequence of nested invariant manifolds and respective reduced systems on those manifolds. Our theoretical results are generally accompanied by rigorous algorithmic descriptions suitable for direct implementation based on existing off-the-shelf software systems, specifically symbolic computation libraries and Satisfiability Modulo Theories solvers. We present computational examples taken from the well-known BioModels database using our own prototypical implementations.


Introduction
Biological network models describing elements in interaction are used in many areas of biology and medicine. Chemical reaction networks are used as models of cellular biochemistry, including gene regulatory networks, metabolic networks and signaling networks. In epidemiology and ecology, compartmental models can be described as networks of interactions between compartments. Both in chemical reaction networks and in compartmental models the probability that two elements interact is assumed proportional to their abundances. This property, called mass action law in biochemistry, leads to polynomial differential equations in the kinetics.
For differential equations that describe the development of such networks over time a crucial question is concerned with reduction of dimension. We illustrate such a reduction and the steps involved for 1 arXiv:2010.10129v1 [q-bio.MN] 20 Oct 2020 the classical Michaelis-Menten system, an archetype of enzymatic reactions. The differential equations for the concentrations of relevant chemical species, which are substrate and enzyme-substrate complex, have the formẏ 1 = −εk 1 y 1 + (k 1 y 1 + k −1 )y 2 y 2 = εk 1 y 1 − (k 1 y 1 + k −1 + k 2 )y 2 , involving a small parameter ε that represents the ratio of the total concentration of the enzyme to the concentration of the substrate. The fact that this ratio is small is an assumption of the model that has to be verified in applications. In a first step toward reduction, a scaling transformation y 1 = x 1 and y 2 = εx 2 yieldsẋ In a second step, one uses singular perturbation theory to obtain the famous Michaelis-Menten equation. It consists of two components: First, we obtain a one dimensional invariant manifold given approximately by the quasi-steady state condition k 1 x 1 − (k 1 x 1 + k −1 + k 2 )x 2 = 0. This considers the fast variable x 2 to be at the steady state and lowers dimension from two to one. Second, we obtain a reduced system for the slow variable:ẋ With our example, we paraphrased the approach in a seminal paper by Heineken et al. [28], which was the first one to rigorously discuss quasi-steady state from the perspective of singular perturbation theory. Realistic network models may have many species and differential equations. Considerable effort has been put into model order reduction, i.e., finding approximate models with a smaller number of species 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 ε → 0 is central in singular perturbation theory. It is rather obvious that arbitrary scaling transformations are unlikely to provide useful information about a given system. Successful scalings, in contrast, are typically related to the existence of nontrivial invariant manifolds. Applications of scaling rely on the observation that, loosely speaking, any result that holds asymptotically for ε → 0 remains valid for sufficiently small positive ε * , provided some technical conditions are satisfied. To determine scalings of polynomial or rational vector fields that model biological networks, tropical equilibration methods were introduced and developed in a series of papers by Noel et al. [44], Radulescu et al. [47], Samal et al. [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 invariant manifolds and reduced equations. Other methods due to Goeke et al. [26], and recently extended to multiple time scales by Kruff and Walcher [32], determine critical parameter values and manifolds for singular perturbation reductions.
The principal purpose of the present paper is to complement scaling with an algorithmic test for the existence of invariant manifolds and the computation of those manifolds along with corresponding reduced systems of differential equations. In the asymptotic limit, methods from singular perturbation theory, principally developed by Tikhonov [58] and Fenichel [21], are available. A recent extension to multiscale systems by Cardin and Teixeira [10] turns out to be a valuable tool for the systematic computation 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 evolves toward quasi-steady state values, and all slower variables are constant. The sets of quasi-steady state conditions for relaxed variables define invariant manifolds, more precisely, they provide the lowest order approximations to the invariant manifolds. As the set of relaxed variables and thus quasi-steady state conditions increases, the respective invariant manifolds get nested so that later manifolds are contained in 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 are still missing. system of the form z k = ε a k f k (z, ε) = ε a k f k (z, 0) + ε a k,2 p k,2 + · · · + ε a k,w k p k,w k = ε a k f k (z, 0) + o(1) , 1 ≤ k ≤ m, (7) where a k , a k,j ∈ Q, 0 = a 1 < a 2 < . . . < a m , 0 < a k,j , and p k,j are multivariate polynomials in z for 1 ≤ k ≤ m and 2 ≤ j ≤ w k . Note that the case m = 1 is not excluded. By substituting δ := ε 1/q , with a sufficiently large positive integer q, one ensures that only nonnegative integer powers of δ appear: 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 in Sect. 3 that the exact number of time scales finally obtained by our overall approach can actually be smaller than m.
Given certain conditions, which will be made explicit in Theorem 1 and with its application in Sect. 3.1, we may formally truncate the right hand sides of (8) and keep only terms of lowest order in δ: In the sequel, we refer to the transformation process from (1) to (8) as scaling. Strictly speaking, this comprises 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 representing system (1) and a choice of 0 < ε * < 1 for (2). For our practical purposes, the polynomial coefficients in S as well as ε * are taken from Q. Our algorithm is furthermore parameterized with a function c mapping suitable indices to rational numbers and a constant function d yielding either a tuple D = (d 1 , . . . , d n ) of rational numbers or ⊥. The black-box functions c and d reflect the mathematical assumptions around (2) and (4) that suitable c k,J and d k exist, respectively. Suitable instantiations for the parameters c and d can be realized, e.g., using tropical geometry, which will be the topic of Sect. 2.2. It will turn out that instantiations of d can fail on the given combination of S and ε * , which is signaled by the return value ⊥ of d, and checked right away in l.1 of Algorithm 1.

Scaling via Tropical Geometry
So far, the above transformations leading to (4) are a formal exercise. No particular strategy was applied for choosing ε * ∈ (0, 1). Early model reduction studies used dimensional analysis to obtain ε * 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 approach starts 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 descriptions of all the quantities occurring in the differential equations (parameters, monomials, time scales), in a similar way to describing curves by continuously varying real parameters. This interpretation is rooted in physics where it unravels scaling laws. It is also natural for any computation with orders of magnitude. In tropical geometry, it is encountered in several places: as Litvinov-Maslov dequantization of real numbers leading to degeneration of complex algebraic varieties into tropical varieties [37,61], or in the theory of Puiseux series in relation to tropical varieties and pre-varieties [4].
This gives the following invariant: Denote S = ( m k=1 T k ⊕ P k ) σ, where (x = g) ⊕ p stands for x = g + p and is applied elementwise. Then S is equal to S up to multiplication of the differential equationẏ i = J γ i,J y J in S with a positive scalar factor 1/ε µ+di * .
Secondly, the orders D = (d 1 , d 2 , . . . , d n ) satisfy certain constraints. These constraints result heuristically from the idea of compensation of dominant monomials [43]. Slow dynamics is possible if for each dominant, i.e., much larger than the other, monomial on the right hand side of (6), there is at least one other monomial of the same order, but with opposite sign. This condition, named tropical equilibration condition [43,44,46,47,51,50], reads 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 tropical scaling. The result depends on which parameters are close and which are very different as dictated by the coarse-graining procedure, i.e., by the choice of ε * . Decreasing ε * destroys details and parameters tend to have the same order of magnitude. Increasing ε * refines details and parameters range over several orders of magnitude. For instance, using (10) and p = 1 parameters k 1 = 0.1 and k 2 = 0.01 have orders c 1 = 1 and c 2 = 2 for ε * = 1/10, but c 1 = c 2 = 1 for ε * = 1/50. This is the perspective taken in [43,44,51].
It is noteworthy that in the context of singular perturbation methods (cf. Sect. 3), which provide asymptotic 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 with tropical versions as given in Algorithm 2 and Algorithm 3, respectively.
We use ·, · to denote the standard scalar product in Q n+1 .

15:
for m := 1 to c do 16: Algorithm 2 explicitly uses, besides the parameters k and J specified for c in Algorithm 1, also the right hand sides of the input system (1) and the choice of ε * . As yet another parameter it takes the desired precision p for rounding in (10). Notice that the use of this extra information is compatible with 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 no parameters at all. In l.1 we use Algorithm 4 as a subalgorithm for tropical equilibration. One obtains a disjunctive normal form Π, which explicitly describes a set P = { p ∈ Q n | Π(p) } as a finite union of convex polyhedra, as known from tropical geometry. Every (d 1 , . . . , d n ) ∈ P satisfies (11). The satisfiability condition in l.2 tests whether P = ∅. 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. Such an exceptional situation is signaled with a return value ⊥ in l.3. In the regular case P = 0, the choice (d 1 , . . . , d n ) in l.5 is provided by the SMT solver. From a practical point of view, the disjunctive normal form 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 1 , . . . , d n ) lying in a specific convex polyhedron P ⊆ P, which technically corresponds to one conjunction in Π. Such choices are subtle and typically require human interaction. For instance, when the chain of reduced dynamical systems ends with a steady state, it is interesting to consider the polyhedron P that is closest to that steady state. Such strategies are not implemented in the current version of our algorithm.
At this stage we have obtained a scaled system as defined in Sect. 2.1, including partitioning. The focus of the next section is to utilize this scaling for analytically substantiated reductions.

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 improved by Hoppensteadt, show the convergence of the solution of system (8) to the solution of a differentialalgebraic system in which the slowest variables z m follow differential equations and the remaining fast variables 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 define a slow invariant manifold that is persistent for δ > 0 [21]. For geometrical singular perturbations, differentiability in δ is needed in system (8).
Samal et al. have noted that Tikhonov's theorem is applicable to tropically scaled systems [51]. For instance, with δ 1 = δ b2 , system (8) may be rewritten as However, this approach comes with certain limitations. To start with, it allows only two time scales. Furthermore, in case b 2 > 1, there may be differentiability issues with respect to δ 1 , and some care has to be taken when one tries to apply to (12) also Fenichel's results [21]. In this section, we are going to generalize geometrical singular perturbations, and compute invariant manifolds and reduced models for more than two time scales, introducing further δ 2 , . . . , δ . Our generalization is based on a recent paper 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 introduce the stronger notion of hyperbolic attractivity. In Sect. 3.2 we describe efficient algorithmic tests for hyperbolic attractivity. Section 3.3 gives sufficient algorithmic criteria addressing the above-mentioned differentiability issues.

Application of a Fenichel Theory for Multiple Time Scales
From now on we consider our system (8) over the positive first orthant U = (0, ∞) n ⊆ R n . A recent paper by Cardin and Teixeira [10] generalizes Fenichel's theory to provide a solid foundation to obtain more than one nontrivial invariant manifold. This allows, in particular, the reduction of multi-time scale systems such as system (8). Technically, the approach considers a multi-parameter system using time scale factors δ 1 , δ 1 δ 2 , . . . instead of increasing powers of one single δ.
We let ∈ {2, . . . , m} and define and furthermore These definitions allow us to express also all δ b k,j occurring in (8) as products of powers of δ 1 , . . . , δ −1 , with nonnegative but possibly non-integer rational exponents, via expressing each b k,j as a nonnegative rational linear combination of β 1 , . . . , β −1 . This yields Moreover, we express which is obtained by writing each b k − b as a nonnegative rational linear combination of β 1 , . . . , β −1 .
In these terms our system (8) translates to In terms of the right hand sides of (16) the application of relevant results in [10] requires that g 1 , . . . , g and η +1 g +1 , . . . , η m g m are smooth on an open neighborhood of U × [0, We are going to tacitly assume such smoothness here and address this issue from an algorithmic point of view in Sect. 3.3.
We are now ready to transform our system into time scales as follows, where possibly > 2: In time scale τ k , with 1 ≤ k ≤ , system (16) then becomes . . .
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 approximated by the elimination of higher order terms. We are now going to introduce a construction required for a justification of this approximation, which also clarifies the greatest possible choice for ≤ m above. Define F 0 = 0 and With system (8) in mind, we are going to use f k (z, 0) in favor of g k (z, 0, . . . , 0). It is easy to see that both are equal. We define furthermore The sets M k are obtained from varieties defined by the lists M k via intersection with the first orthant. Furthermore, U = M 0 ⊇ M 1 ⊇ · · · ⊇ M m establishes a chain of nested subvarieties, again intersected with the first orthant. We say that M 1 is hyperbolically attractive on M 0 , if M 1 = ∅ and for all z ∈ M 1 all eigenvalues of the Jacobian D z1 f 1 (z, 0) have negative real parts. Therefore M 1 is a manifold. For k ∈ {2, . . . , m}, M k is hyperbolically attractive on M k−1 , if M k = ∅ and the following holds. Recall that using the defining polynomials F k−1 of M k−1 , the implicit function theorem yields a unique local resolution of Z k−1 as functions of z k , . . . , z m , and thus we obtain We now require that for all z ∈ M k all eigenvalues of . . , M −1 M , then we simply write M 0 · · · M , and call this a hyperbolically attractive -chain. Such a chain is called maximal if either = m or M M +1 . Let M 0 · · · M be a hyperbolically attractive -chain. Consider for each k ∈ {1, . . . , } the following differential-algebraic system: In the limiting caseδ = 0, this corresponds to system (17). Recall that and equivalently rewrite (19) as a triplet (M k−1 , T k , R k ) with entries as follows: For a given index k, we call (M k−1 , T k , R k ) a reduced system on M k−1 , where the relevant hyperbolic attractivity relation is M k−1 M k . In order to indicate the relevance of M 0 · · · M we write (M 0 , T 1 , R 1 ) · · · (M −1 , T , R ) also for reduced systems, where M is not made explicit but relevant for the last triplet. Slightly abusing language, we speak of a hyperbolically attractive -chain of reduced systems, which is maximal if M 0 · · · M is. The following theorem is a consequence of [10, Theorem A and Corollary A], specialized to the situation at hand.
For k ∈ {1, . . . , }, the M k−1 are critical manifolds, which contain only stationary points. The systems (T k , R k ) of ordinary differential equations on M k−1 approximate invariant manifolds N k−1 in the sense of the theorem. They furthermore approximate solutions in time scale τ k of system (16), which is equivalent to our system (8). In other words, system (8) admits a succession of invariant manifolds,
For k ∈ {1, . . . , }, M k−1 is a list of real constraints defining M k−1 ⊆ R n ; T k is a list of differential equations; 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−1 , T k , R k ) represent reduced systems according to (20).
if not ϕ then 9: break 10: end if 11: 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 f k (z, 0) without the higher order terms enter the reduced systems (M k−1 , T k , R k ).
Algorithm 5 now starts with the output (T 1 , . . . , T m ) of Algorithm 1, which represents the scaled system (9). Notice that each T k already meets the specification in (20). In l.1 we define U to contain defining inequalities of the first orthant U. Starting with k = 1, the for-loop in l.3-14 successively constructs M k and R k such that in combination with T k from the input, (M k−1 , T k , R k ) forms a reduced system as in (20). The loop stops when either k = m + 1 or a test for hyperbolic attractivity in l.7 finds that M k−1 M k . We are going to discuss this test in detail in the next Sect. 3.2. Note that we maintain a matrix A for storing information between the subsequent calls of our test. In either case we arrive at a maximal hyperbolically attractive (k − 1)-chain of reduced systems given as a list ((M 0 , T 1 , R 1 ), . . . , (M k−2 , T k−1 , R k−1 )). Following the notational convention used throughout this section we set to k − 1 in l. 16. The test in l.17-19 reflects the choice of ∈ {2, . . . , m} at the beginning of this section. Finally, l.20 uses the second input (P 1 , . . . , P m ) of the algorithm to address the smoothness requirements for system (16). We are going to discuss the corresponding procedure in detail in Sect. 3.3. It will turn out that this procedure provides only a sufficient test. Therefore we issue in case of failure only a warning, allowing the user to verify smoothness a posteriori, using alternative algorithms or human intelligence. One might mention that it is actually sufficient to consider weaker finite differentiability conditions instead of smoothness, which can be seen by inspection of the proofs in [10].
From an application point of view, attracting invariant manifolds are relevant in the context of biological 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 computational perspective, hyperbolic attractivity can be straightforwardly tested using the Hurwitz criterion, as we are going to make explicit in Sect. 3.2.
The relevant results in [10], in contrast, are based on a series of hyperbolicity conditions, which are somewhat 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], which checks the number of purely imaginary eigenvalues of a real polynomial via the Cauchy index of a related rational function.

Verification of Hyperbolic Attractivity
Our definition of hyperbolic attractivity M k−1 M k refers to the eigenvalues of the Jacobians of the f * k , which cannot be directly obtained from the Jacobians of the f k [10,11]. Generalizing work on systems with three time scales [32], we take in this section a linear algebra approach to obtain the relevant eigenvalues without computing the f * k . To start with, recall the well-known Hurwitz criterion [30]: Then all complex zeros of f have negative real parts if and only if ∆ 1 > 0, . . . , ∆ n > 0. Notice that ∆ n = a n ∆ n−1 , and therefore ∆ n > 0 can be equivalently replaced with a n > 0.
We call H n the Hurwitz matrix and ∆ i the i-th Hurwitz determinant of f . Furthermore, we refer to Γ = (∆ 1 > 0 ∧ · · · ∧ ∆ n−1 > 0 ∧ a n > 0) as the Hurwitz conditions for f .
Our first result generalizes [32, Proposition 1 (ii)]. The proof is straightforward by induction.
In particular, one can choose * 1 = · · · = * k−1 = * with sufficiently small * and consider Let Γ k denote the Hurwitz conditions for the characteristic polynomial of J k . Then Lemma 3 allows to state hyperbolic attractivity M 0 · · · M as a first-order formula over the reals as follows: On these grounds, any real decision procedure [57,13,63] provides an effective test for hyperbolic attractivity. However, our formulation (21) uses a quantifier alternation ∃σ∀ in its second part. From a theoretical point of view, such a bounded alternation does not affect the asymptotic worst-case complexity, which remains single exponential [27]. From a practical point of view, we would like to continue using SMT solving over a quantifier-free logic. Our next result allows a suitable first-order formulation without quantifier alternation. Its proof combines [32, Lemma 3] with our Lemma 3.
In contrast to (21), our first-order characterization in Corollary 5 has no quantifier alternation. Note that the two top-level components of (22) establish two independent decision problems, addressing non-emptiness of the manifold and our requirement on the eigenvalues, respectively. It is easy to see that for all ∈ {1, . . . , m} and all k ∈ {1, . . . , − 1}, ϕ entails ϕ k . Thus (22) can be equivalently rewritten as k=1 (ϕ k ∧ ψ k ), explicitly: Our approach tests the conjunction in (23) using a for-loop over k in Algorithm 5. Technically, this construction ensures with the test for M k−1 M k in Algorithm 6 that M 0 · · · M k−1 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 we construct in 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.20 the validity of ψ k and return a corresponding Boolean value. We additionally return A = A k for reuse with the next iteration. The validity tests for ϕ k and ψ k in l.1 and l.20, respectively, amount to SMT solving, this time using the logic QF_NRA [2] for quantifier-free nonlinear real arithmetic. Recall the positive integer parameter p used for the precision with both Algorithm 2 and Algorithm 3. For p > 1 symbolic computation possibly yields fractional powers of numbers in the defining equations for manifolds as well as in the vector fields of the differential equations. However, such expressions are not covered by QF_NRA. When this happens, we catch the corresponding error from the SMT solver and restart with floats.
In l.1-8 of Algorithm 7 we compute β 1 , . . . , β −1 as defined in (13) and simultaneously obtain b 1 , . . . , b . In l.9-14 we compute the right hand sides of the conditions in (24) or (25), depending on the current index k. For checking those conditions in l. 16 we once more employ SMT solving, this time using the adequate logic QF_LIA [2] for quantifier-free linear integer arithmetic. Since we are aiming at nonnegative integer solutions, we introduce explicit non-negativity conditions r 1 ≥ 0, . . . , r −1 ≥ 0. In We check here a sufficient criterion for smoothness as required for (16). Output: "true" or "failed" in terms of a 3-valued logic; for all p in P k do 12: end for 14: end for 15: for all e ∈ E do 16: if not Z |= ∃r 1 . . . ∃r −1 (r 1 ≥ 0 ∧ · · · ∧ r −1 ≥ 0 ∧ (β 1 , . . . , β −1 ), (r 1 , . . . , r −1 ) = e) then 17: return failed 18: end if 19: end for 20: return true case of unsatisfiability Algorithm 7 returns "failed" in l.17. When this happens, the calling Algorithm 5 issues a warning but continues. This protocol is owed to the fact that our procedure provides only a sufficient test, which could be supplemented with other software or human intuition. In case of satisfiability, in contrast, smoothness is guaranteed, we reach l.20, and return "true." We remark that the computation time spent on E is negligible compared to the SMT solving later on. The construction of the entire set E beforehand avoids duplicate SMT instances.

Algebraic Simplification of Reduced Systems
In the output (M 0 , T 1 , R 1 ), . . . , (M −1 , T , R ) of Algorithm 5, the T k are taken literally from the input, and the M k−1 and R k are obtained via quite straightforward rewriting of the input. As a matter of fact, the computationally hard part of Algorithm 5 consists in the computation of the upper index . We now want to rewrite the triplets (M k−1 , T k , R k ) once more, aiming at less straightforward but simpler and, hopefully, more intuitive representations. The principal idea is to heuristically eliminate on the right hand side of the differential equations in T k those variables whose derivatives have already occurred as left hand sides in one of the T 1 , . . . , T k−1 . Of course, our simplifications will preserve all relevant properties of (M 0 , T 1 , R 1 ), . . . , (M −1 , T , R ), such as hyperbolic attractivity and sufficient differentiability. 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−1 = (z 1 , . . . , z k−1 ). In l.1-5 we construct a block term order ω on all variables {x 1 , . . . , x n } so that variables from Z k−1 are larger than variables from z k . This ensures that all multivariate polynomial reductions modulo ω throughout our algorithm will eliminate variables from Z k−1 in favor of variables : ω := a block term order with z 1 · · · z y 6: for k := 1 to do 7: 10: T k := () 11: end for 14: end for 15: 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 1 , . . . , z , y will heuristically give more efficient computations.
Recall that the radical ideal F of F is the infinite set of all polynomials with the same common complex roots as F . In l.8, we compute a finite reduced Gröbner basis G modulo ω of that radical. If radical computation is not available on the software side, then the algorithm remains correct with a Gröbner basis of the ideal F 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−1 . In l.12, reduction modulo ω, which comes with heuristic elimination of variables, applies once more 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 the following property: if two polynomials g 1 , g 2 coincide on the manifold M k−1 defined by M k−1 , then they reduce to the same normal form h. In particular, if g 1 vanishes on M k−1 , then it reduces to 0. We call the output of Algorithm 8 simplified reduced systems.

Back-Transformation of Reduced Systems
Let ∈ {2, . . . , m} and k ∈ {1, . . . , }. Recall that a triplet (M k−1 , T k , R k ) obtained from Algorithm 5 describes a reduced system according to (20). A corresponding simplified system (M k−1 , T k , R k ) is obtained from Algorithm 8 via an equivalence transformation on the set of equations M k−1 and further equivalence transformations modulo M k−1 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 for both these outputs scaling can be reversed using the substitution obtained with Algorithm 1. For our discussion here, we use names M k−1 , T k , R k as in the unsimplified system.
The application of σ to all components of (M k−1 , T k , R k ) yields a raw back-transformation as follows:

Algorithm 9 TransformBack
for all end for 9: In T k σ, we multiply by ε µ+dj * in order to arrive at differential equations in dyj dt . Furthermore, recall that the explicit factor δ b k in the original T k corresponds to a time scale δ b k τ . The corresponding time scale in t is given by which we make explicit by equivalently rewriting T k σ as Similarly, R k σ can be rewritten as R * k = (ẏ j = 0 | x j ∈ z k+1 ∪ · · · ∪ z m ), and we set M * k = M k σ. We call (M * 0 , T * 1 , R * 1 ), . . . , (M * −1 , T * , R * ) back-transformed reduced systems. In terms of the definitions 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 systems in 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 ε dj * f σ as (y j f /x j )σ. Let us discuss what has been gained in ((M * 0 , T * 1 , R * 1 ), . . . , (M * −1 , 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 ourselves to the positive first orthant U using defining inequalities U = {x 1 > 0, . . . , x n > 0} translates into U σ = {y 1 /ε d1 * > 0, . . . , y n /ε dn * > 0}, which is again the positive first orthant U. The manifolds M k−1 described by M k−1 lead to manifolds M * k−1 described by M * k−1 , preserving the nestedness Moreover, the system (T * k , R * k ) defines differential equations on M * k−1 . We are now faced with a discrepancy. On the one hand, we fix ε = ε * . On the other hand, the requirement thatδ be sufficiently small in Theorem 1 entails that ε be sufficiently small. It is of crucial 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, singular perturbation theory is typically concerned with asymptotic results, which are not helpful here.
In case of persistence, there exist nested invariant manifolds N * k−1 which are Hausdorff-close to M * k−1 for system (1). Moreover, the differential equations T * k associated with M * m−1 correspond to the kth Algorithm 10 TropicalMultiReduce Input: 1. A list S = (ẏ 1 = f 1 , . . . ,ẏ n = f n ) of autonomous first-order ordinary differential equations where f 1 , . . . , f n ∈ Q[y 1 , . . . , y n ]; 2. ε * ∈ (0, 1) ∩ Q; 3. p ∈ N \ {0} Output: A list ((M * 0 , T * 1 , R * 1 ), . . . , (M * −1 , T * , R * )) of triplets where ∈ {2, . . . , m}, or the empty list. For k ∈ {1, . . . , }, M * k−1 is a list of real constraints defining M * k−1 ⊆ R n ; T * k is a list of differential equations; R * k is a list of trivial differential equationsẏ = 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. We have achieved a decomposition of (1) into systems of smaller dimension. At the very least, one obtains a well-educated guess about possible candidates for invariant manifolds and reductions. For the investigation of those candidates one may check the N * k for approximate invariance using, e.g., numerical methods, or by applying criteria proposed 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 presented on the original scale. Figure 5 explains the functional dependencies and principal data flow between our algorithms graphically.

Computational Examples
Based on our explicit algorithms in the present work, we have developed two independent software prototypes realizing all methods described here. The first one is in Python using SymPy [40] for symbolic computation, pySMT [24] as an interface to the SMT solver MathSAT5 [12], and SMTcut for the computation of tropical equilibrations [39]. The second one is a Maple package, which makes use of Maple's built-in SMTLIB package [22] for using the SMT solver z3 [17]. For our computations here we have used our Python code. Computation results are identical with both systems, and timings are similar. We have conducted our computations on a standard desktop computer with an 3.3 GHz 6-core Intel 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 system from the BioModels database. The subsequent subsections showcase several further such examples in 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.

An Epidemic Model of the Bird Flu Virus H5N6
We consider a model related to the transmission dynamics of subtype H5N6 of the influenza A virus in the Philippines in August 2017 [35]. That model is identified as BIOMD0000000716 in the BioModels database, a repository of mathematical models of biological processes [34]. The model specifies four species: 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 1 , y 2 , y 3 , y 4 , respectively. The

Algorithm 6 IsHyperbolicallyAttractive
Notice that our implementations conveniently rewrite equational constraints as monomial equations with numerical right hand sides when possible. This supports readability but is not essential for the simplifications applied here, which are based on Gröbner basis theory. Comparing T 2 with T 2 , we see that the equation for x 1 x 2 in M 1 is plugged in. Similarly, M 2 is simplified to M 2 , which is in turn used to reduce T 3 to T 3 .
The back-transformed reduced systems as computed by Algorithm 9 read as follows: We compare T * 1 , . . . , T * 3 to the input system S: In the equation forẏ 1 , the monomial in y 1 is identified as a higher order term with respect to δ and discarded by Algorithm 1. In the equation forẏ 2 , the monomial in y 1 y 2 has been Gröbner-reduced to a constant modulo the defining equation in M 1 . Similarly, the equation forẏ 3 loses its monomial in y 2 y 3 by truncation of higher order terms, and in the equation foṙ y 4 , the monomial in y 2 y 3 is Gröbner-reduced to a monomial in y 3 .
Notice the explicit constant factors on the right hand sides of the differential equations in T * 1 , . . . , T * 3 . They originate from factors δ b k in the respective scaled systems T 1 , . . . , T 3 , corresponding to (8). They are left explicit to make the time scale of the differential equations apparent. We see that the system T * 2 • R * 2 is 125 times slower than T * 1 • R * 1 , and T * 3 • R * 3 is another 125 times slower. This multiple time scale reduction of the bird flu model emphasizes a cascade of successive relaxations of different model variables. First, the population of susceptible birds relaxes. As explained in the introduction, by relaxation we mean that these variables reach quasi-steady state values. This relaxation is 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 * 3 .

Caspase Activation Pathway
BIOMD0000000102 is a quantitative kinetic model that examines the intrinsic pathway of caspase activation that is essential for apoptosis induction by various stimuli including cytotoxic stress [36]. Species concentrations over time are mapped as follows:

TGF-β Pathway
BIOMD0000000101, is a simple representation of the TGF-β signaling pathway that plays a central role in 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 receptor complex-plasma membrane), lRIRII_endo (ligand receptor complex-endosome), RI_endo (receptor 1 endosome), and RII_endo (receptor 2 endosome), are mapped to differential variables y 1 , y 2 , y 3 , y 4 , y 5 , and y 6 , 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 We choose ε * = 1 5 , p = 1, and select d = (0, −4, −1, −2, −1, −5) from the tropical equilibrium. Our back-transformed reduced systems read as follows: The total computation time was 0.906 s.
The multiple time scale reduction of the TGF-β model emphasizes a cascade of successive relaxations of concentrations of different species. First, the concentration of receptor 1 relaxes rapidly. Then follows the membrane complex, and, even slower, the relaxation of receptor 2.
The multiple time scale reduction of this avian influenza model emphasizes a cascade of successive relaxations of different model variables. First, the susceptible bird population relaxes rapidly. The reduced equation T 1 and manifold M 1 suggest that the bird population dynamics is of the Allee type and evolves toward the stable extinct state. It follows the relaxation of infected human population that also evolves toward the extinct state, the end of the epidemics.

Concluding Remarks
We provided a symbolic method for automatic model reduction of biological networks described by ordinary differential equations with multiple time scales. This method is applicable to systems with two time scales or more, superseding traditional slow-fast reduction methods that can cope with only two time scales. We also proposed, for the first time, the automatic verification of the hyperbolicity conditions required for the validity of the reduction. Our theoretical framework is accompanied by rigorous algorithms and prototypical implementations, which we successfully applied to real-world problems from the BioModels database [34].
We would like to list some open points and possible extensions of our research here. Our reduction algorithm is based on a fixed scaling (8) leading to a fixed ordering of the time scales of different variables. In our reduction scheme, different variables relax hierarchically, firstly the fastest ones, then the second fastest, and lastly the slowest ones, which justifies our geometric picture of nested invariant manifolds. However, there are situations, e.g. in models of relaxation oscillations, when the ordering of time scales changes with time: variables that were fast can become slow at a later time, and vice versa. In order to cope with such situations, one would like to use different scalings for different time segments. 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 reduction may stop early in this hierarchy when hyperbolic attractivity is not satisfied at some stage. One possible reason is the presence of conservation laws, also known as first integrals, at the given reduction stage. Such conservation laws necessarily force an eigenvalue zero for the Jacobian. A theorem by Schneider and Wilhelm [52] can be employed to reduce such a setting to the hyperbolically attractive case. As for the behavior of first integrals when proceeding to the reduced system, see the discussion of the non-standard case in [33] for two time scales; an extension to multiple time scales should be straightforward. Work in progress is concerned with the introduction of novel slow variables, one for every independent conservation law of the fast subsystem, applying this to networks with multiple time scales, and approximate linear and polynomial conservation laws.
More generally, it is of interest to consider cases when hyperbolic attractivity fails but hyperbolicity still holds: In such cases, Cardin and Teixeira show there still exist invariant manifolds [10]. Testing for hyperbolicity is more involved than testing for hyperbolic attractivity, but in theory it is well understood, 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 cycles may 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 sufficient finite differentiability. Therefore, given a differential equation system and a scaling, invariant manifolds and corresponding reduced systems exist for C p functions with fixed p < ∞. Going through the details will involve intricate analysis that is left to future work.
In the introduction we sketched a Michaelis-Menten system abstracting from the known numerical values for the reaction rate constants k 1 , k −1 , k 2 . It would be indeed interesting to work on such parametric data. In the presence of parameters, one would consider effective quantifier elimination over real closed fields [14,63,31,55] as a generalization of SMT solving. Robust implementations are freely available [8,18] and well supported. They have been successfully applied to problems in chemical reaction network theory during the past decade [56,62,7,19]. Such a generalization is not quite straightforward. With the tropical scaling in Sect. 2.2, Algorithm 2 would introduce logarithms of polynomials in the parametric coefficients, which is not compatible with the logical framework used here. Similar tropicalization methods, which are unfortunately not compatible with our abstract view on scaling in Sect. 2.1, require only logarithms of individual parametric coefficients [51]. Such a more special 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 automatic strategies for determining good default values for ε * and for choices of d from the tropical equilibration in Algorithm 3.

A.1. Failure of Tropicalization With an Unbalanced Monomial
BIOMD0000000609 describes the metabolism and the related hepatotoxicity of acetaminophen, a pain killer [48]. The species concentrations over time of Sulphate__PAPS, GSH, NAPQI, Paracetamol_APAP, and Protein_adducts are mapped to differential variables y 1 , y 2 , y 3 , y 4 , and y 5 , respectively. The input system is given by Since there is only one monomial on the right hand side of the equation forẏ 5 , equilibration is impossible. This causes Algorithm 4 to return in l.23 a disjunctive normal form Π equivalent to "false", which describes the empty set. Hence Algorithm 3 returns ⊥, and Algorithm 1 returns the empty list. The total computation time was 0.006 s.

A.2. Failure of Hyperbolic Attractivity due to an Empty Manifold in the First Orthant
BIOMD0000000726 examines the transmission dynamics of rabies between dogs and humans [49]. The mapping of the model variables over time and our differential variables is as follows: Species Differential variable Description S_d y 1 susceptible dogs E_d y 2 exposed dogs I_d y 3 infectious dogs R_d y 4 recovered dogs S_h y 5 susceptible humans E_h y 6 exposed humans I_h y 7 infectious humans R_h y 8 recovered humans The input system is given by
We choose parameters ε * = 1 2 , p = 1, and point d = (2, 1, 1). Algorithm 5 returns an empty list, since the test for hyperbolic attractivity fails in Algorithm 6, l.20, even though the manifold M 1 , defined by M 1 = [− 37 40 x 1 x 2 + x 1 = 0], is not empty in the first orthant, as has been ensured in l.1. Obviously, the simplified and back-translated systems are empty lists as well. The total computation time was 0.453 s.