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

TL;DR

This paper introduces a coordinate-independent reduction method for systems with three time scales, extending existing theories to arbitrary parameter-dependent systems and applying the approach to biochemical models.

Contribution

It develops a novel coordinate-independent reduction technique for three time-scale systems and extends Tikhonov-Fenichel parameter analysis to this setting.

Findings

Applicable to systems with three time scales.

Extended Tikhonov-Fenichel parameter analysis.

Demonstrated on biochemical systems.

Abstract

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.

Equations

$$\begin{array}[]{rcl}\dot{x}_{1}&=&\phantom{\varepsilon_{1}}\phantom{\varepsilon_{2}}f_{1}(x,\varepsilon_{1},\varepsilon_{2})\\ \dot{x}_{2}&=&\phantom{\varepsilon_{1}}\varepsilon_{1}f_{2}(x,\varepsilon_{1},\varepsilon_{2})\\ \dot{x}_{3}&=&\varepsilon_{1}\varepsilon_{2}f_{3}(x,\varepsilon_{1},\varepsilon_{2})\\ \end{array};\quad\text{briefly }\dot{x}=f(x,\varepsilon_{1},\varepsilon_{2}).$$

$${\mathcal{M}}_1:=\left\{x\in U;f_1(x,0,0)=0\right\}$$

$${\mathcal{M}}_2:=\left\{x\in U;f_1(x,0,0)=f_2(x,0,0)=0\right\},$$

$$f_2(x,{\epsilon }_1,{\epsilon }_2)={\tilde{f}}_2(x_2,x_3,{\epsilon }_1,{\epsilon }_2):=f_2(g(x_2,x_3,{\epsilon }_1,{\epsilon }_2),x_2,x_3,{\epsilon }_1,{\epsilon }_2).$$

$$\begin{array}[]{rcl}0&=&\phantom{\varepsilon_{2}}f_{1}(x,0,\varepsilon_{2})\\ \dot{x}_{2}&=&\phantom{\varepsilon_{1}}f_{2}(x,0,\varepsilon_{2})\\ \dot{x}_{3}&=&\varepsilon_{2}f_{3}(x,0,\varepsilon_{2})\\ \end{array}$$

$${\mathcal{M}}_2^{{\epsilon }_2}:=\left\{x\in U;f_1(x,0,{\epsilon }_2)=0\right\},$$

$$\begin{array}[]{rcl}0&=&f_{1}(x,0,0)\\ \dot{x}_{2}&=&f_{2}(x,0,0)\\ \dot{x}_{3}&=&0\\ \end{array}$$

$$\begin{array}[]{rcl}0&=&f_{1}(x,0,0)\\ 0&=&f_{2}(x,0,0)\\ \dot{x}_{3}&=&f_{3}(x,0,0)\\ \end{array}$$

$$B_1(x):=-D_{x_1}{f}_2(x,0,0)D_{x_1}{f}_1(x,0,0{)}^{-1}{D}_{x_2}{f}_1(x,0,0)+D_{x_2}{f}_2(x,0,0)$$

$$B_2(x,\epsilon ):=\left(\begin{array}{cc}{D}_{x_1}{f}_1(x,0,0)& D_{x_2}{f}_1(x,0,0)\\ \epsilon D_{x_1}{f}_2(x,0,0)& \epsilon D_{x_2}{f}_2(x,0,0)\end{array}\right)$$

$$f_1(g(x_2,x_3,{\epsilon }_1,{\epsilon }_2),x_2,x_3,{\epsilon }_1,{\epsilon }_2)=0$$

$$D_{x_2}g(x_2,x_3)=-D_{x_1}{f}_1(x,,{\epsilon }_1,{\epsilon }_2{)}^{-1}{D}_{x_2}{f}_1(x,,{\epsilon }_1,{\epsilon }_2)$$

$$\begin{array}[]{rclll}\dot{x}_{1}&=&\widehat{f}_{1}(x,\varepsilon_{2})&+\varepsilon_{1}\widehat{f}_{1,1}(x,\varepsilon_{1})&+\varepsilon_{1}\varepsilon_{2}\widehat{f}_{1,2}(x,\varepsilon_{1},\varepsilon_{2})\\ \dot{x}_{2}&=&&\phantom{+}\varepsilon_{1}\widehat{f}_{2}(x,\varepsilon_{1})&+\varepsilon_{1}\varepsilon_{2}\widehat{f}_{2,2}(x,\varepsilon_{1},\varepsilon_{2})\\ \dot{x}_{3}&=&&&\phantom{+}\varepsilon_{1}\varepsilon_{2}\widehat{f}_{3}(x,\varepsilon_{1},\varepsilon_{2})\\ \end{array}$$

$$\begin{array}[]{rcl}0&=&\phantom{\varepsilon_{2}}\widehat{f}_{1}(x,\varepsilon_{2})\\ \dot{x}_{2}&=&\phantom{\varepsilon_{1}}\widehat{f}_{2}(x,0)\\ \dot{x}_{3}&=&\varepsilon_{2}f_{3}(x,0,\varepsilon_{2})\\ \end{array}$$

$$\dot{x}=g^{(0,0)}(x,{\epsilon }_1,{\epsilon }_2)+{\epsilon }_1\left(g^{(1,0)}(x,{\epsilon }_1,{\epsilon }_2)+{\epsilon }_2{g}^{(1,1)}(x,{\epsilon }_1,{\epsilon }_2)\right)$$

$$g^{(0,0)}(x,{\epsilon }_1,{\epsilon }_2)+{\epsilon }_1{g}^{(1,0)}(x,{\epsilon }_1,{\epsilon }_2)$$

$$\dot{x}=g^{(0,0)}(x,{\epsilon }_2)+{\epsilon }_1\left(g^{(1,0)}(x,{\epsilon }_1)+{\epsilon }_2{g}^{(1,1)}(x)\right)+O({\epsilon }_2({\epsilon }_1+{\epsilon }_2)),$$

$$g^{(0,0)}(x,{\epsilon }_2)=P_1(x,{\epsilon }_2){\mu }_1(x,{\epsilon }_2),$$

$$A_1(x,{\epsilon }_2):=D{\mu }_1(x,{\epsilon }_2)P_1(x,{\epsilon }_2)$$

$$g^{(0,0)}(x,{\epsilon }_2)+{\epsilon }_1{g}^{(1,0)}(x,{\epsilon }_1)=P_1(x,{\epsilon }_2){\mu }_1(x,{\epsilon }_2)+{\epsilon }_1{P}_2(x,{\epsilon }_1){\mu }_2(x,{\epsilon }_1),$$

$$A_2(x,{\epsilon }_1,{\epsilon }_2):=\left(\begin{array}{c}D{\mu }_1(x,{\epsilon }_2)\\ D{\mu }_2(x,{\epsilon }_1)\end{array}\right)\left(\begin{array}{cc}{P}_1(x,{\epsilon }_2)& {\epsilon }_1{P}_2(x,{\epsilon }_1)\end{array}\right)$$

$$Q_1(x,{\epsilon }_2):=I_n-P_1(x,{\epsilon }_2)A_1(x,{\epsilon }_2{)}^{-1}{D}_x{\mu }_1(x,{\epsilon }_2).$$

$$\dot{x}=Q_1(x,{\epsilon }_2)\left(g^{(1,0)}(x,0)+{\epsilon }_2{g}^{(1,1)}(x)\right)$$

$${\tilde{Q}}_2(x,{\epsilon }_1):=I_n-\left(\begin{array}{cc}{P}_1(x,0),& {\epsilon }_1{P}_2(x,0)\end{array}\right)A_2(x,{\epsilon }_1,0{)}^{-1}\left(\begin{array}{c}{D}_x{\mu }_1(x,0)\\ D_x{\mu }_2(x,0)\end{array}\right).$$

$$\dot{x}=Q_2(x)g^{(1,1)}(x)$$

$$D_x{g}^{(0,0)}(x,0)+{\epsilon }_1{D}_x{g}^{(1,0)}(x,{\epsilon }_1)$$

$$\left(\begin{array}{c}{D}_x{\mu }_1\\ D_x{\mu }_2\end{array}\right)=\left(\begin{array}{cc}{A}_1& A_2\end{array}\right)$$

$$\left(\begin{array}{c}-A_1^{-1}{A}_2\\ I\end{array}\right)$$

$$\begin{array}[]{rcccl}E+S&\overset{k_{1}}{\underset{k_{-1}}{\rightleftharpoons}}&C_{1}&\overset{k_{2}}{\rightharpoonup}&E+P,\\ E+I&\overset{k_{3}}{\underset{k_{-3}}{\rightleftharpoons}}&C_{2}&&\end{array}$$

$$\begin{array}[]{rcl}\dot{s}&=&k_{-1}c_{1}-k_{1}s(e_{0}-c_{1}-c_{2})\\ \dot{c}_{1}&=&k_{1}s(e_{0}-c_{1}-c_{2})-(k_{-1}+k_{2})c_{1}\\ \dot{c}_{2}&=&k_{3}(e_{0}-c_{1}-c_{2})(i_{0}-c_{2})-k_{-3}c_{2}\end{array}$$

Full text

MSC (2010): 92C45, 34E15, 37D10

Key words: Reaction network, dimension reduction, invariant set, multiple time scales.

1 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 $\varepsilon$ (or, in other words, for systems with two time scales).

While Tikhonov’s and Fenichel’s theory is concerned with first order approximations in $\varepsilon$, there exist approaches to include higher order terms in $\varepsilon$, 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 $\varepsilon_{1},\,\varepsilon_{2},\ldots$, 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 coordinate-independent 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.

2 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

[TABLE]

Here $x=(x_{1},x_{2},x_{3})^{\rm tr}\in\mathbb{R}^{n}$ with $x_{1}\in\mathbb{R}^{n_{1}}$, $x_{2}\in\mathbb{R}^{n_{2}}$, and $x_{3}\in\mathbb{R}^{n_{3}}$, and $f$ is smooth on an open neighborhood of $U\times[0,\delta_{1})\times[0,\delta_{2})$, with $U\subseteq\mathbb{R}^{n}$ open and nonempty, and $\delta_{1}>0$, $\delta_{2}>0$.

We define

[TABLE]

and

[TABLE]

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\in{\cal M}_{1}$, all the eigenvalues of $D_{x_{1}}f_{1}(x,0,0)$111For a smooth function $g=g(x,y,\ldots)$ we denote the partial derivatives by $D_{x}g$, $D_{y}g$ etc. have nonzero real parts.

For sufficiently small $\varepsilon_{1},\varepsilon_{2}$ this condition implies local solvability of the implicit equation $f_{1}(x,\varepsilon_{1},\varepsilon_{2})=0$ in the form $x_{1}=g(x_{2},x_{3},\varepsilon_{1},\varepsilon_{2})$, and one may furthermore write

[TABLE]

• •

Second hyperbolicity condition: For every $x\in{\cal M}_{2}$, all the eigenvalues of $D_{x_{2}}\widetilde{f}_{2}(x,0,0)$ have nonzero real parts.222In [3] the second hyperbolicity condition is erroneously written for $f_{2}$ rather than $\widetilde{f}_{2}$. The authors are aware of this and will publish a corrigendum.

By suitable choice of $U$, $\delta_{1}$ and $\delta_{2}$ we may assume that ${\cal M}_{1}$ and ${\cal M}_{2}$ are submanifolds.

Continuing to follow [3] we introduce the auxiliary system

[TABLE]

on

[TABLE]

and the intermediate reduced system

[TABLE]

on ${\cal M}_{1}$. In both equations (4) and (5) above the dot denotes differentiation with respect to $\tau_{2}:=\varepsilon_{1}t$. By suitable choice of $\delta_{2}$ we may also assume that every ${\cal M}_{2}^{\varepsilon_{2}}$ is a submanifold of $\mathbb{R}^{n}$.

Finally we define the completely reduced system

[TABLE]

on ${\cal M}_{2}$, where the dot in (6) denotes differentiation with respect to $\tau_{3}:=\varepsilon_{1}\varepsilon_{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_{x_{1}}f_{1}(x,0,0)$ has only eigenvalues with negative real part on $\mathcal{M}_{1}$ and if furthermore $D_{x_{2}}\widetilde{f}_{2}(x,0,0)$ has only eigenvalues with negative real part on $\mathcal{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.

1. (a)

Let $\mathcal{N}\subseteq\mathcal{M}_{2}$ be a compact submanifold (with nonempty interior in the relative topology, and possibly with boundary). Then for all sufficiently small $\varepsilon_{1},\varepsilon_{2}$ there exists a locally invariant manifold $\mathcal{N}_{\varepsilon_{1},\varepsilon_{2}}$ for system (1) which is $O(\varepsilon_{1}+\varepsilon_{2})$ close to $\mathcal{N}$, diffeomorphic to $\mathcal{N}$ and locally exponentially attracting. Given the appropriate time scales, solutions of (1) on $\mathcal{N}_{\varepsilon_{1},\varepsilon_{2}}$ converge to solutions of (6) on $\mathcal{N}$. 2. (b)

Let $\varepsilon_{2}$ be sufficiently small and let $\mathcal{L}\subseteq\mathcal{M}_{2}^{\varepsilon_{2}}$ be a compact submanifold (with nonempty interior in the relative topology, and possibly with boundary). Then for all sufficiently small $\varepsilon_{1}$ there exists a locally invariant manifold $\mathcal{L}_{\varepsilon_{1},\varepsilon_{2}}$ for system (1) which is $O(\varepsilon_{1}+\varepsilon_{2})$ close to $\mathcal{L}$, diffeomorphic to $\mathcal{L}$ and locally exponentially attracting. Given the appropriate time scales, solutions of (1) on $\mathcal{L}_{\varepsilon_{1},\varepsilon_{2}}$ converge to solutions of (5) on $\mathcal{L}$.

As given, the part regarding $\widetilde{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.

1. (i)

$D_{x_{1}}f_{1}(x,0,0)$* has only eigenvalues with negative real parts on $\mathcal{M}_{1}$, and*

[TABLE]

has only eigenvalues with negative real parts on $\mathcal{M}_{2}$. 2. (ii)

$D_{x_{1}}f_{1}(x,0,0)$* has only eigenvalues with negative real parts on $\mathcal{M}_{1}$, and for all sufficiently small $\varepsilon>0$ the matrix*

[TABLE]

has only eigenvalues with negative real parts on $\mathcal{M}_{2}$.

Proof.

We use the notions introduced with the hyperbolicity condition (H) and Definition 1. From

[TABLE]

one gets by the chain rule

[TABLE]

when $f_{1}(g(x_{2},x_{3},\varepsilon_{1},\varepsilon_{2}),x_{2},x_{3},\varepsilon_{1},\varepsilon_{2})=0$, 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**.**

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

[TABLE]

with only the $\widehat{f}_{i}$ remaining in the subsequent reductions. Thus the auxiliary system becomes

[TABLE]

and there are analogous modifications for the intermediate and the fully reduced system. 2. (b)

The passage from (1) to the completely reduced system (6) can evidently be obtained in the following manner: Fix $\varepsilon_{1}>0$ and reduce (1) with respect to the small parameter $\varepsilon_{2}$ (in time scale $\varepsilon_{2}t$). Then let $\varepsilon_{1}\to 0$, rescaling time once more to $\tau_{3}$. We will use this observation later on.

3 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

[TABLE]

on an open neighborhood of $\widetilde{U}\times[0,\delta_{1})\times[0,\delta_{2})\subseteq\mathbb{R}^{n}\times\mathbb{R}\times\mathbb{R}$ ($\widetilde{U}\subseteq\mathbb{R}^{n}$ open), evidently satisfying the following conditions:

1. (i)

For all sufficiently small $\varepsilon_{1}\geq 0,\,\varepsilon_{2}\geq 0$, the zeros of $g^{(0,0)}(x,\varepsilon_{1},\varepsilon_{2})$ form a submanifold $\widetilde{\mathcal{M}}_{1}\subseteq\widetilde{U}$, of codimension $n_{1}$, $1\leq n_{1}<n$. Given any compact submanifold $\mathcal{P}_{1}\subseteq\widetilde{\mathcal{M}}_{1}$, there exists $\theta_{1}>0$ such that at every $y\in\mathcal{P}_{1}$ the derivative $D_{x}g(y,\varepsilon_{1},\varepsilon_{2})$ admits the eigenvalue zero with algebraic and geometric multiplicity $n-n_{1}$, and the remaining eigenvalues have real parts $\leq-\theta_{1}$. 2. (ii)

For all sufficiently small $\varepsilon_{1}>0,\,\varepsilon_{2}\geq 0$ the zeros of

[TABLE]

form a submanifold $\widetilde{\mathcal{M}_{2}}\subseteq\widetilde{U}$, of codimension $n_{1}+n_{2}$, $1\leq n_{2}<n-n_{1}$. Moreover, for any compact submanifold $\mathcal{P}_{2}\subseteq\widetilde{\mathcal{M}}_{2}$ there exists a $\theta_{2}>0$ with the following property: At every $y\in\mathcal{P}_{2}$ the derivative $D_{x}g^{(0,0)}(y,\varepsilon_{1},\varepsilon_{2})+\varepsilon_{1}D_{x}g^{(1,0)}(y,\varepsilon_{1},\varepsilon_{2})$ admits the eigenvalue zero with algebraic and geometric multiplicity $n-n_{1}-n_{2}$, and the remaining eigenvalues have real parts $\leq-\theta_{2}\varepsilon_{1}$.

By Remark 1 one may assume that system (7) is in the special form

[TABLE]

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.

Lemma 1**.**

1. (a)

There exists a local diffeomorphism transforming system (8) to a system of type (1) with condition (HA) if and only if conditions (i) and (ii) above hold. 2. (b)

Condition (i) for system (8) is equivalent to the following: For any $y\in\widetilde{\mathcal{M}_{1}}$ there exist a neighborhood $U_{1,y}$, a smooth map $P_{1}:\,U_{1,y}\to\mathbb{R}^{n\times n_{1}}$ such that $P_{1}(y,\varepsilon_{2})$ has rank $n_{1}$, and a smooth map $\mu_{1}:\,U_{1,y}\to\mathbb{R}^{n_{1}}$ such that $D_{x}\mu_{1}(y,\varepsilon_{2})$ has rank $n_{1}$, yielding a decomposition

[TABLE]

and moreover there is a $\theta_{1}>0$ such that

[TABLE]

has only eigenvalues with real part $\leq-\theta_{1}$, for all $x\in U_{1,y}$. 3. (c)

In presence of condition (i), condition (ii) for system (8) is equivalent to the following: For every (sufficiently small) $\varepsilon_{1}>0$ and any $y\in\widetilde{\mathcal{M}_{2}}$ there exist a neighborhood $U_{2,y}$, a smooth map $P_{2}:\,U_{2,y}\to\mathbb{R}^{n\times n_{2}}$ such that $\left(P_{1}(y,\varepsilon_{2}),\varepsilon_{1}P_{2}(y,\varepsilon_{1})\right)$ has rank $n_{1}+n_{2}$, and a smooth map $\mu_{2}:\,U_{2,y}\to\mathbb{R}^{n_{2}}$ such that $\left(D_{x}\mu_{1}(y,\varepsilon_{2}),D_{x}\mu_{2}(y,\varepsilon_{2})\right)^{\rm tr}$ has rank $n_{1}+n_{2}$, yielding a decomposition

[TABLE]

and moreover there is a $\theta_{2}>0$ such that

[TABLE]

has only eigenvalues with real part $\leq-\theta_{2}\varepsilon_{1}$, for all $x\in U_{2,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)}+\varepsilon_{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}\,\mu_{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 $\mu_{1}$, and there is an algorithmic approach to its computation. See [7] for details.

Next we use the decompositions to compute reductions.

Proposition 2**.**

1. (a)

In arbitrary coordinates the reduction corresponding to the passage from system (1) to the auxiliary system may be obtained as follows:

Given $\varepsilon_{2}\geq 0$, determine the projection matrix

[TABLE]

The auxiliary system (4) for $\varepsilon_{2}>0$ then corresponds to

[TABLE]

on the local invariant manifold defined by $\mu_{1}(x,\varepsilon_{2})=0$. The equation corresponding to the intermediate reduced system (5) is obtained by setting $\varepsilon_{2}=0$. 2. (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 $\varepsilon_{1}>0$, determine the projection matrix

[TABLE]

Then $\widetilde{Q}_{2}(x,\varepsilon_{1})$ extends smoothly to a matrix valued function $Q_{2}(x)$ at $\varepsilon_{1}=0$. The equation corresponding to the completely reduced system in arbitrary coordinates is given by

[TABLE]

on the local invariant manifold defined by $\mu_{1}(x,0)=\mu_{2}(x,0)=0$.

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 $\widetilde{Q}_{2}$ as $\varepsilon_{1}\to 0$, since $A_{2}(x,0)$ is non-invertible. To resolve this difficulty, recall that $\widetilde{Q}_{2}(x,\varepsilon_{1})$ is the projection map onto the kernel of

[TABLE]

along the image, for $x\in\widetilde{\mathcal{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},\,\varepsilon_{1}P_{2})$, which in turn equals the column space $W_{1}$ of $(P_{1},\,P_{2})$. The latter matrix has full rank at $\varepsilon_{1}=0$, and its entries depend smoothly on $\varepsilon_{1}$ and $x$. Moreover the kernel is equal to the kernel of $(D_{x}\mu_{1},\,D_{x}\mu_{2})^{\rm tr}$, and we may assume w.l.o.g. that

[TABLE]

with invertible $A_{1}$, whence the kernel is equal to the column space $W_{2}$ of the matrix

[TABLE]

with entries depending smoothly on $\varepsilon_{1}$. Thus there remains to verify that the matrix of the projection onto $W_{2}$ along $W_{1}$ depends smoothly on $\varepsilon_{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 $\widetilde{\mathcal{M}}_{1}$ , the corresponding initial value (up to an error of order $\varepsilon_{1}+\varepsilon_{2}$) for the auxiliary system and for the intermediate reduced system is the (locally unique) intersection of $\widetilde{\mathcal{M}}_{1}$ and the level sets of $n-n_{1}$ independent first integrals of $\dot{x}=g^{(0,0)}(x,0)$ which contain $y$.

• •

Assuming that $y$ is sufficiently close to $\widetilde{\mathcal{M}}_{2}$ , the corresponding initial value (up to an error of order $\varepsilon_{1}+\varepsilon_{2}$) for the auxiliary system and for the intermediate reduced system is the (locally unique) intersection of $\widetilde{\mathcal{M}}_{2}$ and the level sets of $n-n_{1}-n_{2}$ independent common first integrals of $\dot{x}=g^{(0,0)}(x,0)$ and $\dot{x}=g^{(1,0)}(x,0,0)$ which contain $y$. (A direct application of Proposition 2 in **[7]** would lead to simultaneous first integrals of $\dot{x}=g^{(0,0)}(x,0)+\varepsilon_{1}g^{(1,0)}(x,\varepsilon_{1},0)$ for all $\varepsilon_{1}$. 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

[TABLE]

which leads (with the usual assumptions of mass action kinetics, spatial homogeneity and constant thermodynamical parameters) to the differential equation system

[TABLE]

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})^{\rm tr}$ and assume $k_{2}=\varepsilon_{1}\varepsilon_{2}k_{2}^{*}$, $k_{3}=\varepsilon_{1}k_{3}^{*}$ and $k_{-3}=\varepsilon_{1}k_{-3}^{*}$. (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), with

[TABLE]

Moreover $\widetilde{\mathcal{M}}_{2}$ is contained in the common zero set of

[TABLE]

$\widetilde{\mathcal{M}}_{1}$ is contained in the zero set of $\mu_{1}$, and we have

[TABLE]

We determine the auxiliary system and the intermediate reduced system. With

[TABLE]

one has

[TABLE]

and furthermore

[TABLE]

Application to

[TABLE]

yields the auxiliary system (in time scale $\varepsilon_{1}t$) on $\widetilde{\mathcal{M}}_{1}$:

[TABLE]

Setting $\varepsilon_{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 $\widetilde{\mathcal{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 $\widetilde{\mathcal{M}_{1}}$, thus the system

[TABLE]

which leads to quadratic equations for $s$ and $c_{1}$.

To find the fully reduced system one first computes

[TABLE]

and

[TABLE]

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 $\varepsilon_{1}\varepsilon_{2}t$). It is given by

[TABLE]

with

[TABLE]

and

[TABLE]

restricted to the invariant curve $\widetilde{\mathcal{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}$, $\mu_{1}=0$ and $\mu_{2}=0$. **

4 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.

4.1 Tikhonov-Fenichel parameter values

Tikhonov-Fenichel parameter values (TFPV) were introduced in [8] for polynomial (or rational) parameter dependent systems

[TABLE]

A TFPV $\widehat{\pi}$ is characterized by the property that small perturbations $\pi=\widehat{\pi}+\varepsilon\rho+\cdots$ along a smooth curve in parameter space $\Pi$ give rise to a singular perturbation reduction for

[TABLE]

with locally exponentially attracting critical manifold. (The definition extends easily to smooth systems but the algorithmic approach relies on the stronger assumption.) There exists an intrinsic characterization of TFPV’s, see [8] Lemmas 1 and 2, for which the characteristic polynomial

[TABLE]

of the Jacobian $D_{x}h(x,\pi)$ is relevant. We recall:

Lemma 2**.**

Given $0<s<n$, a parameter value $\widehat{\pi}$ is a TFPV with locally exponentially attracting critical manifold $Z_{s}$ (depending on $\widehat{\pi}$) of dimension $s$, and $x_{0}\in Z_{s}$, if and only if the following hold:

• •

$h(x_{0},\widehat{\pi})=0$.

• •

The characteristic polynomial $\chi(\tau,x,\pi)$ from (11)) satisfies

1. (i)

$\sigma_{0}(x_{0},\widehat{\pi})=\cdots=\sigma_{s-1}(x_{0},\widehat{\pi})=0$; 2. (ii)

all roots of $\chi(\tau,x_{0},\widehat{\pi})/\tau^{s}$ have negative real parts.

• •

The system $\dot{x}=h(x,\widehat{\pi})$ 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 $\chi(\tau,x_{0},\widehat{\pi})/\tau^{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,\pi)$ on the space $S_{1}+\cdots+S_{d}$ of polynomials in $n$ variables with zero constant term and of degree $\leq 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 $\Pi\subseteq\mathbb{R}_{+}^{m}$ is a semi-algebraic set, and that system (10) admits the positively invariant subset $\mathbb{R}^{n}_{+}$. Then, as was shown in [8], the Tikhonov-Fenichel parameter values for dimension $s$, $1\leq s<n$ form a semi-algebraic subset $\Pi_{s}\subseteq\mathbb{R}^{m}$. We will denote the Zariski closure of $\Pi_{s}$ by $W_{s}$. Thus the elements of $W_{s}$ satisfy all defining equations for $\Pi_{s}$ but not necessarily the defining inequalities.

4.2 Nested Tikhonov-Fenichel parameter values

Generalizing the approach to TFPV in [8]), and taking into account the special form of (7), it seems reasonable to consider surfaces in parameter space. Thus consider a smooth surface of the special form

[TABLE]

defined in some nighborhood of $(0,\,0)$. Substitute $\gamma(\varepsilon_{1},\,\varepsilon_{2})$ for $\pi$ in (10) to get

[TABLE]

with the $g^{(i,j)}$ smooth by Hadamard’s lemma. In order to obtain a system (7) that also satisfies the conditions (i) and (ii) preceding Lemma 1, the following is necessary: There exist $s>0$ and $k>0$ such that $\widehat{\pi}\in\Pi_{s+k}$, and $\widehat{\pi}+\varepsilon_{1}\cdot{\rho_{1}}(\varepsilon_{1})\in\Pi_{s}$ for all sufficiently small $\varepsilon_{1}>0$. (Note that $\varepsilon_{2}$ plays no role in these conditions.) This observation gives rise to:

Definition 2**.**

Given system (10) and $s,\,k>0$ with $s+k<n$, let $\delta>0$ and let

[TABLE]

be a smooth curve such that

1. (i)

$\beta(\varepsilon_{1})\in\Pi_{s}$* for all $\varepsilon_{1}>0$,* 2. (ii)

$\widehat{\pi}:=\beta(0)\in\Pi_{s+k}$.

Then we call $\widehat{\pi}$ a Tikhonov-parameter value (for dimension $s+k$) nested in $\overline{\Pi}_{s}$.

We note some properties of nested TFPV.

Proposition 3**.**

1. (a)

Any TFPV $\widehat{\pi}\in\Pi_{s+k}$ which is nested in $\overline{\Pi}_{s}$ lies in the boundary of $\Pi_{s}$ relative to its Zariski closure $W_{s}$. 2. (b)

Let $\beta$ as in Definition 2, and for $\varepsilon_{1}>0$ consider the decomposition

[TABLE]

according to **[7]**, Theorem 1. Then

[TABLE]

on the critical manifold.

Proof.

Part (a) is a direct consequence of the definition. As for part (b), at $\varepsilon_{1}=0$, with $\widehat{\pi}\in\Pi_{s+k}$ and $x_{0}\in Z_{s+k}$ (using notation from Lemma 2), the coefficient $\sigma_{s}(x_{0},\widehat{\pi})$ of the characteristic polynomial (11) of

[TABLE]

must vanish. This is equivalent to non-invertibility of $D\mu^{*}(x,\,0)\,P^{*}(x,0))$; see e.g. [7], Remark 4.

∎

Remark 4**.**

Proposition 3 opens a starting point for the computation of nested TFPV: Start with system (10) corresponding to “generic” parameter values in $\Pi_{s}$, i.e. parameter values in the intersection of $\Pi_{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 $\overline{\Pi}_{s}$, and one can use part (b) in order to obtain necessary conditions. Practically this may be realized by determining the decomposition $P\cdot\mu$ for generic $\pi\in\Pi_{s}$ and then looking at zeros of $D\mu\cdot P$, with parameters in the boundary. (The boundary may also contain further parameter values in $\Pi_{s}$.)

4.3 Special settings for chemical reaction networks

For chemical reaction networks (CRN) the parameter region is usually given by $\Pi=\mathbb{R}_{+}^{m}$, thus

[TABLE]

and for many such systems and given $s$, the irreducible components of $W_{s}$ are just determined by the vanishing of certain of the $\pi_{i}$; see e.g. [6, 8, 9]. (The underlying reason for this fact is the subject of forthcoming work.) Thus we have, for $\pi$ in a given irreducible component:

1. (i)

Upon relabelling, there is an $\ell$, $0<\ell<m$ such that $\pi_{i}=0$ for all $i\in\{\ell+1,\cdots m\}$; 2. (ii)

the remaining parameters are nonnegative.

In other words, the intersection of $\Pi_{s}$ with the given irreducible component of $W_{s}$ corresponds to some subset of $\overline{\mathbb{R}_{+}^{\ell}}$, with boundary $\overline{\mathbb{R}_{+}^{\ell}}\setminus\mathbb{R}_{+}^{\ell}$. This leads to an obvious case-by-case analysis. Note that boundary points may or may not be contained in $\Pi_{s}$, but there is no loss in starting with “generic” parameter values in the interior $\mathbb{R}^{\ell}_{>0}$. We look at a particular example.

Example 2**.**

We again consider competitive inhibition; see equation (9). Here the parameters are of the form

[TABLE]

From [8], Proposition 8 we have the necessary condition $e_{0}k_{1}k_{2}k_{-3}=0$ for a TFPV in $\Pi_{1}$, with each of the four cases (e.g. $e_{0}=0$ and all other parameter values $\geq 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 by

[TABLE]

where we have used the abbreviations $e=e_{0}-c_{1}-c_{2}$ and $i=i_{0}-c_{2}$; note that $e\geq 0$ and $i\geq 0$ by design of (9). By Remark 4, nested TFPV’s for dimension two and corresponding points in the critical manifold necesarily satisfy $\det\left(D\mu\cdot P\right)=0$, with

[TABLE]

[TABLE]

Proceeding according to Remark 4, we determine the vanishing set of

[TABLE]

Since all the variables and parameters are nonnegative, this sum equals zero if and only if every summand vanishes. In particular, $k_{-1}\cdot k_{-3}$ has to vanish for any nested TFPV. We look at the two ensuing cases.

1. (i)

$k_{-1}=0$: Then the remaining conditions are

[TABLE]

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. 2. (ii)

$k_{-3}=0$: In this case there remains

[TABLE]

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.

The case $k_{3}=k_{-3}=0$ leads to system (9) with $k_{3}=\varepsilon_{1}k_{3}^{*}$, $k_{-3}=\varepsilon_{1}k_{-3}^{*}$ and $k_{2}=\varepsilon_{1}\varepsilon_{2}k_{2}^{*}$, the reduction of which was discussed in Example 1.**

5 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 $\Pi_{1}\subseteq W_{1}$ which themselves belong to $\Pi_{1}$, as well as certain degenerate cases for $\Pi_{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.

5.1 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 $\Pi_{1}$ there are three remaining cases, viz. $e_{0}=0$, $k_{1}=0$ and $k_{-3}=0$.

1. (a)

For $e_{0}=0$, the system is given by

[TABLE]

We only consider the generic case for $\Pi_{1}$, thus all the remaining parameters are $>0$. Then the (only possible) decomposition $P\cdot\mu$ for the right hand side is given by

[TABLE]

For nested TFPV, a simple computation yields the necessary condition

[TABLE]

with all terms positive; thus every summand must vanish, and in particular

[TABLE]

In case $k_{-1}=k_{2}=0$ system (13) admits a two dimensional variety of stationary points, with $k_{1}s\not=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}=\varepsilon 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}=\varepsilon_{1}\varepsilon_{2}e_{0}^{*}$, hence every term $k_{1}e_{0}s$ is of the form $\varepsilon_{1}^{2}\varepsilon_{2}\cdot(\cdots)$, 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. 2. (b)

Next we consider the system with $k_{1}=0$, i.e.

[TABLE]

Because of Example 2 we may assume that $k_{2}\neq 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}=\varepsilon_{1}\varepsilon_{2}k_{1}^{*}$, $k_{3}=\varepsilon_{1}k_{3}^{*}$ and $k_{-3}=\varepsilon_{1}k_{-3}^{*}$, and (9) becomes

[TABLE]

We compute the reductions for this case. For the auxiliary system (on the critical variety defined by $c_{1}=0$) we obtain the decomposition

[TABLE]

and a straightforward computation yields the projection matrix

[TABLE]

and the auxiliary system

[TABLE]

on the variety defined by $c_{1}=0$. The intermediate reduced system is obtained setting $\varepsilon_{2}=0$.

Turning to the complete reduction, the decomposition of the “fast part” of (14) is given by

[TABLE]

One obtains

[TABLE]

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 $\widetilde{c}_{2}$ of the quadratic equation

[TABLE]

The completely reduced system will automatically yield $\dot{c}_{1}=0$ and $\dot{c}_{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 equation

[TABLE]

with the dot denoting differentiation with respect to $\varepsilon_{1}\varepsilon_{2}t$. 3. (c)

Finally, we deal with the case $k_{-3}=0$, which does not automatically yield a one dimensional variety of stationary points. System (9) becomes

[TABLE]

We may assume that $k_{1}\not=0$ and $k_{2}\not=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}\not=0$ then $i=0$ forces $c_{2}=i_{0}$; the corresponding parameter values are not in $\Pi_{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 system

[TABLE]

To obtain a two dimensional variety of stationary points one has to check the boundary of $\Pi_{1}$ for nested TFPV, which splits into four cases. We only discuss the case $k_{1}=0$ here, thus (9) with small parameters becomes

[TABLE]

The computation of the auxiliary system runs similar to the reduction of (14) and yields

[TABLE]

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; yielding

[TABLE]

5.2 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

[TABLE]

yields, with the usual assumptions and stoichiometry, the differential equation

[TABLE]

where all appearing constants are non-negative. According to Goeke [6], Kap. 9.4, necessary conditions for TFPV are given by

[TABLE]

5.2.1 Case $k_{1}=0$

When we substitute $k_{1}=0$ in equation (16) we obtain

[TABLE]

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

[TABLE]

and necessary conditions for nested TFPV from

[TABLE]

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 $k_{1}=\varepsilon_{1}\varepsilon_{2}k_{1}^{*}$, $k_{-3}=\varepsilon_{1}k_{-3}^{*}$, $k_{4}=\varepsilon_{1}k_{4}^{*}$ . 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

[TABLE]

and with $x=(s,\,c_{1},\,c_{2})^{\rm tr}$ we get

[TABLE]

with

[TABLE]

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)}+\varepsilon_{1}g^{(1,0)}$ is given by

[TABLE]

with $\mu_{1}=c_{1}$, $\mu_{2}=c_{2}$, and

[TABLE]

The projection matrix is

[TABLE]

and the fully reduced system in very slow time on the invariant manifold $c_{1}=c_{2}=0$ is given by the equation

[TABLE]

Similarly one computes the intermediate system on the two dimensional variety given by $c_{1}=0$ from the decomposition $P_{1}\cdot\mu_{1}$:

[TABLE]

5.2.2 Case $e_{0}=0$

From $e_{0}=0$ one also obtains a component of $W_{1}$, and system (16) specializes to

[TABLE]

The right hand side has a factorization $P\cdot\mu$ with

[TABLE]

and in order to obtain nested TFPV we examine all variable-parameter configurations that satisfy

[TABLE]

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 $k_{1}=k_{-1}=k_{2}=0$, $k_{1}=k_{-3}=k_{4}=0$ and $k_{3}=k_{-3}=k_{4}=0$.

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}\not=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}=\varepsilon_{1}k_{1}^{*}$ and $e_{0}=\varepsilon_{1}\varepsilon_{2}e_{0}^{*}$, which implies $k_{1}e_{0}=O(\varepsilon_{1}^{2}\varepsilon_{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

[TABLE]

we compute:

[TABLE]

The intermediate reduced system on the invariant variety $c_{1}=0$ is then:

[TABLE]

5.2.3 Case $k_{-3}=k_{4}=0$

These conditions define a component of $W_{1}$, and generically the critical manifold is given by $s=c_{1}=0$. System (16) is given by

[TABLE]

and the product decomposition (which we do not write down here) yields

[TABLE]

as necessary conditions for nested TFPV. One possible case is $k_{1}=0$ with critical manifold $c_{1}=0$. The remaining cases are:

1. (i)

$k_{2}=k_{-1}=0$ with variety $s=0$; 2. (ii)

$k_{2}=k_{3}=0$ with variety $k_{1}(e_{0}-c_{1}-c_{2})s-k_{-1}c_{1}=0$.

Note that the condition $e_{0}-c_{2}=0$ does not yield a two dimensional critical variety.

5.2.4 Case $k_{2}=0$

In this situation system (16) simplifies to

[TABLE]

The condition $k_{2}=0$ by itself does not define an irreducible component of $W_{1}$; in other words it does not guarantee the existence of a one dimensional variety of stationary points. Therefore we first investigate sufficient conditions, using the observation $\dot{s}+\dot{c}_{1}+2\dot{c}_{2}=-k_{4}c_{2}$.

1. (a)

For $k_{4}\not=0$ this observation implies that any stationary point satisfies $c_{2}=0$, and the remaining condition is $k_{3}sc_{1}=0$.

1. (i)

In case $k_{3}\not=0$ we have either $s=0$, with the variety of stationary points given by $s=c_{2}=0$; in turn this yields the parameter configuration

[TABLE] 2. (ii)

Alternatively we have $c_{1}=0$, the variety is given by $c_{1}=c_{2}=0$, and one must have $k_{1}e_{0}=0$. We obtain the possible parameter configurations

[TABLE] 2. (b)

In case $k_{4}=0$ the remaining system is

[TABLE]

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.

1. (i)

Given that $k_{3}\not=0$ and $k_{-3}\not=0$, the variety is given either by $c_{2}=s=0$, which yields the parameter conditions

[TABLE]

or the variety is given by $c_{1}=c_{2}=0$, with parameter conditions

[TABLE]

all of these are special cases from (a). 2. (ii)

In case $k_{3}=0$ we obtain the one dimensional variety given by $c_{2}=0$ and $k_{1}(e_{0}-c_{1})s-k_{-1}c_{1}=0$; thus we have the parameter condition

[TABLE]

which defines a component of $W_{1}$. 3. (iii)

In case $k_{-3}=0$ one obtains the variety $s=c_{1}=0$, with parameter conditions

[TABLE]

For all these parameters the next task is to discuss conditions for embedded TFPV. We will only do so for two cases.

In the case $k_{-1}=k_{2}=0$ one has a decomposition

[TABLE]

which yields

[TABLE]

We take a closer look at the case $k_{1}=k_{3}=0$, with critical variety $c_{2}=0$. The surface in parameter space

[TABLE]

yields

[TABLE]

The intermediate reduced system is as follows:

[TABLE]

and the completely reduced system (on $s=c_{2}=0$) is given by:

[TABLE] 2. 2.

In the case $k_{2}=k_{3}=k_{4}=0$ we have the one dimensional critical manifold

[TABLE]

and the right hand side of the system at $k_{2}=k_{3}=k_{4}=0$ can be decomposed into $P\cdot\mu$, with

[TABLE]

This yields

[TABLE]

We investigate the case $k_{-3}=0$. Additionally setting $k_{-3}=0$ we obtain the two dimensional critical manifold defined by

[TABLE]

Following the usual procedure we consider the surface

[TABLE]

in parameter space, and thus

[TABLE]

Here the intermediate reduced system is given by

[TABLE]

on $\mu_{2}=0$, and the fully reduced system is given by

[TABLE]

6 More time scales

In this section we give a brief outline on extending the coordinate-free approach to more than three time scales. Thus let $N\geq 3$ and first consider a system with $N-1$ small parameters of the form

[TABLE]

with separated variables. By a smooth coordinate transformation this becomes

[TABLE]

with the very last term in the embedded brackets being $\varepsilon_{N-1}g^{(1,\ldots,1)}$. Here all $g^{(i_{1},\ldots,i_{N-1})}$ are functions of $(x,\varepsilon_{1},\ldots,\varepsilon_{N-1})$. Moreover conditions (i), (ii) preceding Lemma 1 generalize in an obvious manner to the vanishing sets of

[TABLE]

and as in Proposition 2 one obtains decompositions

[TABLE]

Likewise, one generalizes the definitions of $A_{1},\,A_{2}$ and the constructions of the projection matrices $Q_{j}$, the latter extending smoothly to $\varepsilon_{1}=\cdots\varepsilon_{N-1}=0$. This yields the various (intermediate) reductions.

Given a general parameter dependent system (10), nested Tikhonov-Fenichel parameter values may be found via the ansatz

[TABLE]

and the ensuing decomposition of $h(x,\gamma(\varepsilon_{1},\ldots,\varepsilon_{N-1}))$ analogous to the one in (12). Thus the problem is to find $s>0$, $0<k_{1}<\cdots<k_{N-1}$ with $s+k_{N-1}<n$ and a smooth map

[TABLE]

defined in some neighborhood of [math], such that

[TABLE]

whenever all $\varepsilon_{j}>0$. Rather obvious generalizations of Proposition 3 hold, and the strategy outlined in Remark 4 remains applicable.

Appendix

For the reader’s convenience we state and prove here two lemmas.

Lemma 3**.**

Let $r_{1},\,r_{2}$ be positive integers and

[TABLE]

Moreover assume that all eigenvalues of $A$ have real part $<0$. Then the following are equivalent.

1. (i)

All eigenvalues of $-CA^{-1}B+D$ have real part $<0$. 2. (ii)

There exists $\delta>0$ such that, for every $\varepsilon\in(0,\,\delta)$, all eigenvalues of

[TABLE]

have real part $<0$.

Proof.

Consider the singularly perturbed linear differential equation

[TABLE]

Introducing $z:=x+A^{-1}By$ one can rewrite this as

[TABLE]

Here the fast system is just given by $\dot{z}=Az$, and the slow system (on the critical manifold defined by $z=0$) is given by

[TABLE]

Using Tikhonov’s theorem (in the form given e.g. in Verhulst [19], Ch. 8), one sees that both conditions (i), (ii) are equivalent to exponential attractivity of the stationary point [math] for the linear system. ∎

Lemma 4**.**

Let $V\subseteq\mathbb{R}^{n}$ be open and nonempty, $0<r<n$, $\delta>0$ and

[TABLE]

be smooth functions (defined in some neighborhood of $V\times[0,\,\delta)$) such that $\mathbb{R}^{n}$ is the sum of the image $W_{1}$ of $B_{1}$ and the image $W_{2}$ of $B_{2}$, for every $(x,\varepsilon)$. Then the entries of the matrix $Q(x,\varepsilon)\in\mathbb{R}^{n\times n}$ which sends $v\in\mathbb{R}^{n}$ to its $W_{2}$-component with respect to the direct sum decomposition $W_{1}\oplus W_{2}$ depend smoothly on $(x,\varepsilon)$.

Proof.

We suppress the arguments $(x,\varepsilon)$ in the notation. By assumption $C:=\begin{pmatrix}B_{1},&B_{2}\end{pmatrix}$ is invertible, and the entries of $C^{-1}$ depend smoothly on $(x,\varepsilon)$. With the projection matrix given by

[TABLE]

the assertion is obvious.

∎

Acknowledgement. The work of both authors has been supported by the bilateral project ANR-17-CE40-0036 and DFG-391322026 SYMBIONT.