Coordinate-independent singular perturbation reduction for systems with three time scales

On the basis of recent work by Cardin and Teixeira on ordinary differential equations with more than two time scales, we devise a coordinate-independent reduction for systems with three time scales; thus no a priori separation of variables into fast, slow etc. is required. Moreover we consider arbitrary parameter dependent systems and extend earlier work on Tikhonov-Fenichel parameter values -- i.e. parameter values from which singularly perturbed systems emanate upon small perturbations -- to the three time-scale setting. We apply our results to two standard systems from biochemistry.


Introduction and overview
Ordinary differential equations involving a small parameter appear frequently in mathematics and in science. Their principal use in chemistry and biochemistry -which is the main topic of the present paper -is to find certain (attracting) invariant sets and to achieve reduction of dimension. The mathematical basis is singular perturbation theory, originally due to Tikhonov [18] and Fenichel [4], for systems with one small parameter ε (or, in other words, for systems with two time scales). While Tikhonov's and Fenichel's theory is concerned with first order approximations in ε, there exist approaches to include higher order terms in ε, e.g. to improve accuracy in the approximation of invariant manifolds; see for instance the critical survey in Kaper and Kaper [10]. More recently, Noel et al. [13], Radulescu et al. [15], Samal et al. [16,17] developed an algorithmic method to compute slow-fast scenarios in chemical reaction networks, using tropical geometry. Concerning the existence (or persistence) of invariant sets obtained by such (a priori formal) calculations one may invoke hyperbolicity properties; for instance Theorem 4.1 in Chicone [1] is very useful in this respect. A direct method for chemical reaction networks involving different orders of a single small parameter, given certain properties of the system, is due to Cappeletti and Wiuf [2]. A different perspective is the consideration of systems with more than two time scales by introducing, cum grano salis, several small parameters ε 1 , ε 2 , . . ., and to obtain invariant manifolds and reduction on this basis. (One has the option to set all parameters equal in the end.) Recently Cardin and Teixeira [3] generalized Fenichel's fundamental theorems, proving results on invariant sets and reductions of systems with more than two time scales. Here, the differential equation systems are assumed to have variables separated into blocks of fast, slow, "very slow" ones, and so on. The present paper is based, on the one hand, on Cardin and Teixeira [3]. On the other hand, we extend earlier work [7,8] that is concerned with coordinateindependent reduction (not requiring an a priori separation of slow and fast variables), as well as with the basic question of finding -in arbitrary parameter dependent systems -critical parameter values from which singular perturbation reductions emanate. We will focus on the three time-scale setting, essentially to keep notation manageable, and will only briefly sketch extensions to more than three time scales. Furthermore we will mostly consider systems that satisfy not only the normal hyperbolicity conditions from [3] but have the stronger feature of exponential attractivity. One reason for this restriction lies in our interest in chemical reaction networks. But beyond this practical consideration, the algorithms to compute critical parameter values for singular perturbation scenarios indeed requires this additional property. The paper is organized as follows. In Section 2 we review the work by Cardin and Teixeira [3]. Section 3 generalizes the coordinate-independent reduction algorithm from [7] to three-timescale systems. In Section 4 we start from a general parameter dependent system and extend the work from [8] on critical parameter values (Tikhonov-Fenichel parameter values) to three time scales (resp. two "small parameters"), and in Section 5 we discuss two classical examples (cooperativity with two complexes, competitive inhibition) from biochemistry in detail. Section 6 contains a few remarks about more than three time scales, and finally, for the reader's convenience, we prove some essentially known facts in an Appendix.

Separated fast and slow variables
In this section we review and specialize results from Cardin and Teixeira [3] for a parameter dependent ordinary differential equation system ; brieflyẋ = f (x, ε 1 , ε 2 ).
Here x = (x 1 , x 2 , x 3 ) tr ∈ R n with x 1 ∈ R n1 , x 2 ∈ R n2 , and x 3 ∈ R n3 , and f is smooth on an open neighborhood of U × [0, δ 1 ) × [0, δ 2 ), with U ⊆ R n open and nonempty, and δ 1 > 0, δ 2 > 0. We define and we will assume throughout that these sets are nonempty. Cardin and Teixeira require some hyperbolicity conditions, which we state here in slightly stronger versions, for the sake of simplicity: • First hyperbolicity condition: For every x ∈ M 1 , all the eigenvalues of D x1 f 1 (x, 0, 0) 1 have nonzero real parts. For sufficiently small ε 1 , ε 2 this condition implies local solvability of the implicit equation f 1 (x, ε 1 , ε 2 ) = 0 in the form x 1 = g(x 2 , x 3 , ε 1 , ε 2 ), and one may furthermore write • Second hyperbolicity condition: For every x ∈ M 2 , all the eigenvalues of D x2 f 2 (x, 0, 0) have nonzero real parts. 2 By suitable choice of U , δ 1 and δ 2 we may assume that M 1 and M 2 are submanifolds. Continuing to follow [3] we introduce the auxiliary system and the intermediate reduced system In both equations (4) and (5) above the dot denotes differentiation with respect to τ 2 := ε 1 t. By suitable choice of δ 2 we may also assume that every M ε2 2 is a submanifold of R n . Finally we define the completely reduced system on M 2 , where the dot in (6) denotes differentiation with respect to τ 3 := ε 1 ε 2 t.
We replace the hyperbolicity conditions from [3] by stronger requirements, since in our applications we focus on attracting invariant manifolds. Definition 1. We say that system (1) satisfies the hyperbolic attractivity condition (HA) if D x1 f 1 (x, 0, 0) has only eigenvalues with negative real part on M 1 and if furthermore D x2 f 2 (x, 0, 0) has only eigenvalues with negative real part on M 2 .
Our starting point is the following theorem, specialized from Cardin and Teixeira [3], Theorems A, B and Corollary A. Some of our statements are informal; for rigorous statements and pertinent definitions we refer to [3]. Theorem 1. Let system (1) be given, with (HA) satisfied.
(a) Let N ⊆ M 2 be a compact submanifold (with nonempty interior in the relative topology, and possibly with boundary). Then for all sufficiently small ε 1 , ε 2 there exists a locally invariant manifold N ε1,ε2 for system (1) which is O(ε 1 + ε 2 ) close to N , diffeomorphic to N and locally exponentially attracting. Given the appropriate time scales, solutions of (1) on N ε1,ε2 converge to solutions of (6) on N .
(b) Let ε 2 be sufficiently small and let L ⊆ M ε2 2 be a compact submanifold (with nonempty interior in the relative topology, and possibly with boundary). Then for all sufficiently small ε 1 there exists a locally invariant manifold L ε1,ε2 for system (1) which is O(ε 1 + ε 2 ) close to L, diffeomorphic to L and locally exponentially attracting. Given the appropriate time scales, solutions of (1) on L ε1,ε2 converge to solutions of (5) on L.
As given, the part regarding f 2 in condition (HA) is not ready to use in applications. We provide two equivalent versions. Proposition 1. Condition (HA) is equivalent to either of the following conditions.
(i) D x1 f 1 (x, 0, 0) has only eigenvalues with negative real parts on M 1 , and has only eigenvalues with negative real parts on M 2 .
(ii) D x1 f 1 (x, 0, 0) has only eigenvalues with negative real parts on M 1 , and for all sufficiently small ε > 0 the matrix has only eigenvalues with negative real parts on M 2 .
Proof. We use the notions introduced with the hyperbolicity condition (H) and Definition 1. From one gets by the chain rule , and a further application of the chain rule shows the equivalence of (HA) and (i). The equivalence of (i) and (ii) follows from Lemma 3 in the Appendix.
Remark 1. (a) One may rewrite systems (1) through (6) to some extent, with no effect on the reductions. Using Hadamard's lemma, one may restate (1) asẋ with only the f i remaining in the subsequent reductions. Thus the auxiliary system becomes and there are analogous modifications for the intermediate and the fully reduced system.
(b) The passage from (1) to the completely reduced system (6) can evidently be obtained in the following manner: Fix ε 1 > 0 and reduce (1) with respect to the small parameter ε 2 (in time scale ε 2 t). Then let ε 1 → 0, rescaling time once more to τ 3 . We will use this observation later on.

Coordinate-free reduction
In the present section we generalize the coordinate-independent reduction procedure from [14,7] to the three-timescale setting. The first task is to intrinsically characterize those systems which admit a transformation to "standard form" (1). Reversing matters, applying a (local) smooth coordinate transformation to equation (1) yields a smooth system , evidently satisfying the following conditions: (i) For all sufficiently small ε 1 ≥ 0, ε 2 ≥ 0, the zeros of g (0,0) (x, ε 1 , ε 2 ) form a submanifold M 1 ⊆ U , of codimension n 1 , 1 ≤ n 1 < n. Given any compact submanifold P 1 ⊆ M 1 , there exists θ 1 > 0 such that at every y ∈ P 1 the derivative D x g(y, ε 1 , ε 2 ) admits the eigenvalue zero with algebraic and geometric multiplicity n − n 1 , and the remaining eigenvalues have real parts ≤ −θ 1 .
By Remark 1 one may assume that system (7) is in the special form adjusting conditions (i) and (ii) accordingly. Conditions (i) and (ii) are certainly necessary for (7) or (8) to be a transformed version of (1). The first part of the next lemma shows sufficiency. (b) Condition (i) for system (8) is equivalent to the following: For any y ∈ M 1 there exist a neighborhood U 1,y , a smooth map P 1 : U 1,y → R n×n1 such that P 1 (y, ε 2 ) has rank n 1 , and a smooth map µ 1 : U 1,y → R n1 such that D x µ 1 (y, ε 2 ) has rank n 1 , yielding a decomposition and moreover there is a θ 1 > 0 such that has only eigenvalues with real part ≤ −θ 1 , for all x ∈ U 1,y .
Proof. The nontrivial assertion of part (a) follows from the existence of n − n 1 independent first integrals of g (0,0) in a neighborhood of y, which was noted by Fenichel [4], Lemma 5.3 for smooth vector fields, and shown in [14], Proposition 2.2 for the analytic setting, and likewise from the existence of n − n 1 − n 2 independent first integrals of g (0,0) + ε 1 g (1,0) in a neighborhood of y. These first integrals determine slow and "very slow" variables. Parts (b) and (c) are straightforward applications of [7], Theorem 1, Remark 4 and Remark 2.
Remark 2. The existence of the decomposition g (0,0) = P 1 µ 1 in part (b) (as well as the decomposition in part (c)) is a consequence of the implicit function theorem in the smooth or analytic case. For polynomial or rational vector fields there exists a decomposition with rational functions as entries of P 1 and µ 1 , and there is an algorithmic approach to its computation. See [7] for details.
Next we use the decompositions to compute reductions.
Proposition 2. (a) In arbitrary coordinates the reduction corresponding to the passage from system (1) to the auxiliary system may be obtained as follows: Given ε 2 ≥ 0, determine the projection matrix The auxiliary system (4) for ε 2 > 0 then corresponds tȯ on the local invariant manifold defined by µ 1 (x, ε 2 ) = 0. The equation corresponding to the intermediate reduced system (5) is obtained by setting ε 2 = 0.
(b) In arbitrary coordinates the reduction corresponding to the passage from system (1) to the completely reduced system (6) may be obtained as follows: Given ε 1 > 0, determine the projection matrix Then Q 2 (x, ε 1 ) extends smoothly to a matrix valued function Q 2 (x) at ε 1 = 0. The equation corresponding to the completely reduced system in arbitrary coordinates is given byẋ Proof. Part (a) is a direct application of [7], Theorem 1. For part (b) this theorem is also applicable, but there is a technical problem involving Q 2 as along the image, for x ∈ M 2 (see [7], Remark 1). With the conditions given in Lemma 1 (c) the image is equal to the column space of (P 1 , ε 1 P 2 ), which in turn equals the column space W 1 of (P 1 , P 2 ). The latter matrix has full rank at ε 1 = 0, and its entries depend smoothly on ε 1 and x. Moreover the kernel is equal to the kernel of (D x µ 1 , D x µ 2 ) tr , and we may assume w.l.o.g. that with entries depending smoothly on ε 1 . Thus there remains to verify that the matrix of the projection onto W 2 along W 1 depends smoothly on ε 1 . For the sake of completeness we give a proof of this fact in Lemma 4, Appendix.
We note that the reduction also works, including convergence properties, under the weaker assumption corresponding to (H) rather than (AH) for A 1 and A 2 in Lemma 1.

Remark 3.
While Proposition 2 provides the reduced equations, one also needs initial values for these, which may be obtained from an initial value y of system 8 with the help of the first integrals noted in the proof of Lemma 1(a); see [7], Proposition 2: • Assuming that y is sufficiently close to M 1 , the corresponding initial value (up to an error of order ε 1 + ε 2 ) for the auxiliary system and for the intermediate reduced system is the (locally unique) intersection of M 1 and the level sets of n − n 1 independent first integrals ofẋ = g (0,0) (x, 0) which contain y.
• Assuming that y is sufficiently close to M 2 , the corresponding initial value (up to an error of order ε 1 + ε 2 ) for the auxiliary system and for the intermediate reduced system is the (locally unique) intersection of M 2 and the level sets of n − n 1 − n 2 independent common first integrals ofẋ = This is equivalent to the condition stated.) To illustrate the procedure with an example, we recall the competitive inhibition network with substrate S, enzyme E, inhibitor I and two complexes C 1 , C 2 ; see for instance Keener and Sneyd [11]. The reaction scheme is given by which leads (with the usual assumptions of mass action kinetics, spatial homogeneity and constant thermodynamical parameters) to the differential equation system for the concentrations. (The original system is five dimensional; the two linear first integrals e + c 1 + c 2 and i + c 2 yield reduction to dimension three.) Example 1. In system (9) set x = (s, c 1 , c 2 ) tr and assume (Coloquially speaking, binding to the inhibitor and degradation from the inhibitor complex are slow, while degradation from the substrate complex to enzyme and product is very slow.) This is of the type (8), Moreover M 2 is contained in the common zero set of M 1 is contained in the zero set of µ 1 , and we have We determine the auxiliary system and the intermediate reduced system. With and furthermore Application to Setting ε 2 = 0 one obtains the intermediate reduced system. When the initial values for system (9) are given by (s 0 , c 1,0 , c 2,0 ), to obtain the approximate initial values (s * 0 , c * 1,0 , c * 2,0 ) on M 1 one uses (according to Remark 3) the two first integrals s + c 1 and c 2 of g (0,0) and the defining equation for M 1 , thus the system which leads to quadratic equations for s and c 1 .
To find the fully reduced system one first computes .
The computation of the projection matrix is straightforward (although a software system is helpful) but the output is sizeable. We just record the fully reduced system (in time scale ε 1 ε 2 t). It is given bẏ restricted to the invariant curve M 2 . Finally, given initial values (s 0 , c 1,0 , c 2,0 ) for system (9), approximate initial values for the fully reduced system may be determined by solving the algebraic equations s + c 1 = s 0 + c 1,0 , µ 1 = 0 and µ 2 = 0.

Critical parameter values
Typically in applications one starts with a general parameter dependent system, rather than a system of type (1) or (7) with pre-assigned "small parameters". Therefore the first task is to determine critical parameter values, for which small perturbations lead to singular perturbation scenarios. Thus we consider Tikhonov-Fenichel parameter values, as defined in [8] for two time scales, and extend the notion to the three time scale setting.
• The systemẋ = h(x, π) admits s independent local analytic first integrals at x 0 .
All the conditions in the lemma can be represented by polynomial equations and inequalities. The condition on the roots of χ(τ, x 0 , π)/τ s is characterized by inequalities: There exist n − s Hurwitz determinants (see e.g. Gantmacher [5], Ch. V, §6, Thm. 4 ff.) which must attain values > 0. Moreover, the existence requirement for s independent first integrals leads to a series of polynomial equations via degree by degree evaluation of Taylor expansions. More precisely, for every d > 0 there is an induced action of D x h(x, π) on the space S 1 + · · · + S d of polynomials in n variables with zero constant term and of degree ≤ d. Extending condition (i), the characteristic polynomial of this action (which coincides with (11) for d = 1) must have vanishing coefficients for all sufficiently small powers of the indeterminate. (No further inequalities appear, due to the structure of the eigenvalues for this action.) Thanks to Hilbert's Basissatz, finitely many of these equations suffice. A full account is given in [8].
For the remainder of this section we assume that Π ⊆ R m + is a semi-algebraic set, and that system (10) admits the positively invariant subset R n + . Then, as was shown in [8], the Tikhonov-Fenichel parameter values for dimension s, 1 ≤ s < n form a semi-algebraic subset Π s ⊆ R m . We will denote the Zariski closure of Π s by W s . Thus the elements of W s satisfy all defining equations for Π s but not necessarily the defining inequalities.
We note some properties of nested TFPV.
Proposition 3. (a) Any TFPV π ∈ Π s+k which is nested in Π s lies in the boundary of Π s relative to its Zariski closure W s .
Remark 4. Proposition 3 opens a starting point for the computation of nested TFPV: Start with system (10) corresponding to "generic" parameter values in Π s , i.e. parameter values in the intersection of Π s with an irreducible component of the Zariski closure W s . In order to find nested parameters for higher dimension one only needs to look at the boundary of Π s , and one can use part (b) in order to obtain necessary conditions. Practically this may be realized by determining the decomposition P ·µ for generic π ∈ Π s and then looking at zeros of Dµ · P , with parameters in the boundary. (The boundary may also contain further parameter values in Π s .)

Special settings for chemical reaction networks
For chemical reaction networks (CRN) the parameter region is usually given by and for many such systems and given s, the irreducible components of W s are just determined by the vanishing of certain of the π i ; see e.g. [6,8,9]. (The underlying reason for this fact is the subject of forthcoming work.) Thus we have, for π in a given irreducible component: (i) Upon relabelling, there is an ℓ, 0 < ℓ < m such that π i = 0 for all i ∈ {ℓ + 1, · · · m}; (ii) the remaining parameters are nonnegative.
In other words, the intersection of Π s with the given irreducible component of W s corresponds to some subset of R ℓ + , with boundary R ℓ + \ R ℓ + . This leads to an obvious case-by-case analysis. Note that boundary points may or may not be contained in Π s , but there is no loss in starting with "generic" parameter values in the interior R ℓ >0 . We look at a particular example. Example 2. We again consider competitive inhibition; see equation (9). Here the parameters are of the form From [8], Proposition 8 we have the necessary condition e 0 k 1 k 2 k −3 = 0 for a TFPV in Π 1 , with each of the four cases (e.g. e 0 = 0 and all other parameter values ≥ 0) yielding a singular perturbation reduction with attracting one dimensional critical manifold. Hence W 1 has four irreducible components. In order to find nested TFPV's for dimension 2 we perform a case-by-case investigation. We only consider one case here; see Section 5 for the remaining ones. For the case k 2 = 0 the system is given bẏ where we have used the abbreviations e = e 0 − c 1 − c 2 and i = i 0 − c 2 ; note that e ≥ 0 and i ≥ 0 by design of (9). By Remark 4, nested TFPV's for dimension two and corresponding points in the critical manifold necesarily satisfy det (Dµ · P ) = 0, with Proceeding according to Remark 4, we determine the vanishing set of Since all the variables and parameters are nonnegative, this sum equals zero if and only if every summand vanishes. In particular, k −1 · k −3 has to vanish for any nested TFPV. We look at the two ensuing cases.
(i) k −1 = 0: Then the remaining conditions are If k 1 = 0 or k 3 = k −3 = 0 we obtain a two dimensional variety of stationary points. Checking the attractivity conditions (HA), one finds that these cases yield nested TFPV. If e = s = 0 holds then we get e 0 −c 1 −c 2 = 0 = s which corresponds to a one dimensional variety. In case e = k −3 = 0 we get c 1 = 0 while c 2 and s are arbitrary, thus we have a two dimensional (attracting) variety of stationary points.
(ii) k −3 = 0: In this case there remains In view of case (i) we only have to check k 3 = 0 or e = i = 0. In both cases we get a variety of dimension two.

Further examples
In this section we continue the discussion of the competitive inhibitor network, to some extent, and furthermore present a fairly complete investigation of a cooperative system with two complexes, following the strategy outlined in Remark 4. Missing from a complete analysis are some cases concerned with boundary points in Π 1 ⊆ W 1 which themselves belong to Π 1 , as well as certain degenerate cases for Π 2 . Moreover we will not generally record routine calculations to verify conditions such as (HA), and for ease of notation we will frequently use the term "critical manifold" for the Zariski closure of this object, without mentioning the inequalities to be satisfied.

Competitive inhibition (cont.)
We continue to investigate the competitive inhibition network; see equation (9), Examples 1 and 2. The analysis of TFPV which was started in Example 2 will be finished here. For Π 1 there are three remaining cases, viz. e 0 = 0, k 1 = 0 and k −3 = 0.
(a) For e 0 = 0, the system is given by We only consider the generic case for Π 1 , thus all the remaining parameters are > 0. Then the (only possible) decomposition P · µ for the right hand side is given by For nested TFPV, a simple computation yields the necessary condition with all terms positive; thus every summand must vanish, and in particular In case k −1 = k 2 = 0 system (13) admits a two dimensional variety of stationary points, with k 1 s = 0 only if c 1 + c 2 = 0. The intersection of this variety with the positive orthant is only one dimensional, thus we do not obtain a two dimensional critical manifold. The cases with k 1 s = 0 translate to k 1 = εk * 1 for the system with small parameters. (Otherwise the critical manifold would be given by s = 0, which does not contain the line given by c 1 = c 2 = 0.) Moreover we have e 0 = ε 1 ε 2 e * 0 , hence every term k 1 e 0 s is of the form ε 2 1 ε 2 · (· · · ), and the completely reduced system is necessarily trivial. Likewise, the case k −3 = 0 in system (13) leads to k 1 s = 0. To summarize, the case e 0 = 0 yields no interesting reductions for the three time scale setting, in marked contrast to the familiar (quasi-steady state) reduction for small initial enzyme concentration with two time scales.
(b) Next we consider the system with k 1 = 0, i.e.
Because of Example 2 we may assume that k 2 = 0, which yields c 1 = 0 for stationary points. In case k −3 = k 3 = 0 we indeed have a two dimensional critical manifold. Turning to small parameters we have k 1 = ε 1 ε 2 k * 1 , k 3 = ε 1 k * 3 and k −3 = ε 1 k * −3 , and (9) becomes We compute the reductions for this case. For the auxiliary system (on the critical variety defined by c 1 = 0) we obtain the decomposition  , and a straightforward computation yields the projection matrix on the variety defined by c 1 = 0. The intermediate reduced system is obtained setting ε 2 = 0. Turning to the complete reduction, the decomposition of the "fast part" of (14) is given by  One obtains and may continue as prescribed by Proposition 2. There are shortcuts, though: First note that the critical manifold is given by c 1 = 0, and c 2 constant and equal to the smaller solution c 2 of the quadratic equation The completely reduced system will automatically yieldċ 1 = 0 andċ 2 = 0, hence only the first row of the projection matrix needs to be computed. As the final result of the reduction procedure we get the equatioṅ with the dot denoting differentiation with respect to ε 1 ε 2 t.
(c) Finally, we deal with the case k −3 = 0, which does not automatically yield a one dimensional variety of stationary points. System (9) becomeṡ We may assume that k 1 = 0 and k 2 = 0, otherwise one would arrive at (non-generic) subcases of previously discussed systems. From this we obtain c 1 = 0 and es = 0 as necessary conditions. Now e = 0 and nonnegativity of varaibles imply e 0 = 0; a previously discussed case, therefore every stationary point satisfies s = 0. If k 3 = 0 then i = 0 forces c 2 = i 0 ; the corresponding parameter values are not in Π 1 . So the only case remaining is k 3 = k −3 = 0 (very slow binding to the inhibitor, very slow degradation of the inhibitor complex), with systeṁ To obtain a two dimensional variety of stationary points one has to check the boundary of Π 1 for nested TFPV, which splits into four cases. We only discuss the case k 1 = 0 here, thus (9) with small parameters becomes The computation of the auxiliary system runs similar to the reduction of (14) and yields  ṡ on the invariant variety given by c 1 = 0. Finally, the completely reduced system lives on the variety defined by c 1 = s = 0 (a coordinate subspace), and therefore by [9], Proposition 5 the reduced system may be directly obtained via "classical" QSS reduction; yieldinġ

A cooperative system
In this subsection we study the standard cooperative system involving substrate S, two complexes C 1 , C 2 , enzyme E and product P . The reaction scheme yields, with the usual assumptions and stoichiometry, the differential equation where all appearing constants are non-negative. According to Goeke [6], Kap. 9.4, necessary conditions for TFPV are given by e 0 k 1 k 2 (k −3 + k 4 ) = 0.

Case k 1 = 0
When we substitute k 1 = 0 in equation (16) we obtaiṅ Hence considering the generic case (all remaining parameters > 0) we obtain an irreducible component of W 1 given by k 1 = 0, and the critical manifold is given by c 1 = c 2 = 0. We get a decomposition with and necessary conditions for nested TFPV from Thie first set of conditions does not, by itself, yield a two dimensional critical manifold, and we will not pursue it further here. The second set, i.e. k −3 = k 4 = 0, yields the two dimensional variety given by c 1 = 0. Considering this setting, we introduce the small parameters in our original system by substituting Ordering the parameters as e 0 , k 1 , k −1 , k 2 , , k 3 , k −3 , k 4 , we thus consider the surface in parameter space given by and with x = (s, c 1 , c 2 ) tr we get (17) h(x, ε 1 , For this system we compute the complete reduction on c 1 = c 2 = 0 and the intermediate reduction on c 1 = 0. In order to compute the completely reduced system, a factorization of g (0,0) + ε 1 g (1,0) is given by with µ 1 = c 1 , µ 2 = c 2 , and The projection matrix is and the fully reduced system in very slow time on the invariant manifold c 1 = c 2 = 0 is given by the equatioṅ Similarly one computes the intermediate system on the two dimensional variety given by c 1 = 0 from the decomposition P 1 · µ 1 : 

Case e 0 = 0
From e 0 = 0 one also obtains a component of W 1 , and system (16) specializes toṡ The right hand side has a factorization P · µ with and in order to obtain nested TFPV we examine all variable-parameter configurations that satisfy The plane given by s = 0 is not a viable candidate for a two dimensional critical manifold since it does not contain the line c 1 = c 2 = 0. This leaves the cases The first of these yields a two dimensional variety (defined by k 3 sc 1 −k −3 c 2 = 0) only under the additional condition k 4 = 0. The second case, whenever k 1 = 0, yields a variety whose intersection with the positive orthant has dimension one, hence is of no relevance. For the third case we obtain a two dimensional variety only if k 1 = 0 or k 2 = 0.
With the exception of this very last case, the completely reduced system will always be trivial, due to k 1 = ε 1 k * 1 and e 0 = ε 1 ε 2 e * 0 , which implies k 1 e 0 = O(ε 2 1 ε 2 ). We consider one spacial case, viz. the intermediate reduction coresponding to the nested TFPV with k 1 = k 3 = k −3 = k 4 = 0; here c 1 = 0 defines the two dimensional critical manifold. Considering we compute:

The intermediate reduced system on the invariant variety
(b) In case k 4 = 0 the remaining system iṡ Adding the first two equations for stationary points shows that k −3 c 2 = 0, and combining this with the third equation yields k 3 sc 1 = 0; in addition one has k 1 es − k −1 c 1 = 0. Thus there are further conditions for the existence of a one dimensional critical variety.
(i) Given that k 3 = 0 and k −3 = 0, the variety is given either by c 2 = s = 0, which yields the parameter conditions or the variety is given by c 1 = c 2 = 0, with parameter conditions For all these parameters the next task is to discuss conditions for embedded TFPV. We will only do so for two cases.
1. In the case k −1 = k 2 = 0 one has a decomposition  We take a closer look at the case k 1 = k 3 = 0, with critical variety c 2 = 0. The surface in parameter space The intermediate reduced system is as follows: and the completely reduced system (on s = c 2 = 0) is given by: 2. In the case k 2 = k 3 = k 4 = 0 we have the one dimensional critical manifold k 1 s · (e 0 − c 1 ) − k −1 c 1 = 0, c 2 = 0 and the right hand side of the system at k 2 = k 3 = k 4 = 0 can be decomposed into P · µ, with This yields det(Dµ · P ) = (k 1 (e 0 − c 1 ) + k 1 s + k −4 ) · k −3 .
Proof. We suppress the arguments (x, ε) in the notation. By assumption C := B 1 , B 2 is invertible, and the entries of C −1 depend smoothly on (x, ε). With the projection matrix given by the assertion is obvious.