On the quasi-steady-state approximation in an open Michaelis-Menten reaction mechanism.

The conditions for the validity of the standard quasi-steady-state approximation in the Michaelis-Menten mechanism in a closed reaction vessel have been well studied, but much less so the conditions for the validity of this approximation for the system with substrate inflow. We analyze quasi-steady-state scenarios for the open system attributable to singular perturbations, as well as less restrictive conditions. For both settings we obtain distinguished invariant manifolds and time scale estimates, and we highlight the special role of singular perturbation parameters in higher order approximations of slow manifolds. We close the paper with a discussion of distinguished invariant manifolds in the global phase portrait.

Introduction

Cellular function involves a large network of transformations of substrates, denoted S, into products, P, which in turn may be further transformed, eliminated, or cycled back into a useful form. While the chemical conversion of S into P can occur spontaneously the rate constant, k, that regulates the speed of the reaction (1.1) will often be very small, so that spontaneous conversion is too slow to sustain life. Moreover, spontaneous conversion allows only the crudest forms of control. Consequently, the reaction must be catalyzed or “sped up.” Enzymes, denoted E, are biochemical catalysts that accelerate the conversion of S into P, and the chemical process by which the conversion of a substrate molecule into a product molecule is accelerated by an enzyme is called an enzymatic reaction.

The simplest description of an enzymatic reaction for a single-substrate, single-product reaction is the Michaelis–Menten mechanism [23, 29, 50], in which the conversion of S into P is achieved via two elementary reactions: the reversible formation of the enzyme-substrate complex, C, and the conversion of S to P in the complex C with (in this simple model) simultaneous dissociation into E and P. Enzymes lower the free-energy barrier separating reactants from products, with the result that (1.2) is generally faster than (1.1) by many orders of magnitude [30, Section 6.2].

The modeling and quantification of enzymatic reaction rates is of particular importance, especially since metabolic disease and dysfunction may arise when these reactions are too slow due, e.g., to a mutation in the corresponding gene. At or near the thermodynamic limit, enzymatic reactions are modeled by nonlinear ordinary differential equations (ODEs), known as rate equations, that obey the law of mass action. While the nonlinear terms in the model equations of enzymatic reactions make the mathematical treatment of the reaction mechanism challenging, avenues for simplification often exist. Specifically, if the rates of the elementary reactions that comprise the catalytic reaction are disproportionate, the ODE model will be multiscale, meaning the complete reaction will consist of disparate slow and fast time scales. Under the influence of distinct fast and slow time scales, the rate of change of c (using lower-case italic letters to represent the concentrations of the corresponding species) is very small relative to the rate of rate of change of s. The exploitation of this almost negligible rate of change warrants a simplification of the form where keff is the effective—but non-elementary—rate function. In the case of the Michaelis–Menten mechanism, keff is a hyperbola in the variable s; in more complicated mechanisms it may adopt the form of, for instance, a Hill-type function. The advantage offered by (1.3) is that the entire reaction is describable in terms of the reactant concentration, s, since the explicit dependence on e and c has been eliminated. The most widely studied example of this kind of reduction is probably the Michaelis–Menten rate law, which can be obtained using the standard quasi-steady-state approximation (sQSSA). More generally, rate laws of the form (1.3) are referred to as quasi-steady-state (QSS) reductions or quasi-steady-state approximations (QSSA). The term QSS speaks to the fact that the concentration of at least one chemical species (typically an intermediate) changes very slowly for the majority of the reaction.* In fact, the rate of change is so small that it is nearly zero (steady-state) but not quite; hence the expression quasi-steady-state.

The principal value of QSS approximations is that they yield a reduction of dimension [10]. In the biochemical arena, the related equilibrium approximation was initially justified via biochemical arguments by Henri [23] and by Michaelis and Menten [29]. Briggs and Haldane [4] later provided a mathematical justification of the sQSSA using an argument that hints at later singular perturbation treatments but lacked formal justification. Only the development of singular perturbation theory some decades later (with seminal contributions by Tikhonov [46], and later Fenichel [9]) laid a solid mathematical foundation, which was used by Heineken et al. [22] to develop criteria for the validity of the sQSSA for the closed Michaelis–Menten system. This history was paralleled in inorganic chemistry, with the initial development of the sQSSA based on ad hoc chemical reasoning [2, 6], followed eventually by more rigorous treatments based on singular perturbation theory [3].

Singular perturbation theory in this context applies to ODEs that depend on a small nonnegative parameter ε, and admit non-isolated stationary points at ε = 0. In practice, e.g. for systems with polynomial or rational right-hand side, the set of stationary points then contains a submanifold of positive dimension, which is called a critical manifold. Given appropriate conditions (see Appendix A for details), one obtains a reduction to a system of smaller dimension constrained to evolve on the critical manifold. The challenge in any application of Fenichel theory resides in finding a small parameter from a given parameter dependent system. Traditional analyses of enzymatic reactions rely heavily on scaling and non-dimensionalization in order to transform the model equations into a standard form, and the utility of scaling analysis is that the small parameter often emerges naturally from the dimensionless equations [41]. A different, more recent approach [16] starts with determining so-called Tikhonov–Fenichel parameter values (TFPV), by searching for parameter combinations at which the system admits non-isolated stationary points, and satisfies further technical conditions (see Section 3.1). From such (dimensional) TFPV one then obtains singular perturbation reductions via small perturbations along a curve in parameter space. In chemical applications, critical manifolds frequently emerge when specific system parameters (such as rate constants) vanish.

While singular perturbation theory provides a very satisfactory toolbox for reduction of chemical reaction networks, examples from the literature indicate that the approach may be too narrow for some applications. Thus in some scenarios, at a certain parameter value there exists a distinguished invariant manifold which is, however, not comprised of stationary points. Formally, this means that a QSS reduction which approximates the system when 0 < ε ≪ 1 is not attributable to Fenichel theory. Nevertheless the QSS reduction is still sometimes a good approximation to the full system when Fenichel theory is inapplicable, and this raises several important questions. First, given the lack of a critical manifold and a fixed reduction procedure, how does one justify a QSS reduction, and how does one go about quantifying its efficacy? Second, if Fenichel theory is not applicable but a QSS reduction still proves to be an accurate approximation, will there be a distinguished invariant manifold that attracts nearby trajectories? In other words, what phase-space structures make the QSS reduction possible in situations where Fenichel theory is extraneous? In the present paper we contribute, on the one hand, to answering these questions for an open Michaelis–Menten system with constant substrate influx. On the other hand, we provide sharper estimates for the accuracy of the sQSSA in singular perturbation scenarios. Finally, we consider distinguished invariant manifolds from a global perspective for the system on the Poincaré sphere.

An open Michaelis–Menten reaction mechanism

The open Michaelis–Menten reaction mechanism we consider here is the classical Michaelis–Menten reaction mechanism with a constant influx of substrate, S, at a rate k0: where k0, k1, k−1 and k2 are rate constants.

Mathematical models for (2.1) come in both deterministic and stochastic forms. Here we consider only the deterministic ODE model that follows the law of mass action near the thermodynamic limit. For a thorough analysis of the chemical master equation corresponding to (2.1), we invite the reader to consult [1, 45].

The mass action model corresponding to (2.1) is given by the following set of nonlinear ODEs: where a dot denotes differentiation with respect to time. Summing equations (2.2b) and (2.2c) reveals the conservation law where e denotes the total enzyme concentration. Employing (2.3) to eliminate (2.2c), and noting that (2.2d) is not coupled to (2.2a) or (2.2b), yields the simplified model from which the time dependence of p and e are readily obtained from (2.2d) and (2.3) once the solution to (2.4) is procured.

In contrast, the mass-action system for the closed Michaelis–Menten reaction mechanism is recovered by setting k0 = 0:

One distinguishing difference between the open and closed systems is that the total substrate concentration, s, is a conserved quantity when the reaction is closed. Therefore, (2.2) with k0 = 0 is equipped with the additional conservation law s = s + c + p, whereas with k0 > 0 one has only one conservation law, (2.3).

It is well known that further simplification of (2.5) is possible via a QSS reduction. The most common reduction is the sQSSA, in which (2.5) is approximated with a differential-algebraic equation consisting of the algebraic equation obtained by setting the right-hand side of equation (2.5b) equal to zero (“ċ = 0”) along with the differential equation (2.5a). This reduces to the single differential equation where K is the Michaelis constant.

The legitimacy of the sQSSA (2.6) for the closed Michaelis–Menten reaction mechanism (2.5) is well-understood. Following an early effort by Briggs and Haldane [4], Heineken, Tsuchiya, and Aris [22] were perhaps the first to prove with some degree of rigor that (2.6) is valid provided e ≪ s0. The qualifier, e ≪ s0, was justified via singular perturbation analysis. Defining , , and T := k1e generates the singularly perturbed dimensionless form of (2.5) where prime denotes differentiation with respect to T, λ := k2/k1 s0, κ := K0, and μ := e/s0. Consequently, the sQSSA (2.6) is justified via Tikhonov’s theorem [46]. Throughout the years, refinements and variations of the condition μ ≪ 1 have been made. Perhaps most famously, Segel [42] and Segel and Slemrod [41] extended the results of Heineken et al. [22] and demonstrated that (2.6) is valid whenever e ≪ K + s0. Embedded in Segel’s estimate is the more restrictive condition, e ≪ K, which is independent of the initial substrate concentration, and is nowadays the almost universally accepted qualifier that justifies (2.6) [7].

While the QSS reductions of the closed Michaelis–Menten reaction are well-studied, analyses pertaining to the validity of the QSSA in open reaction environments are somewhat sparse [1, 18, 44, 45]. The question we address is therefore: when is further reduction of (2.4) possible? The trajectories illustrated in Figure 1 show that there are certainly conditions under which the QSSA estimate of the enzyme-substrate complex, given by Eq (2.6b), which applies equally to the open system, is close to an invariant manifold (here, a trajectory) that attracts nearby trajectories and along which the equilibrium point is eventually approached from almost all initial conditions [11, 20, 39]. We thus ask under what condition is the open sQSSA permissible? At first glance, it seems rather intuitive to postulate that the open sQSSA (2.8) is valid under the same condition that legitimizes the closed sQSSA: e ≪ K. In fact, following the earlier work of Segel and Slemrod [41], Stoleriu et al. [44] introduce the parameter α := k0/(k2e), and suggest that (2.6b) and (2.8) are applicable whenever holds. The inequality (2.9) is less restrictive than the Segel and Slemrod condition, since (2.9) is satisfied as long as k0 is sufficiently close to k2e [Implicitly, Stoleriu et al. assume that α < 1 in Eq (2.9).]

Figure 1.

Trajectories of the open Michaelis–Menten equations (2.4) for (a) k1 = 1, e = 1, k−1 = 1 k2 = 3 and k0 = 2.5 (in arbitrary units), i.e. under conditions where there is an equilibrium point in the first quadrant, marked by a dot; and (b) with parameters as in (a), except k0 = 3.5, under which conditions there is not an equilibrium point in the first quadrant, and the s component of the solution grows without bound. The arrows show the direction of the flow. The dashed curve in both figures is defined by the QSSA equation (2.6b).

The approach used to derive (2.9) was based on the traditional method of comparing time scales: a singular perturbation parameter was recovered through scaling analysis of the mass action equations (2.4). However, it is possible to derive erroneous conclusions regarding the validity of the QSSA, even when great care is taken in scaling and non-dimensionalization methodology (see, for example [13], Section 4). It thus seems prudent to reexamine the basis for the sQSSA in the open Michaelis–Menten mechanism using tools of singular perturbation theory that go beyond scaling arguments.

The Quasi-Steady-State Approximation: Justification from singular perturbation theory

In this section we derive the QSSA directly from Fenichel theory. Details covering projection onto the critical manifold can be found in Appendix A.

The critical manifolds: Tikhonov–Fenichel parameter values

To apply Fenichel theory to the open Michaelis–Menten reaction mechanism, we need a curve of non-isolated equilibrium solutions to form in the first quadrant of ; see [16]. The following Lemma addresses the conditions that ensure the existence of a critical manifold, and records some general qualitative features.

Lemma 1. (a) System (2.4) admits an infinite number of stationary points if and only if one of the following conditions holds.

k0 = k1 = 0;

k0 = e = 0;

k0 = k2 = 0.

If the number of stationary points in the plane is finite then it is equal to zero or one. There exists one stationary point if and only if the genericity conditions are satisfied. In that case the stationary point is equal to

This point lies in the first quadrant if and only if in which case it is an attracting node. The stationary point lies in the second quadrant if and only if k2e – k0 < 0, in which case it is a saddle point.

The first quadrant is positively invariant for system (2.4), and solutions starting in the first quadrant exist for all t ≥ 0. When k−1 + k2 > 0 then every solution that starts in the first quadrant enters the (positively invariant) subset defined by c ≤ e.

System (2.4) admits no non-constant closed trajectory.

Sketch of proof. Parts (a) and (b) are straightforward, as is the first statement in part (c). For the second statement note ṡ + ċ ≤ k0, hence solutions starting in the first quadrant remain in a compact set for all finite t > 0. Finally, when c ≥ e then (2.4) shows that ċ ≤ −(k−1 + k2)e, hence the second statement of part (c) holds. We turn to the proof of part (d): If there exists a non-constant closed trajectory then its interior contains a stationary point. Given a degenerate situation from part (a), the variety of stationary points is unbounded, hence would intersect a closed trajectory if it intersects its interior; a contradiction. This leaves the setting with an isolated stationary point, necessarily of index one, which is only possible when the stationary point (3.2) lies in the first quadrant. By part (c) the closed trajectory must be contained in the strip defined by c ≤ e. But in this strip the divergence of the vector field equals −(k1(e – c) + k1 s + k−1 + k2) < 0, and no closed non-constant trajectory can exist by Bendixson’s criterion. □

Remark 1. The case k0 > k2e, in which the inflow exceeds the enzyme’s clearance capacity, is not physiologically irrelevant since the gene coding for a particular enzyme may suffer a mutation that results in an enzyme with reduced catalytic activity, for example. As a rule, the accumulation of a metabolite will eventually become toxic (or possibly oncogenic) to the cell, and the rate at which S accumulates is therefore of interest. Other situations, e.g. the existence of an alternative but less efficient pathway for eliminating S, or the permeation of S through the cell membrane, would require more elaborate models for their study. Nevertheless, the model under study here would yield useful initial insights into the cellular effects of a mutation to an enzyme.

Lemma 1 ensures the existence of a critical manifold comprised of equilibrium points whenever k0 vanishes along with either e, k1 or k2. We note that in the context of the closed reaction (2.5), parameters with e = 0 (with all remaining parameters > 0), resp. k1 = 0 (remaining parameters > 0), resp. k2 = 0 (remaining parameters > 0), are TFPV. Generally, a TFPV [ ] is characterized by the property that the variety of stationary points has positive dimension, and a generic small perturbation of the parameters results in the formation of a normally hyperbolic invariant manifold, called a slow manifold [17]. Choosing a curve in parameter space which is parameterized by ε and starts at a TFPV, one obtains a system which admits a singular perturbation reduction.

Let denote the parameter vector: π := [k0 e k1 k2 k−1]. By Lemma 1, and upon checking normal hyperbolicity (see Appendix A below), the TFPVs and the critical manifolds, M, are as follows:

Fenichel theory ensures that perturbing π in (3.4) along a curve in parameter space through the TFPV results in the formation of an invariant slow manifold that attracts nearby trajectories at an exponential rate. Formally, the QSSA may be seen as an approximation of the dynamics on the slow manifold, perturbing from a TFPV.

Singular perturbations and the geometry of parameter space

The justification of the QSSA from singular perturbation theory requires us to implicitly equip parameter space with some additional geometric structure. For example, consider the case where both e and k0 vanish in the singular limit. In order to formally apply singular perturbation theory, it must hold that k0 ~ O(e) (in a sense discussed below). Generally speaking, this means that we can apply singular perturbation theory along a parametric curve, Γ, in (e, k0) parameter space, Γ := (e, z(e)), provided z(0) = 0 and

However, a small perturbation suggests that the parameter values will be close to the parameter plane origin located at (e, k0) = (0, 0). In this case z(e) is well-approximated by its tangent line at e = 0, thus it is enough to only consider rays of the form k0 = γe, where γ is a positive constant with dimension t−1. To eliminate the need for a dimensional slope γ, define a ray in parameter space by where the parameters and are nominal values of k0 and e, respectively.

The additional constraint of sampling parameter space along a ray [or in a more general way along a curve satisfying (3.5)] must be imposed in order to justify the open sQSSA from singular perturbation theory. In their analysis of the open Michaelis–Menten reaction (2.4), Stoleriu et al. [44] implicitly performed their analysis along a ray defined by where e(0) > 0 is the initial free enzyme concentration. This ray in parameter space is encoded in their initial conditions, which allow for an arbitrary positive value of e(0), but which specify . The advantage of working along the ray defined by (3.7) is that there is no possibility that the inflow can exceed the clearance capacity of the enzyme, i.e. inequality (3.3) is automatically satisfied.

In order to apply singular perturbation theory, we need to start from a critical manifold, i.e. from one of the cases in the set (3.4). Note that the ray through the (e, k0) parameter plane chosen by Stoleriu et al. [44], Equation (3.7), does not satisfy (3.5) unless e(0) = 0. This leads to difficulties. For example, the condition (2.9) along the ray defined by (3.7) translates to

The inequality (3.8) is satisfied by taking k1 → 0, but this limit alone does not produce a critical manifold. Hence, the singular perturbation machinery is not applicable to legitimizing the open sQSSA (2.8) by this route.

Another issue with the constrained set of initial conditions imposed by (3.7) is that it excludes many initial conditions that are physiologically relevant. For example, a natural initial condition is (s, e, c, p)(0) = (0, e, 0, 0), corresponding to the substrate flow being turned on at time zero (e.g. because the cell is placed in a new environment, or because it has turned on a previously dormant metabolic pathway that produces S), but this initial point is inaccessible if the parametric constraint (3.7) has been imposed. Consequently, it remains an open question whether the results of the analysis apply at arbitrary points in parameter space and for arbitrary initial conditions. In particular, there is no guarantee that the analysis of Stoleriu et al. [44] applies when the inflow exceeds the clearance capacity of the enzyme which, as argued previously, is not an irrelevant case. By contrast, a transformation informed by the basic requirements of singular perturbation theory such as (3.6) allows us to make rigorous statements about the manifold structure of the problem, and imposes no constraints on the initial conditions.

Quasi-steady-state reductions: Projecting onto the critical manifold

Let us now consider the first scenario in which e and k0 vanish in the singular limit. The perturbation of the singular vector field is

The singular limit obtained by setting ε = 0 in (3.9) yields a critical manifold, M, that is identically the s axis:

To compute the corresponding singular perturbation reduction (see Appendix A for specific details), we rewrite the right hand side of (3.9) as h(s, c) + εG(s, c, ε). Furthermore h(s, c) = P(s, c) f (s, c), with

Since DfP = −(k1 s + k−1 + k2) is negative on M, M satisfies the attracting hyperbolicity condition, and Tikhonov-Fenichel reduction is applicable (see [15] and also Appendix A). The singular perturbation reduction is then obtained by projecting G(s, c, 0) onto the tangent space of M at x = (s, c), TM, via the linear operator Π which projects onto the kernel along the image N of the Jacobian Dh(x).

Note that N, which is equal to the range of P(x), is a complementary subspace to TM. For our specific problem (3.9), Π is given by and the corresponding reduction, which agrees with the QSS reduction, is

Equation (3.14) is, of course, the open sQSSA. A similar calculation is easily carried out for the case of small k1 and small k0, as well as small k0 and k2, and we refer the reader to Appendix A for details. The specific QSS reduction that accompanies the perturbation defined by and is which is the linear limiting law obtained in the small-s limit of (3.14).

Accordingly, we have confirmation that the open sQSSA (2.8) is valid under any condition that invokes a scaling of the form and . We further note that a QSS reduction based on Fenichel theory is also possible in case and so that both k0 and k2 vanish in the singular limit. This reduction yields the classical equilibrium approximation, (See Section 5.3 and Appendix A for details.)

Several questions remain. First, what is ε? We have shown that the open sQSSA is valid provided k0 and e are sufficiently small, but what is small when k0 and e are nonzero? Second, from the work of Goeke et al. [17], the QSS may still hold in certain regions of the phase-plane even if Fenichel theory is not applicable. The analysis of Stoleriu et al. [44] is also indirectly suggestive of the idea that the validity of the open sQSSA may not necessarily stem from singular perturbation theory. These observations raise the deeper question: is a scaling of the form , necessary for the validity of the QSSA, or merely sufficient? We address these questions directly in the sections that follow.

Revisiting quasi-steady state for the complex species

The notion of QSS

Singular perturbation theory provides a natural setting for developing conditions under which QSSA holds, but the literature (notably Stoleriu et al. [44] for the open Michaelis–Menten mechanism) suggests that one should consider less restrictive notions as well. In the following we will sketch one such notion. This goes back to Schauer and Heinrich [40], who were the first to note that the minimal requirement for the validity of QSS reduction for some species should be the near-invariance of a corresponding variety, which we call the QSS variety. This variety is defined as the zero set of the rate of change for the species under consideration. The idea of near-invariance was expounded upon by Noethen et al. [32], and further analyzed by Goeke et al. [17]:

A fundamental feature of QSS is that the rate of change of certain sets of species should be close to zero for an extended period of time. (In the Michaelis–Menten reaction, QSS for complex thus means that ċ ≈ 0 for an extended period of time.) In the phase space interpretation, a sizable part of the trajectory should thus be close to the QSS variety which is defined by evaluating the condition ċ = 0. (In the Michaelis–Menten mechanism the QSS variety for complex is thus defined by k1(e – c)s – (k + k2)c = 0; see Eq (2.6b).) The validity of such a condition will depend on the parameters.

According to [17], Section 3.3, the minimal requirement for QSS should therefore be near-invariance of the QSS variety, in the sense that the system parameters are small perturbations of QSS parameter values. By definition, at a QSS parameter value the QSS variety is an invariant set for system (2.4). (In the Michaelis–Menten mechanism one thus requires invariance of the variety defined by Eq (2.6b) for system (2.4) at a QSS parameter value.) The arguments in [17] show that this condition is necessary if one requires arbitrary accuracy of the QSS approximation for suitable parameters. By standard continuous dependence theorems for initial values and parameters (see e.g. Perko [34], Section 2.3), small perturbations of a QSS parameter value yield trajectories that remain close to the QSS variety on compact time intervals; thus the condition is also sufficient. One practical advantage of this notion is that QSS parameter values, similarly to TFPV, are algorithmically accessible for polynomial or rational systems.

The near-invariance condition alone may not be sufficiently strong to satisfy expectations about QSS. One may also require that solutions quickly approach the QSS variety in an initial transient phase. Since the combination of these two features is automatically satisfied in singular perturbation settings, singular perturbations naturally enter the picture. But the singular perturbation scenario is both broader and narrower than QSS for chemical species: It is broader since it also is applicable to settings with slow and fast reactions. On the other hand, we will see below that it is, in a sense, too narrow for sQSS in the open Michaelis–Menten reaction mechanism.

Open Michaelis–Menten: QSS parameter values for complex

The QSS variety for (2.4) is given by

We prefer this to the usual notation w(s) = es/(K + s), which may obscure the role of k1. We first determine all QSS parameter values.

Lemma 2. The QSS parameters of system (2.4) are as follows:

e = 0 with the other parameters arbitrary;

k1 = 0 with the other parameters arbitrary;

k0 = k2 = 0;

k−1 = k2 = 0.

Proof. We proceed along the lines of [17], Section 3.4, using an invariance criterion that employs the Lie derivative, L[·], corresponding to (2.4). The Lie derivative is defined by for any polynomial (more generally, smooth) function φ. For the variety defined by φ = 0 to be invariant it is necessary that

Moreover, the condition is sufficient when φ is irreducible, and it is applicable to the irreducible factors of φ; for details see [17] and the references therein.

Now let ψ(s, c) = 0 define the QSS manifold, thus

The invariance condition for the curve ψ(s, c) = 0 is whenever ψ(s, c) = 0, thus

This product yields three conditions which can be evaluated. Clearly k1 = 0 works and yields (ii). The second condition, e – c = 0, holds on ψ = 0 if and only if (k−1 + k2)e = 0, which yields respectively (iv) and (i). The third condition yields k0 = k2 = 0 when k2 = 0, i.e. (iii). In case k2 ≠ 0 one obtains c = k0/k2, and here the coefficient of s and the constant must vanish. This again leads to conditions already discussed. □

Remark 2. (a) In cases (i) and (ii), the QSS variety is given by c = 0, provided that the other parameters are positive, and the QSS parameter conditions are less restrictive than for singular perturbations, which also require k0 = 0. This is a notable difference to the closed Michaelis–Menten scenario, for which all complex-QSS parameter values are also TFPV. Case (iii) corresponds to a singular perturbation scenario. The dynamics in case (iv) is of some interest in the Michaelis–Menten reaction mechanism without inflow; see [7].

Classical QSS reduction is tantamount to exploiting the fact that if ψ(s, c) = 0 defines a nearly invariant curve, then c ≈ w(s), from which the open sQSSA (2.8) presumably follows. However, a word of caution is in order. When a QSS parameter value is also consistent with a singular perturbation and gives rise to a critical manifold, the classical QSS reduction may differ from the reduction obtained from Fenichel theory (see [17], Section 3.5). For example, ψ(s, c) = 0 is nearly invariant if k0 and k2 are small, but the classical QSS reduction, given by does not agree with the reduction obtained from singular perturbation theory, which is given by (A.18). Convergence to the singular perturbation reduction is guaranteed by Fenichel theory, hence the QSS reduction (4.5) cannot correctly describe the dynamics at lowest order.

For the QSS parameters which do not correspond to singular perturbations, there remains to investigate whether solutions approach this variety, and if so, how fast and how close the approach is. Furthermore, even in the singular perturbation scenario one needs estimates on the initial (boundary layer) behavior, since Fenichel’s theory applies directly only to a neighborhood of the critical variety.

These problems will be addressed via direct estimates, which will also be of help in answering a quantitative question, i.e. how small should e be in order to justify (2.8)? Ultimately, the term small is relative in nature. Therefore, the appropriate question to ask is: For (2.8) to be approximately accurate, e and k0 must be much smaller than what? Before we start this investigation we establish an auxiliary result about the phase plane geometry of (2.4).

Phase plane arguments

From here on we restrict attention to system (2.4) on the positively invariant strip W defined by s ≥ 0 and 0 ≤ c ≤ e. A priori we impose no requirements on the parameters. We look at isoclines, noting that and where denotes the x nullcline. These nullclines define positively invariant sets:

Lemma 3. Consider the “wedge”

Then the following hold:

If the system admits no positive stationary point, thus k0 > k2e, then the c-isocline lies above the s-isocline for all s ≥ 0, and W1 extends to s → ∞. If the system admits the positive stationary point (, ) then the isoclines meet at this point, and , for all points of W1.

W1 is positively invariant for system (2.4), and on W1 one has ṡ ≥ 0.

Proof. Part (a) is straightforward. As for part (b), from (2.4) one sees that on W1, thus

This implies the positive invariance of W1, since the vector field points to the interior of W1 at the boundary (Figure 2). Clearly ṡ ≥ 0 on W1. □

Figure 2.

Sketches of the positively invariant sets W1 in the phase plane for the open Michaelis–Menten reaction mechanism (2.1). The curves are the nullclines, and the arrows show the direction of motion of trajectories as they cross the nullclines. Both nullclines tend asymptotically to c = e as s → ∞. is the s intercept of the s nullcline. Left: k0 > k2e and the two nullclines never meet. Right: k2e > k0 and the nullclines cross at the stationary point (, ). The flow points into the region delimited by the two nullclines, making this region a funnel [24].

Remark 3. Smallness of e and existence of a positive stationary point imply smallness of k0; this leads automatically to the singular perturbation setting. Matters are different when k1 is small.

Everywhere inside the wedge, ċ > 0 and ṡ > 0. Thus, all trajectories inside the wedge have positive slope. Since the flow points into the wedge, the slow manifold must also lie inside the wedge. Thus, the slow manifold has a positive slope for in the first quadrant. Moreover, the slow manifold must enter the first quadrant by crossing through the s axis in the interval (0, ), where is the s intercept of the s nullcline (Figure 2).

In the case that there is a positive equilibrium point, for , ċ < 0 and ṡ < 0 between the two nullclines so that trajectories in this region still have positive slope. The flow is, again, into the region between the two nullclines (Figure 2), so the slow manifold must lie within this region. The slow manifold therefore has positive slope here as well. Moreover, . Thus, the two nullclines pinch together asymptotically. Although we do not pursue this idea here, this property would allow the anti-funnel theorem to be used to prove the existence of a unique slow manifold to the right of the equilibrium point [5, 24]. (See Section 6 for correspondence to the global behavior.)

How small is small: A direct estimate

Given that we are interested in obtaining a condition that ensures phase plane trajectories closely follow the QSS variety corresponding to the c nullcline, we compute an upper bound on the limit supremum (lim sup) of for a solution of (2.4), where w(s) is given by (4.1). To determine such an upper bound, we calculate

The derivative ċ given in (2.4) factors: and substitution of (4.9) into (4.8) yields

Differentiating w(s) with respect to s reveals max ∣w′(s)∣ = k1 e/(k−1 + k2). Denote max ∣ṡ∣ by v and note that v ≤ k0 on W1, due to ṡ ≥ 0.

With

Cauchy’s inequality

Implies which yields

A natural choice for σ is σ := τ0/2 leading to the inequality

Applying Gronwall’s lemma to (4.14) generates an upper estimate for L2:

Proposition 1. (a) For every solution of (2.4) with initial value in W1 one has the estimates with τ := τ0t = (k−1 + k2)t.

Thus with the solution approaches the QSS variety up to an error of ε*2, with time constant λ := (k−1 + k2)−1.

Note that the estimates from the proposition explain the rapid approach of the trajectories in Figure 1 to the QSS variety.

From our analysis of the mathematical energy, L2, we have both a time constant, , as well as a parameter, ε*. The time constant yields a natural dimensional fast time scale, τ, that is equivalent to the fast time scale obtained by Segel [42] for the closed Michaelis–Menten reaction mechanism. Moreover, ε* should in some sense be small for the open sQSSA to be accurate. The difficulty here is that ε* has dimension, and we must scale ε* appropriately to recover a dimensionless parameter. To scale, note that if k0 < k2e, then

Since c ≤ e, we divide the (4.17) through by e, and take the inequality, to be the general qualifier for the validity of open sQSSA (2.8) in W1, when a finite stationary point is located in the first quadrant.

Note that ε vanishes if either k1, e or k2 vanish. However, the use of Fenichel theory also requires k0 to vanish in the singular limit, otherwise the perturbation is non-singular and the accuracy of a specific QSS reduction is attributable only to the near-invariance of the QSS manifold [hence the difference in the justification of the open sQSSA that occurs from the mapping versus the mapping ]. This observation is a definitive difference between our work and that of Stoleriu et al. [44].

Additional insights from solutions of the invariance equation

The sQSSA can be thought of as an attempt to approximate the slow invariant manifold [11]. There are many other methods for approximating the slow manifold, ranging from the method of intrinsic low-dimensional manifolds [28], which is accurate to O(ε) [25], to methods that can be improved order-by-order such as singular-perturbation theory [3, 22, 41], computational singular perturbation theory [27], and Fraser’s iterative method [11, 31]. Here, we study solutions of the invariance equation, the equation that the exact slow manifold satisfies, in order to gain further insights into the role of the TFPV in determining the validity of the sQSSA. The Fraser iterative method will be a major tool, but we will also consider various small-parameter expansions of the iterates.

The invariance equation

Assume that, in accordance with the arguments in Section 4.3, we can represent the slow manifold (at least locally) as the graph of a function c = C(s). If ṡ = ṡ(s, c) and ċ = ċ(s, c), then differentiating the assumed representation of the slow manifold with respect to time, we get the invariance equation [11, 19, 21, 26, 35].

The invariance equation could be solved using a perturbation method. A strategy suggested by the work of the previous sections is to perturb from a TFPV along a curve in parameter space with the TFPV as its endpoint, e.g. the ray (3.6). The scaling parameter ε can then serve as a perturbation parameter, and a perturbation problem of the typical form results, i.e. to compute the i’th term in the perturbation series, we solve an algebraic equation that only depends on the previous terms. However, suppose that we did not know about TFPVs. Then we might try to use the same small parameter as in the closed system, viz. some scaled version of e [4, 22, 42, 44]. In the current framework, we would write , and expand If we implement this program, we find that χ1(s) satisfies the differential equation

Higher-order terms also satisfy differential rather than algebraic equations. These difficulties are linked to the fact that e = 0, of itself, is not a TFPV for the open system. For the TFPVs (3.4a) and (3.4b), since the leading-order term in is O(ε), rescaling the TFPVs balances the terms in the invariance equation such that, to leading order, ċ, ṡ and are all O(ε). As a result, (with slight abuse of notation) dχ/ds first appears to O(ε), and we obtain an algebraic equation for χ The case of TFPV (3.4c) is slightly different. If we rescale and take , the ε0 terms of the invariance equation can be rearranged to

The term in square brackets gives us the critical manifold (3.4c) for ζ0. (The other solution, dζ0/ds = −1, gives the fast foliations of the manifold in the limit ε → 0.) At higher orders, ζ first appears with the O(ε) terms. However, because in the limit ε → 0 for this TFPV, the k0 term in ṡ vanishes, the coefficient of dζ at O(ε) is the term in square brackets in Eq (5.3), which vanishes. Thus, dζ first appears with a non-vanishing coefficient at O(ε), and we again have a perturbation problem involving only algebraic equations.

To recapitulate, rescaling the TFPVs yields a tractable perturbation problem precisely because the TFPVs define critical manifolds. Choosing any path through parameter space that does not reduce to a TFPV as ε → 0 will, by contrast, necessarily yield a troublesome perturbation problem.

Each TFPV yields a different perturbation problem. Rather than studying the perturbation expansions of the slow manifold directly, we turn to Fraser’s method [11, 31], which will allow us to compute a sequence of approximations in a TFPV-agnostic manner. Series expansions of the approximations can then be obtained for any desired TFPV scaling parameter.

In Fraser’s iterative method, we think of the invariance equation as an equation to be solved for in terms of . In this case, we can explicitly rearrange the invariance equation to the functional equation [11]

Observe that if we rescale k0 and e as in (3.6) and let ε → 0 in the functional equation, we recover the critical manifold (3.4a). Similar comments can be made for (k0, k1) and (k0, k2) and the corresponding critical manifolds (3.4b) and (3.4c), respectively. Thus, the critical manifolds are recovered in suitable limits of the functional equation.

We now want to solve Eq (5.4) for the slow manifold. If we knew the derivative of with respect to s along the slow manifold, we could immediately compute from (5.4). Since we do not, we solve the invariance equation by iteration: From some initial guess , we compute the derivative, substitute it into (5.4) to obtain , and iterate. In practice, iterative solution of a functional equation such as (5.4) tends to converge specifically to the slow manifold [11, 36] if it converges at all [37], even though every trajectory is a solution of the invariance equation.

The critical manifolds associated with the TFPVs suggest potential initial functions for iteration. Suppose then that we start iteration from the critical manifold [under either TFPV (3.4a) or (3.4b)] . Then is the sQSSA (2.6b). Figure 3 shows a sequence of iterates calculated from this initial function. Convergence is rapid, although much more so away from the s axis.

Figure 3.

Iterates of Eq (5.4) for the open Michaelis–Menten reaction mechanisms starting from the initial function for the parameters of Figure 1(a). The solid dot marks the location of the equilibrium point. The inset shows an expanded view of the behavior of the iterates near the origin.

As a side note, consider using a vertical initial function, i.e. one for which . The first iterate from such an initial function is the s nullcline, which intercepts the s axis at s = k0/k1e, i.e. at the extreme right end of the possible range of s intercepts of the slow manifold. The sQSSA, on the other hand, is the c nullcline, obtained in one iterative step from the initial function , and it intercepts the s axis at s = 0. The two nullclines thus arise naturally as approximations of the slow manifold by iteration from coordinate axes, and serve as upper and lower bounds for the slow manifold. Similar comments about the relationship of the functional equation to the nullclines have previously been made about closed systems [11, 12, 31].

The TFPVs (k0, e) and the small parameters rediscovered

If we obtain higher iterates using a symbolic algebra system, then make the substitution (3.6), and finally expand in powers of ε, we find that the i’th iterate is consistent with the previous iterate to order ε. In other words, the iterative method builds the perturbation series term-by-term, as was previously observed for various perturbative solutions of the closed system [25, 36]. However, this property does not hold if we, for instance, expand in powers of e, since e = 0 is not, of itself, a TFPV for the open system. These properties parallel those of the direct perturbation calculations (not shown).

The first two non-zero terms of the perturbation series computed along the ray (3.6) can be written as follows:

Division by , the nominal value of the enzyme concentration, has made this expression dimensionless. Thus, the ε2 term represents an error term for the sQSSA. Specifically, the absolute value of the coefficient of ε2, is a dimensionless error parameter such that the error in the sQSSA is small provided this coefficient is small. An elementary calculation shows that δ(s) has a local maximum of in s ∈ (0, ∞) provided . The global maximum of δ(s) for s ≥ 0 is either this local maximum or where the dimensional parameter ε* is defined in Eq (4.16). When the inflow exceeds the enzyme’s clearance capacity, the situation is straightforward, and δ(0) is the correct small parameter. Otherwise, δ will be larger than δ(0) when

We dropped the asterisks here because . This inequality can be solved numerically. It yields k0/k2e < 0.0767. Putting it all together, we have the following:

The sQSSA is a good approximation to the slow manifold globally if k0/k2e < 0.0767 and δ ≪ 1. Comparing Eq (5.7) to (4.18), and noting that in this parameter range, , we conclude that ε ≪ 10 is sufficient for the validity of the sQSSA. This is a somewhat more permissive bound than (4.18).

If k0/k2e > 0.0767, then δ(0) ≪ 1 is the appropriate condition for the validity of the sQSSA in the open system.

Note that this analysis has recovered both of the small parameters identified in Section 4.4, but has also established a sharp boundary for switching from one small parameter to the other. We thus have two complementary methods to obtain small parameters. In any given problem, one or the other method might be unworkable, thus our presentation of both methods here.

The TFPVs (k0, k2) and the equilibrium approximation

We can also expand the iterates using the small parameter implied by (3.4c). If we take , and then expand the second (or higher) iterate in powers of ε, we get where K = k−1/k1. Note that the O(ε0) term is the classical quasi-equilibrium approximation (QEA) for the Michaelis–Menten reaction mechanism. Contrast equations (5.5) and (5.10): The QEA for the open system is only accurate to order ε0, unlike the sQSSA which is accurate to order ε. (This fact is also reflected in Remark 2(b) and in the second example in Appendix A.) This is easily understood given that the QEA lies above the sQSSA at any s > 0, and that the slow manifold, which enters the first quadrant by passing through the positive s semi-axis, lies below the sQSSA for . In the interval , the sQSSA will therefore always be closer to the slow manifold than the QEA. This is unlike the situation in the closed system, where the slow manifold lies between the QEA and sQSSA, and where it is possible to choose parameters such that one or the other approximation is more accurate near the origin. The difference is that the QEA is a nullcline in the closed system, but not in the open system. One implication of this result is that the TFPV (3.4a) is the most natural one to use as a basis for a geometric singular perturbation treatment of the slow manifold. (See also Appendix A for further notes on the expansion from the TFPV (3.4c).)

The TFPVs (k0, k1) and the linear regime

Finally, turning to the TFPV (3.4b), we define a perturbation parameter ε by

A perturbation series based on this small parameter is a polynomial in s due to the appearance of k1 and s together in the rate equations. The first nonzero terms of this series are

Substituting this series along with the parameter definitions (5.11) into ṡ from (2.4), we get, to lowest order in ε, where vmax = k2e and or, restoring the small parameters from (5.11),

This is of course the small-s linear limit of the sQSSA, the previously seen Eq (3.15). An alternative route to this equation is presented in Appendix A.

The open Michaelis–Menten reaction on the Poincaré sphere

It is worthwhile to consider the global behavior of system (2.4) and its distinguished invariant sets to illuminate the role of QSS varieties in a broader context, particularly for systems which do not admit a stationary point in the first quadrant.

It is a standard technique to extend planar polynomial ODE systems to the Poincaré sphere. A good description of the procedure is given in Perko [34], Section 3.10: Given a sphere in , let the phase plane be tangent to its north pole, and consider the bijective central projection from the upper half sphere to the phase plane. Then points on the equator of the sphere may be viewed as points at infinity for the planar system, with each line through the origin corresponding to a pair of antipodal points on the equator.[†] Finally, for the purpose of visualization one applies a parallel projection in the north-south direction from the upper hemisphere to the equatorial plane.

A discussion of the system on the Poincaré sphere thus allows us to understand the behavior of the planar system at infinity. Note that all solutions of the system on the Poincaré sphere, which is compact, exist for all , while this is not necessarily the case for solutions of (2.4) when t ≤ 0 or outside the first quadrant. Any reference to limit sets in the following arguments is to be understood for the system on the sphere. In our analysis we will mostly be interested in the first quadrant.

Stationary points at infinity (i.e. on the equator) for a polynomial planar system are of particular interest. Antipodal pairs of stationary points generally correspond to invariant lines for the homogeneous part of highest degree; see e.g. [48]. For system (2.4) with k1 ≠ 0 we thus need to consider the homogeneous quadratic part

This homogeneous vector field admits three invariant lines, viz.

The stationary points at infinity which are relevant for the first quadrant correspond to the rays and we call the corresponding stationary points at infinity P1, resp. P2. Moreover we denote by P3 the stationary point at infinity which corresponds to .

We present the pertinent results for system (2.4) on the Poincaré sphere (see Appendix B for computations and proofs).

Lemma 4. Assume that the genericity conditions (3.1) are satisfied. Then the following hold for the system on the Poincaré sphere.

The stationary point P1 at infinity is a degenerate saddle when k2e > k0, with the stable manifold contained in the equator. In case k2e < k0 this point is a degenerate attracting node.

The stationary point P2 at infinity is a saddle-node, with a repelling node part on the upper hemisphere.

The stationary point P3 at infinity is a repelling node.

We first describe the behavior of system (2.4) on the relevant part of the Poincaré sphere when there is an isolated stationary point in the first quadrant as illustrated in Figure 4.

Figure 4.

The system on the Poincaré Sphere for the open Michaelis–Menten reaction mechanism in case k2e > k0. The distinguished trajectory is colored green.

Proposition 2. Assume that the genericity conditions (3.1) hold, and let k2e > k0. Then every solution starting in the first quadrant converges toward as t → ∞. There is a unique distinguished trajectory that connects the saddle P1 at infinity to P0. Moreover this trajectory is asymptotic in the phase plane to the line c = e → −∞.

We turn to the case when P0 lies in the second quadrant; see Figure 5. Here, considering the system on the Poincaré sphere is necessary to understand the global dynamics, and moreover a proper understanding requires us to look beyond the first quadrant.

Figure 5.

The system on the Poincaré Sphere for the open Michaelis–Menten reaction mechanism in case k2e < k0. The distinguished trajectory is colored green.

Proposition 3. Assume that the genericity conditions (3.1) hold, and let k2e < k0. Then every solution that starts in the first quadrant converges to P1 as t → ∞, and the corresponding trajectory in the phase plane is asymptotic to the line c = e0 to P1.

Remark 4. In view of Proposition 3, the mathematically distinguished trajectory connecting P0 and P1 may be seen as a natural candidate for a “global” slow manifold in appropriate parameter regimes. We provide a few more details here. From the proof (see also Figure 5) one finds that the two components of the unstable manifold of P0 connect respectively to P1 and to the antipode of P3. Solutions in the open upper hemisphere, unless they start on the stable manifold of P0, converge either to P1 or to the antipode of P3 as t ψ ∞. Moreover one component of the stable manifold of P0 connects to P2 (which is the only available alpha limit point), and the other may connect either to the antipode of P1, or to the antipode of P2, or to P3. (Topological arguments do not yield more precise information, and we will not delve further into this matter.) The stable manifold of P0 separates the regions of attraction for P1 and the antipode of P3 in the open upper hemisphere. In turn, the region of attraction for P1 is separated by the distinguished trajectory into two subregions. For one of these subregions, the alpha limit set of all points is equal to {P2}, thus one may briefly say that all trajectories in this region come down from c = ∞. For the other subregion, a similarly concise statement does not seem possible: The set of alpha limit points certainly includes P3, but it may also include the antipode of P2 or of P1.

Discussion

The open Michaelis–Menten reaction mechanism, although of definitive relevance in biochemistry, has attracted less attention than the classical closed mechanism without influx. We investigated the sQSSA for this system from two perspectives:

We considered QSS as a singular perturbation phenomenon. We determined all parameter combinations (TFPVs) from which singular perturbation reductions emanate via a small perturbation, and in particular we identified the relevant parameter values for sQSSA, which are given by k0 = e = 0, with all other parameters positive. By singular perturbation theory one obtains the familiar QSS reduction.

Motivated by a more general notion of QSS (proposed by Schauer and Heinrich [40]) and by the results of Stoleriu et al. [44], we obtained quasi-steady state reduction by direct estimates from QSS parameter values given by e = 0, all other parameters positive. Thus sQSS reduction for the open MM reaction is applicable to a wider range of parameters than for singular perturbation reduction. Note that such a phenomenon does not appear in the closed Michaelis–Menten system. By these estimates we obtained a justification of central results in [44], and could also determine their range.

However, considering the fine structure of slow manifolds by analysis of higher order approximations revealed the special role (and higher accuracy of approximation) of Tihkonov-Fenichel parameter values in contrast to QSS parameter values.

Finally, we took a global mathematical perspective to investigate scenarios with no positive equilibrium.