Dynamical aspects of the total QSSA in Enzyme Kinetics

TL;DR

This paper provides a mathematical foundation for the quasi-steady state approximation in enzyme kinetics by linking it to the center manifold theory, enhancing understanding of enzyme reaction dynamics.

Contribution

It proves the asymptotic equivalence between QSSA and the center manifold in enzyme kinetics, offering a rigorous basis for mechanistic analysis.

Findings

QSSA is asymptotically equivalent to the center manifold.

Provides a mathematical foundation for enzyme reaction analysis.

Links Tihonov theory with enzyme kinetics approximations.

Abstract

In this paper we prove that the well-known quasi-steady state approximations, commonly used in enzyme kinetics, which can be interpreted as the reduced system of a differential system depending on a perturbative parameter, according to Tihonov theory, are asymptotically equivalent to the center manifold of the system. This allows to give a mathematical foundation for the application of a mechanistic method to determine the center manifold of (at this moment, still simple) enzyme reactions.

Equations

$$X+E\text{raise}0.94722pt\text{ext@arrow}3095\text{arrowfill@}\text{relbar}\text{relbar}\rightharpoonup k_1\text{lower}0.94722pt\text{ext@arrow}0359\text{arrowfill@}\leftharpoondown \text{relbar}\text{relbar}{k}_{-1}C\stackrel{k_2}{\to }{X}_p+E$$

$$X+C+X_p=X_T\text{and}E+C=E_T,$$

$$\dot{X}(t)=-k_{1}X(t)\left(E_T-C(t)\right)+k_{-1}C(t)$$

$$\dot{C}(t)$$

$$=k_1\left[X(t)\left(E_T-C(t)\right)-K_{M}C(t)\right]$$

$$X(0)=X_T,C(0)=0.$$

$$C(t)=\frac{E_{T}X(t)}{X(t)+K_M},$$

$$\dot{X}(t)=-k_{2}C(t)=-\frac{V_{max}X(t)}{X(t)+K_M},X(0)=X_T,$$

$$\dot{\bar{X}}(t)=-k_{2}C(t),$$

$$\dot{C}(t)=k_1\left[\bar{X}(t)E_T-\left(\bar{X}(t)+E_T+K_M\right)C(t)+C^2(t)\right],$$

$$\overline{X}(0)=X_T,C(0)=0.$$

$$\dot{\bar{X}}\approx -k{C}_{-}(\bar{X}),\bar{X}(0)=X_T,$$

$$C_{-}(\bar{X})=\frac{(E_T+K_M+\bar{X})-\sqrt{(E_T+K_M+\bar{X}{)}^2-4E_T\bar{X}}}{2}$$

$$\dot{x}=f(x),x\in {\mathbb{R}}^n.$$

$$\dot{y}=Ay,y\in {\mathbb{R}}^n,$$

$$E^s=\mathrm{span}\{e_1,...,e_s\},$$

$$E^u=\mathrm{span}\{e_{s+1},...,e_{s+u}\},s+u+c=n$$

$$E^c=\mathrm{span}\{e_{s+u+1},...,e_{s+u+c}\},$$

$$\left(\begin{array}[]{c}\dot{u}\\ \dot{v}\\ \dot{w}\end{array}\right)=\left(\begin{array}[]{ccc}A_{s}&0&0\\ 0&A_{u}&0\\ 0&0&A_{c}\\ \end{array}\right)\left(\begin{array}[]{c}u\\ v\\ w\end{array}\right),$$

$$\{\begin{array}{ll}& \dot{u}=A_{s}u+R_s(u,v,w),\\ & \dot{v}=A_{u}v+R_u(u,v,w),\\ & \dot{w}=A_{c}w+R_c(u,v,w),\end{array}$$

$$W^{s}_{loc}(0)=\Bigl{\{}(u,v,w)\in\mathbb{R}^{s}\times\mathbb{R}^{u}\times\mathbb{R}^{c}|v=h^{s}_{v}(u),w=h^{s}_{w}(u);$$

$$Dh^{s}_{v}(0)=Dh^{s}_{w}(0)=0;\ |u|\ \text{sufficiently small}\Bigr{\}}$$

$$W^{u}_{loc}(0)=\Bigl{\{}(u,v,w)\in\mathbb{R}^{s}\times\mathbb{R}^{u}\times\mathbb{R}^{c}|u=h^{u}_{u}(v),w=h^{u}_{w}(v);$$

$$Dh^{u}_{u}(0)=Dh^{u}_{w}(0)=0;\ |v|\ \text{sufficiently small}\Bigr{\}}$$

$$W^{c}_{loc}(0)=\Bigl{\{}(u,v,w)\in\mathbb{R}^{s}\times\mathbb{R}^{u}\times\mathbb{R}^{c}|u=h^{c}_{u}(w),v=h^{c}_{v}(w);$$

$$Dh^{c}_{u}(0)=Dh^{c}_{v}(0)=0;\ |w|\ \text{sufficiently small}\Bigr{\}}$$

$$\dot{x}=Ax+f(x,y),$$

$$\dot{y}=By+g(x,y),(x,y)\in {\mathbb{R}}^c\times {\mathbb{R}}^s,$$

$$f(0,0)=0,Df(0,0)=0,$$

$$g(0,0)=0,Dg(0,0)=0.$$

$$W^c(0)=\{(x,y)\in {\mathbb{R}}^c\times {\mathbb{R}}^s|y=h(x),|x|<\delta ,h(0)=0,Dh(0)=0\}$$

$$\dot{u}=Au+f(u,h(u)),u\in {\mathbb{R}}^c.$$

$$x(t)=u(t)+O\left(e^{-\gamma t}\right),$$

$$y(t)=h(u(t))+O\left(e^{-\gamma t}\right),$$

$$y=h(x).$$

$$\dot{y}=Dh(x)\dot{x}.$$

Full text

1 Introduction

The††footnotetext: Keywords: Michaelis-Menten kinetics, Center manifold, Tihonov’s Theorem, quasi-steady state approximation, singular perturbation. mathematical††footnotetext: 1 Dipartimento di Ingegneria Meccanica e Aerospaziale, Via Eudossiana n. 18, 00184 Roma. ††footnotetext: 2 Laboratorio di Strutture e Materiali Intelligenti - Sapienza University Palazzo Caetani, via San Pasquale snc, 04012 - Cisterna di Latina (LT) Italy. ††footnotetext: 3 Istituto di Analisi dei Sistemi ed Informatica “Antonio Ruberti” (IASI-CNR) Piazzale A. Moro, 7 00185 Rome - Italy. ††footnotetext: 4Dept. of Economics, University of Roma TRE, via Silvio D’Amico 77, 00145 Rome, Italy. study of enzyme kinetics (or Michaelis-Menten kinetics) is mainly based on the standard Quasi-Steady State Approximation (sQSSA) - which has been used in Biochemistry, since the pioneering papers by Bodenstein [7] and Chapman and Underhill [10] - which starts from the observation that enzyme reactions are characterized by a first, short transient phase, where the intermediate complex rapidly grows, and a second, longer quasi-equilibrium phase, where the complex slowly decays in the product, in general the activated substrate. Each phase of the reaction has a time scale ($t_{c}$ and $t_{s}$, respectively).

The quasi-steady state approximation is a very efficient way of simplification for the description of a typical saturation phenomenon, occurring in enzyme kinetics, but present in several other biological systems. Let us cite, just as non exhaustive examples, the Monod-Wyman-Changeux molecular model of cooperativity in allosteric reactions [27], or, more recently, the model of Lekszycki and coworkers concerning the bone regeneration [26, 17] and the model of infarcted cardiac tissue regeneration by means of stem cells [1], where saturation phenomena are observed.

The reaction can be described as follows.

Let us consider an enzyme, $E$, which reacts with a protein, $X$, resulting in an intermediate complex, $C$. In turn, this complex can break down into a product, $X_{p}$, and the enzyme $E$. It is frequently assumed that the formation of $C$ is reversible while its breakup is not. The process is represented by the following sequence of reactions

[TABLE]

where $k_{1},k_{-1},k_{2}$ are the reaction rates.

For notational convenience we will use variable names to denote both a chemical species and its concentration. For example, $E$ denotes both an enzyme and its concentration. Reaction (1) obeys two natural constrains: the total amounts of protein and enzyme remain constant. Therefore,

[TABLE]

for positive constants $X_{T}$ and $E_{T}$. In conjunction with the constraints (2), the following Cauchy Problem for a system of two ordinary differential equations can be used to model reaction (1):

[TABLE]

where $\dot{x}(t)=\frac{dx}{dt}$ and where $E_{T}$, $S_{T}$, $k_{1}$, $k_{2}$, $k_{-1}$ are viewed as fixed positive constants and $K_{M}=(k_{-1}+k_{2})/k_{1}$ is the Michaelis affinity constant. Similarly, we can define the dissociation constant $K_{D}=k_{-1}/k_{1}$ and the Van Slyke-Cullen constant $K=k_{2}/k_{1}$.

Since, after a short transient, where the complex $C$ rapidly grows, reaching its maximal concentration, it slowly decays, the sQSSA consists in supposing that, after the transient phase, the complex can be considered in a quasi-equilibrium, i.e., posing $\displaystyle\frac{dC}{dt}\cong 0$. With this approximation, the system becomes the differential-algebraic system

[TABLE]

($V_{max}=k_{2}E_{T}$) where only the initial condition $X(0)=X_{T}$ can be imposed, because the sQSSA describes only the slow phase, where the initial value of $C(t)$ is its maximal value, instead of [math].

In the Sixties of the last Century mathematicians (see, in particular, [20]) interpreted the sQSSA in terms of leading order term of asymptotic expansions with respect to a perturbation parameter $\varepsilon$, which must be supposed small. Heineken et al. used $\varepsilon_{HTA}=E_{T}/S_{T}$, because in literature it is widely used to impose that the initial concentration of the enzyme $E$ is much less than the concentration of the substrate $X$.

The parameter can also arise by virtue of a biochemical condition imposing the separation between the two timescales $t_{c}$ and $t_{s}$ characterizing the reaction (see also [25, 40, 41, 42, 32]). In this way, Segel-Slemrod [42] showed that the sQSSA can be obtained also as the leading order of an asymptotic expansion in terms of $\varepsilon_{SS}=\frac{E_{T}}{S_{T}+K_{M}}$, enlarging the parameter range of validity of the sQSSA.

Inspired by the papers by Laidler [23], Swoboda [44, 45], and Schauer and Heinrich [38], Borghans et al. [8] introduced a different approximation, called total Quasi-Steady State approximation (tQSSA) which uses the new variable $\overline{X}=X+C$, called total substrate.

Formally, introducing the ”lumped” variable $\bar{X}:=X+C$, problem (1) can be rewritten as

[TABLE]

Also the tQSSA posits that $C$ equilibrates quickly compared to $\bar{X}$.

Imposing also in this case a quasi-steady state approximation ($\displaystyle\frac{dC}{dt}\cong 0$), we obtain

[TABLE]

where

[TABLE]

is the only biologically allowed solution of $\displaystyle\frac{dC}{dt}=0$.

Let us remark that since, thanks to the conservations laws, $P(t)=X_{T}-\bar{X}(t)$, the tQSSA can be viewed as the other side of the coin of Laidler’s theory, though the approach followed in [38, 8] implicitly contains much more information about the reliability of the approximation, as shown in [2].

Also the tQSSA can be seen as the leading order term of an asymptotic expansion in terms of a suitable parameter, $\displaystyle\varepsilon=\frac{E_{T}}{(K_{M}+E_{T}+X_{T})^{2}}$, producing a new approximation, which is valid in a much wider parameter range. The parameter, introduced in [8], appears already in a paper by Palsson [32], where the author determines sufficient conditions for the application of any Quasi-Steady State Approximation, based again on the time scale separation. Taking into account that the perturbation parameter is always less than $1/4$, its introduction in terms of time scale separation appears much more natural than the previous parameters. This result gives a theoretical mathematical foundation of the choice of the parameter in the tQSSA. Moreover, several authors (see, for example, [40, 41]) study the transient phase of the reaction supposing that in this phase $X$ does not change substantially. This hypothesis is not realistic, while, using the total substrate $\bar{X}$, we observe that at time 0, we have $\dot{\bar{X}}(0)=0$, which addresses much more naturally the request of small changes of the total substrate in the initial time of the reaction.

In Figure 1 we show the different efficiency of the two quasi-steady state approximations, when the parameters are stressed in such a way that the sQSSA is no more valid.

In previous literature the different QSSAs are approached by means of two different tools: Tihonov’s Theorem [46, 47, 48], which studies the asymptotic stability of systems of differential equations characterized by the presence of small perturbative parameters and Center Manifold Theory, which is one of the most powerful tools to study the dimensional reduction of differential systems. For example, on the one hand, Heineken et al. [20] and Dvořák and Šiška [16] quote Tihonov’s Theorem in order to justify the sQSSA, while Khho and Hegland [21] refer to this theorem to apply the tQSSA; on the other hand, other authors [28, 22] interpret the sQSSA and the tQSSA, respectively, as the slow manifold of the Michaelis-Menten kinetics.

However, at the best of our knowledge, the two approaches are not yet compared, in order to check whether there exist any equivalences between the so-called singular points, introduced by Tihonov, and the center manifolds, as studied, for example, by Carr [9].

Applying the techniques exposed in [52, 53], we show that the two approximations are asymptotically equivalent, concluding that the sQSSA and the tQSSA can be interpreted both as the singular point of the Michaelis-Menten kinetics and its center manifold.

This means that Tihonov’s Theorem implies that any QSSA can be mathematically interpreted as the study of the reduced system of the original system setting the perturbative parameter $\varepsilon=0$, instead of setting the derivative of the complex $C(t)$ equal to zero.

This fact formally allows the application of a ”mechanistic” passage, consisting in equating to zero the derivatives of the complexes, in the single reaction scheme, as in more complex reactions, because this is the simplest way to reach (an approximate) expression of the center manifold. For example, Wang and Sontag [50] apply this technique for the study of the sQSSA of the double phosphorylation-dephosphorylation cycle.

As shown in [35], however, these approximation are no more applicable to mechanisms where oscillations can appear and an a priori analysis of the application should be performed every time we have to deal with any QSSA.

The paper is organized as follows in Section 2 we recall all the main definitions and properties of Tihonov’s Theory and of the Center Manifold Theory; in Section 3 we show the equivalence of the two approaches in the case of the sQSSA, of the tQSSA and for a class of more general systems of differential equations characterized by the presence of a small, perturbative parameter; in Section 4 we discuss some future applications of this results to more complex enzymatic reactions.

2 Preliminary results and notations on nonlinear dynamical systems.

In this section we introduce the notations and summarize the results we need for the formulation of the problem investigated here. For convenience of the reader we closely follow the notations of the fundamental book by Wiggins [53].

We will investigate systems in the class of general autonomous vector fields

[TABLE]

It is natural to consider the linearized system

[TABLE]

associated to the vector field (8), if $\bar{x}\in\mathbb{R}^{n}$ is one of its fixed points, and the constant Jacobian $n\times n$ matrix $A=Df(\bar{x})$. Then $\mathbb{R}^{n}$ can be represented as the direct sum of three subspaces denoted $E^{s}$, $E^{u}$, and $E^{c}$, which are defined as follows:

[TABLE]

where $\{e_{1},...,e_{s}\}$, $\{e_{s+1},...,e_{s+u}\}$, $\{e_{s+u+1},...,e_{s+u+c}\}$ are the (generalized) eigenvectors of $A$ corresponding to the eigenvalues having negative real part, positive real part and zero real part, respectively. $E^{s}$, $E^{u}$, and $E^{c}$ are referred to as the stable, unstable, and center subspaces, respectively. They are invariant subspaces (or manifolds) since solutions of (9) with initial conditions entirely contained in either $E^{s}$, $E^{u}$, or $E^{c}$ must forever remain in that particular subspace for all time.

It is well known that there exists a linear transformation $T$ which transforms the linear equation (9) into block diagonal form

[TABLE]

where $T^{-1}y\equiv(u,v,w)\in\mathbb{R}^{s}\times\mathbb{R}^{u}\times\mathbb{R}^{c}$, $s+u+c=n$, $A_{s}$ is an $s\times s$ matrix having eigenvalues with negative real part, $A_{u}$ is an $u\times u$ matrix having eigenvalues with positive real part, and $A_{c}$ is an $c\times c$ matrix having eigenvalues with zero real part. The [math] in the block diagonal form (11) indicate appropriately sized block consisting of all zero’s. Using this same linear transformation to transform the coordinates of the nonlinear vector field (8) gives the equation

[TABLE]

where $R_{s}(u,v,w)$, $R_{u}(u,v,w)$, and $R_{c}(u,v,w)$ are the first $s$, $u$, and $c$ components, respectively, of the vector $T^{-1}RTy$.

The following theorem shows how this structure changes when the nonlinear vector field (12) is considered. It is stated without proof (see [52] for details).

Theorem 2.1** (Local Stable, Unstable, and Center Manifolds of Fixed Points).**

Suppose (12) is $C^{r}$, $r\geq 2$. Then the fixed point $(u,v,w)=0$ of (12) possesses a $C^{r}$ s-dimensional local, invariant stable manifold, $W^{s}_{loc}(0)$, a $C^{r}$ u-dimensional local, invariant unstable manifold, $W^{u}_{loc}(0)$, and a $C^{r}$ c-dimensional local, invariant center manifold, $W^{c}_{loc}(0)$, all intersecting in $(u,v,w)=0$. These manifolds are all tangent to the respective invariant subspaces of the linear vector field (11) at the origin and, hence, are locally representable as graphs. In particular, we have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $h^{s}_{v}(u)$, $h^{s}_{w}(u)$, $h^{u}_{u}(v)$, $h^{u}_{w}(v)$, $h^{c}_{u}(w)$, and $h^{c}_{v}(w)$ are $C^{r}$ functions. Moreover, trajectories in $W^{s}_{loc}(0)$ and $W^{u}_{loc}(0)$ have the same asymptotic properties as trajectories in $E^{s}$ and $E^{u}$, respectively. Namely, trajectories of (12) with initial conditions in $W^{s}_{loc}(0)$ (resp., $W^{u}_{loc}(0)$) approach the origin at an exponential rate asymptotically as $t\rightarrow+\infty$ (resp., $t\rightarrow-\infty$).

If the eigenvalues of the center subspace are all precisely zero - rather than having just real part zero - then a center manifold is called a slow manifold.

2.1 Center Manifolds

If $E^{u}=\emptyset$, then any orbit will rapidly decay to $E^{c}$. Thus, in order to investigate the long-time behavior (i.e., stability) we need only to investigate the system restricted to $E^{c}$. This simple reasoning is the foundation of the “reduction principle” applied to the study of the stability of nonhyperbolic fixed points of nonlinear vector fields.

For our purposes, let us consider vector fields of the following form

[TABLE]

where

[TABLE]

In the above, $A$ is a $c\times c$ matrix having eigenvalues with zero real parts, $B$ is an $s\times s$ matrix having eigenvalues with negative real parts, and $f$ and $g$ are $C^{r}$ functions ($r\geq 2$).

For the sake of notation simplicity, let us write the center manifold in the following way:

[TABLE]

with $\delta$ sufficiently small.

Remark 2.2**.**

We remark that the conditions $h(0)=0$ and $Dh(0)=0$ imply that $W^{c}(0)$ is tangent to $E^{c}$ at $(x,y)=(0,0)$.

The following three theorems are taken from the seminal book [9], as reported in [53].

Theorem 2.3** (Existence).**

There exists a $C^{r}$ center manifold for (2.1). The dynamics of (2.1) restricted to the center manifold is, for $u$ sufficiently small, given by the following c-dimensional vector field

[TABLE]

The next result implies that the dynamics of (16) near $u=0$ determines the dynamics of (2.1) near $(x,y)=(0,0)$.

Theorem 2.4** (Stability).**

i) Let the zero solution of (16) be stable (asymptotically stable) (unstable); then the zero solution of (2.1) is also stable (asymptotically stable) (unstable). ii) Let the zero solution of (16) be stable. Then if $(x(t),y(t))$ is a solution of (2.1) with $(x(0),y(0))$ sufficiently small, there is a solution $u(t)$ of (16) such that as $t\rightarrow\infty$

[TABLE]

where $\gamma>0$ is a constant.

It is possible to compute the center manifold so that we can reap the benefits of Theorem 2.4. Using invariance of $W^{c}(0)$ under the dynamics of (2.1), we derive a quasi-linear partial differential equation that $h(x)$ must satisfy, in order for its graph to be a center manifold for (2.1). This is done as follows:

The $(x,y)$ coordinates of any point on $W^{c}(0)$ must satisfy

[TABLE] 2. 2.

Differentiating (18) with respect to time implies that the $(\dot{x},\dot{y})$ coordinates of any point on $W^{c}(0)$ must satisfy

[TABLE] 3. 3.

Any point on $W^{c}(0)$ obeys the dynamics generated by (2.1). Therefore, substituting

[TABLE]

into (19) gives

[TABLE]

or

[TABLE]

Then, to find a center manifold, all we need to do is to solve (22).

Theorem 2.5** (Approximation, [53]).**

Let $\phi:\mathbb{R}^{c}\rightarrow\mathbb{R}^{s}$ be a $C^{1}$ mapping with $\phi(0)=D\phi(0)=0$ such that $\mathcal{N}\left(\phi(x)\right)=O\left(|x|^{q}\right)$ as $x\rightarrow 0$ for some $q>1$. Then

[TABLE]

The theorem gives us a method for computing an approximate solution of (22) to any desired degree of accuracy. So, for this task, we will employ power series expansions (note that by Remark 2.2 power series expansions start from second order).

Suppose now that (2.1) depends on a vector of parameters $\varepsilon\in\mathbb{R}^{p}$:

[TABLE]

where

[TABLE]

with the same assumptions as in (2.1).

Following Wiggins [52, 53], we will handle parameterized systems including the parameter $\varepsilon$ as a new dependent variable as follows

[TABLE]

This system has a fixed point at $(x,y,\varepsilon)=(0,0,0)$. The matrix associated with the linearization of (2.1) around this fixed point has $s$ eigenvalues with negative real part and $c+p$ eigenvalues with zero real part. Let us now apply center manifold theory. Modifying definition given in formula (15), a center manifold will be represented as the graph of $h(x,\varepsilon)$ for $x$ and $\varepsilon$ sufficiently small. Theorem 2.3 still applies, with the vector field reduced to the center manifold given by

[TABLE]

Let us calculate the center manifold. Using invariance of the graph of $h(x,\varepsilon)$ under the dynamics generated by (2.1), we have

[TABLE]

Substituting (2.1) into (27) results in the following quasi-linear partial differential equation that $h(x,\varepsilon)$ must satisfy in order for its graph to be a center manifold.

[TABLE]

Although center manifolds exist, they do not need to be unique. This can be seen from a well-known example due to Anosov (see [43], [53]). It can be proven (see, among others, [9] as reported in [53]) that any couple of center manifolds of a given fixed point differ by (at most) exponentially small terms. Thus, the Taylor series expansions of any two center manifolds agree to all orders.

Moreover, it can be shown that, due to the attractive nature of the center manifold, certain orbits (for example, fixed points, periodic orbits, homoclinic orbits, and heteroclinic orbits) that remain close to the origin for all the time must be on every center manifold of a given fixed point.

2.2 A different viewpoint: Singular Perturbations

For this subsection we will refer to the widespread book by W. Wasow [51], and - in particular - to its relevant section on singular perturbations. A systematic study of the qualitative aspects of such singular perturbation problems can be found in a series of papers by Tihonov ([46], [47] and [48]).

We consider differential systems of the form

[TABLE]

where $x$ is $c$-dimensional vector and $y$ an $s$-dimensional vector. All variables are real, and $\varepsilon$ is positive.

We assume that:

• (A)

The functions $f$ and $g$ in (2.2) are continuous in an open region $\Omega$ of the $(x,y)$-space.

• (B)

There is an $s$-dimensional vector function $\phi(x)$ continuous in $\xi_{1}\leq x\leq\xi_{2}$ such that the points $(x,\phi(x))$, for all $\xi_{1}\leq x\leq\xi_{2}$, are in $\Omega$ and

[TABLE]

• (C)

There exists a number $\eta>0$, independent of $x$, such that the relations

[TABLE]

imply

[TABLE]

The function $\phi(x)$ will be referred to as a root of the equation $g(x,y)=0$. It is not excluded that $g(x,y)=0$ may have other roots besides $\phi(x)$. A root $\phi(x)$ that satisfies condition C will be called isolated in $\xi_{1}\leq x\leq\xi_{2}$.

Definition 2.6**.**

The system of differential equations

[TABLE]

in which $x$ is a parameter, will be called the boundary layer equation belonging to the system (2.2).

To (2.2) there corresponds the reduced (or degenerate) system

[TABLE]

The solutions of (2.2) and (2.6) define trajectories $\left(x(t,\varepsilon),y(t,\varepsilon)\right)$ and $\left(x_{0}(t),y_{0}(t)\right)$ in the $(x,y)$-space.

We also assume:

• (D)

The singular point $y=\phi(x)$ of the boundary layer equation (30) is asymptotically stable for all $\xi_{1}\leq x\leq\xi_{2}$.

The root $\phi(x)$ will be called, briefly, a stable root in $\xi_{1}\leq x\leq\xi_{2}$, if assumption (D) is satisfied.

In accordance with our previous terminology we refer to the problem consisting of equations (2.2) together with the initial condition

[TABLE]

as the full problem. The reduced problem is here defined by

[TABLE]

[TABLE]

The differential equation (2.2) is, of course, obtained by setting $\varepsilon=0$ in (2.2) and determining the root $y=\phi(x)$ of the equation $g(x,y)=0$. Moreover, we assume:

• (E)

The full problem, as well as the reduced one, has a unique solution in an interval $0\leq t\leq T$.

• (F)

The asymptotic stability of the singular point $y=\phi(x)$ is uniform with respect to $x$ in $\xi_{1}\leq x\leq\xi_{2}$.

Let $\mu>0$. The set of points in the $(x,y)$-space for which the inequalities

[TABLE]

hold will be called a “$\mu$-tube”. The set

[TABLE]

constitutes the “lateral boundary” of the $\mu$-tube.

Lemma 2.7**.**

Suppose assumptions (A) to (F) are satisfied. Let $\mu>0$ be arbitrary but so small that the closure of the $\mu$-tube lies in $\Omega$. There exist then two numbers $\gamma(\mu)$ and $\varepsilon(\mu)$ such that for $\varepsilon<\varepsilon(\mu)$ the following is true: Any solution of the full equation that is in the interior of the $\mu$-tube for some value $\tilde{t}$ of $t$, $0\leq\tilde{t}\leq T$, and in the closure of the $\mu$-tube for all $t$ in $\tilde{t}\leq t<T$, does not meet the lateral surface of the $\mu$-tube for $\tilde{t}\leq t<T$.

The lemma states that, for small $\varepsilon$, any solution that comes close to the curve $y=\phi(x)$ in $\xi_{1}\leq x\leq\xi_{2}$ remains close to it, as long as $\xi_{1}\leq x\leq\xi_{2}$.

For a convenient formulation of Tihonov’s Theorem, according to [51], we introduce one more term.

Definition 2.8**.**

A point $(\alpha,\beta)\in\Omega$, $\xi_{1}\leq\alpha\leq\xi_{2}$ is said to lie in the domain of influence of the stable root $y=\phi(x)$ if the solution of the problem

[TABLE]

exists and remains in $\Omega$ for all $\tau>0$, and if it tends to $\phi(\alpha)$, as $\tau\rightarrow+\infty$.

Theorem 2.9**.**

Let Assumptions (A) to (F) be satisfied and let $(\alpha,\beta)$ be a point in the domain of influence of the root $y=\phi(x)$. Then the solution $x(t,\varepsilon)$, $y(t,\varepsilon)$ of the full initial value problem (2.2), (32) is linked with the solution $(x_{0}(t)$, $y_{0}(t)=\phi(x_{0}(t)))$ of the reduced problem (2.2), (34) by the limiting relations

[TABLE]

Here $T_{0}$ is any number such that $y=\phi(x_{0}(t))$ is an isolated stable root of $g\left(x_{0}(t),y\right)=0$ for $0\leq t\leq T_{0}$. The convergence is uniform in $0\leq t\leq T_{0}$, for $x(t,\varepsilon)$, and in any interval $0<t_{1}\leq t\leq T_{0}$ for $y(t,\varepsilon)$.

Tihonov’s Theorem 2.9 is only the first step in the asymptotic solution of initial value problems of the singular perturbation type. The most natural approach to this problems is to attempt a solution (outer solution) in the form of a series in powers of $\varepsilon$:

[TABLE]

and to determine the coefficients $x_{r}(t)$, $y_{r}(t)$ by means of formal substitution and comparison of coefficients.

It is clear that we have to relate the series (36) to the behavior of the solution of (2.2) in the boundary layer, as shown, for example, in [20, 42, 11]. For values of $t$ that are of order $O(\varepsilon)$ the solution to our perturbation problem can be found starting from the stretching transformation $t=\tau\varepsilon$. Hence, the stretched form of the original problem is

[TABLE]

Also in this case we determine the solution of the transient phase (inner solution) in terms of a series in powers of $\varepsilon$.

The developments of the passages is beyond the scope of this paper. For the other accounts and a more detailed discussion, see [51].

3 Main results

3.1 The Heineken-Tsuchiya-Aris system

Let us consider the enzymatic reaction described in (1) and (1). Clearly $(X,C)=(0,0)$ is a fixed point of (1). Following [20], let us first adimensionalize equations (1). Let us observe that we could use different adimensionalization procedures, in particular using the parameter $\varepsilon_{SS}=\frac{E_{T}}{S_{T}+K_{M}}$, as proposed in [42]. However, we follow the simpler scheme shown in [20], just in order to test our theoretical results and compare them with the results shown in [9]:

[TABLE]

where

[TABLE]

and

[TABLE]

This is the Heineken-Tsuchiya-Aris system [20]. Carr ([9], p.8, example 3) uses equations

[TABLE]

To obtain (3.1) from (3.1) just impose $\kappa+\lambda=c$ and $\kappa=1$. We will start from (3.1) for having a more realistic system. By applying the sQSSA (which corresponds to impose $\varepsilon=0$), we have the reduced system (outer solution) of (3.1)

[TABLE]

As above remarked, Heineken et al. [20] and Dvořák and Šiška [16] quote Tihonov’s Theorem in order to justify the sQSSA.

The aim of this subsection is now to determine the center manifold for (3.1), using the techniques described in [52, 53] and to show that it is asymptotically equivalent to the singular points related to Tihonov theory.

To this aim, let us now set $\tau=\varepsilon s$. Equations (3.1) can be rewritten in the equivalent form (inner solution)

[TABLE]

where $\varphi(u,v)=-u+(u+\kappa-\lambda)v$. In order to obtain for (3.1) a block form, of the type (2.1), where the submatrix having eigenvalues with zero real parts is separated from the submatrix having eigenvalues with negative real parts, we operate the substitution $w:=u-\kappa v$, i.e., $v=\frac{u-w}{\kappa}$. Hence,

[TABLE]

and

[TABLE]

Following [53], the way we will handle parametrized systems consists of including the parameter $\varepsilon$ as a new dependent variable as in (3.1), which merely acts to augment the matrix $A$ by adding a new center direction that has no dynamics. In this way, system (3.1) becomes

[TABLE]

where the parameter $\varepsilon$ is a new variable and the system could have also other fixed points.

The associated linearized system has a diagonal form, where the eigenvalues are given by [math] (multiplicity $2$) and $-\kappa$.

To find a center manifold, all we need to do is to solve (22) for system (3.1), employing Theorem 2.5, which gives us a method for computing an approximate solution of (22) to any desired degree of accuracy. Referring to (22) and (2.1), $A=0$, $B=-\kappa$, so we search for a function $w=h(u,\varepsilon)$ such that

[TABLE]

where

[TABLE]

Using Theorem 2.5 we assume

[TABLE]

Substituting (46) into (44), one has:

[TABLE]

where

[TABLE]

Accordingly, substituting (3.1) into (3.1), one has

[TABLE]

Truncating at second order terms, we obtain:

[TABLE]

Equating terms of the same power to zero gives $a_{1}=\frac{1}{\kappa}$, $a_{2}=-\frac{\lambda}{\kappa^{2}}$ and $a_{3}=0$. Hence, the center manifold for system (3.1) is:

[TABLE]

which, for $\kappa=1$, gives the result shown in [9]. Finally, substituting (50) into (3.1), we obtain the vector field reduced to the center manifold, according to equation (16) of Theorem 2.3. In fact, if $a_{1}=\frac{1}{\kappa}$, $a_{2}=-\frac{\lambda}{\kappa^{2}}$ and $a_{3}=0$, formula (3.1) becomes:

[TABLE]

Thus:

[TABLE]

or, in terms of the original time scale,

[TABLE]

Let us now conclude showing that the center manifold obtained following this method is asymptotically sufficiently close to (3.1). We can obtain back the equation in $v$. In fact, since $v=\frac{u-w}{\kappa}$, from (50) and

[TABLE]

we have

[TABLE]

Considering $\varepsilon\ll 1$, one has

[TABLE]

which is the second equation of (3.1). We can conclude that, supposing $\varepsilon\ll 1$, the center manifold determined in this way approximates the solution given by the sQSSA, which coincides with the root related to the application of Tihonov’s Theorem.

In Figure (2) we compare the sQSSA of system (3.1), obtained from (3.1), with the center manifold (53), at the zeroth order and at the first order in $\varepsilon$, respectively. Obviously, the latter gives a better approximation of the numerical solution of (3.1), while the former well approximates the sQSSA curve, which in fact can be considered the zeroth order term of an asymptotic expansion of the solution in terms of $\varepsilon$.

3.2 The total quasi-steady state approximation (tQSSA)

Let us now consider system (1). Let us adimensionalize the system, as in [39, 11]

[TABLE]

where

[TABLE]

and

[TABLE]

with system parameters

[TABLE]

such that $\sigma+\eta+\kappa_{m}=1$.

By applying the tQSSA (which corresponds to imposing $\varepsilon=0$) to (3.2), we have

[TABLE]

and, solving in $v$:

[TABLE]

which represents the singular point (or outer solution) of (3.2), where $\eta$, $\sigma$, $\kappa_{m}$ are viewed as fixed positive constants and $\varepsilon$ is the parameter. Its fixed point is $(u,v)=(0,0)$.

The new parameter $\varepsilon$ appears already in [32] and was used in [39, 15, 11] to determine the asymptotic expansions whose leading term is just the tQSSA. In 2008 Khoo and Hegland [21] applied Tihonov’s Theorem [46, 47] in order to study the tQSSA, obtaining similar results as in [8].

The aim of this subsection is now to determine the center manifold for (3.2), using the techniques described in [52, 53] and to show that they are asymptotically equivalent to the singular points related to Tihonov theory.

To this aim, let us now set $\tau=\varepsilon s$, system (3.2) can be rewritten in the form (inner solution)

[TABLE]

In order to obtain for (57) a block form, of the type (2.1), we make the substitution $w:=u-(\eta+\kappa_{m})v$, i.e., $v=\frac{u-w}{\eta+\kappa_{m}}$. Hence,

[TABLE]

[TABLE]

[TABLE]

Doing so, and introducing the new variable $\varepsilon$, system (57) becomes

[TABLE]

The associated linearized system has a diagonal form and, in fact, the eigenvalues are given by [math] (with multiplicity $2$) and $-(\eta+\kappa_{m})$.

Also in this case the eigenvalue [math] has multiplicity $2$.

We solve (22) for system (58), employing Theorem 2.5 and determine the center manifold. Referring to (22) and (2.1), we have that $A=0$, $B=-(\eta+\kappa_{m})$. Accordingly, we search a function $h(u,\varepsilon)$ such that

[TABLE]

Using Theorem 2.5 we assume

[TABLE]

Substituting (60) into (3.2), one has:

[TABLE]

Truncating at second order terms, we obtain:

[TABLE]

so,

[TABLE]

[TABLE]

from which:

[TABLE]

Hence, the center manifold for system (58) is:

[TABLE]

Finally, substituting (62) into (58) we obtain the vector field reduced to the center manifold, according to equation (16) of Theorem 2.3. Then:

[TABLE]

or, in terms of the original time scale,

[TABLE]

Let us show that the center manifold obtained in (62) is asymptotically close to the root given by (56), in terms of Tihonov’s Theorem. From (62), and since $v=\frac{u-w}{\eta+\kappa_{m}}$, with

[TABLE]

we have

[TABLE]

Since equation (56) is obtained putting $\varepsilon\ll 1$, one has

[TABLE]

while, approximating the square root in (56), one has

[TABLE]

It follows that both (3.2) and (3.2) are asymptotic to $\frac{u}{\eta+\kappa_{m}}$ when $u\rightarrow 0$. This means that both the expressions can be intended as two different approximations of the center manifold.

In Figure (3) we compare the tQSSA of system (3.2), obtained from (56), with the center manifold (65), at the zeroth order and at the first order in $\varepsilon$, respectively. Obviously, the latter gives a better approximation of the numerical solution of (3.2), while the former well approximates the tQSSA curve, which in fact can be considered the zeroth order term of an asymptotic expansion of the solution in terms of $\varepsilon$.

3.3 A more general viewpoint

Let us consider now a more general system of the following form (outer solution)

[TABLE]

and the corresponding inner solution

[TABLE]

(with $\tau=\varepsilon s$) where

[TABLE]

The origin is a fixed point for (69). Heineken-Tsuchiya-Aris system (3.1) and the system obtained by the tQSSA approximation (57), are particular cases of the system (69)-(70).

We are able to state a more general theorem concerning the center manifold, which is the main result of our paper.

Let $w:=au+bv$; hence:

[TABLE]

[TABLE]

Doing so, system (69) becomes, for $a,b\in\mathbb{R}$ and $b<0$,

[TABLE]

The associated linearized system has a block form of type (2.1) and, in fact, the eigenvalues are given by [math] (with multiplicity $2$) and $b<0$. Thus in every system of the form (3.3) we are in presence of a center manifold.

We write equation (22) for system (3.3), employing Theorem 2.5. Referring to (22) and (2.1), we have that $A=0$, $B=b$. Accordingly, we search for a function $w=h(u,\varepsilon)$ such that

[TABLE]

[TABLE]

from which, since $D_{u}h(u,\varepsilon)\varepsilon\varphi\left(u,\frac{w-au}{b}\right)$ is a function at least of third order in $\varepsilon$ and $u$, while we are interested in a second order expression of function $h(u,\varepsilon)$, we can neglect this term and focus on

[TABLE]

Theorem 3.1**.**

The center manifold of (69) and the isolated point of (69) are asymptotically equivalent.

Proof.

Step 1.

Using Theorem 2.5 we assume

[TABLE]

and it is trivial to prove that $h(u,\varepsilon)$ satisfies (72) for $\lambda_{3}=0$. Moreover, from (70),

[TABLE]

where $\Theta(u^{2},uw,w^{2})$ contains the quadratic terms in $u$ and $w$.

Since the terms in $uw$ and $w^{2}$, with $w=h(u,\varepsilon)=\lambda_{1}u^{2}+\lambda_{2}u\varepsilon+\dots$, are at least of third order in $\varepsilon$ and $u$, we consider only term in $u^{2}$. Therefore,

[TABLE]

while for $\varphi$ it is sufficient to consider the first order expansion in $u$ because, otherwise, in (72) we would have third order terms for $\varepsilon\varphi\left(u,\frac{w-au}{b}\right)$ in $\varepsilon$ and $u$. Thus,

[TABLE]

where we recall that $v=\frac{w-au}{b}$. Accordingly, equation (72) becomes:

[TABLE]

Equating to zero terms of the same power gives

[TABLE]

Hence, the center manifold for system (69) is

[TABLE]

Setting in the RHS $\varepsilon=0$, we obtain the center manifold $w=h(u,0)$ of (69).

Step 2. Singular Point Technique [51]

On the other hand,

[TABLE]

Since, setting $\varepsilon=0$, we have that $v=\frac{w-au}{b}=\frac{\lambda_{1}u^{2}-au}{b}$, equation (72) becomes

[TABLE]

which gives an identity up to $O(u^{2})$, if we substitute $\lambda_{1}$ as above and if we operate a Taylor expansion around $(u,v)=(0,0)$.

The vector field reduced to the center manifold, from equation (16) of Theorem 2.3, is:

[TABLE]

or, in terms of the original time scale,

[TABLE]

Moreover:

[TABLE]

where $v=\frac{w-au}{b}$ and all the derivatives are calculated in $(0,0)$. Hence,

[TABLE]

and, since $\lambda_{2}=-\frac{a}{b}\left[\varphi_{u}(0,0)-\frac{a}{b}\varphi_{v}(0,0)\right]$, we have:

[TABLE]

for $v=\frac{w-au}{b}$. From (73), the vector field reduced to the center manifold, in terms of the original time scale, near the origin, becomes:

[TABLE]

for values of $\lambda_{1}$ and $\lambda_{2}$ as above.

Summarizing, we have obtained two relations:

a) From the Center Manifold Theory: considering (72), for $w=h(u,\varepsilon)$, and setting $\varepsilon=0$, we have:

[TABLE]

b) From Singular Perturbation Techniques: by assumption B of section 2.2, and since $g(u,w)=bw+a\varepsilon\varphi\left(u,\frac{w-au}{b}\right)+b\psi\left(u,\frac{w-au}{b}\right)$ in (3.3), we have, putting $\varepsilon=0$:

[TABLE]

Comparing (81) and (82), we observe a relation between $h(u,0)$ and $\phi(u)$ but we cannot infer that $h(u,0)=\phi(u)$, due to the non-uniqueness of center manifold. However, in the above steps we have proven that

[TABLE]

Q.E.D. ∎

This theorem means that the center manifolds obtained by means of (81) and (82) are asymptotically equivalent. This allows us to interpret any QSSA, obtained imposing $\dot{C}=0$, as a manifold which is asymptotically equivalent to the center manifold, as confirmed by equation (54) for Heineken-Tsuchiya-Aris system, and by equations (3.2)-(3.2) for the tQSSA.

This explains why, in order to achieve the center manifold of (3.1) and (3.2), it is sufficient to consider - for $u\rightarrow 0$ - the expression obtained equating to zero the second equation of these systems, (i.e. for $\varepsilon=0$). We recall that in many papers (see, for example, [20, 16, 21], who refer to Tihonov’s Theorem, and [22], according to Singular Perturbation Theory), the center manifold is obtained equating to zero the right hand side of the equation of the form:

[TABLE]

4 Conclusions and Perspectives

The quasi-steady state approximation has been a challenge for applied mathematicians, who had to explain the feasibility of an approximation which imposes to the complex $C$ both to be constant and to depend on $X$. Some Biochemistry texts (see, for example, [24, 54, 37, 19]) mislead the reader, interpreting the QSSA as a true equality, which brings to assert that the ratio $\displaystyle{E(t)\,S(t)}/{C(t)}$ is constant during all the quasi-steady state phase. This is obviously not true. In [5] the authors solve the apparent incongruence, determining the asymptotic value of $\displaystyle{E(t)\,S(t)}/{C(t)}$, showing that, for every choice of the kinetic parameters and of the initial conditions,

[TABLE]

(where $\alpha=\frac{k_{1}}{2}(K_{M}+E_{T})\left[1-\sqrt{1-\frac{4k_{2}E_{T}}{k_{1}(K_{M}+E_{T})^{2}}}\;\right]$ ), differently from what is wrongly stated.

Heineken et al. [20] and successively other authors [42, 39, 11] interpreted sQSSA and tQSSA as leading order expansions of the solutions in terms of a suitable parameter, which has to be considered small.

This interpretation allows us to embed the QSSA theory in a framework which is related to Tihonov’s Theorem [46, 47, 48, 20, 42, 51, 16, 21], where the parameter multiplies the derivative of $C$ and the QSSA can be obtained as the singular point of the original system, setting $\varepsilon=0$.

In this paper we have shown that, at least in the classical simple scheme (1), the approximation obtained applying Tihonov’s Theorem is asymptotically equivalent to the center manifold of the system, which means that reduced system and center manifold are two sides of the same coin.

Once again, the total QSSA has shown to be much more efficient and natural than the standard one, mainly thanks to the fact that the parameter used for the expansions in the total framework is always less than $\displaystyle\frac{1}{4}$.

In our actual researches we are applying the techniques shown in this paper to more complex enzyme reactions, as the fully competitive inhibition [4], the phosphorylation-dephosphorylation cycle (or Goldbeter-Koshland switch [18]), the linear double phosphorylation reaction, the double phosphorylation-dephosphorylation cycle [50] and, more in general, futile cycles [49].

These mechanisms were already studied in terms of tQSSA in previous papers [33, 3, 34, 35, 36, 12, 13, 6].

The techniques here shown will allow to read the tQSSA as the leading term of an asymptotic expansion in terms of a suitable perturbation parameter, in these more complex cases, too.

Acknowledgements

The authors are deeply grateful to Prof. Enzo Orsingher, from Sapienza University (Rome, Italy) and Prof. Jan Andres, from Palacky University (Olomouc, Czech Republic) for their translations of the papers [46, 47, 48] and some precious clarifications concerning some passages of the papers.