and thus, with $\dot{c}_{1}=\dot{c}_{2}=0$,

$\dot{s}=-\cfrac{k_{1}e_{0}k_{2}s}{k_{1}s+k_{2}+k_{-1}},\quad s(0)=s_{0},$
(72)

which is valid asymptotically as $\varepsilon\to 0$. Here, we recover the
familiar Michaelis–Menten equation in the limit when the concentrations of
both enzyme and inhibitor approach zero of order $\varepsilon$. (The same
reduction is obtained for $k_{3}=\varepsilon k_{3}^{*}$ and $e_{0}=\varepsilon
e_{0}^{*}$.)

From a different perspective, when the term $k_{3}i_{0}$ vanishes, the
subspace $W:=\\{(s,c_{1},c_{2})\in\mathbb{R}^{3}:c_{2}=0\\}$ is invariant, and
a slight perturbation (not necessarily of order $\varepsilon$) results in the
near–invariance of $W$. Considering the expressions (69) and (70) for
$\varepsilon^{*}$ and $\mu^{*}$, respectively, one sees that $k_{3}i_{0}\to 0$
has no strong effect on these parameters, and that $\varepsilon_{MM}$ is a
good upper estimate for $\varepsilon_{U}$. One may rewrite (68) as

$\dot{s}=-\cfrac{k_{1}e_{0}k_{2}s}{k_{1}s(1+k_{3}i_{0}/k_{-3})+k_{-1}+k_{2}},$
(73)

thus, when $k_{3}i_{0}/k_{-3}\ll 1$, then the standard Michaelis–Menten
reduction is approximately valid. In this case, the dynamics are effectively
two-dimensional. Hence (for the given initial values) the magnitude of
$\mu_{U}$ is irrelevant, and (68) will hold even if $1<k_{1}e_{0}/k_{-3}$ and
$1<\mu_{U}$ 202020The implicit assumptions in the proof of Proposition 3 are
then not satisfied. since $\varepsilon_{MM}\ll 1$ automatically ensures the
validity of (68) when $W$ is nearly invariant (see, Figure 8).

Figure 8: Uncompetitive inhibition reaction mechanism: Near-invariance may
lead to scenarios in which the reduction 68) is accurate even when
$1<\mu_{U}$. The solid black curve is the numerical solution to (65). The
broken red curve is the numerical solution to (68). Time has been mapped to
the $\tau$ scale: $\tau=t/T$, $\tau\in[0,1]$. The parameters (in arbitrary
units) are: $k_{1}=1.0$, $e_{0}=1.0$, $k_{-3}=10^{-1}$, $k_{3}=10^{-2}$,
$i_{0}=10^{-3}$, $k_{2}=k_{-1}=10^{2}$ with $\varepsilon_{U}\approx 2.5\times
10^{-3}$ and $\mu_{U}\approx 5.0.$ Left panel: $s_{0}=10.0$. Right panel: The
long-time validity is verified with $s_{0}=10^{-1}$. Note that
$k_{1}e_{0}/k_{-3}$ is large, thus $\mu_{U}$ is of order one. However, since
$k_{3}i_{0}/k_{-3}=10^{-2}$ and $c_{1}(0)=c_{2}(0)=0$, the dynamics prior to
reduction are essentially two-dimensional. Consequently, (68) holds since
$\varepsilon_{MM}\ll 1$.

###  6.3 Competitive inhibition reaction mechanism

The irreversible competitive inhibition reaction mechanism

$\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}$
(74)

corresponds (with mass action kinetics and stoichiometric conservation laws)
to the ODE system

$\begin{array}[]{rclclclclcl}\dot{s}=&-&k_{1}(e_{0}-c_{1}-c_{2})s&+&k_{-1}c_{1},&\\\
\dot{c}_{1}=&&k_{1}(e_{0}-c_{1}-c_{2})s&-&(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}$
(75)

The usual initial conditions are $s(0)=s_{0},\,e(0)=e_{0},\,i(0)=i_{0}$ and
$c_{1}(0)=c_{2}(0)=p(0)=0$. We fix a reference value $e_{0}^{*}$, then from
the conservation laws we obtain the compact positively invariant set

$K:=\\{(s,c_{1},c_{2})\in\mathbb{R}_{\geq 0}^{3}:0\leq s\leq s_{0},\;\;0\leq
c_{1}\leq e_{0}^{*},\;\;0\leq c_{2}\leq\min\\{e_{0}^{*},\,i_{0}\\}\\}.$

The parameter space $\Pi=\mathbb{R}_{\geq 0}^{8}$ has elements

$(e_{0},s_{0},i_{0},k_{1},k_{-1},k_{2},k_{3},k_{-3})^{\rm tr},$

and it is known that $e_{0}=0$, with all other parameters positive, defines a
TFPV,

$\widehat{\pi}:=(0,s_{0},i_{0},k_{1},k_{-1},k_{2},k_{3},k_{-3})^{\rm tr}$ (76)

with corresponding critical manifold

$\widetilde{Y}:=\\{(s,c_{1},c_{2})\in\mathbb{R}^{3}_{\geq
0}:c_{1}=c_{2}=0\\}.$

(see, below for nondegeneracy conditions on the remaining parameters.) We
choose the parameter ray direction $\rho=(e_{0}^{*},0,\ldots,0)^{\rm tr}$,
with $e_{0}=\varepsilon e_{0}^{*}$. The singular perturbation reduction (see,
Goeke and Walcher [20, Section 3.2]) yields the equation

$\dot{s}=-\dfrac{k_{1}k_{-3}k_{2}e_{0}s}{(k_{1}s+k_{-1}+k_{2})k_{-3}+k_{3}i_{0}(k_{-1}+k_{2})}.$
(77)

The reduced equation (68) has been previously reported in the literature (see,
e.g. Schnell and Mendoza [37]. Note that the reduction (77) again differs from
the classical QSS reduction (see e.g. Keener and Sneyd [24, Section 1.4.3]).
However, (77) and the classical reduction agree up to a term of order
$\varepsilon^{2}$ and are therefore asymptotically equivalent.

####  6.3.1 Asymptotic small parameters

The coefficients of the characteristic polynomial on the critical manifold are

$\begin{array}[]{rcl}\sigma_{1}&=&k_{1}s+k_{-1}+k_{2}+k_{3}i_{0}+k_{-3};\\\
\sigma_{2}&=&k_{-3}k_{1}s+(k_{3}i_{0}+k_{-3})(k_{-1}+k_{2});\\\
\widehat{\sigma}_{3}&=&k_{2}k_{1}e_{0}^{*}\cdot\left(k_{3}i_{0}+k_{-3}\right).\end{array}$

We generally assume that all parameters are contained in a compact subset of
the positive orthant, hence are bounded above by certain positive constants.
Moreover $\sigma_{1}(\widehat{\pi},s)$ and $\sigma_{2}(\widehat{\pi},s)$
satisfy the TFPV property

$\begin{array}[]{rcl}k_{3}i_{0}+k_{-3}+k_{-1}+k_{2}&=&\min\sigma_{1}>0,\\\
(k_{-1}+k_{2})(k_{3}i_{0}+k_{-3})&=&\min\sigma_{2}>0\end{array}$

if and only if $k_{-1}+k_{2}$ and $k_{3}i_{0}+k_{-3}$ are bounded below by
certain positive constants. More restrictively, we will assume that $i_{0}$ is
bounded below by some positive constant. Finally $k_{2}$ and $k_{-3}$ should
be bounded below by positive constants, lest the reduced equation (77) is too
close to trivial.

With $\sigma_{1},\sigma_{2}$ and $\widehat{\sigma}_{3}$ evaluated at
$\widehat{\pi}$, we obtain the distinguished small parameter

$\begin{array}[]{rcl}\varepsilon^{*}&=&\varepsilon\displaystyle\sup_{\widehat{Y}\cap
K}\cfrac{\widehat{\sigma}_{3}(s,\widehat{\pi},\rho,0)}{\sigma_{1}(s,\widehat{\pi})\sigma_{2}(s,\widehat{\pi})}\\\
&=&\cfrac{k_{2}k_{1}e_{0}}{(k_{-1}+k_{2})^{2}}\cdot\cfrac{k_{-1}+k_{2}}{k_{-1}+k_{2}+k_{3}i_{0}+k_{-3}}\\\
&=&\varepsilon_{MM}\cdot\cfrac{k_{-1}+k_{2}}{k_{-1}+k_{2}+k_{3}i_{0}+k_{-3}}:=\varepsilon_{I}\\\
\end{array}$ (78)

To verify the equalities, note that $\widehat{\sigma}_{3}$ is constant while
$\sigma_{1},\,\sigma_{2}$ are increasing with $s$.

It is straightforward to check that $\sigma_{1}^{2}-4\sigma_{2}\geq 0$, thus
all eigenvalues are real and Proposition 3 is applicable. Determining the
parameter $\mu^{*}$ requires a distinction of cases. The derivative of

$s\mapsto
q(s):=\cfrac{\widehat{\sigma}_{3}(s,\widehat{\pi},\rho,0)\sigma_{1}(s,\widehat{\pi})}{\sigma_{2}(s,\widehat{\pi})^{2}}$

is a rational function in $s$ with numerator of degree one. Both coefficients
are negative if and only if

$2k_{-3}(k_{3}i_{0}+k_{-3})+k_{-3}(k_{-1}+k_{2})\geq(k_{-1}+k_{2})i_{0}k_{3},$
(79)

otherwise they have opposite signs. Note that (79) is satisfied whenever
$(k_{3}i_{0})/(k_{-3})\leq 1$. This inequality admits a direct interpretation
in terms of the reaction mechanism. On the one hand, it places a lower bound
on the allowable size of $k_{-3}$. More importantly, it holds whenever the
inhibitor concentration is not too high, thus it is controllable by
experimental design.

When (79) holds then $s\mapsto q(s)$ is strictly decreasing for $s\geq 0$, and

$\mu^{*}=\mu_{I}^{(1)}:=\varepsilon\cfrac{\widehat{\sigma}_{3}(0,\widehat{\pi},\rho,0)\sigma_{1}(0,\widehat{\pi})}{\sigma_{2}(0,\widehat{\pi})^{2}}=\varepsilon_{MM}\cdot\cfrac{k_{-1}+k_{2}+k_{3}i_{0}+k_{-3}}{k_{3}i_{0}+k_{-3}}.$
(80)

Whenever (79) does not hold then a straightforward calculation shows that the
maximum of $s\mapsto q(s)$ for $0\leq s<\infty$ is given by

$\begin{array}[]{rcl}\mu^{*}=\mu_{I}^{(2)}&=&\dfrac{k_{2}k_{1}e_{0}\cdot(k_{3}i_{0}+k_{-3})}{4k_{-3}\cdot\bigg{(}k_{3}i_{0}(k_{-1}+k_{2})-k_{-3}(k_{3}i_{0}+k_{-3})\bigg{)}}\\\
&=&\varepsilon_{MM}\cdot\dfrac{(k_{-1}+k_{2})^{2}\cdot(k_{3}i_{0}+k_{-3})}{4k_{-3}\cdot\bigg{(}k_{3}i_{0}(k_{-1}+k_{2})-k_{-3}(k_{3}i_{0}+k_{-3})\bigg{)}}.\\\
\end{array}$ (81)

This expression is somewhat unwieldy. But $\mu_{I}^{(2)}$ admits an obvious
lower bound, obtained by discarding the negative term in the denominator:

$\mu_{I}^{(2)}\geq\cfrac{1}{4}\cdot\cfrac{k_{1}k_{2}e_{0}}{k_{-3}(k_{-1}+k_{2})}\cdot\bigg{(}1+\cfrac{k_{-3}}{k_{3}i_{0}}\bigg{)}.$

Moreover, the negation of (79) provides an estimate for the denominator which
yields an upper bound

$\mu_{I}^{(2)}\leq\cfrac{1}{4}\cdot\cfrac{k_{2}k_{1}e_{0}(k_{3}i_{0}+k_{-3})}{k_{-3}^{2}(k_{3}i_{0}+k_{-3}+k_{-1}+k_{2})}.$

The lower bound shows that is it necessary to require
$k_{1}e_{0}\ll\min\\{k_{-3},k_{3}i_{0}\\}$ whenever $k_{-1}$ and $k_{2}$ are
of the same order.

Finally, whenever $s_{0}$ is not too large one may also consider the estimate

$\begin{array}[]{rcl}\mu^{*}&\leq&\varepsilon\cfrac{\widehat{\sigma}_{3}(s_{0},\widehat{\pi},\rho,0)\sigma_{1}(s_{0},\widehat{\pi})}{\sigma_{2}(0,\widehat{\pi})^{2}}\\\
&=&{\cfrac{k_{2}k_{1}e_{0}}{(k_{-1}+k_{2})^{2}}}\cdot\cfrac{k_{1}s_{0}+k_{-1}+k_{2}+k_{3}i_{0}+k_{-3}}{k_{3}i_{0}+k_{-3}}\\\
&=&\varepsilon_{MM}\cdot\cfrac{k_{1}s_{0}+k_{-1}+k_{2}+k_{3}i_{0}+k_{-3}}{k_{3}i_{0}+k_{-3}}=:\widetilde{\mu}_{I},\end{array}$
(82)

which is a direct consequence of monotonicity properties of the $\sigma_{j}$.
This inequality is exact whenever (79) does not hold and $s_{0}$ is smaller
than the argument of $\max q$.

####  6.3.2 Numerical simulations

Generally, by Fenichel theory the accuracy of the reduction (77) improves
along the parameter ray as $\varepsilon_{I}\to 0$ and $\mu_{I}^{(i)}\to 0$,
for $i=1,2$, respectively. Continuing the procedure employed in the previous
case studies, we illustrate the efficacy of the qualifiers $\varepsilon_{I}\ll
1,\mu_{I}^{(i)}\ll 1$ (with appropriate index $i$) with several numerical
simulations:

  1. 1.

For our first example, we once again consider the case in which all parameters
except $e_{0}$ are equal, which is a representative of parameters of the same
magnitude. Numerical simulations confirm that the accuracy of (77) improves as
$\varepsilon_{I}\to 0$ and $\mu_{I}^{(1)}\to 0$ (see, Figure 9).

Figure 9: Competitive inhibition reaction mechanism: With parameters of unit
magnitude, and (79) valid, numerical simulations indicate that the accuracy of
(77) improves along the parameter ray as both $\varepsilon_{I}\to 0$ and
$\mu_{I}^{(1)}\to 0$. In all panels, the parameters (in arbitrary units) are:
$s_{0}=1.0$, $k_{1}=1.0$, $k_{2}=1.0$, $k_{-1}=1.0$, $k_{3}=1.0$, $k_{-3}=1.0$
and $i_{0}=1.0$. The solid black curve is the numerical solution to the mass
action system (75). The broken red curve is the numerical solution to (77).
Time has been mapped to the $\tau$ scale: $\tau=t/T$, $\tau\in[0,1]$ Top Left
panel: $e_{0}=1.0$ with $\varepsilon_{I}=1.25\times 10^{-1}$ and
$\mu_{I}^{(1)}=5\times 10^{-1}$. Top Right panel: $e_{0}=10^{-1}$ with
$\varepsilon_{I}=1.25\times 10^{-2}$ and $\mu_{I}^{(1)}=5\times 10^{-2}$.
Bottom Left panel: $e_{0}=10^{-2}$ with $\varepsilon_{I}=1.25\times 10^{-3}$
and $\mu_{I}^{(1)}=5\times 10^{-3}$. The reduction 77) is nearly
indistinguishable from (75). Bottom Right panel: $e_{0}=10^{-4}$ with
$\varepsilon_{I}=1.25\times 10^{-4}$ and $\mu_{I}^{(1)}=5\times 10^{-4}$. The
QSS reduction (77) is again practically indistinguishable from (75).

  2. 2.

In our second example, we demonstrate the effectiveness of $\varepsilon_{I}$
and $\mu_{I}^{(i)}$ with parameters that have disparate magnitudes. We observe
that $\mu_{I}^{(1)}$ is the definitive indicator of the accuracy of (77) when
(79) holds, while $\mu_{I}^{(2)}$ is the indicator of the accuracy of (77)
whenever (79) fails, reflecting the fact that one eigenvalue must have much
smaller absolute value than the other two (see, Figures 10 and 11).

Figure 10: Competitive inhibition reaction mechanism: When parameter values
are disparate in magnitude, numerical simulations indicate that the accuracy
of (77) improves along the parameter ray as $\mu_{I}^{(1)}\to 0$ when (79)
holds. In all panels, the parameters (in arbitrary units) are: $s_{0}=10.0$,
$k_{1}=1.0$, $k_{2}=k_{-1}=10^{2}$, $k_{3}=k_{-3}=10^{-1}$ and $i_{0}=1.0$.
The solid black curve is the numerical solution for $s$ to the mass action
system (75). The broken red curve is the numerical solution to (77). Time has
been mapped to the $\tau$ scale: $\tau=t/T$, $\tau\in[0,1]$. Top Left panel:
$e_{0}=1.0$, $\varepsilon_{I}\approx 2.5\times 10^{-3}$, and
$\mu_{I}^{(1)}\approx 2.5$. Top Right panel: $e_{0}=10^{-1}$,
$\varepsilon_{I}\approx 2.5\times 10^{-4}$, and $\mu_{I}^{(1)}\approx
2.5\times 10^{-1}$. Bottom Left panel: $e_{0}=10^{-2}$,
$\varepsilon_{I}\approx 2.5\times 10^{-5}$, and $\mu_{I}^{(1)}\approx
2.5\times 10^{-2}$. Bottom Right panel: $e_{0}=10^{-3}$,
$\varepsilon_{I}\approx 2.5\times 10^{-6}$, $\mu_{I}^{(1)}\approx 2.5\times
10^{-3}$ and the QSS reduction (77) is nearly indistinguishable from (75).
Collectively, these simulations indicate that $\mu_{I}^{(1)}\ll 1$ is the
qualifier that ensures the validity of (77) when (79) holds.

Figure 11: Competitive inhibition reaction mechanism: When parameter values
are disparate in magnitude, numerical simulations indicate that the accuracy
of (77) improves along the search direction as $\mu_{I}^{(2)}\to 0$ when (79)
fails. In all panels, the parameters (in arbitrary units) are: $s_{0}=10.0$,
$k_{1}=1.0$, $k_{2}=k_{-1}=10^{2}$, $k_{-3}=10^{-2}$ and $k_{3}=i_{0}=1.0$.
The solid black curve is the numerical solution for $s$ to the mass action
system (75). The broken red curve is the numerical solution to (77). Time has
been mapped to the $\tau$ scale: $\tau=t/T$, $\tau\in[0,1]$. Top Left panel:
$e_{0}=1.0$, $\varepsilon_{I}\approx 2.5\times 10^{-3}$, and
$\mu_{I}^{(2)}\approx 12.6$. Top Right panel: $e_{0}=10^{-1}$,
$\varepsilon_{I}\approx 2.5\times 10^{-4}$, and $\mu_{I}^{(2)}\approx
12.6\times 10^{-1}$. Bottom Left panel: $e_{0}=10^{-2}$,
$\varepsilon_{I}\approx 2.5\times 10^{-5}$, and $\mu_{I}^{(2)}\approx
12.6\times 10^{-2}$. Bottom Right panel: $e_{0}=10^{-3}$,
$\varepsilon_{I}\approx 2.5\times 10^{-6}$, and $\mu_{I}^{(2)}\approx
12.6\times 10^{-3}$. Observe the QSS reduction (77) is nearly
indistinguishable from (75), which indicates that $\mu_{I}^{(2)}\ll 1$ is the
qualifier that ensures the validity of (77) when (79) fails.

####  6.3.3 Near–invariance

As in the previous sections, we now discuss special instances of near-
invariance. The inhibitory mechanism can be turned off by requiring
$k_{3}i_{0}=0$, which implies that, for sufficiently small $k_{3}$ or $i_{0}$,
the subspace $U:=\\{(s,c_{1},c_{2})\in\mathbb{R}^{3}:c_{2}=0\\}$ will be
nearly invariant. One perspective is to define the parameter ray by
$k_{3}=\varepsilon k_{3}^{*}$ and $e_{0}=\varepsilon e_{0}^{*}$. Then, the
perturbation form of the mass action equations is

$\begin{pmatrix}\dot{s}\\\ \dot{c}_{1}\\\
\dot{c}_{2}\end{pmatrix}=\begin{pmatrix}k_{1}s+k_{-1}&\;\;k_{1}s\\\
-k_{1}s-(k_{-1}+k_{2})&-k_{1}s\\\
0&-k_{-3}\end{pmatrix}\begin{pmatrix}c_{1}\\\
c_{2}\end{pmatrix}+\varepsilon\begin{pmatrix}-k_{1}e_{0}^{*}s\\\
k_{1}e_{0}^{*}s\\\
-k_{3}^{*}(c_{1}+c_{2})(i_{0}-c_{2})\end{pmatrix}+\mathcal{O}(\varepsilon^{2}).$

The QSS reduction is obtained by projecting the leading order perturbation
onto the critical manifold, thus

$\begin{pmatrix}\dot{s}\\\ \dot{c}_{1}\\\
\dot{c}_{2}\end{pmatrix}=\varepsilon\begin{pmatrix}1&\cfrac{k_{1}s+k_{-1}}{k_{1}s+k_{-1}+k_{2}}&\cfrac{k_{1}k_{2}s}{(k_{1}s+k_{-1}+k_{2})k_{-3}}\\\
0&0&0\\\ 0&0&0\end{pmatrix}\begin{pmatrix}-k_{1}e_{0}^{*}s\\\
\;\;k_{1}e_{0}^{*}s\\\ 0\end{pmatrix},$

from which we recover

$\dot{s}=-\cfrac{k_{1}k_{2}e_{0}s}{k_{1}s+k_{-1}+k_{2}},$

i.e., the sQSSA of the MM reaction mechanism. This is not surprising. With
initial conditions $s(0)=s_{0},c_{1}(0)=c_{2}(0)=0$, the dynamics are
approximately two-dimensional. From a different perspective (taking
independent limits), we can write (77) as

$\dot{s}=-\cfrac{k_{1}k_{2}e_{0}s}{k_{1}s+(k_{-1}+k_{2})\bigg{(}1+\cfrac{k_{3}i_{0}}{k_{-3}}\bigg{)}},$

from which it is clear by inspection that the sQSSA is recoverable from (77)
whenever $k_{3}i_{0}/k_{-3}\ll 1$. Consequently, whenever
$k_{3}i_{0}/k_{-3}\ll 1$ we need only consider the magnitude of
$\varepsilon_{I}$ to ascertain the accuracy of (77) (see, Figure 12).

Figure 12: Competitive inhibition reaction mechanism: If $k_{3}i_{0}/k_{-3}\ll
1$, the subspace $U$ is nearly invariant, and (77) will be accurate provided
$\varepsilon_{RS}\ll 1$. In both panels, the parameters (in arbitrary units)
are: $e_{0}=10.0$, $k_{1}=1.0$, $k_{2}=k_{-1}=10^{3}$, $k_{3}=10^{-7}$,
$k_{-3}=10^{-1}$ and $i_{0}=10^{0}$. It is straightforward to verify that (79)
is satisfied, and that $\varepsilon_{MM}=2.5\times 10^{-3}$, and
$\mu_{I}^{(1)}\approx 50$. The solid black curve is the numerical solution for
$s$ to the mass action system (75). The broken red curve is the numerical
solution to (77). Time has been mapped to the $\tau$ scale: $\tau=t/T$,
$\tau\in[0,1]$. Left panel: $s_{0}=10.0$ and the reduction (77) is very
accurate. Right panel: Here $s_{0}=10^{-1}$ and the reduction (77) is still
accurate; this confirms the long-time validity of (77) is regulated by the
magnitude of $\varepsilon_{MM}$ whenever $k_{3}i_{0}/k_{-3}\ll 1$.

##  7 Case Studies: Reduction from dimension three to two

In this section, we further discuss the uncompetitive and competitive
inhibition reaction mechanisms, but now we consider exemplary cases of
reduction to dimension two. These scenarios are of less practical relevance
than those in the previous section, but we present them for illustrative
purposes. We will provide less detailed discussions, and will be content to
show the feasibility of the method. The results from Section 5, in particular
Proposition 5 and Section 5.4, will be employed. The determination of
distinguished parameters now amounts to finding (or estimating) the maximum
and minimum of rational functions in two variables on some compact set.

###  7.1 Uncompetitive inhibition reaction mechanism

For the uncompetitive inhibition reaction mechanism, (64) and (65), one sees
that $k_{1}=k_{-3}=0$, with all other parameters contained in some compact
subset of the open positive orthant, defines a TFPV $\widehat{\pi}$, with a
two-dimensional critical manifold $\widetilde{Y}$ given by $c_{1}=0$. The TFPV
conditions mean that both elementary reactions responsible for the formation
of $C_{1}$ are slow.

We consider system (65) with initial values $s(0)=s_{0},\,c_{1}(0)=c_{2}(0)=0$
on the compact positively invariant set

$K:=\\{(s,c_{1},c_{2})\in\mathbb{R}^{3}_{\geq 0}:0\leq s\leq
s_{0},\;\;c_{1}+c_{2}\leq e_{0},\;\;c_{2}\leq\min\\{e_{0},i_{0}\\}\\},$

and take the ray direction

$\rho=(0,0,k_{1}^{*},0,0,0,k_{-3}^{*})^{\rm tr}$

in parameter space, with $k_{i}=\varepsilon k_{i}^{*}$, $k_{i}^{*}>0$, for
$i\in\\{1,\,-3\\}$. Straightforward computations yield the reduced system

$\begin{pmatrix}\dot{s}\\\
\dot{c}_{2}\end{pmatrix}=\dfrac{1}{k_{-1}+k_{2}+{k_{3}}(i_{0}-c_{2})}\begin{pmatrix}-k_{1}(e_{0}-c_{2})(k_{2}+k_{3}(i_{0}-c_{2}))s+k_{-3}k_{-1}c_{2}\\\
k_{1}k_{3}(e_{0}-c_{2})(i_{0}-c_{2})s-k_{-3}(k_{-1}+k_{2})c_{2}\end{pmatrix},$
(83)

with initial conditions $s(0)=s_{0},\,c_{2}(0)=0.$ A straightforward phase
plane analysis of system (83) (respectively, of the orbitally equivalent
system with the common denominator discarded) shows that every solution in the
positive quadrant converges to the stationary point $0$.

####  7.1.1 Asymptotic small parameters

We now determine dimensionless parameters that gauge the accuracy of (83). For
the sake of brevity, we will restrict attention to the case $e_{0}>i_{0}$. The
coefficients of the characteristic polynomial on $\widetilde{Y}$ are given by

$\begin{array}[]{rcl}\widetilde{\sigma}_{1}&=&k_{-1}+k_{2}+k_{3}(i_{0}-c_{2})+\varepsilon\,(\cdots)\\\
\widetilde{\sigma}_{2}&=&k_{1}\left(k_{3}(i_{0}-c_{2})(e_{0}-c_{2}+s)+k_{2}(e_{0}-c_{2})\right)+k_{-3}(k_{-1}+k_{2})+\varepsilon^{2}\,(\cdots)\\\
\widetilde{\sigma}_{3}&=&k_{1}k_{-3}k_{2}(e_{0}-c_{2}).\end{array}$

Thus,

$\begin{array}[]{rcl}\sigma_{1}&=&k_{-1}+k_{2}+k_{3}(i_{0}-c_{2})\\\
\widehat{\sigma}_{2}&=&k_{1}^{*}\left(k_{3}(i_{0}-c_{2})(e_{0}-c_{2}+s)+k_{2}(e_{0}-c_{2})\right)+k_{-3}^{*}(k_{-1}+k_{2})\\\
\widehat{\sigma}_{3}&=&k_{1}^{*}k_{-3}^{*}k_{2}(e_{0}-c_{2}),\end{array}$

and the first nondegeneracy condition from Lemma 7 is satisfied since
$e_{0}>i_{0}$.

According to Propositions 5 and 6 and their proofs, for timescale comparisons
we consider the rational function

$q(s,c_{2})=\dfrac{\widehat{\sigma}_{2}(s,c_{2}.\widehat{\pi},\rho,0)}{\sigma_{1}(s,c_{2},\widehat{\pi})^{2}},\quad
0\leq s\leq s_{0},\quad 0\leq c_{2}\leq i_{0}.$

Since $\widehat{\sigma}_{2}$ decreases with $c_{2}$ and increases with $s$,
while $\sigma_{1}$ decreases with $c_{2}$, we obtain an upper estimate from

$\varepsilon^{*}\leq\varepsilon\dfrac{\max\widehat{\sigma}_{2}}{(\min\sigma_{1})^{2}}=\dfrac{k_{1}(k_{3}i_{0}(e_{0}+s_{0})+k_{2}e_{0})+k_{-3}(k_{-1}+k_{2})}{(k_{-1}+k_{2})^{2}}=:\delta^{*}.$

Moreover from $\varepsilon^{*}\geq q(s_{0},0)$, we find that

$\dfrac{k_{1}(k_{3}i_{0}(e_{0}+s_{0})+k_{2}e_{0})+k_{-3}(k_{-1}+k_{2})}{(k_{-1}+k_{2}+k_{3}i_{0})^{2}}\leq\varepsilon^{*}.$

Likewise, we obtain lower timescale estimates from

$\varepsilon_{*}\geq\varepsilon\dfrac{\min\widehat{\sigma}_{2}}{(\max\sigma_{1})^{2}}=\dfrac{k_{1}k_{2}(e_{0}-i_{0})+k_{-3}(k_{-1}+k_{2})}{(k_{-1}+k_{2}+k_{3}i_{0})^{2}}=:\delta_{*}.$

Thus, for $i_{0}$ not too large, the estimates by $\delta^{*}$ and
$\delta_{*}$ are quite acceptable.

To estimate the disparity of the slow eigenvalues, according to Section 5.4,
we consider

$\kappa^{*}=\max\dfrac{\sigma_{1}\widehat{\sigma}_{3}}{\widehat{\sigma}_{2}^{2}}\leq\dfrac{(k_{-1}+k_{2}+k_{3}i_{0})k_{1}^{*}k_{-3}^{*}k_{2}e_{0}}{(k_{1}^{*}k_{2}e_{0}+k_{-3}^{*}(k_{-1}+k_{2}))^{2}}=:\nu^{*}$

as well as

$\kappa_{*}=\min\dfrac{\sigma_{1}\widehat{\sigma}_{3}}{\widehat{\sigma}_{2}^{2}}\geq\dfrac{(k_{-1}+k_{2})k_{1}^{*}k_{-3}^{*}k_{2}(e_{0}-i_{0})}{(k_{1}^{*}(k_{3}i_{0}(e_{0}+s_{0})+k_{2}e_{0})+k_{-3}^{*}(k_{-1}+k_{2}))^{2}}=:\nu_{*}.$

Whenever $i_{0}$ is not too large, these two parameters are close, and so are
the slow eigenvalues.

####  7.1.2 Numerical Simulations

From Fenichel theory, it is known that the accuracy of the reduction (83)
improves along the perturbation direction as $\varepsilon\to 0$. We include
some numerical simulations to gauge the efficacy of the parameter
$\delta^{*}$:

  1. 1.

Following the outline established in Section 6, we first consider a case when
all parameters are of unit order. Numerical simulations confirm that the
accuracy of (83) improves as $\delta^{*}\to 0$ along the parameter ray
direction (see, Figure 13). We include the values of $\delta_{*}$ to indicate
the variation of timescale ratios.

Figure 13: Uncompetitive inhibition reaction mechanism with reduction to
dimension two. In all panels, the parameters (in arbitrary units) are:
$e_{0}=2.0$, $s_{0}=1.0$, $k_{2}=k_{-1}=1.0$, $k_{3}=1.0$, and $i_{0}=1.0$.
The solid black curve is the numerical solution to the mass action system
(65). The dashed/dotted red curve is the numerical solution to (83). Time has
been mapped to the $\tau$ scale: $\tau=t/T$, $\tau\in[0,1]$. Top panels:
$k_{1}=k_{-3}=10^{-1}$ with $\delta^{*}=1.75\times 10^{-1}$ and
$\delta_{*}=3.33\times 10^{-2}$. Bottom panels: $k_{1}=k_{-3}=10^{-3}$ with
$\delta^{*}=1.75\times 10^{-3}$ and $\delta_{*}=3.33\times 10^{-4}$. As
expected, the accuracy of (83) improves as the perturbation decreases along
the parameter ray.

  2. 2.

In a second set of simulations, we consider the case of parameter values that
are disparate in magnitude. Numerical simulations confirm once again that the
accuracy of (83) improves along the parameter ray direction as $\delta^{*}\to
0$ (see, Figure 14).

Figure 14: Uncompetitive inhibition reaction mechanism with reduction to
dimension two and disparate parameter values. In all panels, the parameters
(in arbitrary units) are: $e_{0}=10.0$, $s_{0}=5.0$, $k_{2}=10^{3}$,
$k_{-1}=10^{2}$, $k_{3}=10.0$, and $i_{0}=1.0$. The solid black curve is the
numerical solution to the mass action system (65). The broken red curve is the
numerical solution to (83). Time has been mapped to the $\tau$ scale:
$\tau=t/T$, $\tau\in[0,1]$. Top panels: $k_{1}=k_{-3}=10^{2}$ with
$\delta^{*}\approx 9.3\times 10^{-1}$ and $\delta_{*}\approx 8.2\times
10^{-1}$. Bottom panels: $k_{1}=k_{-3}=1.0$ with $\delta^{*}\approx 9.3\times
10^{-3}$ and $\delta_{*}\approx 8.3\times 10^{-3}$. Again, it is clear that
the accuracy of (83) improves as the perturbation decreases along the
parameter ray. Notably, the approximation in the second case is very good
although $k_{1}=k_{-3}=1$. (As always, the expression “$\ll 1$” should not be
taken too literally.) As for measuring the discrepancy of the “slow”
eigenvalues, one finds $\kappa_{*}\geq\nu_{*}\approx 7.8\cdot 10^{-2}$ and
$\kappa^{*}\leq\nu^{*}\approx 9.0\cdot 10^{-2}$. This indicates a ratio of
about $10^{-1}$.

###  7.2 Competitive inhibition reaction mechanism

For the competitive inhibition reaction mechanism, (74) and (75), we consider
the case that formation of complex $C_{1}$, and both formation and degradation
of complex $C_{2}$, are slow. Setting $k_{1}=k_{3}=k_{-3}=0$, with all the
other parameters contained in a compact subset of the positive orthant,
defines a TFPV $\widehat{\pi}$ for dimension $s=2$, the critical manifold
$\widetilde{Y}$ being given by $c_{1}=0$ (see, Kruff and Walcher [25]). We
consider the system on the compact positively invariant set $K$ defined by

$K:=\\{(s,c_{1},c_{2})\in\mathbb{R}_{\geq 0}^{3}:0\leq s\leq s_{0},\;\;0\leq
c_{1}\leq e_{0},\;\;0\leq c_{2}\leq\min\\{e_{0},\,i_{0}\\}\\},$

choosing the ray direction

$\rho=(0,0,k_{1}^{*},0,0,k_{3}^{*},k_{-3}^{*})^{\rm tr}$

in parameter space, and $k_{i}=\varepsilon k_{i}^{*}$ with $k_{i}^{*}>0$ for
$i\in\\{1,\,3,\,-3\\}$. Standard computations yield the reduced system

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

with initial conditions $s(0)=s_{0}$, $c_{2}(0)=0$. The qualitative behavior
of this system is easily determined. All solutions in the positive quadrant
converge to a stationary point $(0,c_{2}^{*})$, with
$0<c_{2}^{*}<\min\\{e_{0},i_{0}\\}$.212121Incidentally, it is possible to
compute the solutions to (84) via quadratures. The second equation is
separable, and upon substitution the first equation is non-autonomous linear.

####  7.2.1 Asymptotic small parameters

For the sake of brevity, we will consider only the case $e_{0}>i_{0}$. The
coefficients of the characteristic polynomial on $\widetilde{Y}$ are

$\begin{array}[]{rcl}\widetilde{\sigma}_{1}&=&k_{-1}+k_{2}+\varepsilon\,(\cdots)\\\
\widetilde{\sigma}_{2}&=&k_{1}k_{2}(e_{0}-c_{2})+k_{3}(e_{0}+i_{0}-2c_{2})(k_{-1}+k_{2})+k_{-3}(k_{-1}+k_{2})+\varepsilon^{2}\,(\cdots)\\\
\widetilde{\sigma}_{3}&=&k_{1}k_{3}\cdot\left(k_{2}e_{0}(e_{0}+i_{0}-2c_{2})+k_{2}c_{2}(2c_{2}-e_{0}-i_{0})\right)+k_{1}k_{-3}\cdot
k_{2}(e_{0}-c_{2}),\end{array}$

and we obtain

$\begin{array}[]{rcl}\sigma_{1}&=&k_{-1}+k_{2}\\\
\widehat{\sigma}_{2}&=&k_{1}^{*}k_{2}(e_{0}-c_{2})+k_{3}^{*}(e_{0}+i_{0}-2c_{2})(k_{-1}+k_{2})+k_{-3}^{*}(k_{-1}+k_{2})\\\
\widehat{\sigma}_{3}&=&k_{1}^{*}k_{2}(e_{0}-c_{2})\left(k_{3}^{*}(e_{0}+i_{0}-2c_{2})+k_{-3}^{*}\right).\end{array}$

The nondegeneracy conditions are satisfied (also at $c_{2}=i_{0}$), due to
$e_{0}>i_{0}$. As for timescales, we need to analyze the rational function

$q(s,c_{2})=\dfrac{\widehat{\sigma}_{2}(s,c_{2},\widehat{\pi},\rho,0)}{\sigma_{1}(s,c_{2},\widehat{\pi})^{2}},\quad
0\leq s\leq s_{0},\quad 0\leq c_{1}\leq e_{0},\quad 0\leq c_{2}\leq i_{0}.$

Since $\sigma_{1}$ is constant and $\widehat{\sigma}_{2}$ is decreasing with
$c_{2}$, attaining its maximum at $c_{2}=0$, and its minimum at $c_{2}=i_{0}$,
we find the distinguished parameters

$\varepsilon^{*}=\dfrac{k_{1}k_{2}e_{0}+(k_{3}(e_{0}+i_{0})+k_{-3})(k_{-1}+k_{2})}{(k_{-1}+k_{2})^{2}}$

and

$\varepsilon_{*}=\dfrac{k_{1}k_{2}(e_{0}-i_{0})+(k_{3}(e_{0}-i_{0})+k_{-3})(k_{-1}+k_{2})}{(k_{-1}+k_{2})^{2}}.$

Furthermore, according to Section 5.4, we consider

$\kappa^{*}=\max\dfrac{\sigma_{1}\widehat{\sigma}_{3}}{\widehat{\sigma}_{2}^{2}}\leq\dfrac{k_{1}^{*}k_{2}e_{0}(k_{3}^{*}(e_{0}+i_{0})+k_{-3}^{*})(k_{-1}+k_{2})}{\left(k_{1}^{*}k_{2}(e_{0}-i_{0})+(k_{3}^{*}(e_{0}-i_{0})+k_{-3}^{*})(k_{-1}+k_{2})\right)^{2}}=:\nu^{*}$
(85)

and

$\kappa_{*}=\min\dfrac{\sigma_{1}\widehat{\sigma}_{3}}{\widehat{\sigma}_{2}^{2}}\geq\dfrac{k_{1}^{*}k_{2}(e_{0}-i_{0})(k_{3}^{*}(e_{0}-i_{0})+k_{-3}^{*})(k_{-1}+k_{2})}{\left(k_{1}^{*}k_{2}e_{0}+(k_{3}^{*}(e_{0}+i_{0})+k_{-3}^{*})(k_{-1}+k_{2})\right)^{2}}=:\nu_{*}$
(86)

to measure the disparity between the eigenvalues $\lambda_{2}$ and
$\lambda_{3}$.

####  7.2.2 Numerical simulations

We present numerical examples to gauge the accuracy of the reduction (84) with
decreasing $\varepsilon^{*}$:

  1. 1.

For our first example, we consider the case with
$\pi:=\varepsilon(k_{1}^{*},1.0,1.0,1.0,1.0.,k_{3}^{*},k_{-3}^{*},1.0)^{\text{tr.}}$
(see, Figure 15). We include the values of $\varepsilon^{*},\,\varepsilon_{*}$
to indicate the variation of timescale ratios.

Figure 15: Competitive inhibition reaction mechanism with reduction to
dimension two. The accuracy of (84) improves along the parameter ray. In all
panels, the parameters (in arbitrary units) are: $e_{0}=1.0$, $s_{0}=1.0$,
$k_{2}=1.0$, $k_{-1}=1.0$, and $i_{0}=1.0$. The solid black curve is the
numerical solution to the mass action system (75). The broken red curve is the
numerical solution to (84). Time has been mapped to the $\tau$ scale:
$\tau=t/T$, $\tau\in[0,1]$. Top panels: Simulation performed with
$k_{1}=k_{-3}=k_{3}=10^{-1}$ and $\varepsilon^{*}=1.75\times 10^{-1}$,
$\varepsilon_{*}=5.0\times 10^{-2}$. Bottom panels: Simulation performed with
$k_{1}=k_{-3}=k_{3}=10^{-3}$ and $\varepsilon^{*}=1.75\times 10^{-3}$,
$\varepsilon_{*}=5.0\times 10^{-4}$. The singular perturbation reduction (84)
is practically indistinguishable from (75).

  2. 2.

For our second example, we again consider a case when parameters are of
differing magnitudes. Once more, numerical simulations confirm that the QSS
reduction (84) improves as $\varepsilon^{*}\to 0$ along the parameter ray
(see, Figure 16).

Figure 16: Competitive inhibition reaction mechanism with reduction to
dimension two and disparate parameters. In all panels, the “non-small”
parameters (in arbitrary units) are: $e_{0}=10^{2}$, $s_{0}=10^{2}$,
$k_{2}=10^{2}$, $k_{-1}=1.0$, and $i_{0}=10.0$. The solid black curve is the
numerical solution to the mass action system (75). The broken red curve is the
numerical solution to (84). Time has been mapped to the $\tau$ scale:
$\tau=t/T$, $\tau\in[0,1]$. Top panels: Simulation performed with
$k_{1}=k_{-3}=k_{3}=1.0$ and $\varepsilon^{*}\approx 2.08$,
$\varepsilon_{*}\approx 1.78$. This scenario is outside the range of
applicability for Proposition 6. Bottom panels: Simulation performed with
$k_{1}=k_{-3}=k_{3}=10^{-2}$ and $\varepsilon^{*}\approx 2.08\times 10^{-2}$,
$\varepsilon_{*}\approx 1.78\times 10^{-2}$. The singular perturbation
reduction (84) here is very close to (75).

####  7.2.3 The case of very small $k_{1}$: Three timescales

Finally, we discuss a scenario mentioned in Section 5.4. From equations (85)
and (86), below one sees that both $\kappa^{*}$ and $\kappa_{*}$ approach zero
as $k_{1}^{*}\to 0$. This may indicate three timescales. Moreover, from
equation (75), one sees that the plane defined by $c_{1}=0$ is invariant when
$k_{1}=0$, thus nearly invariant when $k_{1}$ is small. A coordinate-
independent approach to a three-timescale scenario was presented in Kruff and
Walcher [25], based on work of Cardin and Teixeira [6]. We introduce two small
parameters $\varepsilon_{1},\,\varepsilon_{2}$ and

$k_{3}=\varepsilon_{1}k_{3}^{\dagger},\,k_{-3}=\varepsilon_{1}k_{-3}^{\dagger},\,k_{1}=\varepsilon_{1}\varepsilon_{2}k_{1}^{\dagger},$

and rewrite system (75) with three timescales. As detailed in [25], the system
admits a sequence of two reductions, with nested invariant manifolds:

A reduction to slow dynamics on a two-dimensional invariant manifold close to
$c_{1}=0$, with reduced system

$\begin{pmatrix}\dfrac{ds}{d\tau_{1}}\\\
\dfrac{dc_{2}}{d\tau_{1}}\end{pmatrix}=\begin{pmatrix}0\\\
k_{3}^{\dagger}(e_{0}-c_{2})(i_{0}-c_{2})-k_{-3}^{\dagger}c_{2}\end{pmatrix}$
(87)

with $\tau_{1}=\varepsilon_{1}t$.

A subsequent reduction to “very slow” dynamics on a one-dimensional invariant
manifold close to $c_{1}=0,c_{2}=\widetilde{c}_{2}$, with

$\widetilde{c}_{2}=\cfrac{k_{3}(e_{0}+i_{0})+k_{-3}-\sqrt{(k_{3}(e_{0}+i_{0})+k_{-3})^{2}-4e_{0}i_{0}k_{3}^{2}}}{2k_{3}}.$

The fully reduced one-dimensional equation is then

$\dfrac{ds}{d\tau_{2}}=-\cfrac{k_{2}\cdot
k_{1}^{\dagger}(e_{0}-\widetilde{c}_{2})}{k_{-2}+k_{-1}}s$

with $\tau_{2}=\varepsilon_{1}\varepsilon_{2}t$, or, restated in fast time,

$\dot{s}=-\cfrac{k_{2}\cdot k_{1}(e_{0}-\widetilde{c}_{2})}{k_{-2}+k_{-1}}s.$
(88)

Figure 17 illustrates the “slow–very slow” dynamics for a numerical example.

Figure 17: Competitive inhibition reaction mechanism with a three timescale
scenario. In both panels, the parameters (in arbitrary units) are:
$s_{0}=10.0$, $e_{0}=10^{2}$, $k_{1}=1\times 10^{-4}$, $k_{2}=2\times 10^{3}$,
$k_{-1}=1.0$ ,$k_{3}=10^{-1}$, $k_{-3}=10^{-3}$ and $i_{0}=50.0$. The solid
black curve is the numerical solution to the mass action system (75). The
dashed/dotted red curves are the numerical solutions to (87) and (88). Time
has been mapped to the $\tau$ scale: $\tau=t/T$. For the chosen parameter
values, $\varepsilon^{*}\approx 7.5\times 10^{-3}$ and $\varepsilon_{*}\approx
2.5\times 10^{-3}$. Moreover, the two slow eigenvalues are disparate since
$\nu^{*}\approx 6\times 10^{-3}$, this is consistent with a three timescale
scenario. Top Left panel: The initial accumulation of $c_{1}$ occurs on the
fast timescale; the concentrations of $c_{2}$ and $s$ are approximately
constant on the fast timescale. Top Right panel: The reduction (87) is
accurate on the slow timescale as $c_{2}$ approaches its threshold value,
$\widetilde{c}_{2}$. Bottom Left panel: The reduction (88) is accurate on the
the very slow timescale, $\tau_{2}$, on which the depletion of $s$ is
significant. Bottom Right panel: On the very slow timescale, $c_{2}$ is
effectively constant: $c_{2}\approx\widetilde{c}_{2}$.

##  8 Discussion

While the underlying theory and the qualitative analysis concerning the
reduction of biochemical and chemical reaction networks is well understood and
rests on solid ground, there is a sizable gap between available theory and
applications to parameter identification problems in laboratory settings,
where heuristics and ad hoc approaches are (perforce) still prevalent. Closing
the gap requires further, more precise theoretical results. The present paper
contributes towards this goal, by introducing a general consistent method to
obtain perturbation parameters, based on local linear timescales. Note that by
its nature, our approach is focused on and limited to the local behavior.

We briefly recall the context and reviewing the results of the present paper:

  1. 1.

We start from a singular perturbation reduction with a well-defined critical
manifold. This is crucial to ensure appropriateness of linearizability.
Considering the three steps (as outlined in the Introduction) that are
necessary for a global quantitative estimate of the approximation error, our
results amount to an essential part of Step 1. In absence of results
concerning Steps 2 and 3, direct applications are limited. But, our results
permit consistency checks, which show that certain common perturbation
parameters are not feasible.

  2. 2.

Using classical results from algebra to approximate eigenvalue ratios in the
asymptotic limit, we obtained parameters that are computable, palatable, and
admit a biochemical interpretation.

  3. 3.

We first applied our methods to the Michaelis–Menten reaction mechanism. As it
turns out, even for such a familiar system our approach provides new and
elucidating perturbation parameters. Moreover, we included a partial
discussion of Step 3 for the irreversible system with small product formation
rate.

  4. 4.

For two relevant extensions of the Michaelis–Menten reaction mechanism (like
the uncompetitive inhibition and competitive inhibition), and a non-
Michaelis–Menten reaction mechanism (like the cooperative system with two
complexes), we derived perturbation parameters in the spirit of Segel and
Slemrod [40], but without resorting to nonlinear timescales. This stands in
contrast to the practice of using $\varepsilon_{BH}$ or $\varepsilon_{SSl}$,
or ad-hoc modifications of these. We augmented these results by an extensive
discussion of numerical examples to illustrate the efficacy of these
parameters, but also to highlight the importance of the compactness
requirements we impose throughout. We also discussed one case that leads to a
system with three timescales.

  5. 5.

Finally, we discussed exemplary cases of reduction from dimension three to
dimension two for both reaction inhibition scenarios, to verify the
feasibility of our approach. Numerical simulations illustrate the quality and
accuracy of the approximations.

The remaining items (Step 2 and Step 3) as stated in the Introduction need to
be handled on a case-by-case basis. We will provide a complete analysis of the
irreversible Michaelis–Menten reaction mechanism with low enzyme in
forthcoming work.

##  9 Appendix

In this section, we collect some technical matters and proofs, as well as
recalling some known results for which a concise presentation seems
appropriate and useful.

###  9.1 Lyapunov function arguments

Lyapunov functions can be used to estimate the approach to the slow manifold
in a singularly perturbed system, as was mentioned in the Introduction. This
estimate gives rise to a small parameter $\varepsilon_{L}$ which controls the
distance of the solution to the slow manifold. We give an account of the
relevant facts here.

We first state an auxiliary result that goes back to Lyapunov.

######  Lemma 5.

Let $Q$ be a real $n\times n$-matrix, with eigenvalues
$\mu_{1},\ldots,\mu_{n}$, and let $\delta>0$. Then there exists a scalar
product $\left<\cdot,\cdot\right>$ on $\mathbb{R}^{n}$ such that for all $x$
one has

$\left(\min_{1\leq i\leq n}{\rm
Re}\,\mu_{i}-\delta\right)\,\left<x,\,x\right>\leq\left<x,\,Qx\right>\leq\left(\max_{1\leq
i\leq n}{\rm Re}\,\mu_{i}+\delta\right)\,\left<x,\,x\right>,$

and

$\left(\min_{1\leq i\leq
n}|\mu_{i}|^{2}-\delta\right)\,\left<x,\,x\right>\leq\left<Qx,\,Qx\right>\leq\left(\max_{1\leq
i\leq n}|\mu_{i}|^{2}+\delta\right)\,\left<x,\,x\right>.$

This can be proven as in Walter [46, Chapter VII, §30] (see, also Arnold [2,
Chapter 22]). For matrices that are diagonalizable over $\mathbb{C}$ build a
real basis from real and imaginary parts of a complex eigenbasis. For the non-
diagonalizable case, by suitable choice of basis elements the nilpotent part
can be chosen to have norm $<\delta$.

####  9.1.1 Estimates

This presentation follows Berglund and Gentz [3, Section 2.1 ff.], but for
illustrative purposes, we are satisfied with a local version. Consider a
smooth system

$\begin{array}[]{rcl}\dot{x}&=&\varepsilon\widetilde{f}_{1}(x,y,\varepsilon)\\\
\dot{y}&=&f_{2}(x,y,\varepsilon)\end{array},\quad\text{briefly}\begin{pmatrix}\dot{x}\\\
\dot{y}\end{pmatrix}=F(\begin{pmatrix}x\\\ y\end{pmatrix})$ (89)

with $\begin{pmatrix}x\\\ y\end{pmatrix}$ in some open subset of
$\mathbb{R}^{n}$, $x\in\mathbb{R}^{m}$, and a nonnegative parameter
$\varepsilon$. Moreover let $\begin{pmatrix}x_{0}\\\ y_{0}\end{pmatrix}$ be
such that $f_{2}(x_{0},y_{0},0)=0$, and $M$ a suitable compact neighborhood of
this point. (More conditions on $M$ will be implicitly imposed below, by
further assumptions.)

  * •

Assume furthermore that

$f_{2}(x,y,\varepsilon)=0\Longleftrightarrow y=g(x,\varepsilon)$

for $\begin{pmatrix}x\\\ y\end{pmatrix}\in M$ and
$\varepsilon\leq\varepsilon_{\rm max}$, with some positive $\varepsilon_{\rm
max}$ and a smooth function $g$. The zero set $Y_{\varepsilon}$ of
$f_{2}(\cdot,\cdot,\varepsilon)$ in $M$ will be called the slow manifold, or
QSS manifold,222222This is an order $\varepsilon$ approximation of the slow
manifold in a singular perturbation setting. for $\varepsilon$. By Hadamard’s
lemma, after possibly shrinking $M$ there exists a smooth matrix valued
function $A$ such that

$f_{2}(x,y,\varepsilon)=A(x,y,\varepsilon)\cdot(y-g(x,\varepsilon)).$

Thus, we may rewrite system (89) as

$\begin{array}[]{rcl}\dot{x}&=&\varepsilon\widetilde{f}_{1}(x,y,\varepsilon)\\\
\dot{y}&=&A(x,y,\varepsilon)\cdot(y-g(x,\varepsilon)).\end{array}$ (90)

With

$D_{y}f_{2}(x,y,\varepsilon)=A(x,y,\varepsilon)+\left(D_{y}A(x,y,\varepsilon)\,\right)(y-g(x,\varepsilon)),$

one finds in particular

$D_{y}f_{2}(x,y,\varepsilon)=A(x,y,\varepsilon)\text{ on }Y_{\varepsilon}.$

  * •

Now assume that all eigenvalues of $A(x_{0},y_{0},0)$ have negative real
parts. By continuity and suitable choice of $M$ and $\varepsilon_{\rm max}$,
all eigenvalues of $A(x,y,\varepsilon)$ have negative real part for
$\begin{pmatrix}x\\\ y\end{pmatrix}\in M$ and
$0\leq\varepsilon\leq\varepsilon_{\rm max}$. Due to Lemma 5, there exists a
scalar product $\left<\cdot,\cdot\right>$ on $\mathbb{R}^{n-m}$ and some
$\gamma>0$ such that

$\left<z,A(x_{0},y_{0},0)z\right>\leq-2\gamma\left<z,z\right>,\quad\text{ all
}z\in\mathbb{R}^{n-m}.$

(Recall the correspondence between $2\gamma$ and eigenvalues.) Thus, we may
assume that

$\left<z,A(x,y,\varepsilon)z\right>\leq-\gamma\left<z,z\right>,\quad\text{ all
}z\in\mathbb{R}^{n-m},$ (91)

on $M$, with $0\leq\varepsilon\leq\varepsilon_{\rm max}$. Denote by
$\|\cdot\|$ the norm associated with this scalar product.

The following line of arguments is a slight variant of classical reasoning
(which uses Gronwall’s lemma, see e.g. Evans [15, Appendix B for the latter]).
For solutions of (89) we find

$\begin{array}[]{rcl}\dfrac{d}{dt}\left<y-g(x,\varepsilon),y-g(x,\varepsilon)\right>&=&2\left<y-g(x,\varepsilon),\dot{y}-D_{x}g(x,\varepsilon)f_{1}(x,y,\varepsilon)\right>\\\
&=&2\left<y-g(x,\varepsilon),A(x,y,\varepsilon)\left(y-g(x,\varepsilon)\right)\right>\\\
&&\quad-2\left<y-g(x,\varepsilon),D_{x}g(x,\varepsilon)f_{1}(x,y,\varepsilon)\right>.\end{array}$

The first term on the right hand side can be estimated by
$-\gamma\cdot\left<y-g(x,\varepsilon),y-g(x,\varepsilon)\right>$. As for the
second term, by Cauchy-Schwarz one has

$\begin{array}[]{rcl}2\left|\left<y-g(x,\varepsilon),D_{x}g(x,\varepsilon)f_{1}(x,y,\varepsilon)\right>\right|&\leq&2\|y-g(x,\varepsilon)\|_{2}\cdot\|D_{x}g(x,\varepsilon)f_{1}(x,y,\varepsilon)\|_{2}\\\
&\leq&2\|y-g(x,\varepsilon)\|_{2}\cdot\left(\|D_{x}g(x,\varepsilon)\|\cdot\|f_{1}(x,y,\varepsilon)\|\right)\end{array}$

with suitable norms in the second and third factor.

Now, there exists a positive constant $\kappa=\varepsilon\widetilde{\kappa}$
such that

$\|D_{x}g(x,\varepsilon)\|\cdot\|f_{1}(x,y,\varepsilon)\|\leq\kappa.$

So, for $V:=\|y-g(x,\varepsilon)\|^{2}$ one obtains the differential
inequality

$\frac{dV}{dt}\leq-2\gamma V+2\kappa\sqrt{V}.$

Comparison with the solution of the corresponding Bernoulli equation yields

$V\leq V(0)\exp(-\gamma
t)+\left(\frac{\kappa}{\gamma}\right)^{2}\cdot(1-\exp(-\gamma t)),$ (92)

thus $\|y-g(x,\varepsilon)\|=\sqrt{V(t)}$ can be estimated e.g. by
$\sqrt{2}\frac{\kappa}{\gamma}$ as $t\to\infty$.232323One may replace
$\sqrt{2}$ by any smaller constant which is $>1$. Therefore, after a transient
phase the proximity of the solution to the slow manifold is controlled by

$\varepsilon_{L}:=\sqrt{2}\frac{\kappa}{\gamma}=\sqrt{2}\varepsilon\frac{\widetilde{\kappa}}{\gamma}.$
(93)

More precisely, once $V(0)\exp(-\gamma
t)\leq\left(\dfrac{\kappa}{\gamma}\right)^{2}$, the stated estimate holds. The
inequality is satisfied whenever

$t\geq\frac{1}{\gamma}\,\log\left(\dfrac{\gamma^{2}V(0)}{\kappa^{2}}\right)\sim\log\dfrac{1}{\varepsilon_{L}},$

and this indicates that the time span for the approach to the QSS manifold is
of order $|\log\varepsilon_{L}|$ in the fast timescale, and of order
$\varepsilon_{L}|\log\varepsilon_{L}|$ in the slow timescale
$\varepsilon_{L}t.$ (A more detailed analysis will provide a lower estimate by
a variant of (91), and confirm that the asymptotic estimate cannot be
improved.) In particular time spans of order $1$ will not suffice for the
transient.

In reaction network settings, $\varepsilon_{L}$ is a dimensional parameter
(with dimension concentration); a suitable normalization needs to be chosen.

####  9.1.2 A correspondence to eigenvalues

We sketch the relation of the small parameter $\varepsilon_{L}$ to eigenvalues
of the Jacobian. For the sake of simplicity, we only consider the
linearization here, disregarding higher order terms. Given the system

$\begin{array}[]{rcccccc}\dot{x}&=&-\varepsilon\widetilde{U}x&+&\varepsilon\widetilde{V}y&&\\\
\dot{y}&=&Wx&-&Zy&=&-Z\left(y-Z^{-1}W\right)x\\\ \end{array}\quad,\text{
briefly}\begin{pmatrix}\dot{x}\\\ \dot{y}\end{pmatrix}=F(x,y,\varepsilon),$

and keeping the notation from above, we have $A=-Z$, $g(x)=Z^{-1}Wx$,
$D_{x}g=Z^{-1}W$. The slow manifold $Y_{\varepsilon}$ is given by $Wx-Zy=0$,
up to higher order terms. Moreover

$\widetilde{f}_{1}=-\widetilde{U}x+\widetilde{V}y=\left(-\widetilde{U}+\widetilde{V}Z^{-1}W\right)x\quad\text{
on }Y_{\varepsilon}.$

Now consider the eigenvalues of the matrix
$DF=\begin{pmatrix}-\varepsilon\widetilde{U}&\varepsilon\widetilde{V}\\\
W&-Z\end{pmatrix}$; see also Lemma 7. Thus, let
$\alpha_{0}+\varepsilon\alpha_{1}+\cdots$ be an eigenvalue with eigenvector
$\begin{pmatrix}x_{0}+\varepsilon x_{1}+\cdots\\\ y_{0}+\varepsilon
y_{1}+\cdots\end{pmatrix}$; $\begin{pmatrix}x_{0}\\\
y_{0}\end{pmatrix}\not=0$. For $\alpha_{0}\not=0$, comparing lowest order
terms in the eigenvalue condition yields

$x_{0}=0\text{ and }Wx_{0}-Zy_{0}=\alpha_{0}y_{0},$

thus $-\alpha_{0}$ is an eigenvalue of $Z$. By Lemma 5, we see that $2\gamma$
can be chosen near the nonzero eigenvalue of $DF(x,y,0)$ with smallest
absolute real part.

For $\alpha_{0}=0$, thus the eigenvalue has order $\varepsilon$, comparing
lowest orders in the eigenvalue condition yields

$-\widetilde{U}x_{0}+\widetilde{V}y_{0}=\alpha_{1}x_{0}\text{ and
}Wx_{0}-Zy_{0}=0,$

hence $\alpha_{1}$ is an eigenvalue for
$-\widetilde{U}+\widetilde{V}Z^{-1}W=\widetilde{f}_{1}$. An upper estimate for
$\varepsilon\|\widetilde{f}_{1}\|$ can be obtained from Lemma 5: Choose the
order $\varepsilon$ eigenvalue with greatest absolute value, multiplied by
some factor accounting for a coordinate change. Thus, we see that $\kappa$ is
composed of the factor $\|Z^{-1}W\|$ (which reflects the geometry of the slow
manifold), the absolutely largest eigenvalue of order $\varepsilon$ and some
multiplicative constants from coordinate transformations. In our local
setting, all the multiplicative constants mentioned above are of order one.

To summarize, the small parameter $\varepsilon_{L}=\kappa/\gamma$ is
determined by the ratio of the largest absolute eigenvalue of order
$\varepsilon$ to the smallest absolute real part of eigenvalues of order one.
From this perspective, for slow manifolds of dimension one in particular, the
relevance of the parameters $\varepsilon^{*}$ and $\mu^{*}$ is obvious. Their
advantage lies in their (relative) computational accessibility. Likewise,
$\varepsilon^{*}$ is a both relevant and computationally accessible parameter
for three-dimensional systems with two-dimensional slow manifolds.

####  9.1.3 Remarks on Steps 2 and 3

Lyapunov function arguments provide a small parameter $\varepsilon_{L}$ which
characterizes closeness of a solution of (89) to the slow manifold. This takes
care of Step 1 described in the Introduction, and clarifies the role of
eigenvalues up to ($\varepsilon$-independent) factors due to coordinate
changes.

For the ultimate goal of obtaining quantitative estimates for the discrepancy
between the true solution and the singular perturbation approximation, one
needs to go further. In Step 2, an appropriate critical time for the onset of
the slow dynamics, as well as an appropriate initial value for the reduced
system, must be determined. As for Step 3, by a continuity and compactness
argument, the right hand sides of the full and the reduced equation differ by
$\varepsilon_{L}$ times some constant. With this, and an error estimate for
the initial value for the reduced system, continuous dependence provides an
estimate of the approximation error on compact time intervals. Further work
may be required, since one is mostly interested in unbounded time intervals,
so one cannot rely only on standard continuous dependence theorems.

In the present manuscript, we generally did not address the determination of
$\varepsilon_{L}$ in examples and case studies. The only exception is
irreversible Michaelis–Menten with slow product formation (see, Section 4.1),
which also contains partial results for Step 3. For the (more familiar and
more relevant) irreversible Michaelis–Menten system with small enzyme
concentration all three steps can dealt with completely (even if some
complications arise), as will be shown in a forthcoming paper. For any system
of dimension $>2$, even completing Step 1 seems quite demanding.

###  9.2 A proof of Lemma 2

###### Proof.

Part (a) is a special case of Lemma 7 below. To prove part (b), abbreviate
$\sigma^{*}_{i}(x):=\sigma_{i}(x,\widehat{\pi})$ for $x\in\widetilde{Y}\cap
K$, $1\leq i\leq n-1$. Then the nonzero roots of the characteristic polynomial
$\chi$ are the roots of

$\zeta(x,\tau):=\tau^{n-1}+\sigma^{*}_{1}(x)\tau^{n-2}+\cdots+\sigma^{*}_{n-1}(x).$

By the blanket assumptions, the $\sigma^{*}_{i}$ are bounded above and below
by positive constants, hence the absolute values of all zeros of the
$\zeta(x,\cdot)$ are bounded above by some constant. Since $\widehat{\pi}$ is
a TFPV, all zeros have negative real parts. Now assume that for every positive
constant $\delta$, some $\zeta(x,\tau)$ has a zero with real part
$\geq-\delta$. Then there exist sequences $(x_{k})$ in $\widetilde{Y}\cap K$
and $(\mu_{k})$ in $\mathbb{C}$ such that $\zeta(x_{k},\mu_{k})=0$ and ${\rm
Re}\,\mu_{k}\to 0$. Due to boundedness of the sequence $(\mu_{k})$ and
compactness of $\widetilde{Y}\cap K$ we may assume that the $\mu_{k}$ converge
to $\mu^{*}$, ${\rm Re}\,\mu^{*}=0$, and the $x_{k}$ converge to
$x^{*}\in\widetilde{Y}\cap K$. By continuity $\zeta(x^{*},\mu^{*})=0$; a
contradiction. Part (c) follows by continuity and compactness arguments. ∎

###  9.3 Parameter dependence of eigenvalues

Recall from (39) the definition

$\widetilde{\sigma}_{i}(x,\varepsilon):=\sigma_{i}(x,\widehat{\pi}+\varepsilon\rho),\quad
1\leq i\leq n,\quad\widetilde{\sigma}_{0}:=1.$

We first prove (44), concerning the orders of the $\widetilde{\sigma}_{i}$
whenever $s>1$.

######  Lemma 6.

Let $\widehat{\pi}$ be a TFPV for dimension $s$, with critical manifold
$\widetilde{Y}$. Then for all $x\in\widetilde{Y}\cap K$ one has

$\widetilde{\sigma}_{i}(x,\varepsilon)=\varepsilon^{i-n+s}\widehat{\sigma}_{i}(x,\varepsilon)\text{
for all }x\in\widetilde{Y}\cap K,\quad n-s\leq i\leq n,$

with polynomial $\widehat{\sigma}_{i}$.

###### Proof.

The arguments we will use are similar to those in the proof of Goeke et al.
[21, Proposition 3]. We set

$\widetilde{h}(x,\varepsilon):=h(x,\widehat{\pi}+\varepsilon\rho)\text{ for
}x\in\widetilde{Y}\cap K.$

There exists a local transformation of $\widetilde{h}$ into Tikhonov standard
form. Thus, there exists a local analytic diffeomorphism $\Phi$ and a vector
field $\widetilde{q}$ such that

$D\Phi(x)\widetilde{h}(x,\varepsilon)=\widetilde{q}(\Phi(x),\varepsilon)$

and consequently

$D\Phi(x)D\widetilde{h}(x,\varepsilon)=D\widetilde{q}(\Phi(x),\varepsilon)D\Phi(x),\quad
x\in\widetilde{Y}\cap K.$

Therefore the Jacobian of $\widetilde{h}$ at $x$ and the Jacobian of
$\widetilde{q}$ at $\Phi(x)$ are conjugate; in particular they have the same
characteristic polynomial. Denoting by $\widetilde{\nu}_{i}(y)$ the
coefficients of the characteristic polynomial of $D\widetilde{q}(y)$, this
means

$\widetilde{\sigma}_{i}(x,\varepsilon)=\widetilde{\nu}_{i}(\Phi(x),\varepsilon)\text{
for all }x\in\widetilde{Y}\cap K.$

Since $\widetilde{q}$ is in Tikhonov standard form, we have

$\widetilde{q}(y,\varepsilon)=\begin{pmatrix}\varepsilon\widetilde{q}_{1}(y,\varepsilon)\\\
\widetilde{q}_{2}(y,\varepsilon)\end{pmatrix},$

with $q_{1}$ having $s$ entries, and

$D\widetilde{q}(y,\varepsilon)=\begin{pmatrix}\varepsilon
D\widetilde{q}_{1}(y,\varepsilon)\\\
D\widetilde{q}_{2}(y,\varepsilon)\end{pmatrix}.$

Thus, every entry of the first $s$ rows of the Jacobian is a multiple of
$\varepsilon$, and with the Laplace expansion of the determinant this implies

$\widetilde{\nu}_{i}(x,\varepsilon)=\varepsilon^{i-n+s}\,\widehat{\nu}_{i}(x,\varepsilon),\quad
n-s<i\leq n,$

and finally (44). ∎

Now we turn to determining the orders of the eigenvalues.

######  Lemma 7.

With objects and notation as in Lemma 6, let (44) hold, and furthermore
consider the nondegeneracy conditions:

  1. (i)

$\widehat{\sigma}_{n-s}(x,0)\not=0$ and $\widehat{\sigma}_{n}(x,0)\not=0$ on
$\widetilde{Y}\cap K$.

  2. (ii)

The polynomials

$\widehat{\sigma}_{n-s}(x,0)\tau^{s}+\widehat{\sigma}_{n-s+1}(x,0)\tau^{s-1}+\cdots+\widehat{\sigma}_{n}(x,0)$
(94)

admit only simple zeros, for all $x\in\widetilde{Y}\cap K$.

  1. (a)

Whenever (i) holds then the zeros $\lambda_{i}(x,\varepsilon)$ of the
characteristic polynomial can be labeled such that

$\lambda_{1}(x,0)\not=0,\ldots,\lambda_{n-s}(x,0)\not=0\quad\text{ on
}\widetilde{Y}\cap K,$

and

$\lambda_{i}(x,\varepsilon)=\varepsilon\widehat{\lambda}_{i}(x,\varepsilon),\quad
x\in\widetilde{Y}\cap K,\quad i>n-s,$

with continuous $\widehat{\lambda}_{i}$ such that
$\widehat{\lambda}_{i}(x,0)\not=0$ on $\widetilde{Y}\cap K$, $n-s+1\leq i\leq
n$.242424The $\widehat{\lambda}_{i}$ can be represented as convergent power
series in $(x,\varepsilon^{1/m})$ for some positive integer $m$.

  2. (b)

Whenever (ii) holds in addition to (i) then all $\widehat{\lambda}_{i}$,
$n-s+1\leq i\leq n$, are analytic in $(x,\varepsilon)$.

###### Proof.

The proof rests on the Newton-Puiseux theorem and on Hensel’s lemma; we refer
specifically to Abhyankar [1, Lectures 12 and 13]. According to Newton-Puiseux
the equation $\lambda^{n}+\sum\widetilde{\sigma}_{i}\lambda^{n-i}=0$ admits
series solutions

$\lambda=\alpha\varepsilon^{\gamma}+\cdots$

in rational exponents of $\varepsilon$, with a positive rational number
$\gamma$ and $\alpha\not=0$. For such an expansion to hold with some $\gamma$
and $\alpha\not=0$, cancellation of lowest order terms in (10) is necessary.
The lowest orders of the terms in the monomials are

$(n-i)\gamma\text{ for }0\leq i\leq s,\quad\text{ and
}\quad(n-j)\gamma+j-n+s\text{ for }s+1\leq j\leq n,$

and for cancellation one must have equality between two of these orders.
Clearly two orders in the first block cannot be equal. Assuming that an order
from the first block equals an order in the second block, we get

$(n-i)\gamma=(n-j)\gamma+j-n+s\Rightarrow\gamma=\frac{j-(n-s)}{j-i}<1\text{
unless }i=n-s.$

But in case $\gamma<1$ the lowest order equals $s\gamma$, with no
cancellation; so only $\gamma=1$ remains. Finally, if two orders in the second
block are equal then one directly sees $\gamma=1$. This shows part (a).

Continuing the argument, $\gamma=1$ implies that precisely the monomials of
degree $\leq n-s$ contribute to the lowest order, and the ansatz yields

$\widehat{\sigma}_{n-s}(x,0)\alpha^{s}+\widehat{\sigma}_{n-s+1}(x,0)\alpha^{s-1}+\cdots+\widehat{\sigma}_{n}(x,0)=0,$

thus $s$ distinct choices for $\alpha$ by condition (ii), and $\alpha\not=0$.
By Hensel’s lemma, each choice for $\alpha$ yields a series
$\lambda=\alpha\varepsilon+\cdots$, in positive integer powers of
$\varepsilon$. This shows part (b). ∎

######  Remark 6.

In case $s=1$ the second condition is automatic. Therefore Lemma 2 (a) is also
proven.

## References

  * [1] S.S. Abhyankar: Algebraic geometry for scientists and engineers. American Mathematical Society, Providence, RI (1990).
  * [2] V.I. Arnold: Ordinary differential equations. Springer, Berlin (1992).
  * [3] N. Berglund, B. Gentz: Noise-induced phenomena in slow-fast dynamical systems. A sample-paths approach. Springer, London (2006).
  * [4] G. E. Briggs, J. B. S. Haldane: A note on the kinetics of enzyme action. Biochem. J., 19, 338 - 339 (1925).
  * [5] M.S. Calder, D. Siegel: Properties of the Michaelis-–Menten mechanism in phase space. J. Math. Anal. Appl. 339, 1044–1064 (2008).
  * [6] P.T. Cardin, M.A. Texeira: Fenichel theory for multiple time scale singular perturbation problems. SIAM J. Appl. Dyn. Sys. 16, 1452-1452 (2017)
  * [7] C. Chicone: Ordinary differential equations with applications.  Second edition. Springer, New York (2006).
  * [8] B. Choi, G. A. Rempala, J. K. Kim: Beyond the Michaelis–Menten equation: Accurate and efficient estimation of enzyme kinetic parameters. Sci. Rep. 7, 1-11 (2017).
  * [9] J. Eilertsen, M.R. Roussel, S. Schnell, S. Walcher: On the quasi-steady state approximation in an open Michaelis–Menten reaction mechanism. AIMS Math. 6, no. 7, 6781 – 6814 (2021).
  * [10] J. Eilertsen, S. Schnell: A kinetic analysis of coupled (or auxiliary) enzyme reactions, Bull. Math Biol. 80, 3154-3183 (2018).
  * [11] J. Eilertsen, S. Schnell: The quasi-steady-state approximations revisited: Timescales, small parameters, singularities, and normal forms in enzyme kinetics, Math. Biosci. 325, 108339 (2020).
  * [12] J. Eilertsen, S. Schnell, S. Walcher: On the anti-quasi-steady-state conditions of enzyme kinetics. Mathematical Biosciences 350, 108870 (2022).
  * [13] J. Eilertsen, W. Stroberg, S. Schnell: Phase-plane geometries in coupled enzyme assays, Math. Biosci. 306, 126-135 (2018).
  * [14] J. Eilertsen, M.A. Tyczynska, S. Schnell: Hunting $\epsilon$: The origin and validity of quasi-steady-state reductions in enzyme kinetics. SIAM J. Appl. Dyn. Syst. 20, 2450-2481 (2021).
  * [15] L.C. Evans: Partial Differential Equations, 2nd ed. American Mathematical Society, Graduate Studies in Mathematics, v. 19 (2014).
  * [16] P.A. Eyers, J.M. Murphy: The evolving world of pseudoenzymes: Proteins, prejudice and zombies. BMC Biol 14, 98 (2016).
  * [17] N. Fenichel: Geometric singular perturbation theory for ordinary differential equations. J. Differential Equations 31(1), 53–98 (1979).
  * [18] F.R. Gantmacher: Applications of the theory of matrices. Dover, Mineola (2005).
  * [19] A. Goeke, S. Walcher: Quasi-steady state: Searching for and utilizing small parameters, In: A. Johann, H.-P. Kruse, F. Rupp, S. Schmitz (eds.): Recent Trends in Dynamical Systems, Springer-Verlag, New York, pp. 153 - 178 (2013).
  * [20] A. Goeke, S. Walcher: A constructive approach to quasi-steady state reduction. J. Math. Chem. 52, 2596 - 2626 (2014).
  * [21] A. Goeke, S. Walcher, E. Zerz: Determining “small parameters” for quasi-steady state. J. Diff. Equations 259, 1149–1180 (2015).
  * [22] A. Goeke, S. Walcher, E. Zerz: Classical quasi-steady state reduction – A mathematical characterization. Physica D 345, 11 - 26 (2017).
  * [23] F.G. Heineken, H.M. Tsuchiya, R. Aris: On the mathematical status of the pseudo-steady hypothesis of biochemical kinetics. Math. Biosci. 1, 95 - 113 (1967).
  * [24] J. Keener, J. Sneyd: Mathematical physiology I: Cellular physiology, Second Ed. Springer-Verlag, New York (2009).
  * [25] N. Kruff, S. Walcher: Coordinate-independent singular perturbation reduction for systems with three time scales. Math. Biosci. Eng. 16, 5062 – 5091 (2019).
  * [26] M. Krupa, P. Szmolyan: Extending slow manifolds near transcritical and pitchfork singularities. Nonlinearity 14, 1473 – 1491 (2001).
  * [27] S.H. Lam, D.A. Goussis: The CSP method for simplifying kinetics. Int. J. Chem. Kinet. 26, 461 - 486 (1994)
  * [28] L. Michaelis, M.L. Menten: Die Kinetik der Invertinwirkung. Biochem. Z. 49, 333 - 369 (1913).
  * [29] W. G. Miller, R. A. Alberty: Kinetics of the reversible Michaelis–Menten mechanism and the applicability of the steady-state approximation. J. Am. Chem. Soc. 80 5146-5151 (1958).
  * [30] J.M. Murphy, H. Farhan, P.A. Eyers: Bio-Zombie: The rise of pseudoenzymes in biology. Biochem. Soc. Trans. 45 537-544 (2017).
  * [31] L. Noethen, S. Walcher: Quasi-steady state in the Michaelis–Menten system. Nonlinear Anal. Real World Appl. 8, 1512 - 1535 (2007).
  * [32] L. Noethen, S. Walcher: Tikhonov’s theorem and quasi-steady state. Discrete Contin. Dyn. Syst. Ser. B 16(3), 945–961 (2011).
  * [33] B.O. Palsson, E.N. Lightfoot: Mathematical modelling of dynamics and control in metabolic networks. I. On Michaelis–Menten kinetics. J. Theor. Biol. 111, 273 - 302 (1984).
  * [34] D.G. Patsatzis, D.A. Goussis: A new Michaelis–Menten equation valid everywhere multi-scale dynamics prevails. Math. Biosci. 315, 108220, 13 pp. (2019).
  * [35] J.G. Reich, E.E. Selkov: Mathematical analysis of metabolic networks. FEBS Letters 40, Suppl. 1, S119 - S127 (1974).
  * [36] S. Schnell, P.K. Maini: Enzyme kinetics at high enzyme concentration. Bull. Math. Biol. 62, 483–499 (2000).
  * [37] S. Schnell, C. Mendoza: A fast method to estimate kinetic constants for enzyme inhibitors. Acta Biotheor. 49, 109–113 (2001).
  * [38] S. Schnell: Validity of the Michaelis-‐Menten equation–steady‐state, or reactant stationary assumption: that is the question. FEBS Journal 281, 464–-472 (2014).
  * [39] L.A. Segel: On the validity of the steady state assumption of enzyme kinetics. Bull. Math. Biol. 50, 579–593 (1988).
  * [40] L.A. Segel, M. Slemrod: The quasi-steady-state assumption: A case study in perturbation. SIAM Review 31, 446 - 477 (1989).
  * [41] M. Seshadri, G. Fritzsch: Analytical solutions of a simple enzyme kinetic problem by a perturbative procedure. Biophys. Struct. Mech. 6, 111 - 123 (1980).
  * [42] W. Stroberg, S. Schnell: On the estimation errors of $K_{M}$ and $V$ from time-course experiments using the Michaelis–Menten equation. Biophys. Chem. 219, 17-27 (2016).
  * [43] A.R. Tzafriri: Michaelis–Menten kinetics at high enzyme concentrations. Bull. Math. Biol. 65, 1111–1129 (2003).
  * [44] A.N. Tikhonov: Systems of differential equations containing a small parameter multiplying the derivative (in Russian). Math. Sb. 31, 575–586 (1952).
  * [45] F. Verhulst: Methods and applications of singular perturbations. Boundary layers and multiple timescale dynamics, Springer, New York (2005).
  * [46] W. Walter: Ordinary differential equations. Graduate Texts in Mathematics 182. Springer, New York (1998).
  * [47] F.Y.M. Wan: Mathematical models and their analysis. SIAM, Philadelphia (2018).
  * [48] A. Zagaris, H.G. Kaper, T.J. Kaper: Analysis of the computational singular perturbation reduction method for chemical kinetics. J. Nonlinear Sci. 14, 59–91 (2004).