A Survey of Some Methods for Real Quantifier Elimination, Decision, and Satisfiability and Their Applications

Effective quantifier elimination procedures for first-order theories provide a powerful tool for generically solving a wide range of problems based on logical specifications. In contrast to general first-order provers, quantifier elimination procedures are based on a fixed set of admissible logical symbols with an implicitly fixed semantics. This admits the use of sub-algorithms from symbolic computation. We are going to focus on quantifier elimination for the reals and its applications giving examples from geometry, verification, and the life sciences. Beyond quantifier elimination we are going to discuss recent results with a subtropical procedure for an existential fragment of the reals. This incomplete decision procedure has been successfully applied to the analysis of reaction systems in chemistry and in the life sciences.


T. Sturm
both of which are excellent. In Sect. 3 we discuss the application of real quantifier elimination and variants thereof to geometric theorem proving, verification, and bifurcation analysis in the life sciences. In Sect. 4 we introduce an incomplete heuristic satisfiability checking procedure for the reals in order to consider in Sect. 5 larger examples from the life sciences and also from chemistry. The applications here are very similar to the life science applications discussed in detail in Sect. 3. Therefore 2 we only summarize the models and focus on the current limit of problem sizes tractable by symbolic methods in that area. We end with some concluding remarks in Sect. 6.
This article is a survey of my own work, partly with collaborators. There is no claim that it contains unpublished original results. I allowed myself to take over single sentences or short passages from my own publications without explicitly mentioning this.

Introduction
The following formal statement ϕ over the reals asks whether or not one can find for all x ∈ R some y ∈ R such that a certain polynomial p ∈ Z[a, b, x, y] is strictly positive while another such polynomial q is zero or negative: 3 ϕ= ∀x∃y( p > 0 ∧ q ≤ 0) where p= x 2 + x y + b, q= x + ay 2 + b.
We have to expect that the validity of ϕ depends on the choices of real values for the parameters a and b. Even for readers familiar with real algebraic geometry, a solution is probably not easy to see. 4 It gets easier when considering 1 https://www.usna.edu/Users/cs/wcbrown/research/ISSAC04/Tutorial.html. 2 And also due to reasonable page limits imposed by the editors. 3 This instructive example has been used by Hoon Hong in a colloquium talk at the University of Passau in 2005. 4 The free tool GeoGebra greatly supports intuition by plotting the ranges of the two inequalities with respect to (x, y) ∈ R 2 and allowing transformations of that plot via sliders for the values of a and b; our example can be found at https://www.geogebra.org/m/ f7UG3zpq. ¬ϕ, which is equivalent to ∃x∀y( p ≤ 0 ∨ q > 0). When a ≥ 0, we can choose x = −b + 1 to satisfy q > 0. When b ≤ 0, we can choose x = 0 to satisfy p ≤ 0. Thusφ= a ≥ 0 ∨ b ≤ 0 implies ¬ϕ. Equivalently, ϕ implies ϕ = a < 0 ∧ b > 0. Vice versa, it is not hard to see that ϕ implies ϕ: Given x ∈ R choose y ∈ R with sgn(y) = sgn(x) and |y| = (x + b)/a . Formally, we are considering interpreted first-order logic with equality over a finite language 5 L = (0, 1, +, −, ·, ≤) where all symbols have their usual interpretations over the reals. We admit ourselves to use relations "<", ">", "≥", " =" and integers as short notations for corresponding quantifier-free formulas and terms, respectively. Given a first-order L-formula ϕ, our goal is to compute a quantifier-free L-formula ϕ such that R | ∀(ϕ ←→ ϕ ), where "∀" denotes the universal closure. This is called real quantifier elimination. We call ϕ an L-sentence, if all occurring variables are quantified, in other words, if there are no parameters. When applying quantifier-elimination to a sentence ϕ, the obtained quantifier-free formula ϕ will not contain any variables and can be straightforwardly simplified to either "true" or "false". This way, real quantifier elimination establishes in particular a decision procedure. In the special case that ϕ is of the form ∃ψ, where "∃" denotes the existential closure and ψ is quantifierfree, quantifier elimination establishes a satisfiability checking procedure. In Sect. 4 we will study a satisfiability checking approach that is not based on quantifier elimination.
In symbolic computation our formal framework is often called the Tarski Algebra, while in satisfiability modulo theories solving (SMT) it is referred to as nonlinear real arithmetic. One might mention that from a point of view of algebraic model theory we are considering the model class RCF = Mod(Th(R)) of real closed fields, containing all L-structures in which exactly the same first-order L-sentences hold as in R. Natural enumerable axiomatizations of RCF can be found in [46]. By the Löwenheim-Skolem Theorem, RCF is a proper class containing models of arbitrary infinite cardinality. Among those models are the real numbers R, the countable model of all real algebraic numbers, or the non-Archimedian extension field R(∞), in which ∞ is algebraically free over R but strictly greater than all real numbers. 6 Since from the point of view of first-order L-formulas, the models in RCF cannot be distinguished, it is safe to assume with the application of real quantifier elimination that the domain of the computation is R, or any other model in RCF. We shall see that under the hood of real quantifier elimination and decision procedures, non-standard models, in particular the ones mentioned above, play an important role.

History
The first real quantifier elimination procedure was published by Tarski at the end of the 1940s [65]. 7 As early as 1954, concluding remarks in a technical report by Davis to the US Army on an implementation of a similar procedure for Presburger arithmetic point at a very early interest in software implementations of real quantifier elimination [16].
During the 1970s Collins developed the first elementary recursive real quantifier elimination procedure [11,12], which was based on cylindrical algebraic decomposition (CAD). An implementation by Arnon was available around 1980 [1]. CAD has undergone many improvements since and establishes an active research area until today [7,14,40]. Hong reimplemented CAD in an interactive software Qepcad [13]. Qepcad developed into Qepcad B, which is now maintained by Brown [6].
From the mid 1980s to the early 1990s there was a strong interest in the asymptotic worst-case time complexity of the decision problem for RCF. In 1988 Davenport-Heintz [15] and Weispfenning [72] independently showed that it is doubly exponential. Weispfenning's article actually brought even stronger results: first, it showed that the decision problem is doubly exponential already for linear formulas, where the total degree in all quantified variables does not exceed 1. Next, considering finer complexity parameters than the input word length it showed 5 Alternatively called a signature. 6 Certainly all these models are unique only up to isomorphisms. 7 Tarski's results are actually around 10 years older, but publication was delayed due to World War II. that the problem for linear formulas is doubly exponential in the number of quantifier alternations but only singly exponential in the number of quantifiers (when bounding the number of alternations) and only polynomial in the number of variables (when bounding the number of quantifiers). Finally, it came with a corresponding effective procedure called virtual substitution (VS), where the idea is for the elimination of an existential quantifier to formally substitute in a generalized sense sufficiently many parametric zeros of certain polynomials contained in the input formula. Weispfenning's complexity results were followed by further research by Grigoriev, Renegar, Basu-Pollack-Roy, and others developing entirely new quantifier elimination procedures with strong theoretical results taking into consideration even finer complexity parameters like polynomial degrees or coefficient sizes [2,27,48]. 8 Those asymptotically fast procedures remained unimplemented, and there are in fact substantial arguments that the trade-off points are too large for the procedures to be relevant in practice [30].
VS, in contrast, turned out to be practical. It was implemented by Weispfenning's students in the computer logic system Redlog [17,19,61], which has developed since into an established tool with more than 350 citations in the scientific literature for successful applications mostly in the sciences. 9 The focus is on low-degree problems with few quantifier alternations and with a significant number of parameters, as we know today that CAD-and thus Qepcad B-is doubly exponential in all variables. All examples discussed in this article have been computed with an implementation of VS in Redlog that is limited to formulas where for the elimination of a quantifier the total degree of the corresponding quantified variable does not exceed 2. Nevertheless we will encounter significantly higher degrees, because VS implementations are supplemented with powerful simplification strategies; see, e.g., [20] and [36,Section 5]. The nonlinear bound of 2 is based on work by Weispfenning in the early 1990s, which also theoretically discussed arbitrary degree bounds [73]. Only recently, Košta has theoretically improved Weispfenning's framework for higher degrees and provided a generic implementation in Redlog, which can be instantiated with substitution tables up to a chosen degree bound [35]. Highly optimized such tables are available up to degree 3; providing practically useful tables for higher degrees is a realistic but challenging task.

Degree Bounds
It is important to understand that the degree bounds for VS have to be satisfied for the quantified variable with the elimination of every single quantifier. Quantifiers are essentially eliminated one by one. Except for Weispfenning's original linear case [39,72], where there are no products of quantified variables at all, the elimination of a quantifier will in general increase the degrees of other quantified variables. As a consequence, one cannot tell by inspection of the input whether or not quantifier elimination will succeed subject to the current degree bound. This in turn makes it hard to give meaningful upper complexity bounds. In fact, Weispfenning's bounds discussed above have been proved only the linear case. From a practical point of view, hundreds of meaningful problems have been solved during the past 25 years with a degree bound of only 2. Those computations have given the following impression: Firstly, input formulas modelling real world situations, in contrast to randomly generated ones, are surprisingly well-behaved concerning the increase of degrees during elimination. Secondly, the observed computation times appear compatible with the theoretical results for the linear case concerning improved performance with bounded quantifier alternation and significantly less sensitivity to numbers of parameters than to numbers of quantified variables.

Virtual Substitution in a Nutshell
Every L-formula can be equivalently transformed into a prenex L-formula consisting of a sequence of quantifiers followed by a quantifier-free formula. Without loss of generality, all right hand sides of equations and inequalities are 0. Given such a prenex formula, VS successively eliminates the quantifiers starting with the innermost one. Using the equivalence ∀x 1 ψ ←→ ¬∃x 1 ¬ψ we may w.l.o.g. assume that the innermost quantifier is an existential one. More generally, there is an innermost block of finitely many existential quantifiers to be eliminated, possibly with further quantifier blocks outside. The idea is to replace ∃x 1 with a finite disjunction such that where E = {. . . , (γ, t), . . .} is a finite elimination set containing sufficiently many test points t. Notice that before the next elimination step ∃x 2 can be moved into the disjunction. This explains the better performance of the procedure with bounded quantifier alternation.
Fixing an interpretation of all variables except x 1 it is not hard to see that the set of satisfying values for ψ with respect to x 1 is a finite union of intervals. Furthermore, the upper and lower bounds of those intervals are −∞, ∞, or zeros of the left hand side polynomials in ψ. On these grounds, the test points t in E are derived as formal solutions of the polynomials in ψ, and the corresponding guards γ guarantee the existence of those solutions. For instance, elimination set elements produced by ax 2 − 3x + 7 ≤ 0 would include , a = 0, 7 3 .
However, standard term substitution is not sufficient here: For instance, we must substitute quotients and square roots, both of which are not in our language L. Furthermore, with strict inequalities we need, e.g., t − ε with infinitesimal ε, and similarly we need ∞. Those substitutions must be carried out in such a way that the substitution results contain only symbols from L. The key idea is to use a virtual substitution, which does not map terms to terms but atomic L-formulas to quantifier free L-formulas.
: Quotients: Infinity: Positive infinitesimals: Formal solutions of quadratic equations: It turns out that virtual substitution is actually strong enough to substitute abstract formal representations of roots of polynomials of arbitrary degrees with suitable guards. This explains the role of the substitution tables with the new generic implementation mentioned at the end of Sect. 2.2 [35].

Variants of Real Quantifier Elimination
There are numerous variants and extensions of the concept of quantifier elimination. We are going to discuss here and use with our applications in the following section three of them.

Extended Quantifier Elimination
Extended quantifier elimination [73] does not form a disjunction over the virtual substitution results. Instead, they are kept separate and associated with formal assignments of the substituted test points: There is a precise semantics as follows: For fixed choices of all parameters, whenever some left hand side condition holds, then ∃x 1 ψ holds, and the corresponding right hand side assignment is one sample solution. As a simple example consider Notice that our specification of the semantics above would allow the simplification component of virtual substitution to delete one of the first two rows from the result. In the last row we can see that an infinite element has been introduced as a test point. In the course of the elimination of further variables this can happen several times, which corresponds to a stack of extensions of real closed fields, where every newly introduced infinite element is infinite not only compared to standard real numbers but also compared to the already existing infinite elements. Therefore all such elements are carefully distinguished using indices.
With the elimination of several existential quantifiers one obtains assignments for all corresponding variables, on which one can perform back-substitution in the style of linear programming. For fixed choices of parameters, non-standard elements can be efficiently eliminated from the sample solutions in a post-processing step [36]. Finally, with universal quantifiers we have the following dual semantics: For fixed choices of all parameters, whenever some left hand side condition does not hold, then ∀x 1 ψ does not hold, and the corresponding right hand side assignment is one counterexample.

Generic Quantifier Elimination
The combinatorial complexity of virtual substitution mostly arises from case distinctions on the vanishing of certain coefficients. We have seen an example with the test points derived from ax 2 − 3x + 7 ≤ 0 in Sect. 2.4. The idea of generic quantifier elimination [18,54,62] is to avoid those case distinctions when the corresponding coefficient does not contain any quantified variables. It assumes the generic case that the coefficient does not vanish. Those assumptions are collected in a global theory , which establishes part of the output: Again, there is a precise semantics stating that the elimination result obtained this way is correct for all choices of parameters that satisfy , formally: Consider again our simple example from the previous section: Notice that we only assume negations of equations but no ordering inequalities. Hence cannot become inconsistent. Furthermore, the set of choices for the parameters that does not satisfy the assumptions in has a lower dimension than the parameter space.

Positive Real Quantifier Elimination
With many natural problems, variables may be assumed to be positive. The idea of positive quantifier elimination is to have a variant of the elimination procedure which systematically exploits that assumption [57,64]. First implementations of that idea, as the one used with our life science application in Sect. 3.3, implicitly assumed that all occurring variables are strictly positive. With a more modern approach now taken in Redlog relevant positivity conditions are explicitly encoded in the input formula and extracted by the VS procedure. This way positive quantifier elimination turns into an optimization of the procedure which is transparent for the user.

Geometry
We focus here on geometric theorem proving. Related applications of real quantifier elimination methods include computational geometry [58] and solid modeling [59,60]. Theorems of elementary geometry have traditionally been considered an important test case for the scope of methods in automatic theorem proving. In particular, they have stimulated a variety of algebraic techniques for their solution. There is an established style of algebraic modeling with parameters u and dependent variables x, which can be found already in Hilbert's Grundlagen der Geometrie 10 [29], originally published in 1899. Using a suitably positioned coordinate system, a geometric configuration and a claimed geometric property of that configuration are translated into a formula x] describe the configuration and the claim, respectively. The parameters u describe geometric entities that can be freely chosen, e.g., three corners of a triangle, while the dependent variables x correspond to entities that are uniquely determined by the existing construction, e.g., the center of the circumcircle of that triangle. Around 1977, Wu discovered the Wu-Ritt method 11 [74,75] 12 adapting a method for differential algebra developed by Ritt [49,50]. Chou made important contributions to the development and application of the Wu-Ritt method. Most of Chou's work has been summarized in an impressive monograph including more than 500 geometric theorems automatically proved by this method [9]. Kapur [34] developed a prover based on a radical membership test via Gröbner Basis [8] computations. Kutzler and Stifter [37] developed another prover also based on Gröbner techniques. They modified the notion of reduction to a pseudo reduction, which is somewhat similar to the notion of reduction used with the Wu-Ritt method. Wang [69,70] developed complex elimination methods based on ideas by Seidenberg [51][52][53] but also inspired by Wu-Ritt. All those methods, which have turned out to be quite successful, try to prove ϕ as a statement about complex numbers. Since ϕ is a universal assertion, the validity of ϕ over the complex numbers entails the validity of ϕ over the reals and thus an automatic proof of the original geometric assertion. If, in contrast, ϕ turns out to be false over the complex numbers, no decision on the validity of the geometric statement can be made. It is due to the historical focus on complex methods that the algebraic translations introduced above are purely equational. From a geometric point of view it makes sense to also have inequalities comparing, e.g., lengths. Employing real quantifier elimination will allow us to do so.
Before studying an example in detail, we want to have a closer look at the available freedom with the parameters u. We had mentioned earlier than one would, e.g., introduce (u 1 , u 2 ), (u 3 , u 4 ), (u 5 , u 6 ) as corners of a triangle in the 10 Foundations of Geometry. 11 Also known as Wu-Ritt reduction or Wu's method. 12 This article appeared in two different journals in 1984 and 1986. The reason was limited availability of the first one. Both these publications are in English. We thank Xiao-Shan Gao for clarifying this. plane. However, this requires certain non-degeneracy conditions stating, e.g., that those three points be not collinear. When they are, they do not form a triangle. This breaks the geometric construction and in general also the validity of ϕ. Furthermore, for a geometric theorem to hold it might be necessary to exclude certain valid configurations like parallel lines or rectangular triangles. Since such cases are not really degenerate, non-degeneracy conditions have also been called subsidiary conditions. It is in general by no means easy to find and add all necessary subsidiary conditions for the validity of a theorem. Illustrating examples for hidden conditions have been given by Chou [9, pp. 43-44].

T. Sturm
As a matter of fact, the results produced with the above-mentioned algebraic methods are correct also only under certain assumptions in the shape of negated equations. Those assumptions would state, e.g., that during the algebraic computations certain parametric polynomial coefficients c(u) do not vanish. They are collected and explicitly provided with the output. We are going to apply here the generic real quantifier elimination from Sect. 2.5.2, which also produces such assumptions. It is a most amazing fact that the assumptions made by the proof methods mostly coincide with the subsidiary conditions necessary for the validity of the corresponding geometric theorem [9,18,62]. That is, an apparent weakness of the algebraic procedures turns out as an adequate tool for actively discovering missing hypotheses. The presence of this phenomenon with many different algebraic methods suggests that the approach to translate geometry to algebra is very natural.
We are going to study a variant of the Steiner-Lehmus Theorem taken from McPhee et al. [41]. In its original form the theorem states that any triangle with two equal internal bisectors is isosceles. It had been firstly mentioned in a letter by Christian Ludolf Lehmus to Charles-François Sturm in 1840, asking for a purely geometric proof. Sturm passed the question on to several mathematicians, and Jakob Steiner provided such a proof. The contrapositive of the original Steiner-Lehmus Theorem follows immediately from the following, stronger, variant; compare Fig. 2b: (variant). Assume that ABC is a triangle such that AB > AC. Then the angle bisector from B to AC is longer than the angle bisector from C to AB. That is, the longer bisector goes to the shorter side.

Steiner-Lehmus Theorem
For our algebraic translation, we take coordinates A = (−1, 0), B = (1, 0), and C = (u 1 , u 2 ) for the triangle. We assume w.l.o.g. that both C and the center M = (0, x 1 ) of the circumcircle of ABC lie above AB: The center M = (0, x 1 ) of the circumcircle is determined by its radius: We now construct the bisector of ACB: The point V = (0, x 2 ) is the lower extremity of the circumcircle: and X = (x 3 , 0) is the intersection between CV and AB. That is, the points CXV are collinear: The line segment CX is the bisector on the side AB. For constructing the bisector on AC using the above technique, we would have to guarantee that both M and B lie on the same side of AC within the circumcircle, which would introduce tedious case distinctions. We thus construct the bisector on AC more straightforwardly: Let W = (x 4 , 0) be the point on the line AB that lies left of B with distance BC: x 6 ) of the bisector on AC is the unique point with equal distance to both W and C and with AYC collinear: Finally we add to our hypothesises the fact that AC is shorter than AB: 2 < 2 2 and state the conclusion that CX is shorter than BY : In terms of these definitions of h 1 , …, h 7 and g, our algebraic translation of the Steiner-Lehmus Theorem reads as follows: Application of generic quantifier elimination to ϕ yields a quantifier-free formula ϕ containing 231 atomic formulas plus the following assumptions: The first assumption can be rewritten as (u 1 −1) 2 +u 2 2 = 4, i.e., the distance between B and C is different from 2. In other words, the triangle ABC is not isosceles. The assumption u 1 = 0 also states that the triangle is not isosceles, and u 2 = 0 is a natural non-degeneracy condition for the triangle. While with other examples one frequently obtains ϕ = true, the situation here is more complicated. However, applying cylindrical algebraic decomposition to ∀u 1 ∀u 2 −→ ϕ we obtain "true", which means that ϕ is indeed true under the assumptions already made. The overall computation time for this example is around 2 s.

Verification
Given a system, typically modeled by differential equations, and given a property, the verification problem asks whether or not the system satisfies the property. We are going to prove collision avoidance for two cars under cruise control [56]. The rear car uses a cruise control law that actively adjusts its acceleration based on its own velocity, acceleration, the relative velocity of the leading car, and the gap between the two cars. Proving collision avoidance means showing that the two cars will not collide assuming that the cruise control is activated in a safe initial configuration. The cruise control law is taken from the leader control developed in [26] and also discussed in [47,66]; see also the Berkeley Path project. 14 The system dynamics and the control law are given as a system of ordinary differential equations in Fig. 3, whereḟ := d f / d t is the common short notation for the derivative with respect to time.  For the initial state we choose a gap, assume that the rear car is driving at constant speed, and introduce parameters c 1 and c 2 for the speeds of the leading and the rear car, respectively: A state is defined to be safe whenever there is a positive gap between the two cars. Of course, there is an additional safety distance, which is not made explicit in our model: We use the certificate-based approach for verification [28,45,66], where the idea is to discover an invariant set Inv, which is defined by the following properties: 15 (I 1 ) Init ⊆ Inv (I 2 ) Inv ⊆ Safe (I 3 ) The system dynamics does not admit the system to leave the set Inv.
Such a set Inv serves as a certificate for safety. We make an ansatz assuming Inv can be defined in terms of the non-negativity of a suitable linear combination of entities in our model. We may normalize one of the linear factors to 1: The system dynamics can be over-approximated in Tarski algebra as follows: In terms of these definitions we formally write down the certificate conditions (I 1 )-(I 3 ): The derivative of p in ϕ 3 can be symbolically computed using the system dynamics, which specifies that accelerations are derivatives of velocities, and the control law, which explicitly yields the derivative of the acceleration a of the rear car: On the basis of these preparations the following formula with parameters c 1 and c 2 corresponding to initial speeds of the two cars asks for the existence of linear factors λ 1 , …, λ 4 such that Inv is in fact an invariant: Redlog's virtual substitution eliminates all variables except λ 3 , where it finally fails due to the degree bound. 16 The result is a disjunction ϕ 584 of 584 subformulas altogether containing 33,365 atomic formulas in our parameters 15 Notice that there is natural correspondence between our formulas and their respective sets of satisfying values in a real space of suitable dimension. For the sake of a concise description we do not make this explicit. 16 Recall that with the elimination of a block of either existential or universal quantifiers the order of the variables can be freely chosen, and Redlog uses heuristics to do so. Using CAD-based techniques [5] and the assumption that both c 1 and c 2 are positive, we manage to equivalently simplify the 33 out of the 584 subformulas to the following condition, which is quantifier-free and contains only the initial velocity c 2 of the rear car: It is easy to see that for positive c 2 our ϕ 33 can be equivalently simplified to ϕ 33= c 2 ≤ 10 √ 3 + 15. Since R | ϕ 33 ←→ ϕ 33 , R | ϕ 33 −→ ϕ 584 and R | ϕ 584 ←→ ϕ, our condition ϕ 33 is sufficient for ϕ. We have established collision freedom for initial v ≤ 10 √ 3 + 15 ≈ 32.32 of the rear car without imposing any restriction on the dynamics of the leading car. The entire computation takes about 1 min of CPU time. Figure 4a shows a simulation of the system for admissible initial velocities v f = 20 and v = 15.
However, there are limitations with our model. We have assumed that the system dynamics in Fig. 3 remains always satisfied. This appears natural, because it is about physical limitations, viz. maximal possible acceleration and braking, of the cars. Let us consider a scenario in which the leading car starts at velocity v f = 15 decelerating at a f = −3, and the rear car starts at velocity 30 with zero acceleration and an initial gap of 10. This satisfies both the system dynamics and our computed safety condition. Surprisingly, a simulation plot in Fig. 4b shows that the cars collide. This is compatible with our model, because the system dynamics is violated: The control law causes the rear car to brake harder than −5. There are several ways out. One could switch off cruise control and alert the driver when its computed acceleration values cannot be realized in reality. One could also refine our model to assume that in such cases acceleration goes exactly to the physical limits.

Life Sciences
In chemistry and systems biology there is a considerable interest in Hopf bifurcations due to their relationship with the occurrence of oscillations. Although the nature of that relationship is very subtle [3,43] there is a consensus that relevant oscillations typically occur in the presence of Hopf bifurcations, and considerable work has been done to investigate various chemical and biological systems with respect to Hopf bifurcation fixed points [24,42,44,55,67].
Numerical packages like AUTO 17 or XPPAUT 18 can locate Hopf bifurcations but cannot prove their absence. Furthermore, they might miss bifurcations that exist only for small ranges of parameters, which is not only of theoretical interest in the context of chemical and biological reaction systems. Whereas existence of Hopf bifurcations is known to be decidable [25,33,38,43] early symbolic investigations carried out for specific parameterized polynomial vector fields arising from larger examples [3,25] were not fully algorithmic but required a sequence of symbolic computations intervened with ad hoc insights and decisions made by a human, and sometimes with sophisticated coordinate transformations. Our example discussed here dates back to 2008. It is to our knowledge the first fully automatic symbolic investigation of Hopf bifurcations for a biologically relevant model [57,64]. 19 Boulier et al. have studied a model related to the gene regulatory network controlling the circadian clock of a certain unicellular green alga [3]. Figure 5 gives an overview. We are going to summarize the development of a simplified system of differential equations in [3], which will serve as a starting point for our methods discussed here. The following system describes the dynamics of the reaction system in Fig. 5a: Lowercase greek letters are constants, and uppercase roman letters are functions in time. All constants except γ 0 correspond to the reaction rates in Fig. 5a; G ∈ [0, . . . , γ 0 ] is the averaged state of the gene, where G = 0 denotes that a polymer is bound to the gene promoter, and G = γ 0 denotes that this is not the case; M is the corresponding concentration of active messenger RNA; P, P 2 , …, P n are concentrations of the protein and its polymers; A i is defined as 1 ε κ − i+1 P i+1 − κ + i+1 P i P , where the factor 1 ε with small positive ε expresses the assumption that the polymerization reactions are fast compared to the other reactions. Note that all occurring constants and functions are strictly positive.
Using quasi-steady-state assumptions thatṖ 2 , …,Ṗ n are small and furthermore neglecting terms with a factor ε one approximateṡ P = β M − δ P P + n ϑ(γ 0 − G) − αG P n and P n =ᾱ P n , whereᾱ = κ + 1 · · · κ + n−1 This yields the following simpler system: At that point Boulier et al. apply various transformations, which we only sketch to the extent necessary to understand the variables occurring in their following further simplified system: The product αᾱ has been replaced with α, P has been replaced with ϑ α 1 n P, and new variables λ and μ have been introduced for terms f − b and b γ 0 , respectively. The factor δ M of M in the second equation has been eliminated 17 http://indy.cs.concordia.ca/auto/. 18 http://www.math.pitt.edu/~bard/xpp/xpp.html. 19 More than 10 years earlier, Hong-Liska-Steinberg had applied quantifier elimination for testing stability of systems of differential equations [31]. In that course α n ϑ n−1 has been replaced with α. Finally, M has been rescaled so that δ P = β, which is now called δ; again other constants had to be adapted. 21 Notice that the new constant λ, in contrast to all other constants and functions, is not necessarily positive. In fact, for λ > 0 messenger RNA transcription is enhanced when polymer is bound to the gene promoter, while for λ < 0 it is reduced. This corresponds to the postive and negative feedback loops mentioned in Fig. 5, respectively.
From the point of view of our quantifier elimination methods that last system can be analyzed for fixed choices of n ∈ N. We use a method by El Kahoui and Weber to generate from that system first-order formulas over the reals [33]. The method is applicable whenever the vector field consists of polynomials in the variables and parameters. For n = 9 we obtain We have introduced real variables v 1 , v 2 , and v 3 for the functions G, M, and P, respectively. The first line of our formula equates the vector field of the differential equations to zero thus expressing that the system is in steady state. In the third line, 2 and 1 are the corresponding Hurwitz determinants of the characteristic polynomial of the Jacobian of the vector field. In general, the condition n−1 = 0 ∧ n−2 > 0 ∧ n−3 > 0 characterizes a Hopf bifurcation, and n−4 > 0 ∧ · · · ∧ 1 > 0 in addition guarantees an empty unstable manifold; in addition, the absolute summand of the characteristic polynomial used for the Hurwitz determinants must be positive, which is not an issue with our model [33].
Hence our ϕ 9 is equivalent over the reals to the existence of a Hopf bifurcation with empty unstable manifold, and quantifier elimination can in principle yield necessary and sufficient conditions in terms of the parameters. In practice this turns out to be not feasible with existing implementations of real quantifier elimination. Existing CAD-based implementations cannot decompose the corresponding 9-dimensional space in reasonable time, while VS-based implementations cannot cope with the high degrees occurring here. A bit surprisingly, VS does succeed, in contrast, on the existential closuresφ n= ∃ϕ n , although there are even more variables to eliminate. The reason is that there is more freedom with the heuristic choice of the elimination order, and powerful simplification techniques [20] help to get around solving with high degrees.
Our results for n ∈ {2, . . . , 10} are collected in a table in Fig. 6a. 22 We use extended positive quantifier elimination, which assumes that all occurring variables are strictly positive, in combination with a case distinction on the sign of λ. In the satisfiable case we obtain exact symbolic sample solutions. Figure 6b gives float approximations of such a sample solution forφ 9 .
Alternative reductions to real quantifier elimination of our model have been discussed in [71]. The vast majority of chemical and biological models as they can be found, e.g., in the BioModels database 23 is considerably larger with respect to both dimension and degrees than our example here. in Sect. 5 we are going to discuss such models based on a different translation into logic. For deciding satisfiability we will use a heuristic method, which is not based on quantifier elimination. We are going to introduce this method in the following section.

Subtropical Methods
We are interested in heuristically finding zeros (z 1 , . . . , z n ) ∈ R n of large polynomials f ∈ Z[x 1 , . . . , x n ] with z 1 > 0, …, z n > 0. In other words, our zeros should be located strictly in the first hyper-quadrant.
Plugging in the point (1, . . . , 1) yields f (1, . . . , 1) = y, with three possibilities for the sign of y. If y = 0, then we are done. If y > 0, then we consider instead − f , which has the same zeros as f . Therefore we may assume in 22 For independent reasons one knows with the first Hopf bifurcation at n = 9 that there will be Hopf bifurcations for all n > 9. 23 http://www.ebi.ac.uk/biomodels-main/. the sequel that f (1, . . . , 1) = y < 0. In Sect. 4.1 we are going to find p ∈ R n with all positive coordinates such that f (p) > 0. In Sect. 4.2 we are going to subsequently apply the intermediate value theorem to construct a zero. A detailed formal discussion of our method summarized here can be found in [63].

Heuristically Finding a Positive Point
To get a first idea consider g = −128x 10 1 x 2 2 + 2x 3 1 x 9 2 − 4096. For g to become positive we need a point p ∈ R 2 where the second term x 3 1 x 9 2 , which has a coefficient with positive sign, dominates in the sense that it has a large absolute value compared to x 10 1 x 2 2 as well as compared to the coefficients −128 and −4096. This can be achieved by making x 2 larger than x 1 by a sufficient order of magnitude, like p = (2, 2 2 ) with g(p) = 2,093,056 > 0.
We want to generalize this idea. Our discussion will be accompanied by the slightly more complicated polynomial Notice that f (1, 1) = −3 < 0 as required by our overall framework. The slope pictured in Fig. 7a is the real variety of f , i.e., the set of all its zeros. We switch to a more abstract view of f taking into consideration only its terms and the sign of the corresponding coefficients. The terms are determined by the respective exponent vectors, and we remember the signs of the coefficients by forming a corresponding partition. Formally this gives us the positive and the negative support of f : 24 The support of f is pictured in Fig. 7b together with its convex hull, the Newton polytope of f . For our purposes here, the Newton polytope is the set of vertices of the convex hull, denoted by newton( f ) ⊆ supp( f ). We have (0, 2) ∈ newton( f ) ∩ supp + ( f ). Since that point is in the convex hull there exists a separating hyperplane H with a normal vector n = (−3, −2) pointing outwards. In the real world the coordinates of n give a suitable ratio such that for (t −3 , t −2 ) with sufficiently large t the term x 2 2 corresponding to (0, 2) dominates in the sense of the introductory example. Plugging in increasing powers of 2 for t we find that already for p = (2 −3 , 2 −2 ) = 1 8 , 1 4 we obtain f (p) = 1087 16,384 > 0. In Fig. 7a we see the segment (t −3 , t −2 ) ∈ R 2 | t ∈ [1,2] of the corresponding moment curve, which converges to the origin as t goes to infinity.
The separating hyperplane H and its normal vector n can be efficiently found using linear programming techniques; see [63] for details. The method is incomplete. It fails when the positive support is non-empty but the Newton polytope contains exclusively points from the negative support; formally, supp + ( f ) = ∅ but newton( f ) ⊆ supp − ( f ) and thus supp + ( f ) ∩ newton( f ) = ∅. If, in contrast, supp + ( f ) = ∅, then f is negative definite on the first hyper-quadrant, and we can be certain that there is no zero with positive coordinates.

Constructing a Zero
We have found f (1, 1) < 0 < f 1 8 , 1 4 . By the intermediate value theorem there is at least one zero on the connecting line segment from (1, 1) to 1 8 , 1 4 . We straightforwardly write down that (x 1 , x 2 ) is a zero of f lying on that line segment: Plugging the equations for x 1 and x 2 into the first equation we obtain a univariate equation in y, which has at least one solution for y ∈ ] 0, 1[:

Further Applications in Chemistry and the Life Sciences
In Sect. 3.3 we have studied Hopf bifurcations for a biological reaction system using results on Hurwitz determinants. The variables in our real first-order model corresponded to concentrations of species and reaction rates. Using ideas from stoichiometric network analysis [10], it is possible to analyze the system dynamics in flux space instead of concentration space and to represent the space of steady states with a combination of subnetworks. The subnetworks form a convex cone in flux space [68]. There are techniques for decomposing that cone and limiting the search for Hopf bifurcations to lower dimensional facets [22]. For our purposes here it is sufficient to understand that we are going to obtain a whole family of formulas for a single model corresponding to those facets. Methods for detecting Hopf bifurcations using similar approaches have been used in several hand computations in a semi-algorithmic way for parametric systems, the most elaborate of which is described in [25].
The formulation of a reaction system using convex coordinates 25 in flux space implicitly limits that system to steady states. Technically, the obtained first-order formulas need not equate the vector field to zero anymore (as in the second line the equation of ϕ 9 in Sect. 3.3). Instead they exclusively contain positivity conditions on all variables and conditions n−1 = 0, n−2 > 0, …, 1 > 0 on Hurwitz determinants. Again, all variables are existentially 25 In some of the cited publications convex coordinates are called reaction coordinates.
Here we obtain a collection of 21 polynomial equations of the form n−1 = 0. 29 The largest polynomial n−1 = 0 in that collection has 10 variables with degrees between 5 and 12 and around 863,000 summands. 30 We discovered a suitable zero in 5 cases, including that largest polynomial, and we could exclude the existence of a zero in all remaining 16 cases.
The overall computation time was 16 s, 15 s of which were spent on the above-mentioned largest polynomial, where we found a negative and a positive point in 10 s and computed the zero in another 5 s. To give an impression of the speed of our methods, we mention that reading and parsing that polynomial from a file takes 25 s, and, e.g., determining the degrees of its 10 variables would take 2 s.
In the positive case, our subtropical approach is capable of finding many suitable zeros in a model: On the one hand, supp + ( f ) ∩ newton( f ) typically contains more than one point. On the other hand, with each point in supp + ( f ) ∩ newton( f ) we can follow the corresponding moment curve to obtain infinitely many solutions. Unfortunately, with MAPK we did not succeed in finding any zero that also satisfies the necessary Hopf conditions n−2 > 0 and n−3 > 0 so far.

Concluding Remarks
The quantifier elimination methods used throughout this article offer the expressivity and flexibility of first-order logic. On the other hand, these methods are deeply rooted in symbolic computation, specifically in real algebra. This allows to make use of the rich algorithmic infrastructure coming with state-of-the-art computer algebra systems. We have seen applications of our methods ranging from geometry via verification to the natural sciences. It is noteworthy that with all our applications we never used regular real quantifier elimination but always certain variants. However, those variants have not been designed for our particular applications. Instead they belong to a collection of established tools, applicable to many more domains than we have discussed here. As a counterpart to the very general concept of quantifier elimination we have seen quite specialized subtropical methods. The development of those methods was triggered by the success of quantifier elimination-based methods in chemistry and in the life sciences. We believe that, similarly to real quantifier elimination and its variants, subtropical methods have the potential to develop into a way more general tool than it might appear now. An important next step is the subtropical treatment of input with one or several inequalities as side conditions in both theory and software.