TIKHONOV’S THEOREM AND QUASI-STEADY STATE

. There exists a systematic approach to asymptotic properties for quasi-steady state phenomena via the classical theory of Tikhonov and Fenichel. This observation allows, on the one hand, to settle convergence issues, which are far from trivial in asymptotic expansions. On the other hand, even if one takes convergence for granted, the approach yields a natural way to compute a reduced system on the slow manifold, with a reduced equation that is fre- quently simpler than the one obtained by the ad hoc approach. In particular, the reduced system is always rational. The paper includes a discussion of nec- essary and suﬃcient conditions for applicability of Tikhonov’s and Fenichel’s theorems, computational issues and a direct determination of the reduced sys- tem. The results are applied to several relevant examples.

1. Introduction. The mathematical description and analysis of reacting systems in chemistry and biochemistry frequently leads to a slow-fast separation for the associated differential equations. From a mathematical perspective, methods and results from singular perturbation theory should provide a natural approach. But, quite remarkably, these are rarely used in a systematic manner (albeit frequently invoked as a guiding principle) in the discussion of such differential equations. The principal purpose of the present paper is to show that, and how, a systematic employment of the standard theorems is possible in the context of chemistry and biochemistry. In particular there is a straightforward procedure to compute the reduced system corresponding to Tikhonov's theorem [26]. The most important basis for our results is Fenichel's classical paper [7].
We will consider homogeneously mixed chemically reacting systems with constant thermodynamic parameters, and will therefore discuss (parameter-dependent) ordinary differential equation systems. We will assume these to be derived from mass action kinetics; thus we will deal with polynomial systems. One may distinguish two common incarnations of slow-fast separation: • Slow and fast reactions (or rather, pairs of forward-backward reactions); • Slow and fast variables (QSS; quasi-steady state phenomena).
Informally speaking, one hopes for the fast reactions or variables to run their course quickly, leaving a reduced system moving on some "slow" manifold. In much of the chemistry literature (see e.g. Atkins [1] or Stryer [2]), for QSS the following ad hoc reduction procedure is employed: Set the rates of change for the quasistationary variables equal to zero (thus certain entries of the right-hand side of the differential equation are set equal to zero), and use the ensuing algebraic equations to obtain a system of smaller dimension.
From a mathematical perspective, such reductions pose a nontrivial task. Singular perturbation theory, in particular Tikhonov's [26] and Fenichel's [7] theorems (see the lecture notes by Jones in [11], Verhulst's monograph [29], and also Hoppensteadt's extension [10] to infinite time intervals) is frequently chosen as a framework. For instance, Kaper and Kaper [12] discuss various reduction methods and asymptotic expansions in the light of Fenichel's invariant manifold theory. But a rigorous application of singular perturbation theory is not commonly found when QSS or general slow-fast phenomena are discussed for particular chemically reacting systems. There is a notable exception: For the setting of slow and fast reactions, Schauer and Heinrich [21] provided a systematic discussion, including a transfer to the setting of Tikhonov's theorem, based on Vasil'eva [28]. (Stiefenhofer [25] followed up on the work of Schauer and Heinrich and considered a more general setting.) But in the QSS case, a systematic approach is not widely used for a transformation to standard scenarios of singular perturbation theory. It seems that only Duchêne and Rouchon [6] clearly identified the connection with and the transfer to Tikhonov. On the other hand, Segel and Slemrod, in their seminal paper [23] on the Michaelis-Menten reaction, included a convergence proof but proceeded directly without reference to existing theory.
One may identify three problems concerning the reduction of equations for reacting system, independent of the approach one chooses. To start with, one frequently first needs to identify a small parameter from the system and perhaps from additional hypotheses. This is straightforward for slow and fast reactions (see Schauer and Heinrich [21]), but more delicate for slow and fast variables. Segel and Slemrod [23] employed time scale arguments, as did many others (see Borghans et al. [4], Murray [16], Schnell and Maini [22]). This first problem is discussed in greater detail by the authors in [19], where a different heuristics is proposed. The second problem concerns the range of validity of the asymptotic reduction and convergence questions. In particular, this problem is not settled for the ad hoc reduction procedure used in chemistry. The third problem (assuming a positive resolution of the first two) is the computation of a reduced equation. Except for relatively simple settings, the ad hoc approach leads to complicated nonlinear and parameter-dependent systems of algebraic equations, and the reduced equations may be rather unpleasant. As a remark in Keener and Sneyd ([13], Section 1.2) indicates, reversible reactions are frequently replaced by irreversible ones to simplify the analysis, although from a physiological point of view such a simplification may not be justified. In particular, reversible versions of the reductions for a number of standard reaction equations lead to rather cumbersome equations, if they can be obtained at all.
Fundamental results of singular perturbation theory, such as Tikhonov's theorem, help settle the convergence problem if they are applicable. (We will always refer to the version of Tikhonov's theorem given in Verhulst [29], Thm. 8.1. We will focus on the lowest-order approximation, including convergence questions; higherorder expansions will not be discussed.) But there is a nontrivial obstacle to the application of Tikhonov to parameter-dependent differential equations for chemically reacting systems: Generally, the differential equation is not given in the normal form required for straightforward application of this theorem, and a priori it may be uncertain whether a transformation to such a normal form exists. We will provide necessary, and locally sufficient, criteria for transformability to "Tikhonov normal form", thus taking a necessary step towards application of the theory to specific problems. The crucial question is concerned with the existence of first integrals for the degenerate system (at parameter value zero). In this sense the approach via Tikhonov's and Fenichel's theory may be more technically involved than other methods.
But on the other hand, and perhaps surprisingly, it is possible to directly compute the reduced system corresponding to Tikhonov's theorem, and this computation involves only basic algebraic operations in the case of differential equations with polynomial, or rational, right-hand side. (Essentially this observation is due to Fenichel [7]; see also Duchêne and Rouchon [6].) Thus, Tikhonov and Fenichel provide a straightforward solution to the third problem noted above; this does not seem to have been noticed widely.
We discuss a number of familiar examples for QSS, and derive reduced systems. While in some cases the reduction via the ad hoc procedure used in chemistry leads to the same result as the procedure corresponding to Tikhonov, the latter generally yields reduced systems that are different from (and actually less cumbersome than) those obtained by the standard approach. Convergence problems will also be settled for a number of these examples; this indicates another advantage of our approach.
The article is based on the first author's doctoral dissertation [17], and should be seen as the second in a series beginning with [19].
2. Singular perturbation theory and QSS. Let us introduce the setting and notation first: A differential equation for a chemically reacting system is given in the formẋ with a small parameter ε ≥ 0. (Both n and m are positive integers.) We will assume mass action kinetics; in particular h is analytic in x and ε. There may be additional parameters that h depends on; these will usually be suppressed in the notation. The subset U of R n+m is assumed to have nonempty interior; frequently U will be a compact set, and h will be defined in some neighborhood of U × [0, ε 0 ], for some positive ε 0 .
2.1. Transformation to a normal form. A principal obstacle to the direct application of Tikhonov's and Fenichel's theory lies in the fact that the main theorems are stated for systems in what we will call Tikhonov normal form: One may also view this system in "slow time" by setting τ = ε · t. Then In our setting f and g will be analytic. For easy reference, we recall Tikhonov's theorem, as given in Verhulst [29], Thm. 8.1 (slightly paraphrased and specialized).
(c) The initial value y 2 (0) is contained in an interior subset of the domain of attraction of y 2 , for parameter value y 1,0 . Given these hypotheses, there exists a τ 1 > 0 such that for every τ 0 > 0 the solution of the initial value problem (3) converges on [τ 0 , τ 1 ] to the solution of (4) with given initial value, as → 0.
Verification of the hypotheses for Tikhonov's theorem is facilitated by an additional condition on the eigenvalues of the derivative of g (see e.g. Fenichel [7]), viz.
Moreover we will assume that U contains points with the property g(y 1 , y 2 , 0) = 0. The implicit function theorem then guarantees that g(y 1 , y 2 , 0) = 0 defines a submanifold of (some neighborhood of) D × G, with dimension n. (This defines the integer n in equation (2).) Together with the criterion for asymptotic stability by linearization, we then see that hypotheses (a) and (b) are satisfied. Hypothesis (c) may hold only in some neighborhood of y 2 .
The existence question for a transformation to Tikhonov normal form, and the determination of such a transformation, is a central problem here. Fenichel noted this, as well as the existence of obstacles (see [7], Lemma 5.1 ff.); see also Duchêne and Rouchon [6], Sections 1 and 2. Thus assume that for every small ε ≥ 0 there exists a system of the form (2) and a continuously differentiable map Φ = Φ(x, ε) for x in some neighborhood of U , which is a diffeomorphism in x and sends the solutions of (1) to solutions of (2). A necessary and locally sufficient condition for this property is the identity We abbreviate Ψ(x) := Φ(x, 0), and denote the inverse of Ψ by Γ. For ε = 0 equation (6) then reduces to The existence of an invertible Ψ satisfying identity (7) is necessary and sufficient for the existence of a transformation to a system as desired: Indeed, one may define Φ(x, ε) := Ψ(x); then the higher order terms of f and g are uniquely determined by h. In any case, the lowest-order terms εf (y 1 , y 2 , 0) and g(y 1 , y 2 , 0) are determined by Φ(·, 0) and h alone. Choosing a transformation independent of ε also brings the advantage that it works both for slow and fast time scales.
Proposition 2.1. Given system (1), there exists a solution-preserving map Φ to (2) with the additional property (5) only if the following hold: There exists µ > 0 such that for every x 0 in the zero set M 0 of h(·, 0), the derivative Dh(x 0 , 0) admits the eigenvalue 0 with geometric multiplicity n, and the remaining eigenvalues have real part ≤ −µ. In particular M 0 is an ndimensional submanifold.
Proof. To prove part (a), evaluate the identity (7). As for part (b), differentiate this identity to see that the matrices Dh(x 0 , 0) and are conjugate when h(x 0 , 0) = 0. The implicit function theorem shows the remaining assertion.
This result indicates how to determine a map Ψ such that (7) holds: Find n independent first integrals for h(·, 0) and extend these to a local diffeomorphism. From a computational point of view, finding first integrals explicitly may pose a problem. Moreover, even local existence of nonconstant first integrals in a neighborhood of a stationary point is, a priori, a nontrivial requirement. But the existence problem can be settled.
Proposition 2.2. Assume that condition (b) in Proposition 2.1 is satisfied. Then for every x 0 ∈ M 0 there exists a neighborhood on which n independent analytic first integrals of h(·, 0) are defined, and locally there exists an analytic diffeomorphism Ψ such that (6) is satisfied.
Proof. In this proof we require analyticity. By condition (b) in Proposition 2.1, the Jacobian at a stationary point of h(·, 0) has eigenvalue 0 with multiplicity n, and the remaining eigenvalues have negative real part. According to Bibikov [3], Theorem 10.1 and Theorem 12.2 (and its proof) there exists an invertible analytic transformation Ψ of the system to a quasi-normal form (QNF) of the forṁ The crucial point in Bibikov's argument is convergence of such a transformation: Generally, for every degree there exists an analytic transformation to partial quasinormal formẏ with the r i, of degree at most and . . . standing for terms of higher degree. Because the system admits an n-dimensional manifold of stationary points, one may conclude that all r 1, = 0. This implies the existence of a formal transformation to a formal power series QNF of the form (8), and this, in turn, implies the existence of a convergent transformation. System (8) obviously admits n linear first integrals.
Fenichel also noticed this fact and sketched a proof of Proposition 2.2 using a geometric version of the argument; see [7], Lemma 5.3. (There seems to be no discussion of such matters in Stiefenhofer [25].) Variants of the proof given above have been used before; see e.g. Theorem 1 in [5]. Bibikov's approach also provides an algorithmic method to determine power series expansions for the first integrals, up to any desired degree. However, such approximations may be of little interest for practical purposes.
Remark 1. In the setting of slow and fast reactions one may frequently find sufficiently many linear first integrals for the degenerate system from stoichiometry. This was -at least implicitly -noted and used by Schauer and Heinrich [21]. Thus in the slow and fast reaction scenario an explicit transformation to Tikhonov normal form may be easy to find. For the slow and fast variable scenario the issue seems more complicated.

2.2.
Reduction. There remains to determine a reduced system corresponding to (1). Given the expansion one may use the condition h (0) (x) = 0 to obtain the "slow manifold". Buṫ is not the reduced differential equation: As noted by Stiefenhofer [25], one has to modify this equation to achieve invariance of the zero set M 0 of h (0) . However, there are many possible modifications. To determine the reduction based on (2), assume that the transformation Φ = Ψ is independent of ε, and start with system (4). While this system cannot be transported back directly, one can embed the differential equation into a system on R n+m such that g(y 1 , y 2 , 0) = 0 defines an invariant set, and both systems are equivalent on the zero set of g(·, 0): Consider and g(y 1 , y 2 , 0) = 0. (Actually, the entries of g(·, ·, 0) are first integrals of this system.) Tikhonov's convergence results remain valid -with obvious modifications -for this reduced system: According to [29], Thm. 8.1 (see Subsection 2.1) there exists τ 1 > 0 (independent of ε) such that on every interval [τ 0 , τ 1 ] with τ 0 > 0 the solution (z 1,ε , z 1,ε ) of (3) converges to the solution (z 1 , z 2 ) of (4) with corresponding initial value for z 1 . Due to this differential equation, the invertibility condition on D 2 g, and the implicit function theorem, (z 1 , z 2 ) is differentiable, and solves the first equation of (10) due to (4). The second equation of (10) follows from differentiating g(z 1 , z 2 ) = 0 with respect to τ . Proposition 2.3. A reduced system for (1) corresponding to Tikhonov's theorem is given by with p from equation (10), up to rescaling time. More precisely, this means that the zero set M 0 of g(·, ·, 0) is an invariant set for x = q(x) and that there exists t 1 > 0 such that for each 0 < t 0 < t 1 , every solution of the time scaled version of (1) starting in a suitable neighborhood of M 0 converges to a solution of (11) for τ ∈ [t 0 , t 1 ] as ε → 0.
Proof. By construction, Γ sends solutions of y = p(y) to solutions of x = q(x), and solutions of (2) to solutions of (1).
From Fenichel [7], Lemma 5.4 and the arguments in its proof, one may obtain a direct way to compute a reduced system. Fenichel's argument starts from the observation that, for every x ∈ M 0 , R n+m is the direct sum of the kernel and the image of Dh (0) (x), because algebraic and geometric multiplicity of the eigenvalue 0 coincide. We will take a different path to determine this decomposition; the proof is constructive.
Lemma 2.4. Assume that condition (b) in Proposition 2.1 is satisfied for every zero of h (0) . Then for each x ∈ M 0 there exists a polynomial σ x (τ ) in the indeterminate τ such that σ x (0) = 0 and that τ ·σ x (τ ) annihilates the linear map Dh (0) (x). Moreover there exists a polynomial α x such that For every v ∈ R n+m the projection onto the kernel of Dh (0) (x) which corresponds to the image-kernel direct sum decomposition of R n+m is given by Proof. We use some standard (linear) algebra; see for instance Lang [14], Ch. XV. The characteristic polynomial of T := Dh (0) (x) is of the form (suppressing the subscripts x) with σ(0) = 0. By the Hamilton-Cayley theorem, and because algebraic and geometric multiplicity of the eigenvalue 0 coincide, the minimum polynomial of T divides τ · σ(τ ). There exists a polynomial α such that (in the present case this is obvious), and for every v ∈ R n+m the decomposition corresponds to the image-kernel direct sum decomposition of R n+m with respect to T .
Proposition 2.5. Assume that condition (b) in Proposition 2.1 is satisfied for every zero of h (0) . For each x ∈ M 0 let the polynomial σ x (τ ) be as defined in Lemma 2.4. Then a reduced system corresponding to Tikhonov's theorem is given by Proof. Let Ψ denote a local transformation of h to Tikhonov normal form. As in the proof of Proposition 2.1, one has that for every x 0 ∈ M 0 the matrices Dh (0) (x 0 ) and are conjugate by DΨ(x 0 ). For this system in Tikhonov normal form a reduced equation is given by (10).
One may see this from a different perspective: Letting the reduced equation from (10) corresponds precisely to the projection of f (1) g (1) onto the kernel of T with respect to the image-kernel decomposition. The general assertion now follows from properties of solution-preserving maps and the fact that conjugation of linear maps carries over to the kernel-image decomposition.
We emphasize that this result is basically due to Fenichel; the crucial argument can be found in [7], Section V. Our approach to compute the projection works in more general settings than Fenichel's Lemma 5.4 which, like Stiefenhofer's [25] expression for the reduced equation, requires an explicit equation for the slow manifold. Duchêne and Rouchon [6] take a different path, which yields the same reduced equation up to higher-order terms in ε. With this approach and a numerical approximation of the slow manifold M 0 they determine and numerically integrate a reduced system of a high-dimensional system from combustion kinetics.
If one starts with a vector field h that is polynomial or rational in x then, remarkably, the reduced system (12) will always be rational. The algebraic variety M 0 will automatically be invariant for this system.

Determining transformations in practice.
From a practical perspective, one needs to explore to what extent explicit knowledge of the maps Ψ and Γ is necessary. Again we refer to the version of Tikhonov's theorem given in Verhulst [29]. One has to check whether there exists a transformation to Tikhonov normal form and whether conditions (a), (b) and (c) in [29], Thm. 8.1 are satisfied.
Thus assume that we are given system (1), and that condition (b) in Proposition 2.1 is satisfied for each zero x 0 of h(·, 0). (The latter can be verified algorithmically, due to the Hurwitz-Routh criterion.) Then Proposition 2.2 shows the local existence of a transformation to Tikhonov normal form, and furthermore conditions (a) and (b) in [29], Thm. 8.1 are satisfied. Moreover, one need not know Ψ or Γ explicitly to find M 0 , or to find the reduced system via Proposition 2.5.
Condition (c), on the other hand, may restrict the position of the initial value. This initial value is often fixed by the experimental setting, hence local properties may not suffice, and explicit knowledge of Ψ and its inverse Γ may become necessary. Moreover, in order to verify Hoppensteadt's [10] conditions for the extension of Proposition 2.3 to infinite time intervals, knowledge of the Tikhonov normal form, and thus of Ψ and Γ, seems to be necessary.
This poses the question of explicitly determining first integrals for the degenerate system in the QSS scenario. We are unable to provide a general solution, but the following observation is of some use. Lemma 2.6. If the reacting system described by (1) consists only of first and second order reactions, with mass action kinetics, and one has m = 1, then the first integrals ofẋ = h(x, 0) can be determined from first integrals of a differential equation with right hand side at most linear.
Proof. In this case the polynomial systemẋ = h(x, 0) is quadratic and there exists a one-codimensional submanifold of zeros of h(·, 0). Hence there exist a scalar polynomial function ρ and a polynomial vector field h * of degree < 2 such that h(·, 0) = ρ · h * . The first integrals of h * are first integrals of h(·, 0), and vice versa.
3. Application: The Michaelis-Menten system. In this section we apply our results to the well-known Michaelis-Menten system, see [15]. This system models an enzyme-catalyzed reaction where substrate S and enzyme E combine to form a complex C, which degrades back to substrate and enzyme, or to product P and enzyme. Initially, no complex or product are present. The four-dimensional differential equation for the concentrations admits two linear first integrals and therefore a two-dimensional problem remains. As in any QSS scenario, one first has to identify a suitable "small parameter", but this was taken care of by Segel and Slemrod [23], and in [18]; see also the general heuristics to determine "small parameters" proposed by the authors in [19]. The reader may notice that the particular choice of an "appropriate small parameter" ε seems not very important in the following examples, and that replacing it by C · ε for some positive constant C does not affect the procedure. This is due to the fact that only the limit ε → 0 is being considered here. But the quality of the approximation, and the range where the approximation works, will strongly depend on a good parameter choice.
We will first discuss the reversible Michaelis-Menten reaction with QSS assumption for complex. While the irreversible system is covered in almost any introductory biochemistry text, the reversible system is not discussed as extensively in the literature. Moreover we will consider the reverse QSS assumption for substrate, in the case of irreversible product formation. For both settings it turns out that the ad hoc procedure used in the chemistry literature does not lead to the reduced system found via singular perturbation theory, and that the latter actually are less complicated.
In the examples the variables are chosen and ordered to make the transfer to Tikhonov normal form as convenient as possible. To make the procedure more accessible to less theory-inclined readers, some requisite computations are presented in a relatively detailed manner.
3.1. The reversible system. We first discuss the Michaelis-Menten reaction with reversible product formation, given the classical quasi-steady state assumption for the complex. A convenient version of the differential equation iṡ with initial values c(0) = 0 and s(0) = s 0 . The k i are rate constants, all of which will be assumed > 0, and s 0 resp. e 0 denote initial concentrations for substrate, resp. enzyme. If one considers QSS for the complex C then the small parameter determined by Segel and Slemrod [23] is given by The discussion in [19] confirms that this is an appropriate choice when k −2 is small. Using the abbreviations one obtains the equivalent version

LENA NOETHEN AND SEBASTIAN WALCHER
The nonzero eigenvalue of Dh (0) (0, s) is equal to due to 0 ≤ s ≤ s 0 , and therefore Tikhonov is applicable. To determine a first integral for the system with ε = 0, one only needs to consider the linear systeṁ We follow the procedure outlined in Gröbner and Knapp [9], Ch. III, Section 2. A preliminary step is to transform to a homogeneous linear system. The stationary point is given by For c * := c −ĉ and s * := s −ŝ, one obtains From the obvious identities denotes the Lie derivative of a scalarvalued function with respect to A) one obtains a first integral Choosing ψ 2 := s * and returning to the original coordinates one has the transformation with inverse Γ(y 1 , y 2 ) = −y 2 + y 1/k1 1 y k−2/k1 2 +ĉ y 2 +ŝ .
We next determine g in the Tikhonov normal form (2). One computes where the first entry is equal to zero by construction, and the second is obvious. Expressing c and s via Γ, one finds g(y 1 , y 2 , 0) = ĉ − y 2 + y To verify the conditions for Tikhonov's theorem, consider g as a function of y 2 with parameter y 1 . The only zero of g in the interesting interval y 2 +ŝ > 0 (which is the only relevant stationary point forẏ 2 = g(y 1 , y 2 , 0)) is given by the unique root of the first factor. Since the nonzero eigenvalue of Dh (0) (0, s) is ≤ −(k −1 + k 2 ), this stationary point is linearly asymptotically stable, which implies global asymptotic stability in dimension one. To summarize, Tikhonov's theorem as given in [29] is applicable.
To compute the reduced system, we use Proposition 2.5. The stationary points of h (0) in the positive orthant are given by c = 0, and one has The polynomial σ according to Lemma 2.4 is obviously given by It is straightforward to find the reduced system on the invariant set given by c = 0: In addition toċ = 0, one obtainṡ This version is different from (and actually less complicated than) the one from the ad hoc reduction procedure. Indeed the ad hoc method, starting from "ċ = 0", or rather from yields a quadratic equation for c as a function of s, and a rather cumbersome reduced differential equation for s(t) for which, in turn, various simplifications have been suggested; see Fraser and Roussel [8], Schauer and Heinrich [20], Seshadri and Fritzsch [24], Tzafriri and Edelman [27]. The approach taken here provides more satisfactory results. Determining the reduced system via the Tikhonov normal form and Proposition 2.3 will yield the same result, with a higher computational effort. But this extra effort provides an additional benefit: One can verify Hoppensteadt's convergence conditions for the transformed system, although working out the details requires some tenacity. (The crucial point which facilitates the analysis is that the normalization following Hoppensteadt's Condition (III) is almost automatic.) Therefore we have convergence on every time interval [t 0 , ∞) with t 0 > 0.
As for the classical Michaelis-Menten system with irreversible product formation, thus k −2 = 0, a variant is convenient. The system for complex and substrate can be written as follows: The orbit equation for h (0) . With ψ 2 (c, s) := s one obtains a transformation Ψ. From this point on the procedure runs analogous to the reversible case. One obtains the classical reduced system, which is just (15) with k −2 = 0.
This seems to be the first instance that the approach via Tikhonov has been employed for the asymptotic reduction of the Michaelis-Menten system with quasistationary complex.
3.2. Reverse QSS. The reverse QSS assumption (RQSSA) for the irreversible Michaelis-Menten reaction was discussed briefly by Segel and Slemrod [23], and more extensively by Borghans et al. [4]. In this variant, substrate S is assumed to be in quasi-steady state. A convenient version of the Michaelis-Menten equations for further discussion is given by: with p(0) = 0 and s(0) = s 0 . Segel and Slemrod [23] proposed the small parameter ε * := k −1 k 1 e 0 to ensure validity of RQSSA. One sees that setting k −1 = 0 in (16) yields a system with only one stationary point, hence Tikhonov cannot be applied. In this case the appropriate choice and interpretation of a small parameter is relevant: From the experimental perspective ε → 0 corresponds to very high e 0 in the reactor. This observation and the arguments in Borghans et al. [4], as well as those in [18], suggest the choice A preliminary scaling seems appropriate to avoid "infinitely large" terms on the right hand side of the equation. Setting This system is already in Tikhonov normal form, with g(p, s, 0) = −(k 1 s 0 + k −1 )s.
Viewing p as a parameter, the differential equation ds dT = g(p, s, 0) admits the unique stationary point s = 0, which is globally asymptotically stable, with eigenvalue −(k 1 s 0 + k −1 ) for the linearization. The reduced system (in fast time) is therefore given by dp dT = εk 2 (s 0 − p) or, returning to the original time scale, The ad hoc reduction method used in chemistry would lead to a different reduced equation (with no a priori guarantee of convergence): Setting the right hand side of the second entry in (16) equal to zero, solving the ensuing quadratic equation for s and substituting in the first line yields a one-dimensional equation with a more complicated right-hand side. See the discussion in Schnell and Maini [22]. For reverse QSS it is straightforward to verify the conditions given by Hoppensteadt in [10], pp. 522 -523: The constant solution (s 0 , 0) of (17) guarantees that Hoppensteadt's condition (I) holds. Condition (II) is unproblematic, and condition (III), as well as the normalization that g = 0 when ε = 0 and s = 0, hold automatically. The continuity and boundedness conditions (IV) and (V) are unproblematic, and the crucial conditions (VI) and (VII) pertain to uniform asymptotic stability of the degenerate system dp dτ = k 2 (s 0 − p) and of the boundary layer system which are both obvious. Therefore one has convergence on every interval [t 0 , ∞) with t 0 > 0. The discrepancy between our reduced system and the one found by the standard procedure is of order o(ε), and in this particular case one can show convergence for the standard reduction by invoking our convergence result and direct estimates. But the point is that a convergence proof is built in for the Tikhonov approach. 4. A cooperative system. In this section we will discuss a well-known threedimensional system which models a cooperative enzyme reaction, taken from Keener and Sneyd [13]. No explicit first integrals seem to be available here for the system at ε = 0, therefore our results are only local. Within this limitation the theory is readily applicable.
The cooperative reaction (see [13], Subsection 1.2.4) involves enzyme, substrate, product and two complexes C 1 and C 2 . As in the Michaelis-Menten reaction, substrate and enzyme combine reversibly to C 1 , which may degrade irreversibly to enzyme and product. Moreover C 1 and substrate may combine to a second complex C 2 , which in turn may degrade irreversibly to product and C 1 . Taking stoichiometry into account, the following three-dimensional system describes this cooperative enzyme-catalyzed reaction: The relevant initial conditions are The standard choice for the small parameter used by Keener and Sneyd [13] and also by Murray, [16], Section 5.3, is ε = e 0 K , with K := s 0 .
According to the arguments in [19], when one is interested in QSS over the full course of the reaction then will be a more appropriate choice. For the consideration of the limit → 0 this will not make a difference. Rewriting the system with ε = e 0 /K yields: Thus and on the positive orthant we have Moreover, from The conditions for Propositions 2.2 and 2.5 are satisfied. To see this, note that σ(0) > 0, and that the minimum of σ is located at k * := −(k −1 + k 2 + k 1 s + k 3 s + k −3 + k 4 )/2 < 0 (thus, in any case, the real parts of the roots are negative), and that s ranges in the compact interval [0, s 0 ]. According to Lemma 2.4 one has to determine and a straightforward computation yields the result Employing one computes the reduced system, according to Proposition 2.5: Changing time scales, one finds the reduced equation: This result is the same as found by the ad hoc procedure, see Keener and Sneyd [13], Eq. (1.42).
But matters are different when one considers reversible product formation; thus product and enzyme may react to form complex C 1 with rate constant k −2 , and product and C 1 may react to form complex C 2 with rate constant k −4 .
In this fully reversible setting the ad hoc approach (solveċ 1 = 0 andċ 2 = 0 for c 1 and c 2 as functions of s) leads to a system of two coupled parameter-dependent quadratic equations for c 1 and c 2 . Using standard algorithmic methods (e.g. Groebner bases), this system will provide a quartic parameter-dependent equation for each of the unknowns c 1 and c 2 . While in principle there still exist the Cardano formulas to solve these equations, the resulting one-dimensional differential equation becomes practically untractable.
The approach via Proposition 2.5, on the other hand, involves larger expressions compared to the irreversible setting, but they are rational. The basic procedure remains the same as in the irreversible case, with straightforward computations, albeit writing them all down would amount to a serious waste of space, and using a (standard) computer algebra system is convenient. The reduced equation turns out to beṡ = −e 0p (s)/q(s) with p(s) = s 2 (k −3 k −4 + k 3 k 4 )(k 1 − k −2 ) + s ((k 1 k 2 + k −1 k −2 )(k −3 + k 4 ) + k −3 k −4 (2k −2 − k 1 )s 0 + k −2 k 3 k 4 s 0 ) − k −1 k −2 (k −3 + k 4 )s 0 − k −2 k −3 k −4 s 2 0 andq (s) = s 2 ((k 1 − k −2 )(k 3 − k −4 )) + s (( This opens a path to discussing the effects of reversible product formation in this cooperative system, with the help of a relatively simple differential equation. If k 1 > k −2 and k 3 > k −4 then the polynomialq will attain only positive values for s ≥ 0, whilep will have exactly one positive root, corresponding to an asymptotically stable equilibrium. As one would hope, this reduced system turns into the one found for the irreversible case when one sets k −2 = k −4 = 0.
5. Concluding remarks. The principal purpose of this paper is to present a systematic approach to the mathematical analysis of QSS via Tikhonov's and Fenichel's theory. To a considerable extent the problem can be reduced to computational issues: Propositions 2.1 and 2.2 provide necessary and locally sufficient criteria, and the computation of the reduced system from Proposition 2.5 leaves only algebraic problems.
The examples in sections 3 and 4 were chosen in view of their relevance, and of the necessity to discuss the new approach for some benchmark problems. The examples may be somewhat misleading in one respect: Generally the zero set of h (0) will form an algebraic variety more complicated than a linear subspace, and the algebra becomes more involved. But this observation, in conjunction with Proposition 2.5, opens up an interesting and relevant field for applied algorithmic algebra. In future papers we will present and discuss reductions of a number of reaction equations, including reversible versions of several standard systems.
The examples do point toward problems posed by global issues; the main focus here is on the explicit determination of first integrals for degenerate systems. There may exist some theoretical basis for the explicit computation of such first integrals, since the degenerate systems are not arbitrary but are derived from a full system for chemical reactions by a well-defined procedure. The problem remains unresolved here, but it will be interesting to pursue it further.