Coordinate-independent singular perturbation reduction for systems with three 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 ^[1] and Fenichel ^[2], 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 ^[3]. More recently, Noel et al. ^[4], Radulescu et al. ^[5], Samal et al. ^[6,7] 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 ^[8] 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 ^[9].

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 ^[10] 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 ^[10]. On the other hand, we extend earlier work ^[11,12] 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 ^[10] 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 algorithm 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 ^[10]. Section 3 generalizes the coordinate-independent reduction algorithm from ^[11] to three-time scale systems. In Section 4 we start from a general parameter dependent system and extend the work from ^[12] 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 ^[10] for a parameter dependent ordinary differential equation system

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

and

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)$^* have nonzero real parts.

^*For a smooth function $g = g(x, y, \ldots)$ we denote the partial derivatives by $D_xg$, $D_yg$ etc.

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

● 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.^†

^†In ^[10] 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 ^[10] we introduce the auxiliary system

on

and the intermediate reduced system

on ${\cal M}_1$. Thus in both equations (2.4) and (2.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

on ${\cal M}_2$, hence the dot in (2.6) denotes differentiation with respect to $\tau_3: = \varepsilon_1 \varepsilon_2 t$. From here on we will follow ^[10] by just using dots for differentiations, with the appropriate time scale evident from the context.

We replace the hyperbolicity conditions from ^[10] by stronger requirements, since in our applications we focus on attracting invariant manifolds.

Definition 1. We say that system (2.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 ^[10], Theorems A, B and Corollary A. Some of our statements are informal; for rigorous statements and pertinent definitions we refer to ^[10].

Theorem 1. Let system (2.1) be given, with (HA) satisfied.

(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 (2.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 (2.1) on $\mathcal N_{\varepsilon_1, \varepsilon_2}$ converge to solutions of (2.6) on $\mathcal N$.

(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 (2.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 (2.1) on $\mathcal L_{\varepsilon_1, \varepsilon_2}$ converge to solutions of (2.5) on $\mathcal L$.

Note that the passage from (2.1) to (2.6) may be seen as a singular perturbation reduction for two time scales. From this perspective, system (2.5) indeed represents an intermediate step. For systems with two time scales this leads to the question whether intermediates exist.

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.

(i) $D_{x_1}f_1(x, 0, 0)$ has only eigenvalues with negative real parts on $\mathcal M_1$, and

has only eigenvalues with negative real parts on $\mathcal M_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

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

one gets by the chain rule

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 (ⅰ). The equivalence of (ⅰ) and (ⅱ) follows from Lemma 3 in the Appendix.

Remark 1. (a) One may rewrite systems (2.1) through (2.6) to some extent, with no effect on the reductions. Using Hadamard's lemma, one may restate (2.1) as

with only the $\widehat 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 (2.1) to the completely reduced system (2.6) can evidently be obtained in the following manner: Fix $\varepsilon_1>0$ and reduce (2.1) with respect to the small parameter $\varepsilon_2$ (in time scale $\varepsilon_2t$). 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 ^[11,13] to the three-time scale setting. The first task is to intrinsically characterize those systems which admit a transformation to "standard form" (2.1). Reversing matters, applying a (local) smooth coordinate transformation to Eq (2.1) yields a smooth system

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:

(ⅰ) 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$.

(ⅱ) For all sufficiently small $\varepsilon_1>0, \, \varepsilon_2\geq 0$ the zeros of

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_xg^{(0, 0)}(y, \varepsilon_1, \varepsilon_2)+\varepsilon_1D_xg^{(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 (3.1) is in the special form

adjusting conditions (ⅰ) and (ⅱ) accordingly. Conditions (ⅰ) and (ⅱ) are certainly necessary for (3.1) or (3.2) to be a transformed version of (2.1). The first part of the next lemma shows sufficiency.

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

(b) Condition (i) for system (3.2) 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

and moreover there is a $\theta_1>0$ such that

has only eigenvalues with real part $\leq -\theta_1$, for all $x\in U_{1, y}$.

(c) In presence of condition (i), condition (ii) for system (3.2) 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

and moreover there is a $\theta_2>0$ such that

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 ^[2], Lemma 5.3 for smooth vector fields, and shown in ^[13], 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_1g^{(1, 0)}$ in a neighborhood of $y$. These first integrals determine slow and "very slow" variables. Parts (b) and (c) are straightforward applications of ^[11], 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 ^[11] for details.

Next we use the decompositions to compute reductions.

Proposition 2. (a) In arbitrary coordinates the reduction corresponding to the passage from system (2.1) to the auxiliary system may be obtained as follows:

Given $\varepsilon_2\geq 0$, determine the projection matrix

The auxiliary system (2.4) for $\varepsilon_2>0$ then corresponds to

on the local invariant manifold defined by $\mu_1(x, \varepsilon_2) = 0$. The equation corresponding to the intermediate reduced system (2.5) is obtained by setting $\varepsilon_2 = 0$.

(b) In arbitrary coordinates the reduction corresponding to the passage from system (2.1) to the completely reduced system (2.6) may be obtained as follows:

Given $\varepsilon_1>0$, determine the projection matrix

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

on the local invariant manifold defined by $\mu_1(x, 0) = \mu_2(x, 0) = 0$.

Proof. Part (a) is a direct application of ^[11], 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

along the image, for $x\in \widetilde {\mathcal M}_2$ (see ^[11], 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

with invertible $M_1$, whence the kernel is equal to the column space $W_2$ of the matrix

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 3.2 with the help of the first integrals noted in the proof of Lemma 1(a); see ^[11], 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 ^[11] 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 ^[14]. 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.) There is a familiar two time scale reduction for this system, viz. the quasi-steady state (QSS) reduction for both complexes; see Eqs (1.51–1.52) in Keener and Sneyd ^[14]. As was shown in ^[15], section 3.2 this QSS reduction is (up to irrelevant higher order terms) identical to the singular perturbation reduction with small parameter $e_0 = \varepsilon e_0^*$, hence Tikhonov guarantees its validity. (It is not generally true that QSS reductions amount to Tikhonov-Fenichel reductions. A sufficient condition is stated in ^[16], Prop. 5, and a list of possible QSS reductions for competitive inhibition - including a characterization which correspond to singular perturbation reductions - is given in section 5.2 of that paper.)

Example 1. In system (3.3) set $x = (s, \, c_1, \, c_2)^{\rm tr}$ and assume $k_2 = \varepsilon_1\varepsilon_2k_2^*$, $k_3 = \varepsilon_1k_3^*$ and $k_{-3} = \varepsilon_1k_{-3}^*$. (Colloquially 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. The related two time scale reduction to a one dimensional equation corresponds to slow formation of enzyme and complex.) This is of the type (3.2), with

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

$\widetilde{\mathcal{M}}_1$ is contained in the zero set of $\mu_1$, and we have

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

one has

and furthermore

Application to

yields the auxiliary system (in time scale $\varepsilon_1 t$) on $\widetilde{\mathcal{M}}_1$:

Setting $\varepsilon_2 = 0$ one obtains the intermediate reduced system.

When the initial values for system (3.3) 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

which leads to quadratic equations for $s$ and $c_1$.

To find the fully reduced system one first computes

and

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

with

and

restricted to the invariant curve $\widetilde{\mathcal M}_2$.

Finally, given initial values $(s_0, c_{1, 0}, c_{2, 0})$ for system (3.3), 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 (2.1) or (3.1) 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 ^[12] 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 ^[12] for polynomial (or rational) parameter dependent systems

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

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 ^[12] Lemmas 1 and 2, for which the characteristic polynomial

of the Jacobian $D_xh(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 (4.2)) satisfies

(i) $\sigma_0(x_0, \widehat\pi) = \cdots = \sigma_{s-1}(x_0, \widehat\pi) = 0$;

(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 ^[17], 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. While there are notable exceptions for chemical reaction networks, with possible first integrals from stoichiometry, this condition is the most problematic to verify in general. One may proceed as follows: For every $d>0$ there is an induced action of $D_xh(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 (ⅰ), the characteristic polynomial of this action (which coincides with (4.2) 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, and thanks to Hilbert's Basissatz, finitely many of these equations suffice. With increasing dimension, one has to deal with feasibility problems. In applications the common procedure is step by step, looking for simplifications at every stage. A full account is given in ^[12].

For the remainder of this section we assume that $\Pi\subseteq \mathbb R_+^m$ is a semi-algebraic set, and that system (4.1) admits the positively invariant subset $\mathbb R^n_+$. Then, as was shown in ^[12], 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 ^[12]), and taking into account the special form of (3.1), it seems reasonable to consider surfaces in parameter space. Thus consider a smooth surface of the special form

defined in some nighborhood of $(0, \, 0)$. Substitute $\gamma(\varepsilon_1, \, \varepsilon_2)$ for $\pi$ in (4.1) to get

with the $g^{(i, j)}$ smooth by Hadamard's lemma. In order to obtain a system (3.1) that also satisfies the conditions (ⅰ) and (ⅱ) 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 (4.1) and $s, \, k>0$ with $s+k < n$, let $\delta>0$ and let

be a smooth curve such that

(i) $\beta(\varepsilon_1)\in \Pi_s$ for all $\varepsilon_1>0$,

(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. (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$.

(b) Let $\beta$ as in Definition 2, and for $\varepsilon_1>0$ consider the decomposition

according to ^[11], Theorem 1. Then

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 (4.2) of

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

Remark 4. Proposition 3 opens a starting point for the computation of nested TFPV: Start with system (4.1) 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

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. ^[12,16,18]. (The underlying reason for this fact is the subject of forthcoming work.) Thus we have, for $\pi$ in a given irreducible component:

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

(ⅱ) 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 Eq (3.3). Here the parameters are of the form

From ^[12], Proposition 8 (see also Goeke's dissertation ^[18], section 9.3) we have the necessary condition $e_0k_1k_{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. For the cases with some $k_i = 0$, which amount to a distinction of slow and fast reactions, the two time scale reductions may be obtained using a computational procedure outlined in Heinrich and Schauer ^[19]; see also Stiefenhofer ^[22]. 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

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 (3.3). 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

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}\cdot k_{-3}$ has to vanish for any nested TFPV. We look at the two ensuing cases.

(ⅰ) $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.

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

In view of case (ⅰ) 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 (3.3) with $k_3 = \varepsilon_1k_3^*$, $k_{-3} = \varepsilon_1k_{-3}^*$ and $k_2 = \varepsilon_1\varepsilon_2k_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. The mathematical analysis yields a rather large number of possible reductions, with parameter conditions that allow an interpretation in a biochemical context. It seems that most of these reductions have not been discussed before in the literature. Here one should note that the principal interest in biochemistry (as in the monograph by Keener and Sneyd) lies in identifying and discussing reaction velocities for given parameter constellations; thus the focus is different from the one in the present paper. (Noel et al. ^[4], Radulescu et al. ^[5], Samal et al. ^[6,7] investigate questions related to those in the present paper.) 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 Eq (3.3), 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$.

(a) For $e_0 = 0$, the system is given by

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

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 (5.1) admits a two dimensional variety of stationary points, with $k_1s\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_1s = 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_2e_0^*$, hence every term $k_1e_0s$ 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 (5.1) leads to $k_1s = 0$.

To summarize, for the case $e_0 = 0$, which leads to the familiar QSS reduction in the two time scale setting, there are no interesting reductions for the three time scale scenario: No parameter value leads to an intermediate reduction to dimension two with nontrivial complete reduction.

(b) Next we consider the system with $k_1 = 0$, i.e.

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_2k_1^*$, $k_3 = \varepsilon_1 k_3^*$ and $k_{-3} = \varepsilon_1k_{-3}^*$, and (3.3) 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

and the auxiliary system

on the variety defined by $c_1 = 0$. (This is a QSS reduction for complex $C_1$, according to ^[16], Proposition 5.) The intermediate reduced system is obtained setting $\varepsilon_2 = 0$.

Turning to the complete reduction, the decomposition of the "fast part" of (5.2) 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 $\widetilde c_2$ of the quadratic equation

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

with the dot denoting differentiation with respect to $\varepsilon_1\varepsilon_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 (3.3) becomes

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

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 (3.3) with small parameters becomes

The computation of the auxiliary system runs similar to the reduction of (5.2) 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 ^[16], Proposition 5 the reduced system may be directly obtained via "classical" QSS reduction; yielding

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

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

where all appearing constants are non-negative; see Keener and Sneyd ^[14], equations (1.62–64), where the QSS reduction for complexes $C_1$ and $C_2$ is discussed. A direct proof of validity via singular perturbation theory - for two time scales - is given in 21ray ^[21], Section 6.5 when $e_0 = \varepsilon e_0^*$; alternatively ^[16], Proposition 5 is applicable. In ^[16], subsection 5.4 a list of parameter conditions for various types of QSS is given, and their relation to Tikhonov-Fenichel reductions is discussed. Keener and Sneyd also mention the (partial) equilibrium approximation, by which one generally determines the product forming velocity given equilibrium in a certain subnetwork. In the present context ^[14], section 1.7, Exercise 10 this means that $k_2 = \varepsilon k_2^*$ and $k_4 = \varepsilon k_4^*$ are small. The discussion in 5.2.4 shows, among other things, that these conditions are not compatible with singular perturbation reductions. (One needs to add further "small parameters", and the resulting reduction formulas are more complicated.) Cooperative systems with three complexes are mentioned in ^[14], and the singular perturbation reduction in the case of arbitrary many complexes, with small $e_0 = \varepsilon e_0^*$, is derived in ^[16], section 5.3. Now we turn to the three time scale setting.

According to Goeke ^[18], Kap. 9.4, necessary conditions for TFPV are given by

5.2.1.   Case $k_1=0$

When we substitute $k_1 = 0$ in Eq (5.4) we obtain

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 $k_1 = \varepsilon_1\varepsilon_2 k_1^*$, $k_{-3} = \varepsilon_1k_{-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

and with $x = (s, \, c_1, \, c_2)^{\rm tr}$ we get

with

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

with $\mu_1 = c_1$, $\mu_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 equation

Similarly one computes the intermediate system on the two dimensional variety given by $c_1 = 0$ from the decomposition $P_1\cdot\mu_1$:

In this context, a brief digression may be appropriate: In their discussion of reaction velocities, Keener and Sneyd ^[14] also derive a Hill function in Eq (1.69)) with arguments that could a priori be interpreted as involving three time scales: They consider "$k_1\to 0$ and $k_3\to\infty$, with constant product $k_1k_3$". Rescaling time, one might think of this as a case of very small $k_1$, with $k_3$ of order one. One system satisfying these conditions is (5.5), but the completely reduced equation manifestly does not follow Hill kinetics. One might also consider the case when $k_3$ is of order one, $k_1$ is very small and all the other $k_i$ are small. This corresponds to a degenerate (but a priori legitimate) subcase of "$k_{-3} = k_4 = 0$" discussed above, with the surface in parameter space now given by

A straightforward analysis shows, however, that the resulting system does not satisfy the rank (or eigenvalue) conditions from Lemma 1(c), hence Tikhonov-Fenichel reduction is not applicable. To summarize: For the given cooperative network, a reaction velocity of Hill type cannot be obtained from time scale arguments.

5.2.2.   Case $e_0=0$

From $e_0 = 0$ one also obtains a component of $W_1$, and system (5.4) specializes to

The right hand side has a factorization $P\cdot \mu$ 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 $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_3sc_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_1k_1^*$ and $e_0 = \varepsilon_1\varepsilon_2 e_0^*$, which implies $k_1e_0 = O(\varepsilon_1^2\varepsilon_2)$. We consider one special 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 $c_1 = 0$ is then:

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 (5.4) is given by

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

as necessary conditions for nested TFPV. One possible case is $k_1 = 0$ with critical manifold $c_1 = 0$. The remaining cases are:

(ⅰ) $k_2 = k_{-1} = 0$ with variety $s = 0$;

(ⅱ) $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 (5.4) simplifies to

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_4c_2$.

(a) For $k_4\not = 0$ this observation implies that any stationary point satisfies $c_2 = 0$, and the remaining condition is $k_3sc_1 = 0$.

(ⅰ) 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

(ⅱ) Alternatively we have $c_1 = 0$, the variety is given by $c_1 = c_2 = 0$, and one must have $k_1e_0 = 0$. We obtain the possible parameter configurations

(b) In case $k_4 = 0$ the remaining system is

Adding the first two equations for stationary points shows that $k_{-3}c_2 = 0$, and combining this with the third equation yields $k_3sc_1 = 0$; in addition one has $k_1es-k_{-1}c_1 = 0$. Thus there are further conditions for the existence of a one dimensional critical variety.

(ⅰ) 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

or the variety is given by $c_1 = c_2 = 0$, with parameter conditions

all of these are special cases from (a).

(ⅱ) 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

which defines a component of $W_1$.

(ⅲ) In case $k_{-3} = 0$ one obtains the variety $s = c_1 = 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

which yields

We take a closer look at the case $k_1 = k_3 = 0$, with critical variety $c_2 = 0$. The surface in parameter space

yields

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

and the right hand side of the system at $k_2 = k_3 = k_4 = 0$ can be decomposed into $P\cdot \mu$, with

This yields

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

Following the usual procedure we consider the surface

in parameter space, and thus

Here the intermediate reduced system is given by

on $\mu_2 = 0$, and the fully reduced system is given by

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

with separated variables. By a smooth coordinate transformation this becomes

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 (ⅰ), (ⅱ) preceding Lemma 1 generalize in an obvious manner to the vanishing sets of

and as in Proposition 2 one obtains decompositions

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 (4.1), nested Tikhonov-Fenichel parameter values may be found via the ansatz

and the ensuing decomposition of $h(x, \gamma(\varepsilon_1, \ldots, \varepsilon_{N-1}))$ analogous to the one in (4.3). 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

defined in some neighborhood of $0$, such that

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

Acknowledgments

The work of both authors has been supported by the bilateral project ANR-17-CE40-0036 and DFG-391322026 SYMBIONT. The authors thank two anonymous reviewers for helpful comments and suggestions.

Conflict of interest

Both authors declare no conflicts of interest in this paper.

Appendix

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

Lemma 3. Let $r_1, \, r_2$ be positive integers and

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

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

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

have real part $< 0$.

Proof. Consider the singularly perturbed linear differential equation

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

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

Using Tikhonov's theorem (in the form given e.g. in Verhulst ^[23], Ch. 8), one sees that both conditions (ⅰ), (ⅱ) are equivalent to exponential attractivity of the stationary point $0$ for the linear system.

Lemma 4. Let $V\subseteq \mathbb R^n$ be open and nonempty, $0 < r < n$, $\delta>0$ and

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

the assertion is obvious.