# The Quasi-Steady-State Approximations revisited: Timescales, small
parameters, singularities, and normal forms in enzyme kinetics

Justin Eilertsen  Santiago Schnell111Corresponding author<EMAIL_ADDRESS>Department of Molecular & Integrative Physiology, University of Michigan
Medical School, Ann Arbor, MI 48109, USA  Department of Computational Medicine
& Bioinformatics, University of Michigan Medical School, Ann Arbor, MI 48109,
USA

###### Abstract

In this work, we revisit the scaling analysis and commonly accepted conditions
for the validity of the standard, reverse and total quasi-steady-state
approximations through the lens of dimensional Tikhonov–Fenichel parameters
and their respective critical manifolds. By combining Tikhonov–Fenichel
parameters with scaling analysis and energy methods, we derive improved upper
bounds on the approximation error for the standard, reverse and total quasi-
steady-state approximations. Furthermore, previous analyses suggest that the
reverse quasi-steady-state approximation is only valid when initial enzyme
concentrations greatly exceed initial substrate concentrations. However, our
results indicate that this approximation can be valid when initial enzyme and
substrate concentrations are of equal magnitude. Using energy methods, we find
that the condition for the validity of the reverse quasi-steady-state
approximation is far less restrictive than was previously assumed, and we
derive a new “small” parameter that determines the validity of this
approximation. In doing so, we extend the established domain of validity for
the reverse quasi-steady-state approximation. Consequently, this opens up the
possibility of utilizing the reverse quasi-steady-state approximation to model
enzyme catalyzed reactions and estimate kinetic parameters in enzymatic assays
at much lower enzyme to substrate ratios than was previously thought.
Moreover, we show for the first time that the critical manifold of the reverse
quasi-steady-state approximation contains a singular point where normal
hyperbolicity is lost. Associated with this singularity is a transcritical
bifurcation, and the corresponding normal form of this bifurcation is
recovered through scaling analysis.

###### keywords:

Enzyme Kinetics , Singular Perturbation , Quasi-Steady-State , Singularity ,
Normal Hyperbolicity , Transcritical Bifurcation

##  1 Introduction

Perhaps the most well-known reaction in biochemistry is the Michaelis-Menten
(MM) reaction mechanism, (1), which describes the catalytic conversation of a
substrate, $S$, into a product, $P$. Catalysis is achieved through means of an
enzyme, $E$, that reversibly binds with the substrate, forming an intermediate
complex, $C$. In turn, $C$ irreversibly disassociates into enzyme and product
molecules:

$\ce{S+E<=>[k_{1}][k_{-1}]C->[k_{2}]E+P}.$ (1)

For theorists and applied mathematicians, whose aim is to mathematically
describe the dynamics of chemical networks and metabolic pathways [1, 2], the
MM reaction mechanism serves as a building block to these more complex
systems. In general, there are two ways to mathematically model enzyme
catalyzed reactions. At low concentrations of chemical species stochastic
models are generally favorable since they describe the seemingly random
collisions of reactant molecules in intracellular environments [3, 4]. In
contrast, when concentrations are high and the chemical species are well-
mixed, the MM mechanism can be appropriately modeled as a set of nonlinear
ordinary differential equations. Although one can argue that deterministic
models are in some sense more manageable than stochastic models, nonlinear
deterministic models rarely admit closed form solutions, and therefore model
reduction techniques must be employed in order to simplify the model equations
so that approximate solutions can obtained and analyzed. Typically, model
reduction is synonymous with approximating the model dynamics on an invariant
manifold, and the advent of powerful computer algebra systems aids in the
systematic reduction of high-dimensional or even infinite-dimensional (i.e.,
partial differential equations and delay differential equations) dynamical
systems. Not surprisingly, there is a large body of literature that clearly
illustrates how model reduction can uncover, quantify, and explain the various
nonlinear phenomena arising not only in chemistry, but also in biology [5, 6,
7] and physics [8, 9, 10, 11].

Historically, the most widely-utilized reduction technique in deterministic
enzyme kinetics has been the singular perturbation method [12], in which a
reduced model is constructed by approximating the flow of the model equations
on a slow invariant manifold (SIM). Since the dimension of the SIM is less
than the dimension of the phase-space, approximation of the dynamics on the
SIM permits a reduction in the dimension of the problem. The singular
perturbation method exploits the presence of disparate fast and slow
timescales; the work of Tikhonov [13] and Gradshtein [14] provides the
necessary rigorous foundation for the construction of a reduced model.
Briefly, when fast and slow timescales are present, the differential equations
that model the MM reaction with a timescale separation can be expressed in the
form

$\displaystyle\dot{x}$ $\displaystyle=f(x,y;\varepsilon),$ (2a)
$\displaystyle\varepsilon\dot{y}$ $\displaystyle=g(x,y;\varepsilon),$ (2b)

where $0<\varepsilon\ll 1$. Setting $\varepsilon=0$ yields a differential-
algebraic-equation

$\displaystyle\dot{x}$ $\displaystyle=f(x,y;0),$ (3a) $\displaystyle 0$
$\displaystyle=g(x,y;0),$ (3b)

and, according to Tikhonov’s theorem, equations (3a)–(3b) provide a very good
approximation to the dynamics of (2a)–(2b) when $\varepsilon$ is sufficiently
small. The reduced system (3a)–(3b) provides a simpler, and often times more
tractable, reduced mathematical model that is commonly referred to as a quasi-
steady-state approximation (QSSA).

The hope is that the condition that supports the validity of (3a)–(3b) (i.e.,
$\varepsilon\ll 1$) is easy to implement in enzymatic assays, so that precise
and accurate measurements of the kinetic parameters pertinent to the reaction
can be made by fitting the QSSA model to experimental time course data [15,
16]. Due to the necessity of disparate timescales, slow manifold reduction may
not be possible in every physical scenario. Therefore, the challenge for
theorists is not only to derive a reduced model that has suitable utility, but
also to determine the unique physical and chemical conditions that permit the
validity of the associated reduction. Thus, an important task of the
theoretician is determine the exact conditions for which the reduced model is
valid [17]. Mathematically, this translates to determining “$\varepsilon$,” a
(typically) dimensionless parameter that may not be unique. The non-uniqueness
of “$\varepsilon$” adds complication, since some “epsilons” are better than
others. The best-known example of the “non-uniqueness dilemma” resides the
history of the derivation of the Michaelis–Menten (MM) equation

$v=\frac{Vs}{K_{M}+s},$ (4)

which is obtained by applying the QSSA to the MM reaction mechanism (1). In
(4), $v$ is the velocity of product formation in the reaction, $V$ is the
limiting rate of the reaction, $K_{M}$ is the Michaelis constant, and $s$ is
the free substrate concentration for the reaction. Alternatively, the MM
equation is often referred to as the standard quasi-steady-state approximation
(sQSSA). In 1967, Heineken et al. [18] formally applied, for the first time,
the standard QSSA to the nonlinear differential equations governing the MM
reaction mechanism (1) via singular perturbation analysis. Based on the
findings of Laidler [19], a pioneer in chemical kinetics and authority on the
physical chemistry of enzymes, Heineken et al. [18] introduced a consistent
“$\varepsilon$” for the MM reaction mechanism through scaling analysis. A more
rigorous analysis of the sQSSA was introduced by Reich and Sel’kov [20] and
Schauer and Heinrich [21], from which other “epsilons” where determined by
proposing conditions that minimized the errors in the implementation of the
sQSSA. In 1988, Lee A. Segel [22] derived the widely accepted criterion for
the validity of the sQSSA and the derivation of the MM equation (4) by
estimating the slow and fast timescales of the reaction. Segel [22]
illustrated that prior physico-chemical knowledge about the reaction dynamic
is instrumental in deriving the fast and slow timescales and uncovering the
most general criteria for the validity of the sQSSA. As a direct result of
Segel’s scaling method, the conditions for the validity of the sQSSA were
derived for suicide substrates [23, 24], alternative substrates [25], fully
competitive enzyme reactions [26], zymogen activation [27], and coupled enzyme
assays [28, 29]. Segel’s scaling approach was also applied to the analysis of
the MM reaction mechanism (1) to extend the validity of the QSSA in different
regions of the initial enzyme and substrate concentration parameter space via
the reverse QSSA (rQSSA) [30, 31], and the total QSSA (tQSSA) [32, 33, 34].

Despite the power of Segel’s scaling and simplification analysis for the MM
reaction mechanism (1), there is still a fundamental challenge with its
implementation via the rQSSA: there has never been a small parameter (i.e., a
specific “$\varepsilon$”), analogous to the one obtained by Segel [22] for the
sQSSA, that is as effective at determining when the rQSSA is valid.
Unfortunately, estimating the fast timescale associated with the rQSSA is
difficult, and there have been fundamental disagreements in the reported
estimates [31]. Thus, timescale estimation continues to be the “Achilles’
heel” of the rQSSA. This raises the obvious question of whether or not
timescale estimation is truly the best approach towards resolving this
problem.

Recently, Walcher and his collaborators [35, 36, 37] demonstrated that
identifying a Tikhonov–Fenichel parameter (TFP) is an effective way to
determine a priori conditions that suggest the validity of the QSSA.
Essentially, a TFP is a dimensional parameter – such as a rate constant or
initial concentration of a species – that, when identically zero, ensures the
existence of a manifold of equilibrium points. Such manifolds are central to
geometric singular perturbation theory (GSPT) and, as a result of Fenichel
[38], it is well-understood that when certain conditions hold, the existence
of a critical manifold of equilibria ensures the existence of a SIM once the
TFP is small but non-zero. In this sense the origin of the SIM can be linked
to the vanishing of a specific dimensional parameter, and a sufficiently small
TFP ensures the existence of a SIM and the corresponding validity of a QSSA.
However, the identification of a TFP does not diminish the importance of the
asymptotic small parameter $\varepsilon$, since it is essential to define what
physically constitutes “small” when a parameter is non-zero. In other words,
we must still answer the question: how small should Tikhonov–Fenichel
parameters be in comparison to other dimensional parameters in order to yield
an accurate reduced model?

In this paper, our primary objective is to determine a specific small
parameter that determines the validity of the rQSSA to the MM reaction
mechanism (1), but also to convey subject matter that can be quite technical
in a language capable over reaching a wider audience that is not limited to
applied mathematicians and physical chemists. In Section 2, we review the
conditions for the validity of the various quasi-steady-state approximations
that are commonly employed to approximate the long-time dynamics of the MM
reaction mechanism (1): namely, the sQSSA, the rQSSA, and tQSSA. In Section 3,
we assess the validity of specific QSS reductions by employing geometric
singular perturbation theory, and illustrate how assumptions about the
validity of the sQSSA based on Segel’s timescale separation can lead to
erroneous conclusions. In Section 4, we introduce two methods that do not rely
on scaling or timescale separation: Tikhonov–Fenichel parameters and energy
methods, both of which can be employed to determine the validity of the QSSA.
In Section 5, we analyze the validity of the rQSSA using the methods discussed
in Section 4, and in Section 6, we discuss timescale separation and the
hierarchy of small parameters that support the justification of the sQSSA,
rQSSA, and the tQSSA. Finally, in our discussion (Section 7), we summarize our
results and critique some of the conclusions drawn about the validity of the
sQSSA and the rQSSA in the previous analyses of Segel and Slemrod [30], and
Schnell and Maini [31], respectively. We also discuss the impact our results
will have on experimental assays, and how the methods we utilize can be
employed to analyze more complicated reactions and experimental assays.

##  2 Applying the Quasi-Steady-State Approximations to the Michaelis–Menten
reaction mechanism: Scaling and simplification approaches

In this section we review the application of the different versions of the
QSSA for the MM reaction mechanism (1), and the mathematical justification for
the application of each approximation. We also present the timescales and
criterion for the validity of the QSSA originally derived by Segel [22] and
Segel and Slemrod [30].

###  2.1 Heuristic estimation of fast and slow timescales using the standard
Quasi-Steady-State Approximation

The mathematical model that describes (deterministically) the reaction
mechanism (1) is a set of nonlinear ordinary differential equations,

$\displaystyle\dot{s}$ $\displaystyle=-k_{1}(e_{0}-c)s+k_{-1}c,$ (5a)
$\displaystyle\dot{c}$ $\displaystyle=k_{1}(e_{0}-c)s-(k_{-1}+k_{2})c,$ (5b)
$\displaystyle\dot{p}$ $\displaystyle=k_{2}c.$ (5c)

The lowercase $s$, $c$ and $p$ denote concentrations of $S$, $C$ and $P$,
respectively, and “$\dot{\phantom{x}}$” denotes differentiation with respect
to time. The equation that models the time evolution of the enzyme
concentration, $e$, has been eliminated via the enzyme conservation law,
$e+c=e_{0}$, and note that adding together (5a)–(5c) yields the substrate
conservation law, $s+c+p=s_{0}$. From the model equations (5a)–(5c), we see
that the MM reaction mechanism (1) is parametrically controlled by the initial
substrate concentration, $s_{0}$, the initial enzyme concentration, $e_{0}$,
and the magnitudes of the rate constants $k_{1}$, $k_{-1}$ and $k_{2}$.

Since the model equations (5a)–(5c) are nonlinear, closed form solutions are
intractable. However, it is well-established that whenever the quantity,
$\varepsilon_{SS}$, is very small,

$\varepsilon_{SS}=\cfrac{e_{0}}{K_{M}+s_{0}}\ll 1,$ (6)

equations (5a)–(5c) can be approximated by the system of differential-
algebraic equations:

$\dot{s}\sim-\cfrac{Vs}{K_{M}+s}\sim-\dot{p},\quad
c\sim\cfrac{e_{0}s}{K_{M}+s}.$ (7)

Here the Michaelis constant is defined as $K_{M}\equiv(k_{-1}+k_{2})/k_{1}$,
and the limiting rate is $V\equiv k_{2}e_{0}$. Formally, the equations (7) are
collectively referred to as the sQSSA, and condition (6) is known as the
reactant-stationary-assumption (RSA) [39]. When the RSA holds, and the sQSSA
(7) is valid, the intermediate complex reaches its maximum value very quickly.
By comparison, the conversion of $S$ to $P$ is slow, and the MM reaction
mechanism (1) is characterized (temporally) by the two disparate timescales,
$t_{C}$ and $t_{\mathcal{D}}$, that respectively quantify the approximate
amount of time it takes $c$ to become maximal, and for $s$ to significantly
deplete. Thus, when the RSA holds, the reaction consists of a fast timescale,
$t_{C}$, and a slow timescale, $t_{\mathcal{D}}$. Both timescales depend on
the rate constants, as well as the initial substrate and enzyme
concentrations:

$\displaystyle t_{C}$ $\displaystyle=\cfrac{1}{k_{1}(s_{0}+K_{M})},$ (8a)
$\displaystyle t_{\mathcal{D}}$ $\displaystyle=\cfrac{K_{M}+s_{0}}{V}.$ (8b)

Segel [22] first derived $t_{C}$ and $t_{\mathcal{D}}$ using heuristic
methods. To estimate $t_{C}$, Segel noted, no doubt from the earlier work of
Heineken, Tsuchiya, and Aris [18], that if $e_{0}\ll s_{0}$, then $c$ should
reach its maximum value very quickly, and there should be correspondingly very
little substrate depletion during the initial accumulation of $c$.
Consequently, one can assume $s\approx s_{0}$ during the initial accumulation
of $c$. To obtain an equation that models the initial increase in $c$, we set
$s=s_{0}$ in (5b), which gives rise to a linear ordinary differential equation
for $c$:

$\dot{c}\sim k_{1}(e_{0}-c)s_{0}-(k_{-1}+k_{2})c.$ (9)

Since (9) is linear, its solution is easily attainable, and is given by:

$c(t)\sim\varepsilon_{SS}s_{0}(1-\displaystyle e^{\displaystyle-t/t_{C}}).$
(10)

Thus, $t_{C}$ is the characteristic timescale of the initial fast transient,
and it is valid as long as $s$ is approximately constant during the initial
accumulation of $c$. The reader should be aware of the fact that (9) is only
valid for timescales on the order of $t_{C}$. Hence, the solution (10) is only
an initial, temporary approximation to the kinetics of (1) when the RSA is
valid. As such, it is often referred to as an initial layer.

Since $t_{C}$ was derived under the presupposition of a negligible depletion
of $s$ during the initial accumulation if $c$, how then, can we determine a
dimensionless parameter that must be small in order to ensure the depletion of
$s$ is in fact negligible over $t_{C}$? Segel [22] reasoned that since the
maximum rate of depletion of $s$ is identically $-k_{1}e_{0}s_{0}$, it
suffices to demand that

$\\{\max|\dot{s}|\cdot t_{C}\ll s_{0}\\}\equiv\\{\varepsilon_{SS}\ll 1\\}$
(11)

hold in order to justify the use of (9). Thus, the RSA (6) ensures that the
approximation (9) is valid.

After $c$ reaches its maximum value, the QSS phase of the reaction ensues, and
$c$ starts to slowly deplete. During the depletion phase $\dot{c}$ is very
small but non-zero, and thus we say that $c$ evolves in a QSS or is “slaved”
by $s$ since, according to the system (7), the approximate concentration of
$c$ is parametrically determined by the concentration of $s$. To estimate the
depletion timescale, Segel [22] assumed that the sQSSA (7) was valid once
$t_{C}\lesssim t$, and that the rate of depletion of substrate was
approximately

$\dot{s}\sim-\cfrac{Vs}{K_{M}+s}\leq 0.$ (12)

Since $s\approx s_{0}$ for $t\lesssim t_{C}$, it stands to reason that the
sQSSA (12) can be supplied with the boundary condition $s(t_{C})=s_{0}$, and
thus the maximum rate of depletion in the QSS phase of the reaction is
approximately

$\max|\dot{s}|\approx\cfrac{Vs_{0}}{K_{M}+s_{0}}.$ (13)

Finally, to calculate the depletion timescale, Segel [22] divided the total
change in substrate over the course of the reaction, $\Delta s=s_{0}$, by the
maximum rate of depletion in the QSS phase (13):

$t_{\mathcal{D}}=\cfrac{|\Delta
s|}{\max|\dot{s}|}=\cfrac{s_{0}}{\cfrac{Vs_{0}}{K_{M}+s_{0}}}=\cfrac{K_{M}+s_{0}}{V}.$
(14)

As noted in the introduction, the effective use of the singular perturbation
method relies on separation of fast and slow timescales, and the question that
naturally arises from the heuristic line of reasoning arises is: Does the RSA
(6) ensure separation of timescales? That is, does it hold that $t_{C}\ll
t_{\mathcal{D}}$ whenever $\varepsilon_{SS}\ll 1$ ? The answer is yes: the RSA
not only ensures that there is negligible formation of product (or depletion
of substrate) for $t\lesssim t_{C}$, it also ensures that the timescales,
$t_{C}$ and $t_{\mathcal{D}}$, are widely separated:

$\cfrac{t_{C}}{t_{\mathcal{D}}}\equiv\epsilon=\cfrac{\varepsilon_{SS}}{\bigg{(}1+\cfrac{k_{-1}}{k_{2}}\bigg{)}\bigg{(}1+\cfrac{s_{0}}{K_{M}}\bigg{)}}\leq\varepsilon_{SS}.$
(15)

Thus, the RSA (i.e., the condition that $\varepsilon_{SS}\ll 1$) induces the
separation of fast and slow timescales. Consequently, the RSA is understood to
be a sufficient condition for the validity of the sQSSA (7), and the
separation of fast and slow timescales (i.e., $\epsilon\ll 1$) is a necessary
(but not sufficient) condition for the sQSSA.

When the RSA is valid, the trajectories in the $(s,c)$-plane have a very
recognizable form: during fast transient phase, the trajectories move almost
vertically towards the $c$-nullcline, given by the curve
$\\{(s,c)\in\mathbb{R}^{2}:c-e_{0}s/(K_{M}+s)=0\\}$, on which $\dot{c}=0$. The
intermediate complex $c$ reaches its maximum value once the trajectory reaches
the $c$-nullcline, at which point the QSS phase begins and the sQSSA is valid.
After $c$ reaches its maximum value, the trajectory closely follows the
$c$-nullcline towards the equilibrium point $(s,c)=(0,0)$ (see figure 1).

The fact that the phase-plane trajectory follows the $c$-nullcline on the slow
timescale, and approaches the $c$-nullcline in vertical fashion over the fast
timescale is the hallmark of the sQSSA. While the trajectory in figure 1
appears as though it is on the $c$-nullcline, it is important to note that it
is not. In fact, the trajectory approaches an invariant manifold,
$\mathcal{M}^{\varepsilon}$, that lies just above the $c$-nullcline. However,
since $\mathcal{M}^{\varepsilon}$ and the $c$-nullcline are very close
together when (6) holds, the $c$-nullcline is used to approximate
$\mathcal{M}^{\varepsilon}$ (again, see figure 1 for a detailed explanation
and illustration). The accuracy of using the $c$-nullcline as an approximation
is assessed asymptotically through scaling and non-dimensionalization methods,
and we review these methods in the subsection that follows.

Figure 1: The RSA ensures the validity of the sQSSA: Fast/slow dynamics occur
when $\varepsilon_{SS}\ll 1$ for the MM reaction mechanism (1). This figure is
an illustration of the phase-plane dynamics when $\varepsilon_{SS}\ll 1$ and
the sQSSA is valid; it serves to convey conceptual aspects of the phase-plane
when the RSA is valid. The red dashed curve is the $c$-nullcline. The black
curve the lies just above the $c$-nullcline is the invariant slow manifold.
The thin lines with double arrows are illustrative of typical phase-plane
trajectories when $\varepsilon_{SS}\ll 1$: trajectories starting on the
$s$-axis approach the $c$-nullcline almost vertically. Once the trajectory
reaches the $c$-nullcline it closely follows $\mathcal{M}^{\varepsilon}$
during the QSS phase of the reaction. Double arrows indicate fast dynamics,
and single arrows represent slow dynamics. The reader should bear in mind
that, slightly contrary to what is illustrated, the tangent vector to the
trajectory is perfectly horizontal when it intercepts the $c$-nullcline, since
its vertical velocity component vanishes at the point of interception.

###  2.2 Asymptotic justification for the sQSSA: Scaling and non-
dimensionalization

The identification of the small parameter (6) as an appropriate condition for
the validity of the standard QSSA is justified through scaling and non-
dimensionalization. Introducing the dimensionless variables $\tau\equiv
t/t_{C},\bar{s}\equiv s/s_{0}$ and $\bar{c}\equiv c/\varepsilon_{SS}s_{0}$
yields

$\displaystyle\cfrac{\text{d}\bar{s}}{\text{d}\tau}$
$\displaystyle=\varepsilon_{SS}[\mu\bar{c}\bar{s}-\bar{s}+\alpha\beta\bar{c}],$
(16a) $\displaystyle\cfrac{\text{d}\bar{c}}{\text{d}\tau}$
$\displaystyle=\bar{s}-\mu\bar{c}\bar{s}-\beta\bar{c},$ (16b)

where $\beta\equiv 1/(1+\sigma)$, $\sigma\equiv s_{0}/K_{M}$, $\kappa\equiv
k_{-1}/k_{2}$, $\mu\equiv 1-\beta$ and $\alpha\equiv\kappa/(1+\kappa)$. Thus,
we see from (16a) and (16b) that if $\varepsilon_{SS}\ll 1$, then $s\approx
s_{0}$ during the transient phase. Over the slow timescale,
$T\equiv\varepsilon_{SS}\tau$, we obtain

$\displaystyle\cfrac{\text{d}\bar{s}}{\text{d}T}$
$\displaystyle=\mu\bar{c}\bar{s}-\bar{s}+\alpha\beta\bar{c},$ (17a)
$\displaystyle\varepsilon_{SS}\cfrac{\text{d}\bar{c}}{\text{d}T}$
$\displaystyle=\bar{s}-\mu\bar{c}\bar{s}-\beta\bar{c}.$ (17b)

Formally, the standard QSSA is an asymptotic approximation of the dynamics on
the slow timescale obtained by setting $\varepsilon_{SS}=0$ in (17b) and
solving for $c$ in terms of $s$. Thus, the standard QSSA is the zeroth-order
approximation to (17a)–(17b) and, when $\varepsilon_{SS}\ll 1$, we refer to
$s$ as a slow variable, and $c$ as a fast variable.

In summary, we have two takeaways from this subsection:

  1. (i)

The RSA is synonymous with $\varepsilon_{SS}\ll 1$, and guarantees that there
is a negligible loss of substrate during the transient (fast) phase of the
reaction. When the RSA is valid, the QSSA system (7) is an appropriate
asymptotic approximation to (5a)–(5c).

  2. (ii)

Timescale separation is synonymous with $\epsilon\equiv
t_{C}/t_{\mathcal{D}}\ll 1$.

While (i) certainly implies (ii), the converse does not hold, which begs the
question: What happens when $\epsilon\ll 1$ but $\varepsilon_{SS}\sim 1$? We
discuss this special case in the subsection that follows.

###  2.3 The extended Quasi-Steady-State Approximation

The RSA and timescale separation provide two different small parameters:
$\varepsilon_{SS}$ and $\epsilon\equiv t_{C}/t_{\mathcal{D}}$. While
$\varepsilon_{SS}\ll 1$ ensures that $\epsilon\ll 1$, the converse is not
true. The obvious question is: Do QSS dynamics still prevail when
$\varepsilon_{SS}\sim 1$ but $t_{C}/t_{\mathcal{D}}\ll 1$? Segel and Slemrod
[30] proposed that the sQSSA (7) should still be valid at some point in the
time course of the reaction as long as $t_{C}\ll t_{\mathcal{D}}$, even if
$\varepsilon_{SS}\sim 1$. Rescaling the mass action equations (5a)–(5b) with
respect to $\bar{T}=t/t_{\mathcal{D}}$ yields

$\displaystyle\nu\beta\cfrac{\text{d}\bar{s}}{\text{d}\bar{T}}$
$\displaystyle=\mu\bar{c}\bar{s}-\bar{s}+\alpha\beta\bar{c},$ (18a)
$\displaystyle\nu\beta\varepsilon_{SS}\cfrac{\text{d}\bar{c}}{\text{d}\bar{T}}$
$\displaystyle=\bar{s}-\mu\bar{c}\bar{s}-\beta\bar{c},$ (18b)

where again, $\nu=k_{2}/(k_{-1}+k_{2})$, and $\beta=K_{M}/(K_{M}+s_{0})$. If
$\varepsilon_{SS}=1$, then setting $\nu=0$ (which implies $\alpha=1$) in
(18a)–(18b) yields a (dimensional) manifold:

$c=\cfrac{e_{0}s}{K_{S}+s},$ (19)

where $K_{S}=k_{-1}/k_{1}$ is the enzyme-substrate dissociation constant.
Interestingly, the manifold (19) is identical to the $s$-nullcline. However,
since the $c$-nullcline and the $s$-nullcline are indistinguishable in the
limiting case that corresponds to $k_{2}=0$, it stand to reason that
trajectories should closely follow both the $s$\- and $c$-nullcline when
$k_{2}$ is incredibly small, but non-zero. Numerical simulations confirm that
phase-plane trajectories closely follow the $c$-nullcline after a brief
transient (see, Figure 2).

Figure 2: The extended QSSA: Fast/slow dynamics still occur when
$\varepsilon_{SS}\sim 1$, but $\nu\ll 1$ for the MM reaction mechanism (1).
The solid black curves are the numerical solutions to the mass action
equations (5a)–(5c), and the dashed/dotted red curve is the ${c}$-nullcline.
In this simulation, $k_{1}=k_{-1}=10$, $k_{2}=0.01$ with $\kappa\equiv
k_{-1}/k_{2}=100$, $\nu=1/101$, and $\varepsilon_{SS}=1$ in arbitrary units
for illustrative purposes. Note that the trajectory is not a vertical straight
line during the initial phase of the reaction. As a consequence, the RSA does
not hold [30, 39], and there is noticeable depletion of substrate during the
transient phase. However, the trajectory still follows the $c$-nullcline after
a brief transient and achieves the slow QSS dynamics. The double arrows
indicate fast dynamics, and single arrows indicate slow dynamics.

The $c$-nullcline “following” that occurs when $\nu\ll 1$ and
$\varepsilon_{SS},\beta\sim 1$ is defined by Segel and Slemrod [30] as the
extended QSSA. One obvious difference between the sQSSA and the extended QSSA
is that there is a noticeable amount of substrate depletion when $\nu\ll 1$
but $\varepsilon_{SS}\sim 1$. Segel and Slemrod [30] proposed that the QSSA
system (7) should still hold after the transient, but that the primary
difference between the sQSSA and the extended QSSA was that one could not take
$s=s_{0}$ at the onset of the QSS phase in the case of the extended QSSA.
Thus, Segel and Slemrod [30] proposed that, when the extended QSSA is valid
but the RSA is invalid, the approximation

$\dot{s}\sim-\cfrac{Vs}{K_{M}+s},$ (20)

must be supplied with a boundary condition “$s_{0}^{e}$” that is less than
$s_{0}$. They employed a graphical method to estimate $s_{0}^{e}$ in their
original manuscript [30], and found that when $\beta=\varepsilon=1$,
$s_{0}^{e}\approx(\sqrt{2}-1)s_{0}$. Thus, they concluded that the QSS phase
of the reaction could be approximated with the initial value problem:

$\dot{s}\sim-\cfrac{Vs}{K_{M}+s},\quad s(0)=(\sqrt{2}-1)s_{0}.$ (21)

Moreover, they found that $s_{0}^{e}\to 0$ as $\beta$ increases, and
determined that the parameter regime where $\nu\beta\ll 1$ but
$\varepsilon_{SS}\sim 1$ results in a rather “uninteresting” domain of
applicability for extended QSSA222 Segel and Slemrod [30] did not take into
account product formation, and assumed that the complete depletion of
substrate was synonymous with the completion of the reaction. Unfortunately,
this assumption is false. For more details, we invite the reader to consult
[29].. Hence, the extended QSSA is formally defined to be the case when
$\nu\ll 1$, but $\varepsilon_{SS},\beta\sim 1$. Up to this point, we have two
dynamical regimes of interest:

  1. (i)

The RSA, in which $\varepsilon_{SS}\ll 1$ and (7) holds.

  2. (ii)

The extended QSSA, in which $\nu\ll 1$, but $\varepsilon_{SS},\beta\sim 1$ and
possibly (21) holds.

As a final remark of this subsection, we note that while the extended QSSA
seems reasonable, there are a few technical deficiencies that engulf the
validity of the approximation. First, setting $\nu=0$ results in the
coalescence of the $s$\- and $c$-nullclines, which implies that $s$ and $c$
are in some sense both in a QSS after the fast transient. Second, there is no
scaling justification for (21). Consequently, other than the observation that
the phase–plane trajectory follows the $c$-nullcline after the transient
phase, there is really no rigorous justification for (21). We will return to
these observations and discuss them in more detail in a later section.

###  2.4 The reverse Quasi-Steady-State Approximation

In contrast to the intermediate complex evolving in a QSS after the fast
transient, the rQSSA occurs when the substrate evolves in QSS after the
initial fast transient. Nguyen and Fraser [40] observed that the rQSSA
corresponds to trajectories closely following the $s$-nullcline instead of the
$c$-nullcline in the $s$–$c$ phase-plane (see Figure 3).

Figure 3: In the rQSSA, the substrate concentration is in a QSS phase in the
MM reaction mechanism (1). This figure is an illustration of the rQSSA. The
thick black curve is an illustration of the slow manifold; the dashed red
curve in the $c$-nullcline, and the dotted blue curve is the $s$-nullcline.
The thin black curve illustrates a typical trajectory when the rQSSA is valid.
During the transient phase, the complex concentration rapidly reaches it
maximum value, and the initial substrate concentration nearly vanishes. After
the transient phase, the trajectory follows the slow manifold,
$\mathcal{M}^{\varepsilon}$, during which time $s$ remains in a QSS phase.
Double arrows represent fast dynamics, and single arrows represent slow
dynamics.

To formulate an approximation to (5a)–(5c), it is assumed that virtually all
of the initial substrate will deplete during the initial transient phase, in
which case

$\cfrac{\text{d}c}{\text{d}t}+\cfrac{\text{d}p}{\text{d}t}=0$ (22)

approximately holds during the QSS phase. It also assumed that there is
negligible formation of product during the transient phase. Thus, in the
rQSSA, the underlying assumption is that there is a preliminary transient
phase in which nearly all of the substrate is depleted, but a negligible
amount of product is generated. As a consequence of this assumption, (22)
admits an approximate conservation law, $s_{0}\sim c+p$, that implies

$\cfrac{\text{d}p}{\text{d}t}=k_{2}(s_{0}-p)$ (23)

is the leading order approximation to $\dot{p}$ in the QSS phase of the
reaction. The above expression governs the production formation under the
rQSSA or rapid equilibrium approximation [41, 42]. Equation (23) is linear and
admits a closed-form solution

$p(t)=s_{0}(1-\displaystyle e^{\displaystyle-k_{2}t}),$ (24)

as well as a corresponding slow timescale: $\tilde{T}=k_{2}t$.

Is there a small parameter homologous to $\varepsilon_{SS}$ for the rQSSA that
justifies the validity of (23)? Typically, less-than-rigorous heuristic
arguments are employed to determine conditions for the validity of (23). By
utilizing timescale and geometrical arguments, Segel and Slemrod [30]
indirectly proposed333 They directly proposed that $K_{S}/e_{0}\ll 1$, but
noted that $k_{2}\approx k_{-1}$, which implies $K_{M}/e_{0}\ll 1$.

$\varepsilon^{*}\equiv\cfrac{K_{M}}{e_{0}}=\eta^{-1}\ll 1$ (25)

as a necessary condition for the validity of the rQSSA. However, there are
several problems with labeling (25) as a sufficient condition. First, as
Schnell and Maini [31] point out, $\varepsilon^{*}\ll 1$ does not ensure that
$\text{d}s/\text{d}t$ vanishes after the initial transient, since $K_{M}\ll
e_{0}\ll s_{0}$ would seem to imply the validity of the QSSA, and Segel and
Slemrod [30] only considered cases where $e_{0}\gg s_{0}$ in their analysis of
the rQSSA. Schnell and Maini [31] suggested that there should be both a
negligible loss of enzyme concentration and a maximal depletion of free
substrate over the fast transient; they proposed

$\varepsilon_{SM}\equiv\eta^{-1}+\cfrac{s_{0}}{e_{0}}=\cfrac{1}{\varepsilon_{SS}}\ll
1$ (26)

as an appropriate qualifier for the validity of the rQSSA. The condition (26)
suggests that since $e$ is in excess with respect to $s$, the majority of the
substrate should vanish over the fast transient. Second, from a more technical
perspective, Goeke et al. [43] remark that $K_{S}/e_{0}\ll 1$ is not
sufficient for the validity of the rQSSA; we will say more about Goeke’s
observation in the section that follows. For now we note that, heuristically,
if $1/k_{2}$ is a slow timescale, then it stands to reason that $k_{2}$ should
also be small in some sense.

In short, most of the conditions in the literature basically imply that
$\varepsilon_{SS}\gg 1$ in order for the rQSSA to be valid. In this work, the
question we want to address through our analysis is: How small should $k_{-1}$
and $k_{2}$ be in comparison to other variables to ensure the validity of the
rQSSA (i.e., what is the equivalent of $\varepsilon_{SS}$ in this case)?
Clearly, it must hold that $\varepsilon^{*}\ll 1$, but how does the initial
concentration of substrate, $s_{0}$, impact the accuracy of the approximation
(23)? This is still an open problem and, using non-scaling methods, we will
close it in a later section.

###  2.5 The total Quasi-Steady-State Approximation

The tQSSA provides an approximation that is valid over a broader parameter
domain. Utilizing the conservation law, $s_{0}=s+c+p$, and replacing
$s=s_{0}-(c+p)$ in (5b) yields

$\displaystyle\dot{c}$ $\displaystyle=k_{1}(e_{0}-c)(s_{0}-c-p)-k_{1}K_{M}c,$
(27a) $\displaystyle\dot{p}$ $\displaystyle=k_{2}c.$ (27b)

Setting the left hand side of (27a) to zero and solving for $c$ gives

$c=\cfrac{1}{2}(e_{0}+K_{M}+s_{0}-p)\pm\cfrac{1}{2}\sqrt{(e_{0}+K_{M}+s_{0}-p)^{2}-4e_{0}(s_{0}-p)}.$
(28)

Of the two possible roots given in (28), only the “$-$” root is physical.
Taking

$c\approx\cfrac{1}{2}(e_{0}+K_{M}+s_{0}-p)-\cfrac{1}{2}\sqrt{(e_{0}+K_{M}+s_{0}-p)^{2}-4e_{0}(s_{0}-p)}\equiv
c_{0}(p),$ (29)

and inserting it into (29) yields the tQSSA:

$\dot{p}\approx k_{2}c_{0}(p).$ (30)

Thus, the complex is taken to be the fast variable that evolves in a QSS after
a brief transient, and $p$, rather than $s$, is the corresponding slow
variable. Due to the fact that $p$ is the designated slow variable, the tQSSA
is valid over much larger domain than either the sQSSA or the rQSSA. As a
result of the more recent work of Borghans et al. [32], the term $s_{0}-p$ is
usually replaced by $s_{T}=s+c$ in the literature (see, for example, [33, 44,
45, 46, 47, 48, 16]), where $s_{T}$ is a lumped variable that denotes the
“total” substrate. Consequently, the equations

$\displaystyle\dot{c}$ $\displaystyle=k_{1}(e_{0}-c)(s_{T}-c)-k_{1}K_{M}c,$
(31a) $\displaystyle\dot{s}_{T}$ $\displaystyle=k_{2}c,$ (31b)

are more commonly employed than (27a)–(27b). Either way, as Pedersen [49]
points out, the coordinate systems $(s_{T},c)$ and $(p,c)$ are merely “two
sides of the same coin”. We will employ the $(p,c)$ coordinate system, and
avoid the lumped variable “$s_{T}$” in this paper as it is not a uniquely
distinguishable physical variable in experimental assays.

The validity of the tQSSA was first analyzed in the $(p,c)$ coordinate system
by Laidler [19], and later by Borghans et al. [32]. However, the accepted
condition for the validity of the tQSSA was derived by Tzafriri [34]. Briefly,
Tzafriri noted that if there is zero formation of product during the initial
transient, then the formation of $c$ during the transient phase can be
approximated by the Riccati equation

$\dot{c}=k_{1}(e_{0}-c)(s_{0}-c)-k_{1}K_{M}c.$ (32)

The solution to (32) admits a natural fast timescale:

$t_{C}^{*}=\cfrac{1}{k_{1}\sqrt{(e_{0}+K_{M}+s_{0})^{2}-4e_{0}s_{0}}}.$ (33)

The corresponding slow timescale, $t_{P}$, is calculated from the method of
Segel [22]

$t_{P}=\cfrac{|\Delta
p|}{\max|k_{2}c_{0}(p)|}=\cfrac{2s_{0}}{(e_{0}+K_{M}+s_{0})-\sqrt{(e_{0}+K_{M}+s_{0})^{2}-4e_{0}s_{0}}},$
(34)

and the separation of $t_{C}^{*}$ and $t_{P}$ is the accepted condition for
validity of the tQSSA:

$\varepsilon_{T}\equiv\cfrac{t_{C}^{*}}{t_{P}}=\cfrac{(e_{0}+K_{M}+s_{0})-\sqrt{(e_{0}+K_{M}+s_{0})^{2}-4e_{0}s_{0}}}{2s_{0}k_{1}\sqrt{(e_{0}+K_{M}+s_{0})^{2}-4e_{0}s_{0}}}\ll
1.$ (35)

Several scaling analyses have been carried out that justify the asymptotic
validity of the tQSSA, and we refer the reader to Schnell and Maini [33] and
Dell‘Acqua and Bersani [47] for detailed perturbation analyses, and Bersani et
al. [50] for a thorough review of the topic. Since the tQSSA is valid whenever
there is negligible formation of product during the fast transient, the
validity of the sQSSA, as well as the rQSSA and extended QSSA, implies the
validity of the tQSSA. In this sense, the tQSSA contains both the standard
QSSA and the rQSSA. Thus, $\varepsilon_{SS},\nu\ll 1\implies\varepsilon_{T}\ll
1$, but $1\ll\varepsilon_{SS}$ also implies $\varepsilon_{T}\ll 1$ (again, see
[32]).

It is worth noting that the fast and slow timescales derived by Tzafriri are
valid whenever there is negligible product formation during the fast
transient. Thus, they are valid in both the sQSSA and rQSSA regimes.
Additionally, the upper bound on $c$ is easily obtained:

$\sup
c\equiv\lambda=c_{0}(0)=\cfrac{1}{2}(e_{0}+K_{M}+s_{0})-\cfrac{1}{2}\sqrt{(e_{0}+K_{M}+s_{0})^{2}-4e_{0}s_{0}}.$
(36)

Consequently, we obtain universal timescales as well an upper bound on $c$
from the tQSSA formulation. At first glance the timescales $t_{C}^{*}$ and
$t_{P}$, as well as the upper bound $\lambda$, appear quite complicated, and
their utility in the context of scaling analysis appears limited. However,
when we introduce non-scaling methods, we will show that in fact these
“ingredients” are quite useful.

The tQSSA wraps up our introduction to each “QSSA” (i.e., the sQSSA, rQSSA,
and tQSSA) that is employed as a reduced model in enzyme kinetics. In summary,
we have:

  1. (i)

The RSA is valid when $\varepsilon_{SS}\ll 1$, and the corresponding QSSA is
given by (7).

  2. (ii)

The extended QSSA is valid when when $\nu\ll 1$, but ($\varepsilon_{SS}\sim 1$
and $\beta\sim 1$). The corresponding QSSA is possibly given by equation (21).

  3. (iii)

The rQSSA is at least valid when $1\ll\varepsilon_{SS}$, and the corresponding
QSSA is given by (23).

  4. (iv)

The tQSSA is valid whenever $\varepsilon_{T}\ll 1$, and the corresponding QSSA
is given by (30).

##  3 Applying the Quasi-Steady State Approximations with geometric singular
perturbation theory

The work of Fenichel [38] consists of a group of theorems that warrant the
existence of a slow invariant manifold, $\mathcal{M}$, in fast/slow systems of
the form of (2a)–(2b). In Section 3.1, we give a basic introduction to the
results obtained by Fenichel and introduce relevant terminology for our paper.
For a more detailed introduction, we refer the reader to [51, 52] and [53]. In
Section 3, we apply Geometric Singular Perturbation Theory (GSPT) to (5a)–(5c)
under the assumption that the RSA is valid. In Section 3.3, we apply GSPT to
(5a)–(5c) in the parameter-space region where the extended QSSA is valid, and
we show how GSPT can be used to prove that Segel and Slemrod’s conclusion that
the approximation given (21) is valid after the initial transient is
incorrect.

###  3.1 The invariance equation

The sQSSA, rQSSA, and the tQSSA are, formally, zeroth-order approximations to
the reaction dynamics on the corresponding slow timescale. However, there is
much more to the story. Whenever the tQSSA is valid, trajectories are
attracted to a slow invariant manifold (SIM). The rigorous study of SIMs is
referred to as geometric singular perturbation theory (GSPT). Briefly, given a
system of the form

$\displaystyle\dot{x}$ $\displaystyle=f(x,y;\varepsilon),$ (37a)
$\displaystyle\varepsilon\dot{y}$ $\displaystyle=g(x,y;\varepsilon),$ (37b)

where $x$ is a slow variable and $y$ is a fast variable, one can approximate
the slow manifold by first assuming that the fast variable is expressible in
terms of the slow variable as $y=h(x)$444The underlying assumption here is
that the slow manifold can be expressed in terms of a global coordinate chart.
More complicated manifolds, such as an n-dimensional torus, may require an
atlas and a collection of transition maps., where $h$ is a function that is to
be determined. Typically, the curve $y=h(x)$ is denoted as
$\mathcal{M}^{\varepsilon}$, and we will use this notation from here on. Since
the slow manifold is invariant, it must satisfy the differential equation
(37a)–(37b). Hence, $y=h(x)$ must satisfy

$\cfrac{\text{d}h}{\text{d}x}\cdot
f(x,h(x);\varepsilon)=\cfrac{1}{\varepsilon}g(x,h(x);\varepsilon),$ (38)

which is known as the invariance equation. The power of the invariance
equation is that, given a scaled equation and an appropriate small parameter,
it allows us to systematically determine the higher-order asymptotic
approximations to the SIM. The function $h(x)$ admits a unique asymptotic
expansion in terms of $\varepsilon$:

$h(x)=h_{0}(x)+\varepsilon
h_{1}(x)+\varepsilon^{2}h_{2}(x)+\mathcal{O}(\varepsilon^{3}),$ (39)

where $h_{0}(x)$ satisfies $g(y,h_{0}(x);0)=0$. Thus, when $g(x,y)$ is not
explicitly dependent on $\varepsilon$, the nullcline associated with the fast
variable commonly serves as a zeroth-order approximation to the slow manifold,
$\mathcal{M}^{\varepsilon}_{0}=\\{(x,y)\in\mathbb{R}^{2}:y=h_{0}(x)\\}$, and
formally the asymptotic approximation to the slow manifold can be expressed
as:

$\mathcal{M}^{\displaystyle\varepsilon}=\mathcal{M}_{0}^{\displaystyle\varepsilon}+\varepsilon\mathcal{M}_{1}^{\displaystyle\varepsilon}+\varepsilon^{2}\mathcal{M}_{2}^{\displaystyle\varepsilon}+...+\mathcal{O}(\varepsilon^{N})+...$
(40)

The terms $\varepsilon^{0},\varepsilon^{1},\varepsilon^{2},..$ in (40) are
called gauge functions, and the terms
$\mathcal{M}_{0}^{\displaystyle\varepsilon},\mathcal{M}_{1}^{\displaystyle\varepsilon},...$
are the coefficients that accompany the gauge functions. Generally speaking,
the slow manifold is not necessarily unique, meaning there can be slow
manifolds. However, they are separated by a distance that is
$\mathcal{O}(e^{-K/\displaystyle\varepsilon})$ where $K$ is $\mathcal{O}(1)$
[53].

In addition to the slow manifold, the work of Fenichel also guarantees the
existence of fast fibers, $\mathcal{F}^{\displaystyle\varepsilon,p}$, which
can be utilized to approximate the dynamics during the transient stage as the
trajectory rapidly approaches the slow manifold.555 In contrast to the
critical manifold, the individual zeroth-order fast fibers to not perturb to
individual invariant fibers when $\varepsilon$ is non-zero. However, the
collective family of fast fibers,
$\displaystyle\bigcup_{p}\mathcal{F}^{\displaystyle\varepsilon,p}$, is
invariant in a certain way. For more details, we invite the reader to consult
[54]. The zeroth-order fiber, $\mathcal{F}_{0}^{\displaystyle\varepsilon,p}$,
which is the zeroth-order approximation to the fibers
$\mathcal{F}^{\displaystyle\varepsilon,p}$, is the solution to the fast
subsystem:

$\displaystyle x^{\prime}$ $\displaystyle=0,$ (41a) $\displaystyle y^{\prime}$
$\displaystyle=f(x,y;0).$ (41b)

The fast subsystem contains a continuous branch of equilibrium points given by
the curve $y=h_{0}(x)$, which is referred to in this context as the critical
manifold. Thus, each zeroth-order fiber links an initial condition
$(x_{0},y_{0})$ to a unique base point (equilibrium point)
$p\in\mathcal{M}_{0}^{\displaystyle\varepsilon}$, and the superscript $p$ in
$\mathcal{F}^{\displaystyle\varepsilon,p}$ is a reference to the unique base
point, $p$. Since we are only interested in the zeroth-order approximation to
$\mathcal{F}^{\displaystyle\varepsilon,p}$, we will henceforth drop the
superscript $p$ for convenience, and use
“$\mathcal{F}_{0}^{\displaystyle\varepsilon}$” to denote an individual zeroth-
order fiber. We will say more about fast fibers and the perspective of the
slow manifold from the context of the fast subsystem when we introduce the
idea of Tikhonov-Fenichel parameters in Section 4.2. However, for now we will
just associate the zeroth-order fast fiber as being the solution to
(41a)–(41b).

###  3.2 Standard Quasi-Steady-State dynamics: Geometric singular perturbation
theory when $\varepsilon_{SS}\ll 1$.

The RSA ensures negligible depletion of substrate during the fast transient,
and suggests666 This is, of course, a violation of the conservation law for
substrate in the MM reaction mechanism (1), but we will proceed as if it were
permissible, and discuss this observation in more detail in section 4. that
the zeroth-order fast fiber,
$\mathcal{F}_{0}^{\displaystyle\varepsilon_{SS}}$, is a straight line that
connects the initial condition $(\bar{s},\bar{c})=(1,0)$ to the base point
$(\bar{s},\bar{c})=(1,1)$:

$\mathcal{F}_{0}^{\displaystyle\varepsilon_{SS}}=\\{(\bar{s},\bar{c})\in\mathbb{R}^{2}:\bar{s}=1,0\leq\bar{c}\leq
1\\}.$ (42)

If $0<\varepsilon_{SS}\ll 1$, then the time it takes the phase-plane
trajectory to reach a $\mathcal{O}(\varepsilon_{SS})$-neighborhood of the slow
manifold is $\tau_{\varepsilon_{SS}}\sim|\ln\varepsilon_{SS}|$ or, in
dimensional time $t\sim t_{C}|\ln\varepsilon_{SS}|$. Moreover, there exists an
invariant manifold, $\mathcal{M}^{\displaystyle\varepsilon_{SS}}$, that is
$\mathcal{O}(\varepsilon_{SS})$ from the critical manifold, and admits a
unique expansion in powers of $\varepsilon_{SS}$ in the form of (40), where
$\mathcal{M}_{0}^{\varepsilon_{SS}}$ is identically the $\bar{c}$-nullcline.
These terms are straightforward to compute by constructing an asymptotic
solution to the invariance equation:

$\cfrac{1}{\displaystyle\varepsilon_{SS}}g(\mathcal{M}^{\displaystyle\varepsilon_{SS}},\bar{s})=D_{\bar{s}}\mathcal{M}^{\displaystyle\varepsilon_{SS}}f(\mathcal{M}^{\displaystyle\varepsilon_{SS}},\bar{s}),$
(43)

and by inserting (40) into (43), it is easy to recover

$\mathcal{M}_{0}^{\varepsilon_{SS}}=\\{(\bar{s},\bar{c})\in\mathbb{R}^{2}_{+}:\bar{c}-\bar{s}(1+\sigma)/(1+\bar{s}\sigma)=0\\}.$
(44)

Finally, inserting the expression for $\mathcal{M}_{0}^{\varepsilon_{SS}}$
into (17a) yields, in dimensional variables, the MM equation

$\cfrac{{\rm d}s}{{\rm d}t}\sim-\cfrac{Vs}{K_{M}+s},$ (45)

which approximates the long-time dynamics of the mass action equations when
the RSA is valid.

###  3.3 Extended Quasi-Steady-State dynamics: Geometric singular perturbation
theory when $\nu\ll 1$, $\beta\sim 1$ and $\varepsilon\sim 1$.

Now let us consider the case of the extended QSSA. The problem can be
reformulated in the $(p,c)$ coordinate system so that it is of the form
(37a)–(37b) (see, for example, [18, 55]). However, we will work out the
solution in the $(s,c)$ coordinate system. In this case, we can guess the
asymptotic expansion of the slow manifold in this coordinate system, and
assume it takes the form

$c=\mathcal{M}^{\displaystyle\tilde{\nu}}=\cfrac{\bar{s}}{\alpha\beta+\mu\bar{s}}+\nu\mathcal{M}^{\displaystyle\tilde{\nu}}_{1}+\mathcal{O}(\displaystyle\tilde{\nu}^{2}),\quad\tilde{\nu}\equiv
1/\kappa=k_{2}/k_{-1},$ (46)

where the curve

$\bar{c}=\cfrac{\bar{s}}{\alpha\beta+\mu\bar{s}},$ (47)

is the dimensionless form of the $s$-nullcline (19). Taking just the zeroth-
order expansion of (46) (i.e., (47)) and inserting it into equation (18a)
while zeroing all terms of order $\nu$ or higher yields

$\cfrac{\text{d}\bar{s}}{\text{d}T}=0.$ (48)

Since the $\bar{c}$\- and $\bar{s}$-nullclines coalesce when $k_{2}=0$, the
zeroth-order approximation to the dynamics on the slow manifold results in the
flow being infinitely slow (i.e., an equilibrium solution). Consequently, we
must go to first-order in $\tilde{\nu}$ to determine the QSS approximation to
the dynamics on the slow manifold. Solving for
$\mathcal{M}^{\displaystyle\tilde{\nu}}_{1}$ by approximating the solution to
the invariance equation777$D_{\bar{s}}$ denotes differentiation with respect
to $\bar{s}$.

$g(\mathcal{M}^{\displaystyle\tilde{\nu}},\bar{s})=D_{\bar{s}}\mathcal{M}^{\displaystyle\tilde{\nu}}f(\mathcal{M}^{\displaystyle\tilde{\nu}},\bar{s})$
(49)

in powers of $\tilde{\nu}$, we recover

$\cfrac{\text{d}\bar{s}}{\text{d}T}=\tilde{\nu}\cfrac{((\bar{s}-1)\theta+1)(\theta-1)\bar{s}}{(\bar{s}-1)^{2}\theta^{2}+(2\bar{s}-2\varepsilon-1)\theta+\varepsilon+1}+\mathcal{O}(\tilde{\nu}^{2}),\quad\theta=\cfrac{\sigma}{\alpha+\sigma}<1,$
(50)

by inserting $\mathcal{M}^{\displaystyle\tilde{\nu}}_{1}$ into (18a), which is
the leading-order non-trivial solution to the dynamics on the slow manifold.
Thus, the QSS approximation for $s$ (50) is different when $\nu\ll 1$ versus
when $\varepsilon_{SS}\ll 1$. In dimensional variables, (50) translates to,

$\cfrac{\text{d}s}{\text{d}t}\sim-\cfrac{Vs(s+K_{S})}{e_{0}K_{S}+(s+K_{S})^{2}},\quad
K_{S}\equiv k_{-1}/k_{1}.$ (51)

The result (51) is not new, and the approximation that can be found in several
papers [56, 57, 58, 43, 55]. The point is that by utilizing the invariance
equation, we have shown that in fact Segel and Slemrod’s extended QSS
approximation (21) is invalid. The difference between the QSS approximations
for $s$ when $\nu\ll 1$ versus when $\varepsilon_{SS}\ll 1$ has a nice
history, and Roussel [55] recently published a thorough explanation of the
different QSS approximations through the lens of Tikhonov’s Theorem and
geometric singular perturbation theory.

In addition to the change in leading-order dynamics on the slow manifold when
$\nu\ll 1$ instead of $\varepsilon_{SS}\ll 1$, the zeroth-order approximation
to the fast fiber, $\mathcal{F}_{0}^{\nu}$, will also change. One major
difference between the transient dynamics of the sQSSA and the extended QSSA
is that we can no longer take $s\approx s_{0}$ during the trajectory’s initial
approach to the slow manifold. Hence, we must determine an appropriate
starting point in order to “match” the inner and outer solutions. Locating an
appropriate starting position on a slow manifold is not necessarily a trivial
exercise. However, mathematicians have developed methods for such a situation
(in [59], A. J. Roberts develops a clever method for estimating starting
positions on centre manifolds). Luckily, in our case, the conservation of
total substrate, $s_{0}=s+c+p$, can be utilized to estimate a starting point.
Over the $\tau$-timescale, the total system is given by:

$\displaystyle\cfrac{\text{d}\bar{s}}{\text{d}\tau}$
$\displaystyle=\varepsilon_{SS}[\mu\bar{c}\bar{s}-\bar{s}+\alpha\beta\bar{c}],$
(52a) $\displaystyle\cfrac{\text{d}\bar{c}}{\text{d}\tau}$
$\displaystyle=\bar{s}-\mu\bar{c}\bar{s}-\beta\bar{c},$ (52b)
$\displaystyle\cfrac{\text{d}\bar{p}}{\text{d}\tau}$
$\displaystyle=\epsilon\bar{c}.$ (52c)

Since $\nu\to 0$ implies $\epsilon\to 0$, we set $\epsilon=0$ and $\alpha=1$
in (52c) and (52a), respectively. Setting the right hand side of (52c) to zero
implies that the total substrate, “$\bar{s}_{T}=\bar{s}+\bar{c},$” is
conserved over the transient timescale when $\nu=0$,

$\bar{s}+\bar{c}=1,$ (53)

and insertion of (53) into the scaled equations yields a Riccati equation for
$\bar{c}$:

$\cfrac{\text{d}\bar{c}}{\text{d}\tau}=1-2\bar{c}+\mu\bar{c}^{2}.$ (54)

Two observations can be made. First, the correction to the duration of the
fast transient when $\nu\ll 1$ (as opposed to $\varepsilon_{SS}$) can be found
by solving the Riccati equation (54). In dimensional form, the solution to
(54) yields (33). Again, this timescale was originally introduced by Tzafriri
[34], and can be taken to be characteristic of the fast transient whenever
there is minimal product formation during the time it takes $c$ to reach its
maximum value (i.e., even when $\varepsilon_{SS}\ll 1$).888The usage of
“characteristic” is slightly abused in this context since the Riccati equation
is nonlinear. Second, the base point $(\bar{c}^{*},\bar{s}^{*})$ of the
zeroth-order fast fiber is given by999 The zeroth-order fast fiber is not a
vertical line in the $(s,c)$ phase–plane. However, it is a vertical line in
the $(p,c)$ phase–plane, since $p$ is the slow variable in this case and not
$s$.

$\bar{c}^{*}=\cfrac{2-\sqrt{4(1-\mu)}}{2\mu},\quad\bar{s}^{*}=1-\bar{c}^{*}.$
(55)

Consequently, if $\varepsilon_{SS}=\sigma=1$ then, after the fast transient,
we expect $\bar{c}\approx 2-\sqrt{2}$ and $\bar{s}\approx\sqrt{2}-1$ (see
figure 4). In dimensional variables, this corresponds to
$c\approx(2-\sqrt{2})\,s_{0}$ and $s=(2-\sqrt{2})\,s_{0}$, and we recover the
original estimate given by Segel and Slemrod [30] (see figure 4).

Figure 4: The extended QSSA equation (21) proposed by Segel and Slemrod [30]
is invalid when $\nu\ll 1$ but $\varepsilon_{SS}=\sigma=1$. In both panels,
the solid black curve is the numerical solutions to the mass action equations
(5a)–(5c), and the dashed/dotted red curve is the $c$-nullcline. left: In this
panel, $s_{0}=20,e_{0}=1,k_{2}=k_{-1}=10,k_{1}=0.1$, and the sQSSA is valid,
and the trajectory moves almost vertically towards the $c$-nullcline, which is
the hallmark of the sQSSA in the $(s,c)$ phase-plane. right: This panel is
illustrates the path of a trajectory when $k_{2}\ll k_{-1}$ but
$\varepsilon_{SS}=1$. The parameters utilized in the numerical simulation are:
$k_{1}=k_{-1}=1,k_{2}=0.01$ with $\kappa=100$ and $\sigma=1$. Notice the
trajectory is no longer vertical during the initial phase of the reaction as
there is substantial loss of substrate during the fast transient. Once the
trajectory reaches the $c$-nullcline, the approximation given by (51) is
valid, and can be equipped with the initial condition $(\sqrt{2}-1)s_{0}$ if
$\varepsilon_{SS}=\sigma=1$. Although the sQSSA looks valid in both panels
based on the proximity of the trajectory to the $c$-nullcline, GSPT indicates
it is only valid in the left panel. In this figure, double arrows represent
fast dynamics, and single arrows represent slow dynamics. All units are
arbitrary.

GSPT, or really just the invariance equation, provides a powerful framework
from which QSS approximations can be carefully derived. Moreover, it can help
to rule out erroneous conclusions about specific QSS approximations such as
(21). Note that the utility of the invariance equation still relies on direct
knowledge of a suitable small parameter in order to construct an asymptotic
expansion to the slow manifold. However, the appropriate small parameter that
warrants the validity of the rQSSA is still an open problem. This raises the
question: are there suitable non-scaling methods that can be employed to
determine the validity of the rQSSA? We introduce two such methods in the
section that follows.

##  4 Applying the Quasi-Steady-State Approximations to the Michaelis–Menten
reaction mechanism: Non-scaled approaches

In this section, we introduce non-scaling approaches to finding parameters or
combinations of parameters that, when made very small, justify the validity of
the sQSSA, the extended QSSA or the rQSSA. We also remark on the origins of
slow manifolds, as this will be critical to uncovering sufficient conditions
for the validity of the rQSSA in Section 5.

###  4.1 Geometric singular perturbation theory: The fast subsystem

Let us consider once more a fast/slow system of the form

$\displaystyle x^{\prime}$ $\displaystyle=\varepsilon f(x,y;\varepsilon),$
(56a) $\displaystyle y^{\prime}$ $\displaystyle=g(x,y;\varepsilon),$ (56b)

which we will refer to as being in Tikhonov standard form. The associated fast
subsystem,

$\displaystyle x^{\prime}$ $\displaystyle=0,$ (57a) $\displaystyle y^{\prime}$
$\displaystyle=g(x,y;0),$ (57b)

contains a branch of fixed points $y=h_{0}(x)$, where $g(x,h_{0}(x);0)=0$.
Again, each fixed point in the set $y=h_{0}(x)$ is called a base point, and
the complete set $y=h_{0}(x)$, denoted hereby as $\mathcal{M}$, is said to be
normally-hyperbolic provided

$\mathfrak{Re}\bigg{(}\cfrac{\partial g(x,y;0)}{\partial
y}\bigg{|}_{y=h_{0}(x)}\bigg{)}\neq 0.$ (58)

In this context, normal hyperbolicity is partially101010 Normally
hyperbolicity provides information about the dominant directions of the flow
near the critical manifold. a statement about the linear stability along
directions that are transverse to the critical manifold. As previously stated,
the zeroth-order fast fibers are solutions to (56a)–(56b), but they are also
individual manifolds themselves. To see this, note that if (58) is less than
zero, then $\mathcal{M}$ is attracting, and trajectories that start close to
$\mathcal{M}$ move closer to $\mathcal{M}$ as $t\to\infty$ (see Figure 5). In
contrast, if (58) is positive, then $\mathcal{M}$ is repelling, and
trajectories that start close to $\mathcal{M}$ move away from $\mathcal{M}$ as
$t\to\infty$. Thus, when the critical manifold is normally hyperbolic and
attracting, the zeroth-order fast fiber is the stable manifold corresponding
to the base point “$p$” of $\mathcal{M}$, and the stable manifold of
$\mathcal{M}$ is comprised of (foliated by) the entire union of zeroth-order
fast fibers.

Figure 5: Illustration of a normally hyperbolic and invariant critical
manifold, $\mathcal{M}$. The black lines with arrows and fixed points
represent trajectories within the unperturbed vector field field; the red
dashed lines represent trajectories in the perturbed vector field. The
straight black lines with double arrowheads represent trajectories moving
towards the fixed points (filled black circles) that lie on the normally
hyperbolic, attracting and invariant manifold $\mathcal{M}$, depicted by the
thick black line. The flow on $\mathcal{M}$ is trivial. However, once
$\varepsilon$ is non-zero and the vector field is smoothly perturbed, an
invariant manifold $\mathcal{M}^{\varepsilon}$ emerges (red), on which the
flow is slow, but no longer trivial. Trajectories that start off of
$\mathcal{M}^{\varepsilon}$ quickly approach it, then follow it tangentially.
Obviously, the fixed points drawn on this manifold are isolated, so the reader
should bear in mind that $\mathcal{M}$ is filled with non-isolated fixed
points.

Why are we interested in manifolds that are comprised of fixed points?
Fenichel’s Theorems describe what happens to normally hyperbolic critical
manifolds of fixed points once the vector field is perturbed. Formally, we
introduce the following theorem due to Fenichel [38], which can be found in
[53]:

Theorem 1: Let $F$ be a $C^{r}$ vector field on $\mathbb{R}^{n}$ with $r\geq
1$. Let $\mathcal{M}$ be a compact and connected $C^{r}$ manifold embedded in
$\mathbb{R}^{n}$. Suppose that $\mathcal{M}$ is normally hyperbolic and
invariant under the flow of $F$. Then, given any $C^{r}$ vector field
$F^{\varepsilon}$ that is sufficiently $C^{1}$–close to $F$, there exists a
$C^{r}$–manifold $\mathcal{M}^{\varepsilon}$ that is invariant under the flow
of $F^{\varepsilon}$ and diffeomorphic to $\mathcal{M}$.

The notation “$C^{r}$” means the vector field has at least $r$ derivatives,
all of which are continuous. The terms “compact” and “connected” are referring
to specific topological properties that $\mathcal{M}$ must be equipped with.
Connected means that $\mathcal{M}$ cannot be expressed as the union of two
disjoint and non-empty sets, and compact means every open cover has a finite
sub-cover. The term $C^{1}$-close is a statement about the distance between
$F$ and $F^{\varepsilon}$:

Definition: Let $F$ and $F^{\varepsilon}$ be two $C^{1}$ vectors fields on
$\mathbb{R}^{n}$, and $\Lambda$ be a compact set. $F$ is said to be
$\vartheta$-close to $F^{\varepsilon}$ on $\Lambda$ if:

$\displaystyle\sup_{\boldsymbol{x}\in\Lambda}$
$\displaystyle||F-F^{\varepsilon}||\leq\vartheta,$ (59a)
$\displaystyle\sup_{\boldsymbol{x}\in\Lambda}$ $\displaystyle||DF-
DF^{\varepsilon}||\leq\vartheta,$ (59b)

where $D$ denotes differentiation. Theorem 1 is one of several theorems that
provides the foundation for GSPT (see [54] for a more technical survey). Later
theorems address compact manifolds with a boundary, in which case invariance
is replaced with local invariance. In a nutshell, one starts with a vector
field that contains a slow, compact critical manifold that consists of fixed
points. Fenichel’s theorems tells us that if a compact subset of the critical
manifold is normally hyperbolic and attracts nearby trajectories, then
perturbing the vector field (i.e., allowing $\varepsilon$ in (57a)–(57b) to be
greater than zero) in an appropriate way results in an attracting slow
manifold.

Let us now discuss a specific point, “$p\in\mathcal{M}$”, where (58) fails and

$\cfrac{\partial g(x,y;0)}{\partial y}\bigg{|}_{p\in\mathcal{M}}=0.$ (60)

A point that satisfies (58) is called a singular point. Singular points can
indicate where a change in the stability of the critical manifold
$\mathcal{M}$ may occur, and can correspond to dynamic bifurcation points. In
particular, a transcritical bifurcation occurs at a singular point where two
normally hyperbolic submanifolds cross and undergo an exchange of stability
(see Figure 6). This specific bifurcation will be of interest when we analyze
the rQSSA.

Figure 6: The transcritical bifurcation point corresponds to a loss of normal
hyperbolicity and an exchange of stability. Here we illustrate a transcritical
bifurcation. Thick dashed lines corresponds to unstable fixed points in which
(58) is greater than zero. Thick solid lines represent stable branches of
fixed points for which (58) is negative. Two manifolds, $\mathcal{M}_{1}$ and
$\mathcal{M}_{2}$, cross at a singular point. At the crossing point (i.e, the
point of intersection), normal hyperbolicity is lost, and (58) fails to hold.
Each manifold undergoes a change in stability at the point of intersection
(i.e., the singular point delimits the point at which the stability of each of
the manifolds changes from attractive to repulsive or vice-versa).

###  4.2 Tikhonov-Fenichel parameters

We now want to take the idea of a critical manifold that is composed of fixed
points one step further. Revisiting the mass action equations (5a)–(5c), we
can ask the following question: Are there specific dimensional parameters
that, when identically zero, result in the formation of a normally hyperbolic
invariant manifold that is composed of fixed points? If so, then it seems
reasonable to assume that these parameters should in some sense be small
whenever a specific QSS reduction is valid. For example, let us first consider
the rQSSA in the $(s,p)$ phase-plane, for which the mass action equations are:

$\displaystyle\dot{s}$
$\displaystyle=-k_{1}(e_{0}-(s_{0}-s-p))s+k_{-1}(s_{0}-s-p),$ (61a)
$\displaystyle\dot{p}$ $\displaystyle=k_{2}(s_{0}-s-p).$ (61b)

We know that the rQSSA implies that $s\approx 0$ during the QSS phase of the
reaction. Define $k_{2}=\varepsilon_{2}k_{2}^{\star}$ and
$k_{-1}=\varepsilon_{1}k_{-1}^{\star}$, where $\varepsilon_{1}$ and
$\varepsilon_{2}$ are dimensionless and small,
$0<\varepsilon_{1},\varepsilon_{2}\ll 1$, and $k_{-1}^{\star},k_{2}^{\star}$
are identically one with the same units as $k_{-1}$ and $k_{2}$ respectively.
Expressing (61a)–(61b) as

$\displaystyle\dot{s}$
$\displaystyle=-k_{1}(e_{0}-(s_{0}-s-p))s+\varepsilon_{1}k_{-1}^{\star}(s_{0}-s-p),$
(62a) $\displaystyle\dot{p}$
$\displaystyle=\varepsilon_{2}k_{2}^{\star}(s_{0}-s-p),$ (62b)

it is clear that the curve $s=0$ is filled with equilibrium points when
$\varepsilon_{1}=\varepsilon_{2}=0$ and $e_{0}>s_{0}$. Furthermore, it is
straightforward to show that the curve $s=0$ is normally hyperbolic and
attracting when $s_{0}<e_{0}$, since

$\cfrac{\partial}{\partial
s}\bigg{[}-k_{1}(e_{0}-(s_{0}-s-p))s\bigg{]}\bigg{|}_{s=0}<0.$ (63)

Consequently, we anticipate

$\dot{p}\sim k_{2}(s_{0}-p)$ (64)

during the slow QSS phase of the reaction when
$0<\varepsilon_{1},\varepsilon_{2}\ll 1$. The takeaway message from this
example is two-fold. First, without scaling and non-dimensionalization, we
have predicted the validity of the rQSSA by choosing dimensional parameters
(i.e., $k_{-1}$ and $k_{2}$) that, upon vanishing, give rise to an invariant
manifold of equilibrium points. This is consistent with the notion of a
normally hyperbolic critical manifold from fast/slow analysis and GSPT.
Second, any dimensionless parameter homologous to $\varepsilon_{SS}$ that acts
as a qualifier for the validity of the rQSSA should vanish whenever $k_{-1}$
and $k_{2}$ vanish. In planar systems, a parameter that, upon vanishing, gives
rise to a manifold of fixed points that is normally hyperbolic, is called a
Tikhonov-Fenichel parameter (TFP); Walcher and collaborators were the first to
conceptualize TFPs [37, 58, 35, 43]. A manifold of equilibrium points is a
special case of a normally hyperbolic manifold; the formal definition of a
normally hyperbolic invariant manifold can by found in [53]. The informal
definition of a normally hyperbolic invariant manifold is that the linearized
flow that is tangent to the manifold is dominated by the flow that is normal
(transversal) to the manifold. This is more or less intuitive when it comes to
manifolds comprised of equilibria, since the flow on such a manifold is
trivial. Nevertheless, manifolds of equilibria can lose normal hyperbolicity
and, as we will see, the loss can play a significant role in determining when
QSS approximations are valid.

The challenging part in the direct utilization of Tikhonov-Fenichel parameters
is that often times the perturbed mass action equations obtained by regulating
a Tikhonov-Fenichel parameter do not necessarily result in a set of equations
that are in Tikhonov standard form. This can make computing a QSS reduction
(without resorting to scaling) less than obvious. For example, the parameters
$e_{0}$ and $k_{1}$ are Tikhonov-Fenichel parameters that warrant the validity
of the sQSSA [35]. Treating $k_{1}$ as a small parameter, we can express
(5a)–(5b) as

$\displaystyle\dot{s}$ $\displaystyle=-\varepsilon
k_{1}^{\star}(e_{0}-c)s-k_{-1}c,$ (65a) $\displaystyle\dot{c}$
$\displaystyle=\varepsilon k_{1}^{\star}(e_{0}-c)s+(k_{-1}+k_{2})c.$ (65b)

The $(s,c)$ phase-plane is equipped with an invariant manifold of equilibria,
$c=0$, when $k_{1}=0$. If $0<k_{1}\ll 1$, then one would expect that
$s_{0}\sim s+p$ on the slow timescale, provided the curve $c=0$ is normally
hyperbolic and attracts nearby trajectories. However, determining a QSS
approximation for $s$ or $c$ from (65a)–(65b) is not obvious without scaling,
since the system is not in Tikhonov standard form. Nonetheless, using
algebraic methods, Walcher and collaborators [43, 36, 57] developed a
surprisingly straightforward method for obtaining QSS approximations from non-
standard systems such as (65a)–(65b). In so doing, they have circumnavigated
the need for scaling analysis to justify the QSSA.

For simple planar systems, Tikhonov-Fenichel parameters provide direct insight
into the conditions that contribute to QSS dynamics, as well the origin (i.e.,
the critical manifold) of the corresponding slow manifold, and in this sense
they are invaluable. However, in order to determine the validity of a
particular QSS reduction so that experiments can be prepared with the
intention of determining enzyme activity, the necessary “smallness” of
Tikhonov-Fenichel parameters must be defined. For example, if $s_{0}\geq
e_{0}$, then the critical set obtained by setting
$\varepsilon_{1}=\varepsilon_{2}=0$ in (62a)–(62b) yields

$e_{0}s-(s_{0}-s-p)s=0,$ (66)

since the term $e_{0}-(s_{0}-s-p)$ can vanish if $e_{0}\leq s_{0}$. In this
situation, it is not entirely clear how to go about properly defining the
critical manifold when $s_{0}\geq e_{0}$. This raises the question as to
whether or not the rQSSA is still valid when $e_{0}\lesssim s_{0}$ and, if so,
then how small must the parameters $k_{-1}$ and $k_{2}$ be in order for the
approximation (64) to hold over the slow timescale? Obviously we expect the
standard QSSA to be valid at extremely large values of $s_{0}$, but what
happens in parameter regions where $s_{0}$ and $e_{0}$ are the same order of
magnitude? This question is central not only to the validity of the rQSSA, but
is also central to understanding the validity of the sQSSA. Clearly, it seems
that something more than scaling analysis and the identification of Tikhonov-
Fenichel parameters is needed in order to address the validity of the rQSSA.
We will introduce a method that permits a possible solution to the validity of
the rQSSA in the subsection that follows.

###  4.3 Energy methods: Revisiting the sQSSA and Michaelis–Menten equation

To obtain a fundamental small parameter that ensures the validity of the
sQSSA, we need to understand where the corresponding slow manifold comes from
and, to do this, we must determine the associated critical manifold. One
possible caveat with the traditional analysis of the sQSSA in the $(s,c)$
phase–plane is that, if we are only interested in bounded trajectories (again,
this rules out the particular case of taking $s_{0}\to\infty$), then taking
$\varepsilon\to 0$ implies that either $k_{1}$ or $e_{0}$ must vanish, as
Walcher and his collaborators have carefully pointed out [43, 58, 35, 37]. In
both of these cases $c$ must vanish, but the dimensionless variable
$\bar{c}=c/\varepsilon s_{0}$ contains a zero denominator in the singular
limit. According to Walcher and Lax [60] this is not a problem. However, in
order to avoid having to deal with this, it is useful to write down the
dimensionless equations in the $(\bar{s},\bar{p})$ coordinate system:

$\displaystyle\cfrac{\text{d}\bar{s}}{\text{d}\tau}$
$\displaystyle=-\varepsilon_{SS}\bar{s}+\mu\bar{s}(1-\bar{s}-\bar{p})+\alpha\beta(1-\bar{s}-\bar{p})$
(67a) $\displaystyle\cfrac{\text{d}\bar{p}}{\text{d}\tau}$
$\displaystyle=\beta(1-\alpha)(1-\bar{s}-\bar{p}).$ (67b)

Setting $\varepsilon_{SS}=0$ in (67a)–(67b) reveals the critical manifold to
be $\bar{p}=1-\bar{s}$, which corresponds to the critical manifold of the
dimensional system (i.e., $c=0$) obtained by setting $k_{1}$ or $e_{0}$ to
zero in (65a)–(65b). Since the curve $s=s_{0}-p$ implies $c=0$, we avoid the
other caveat encountered in the analysis in the $(s,c)$ coordinate system in
section 3: there is no violation of the substrate conservation law
$s_{0}=s+c+p$ along a zeroth-order fast fiber, since initial conditions lie on
the critical manifold, $c=0$. If a non-trivial $c$-nullcline is defined to be
the critical manifold, then there is zero loss of substrate along the one-
dimensional zeroth-order fibers when $\varepsilon_{SS}=0$, which is a
violation of the substrate conservation law, since $s+c>s_{0}$ as the
trajectory approaches a base point that lies above the $s$-axis.

It is expected that trajectories will closely follow the curve $p=s_{0}-s$
once the perturbation is “turned on” and $0<\varepsilon_{SS}\ll 1$. To get an
idea of how “good” the approximation $p\sim s_{0}-s$ is for all time, we need
to compute an upper bound “$M$” such that111111We do not need to write
$|\mathcal{E}|$ since $0\leq s_{0}-s-p$, but later on we will consider cases
where the quantity of interest can vary in sign, so it helps to keep the
notation consistent through each procedure.

$(1-\bar{s}-\bar{p})\leq M.$ (68)

We will start with a differential equation for the energy121212This is the
mathematical energy of a function, and should not be confused with the
physical energy of a particle. $\mathcal{E}^{2}\equiv(s_{0}-s-p)^{2}$.
Multiplying both sides of the dimensional form of (67a) by
“$(s_{0}-s-p)\equiv\mathcal{E}$” yields:

$\displaystyle\cfrac{1}{2}\cfrac{\text{d}\mathcal{E}^{2}}{\text{d}t}$
$\displaystyle=k_{1}e_{0}s\mathcal{E}-k_{1}(s+K_{M})\mathcal{E}^{2},$ (69a)
$\displaystyle\leq k_{1}e_{0}s_{0}|\mathcal{E}|-k_{1}K_{M}\mathcal{E}^{2}.$
(69b)

Now we can can use Cauchy’s inequality with “$\delta$,”

$ab\leq\delta a^{2}+\cfrac{b^{2}}{4\delta},$ (70)

and expand the term $k_{1}e_{0}s_{0}\mathcal{E}$ in (69b):

$k_{1}e_{0}s_{0}|\mathcal{E}|\leq\cfrac{(k_{1}e_{0}s_{0})^{2}}{4\delta}+\delta\mathcal{E}^{2}.$
(71)

Choosing $\delta=k_{1}K_{M}/2$ yields

$k_{1}e_{0}s_{0}|\mathcal{E}|\leq\cfrac{(k_{1}e_{0}s_{0})^{2}}{2k_{1}K_{M}}+\cfrac{k_{1}K_{M}}{2}\mathcal{E}^{2}.$
(72)

After combining (72) with (69b) we have

$\cfrac{\text{d}\mathcal{E}^{2}}{\text{d}t}\leq\cfrac{(k_{1}e_{0}s_{0})^{2}}{k_{1}K_{M}}-k_{1}K_{M}\mathcal{E}^{2},$
(73)

and integrating both sides of (73) yields

$\mathcal{E}^{2}(t)\leq\mathcal{E}^{2}(0)e^{\displaystyle-(k_{1}K_{M})t}+s_{0}^{2}\eta^{2}(1-e^{-\displaystyle
k_{1}K_{M}t}).$ (74)

Replacing “$(1-e^{-k_{1}K_{M}t})$” with “$1$” and dividing both sides of (74)
by $s_{0}^{2}$ reveals

$\displaystyle\limsup_{t\to\infty}\mathcal{\bar{E}}^{2}\leq\eta^{2},\quad\mathcal{\bar{E}}\equiv\mathcal{E}/s_{0},$
(75)

where the limit supremum or “$\limsup$” of a bounded function $w(t)$ is given
by:

$\limsup_{t\to\infty}w(t)=\displaystyle\lim_{T\to\infty}\bigg{(}\sup_{t>T}w(t)\bigg{)}.$
(76)

It is clear from (75) that the long-time error in the approximation is
controlled by $\eta$. We can simplify things further: utilizing the fact that
$\mathcal{E}(0)=0$, and taking the square root of both sides yields

$(1-\bar{s}-\bar{p})\leq\eta.$ (77)

The bound (77) obtained from the energy analysis is promising, but what we are
ultimately after is an expression for $p(t)$ that can be utilized for
parameter estimation. In practice, the approximation131313This approximation
is easily justified through scaling analysis. See [16] for details.

$\cfrac{\text{d}p}{\text{d}t}\sim\cfrac{k_{2}e_{0}(s_{0}-p)}{K_{M}+(s_{0}-p)},$
(78)

is often used to estimate $K_{M}$ and $k_{2}$; the assumption in (78) is that
$c$ is negligible and is well approximated by

$c\approx\cfrac{e_{0}(s_{0}-p)}{K_{M}+s_{0}-p}.$ (79)

The question that immediately follows is: how good is the approximation (79)?
In order to assess the validity of (79), we will derive an upper bound on
$|\mathcal{E}_{c}|$,

$\mathcal{E}_{c}\equiv\bigg{(}c-\cfrac{e_{0}(s_{0}-p)}{K_{M}+s_{0}-p}\bigg{)}.$
(80)

To begin,

$\displaystyle\cfrac{1}{2}\cfrac{\text{d}\mathcal{E}_{c}^{2}}{\text{d}t}$
$\displaystyle=\bigg{(}c-\cfrac{e_{0}(s_{0}-p)}{K_{M}+s_{0}-p}\bigg{)}\bigg{(}\cfrac{\text{d}c}{\text{d}t}+\cfrac{e_{0}K_{M}}{(K_{M}+s_{0}-p)^{2}}\cfrac{\text{d}p}{\text{d}t}\bigg{)},$
(81a)
$\displaystyle=\mathcal{E}_{c}\bigg{(}k_{1}e_{0}(s_{0}-p)-k_{1}c(e_{0}+K_{M}+s_{0}-p)+k_{1}c^{2}+\cfrac{e_{0}K_{M}k_{2}c}{(K_{M}+s_{0}-p)^{2}}\bigg{)},$
(81b)
$\displaystyle=\mathcal{E}_{c}\bigg{(}-k_{1}(K_{M}+s_{0}-p)\mathcal{E}_{c}+k_{1}c(c-e_{0})+\cfrac{e_{0}K_{M}k_{2}c}{(K_{M}+s_{0}-p)^{2}}\bigg{)},$
(81c)
$\displaystyle=-k_{1}(K_{M}+s_{0}-p)\mathcal{E}_{c}^{2}+\mathcal{E}_{c}\bigg{(}k_{1}c(c-e_{0})+\cfrac{e_{0}K_{M}k_{2}c}{(K_{M}+s_{0}-p)^{2}}\bigg{)}.$
(81d)

Next, we need to find an upper bound on the right and side of (81d):

$\displaystyle\cfrac{1}{2}\cfrac{\text{d}\mathcal{E}_{c}^{2}}{\text{d}t}$
$\displaystyle\leq-
k_{1}K_{M}\mathcal{E}_{c}^{2}+|\mathcal{E}_{c}|\max\bigg{|}k_{1}c(c-e_{0})+\cfrac{e_{0}K_{M}k_{2}c}{(K_{M}+s_{0}-p)^{2}}\bigg{|},$
(82a) $\displaystyle\leq-
k_{1}K_{M}\mathcal{E}_{c}^{2}+|\mathcal{E}_{c}|\max|k_{1}c(c-e_{0})|+|\mathcal{E}_{c}|\max\bigg{|}\cfrac{e_{0}K_{M}k_{2}c}{(K_{M}+s_{0}-p)^{2}}\bigg{|},$
(82b) $\displaystyle\leq-
k_{1}K_{M}\mathcal{E}_{c}^{2}+|\mathcal{E}_{c}|\bigg{(}\cfrac{k_{1}e_{0}^{2}}{4}+\cfrac{e_{0}k_{2}\lambda}{K_{M}}\bigg{)}.$
(82c)

Using Cauchy’s inequality with $\delta=k_{1}K_{M}/2$, it follows that

$\cfrac{\text{d}\mathcal{E}_{c}^{2}}{\text{d}t}\leq-
k_{1}K_{M}\mathcal{E}_{c}^{2}+\cfrac{1}{k_{1}K_{M}}\bigg{(}\cfrac{k_{1}e_{0}^{2}}{4}+\cfrac{e_{0}k_{2}\lambda}{K_{M}}\bigg{)}^{2}.$
(83)

Applying Gronwall’s lemma and integrating both sides of (83) yields

$\mathcal{E}_{c}^{2}\leq\mathcal{E}_{c}^{2}(0)\displaystyle e^{\displaystyle-
k_{1}K_{M}t}+\cfrac{1}{k_{1}^{2}K_{M}^{2}}\bigg{(}\cfrac{k_{1}e_{0}^{2}}{4}+\cfrac{e_{0}k_{2}\lambda}{K_{M}}\bigg{)}^{2}.$
(84)

Although the exponential term (84) is in terms of the dimensional time, $t$,
we expect that if (84) is written in terms of the slow time,
$\bar{T}=tk_{2}\varepsilon_{SS}$, then it should decay relatively quickly in
$\bar{T}$. To do this we need to solve

$k_{1}K_{M}t=C_{\bar{T}}k_{2}\varepsilon_{SS}t$ (85)

for the constant, $C_{\bar{T}}$. Solving (85) reveals
$C_{\bar{T}}=1/\nu\varepsilon$, and it clear that the exponential term in (84)
will vanish rapidly in $\bar{T}$ whenever $\varepsilon\nu\ll 1$. Finally,
taking the square root of both sides while noting that
$\sqrt{x+y}\leq\sqrt{x}+\sqrt{y}\;\;\forall x,y\geq 0$, and dividing through
by $e_{0}$ produces the following upper bound:

$|\tilde{c}-h(\bar{p})|\leq\mu\displaystyle
e^{\displaystyle-\bar{T}/2\varepsilon_{SS}\nu}+\cfrac{\eta}{4}+\nu\cfrac{\lambda}{K_{M}},\quad
h(\bar{p})\equiv\cfrac{\sigma(1-\bar{p})}{1+\sigma(1-\bar{p})},$ (86)

where $\tilde{c}\equiv c/e_{0}$. The bound (86) is actually quite revealing.
First, we immediately see that taking $\nu\to 0$ does not support the validity
of the sQSSA (i.e., the MM equation), even though the exponential term decays
rapidly in the slow time whenever $\varepsilon_{SS}\nu\ll 1$, and “nullcline
following” in $(s,c)$ phase-plane is observed when $k_{2}\ll k_{-1}$. Thus, we
do not encounter any sort of dilemma between the sQSSA and the extended QSSA
when energy methods are appropriately employed. Second, if $s_{0}<e_{0}$, then
$\lambda<s_{0}$, and we have

$|\tilde{c}-h(\bar{p})|\leq\mu\displaystyle
e^{\displaystyle-\bar{T}/2\varepsilon_{SS}\nu}+\cfrac{\eta}{4}+\nu\sigma.$
(87)

Since $\sigma<\eta$ when $s_{0}<e_{0}$, it follows that $\eta\to 0$ in order
for the sQSSA to be valid. Third, note that making $s_{0}$ arbitrarily large
only makes the exponential term decay faster, since $\varepsilon_{SS}\to 0$ as
$s_{0}\to\infty$. However, large $s_{0}$ has no influence on the long-time
bound: only $K_{M}$ and $e_{0}$ influence the long-time bound. Finally, if
$e_{0}<s_{0}$, then $\lambda\leq e_{0}$ and it holds that

$|\tilde{c}-h(\bar{p})|\leq\mu\displaystyle
e^{\displaystyle-\bar{T}/2\varepsilon_{SS}\nu}+\cfrac{5}{4}\eta.\phantom{+}\phantom{\nu\sigma}$
(88)

In all cases the natural “small parameter” that arises from the energy
analysis is $\eta$. What can happen if $\varepsilon_{SS}\ll 1$ but $\eta\gg
1$? The answer is that a curious situation arises: the $(s,c)$ phase-plane
trajectory will rapidly approach the $c$-nullcline, and proceed to follow it
closely for a finite amount of time. However, as the reaction progresses, the
phase-plane trajectory, which follows the invariant manifold
$\mathcal{M}^{\varepsilon}$, will start to move away from the $c$-nullcline,
and eventually follow the $s$-nullcline (see figure 7).

Figure 7: Near the origin, the slow manifold lies close to the $s$-nullcline
in the $(s,c)$ phase–plane when $1\ll\eta$. This figure is an illustration and
is meant to convey conceptual features that occur in the phase–plane when
$\eta\gg 1$. The solid black curve is the slow manifold; the dotted blue curve
is the $s$-nullcline, and the dashed red curve is the $c$-nullcline. When
$\eta$ is sufficiently large, the slow manifold lies very close to the
$s$-nullcline at sufficiently small $s$. Comparing with figure 3, this is
exactly the condition we need in order to apply the rQSSA

In conclusion of this section, our analysis suggests that $\eta$ is the
fundamental small parameter that justifies the validity of the sQSSA. The
condition $\eta\ll 1$ can be viewed as a special case of $\varepsilon_{SS}\ll
1$. However, if $\varepsilon_{SS}\ll 1$ and $1\ll\eta$, then we observe a
scenario in the $(s,c)$ phase-plane in which the sQSSA is initially valid
after the fast transient, but steadily loses validity until the rQSSA becomes
valid for the remainder of the reaction. It seems as though the situation is
quite unique, since the sQSSA appears to valid over the slow timescale, but
the rQSSA is valid over a super-slow timescale. If we can understand the
origin of the slow manifold in this case (i.e., its associated critical
manifold), then perhaps it will be possible to explain why this scenario
occurs. Of course, this raises the question: “Can a unique critical manifold
be determined in this case?” The answer to this question resides in the
analysis of the rQSSA, and we will analyze this particular approximation in
great detail in section 5.

##  5 The reverse Quasi-Steady-State Approximation

In this section, we derive the validity of the rQSSA using a combination of
scaling and energy methods.

###  5.1 Justification of the reverse Quasi-Steady-State Approximation in the
$(p,c)$ phase-plane: Scaling and bifurcation analysis

It is straightforward to verify the validity of the rQSSA through scaling
analysis. Assuming that $s_{0}<e_{0}$, it is obvious that $c$ is bounded by
$s_{0}$, which suggests that $\hat{c}=c/s_{0}$ is the appropriate scaled
variable when $s_{0}<e_{0}$. Rescaling the mass action equations (27a)–(27b)
in the $(p,c)$ phase–plane yields

$\displaystyle\varepsilon^{*}\nu\cfrac{\text{d}\hat{c}}{\text{d}\tilde{T}}$
$\displaystyle=(1-\ell\hat{c})(1-\hat{c}-\bar{p})-\varepsilon^{*}\hat{c},$
(89a) $\displaystyle\cfrac{\text{d}\bar{p}}{\text{d}\tilde{T}}$
$\displaystyle=\hat{c},$ (89b)

and the corresponding critical manifold, $\mathcal{M}_{0}$, is obtained by
setting $\varepsilon^{*}=0$,

$\mathcal{M}_{0}=\\{(\bar{p},\hat{c})\in\mathbb{R}^{2}_{+}|(1-\ell\hat{c})(1-\hat{c}-\bar{p})=0\\}.$
(90)

If $s_{0}<e_{0}$, then $(1-\ell\hat{c})>0$, and we can take the critical
manifold to be $1-\hat{c}-\bar{p}=0$. Furthermore, the critical manifold is
normally hyperbolic and attracting, since

$\cfrac{\partial}{\partial\hat{c}}(1-\ell\hat{c})(1-\hat{c}-\bar{p})<-.$ (91)

Thus, due to the fact that the critical manifold is both normally hyperbolic
and attracting, we expect trajectories to approach the curve $c=s_{0}-p$ when
$s_{0}<e_{0}$ and $0<\varepsilon^{*}\ll 1$. But what happens when $e_{0}\leq
s_{0}$? Let us start by setting $e_{0}=s_{0}$, in which case $\ell=1$. The
dimensionless equations are given by

$\displaystyle\varepsilon^{*}\nu\cfrac{\text{d}\hat{c}}{\text{d}\tilde{T}}$
$\displaystyle=(1-\hat{c})(1-\hat{c}-\bar{p})-\varepsilon^{*}\hat{c},$ (92a)
$\displaystyle\cfrac{\text{d}\bar{p}}{\text{d}\tilde{T}}$
$\displaystyle=\hat{c},$ (92b)

and the associated critical set is

$\mathcal{M}_{0}=\\{(\bar{p},\hat{c})\in\mathbb{R}^{2}_{+}|(1-\hat{c})(1-\hat{c}-\bar{p})=0\\}.$
(93)

The set (93) fails to be a manifold at the point where the curves $\hat{c}=1$
and $\hat{c}=1-\bar{p}$ intersect, which is precisely the point
$(\bar{p},\hat{c})=(0,1)$. Moreover, there is a loss of normal hyperbolicity
at this point, since

$\cfrac{\partial}{\partial\hat{c}}\bigg{[}(1-\hat{c})(1-\hat{c}-\bar{p})\bigg{]}\bigg{|}_{(\bar{p},\hat{c})=(0,1)}=0.$
(94)

What is happening at the point where normal hyperbolicity is lost? The sub-
manifolds $\hat{c}=1$ and $\bar{p}=1-\hat{c}$ intersect and exchange stability
at the point $(\bar{p},\hat{c})=(0,1)$ when $K_{M}=0$. To see this, rescale
time and define $\tau^{*}=\tilde{T}/\varepsilon^{*}\nu$:

$\displaystyle\cfrac{\text{d}\hat{c}}{\phantom{i}\text{d}\tau^{*}}$
$\displaystyle=(1-\hat{c})(1-\hat{c}-\bar{p})-\varepsilon^{*}\hat{c},$ (95a)
$\displaystyle\cfrac{\text{d}\bar{p}}{\phantom{i}\text{d}\tau^{*}}$
$\displaystyle=\varepsilon^{*}\nu\hat{c}.$ (95b)

Since $e_{0}=s_{0}$, setting $\varepsilon^{*}=0$ implies that $K_{M}=0$, and
we recover

$\displaystyle\cfrac{\text{d}\hat{c}}{\phantom{i}\text{d}\tau^{*}}$
$\displaystyle=(1-\hat{c})(1-\hat{c}-\bar{p}),$ (96a)
$\displaystyle\cfrac{\text{d}\bar{p}}{\phantom{i}\text{d}\tau^{*}}$
$\displaystyle=0.$ (96b)

Treating $\bar{p}$ as a slowly-varying parameter in the fast subsystem given
by (96a)–(96b), the normal form of the transcritical bifurcation is recovered
by making the change of variables $u=1-\hat{c}$

$\frac{\text{d}u}{\phantom{i}\text{d}\tau^{*}}=\bar{p}u-u^{2}.$ (97)

The normal form equation (97) indicates that the critical
submanifolds141414The transcritical bifurcation is also easily recoverable in
$(s,p)$ coordinates. undergo an exchange of stability at the transcritical
singularity. Clearly, this transcritical singularity influences the dynamics
whenever $K_{M}=0$ and $e_{0}\leq s_{0}$ (see Figure 8).

Figure 8: The critical set of the rQSSA when $e_{0}<s_{0}$ and $K_{M}=0$. This
figure provides a visualization of the invariant set recovered by setting
$e_{0}<s_{0}$ and $k_{-1}=k_{2}=0$. Left: The critical set contains two
transversely intersecting branches of fixed points. Thick solid curves
correspond to non-isolated stable fixed points, and thick dashed lines
correspond to unstable fixed points. The horizontal curve corresponds to the
critical set $\hat{c}=1$, and the diagonal curve corresponds to the critical
set $\bar{p}+\hat{c}=1$. The trajectory rapidly approaches the curve
$\hat{c}=1$, then reaches the curve $1=\hat{c}+\bar{p}$ once $t\sim t_{\ell}$,
and begins descending towards the origin. Clearly, each set constitutes a
normally hyperbolic manifold everywhere except where the branches intersect,
which corresponds to a transcritical singularity. Right: A typical trajectory
(red solid curve) closely follows the attracting critical submanifolds once
the perturbation is turned on.

###  5.2 A heuristic analysis based on timescale separation

The bifurcation analysis from the previous section gave us a global picture of
what the phase-plane dynamics look like when $K_{M}=0$. Moreover, we can use
the bifurcation structure to our advantage and attempt to determine a suitable
“small parameter” that warrants the validity of the rQSSA. To begin, if
$e_{0}<s_{0}$, then there are essentially three stages to the reaction (see
figure 9).

Figure 9: There are three stages in the dynamics of the MM reaction mechanism
(1) when $e_{0}<s_{0}$ and $K_{M}\approx 0$. The red arrows demarcate the
approximate path of a typical trajectory during each stage. Stage one is the
fast transient: the trajectory rapidly approaches the curve $c=e_{0}$, and the
duration of this stage is roughly $t_{C}^{*}$. During the second stage, the
trajectory closely follows the curve $c=e_{0}$; the duration of this stage is
roughly $t_{\ell}$. The final stage corresponds to the the rQSSA, as the
trajectory follows the curve $c=s_{0}-p$. The timescale corresponding to stage
three is $t_{\text{slow}}=1/k_{2}$. If $t_{C}^{*}\ll t_{\ell}\ll 1/k_{2}$,
then the rQSSA is at least heuristically valid.

The first stage is the transient stage, and the phase-plane trajectory moves
almost vertically towards the curve $c=e_{0}$. The approximate duration of
this stage is given by Tzafriri’s fast timescale, $t_{C}^{*}$. After the
initial fast transient the differential equation for product formation is
roughly

$\dot{p}\sim k_{2}e_{0}.$ (98)

Integrating the above expression yields $p\sim k_{2}e_{0}t$. Since $c\approx
e_{0}$ in this regime, the trajectory will reach the vicinity of the $c$-axis
once $p=s_{0}-e_{0}$, which, based on the substrate conservation law, would
imply that $s\approx 0$. Thus, the time it takes the trajectory to reach the
curve $s_{0}=c+p$ is given by $t_{\ell}$:

$t_{\ell}=\cfrac{s_{0}-e_{0}}{k_{2}e_{0}}.$ (99)

Next, the trajectory can be assumed to be very close to the curve $c=s_{0}-p$
and, at this point, the time it takes for the reaction to effectively complete
is $t_{\text{slow}}=1/k_{2}$. If the ratio of $t_{\ell}$ to $t_{\text{slow}}$
is small, i.e.,

$0<\cfrac{t_{\ell}}{t_{\text{slow}}}=\cfrac{s_{0}-e_{0}}{e_{0}}\ll 1,$ (100)

then, from a practical point of view (as opposed to a rigorous point of view),
the rQSSA can be considered valid if $s_{0}>e_{0}$ and $K_{M}$ is incredibly
small.

The heuristic analysis at least partially supports the idea that parameter
domain over which the rQSSA is valid may in fact be larger than previously
assumed. However, the combined conditions, $\varepsilon^{*}\ll 1$ and
$(s_{0}-e_{0})/e_{0}\ll 1$, still do not explain how increasing $s_{0}$
impacts the validity of the approximation. This is a reasonable question that
has plagued theorists for some time now (see, for example [30] and [31]). The
usual methods of scaling analysis and timescale estimation have not provided
very satisfying answers. Consequently, this suggests that a new method needs
to be employed to tackle this problem.

###  5.3 The reverse quasi-steady-state approximation: energy methods

In order to determine a small parameter that ensures the validity of the
rQSSA, it helps to first consider what we expect to happen chemically when
such a parameter is made arbitrarily small. Historically, thanks to the
careful work of Nguyen and Fraser [40], the geometric interpretation of the
rQSSA’s validity has been associated with a phase-plane trajectory closely
following the $s$-nullcline in the $(s,c)$ phase–plane. However, instead of
thinking about the validity of the rQSSA in terms of a trajectory’s propensity
for following the $s$-nullcline, it is perhaps better to think of it in purely
chemical terms: when the rQSSA is valid, the initial substrate concentration
vanishes rapidly with respect to the timescale $\tilde{T}\equiv k_{2}t$.
Translated mathematically, this means that $s$ must dissipate rapidly with
respect to $\tilde{T}$. Thus, we want an equation that describes the
dissipation of $s^{2}/2$:

$\displaystyle\cfrac{1}{2}\cfrac{\text{d}s^{2}}{\text{d}t}$
$\displaystyle=-k_{1}(e_{0}-c)s^{2}+k_{-1}cs,$ (101a) $\displaystyle\leq-
k_{1}(e_{0}-\lambda)s^{2}+k_{-1}s\lambda.$ (101b)

The term $k_{-1}s\lambda$ can be expanded thanks to Cauchy’s inequality

$k_{-1}s\lambda\leq\delta s^{2}+\cfrac{k_{-1}^{2}\lambda^{2}}{4\delta},$ (102)

and choosing $\delta=k_{1}(e_{0}-\lambda)/2$ yields

$\cfrac{1}{2}\cfrac{\text{d}s^{2}}{\text{d}t}\leq-\cfrac{1}{2}k_{1}(e_{0}-c)s^{2}+\cfrac{k_{-1}^{2}\lambda^{2}}{2k_{1}(e_{0}-\lambda)}.$
(103)

From Gronwall’s lemma we have151515From this point forward, we will drop
exponential terms of the form “$(1-e^{-\lambda t})$” that converge to $1$ as
$t\to\infty$ and just replace it by “$1$.”

$\displaystyle s^{2}(t)$ $\displaystyle\leq s_{0}^{2}e^{\displaystyle-
k_{1}(e_{0}-\lambda)t}+\cfrac{k_{-1}^{2}\lambda^{2}}{k_{1}^{2}(e_{0}-\lambda)^{2}}\bigg{(}1-e^{\displaystyle-
k_{1}(e_{0}-\lambda)t}\bigg{)},$ (104a) $\displaystyle\leq
s_{0}^{2}e^{\displaystyle-
k_{1}(e_{0}-\lambda)t}+\cfrac{k_{-1}^{2}\lambda^{2}}{k_{1}^{2}(e_{0}-\lambda)^{2}}.$
(104b)

Diving both sides by $s_{0}^{2}$, taking the square root of both sides, and
invoking the triangle inequality yields:

$s(t)\leq s_{0}e^{\displaystyle-
k_{1}(e_{0}-\lambda)t/2}+\cfrac{k_{-1}\lambda}{k_{1}(e_{0}-\lambda)}.$ (105)

Next, re-writing (105) in terms of $\tilde{T}$,

$\displaystyle\bar{s}(\tilde{T})$ $\displaystyle\leq
e^{\displaystyle-\tilde{T}/2\underline{\varepsilon}\nu}+\tilde{\varepsilon}\cfrac{\lambda}{s_{0}},$
(106a) $\displaystyle\leq
e^{\displaystyle-\tilde{T}/2\underline{\varepsilon}\nu}+\tilde{\varepsilon},$
(106b)

reveals two small parameters:

$\underline{\varepsilon}\equiv\cfrac{K_{M}}{e_{0}-\lambda},\qquad\tilde{\varepsilon}\equiv\cfrac{K_{S}}{e_{0}-\lambda}.$
(107)

Now we have a choice to make: which small parameter corresponds to the
potential validity of the rQSSA? The answer to this question resides in
Fenichel theory: Setting $\underline{\varepsilon}$ to zero results in a
critical manifold of fixed points, whereas making $\tilde{\varepsilon}=0$ by
setting $k_{-1}=0$ results in the manifold $s_{0}=c+p$ being invariant, but it
is not a manifold of equilibrium points, since $k_{-1}$ is not a TFP (see [58]
for details). Moreover, it holds that

$\bar{s}(\tilde{T})\leq
e^{\displaystyle-\tilde{T}/2\underline{\varepsilon}}+\underline{\varepsilon},$
(108)

since $\nu\leq 1$ and $\tilde{\varepsilon}\leq\underline{\varepsilon}$.
Consequently, we expect $s$ to vanish rapidly in $\tilde{T}$ whenever
$\underline{\varepsilon}\ll 1.$ Furthermore, note that

$\lim_{K_{M}\to
0}\underline{\varepsilon}=0\;\;\text{if}\;\;s_{0}<e_{0},\qquad\lim_{K_{M}\to
0}\underline{\varepsilon}=\cfrac{s_{0}-e_{0}}{e_{0}}\;\;\text{if}\;\;s_{0}>e_{0},$
(109)

and we recover the exact small parameter (100) that was obtained heuristically
in Section 5.1.

As a final check, we observe the influence of $\underline{\varepsilon}$ on the
dynamics in the $(p,c)$ phase-plane when $0<\underline{\varepsilon}\ll 1$.
First note that the $c$-nullcline can be approximated by the curve $s_{0}=c+p$
as $K_{M}\to 0$ in any region where $s_{0}<e_{0}$. Thus, if trajectories
closely follow the $c$-nullcline when $K_{M}\ll 1$ and $s_{0}<e_{0}$, then the
rQSSA should be valid. Second, one can show that

$\displaystyle|c-c_{0}(p)|\leq\lambda
e^{\displaystyle-t\zeta_{T}/2}+\cfrac{\eta\lambda}{(1+\eta)(1+\kappa)}\bigg{(}\cfrac{K_{M}}{e_{0}-\lambda+K_{M}}\bigg{)},$
(110)

where $\zeta_{T}\equiv k_{1}(e_{0}-\lambda+K_{M})$ (we will derive this in the
Section 6 using energy methods, but also in Section A using a method pioneered
by Gradshtein [14]). Dividing both sides by $\lambda$ and expressing (110) in
terms of $\tilde{T}$ yields

$\displaystyle|\bar{c}-\bar{c}_{0}(\bar{p})|\leq
e^{\displaystyle-\tilde{T}/2\underline{\varepsilon}\nu}+\cfrac{\eta\nu\underline{\varepsilon}}{(1+\eta)(1+\underline{\varepsilon})},$
(111)

from which it follows that

$\displaystyle|\bar{c}-\bar{c}_{0}(\bar{p})|\leq
e^{\displaystyle-\tilde{T}/2\underline{\varepsilon}}+\underline{\varepsilon}.$
(112)

The upper bound (112) reveals that the phase-plane trajectory will closely
follow the $c$-nullcline after a brief transient. Since the $c$-nullcline is
approximately $s_{0}=c+p$ whenever $s_{0}<e_{0}$ and
$\underline{\varepsilon}\ll 1$, it holds that the rQSSA is valid under these
conditions, and we take $\underline{\varepsilon}$ to be the small parameter
that ensures its validity.

In summary of this section, we have shown that, while scaling methods or
Tikhonov-Fenichel parameters can be used to predict the presence of a singular
point, more analysis is needed to assess the validity of the rQSSA in
neighborhoods impacted by the critical singularity. Hence, the presence of the
singular point in the critical set is what makes the analysis of the rQSSA
challenging. As we have shown, energy methods can be easily employed to assess
the validity of the rQSSA when initial substrate and enzyme concentrations are
of similar magnitude. Moreover, the bound

$\displaystyle|c-c_{0}(p)|\leq\lambda
e^{\displaystyle-t\zeta_{T}/2}+\cfrac{\eta\lambda}{(1+\eta)(1+\kappa)}\bigg{(}\cfrac{K_{M}}{e_{0}-\lambda+K_{M}}\bigg{)},$
(113)

provides direct insight into the origin of each critical manifold, and
suggests that there are three parameters, $\eta$, $\underline{\varepsilon}$
and $\nu$, that, when adequately small, guarantee the long-time validity of
the QSSA, the rQSSA and the tQSSA, respectively. We will discuss this
observation in more detail in section 6.

##  6 Final remarks on timescales and small parameters

In this section we explore the relationship between the fast and slow
timescales of the tQSSA, Tikhonov-Fenichel parameters, and the small
parameters obtained from the energy analysis.

###  6.1 Timescale separation and small parameters: local versus global
conditions for the accuracy of QSS approximations

As mentioned earlier, because the tQSSA encompasses the rQSSA and the sQSSA,
the accepted criterion for the validity of a QSSA is separation of Tzafriri’s
timescales, $\varepsilon_{T}\ll 1$ (see 35). We can express the bound on the
$\limsup$ of $\underline{c}$ and $\overline{c}_{0}(\bar{p})$ in terms of
$T_{z}=t/t_{P}$,

$|\underline{c}-\bar{c}_{0}(\bar{p})|\leq e^{\displaystyle-
T_{z}/2\varepsilon_{D}}+\varepsilon_{L}.$ (114)

To compute (114), we start by computing,

$\displaystyle\limsup_{t\to\infty}(c-h^{-}(p))^{2},$ (115)

where we will use the notation

$\displaystyle h^{-}(p)$
$\displaystyle=\cfrac{1}{2}(e_{0}+s_{0}+K_{M}-p)-\cfrac{1}{2}\sqrt{(e_{0}+s_{0}+K_{M}-p)^{2}-4e_{0}(s_{0}-p)},$
(116a) $\displaystyle h^{+}(p)$
$\displaystyle=\cfrac{1}{2}(e_{0}+s_{0}+K_{M}-p)+\cfrac{1}{2}\sqrt{(e_{0}+s_{0}+K_{M}-p)^{2}-4e_{0}(s_{0}-p)},$
(116b)

for convenience. Thus, the reader should identify $\lambda=h^{-}(0)$ and note
that $c_{0}(p)\equiv h^{-}(p)$. The differential equation for the energy,
$\mathcal{E}^{2}\equiv(c-h^{-}(p))^{2}$, is

$\cfrac{1}{2}\cfrac{\text{d}\mathcal{E}^{2}}{\text{d}t}=\bigg{[}\cfrac{\text{d}c}{\text{d}t}-\cfrac{\text{d}h^{-}}{\text{d}p}\cfrac{\text{d}p}{\text{d}t}\bigg{]}\mathcal{E}.$
(117)

Carefully note that the derivative for $c$ can be factored,
$\dot{c}=k_{1}(c-h^{-}(p))(c-h^{+}(p))$, and thus

$\displaystyle\cfrac{1}{2}\cfrac{\text{d}\mathcal{E}^{2}}{\text{d}t}$
$\displaystyle=\bigg{[}k_{1}(c-h^{-}(p))(c-h^{+}(p))-\cfrac{\text{d}h^{-}}{\text{d}p}\cfrac{\text{d}p}{\text{d}t}\bigg{]}\mathcal{E}$
(118a)
$\displaystyle=\bigg{[}k_{1}\mathcal{E}(c-h^{+}(p))-\cfrac{\text{d}h^{-}}{\text{d}p}\cfrac{\text{d}p}{\text{d}t}\bigg{]}\mathcal{E}$
(118b)
$\displaystyle\leq\mathcal{E}^{2}k_{1}(c-h^{+}(p))+\max\bigg{|}\cfrac{\text{d}h^{-}}{\text{d}p}\bigg{|}\sup\bigg{|}\cfrac{\text{d}p}{\text{d}t}\bigg{|}|\mathcal{E}|.$
(118c)

Next, the term “$c-h^{+}(p)$” can be bounded above:

$\max c-h^{+}(p)\leq\lambda-\min h^{+}(p)=-k_{1}(e_{0}+K_{M}-\lambda).$ (119)

Using the Cauchy-$\delta$ inequality with
$\delta=k_{1}(e_{0}+K_{M}-\lambda)/2$, we obtain

$\max\bigg{|}\cfrac{\text{d}h^{-}}{\text{d}p}\bigg{|}\sup\bigg{|}\cfrac{\text{d}p}{\text{d}t}\bigg{|}|\mathcal{E}|\leq\cfrac{k_{1}(e_{0}+K_{M}-\lambda)}{2}\mathcal{E}^{2}+\cfrac{\bigg{(}\max\bigg{|}\cfrac{\text{d}h^{-}}{\text{d}p}\bigg{|}\bigg{)}^{2}\bigg{(}\sup\bigg{|}\cfrac{\text{d}p}{\text{d}t}\bigg{|}\bigg{)}^{2}}{2k_{1}(e_{0}+K_{M}-\lambda)}.$
(120)

Inserting (120) into (118c) yields

$\cfrac{\text{d}\mathcal{E}^{2}}{\text{d}t}\leq-\mathcal{E}^{2}k_{1}(e_{0}+K_{M}-\lambda)+\cfrac{\bigg{(}\max\bigg{|}\cfrac{\text{d}h^{-}}{\text{d}p}\bigg{|}\bigg{)}^{2}\bigg{(}\sup\bigg{|}\cfrac{\text{d}p}{\text{d}t}\bigg{|}\bigg{)}^{2}}{k_{1}(e_{0}+K_{M}-\lambda)},$
(121)

and from Gronwall’s Lemma and the triangle inequality we recover:

$|\mathcal{E}(t)|\leq|\mathcal{E}(0)|e^{-\displaystyle
k_{1}(e_{0}+K_{M}-\lambda)t}+\cfrac{\max\bigg{|}\cfrac{\text{d}h^{-}}{\text{d}p}\bigg{|}\sup\bigg{|}\cfrac{\text{d}p}{\text{d}t}\bigg{|}}{k_{1}(e_{0}+K_{M}-\lambda)}.$
(122)

Calculating the “$\max$” and “$\sup$” on the right hand side of (122),
dividing both sides by $\lambda$, and expressing the exponential in terms of
$T_{z}=t/t_{P}$ yields (114), where $\varepsilon_{D}$ denotes the “exponential
$\varepsilon$,” and $\varepsilon_{L}$ denotes the “long-time $\varepsilon$:”

$\varepsilon_{D}\equiv\bigg{(}\cfrac{\lambda}{s_{0}}\bigg{)}\bigg{(}\cfrac{1}{1+\kappa}\bigg{)}\bigg{(}\cfrac{\underline{\varepsilon}}{1+\underline{\varepsilon}}\bigg{)},\qquad\varepsilon_{L}\equiv\bigg{(}\cfrac{\eta}{1+\eta}\bigg{)}\bigg{(}\cfrac{1}{1+\kappa}\bigg{)}\bigg{(}\cfrac{\underline{\varepsilon}}{1+\underline{\varepsilon}}\bigg{)}.$
(123)

As expected, $\varepsilon_{D}$ is dependent on initial data, meaning its
magnitude is influenced by the starting location of the trajectory within the
vector field. Moreover, $\varepsilon_{D}$ can be “factored” into three small
parameters:

  1. (i)

$\lambda/s_{0}$: vanishes as a result of zero enzyme to substrate ratio.

  2. (ii)

$1/(1+\kappa)\equiv\nu$: vanishes in the limit of infinitely slow product
formation.

  3. (iii)

$\underline{\varepsilon}/(1+\underline{\varepsilon})$: vanishes at infinitely
high enzyme concentration or, zero $K_{M}$ whenever $s_{0}\leq e_{0}$.

What is the relationship between Tzafriri’s timescales and $\varepsilon_{D}$
and $\varepsilon_{L}$? Tzafriri’s fast timescale, $t_{C}^{*}$, can be written
as

$t_{C}^{*}=\cfrac{1}{k_{1}(K_{M}+e_{0}-\lambda+s_{0}-\lambda)}.$ (124)

Using the identity (124), a simple calculation reveals
$\varepsilon_{T}\leq\varepsilon_{D}$. In fact,
$\varepsilon_{T}\sim\varepsilon_{D}$ at very high enzyme concentration, and
$\varepsilon_{T}\to 0$ if $\varepsilon_{D}\to 0.$ Thus, we can control
timescale separation by controlling $\varepsilon_{D}$. However,
$\varepsilon_{D}\ll 1$ only implies that the phase-plane trajectory moves
rapidly towards the $c$-nullcline; it does not ensure that the trajectory
stays close to the $c$-nullcline as $T_{z}\to\infty$. The long-time validity
of the zeroth-order approximation is regulated by $\varepsilon_{L}$, and
$\varepsilon_{D}\leq\varepsilon_{L}$. Hence, $\varepsilon_{L}$ is the
fundamental small parameter, since it is independent of a particular choice of
initial data and provides an upper bound on the $\limsup$:

$\displaystyle\limsup_{t\to\infty}(c-c_{0}(p))^{2}\leq\lambda^{2}\varepsilon_{L}^{2}.$
(125)

For example, $\varepsilon_{D}$ can be made arbitrarily small by taking $s_{0}$
to be arbitrarily large, but the accuracy of the approximation in regimes were
$s_{0}$ is large may only be temporary at best (see figure 10). Finally, the
magnitude of $\varepsilon_{L}$ is regulated by the TFPs. The parameter
$\varepsilon_{L}$ determines “how small” a TFP must be in comparison to other
parameters in order to ensure the validity of the tQSSA. Thus, there is a
fundamental “hierarchy” or “ordering” of parameters:

$\text{TFP}\prec\varepsilon_{L}\prec\varepsilon_{D}\prec\varepsilon_{T}.$
(126)

From the perspective of the ordering given in (126), timescale separation
should be seen as an effect, not a cause. Segel and Slemrod [30] understood
this; they were very clear in issuing their warning of sufficient versus
necessary conditions for the validity of QSS reductions, and stressed the fact
that timescale separation is merely a necessary condition.

The takeaway from our analysis is that timescale separation is really a local
metric when it comes to assessing the validity of a QSS reduction. This is
primarily because timescales are dependent on initial conditions, and the
vector field can change dramatically once the trajectory leaves a small
neighborhood containing the initial conditions.

Figure 10: The tQSSA can worsen as $t\to\infty$ if $\varepsilon_{LT}$ is not
sufficiently small. In both panels, the numerical solution is obtained with
parameter values: $e_{0}=10,s_{0}=1000,k_{1}=20,k_{-1}=10,k_{2}=10$. Thus,
$\varepsilon_{T}\sim 10^{-6}$, but $\varepsilon_{LT}\approx 0.1$. top panel:
The $\log_{10}$ of the relative error between the numerical solution to the
mass action equations (27a)-(27b) for $c$, and the numerical solution to the
asymptotic approximation (29), denotes $c^{\varepsilon_{T}}$, is demarcated by
the dashed-dotted line. bottom Panel: The dashed-dotted line demarcates the
$\log_{10}$ of the relative error between the numerical solution to the mass
action equations (27a)-(27b) for $p$, and the numerical solution to the
asymptotic approximation (30), denoted $p^{\varepsilon_{T}}$. Even though
timescale separation is quite large, there are regions in the numerically-
obtained time course data where $c^{\varepsilon_{T}}$ approximates $c$ only to
within one digit of accuracy. While the loss in accuracy may not be outright
disastrous in terms of parameter estimation, it is an observation that is
worth pointing out in order to gain a deeper quantitative understanding of the
various QSS approximations utilized for Michaelis–Menten reaction mechanism
(1). The parameters have been assigned arbitrary units for illustrative
purposes.

Roussel and Fraser [41] seem to have been the first to really understand this
in the context of the Michaelis–Menten reaction mechanism (1), as their work
focused on approximating the actual invariant slow manifold. In fact, to the
best of the authors’ knowledge, Roussel and Fraser (see [61], pp. 42-47) were
the first to provide a solid geometrical argument that $\eta$ is in fact the
“global” qualifier for the validity of the sQSSA instead of
$\varepsilon_{SS}$. By approximating the slow manifold to a high degree of
accuracy, Roussel and Fraser [62, 42] were able to demonstrate that $\eta\ll
1$ is necessary for the long-time validity of the sQSSA.

###  6.2 Tightening the limit supremum

It is worth pointing out that a long-time bound smaller than $\varepsilon_{L}$
can be obtained when $e_{0}<s_{0}$. This is to be expected, since
$\underline{\varepsilon}$ does not vanish if $s_{0}>e_{0}$ as $K_{M}\to 0$. To
obtain a tighter bound, note that the upper bound of $c-h^{+}(p)$ given in
(119) is rather liberal. What we want to do is compute $\max(c-h^{+}(p))$ with
more care. First, since $c<h^{+}(p)$ and $0<h^{+}(p)$, it follows that
$c-h^{+}(p)<0$. Second, $h^{+}(p)$ is a decreasing function of $p$, which
means $\min h^{+}(p)$ occurs when $p=s_{0}$. However, if $c\neq 0$, then
$p\leq s_{0}-c$. Thus, we want to maximize:

$c-h^{+}(s_{0}-c)\equiv\theta(c)<0.$ (127)

Since $\theta^{\prime}(c)>0$ and $\theta(c)<0$, the maximum of $\theta(c)$
occurs when $c=\lambda$, from which it is straightforward to show

$\cfrac{\text{d}\mathcal{E}^{2}}{\text{d}t}\leq-\mathcal{E}^{2}k_{1}|\theta(\lambda)|+\cfrac{\bigg{(}\max\bigg{|}\cfrac{\text{d}h^{-}}{\text{d}p}\bigg{|}\bigg{)}^{2}\bigg{(}\sup\bigg{|}\cfrac{\text{d}p}{\text{d}t}\bigg{|}\bigg{)}^{2}}{k_{1}|\theta(\lambda)|},$
(128)

where $\theta(\lambda)$ is given by:

$\theta(\lambda)=\cfrac{\lambda}{2}-\cfrac{1}{2}\bigg{(}e_{0}+K_{M}+\sqrt{(e_{0}+K_{M}+\lambda)^{2}-4e_{0}\lambda}\bigg{)}.$
(129)

An even tighter bound on $\theta(c)$ can be obtained by noting that
$\lambda\leq e_{0}$, and thus $\theta(\lambda)\leq\theta(e_{0})$:

$\theta(c)\leq\theta(e_{0})=-\cfrac{1}{2}\bigg{(}K_{M}+\sqrt{4e_{0}K_{M}+K_{M}^{2}}\bigg{)}.$
(130)

The application of Gronwall’s Lemma reveals

$\limsup_{t\to\infty}(c-h^{-}(p))^{2}\leq\lambda^{2}\bigg{(}\cfrac{\eta}{1+\eta}\bigg{)}^{2}\bigg{(}\cfrac{1}{1+\kappa}\bigg{)}^{2}\bigg{(}\cfrac{2K_{M}}{K_{M}+\sqrt{K_{M}^{2}+4e_{0}K_{M}}}\bigg{)}^{2}\equiv\lambda^{2}\varepsilon_{LT}^{2},$
(131)

which vanishes as $K_{M}\to 0$. Thus, the long-time bound “$\varepsilon_{LT}$”
vanishes if $e_{0}\to 0$, $K_{M}\to 0$, $e_{0}\to\infty$, $K_{M}\to\infty$, or
$k_{2}\to 0$. Finally, note that

$\cfrac{K_{M}}{e_{0}-\lambda+K_{M}}=\cfrac{K_{M}}{-\theta(e_{0})}=\cfrac{K_{M}}{|\theta(e_{0})|}$
(132)

when $e_{0}=s_{0}$, so there is a nice “overlap” of small parameters. Note
also that the rate of convergence as $K_{M}\to 0$ is
$\mathcal{O}(\sqrt{K_{M}})$ when $e_{0}\leq s_{0}$. This is as expected. In
fast/slow systems that contain a transcritical point, trajectories remain at a
distance that is $\mathcal{O}(\sqrt{\varepsilon})$ from the transcritical
point once $0<\varepsilon\ll 1$ (again, see [53] for details).

##  7 Discussion

The primary goal of this work was to illustrate how a combination of scaling
and non-scaling methods can be employed to determine the validity of the
sQSSA, extended QSSA, rQSSA, and tQSSA. Collectively, our work provides an
analysis that admits a considerably clearer understanding of each QSSA, its
corresponding origin, and its validity. Using energy methods, we have
recovered a clear set of associated “small parameters” that correspond to the
three types of fast/slow dynamics considered under the QSSA: the sQSSA, the
rQSSA, and the tQSSA (see Table 1).

Parameter | Approximation | Critical Manifold (Set)
---|---|---
$\eta$ | sQSSA | $1-\bar{s}-\bar{p}=0$
$\underline{\varepsilon}$ | rQSSA | $\\{(\hat{c},\bar{p})\in\mathbb{R}^{+}|\;(1-\hat{c})(1-\hat{c}-\bar{p})=0\\}$
$\nu$ | tQSSA | $\bar{c}_{0}(\bar{p};K_{S})$
Table 1: Small parameters and critical sets for the standard, reverse, and
total QSSA. Each of the critical manifolds that appear when either
$\eta,\underline{\varepsilon}$ or $\nu$ are identically zero correspond to
zero product formation. Hence, the tQSSA is valid in all three cases. The
notation $\bar{c}_{0}(\bar{p};K_{S})$ denotes the $c$-nullcline with $K_{M}$
replaced with $K_{S}$, which occurs when $k_{2}=0.$

The closure of the rQSSA problem for the MM reaction mechanism (1) is a
testament to the power of using a combination of scaling and non-scaling
methods, as we were able to derive a “small parameter” that ensures the
validity of the rQSSA, $\underline{\varepsilon}\equiv K_{M}/(e_{0}-\lambda)$.
In previous studies [30, 31], it was assumed that the rQSSA is valid only when
$s_{0}\ll e_{0}$. The small parameter, $\underline{\varepsilon}$, is a major
improvement on the small parameters previously reported in the literature,
because it quantifies and illustrates the validity of the rQSSA in regimes
where $e_{0}\lesssim s_{0}$. Moreover, for the first time, we have shown that
the critical manifold associated with the rQSSA at finite initial enzyme
concentration contains a transcritical singularity that significantly
influences the validity of the rQSSA whenever $e_{0}\leq s_{0}$ and $K_{M}\ll
e_{0}$, since the critical (sub) manifolds undergo an exchange of stability at
the singular point. On the other hand, this singularity is removed from the
physical domain of interest if one takes $e_{0}\to\infty$. This is the
distinguishing difference between our analysis and the previous analysis of
Schnell and Maini [31].

Our analysis has practical implications in the laboratory for the use of
mathematical approximations under the rQSSA. Experimental assays are typically
implemented to quantify the enzyme activity of a specific reaction. This is
done by fitting experimental data to a particular model equation derived from
a QSS reduction (i.e., either the sQSSA, the rQSSA, or the tQSSA). When the
rQSSA is valid, progress curve experiments for $p$ can be fitted to the
product formation rate expression (23), which will allow the estimation of the
catalytic rate constant, $k_{2}$. In mathematical terms, the procedure of
fitting experimental data to a mathematical model so as to estimate kinetic
parameters is an inverse problem [15, 16]. Contrary to what has been reported,
experimental assays can be carried out when initial substrate and enzyme
ratios are of similar magnitude, and the rQSSA can be used to estimate
$k_{2}$, as long as $\underline{\varepsilon}\ll 1.$

The mathematical analysis of the Michaelis-Menten reaction mechanism (1) has a
long and rich history. We believe that our analysis has deepened the overall
understanding of the QSSA, and that our approach provides an important
improvement to the implementation of this approximation in enzyme kinetics.
The methods we employed can be adapted to analyze complex enzyme catalyzed
reactions, or multiscale dynamical systems. This will be the focus of our
future work.

## Acknowledgement

We are grateful to Sebastian Walcher (RWTH Aachen University) and Wylie
Stroberg (University of Michigan Medical School) for reading and providing
comments to previous versions of this manuscript. Justin Eilertsen is
supported by the University of Michigan Postdoctoral Pediatric Endocrinology
and Diabetes Training Program “Developmental Origins of Metabolic Disorder”
(NIH/NIDDK Grant: K12 DK071212).

##  Appendix A Appendix

Here we present a detailed outline of the the method of slowly-varying
Lyapunov functions applied to the governing equations of the Michaelis–Menten
reaction mechanism as discussed in Section 4.3

###  A.1 Gronwall’s lemma

This differential version can be found in [63]. Suppose that $u(t)\in
C^{1}([0,T];\mathbb{R}_{0}^{+})$ satisfies161616The notation
$L^{1}(\Omega;\mathbb{R}_{0}^{+})$ refers to the space of functions that are
Lebesgue-integrable on $\Omega$ that map $\Omega$ to $\mathbb{R}_{0}^{+}$.

$\cfrac{\text{d}u}{\text{d}t}\leq\lambda u+f(t),\qquad u(0)=u_{0}.$ (133)

Then, for $\lambda\in\mathbb{R}$ and $f(t)\in
L^{1}([0,T];\mathbb{R}_{0}^{+})$, it holds that:

$u(t)\leq\displaystyle e^{\displaystyle\lambda t}u_{0}+e^{\displaystyle\lambda
t}\int_{0}^{t}e^{\displaystyle-\lambda\tau}f(\tau)\text{d}\tau,\quad\forall
t\in[0,T].$ (134)

###  A.2 The “contractivity” condition: Existence of a Lyapunov function

It will help to briefly review the concept of a Lyapunov function before we
introduce the idea of a “contractivity” condition. Consider a simple dynamical
system of the form

$x^{\prime}=f(x),$ (135)

where “′” denotes differentiation with respect to time, $t$. The vector field
$f:\mathbb{R}^{n}\mapsto\mathbb{R}^{n}$ generates a flow, and any scalar-
valued function of $x$ (we will call this function
$L(x):\mathbb{R}^{n}\mapsto\mathbb{R}$) will change with respect to time
according to:

$L^{\prime}=f(x)\cdot\nabla L\equiv\mathcal{G}L,$ (136)

where the operator “$\mathcal{G}$” is known as a “generator.” Now suppose that
the dynamical system (135) has an equilibrium point at $x=x^{*}$, so that
$f(x^{*})=0$. Linearization is typically the preferred method of choice to
determine the stability of $x^{*}$. However, it is well-known that if

$0<L,\;\;\&\;\;\mathcal{G}L<0,\quad\forall
x\in\mathcal{B}\backslash\\{x^{*}\\},$ (137)

where $x^{*}$ is an interior point of the neighborhood (ball) $\mathcal{B}$,
then the equilibrium point $x^{*}$ is asymptotically stable, and “$L$” is
called a Lyapunov function.

So what does this have to do with a contraction condition? Keeping the idea of
a Lyapunov function in mind, consider a fast/slow system:

$\displaystyle x^{\prime}$ $\displaystyle=\varepsilon f(x,y),\quad x(0)=x_{0}$
(138a) $\displaystyle y^{\prime}$ $\displaystyle=g(x,y),\quad y(0)=y_{0},$
(138b)

where the slow variable, $x\in\mathbb{R}^{n}$, the fast variable
$y\in\mathbb{R}^{m}$, and thus
$f(x,y):\mathbb{R}^{n}\times\mathbb{R}^{m}\mapsto\mathbb{R}^{n}$ and
$g(x,y):\mathbb{R}^{n}\times\mathbb{R}^{m}\mapsto\mathbb{R}^{m}$. The “fast
subsystem” associated with (138a)–(138b) is

$y^{\prime}=g(x_{0},y).$ (139)

Assume the fast subsystem has an equilibrium point $y^{*}=h_{0}(x_{0})$, where
$h:\mathbb{R}^{n}\mapsto\mathbb{R}^{m}$, and $g(x_{0},h_{0}(x_{0}))=0$. Since
$x_{0}$ is a constant with respect to the fast subsystem, $g$ can be treated
as a map from $\mathbb{R}^{m}$ to $\mathbb{R}^{m}$. Once again, the stability
of $y^{*}$ can be found via linearization. However, suppose we want to
determine an appropriate Lyapunov function. Generally, finding a Lyapunov
function can be difficult, and so it is common practice to write down a simple
polynomial function and “check” to see if it satisfies the requirements given
in $(\ref{eq:req})$. The simplest choice is a quadratic function,
$L(y)=||y-h_{0}(x_{0})||^{2}$, which is locally positive-definite since
$L(y)>0,\;\;\forall y\neq y^{*}=h_{0}(x_{0})$. The second property we need to
check is whether or not $L^{\prime}$ is locally negative-definite; thus, we
need to show that

$L^{\prime}=2\langle(y-y^{*}),g(x_{0},y)\rangle<0,\quad\forall
y\in\mathcal{B}\backslash\\{y^{*}\\},$ (140)

where $\langle,\rangle$ denotes the usual scalar product between two vectors,
and $||v||=\sqrt{\langle v,v\rangle}$. It is not entirely clear how go about
showing that $L^{\prime}$ is negative definite, since $g(x_{0},y)$ and
$(y-y^{*})$ can change signs. However, a little “trick” can go a long way.
First, rewrite (140) as

$2\langle(y-y^{*}),g(x_{0},y)\rangle=2\langle(y-y^{*}),g(x_{0},y)-g(x_{0},y^{*})\rangle.$
(141)

Now suppose we know that there is some positive number “$\zeta$” such that

$\langle(y-y^{*}),g(x_{0},y)-g(x_{0},y^{*})\rangle\leq-\zeta||y-y^{*}||^{2}.$
(142)

If such a number “$\zeta$” exists, then we can bound $L^{\prime}$

$L^{\prime}=2\langle(y-y^{*}),g(x_{0},y)\rangle\leq-2\zeta||y-y^{*}||^{2}\leq
0,$ (143)

and thus we have shown that $L^{\prime}$ is locally negative-definite.
Consequently, $L=||y-y^{*}||^{2}\equiv||y-h_{0}(x_{0})||^{2}$ is a Lyapunov
function. Why then, do we refer to (142) as a “contractivity” condition? Well,
notice that (143) defines a differential inequality

$L^{\prime}\leq-2\zeta L,$ (144)

and integrating this inequality (courtesy of Gronwall’s lemma) yields

$L\leq L(0)e^{\displaystyle-2\zeta t}.$ (145)

Consequently, the distance between $y_{0}$ and $h_{0}(x_{0})$ decays to zero
exponentially; hence, the distance “contracts.”

In the subsection that follows, we will investigate the growth of
$L=||y-h_{0}(x)||^{2}$, but we will allow $x$ to change in time, since $x$
will vary slowly in comparison to $y$ when $0<\varepsilon\ll 1$. Thus, we will
be interested in finding an upper bound on

$\cfrac{\text{d}L}{\text{d}T}\equiv\dot{L},$ (146)

where $T=\varepsilon t$ is the slow timescale. The motivation for doing this
should be obvious: we want to compute an estimate on how well the reduced
system, $\dot{x}=f(x,h_{0}(x))$, approximates the full system. Because we will
take into account the temporal variation of $x$, $L$ is no longer a Lyapunov
function. However, because $x$ varies slowly, we have chosen to refer to $L$
in this case as a slowly-varying Lyapunov function. One could also refer to
this procedure as simply an “energy method,” since we are ultimately bounding
the energy (norm) of the error $y-h_{0}(x)$.

###  A.3 General procedure: Estimation of bounds via energy methods

Let us start by estimating bounds on the enslavement of the fast variable. The
calculation outlined below follows from the original work of Tikhonov [13] and
Gradshtein [14], and the particular calculation we outline here can be found
in [63]. This procedure will then be used to develop bounds for enslavement in
the $(p,c)$ coordinate system. We will express the procedure in terms with
respect to a the fast/slow system (138a)–(138b) and, since our specific
problem of interest is two-dimensional, we will assume $x\in\mathbb{R}$ and
$y\in\mathbb{R}$. However, the reader should keep in mind that this procedure
extends to higher-dimensional problems.

With respect to the slow time, the system (138a)–(138b) is expressed in form

$\displaystyle\dot{x}$ $\displaystyle=f(x,y),$ (147a)
$\displaystyle\varepsilon\dot{y}$ $\displaystyle=g(x,y),$ (147b)

where “$\dot{\phantom{x}}$” denotes differentiation with respect to the slow
time, $T$. Let $g(\cdot,y)$ be contractive in the sense that

$(g(x,y_{1})-g(x,y_{2}))\cdot(y_{1}-y_{2})\leq-\zeta\cdot(y_{1}-y_{2})^{2},$
(148)

where the “$\cdot$” in (148) and in what follows denotes scalar
multiplication. Next, suppose that $h_{0}(x)$ satisfies $g(x,h_{0}(x))=0$.
Defining $y=h_{0}(x)+z$, where $z=y-h_{0}(x)$, it follows that

$\dot{z}=\dot{y}-D_{x}h_{0}(x)\cdot\dot{x}.$

Then, the differential equation for the energy, $L\equiv
z^{2}=(y-h_{0}(x))^{2}$, is

$\cfrac{1}{2}\dot{z^{2}}=\cfrac{1}{\varepsilon}z\cdot g(x,z+h_{0}(x))-z\cdot
D_{x}h_{0}(x)\cdot\dot{x}.$ (149)

Next, we want to get this expression into a useful form so that we can apply
Gronwall’s lemma. There are essentially two “tricks” that can be utilized to
manipulate this expression: (i.) we can “add/subtract” zero, or (ii.) we can
“multiply” by one. We will choose (i). The term $\varepsilon^{-1}z\cdot
g(x,z+h_{0}(x))$ can be expanded by “subtracting” zero,

$\cfrac{1}{\varepsilon}z\cdot
g(x,z+h_{0}(x))=\cfrac{1}{\varepsilon}\bigg{[}z\cdot g(x,z+h_{0}(x))-z\cdot
g(x,h_{0}(x))\bigg{]},$ (150)

which holds since $g(x,h_{0}(x))=0$. Thus, we obtain

$\cfrac{1}{2}\dot{z^{2}}=\cfrac{1}{\varepsilon}\bigg{[}z\cdot
g(x,z+h_{0}(x))-z\cdot g(x,h_{0}(x))\bigg{]}-z\cdot
D_{x}h_{0}(x)\cdot\dot{x}.$ (151)

Our primary interest is not in finding an actual solution to (151). Instead,
we want to compute an upper bound on the growth of $z^{2}$, which means we
need to convert the differential equation (151) into a differential
inequality. To do this, we will take advantage of the “contractivity”
condition. First rewrite the first term on the right hand side of (151) as:

$\bigg{[}z\cdot g(x,z+h_{0}(x))-z\cdot
g(x,h_{0}(x))\bigg{]}=\bigg{[}g(x,z+h_{0}(x))-g(x,h_{0}(x))\bigg{]}\cdot\bigg{[}z+h_{0}(x)-h_{0}(x)\bigg{]}.$
(152)

Let $z+h_{0}(x)$ denote “$y_{1}$” in (148), and let $h_{0}(x)$ denote
“$y_{2}$” in (148). Then,

$(g(x,z+h_{0}(x))-g(x,h_{0}(x)))\cdot(z+h_{0}(x)-h_{0}(x))\leq-\zeta\cdot
z^{2}$ (153)

It follows from (153) that

$\cfrac{1}{2}\dot{z^{2}}\leq-\cfrac{\zeta}{\varepsilon}z^{2}-z\cdot
D_{x}h_{0}(x)\cdot\dot{x}.$ (154)

Maximizing the inequality (154) yields

$\cfrac{1}{2}\dot{z^{2}}\leq-\cfrac{\zeta}{\varepsilon}z^{2}+|z|\max|D_{x}h_{0}(x)|\max|\dot{x}|.$
(155)

In order to use Gronwall’s lemma, we want the right hand side of (155) to be
in terms of $z^{2}$ instead of $|z|$. To do this, let $a\equiv|z|$ and
$b\equiv\max|D_{x}h_{0}(x)|\max|\dot{x}|$. Then, from Cauchy’s inequality,
($a\cdot b\leq\delta a^{2}+b^{2}/4\delta$), it holds that

$|z|\max|D_{x}h_{0}(x)|\max|\dot{x}|\leq\delta
z^{2}+\cfrac{(\max|D_{x}h_{0}(x)|\max|\dot{x}|)^{2}}{4\delta},\quad\forall\delta>0.$
(156)

Choosing $\delta=\zeta/2\varepsilon$ yields

$|z|\max|D_{x}h_{0}(x)|\max|\dot{x}|\leq\cfrac{\zeta}{2\varepsilon}z^{2}+\varepsilon\cfrac{(\max|D_{x}h_{0}(x)|\max|\dot{x}|)^{2}}{2\zeta}.$
(157)

Substitution of (157) into (155) reveals the differential inequality

$\dot{z^{2}}\leq-\cfrac{\zeta}{\varepsilon}z^{2}+\varepsilon\cfrac{(\max|D_{x}h_{0}(x)|\max|\dot{x}|)^{2}}{2\zeta}.$
(158)

The expression (158) is now in a form that be integrated directly, and it
follows from Gronwall’s lemma that,

$\displaystyle z^{2}$ $\displaystyle\leq z^{2}(0)\displaystyle
e^{\displaystyle-\zeta
T/\varepsilon}+\varepsilon^{2}\cfrac{(\max|D_{x}h_{0}(x)|\max|\dot{x}|)^{2}}{\zeta^{2}}\bigg{(}1-\displaystyle
e^{\displaystyle-\zeta T/\varepsilon}\bigg{)},$ (159a) $\displaystyle\leq
z^{2}(0)\displaystyle e^{\displaystyle-\zeta
T/\varepsilon}+\varepsilon^{2}\cfrac{(\max|D_{x}h_{0}(x)|\max|\dot{x}|)^{2}}{\zeta^{2}}.$
(159b)

Finally, from the triangle inequality, we obtain:

$|z|\equiv|y-h_{0}(x)|\leq|y_{0}-h_{0}(x_{0})|\displaystyle
e^{\displaystyle-\zeta
T/2\varepsilon}+\varepsilon\cfrac{\max|D_{x}h_{0}(x)|\max|\dot{x}|}{\zeta}.$
(160)

###  A.4 Computation of the upper bound for enslavement of $c$ in $(p,c)$
coordinates

Here we will directly calculate the constants $\zeta,\max|D_{x}h(x)|$ and
$\max|\dot{x}|$ need for the upper bound. Note that the computation of $\zeta$
is calculated over the domain
$(p,c)\in\Omega\equiv[0,s_{0}]\times[0,\lambda]$, since trajectories that
start on $(p,c)(0)=(0,0)$ are confined to this domain due to conservation.

In the $(p,c)$ coordinate system, the $c$-nullcline is given by:

$h(p)=\cfrac{1}{2}\bigg{(}e_{0}+K_{M}+s_{0}-p-\sqrt{(e_{0}+K_{M}+s_{0}-p)^{2}-4e_{0}(s_{0}-p)}\bigg{)}.$
(161)

Finding $\max|h^{\prime}(p)|$ is relatively straightforward. First,

$\cfrac{\text{d}^{2}h}{\text{d}p^{2}}=-\cfrac{2e_{0}K_{M}}{((e_{0}+K_{M}+s_{0}-p)^{2}-4e_{0}(s_{0}-p))^{3/2}}<0,$
(162)

and therefore the derivative is monotonically decreasing on the interval
$p\in[0,s_{0}]$. Thus, we obtain:

$\max\bigg{|}\cfrac{\text{d}h}{\text{d}p}\bigg{|}=-\cfrac{\text{d}h}{\text{d}p}\bigg{|}_{p=s_{0}}=\cfrac{e_{0}}{K_{M}+e_{0}}.$
(163)

For contraction, we define

$G(c,p)\equiv\dot{c}=k_{1}c^{2}-k_{1}(K_{M}+e_{0}+s_{0}-p)c+k_{1}e_{0}(s_{0}-p).$
(164)

We will use the Mean Value Theorem to find $\zeta$. The derivative of $G(c,p)$
is given by

$\cfrac{\partial G}{\partial c}=-k_{1}(e_{0}+s_{0}+K_{M}-2c-p)\leq-
k_{1}(e_{0}+K_{M}+s_{0})+\sup k_{1}(2c+p).$ (165)

We need to find $\sup(2c+p)$. A straightforward bound is $2\lambda+s_{0}$.
However, we want to find the supremum subject to the constraints:
$s_{0}=s+c+p$ and $c\leq\lambda$. Employing the conservation law
“$s_{0}=s+c+p$ yields

$2c+p=s_{0}-s+c\leq s_{0}+\lambda.$ (166)

Consequently, since $2c+p=s_{0}+c-s\leq s_{0}+\lambda$, it follows that

$|c-c_{0}(p)|\leq\lambda
e^{\displaystyle-t\zeta_{T}/2}+\cfrac{\eta\lambda}{(1+\eta)(1+\kappa)}\bigg{(}\cfrac{K_{M}}{e_{0}-\lambda+K_{M}}\bigg{)},$
(167)

where $\zeta_{T}$ is given by

$\zeta_{T}=k_{1}(e_{0}-\lambda+K_{M}).$ (168)

Diving both sides of (167) by $\lambda$ and expressing $p=s_{0}\bar{p}$ yields
the dimensionless form of the inequality

$|\underline{c}-\bar{c_{0}}(\bar{p})|\leq
e^{\displaystyle-\tilde{T}/2\underline{\varepsilon}\nu}+\cfrac{\underline{\varepsilon}}{(1+\underline{\varepsilon})}\cfrac{\eta}{(1+\eta)}\cfrac{1}{(1+\kappa)}.$
(169)

Since $\nu\leq 1$, we obtain

$|\underline{c}-\bar{c_{0}}(\bar{p})|\leq
e^{\displaystyle-\tilde{T}/2\underline{\varepsilon}}+\underline{\varepsilon}.$
(170)

Furthermore, since $\zeta_{T}\geq k_{1}K_{M}$, it holds that

$|\underline{c}-\bar{c_{0}}(\bar{p})|\leq
e^{\displaystyle-\bar{T}/2\varepsilon_{SS}}+\eta.$ (171)

Let us define a “long-time epsilon, $\varepsilon_{L}$,”

$\varepsilon_{L}\equiv\cfrac{\underline{\varepsilon}}{(1+\underline{\varepsilon})}\cfrac{\eta}{(1+\eta)}\cfrac{1}{(1+\kappa)}.$
(172)

Expressing the exponential in terms of Tzafriri’s timescale,
$T_{z}=k_{2}\lambda t/s_{0}$, we obtain:

$|\underline{c}-\bar{c_{0}}(\bar{p})|\leq e^{\displaystyle-
s_{0}T_{z}/2\lambda\underline{\varepsilon}\nu}+\varepsilon_{L}.$ (173)

Moreover, we can define a “decay epsilon,” $\varepsilon_{D}\equiv
s_{0}^{-1}\lambda\underline{\varepsilon}\nu$, and write

$|\underline{c}-\bar{c_{0}}(\bar{p})|\leq e^{\displaystyle-
T_{z}/2\varepsilon_{D}}+\varepsilon_{L},$ (174)

where $\varepsilon_{D}\leq\varepsilon_{L}.$

###  A.5 The tQSSA as used in practice

The $c$-nullcline is rarely used for parameter estimation in progress curve
experiments. Typically, at high enzyme concentrations, one uses

$c\approx\cfrac{e_{0}(s_{0}-p)}{e_{0}+K_{M}+s_{0}-p}$ (175)

in order to employ

$\cfrac{\text{d}p}{\text{d}t}\approx\cfrac{k_{2}e_{0}(s_{0}-p)}{e_{0}+K_{M}+s_{0}-p}$
(176)

for parameter estimation (see [16] for a good discussion on the tQSSA and
parameter estimation). Calculating an upper bound on

$|\mathcal{E}_{T}|\equiv\bigg{|}c-\cfrac{e_{0}(s_{0}-p)}{e_{0}+K_{M}+s_{0}-p}\bigg{|}$
(177)

is done using the methods previously. First,

$\displaystyle\cfrac{1}{2}\cfrac{\text{d}\mathcal{E}_{T}^{2}}{\text{d}t}$
$\displaystyle=-k_{1}(e_{0}+K_{M}+s_{0}-p)\mathcal{E}_{T}^{2}+\bigg{[}k_{1}c^{2}+\cfrac{k_{2}ce_{0}(e_{0}+K_{M})}{(e_{0}+K_{M}+s_{0}-p)^{2}}\bigg{]}\mathcal{E}_{T},$
(178a) $\displaystyle\leq-
k_{1}(e_{0}+K_{M})\mathcal{E}_{T}^{2}+\bigg{[}k_{1}\lambda^{2}+\cfrac{k_{2}\lambda
e_{0}}{e_{0}+K_{M}}\bigg{]}|\mathcal{E}_{T}|.$ (178b)