Dynamical aspects of the total QSSA in Enzyme Kinetics
Alberto Maria Bersani, Enrico Bersani, Alessandro Borri, Pierluigi Vellucci
DDynamical aspects of the total QSSA in EnzymeKinetics.
Alberto Maria Bersani Enrico Bersani Alessandro Borri Pierluigi Vellucci [email protected] 10, 2018 Abstract
In this paper we prove that the well-known quasi-steady state ap-proximations, commonly used in enzyme kinetics, which can be inter-preted as the reduced system of a differential system depending on aperturbative parameter, according to Tihonov theory, are asymptot-ically equivalent to the center manifold of the system. This allowsto give a mathematical foundation for the application of a mechanis-tic method to determine the center manifold of (at this moment, stillsimple) enzyme reactions.
The mathematical study of enzyme kinetics (or Michaelis-Menten kinet-ics) is mainly based on the standard Quasi-Steady State Approximation(sQSSA) - which has been used in Biochemistry, since the pioneering pa-pers by Bodenstein [7] and Chapman and Underhill [10] - which starts fromthe observation that enzyme reactions are characterized by a first, shorttransient phase, where the intermediate complex rapidly grows, and a sec-ond, longer quasi-equilibrium phase, where the complex slowly decays in theproduct, in general the activated substrate. Each phase of the reaction hasa time scale ( t c and t s , respectively). Keywords: Michaelis-Menten kinetics, Center manifold, Tihonov’s Theorem, quasi-steady state approximation, singular perturbation. Dipartimento di Ingegneria Meccanica e Aerospaziale, Via Eudossiana n. 18, 00184Roma. Laboratorio di Strutture e Materiali Intelligenti - Sapienza University Palazzo Cae-tani, via San Pasquale snc, 04012 - Cisterna di Latina (LT) Italy. Istituto di Analisi dei Sistemi ed Informatica “Antonio Ruberti” (IASI-CNR) Piaz-zale A. Moro, 7 00185 Rome - Italy. Dept. of Economics, University of Roma TRE, via Silvio D’Amico 77, 00145 Rome,Italy. a r X i v : . [ m a t h . D S ] F e b he quasi-steady state approximation is a very efficient way of simpli-fication for the description of a typical saturation phenomenon, occurringin enzyme kinetics, but present in several other biological systems. Let uscite, just as non exhaustive examples, the Monod-Wyman-Changeux molec-ular model of cooperativity in allosteric reactions [27], or, more recently, themodel 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 , resultingin an intermediate complex, C . In turn, this complex can break down into aproduct, X p , and the enzyme E . It is frequently assumed that the formationof C is reversible while its breakup is not. The process is represented by thefollowing sequence of reactions X + E k −− (cid:42)(cid:41) −− k − C k −→ X p + E (1)where k , k − , k are the reaction rates.For notational convenience we will use variable names to denote botha chemical species and its concentration. For example, E denotes both anenzyme and its concentration. Reaction (1) obeys two natural constrains:the total amounts of protein and enzyme remain constant. Therefore, X + C + X p = X T and E + C = E T , (2)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 equa-tions can be used to model reaction (1):˙ X ( t ) = − k X ( t ) ( E T − C ( t )) + k − C ( t )˙ C ( t ) = k X ( t ) ( E T − C ( t )) − ( k − + k ) C ( t )= k [ X ( t ) ( E T − C ( t )) − K M C ( t )] X (0) = X T , C (0) = 0 . (3)where ˙ x ( t ) = dxdt and where E T , S T , k , k , k − are viewed as fixed positiveconstants and K M = ( k − + k ) /k is the Michaelis affinity constant. Sim-ilarly, we can define the dissociation constant K D = k − /k and the VanSlyke-Cullen constant K = k /k .Since, after a short transient, where the complex C rapidly grows, reach-ing its maximal concentration, it slowly decays, the sQSSA consists in sup-posing that, after the transient phase, the complex can be considered in aquasi-equilibrium, i.e., posing dCdt ∼ = 0. With this approximation, the systembecomes the differential-algebraic system2 ( t ) = E T X ( t ) X ( t ) + K M , ˙ X ( t ) = − k C ( t ) = − V max X ( t ) X ( t ) + K M , X (0) = X T , (4)( V max = k 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 0.In the Sixties of the last Century mathematicians (see, in particular,[20]) interpreted the sQSSA in terms of leading order term of asymptoticexpansions with respect to a perturbation parameter ε , which must be sup-posed small. Heineken et al. used ε HT A = E T /S T , because in literature itis widely used to impose that the initial concentration of the enzyme E ismuch less than the concentration of the substrate X .The parameter can also arise by virtue of a biochemical condition im-posing the separation between the two timescales t c and t s characterizingthe 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 anasymptotic expansion in terms of ε SS = E T S T + K M , enlarging the parameterrange of validity of the sQSSA.Inspired by the papers by Laidler [23], Swoboda [44, 45], and Schauerand Heinrich [38], Borghans et al. [8] introduced a different approximation,called total Quasi-Steady State approximation (tQSSA) which uses the newvariable X = X + C , called total substrate.Formally, introducing the ”lumped” variable ¯ X := X + C , problem (3)can be rewritten as˙¯ X ( t ) = − k C ( t ) , ˙ C ( t ) = k (cid:2) ¯ X ( t ) E T − (cid:0) ¯ X ( t ) + E T + K M (cid:1) C ( t ) + C ( t ) (cid:3) ,X (0) = X T , C (0) = 0 . (5)Also the tQSSA posits that C equilibrates quickly compared to ¯ X .Imposing also in this case a quasi-steady state approximation ( dCdt ∼ = 0),we obtain ˙¯ X ≈ − k C − ( ¯ X ) , ¯ X (0) = X T , (6)where C − ( ¯ X ) = ( E T + K M + ¯ X ) − (cid:112) ( E T + K M + ¯ X ) − E T ¯ X dCdt = 0.3et us remark that since, thanks to the conservations laws, P ( t ) = X T − ¯ X ( t ), the tQSSA can be viewed as the other side of the coin of Laidler’stheory, though the approach followed in [38, 8] implicitly contains muchmore information about the reliability of the approximation, as shown in[2]. Also the tQSSA can be seen as the leading order term of an asymptoticexpansion in terms of a suitable parameter, ε = E T ( K M + E T + X T ) , produc-ing 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 ofany Quasi-Steady State Approximation, based again on the time scale sep-aration. Taking into account that the perturbation parameter is always lessthan 1 /
4, its introduction in terms of time scale separation appears muchmore natural than the previous parameters. This result gives a theoreti-cal mathematical foundation of the choice of the parameter in the tQSSA.Moreover, several authors (see, for example, [40, 41]) study the transientphase of the reaction supposing that in this phase X does not change sub-stantially. This hypothesis is not realistic, while, using the total substrate¯ X , we observe that at time 0, we have ˙¯ X (0) = 0, which addresses muchmore naturally the request of small changes of the total substrate in theinitial time of the reaction.In Figure 1 we show the different efficiency of the two quasi-steady stateapproximations, when the parameters are stressed in such a way that thesQSSA is no more valid.Figure 1: Comparison of the complexes (left) and of the substrates (right),solution of the system (3), with their sQSSA (4) and tQSSA (6). Theparameter set is the following: k = k = 1; k − = 4; E T = 89; X T =100; K M = 5; K = 4; ε SS = 0 . ε = 0 .
01. The inadequacy of the sQSSA,mainly in the first part of the reaction, is evident, while the tQSSA inindistinguishable from the numerical solution of the system.4n previous literature the different QSSAs are approached by means oftwo different tools: Tihonov’s Theorem [46, 47, 48], which studies the asymp-totic stability of systems of differential equations characterized by the pres-ence of small perturbative parameters and Center Manifold Theory, which isone of the most powerful tools to study the dimensional reduction of differen-tial systems. For example, on the one hand, Heineken et al. [20] and Dvoˇr´akand ˇSiˇska [16] quote Tihonov’s Theorem in order to justify the sQSSA, whileKhho and Hegland [21] refer to this theorem to apply the tQSSA; on theother 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 yetcompared, in order to check whether there exist any equivalences betweenthe so-called singular points, introduced by Tihonov, and the center mani-folds, as studied, for example, by Carr [9].Applying the techniques exposed in [52, 53], we show that the two ap-proximations are asymptotically equivalent, concluding that the sQSSA andthe 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 bemathematically interpreted as the study of the reduced system of the origi-nal system setting the perturbative parameter ε = 0, instead of setting thederivative 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 singlereaction scheme, as in more complex reactions, because this is the simplestway to reach (an approximate) expression of the center manifold. For exam-ple, Wang and Sontag [50] apply this technique for the study of the sQSSAof the double phosphorylation-dephosphorylation cycle.As shown in [35], however, these approximation are no more applicableto mechanisms where oscillations can appear and an a priori analysis of theapplication 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 maindefinitions and properties of Tihonov’s Theory and of the Center ManifoldTheory; in Section 3 we show the equivalence of the two approaches in thecase of the sQSSA, of the tQSSA and for a class of more general systems ofdifferential equations characterized by the presence of a small, perturbativeparameter; in Section 4 we discuss some future applications of this resultsto more complex enzymatic reactions.5 Preliminary results and notations on nonlineardynamical systems.
In this section we introduce the notations and summarize the results weneed for the formulation of the problem investigated here. For convenienceof the reader we closely follow the notations of the fundamental book byWiggins [53].We will investigate systems in the class of general autonomous vectorfields ˙ x = f ( x ) , x ∈ R n . (8)It is natural to consider the linearized system˙ y = Ay, y ∈ R n , (9)associated to the vector field (8), if ¯ x ∈ R n is one of its fixed points, and theconstant Jacobian n × n matrix A = Df (¯ x ). Then R n can be representedas the direct sum of three subspaces denoted E s , E u , and E c , which aredefined as follows: E s = span { e , ..., e s } ,E u = span { e s +1 , ..., e s + u } , s + u + c = nE c = span { e s + u +1 , ..., e s + u + c } , (10)where { e , ..., 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 re-ferred to as the stable, unstable, and center subspaces, respectively. Theyare invariant subspaces (or manifolds) since solutions of (9) with initial con-ditions entirely contained in either E s , E u , or E c must forever remain inthat particular subspace for all time.It is well known that there exists a linear transformation T which trans-forms the linear equation (9) into block diagonal form ˙ u ˙ v ˙ w = A s A u
00 0 A c uvw , (11)where T − y ≡ ( u, v, w ) ∈ R s × R u × R c , s + u + c = n , A s is an s × s matrix having eigenvalues with negative real part, A u is an u × u matrixhaving eigenvalues with positive real part, and A c is an c × c matrix havingeigenvalues with zero real part. The 0 in the block diagonal form (11)indicate appropriately sized block consisting of all zero’s. Using this samelinear transformation to transform the coordinates of the nonlinear vector6eld (8) gives the equation ˙ u = A s u + R s ( u, v, w ) , ˙ v = A u v + R u ( u, v, w ) , ˙ w = A c w + R c ( u, v, w ) , (12)where R s ( u, v, w ), R u ( u, v, w ), and R c ( u, v, w ) are the first s , u , and c com-ponents, respectively, of the vector T − RT y .The following theorem shows how this structure changes when the non-linear vector field (12) is considered. It is stated without proof (see [52] fordetails).
Theorem 2.1 (Local Stable, Unstable, and Center Manifolds of FixedPoints) . Suppose (12) is C r , r ≥ . Then the fixed point ( u, v, w ) = 0 of(12) possesses a C r s-dimensional local, invariant stable manifold, W sloc (0) ,a C r u-dimensional local, invariant unstable manifold, W uloc (0) , and a C r c-dimensional local, invariant center manifold, W cloc (0) , all intersecting in ( u, v, w ) = 0 . These manifolds are all tangent to the respective invariantsubspaces of the linear vector field (11) at the origin and, hence, are locallyrepresentable as graphs. In particular, we have W sloc (0) = (cid:110) ( u, v, w ) ∈ R s × R u × R c | v = h sv ( u ) , w = h sw ( u ); Dh sv (0) = Dh sw (0) = 0; | u | sufficiently small (cid:111) W uloc (0) = (cid:110) ( u, v, w ) ∈ R s × R u × R c | u = h uu ( v ) , w = h uw ( v ); Dh uu (0) = Dh uw (0) = 0; | v | sufficiently small (cid:111) W cloc (0) = (cid:110) ( u, v, w ) ∈ R s × R u × R c | u = h cu ( w ) , v = h cv ( w ); Dh cu (0) = Dh cv (0) = 0; | w | sufficiently small (cid:111) where h sv ( u ) , h sw ( u ) , h uu ( v ) , h uw ( v ) , h cu ( w ) , and h cv ( w ) are C r functions.Moreover, trajectories in W sloc (0) and W uloc (0) have the same asymptoticproperties as trajectories in E s and E u , respectively. Namely, trajectories of(12) with initial conditions in W sloc (0) (resp., W uloc (0) ) approach the originat an exponential rate asymptotically as t → + ∞ (resp., t → −∞ ). If the eigenvalues of the center subspace are all precisely zero - ratherthan having just real part zero - then a center manifold is called a slowmanifold . 7 .1 Center Manifolds If E u = ∅ , then any orbit will rapidly decay to E c . Thus, in order toinvestigate the long-time behavior (i.e., stability) we need only to investigatethe system restricted to E c . This simple reasoning is the foundation of the“reduction principle” applied to the study of the stability of nonhyperbolicfixed points of nonlinear vector fields.For our purposes, let us consider vector fields of the following form˙ x = Ax + f ( x, y ) , ˙ y = By + g ( x, y ) , ( x, y ) ∈ R c × R s , (13)where f (0 ,
0) = 0 , Df (0 ,
0) = 0 ,g (0 ,
0) = 0 , Dg (0 ,
0) = 0 . (14)In the above, A is a c × c matrix having eigenvalues with zero real parts, B is an s × s matrix having eigenvalues with negative real parts, and f and g are C r functions ( r ≥ W c (0) = { ( x, y ) ∈ R c × R s | y = h ( x ) , | x | < δ, h (0) = 0 , Dh (0) = 0 } (15)with δ sufficiently small. Remark 2.2.
We remark that the conditions h (0) = 0 and Dh (0) = 0 implythat W c (0) is tangent to E c at ( x, y ) = (0 , . The following three theorems are taken from the seminal book [9], asreported in [53].
Theorem 2.3 (Existence) . There exists a C r center manifold for (13). Thedynamics of (13) restricted to the center manifold is, for u sufficiently small,given by the following c-dimensional vector field ˙ u = Au + f ( u, h ( u )) , u ∈ R c . (16)The next result implies that the dynamics of (16) near u = 0 determinesthe dynamics of (13) near ( x, y ) = (0 , Theorem 2.4 (Stability) . i) Let the zero solution of (16) be stable (asymp-totically stable) (unstable); then the zero solution of (13) 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 (13) with ( x (0) , y (0)) sufficiently small,there is a solution u ( t ) of (16) such that as t → ∞ x ( t ) = u ( t ) + O (cid:0) e − γt (cid:1) ,y ( t ) = h ( u ( t )) + O (cid:0) e − γt (cid:1) , (17) where γ > is a constant.
8t is possible to compute the center manifold so that we can reap thebenefits of Theorem 2.4. Using invariance of W c (0) under the dynamics of(13), we derive a quasi-linear partial differential equation that h ( x ) mustsatisfy, in order for its graph to be a center manifold for (13). This is doneas follows:1. The ( x, y ) coordinates of any point on W c (0) must satisfy y = h ( x ) . (18)2. Differentiating (18) with respect to time implies that the ( ˙ x, ˙ y ) coor-dinates of any point on W c (0) must satisfy˙ y = Dh ( x ) ˙ x. (19)3. Any point on W c (0) obeys the dynamics generated by (13). Therefore,substituting ˙ x = Ax + f ( x, h ( x )) , ˙ y = Bh ( x ) + g ( x, h ( x )) , ( x, y ) ∈ R c × R s , (20)into (19) gives Bh ( x ) + g ( x, h ( x )) = Dh ( x ) ( Ax + f ( x, h ( x ))) . (21)or N ( h ( x )) ≡ Dh ( x ) ( Ax + f ( x, h ( x ))) − Bh ( x ) − g ( x, h ( x )) = 0 (22)Then, to find a center manifold, all we need to do is to solve (22). Theorem 2.5 (Approximation, [53]) . Let φ : R c → R s be a C mappingwith φ (0) = Dφ (0) = 0 such that N ( φ ( x )) = O ( | x | q ) as x → for some q > . Then | h ( x ) − φ ( x ) | = O ( | x | q ) , as x → . The theorem gives us a method for computing an approximate solutionof (22) to any desired degree of accuracy. So, for this task, we will employpower series expansions (note that by Remark 2.2 power series expansionsstart from second order).Suppose now that (13) depends on a vector of parameters ε ∈ R p :˙ x = Ax + f ( x, y, ε ) , ˙ y = By + g ( x, y, ε ) , ( x, y, ε ) ∈ R c × R s × R p , (23)9here f (0 , ,
0) = 0 , Df (0 , ,
0) = 0 ,g (0 , ,
0) = 0 , Dg (0 , ,
0) = 0 . (24)with the same assumptions as in (13).Following Wiggins [52, 53], we will handle parameterized systems includ-ing the parameter ε as a new dependent variable as follows˙ x = Ax + f ( x, y, ε ) , ˙ ε = 0 , ˙ y = By + g ( x, y, ε ) , ( x, y, ε ) ∈ R c × R s × R p , (25)This system has a fixed point at ( x, y, ε ) = (0 , , s eigenvalues withnegative real part and c + p eigenvalues with zero real part. Let us now applycenter manifold theory. Modifying definition given in formula (15), a centermanifold will be represented as the graph of h ( x, ε ) for x and ε sufficientlysmall. Theorem 2.3 still applies, with the vector field reduced to the centermanifold given by˙ u = Au + f ( u, h ( u, ε ) , ε ) , ˙ ε = 0 , ( u, ε ) ∈ R c × R p . (26)Let us calculate the center manifold. Using invariance of the graph of h ( x, ε )under the dynamics generated by (25), we have˙ y = D x h ( x, ε ) ˙ x + D ε h ( x, ε ) ˙ ε = Bh ( x, ε ) + g ( x, h ( x, ε ) , ε ) . (27)Substituting (25) into (27) results in the following quasi-linear partialdifferential equation that h ( x, ε ) must satisfy in order for its graph to be acenter manifold. N ( h ( x, ε )) ≡ D x h ( x, ε ) [ Ax + f ( x, h ( x, ε ) , ε )] + − Bh ( x, ε ) − g ( x, h ( x, ε ) , ε ) = 0 (28)Although center manifolds exist, they do not need to be unique. Thiscan be seen from a well-known example due to Anosov (see [43], [53]). Itcan be proven (see, among others, [9] as reported in [53]) that any coupleof center manifolds of a given fixed point differ by (at most) exponentiallysmall terms. Thus, the Taylor series expansions of any two center manifoldsagree to all orders.Moreover, it can be shown that, due to the attractive nature of thecenter manifold, certain orbits (for example, fixed points, periodic orbits,homoclinic orbits, and heteroclinic orbits) that remain close to the originfor all the time must be on every center manifold of a given fixed point.10 .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. Asystematic study of the qualitative aspects of such singular perturbationproblems can be found in a series of papers by Tihonov ([46], [47] and [48]).We consider differential systems of the form dxdt = f ( x, y ) ε dydt = g ( x, y ) , (29)where x is c -dimensional vector and y an s -dimensional vector. All variablesare real, and ε is positive.We assume that:(A) The functions f and g in (29) are continuous in an open region Ω ofthe ( x, y )-space.(B) There is an s -dimensional vector function φ ( x ) continuous in ξ ≤ x ≤ ξ such that the points ( x, φ ( x )), for all ξ ≤ x ≤ ξ , are in Ω and g ( x, φ ( x )) ≡ . (C) There exists a number η >
0, independent of x , such that the relations (cid:107) y − φ ( x ) (cid:107) < η, y (cid:54) = φ ( x ) in ξ ≤ x ≤ ξ imply g ( x, y ) (cid:54) = 0 , in ξ ≤ x ≤ ξ . The function φ ( 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 φ ( x ). Aroot φ ( x ) that satisfies condition C will be called isolated in ξ ≤ x ≤ ξ . Definition 2.6.
The system of differential equations ε dydt = g ( x, y ) (30) in which x is a parameter, will be called the boundary layer equation belong-ing to the system (29).To (29) there corresponds the reduced (or degenerate) system dx dt = f ( x , y )0 = g ( x , y ) . (31) The solutions of (29) and (31) define trajectories ( x ( t, ε ) , y ( t, ε )) and ( x ( t ) , y ( t )) in the ( x, y ) -space.We also assume: y = φ ( x ) of the boundary layer equation (30) isasymptotically stable for all ξ ≤ x ≤ ξ .The root φ ( x ) will be called, briefly, a stable root in ξ ≤ x ≤ ξ , if assump-tion (D) is satisfied.In accordance with our previous terminology we refer to the problemconsisting of equations (29) together with the initial condition x = α, y = β, for t = 0 (32)as the full problem. The reduced problem is here defined by dxdt = f ( x, φ ( x )) y = φ ( x ) , (33) x = α, for t = 0 (34)The differential equation (33) is, of course, obtained by setting ε = 0 in (29)and determining the root y = φ ( 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 inan interval 0 ≤ t ≤ T .(F) The asymptotic stability of the singular point y = φ ( x ) is uniformwith respect to x in ξ ≤ x ≤ ξ .Let µ >
0. The set of points in the ( x, y )-space for which the inequalities (cid:107) y − φ ( x ) (cid:107) < µ, ξ ≤ x ≤ ξ hold will be called a “ µ -tube”. The set (cid:107) y − φ ( x ) (cid:107) = µ, ξ ≤ x ≤ ξ constitutes the “lateral boundary” of the µ -tube. Lemma 2.7.
Suppose assumptions (A) to (F) are satisfied. Let µ > bearbitrary but so small that the closure of the µ -tube lies in Ω . There existthen two numbers γ ( µ ) and ε ( µ ) such that for ε < ε ( µ ) the following is true:Any solution of the full equation that is in the interior of the µ -tube forsome value ˜ t of t , ≤ ˜ t ≤ T , and in the closure of the µ -tube for all t in ˜ t ≤ t < T , does not meet the lateral surface of the µ -tube for ˜ t ≤ t < T . The lemma states that, for small ε , any solution that comes close to thecurve y = φ ( x ) in ξ ≤ x ≤ ξ remains close to it, as long as ξ ≤ x ≤ ξ .For a convenient formulation of Tihonov’s Theorem, according to [51],we introduce one more term. 12 efinition 2.8. A point ( α, β ) ∈ Ω , ξ ≤ α ≤ ξ is said to lie in the domainof influence of the stable root y = φ ( x ) if the solution of the problem dy/dτ = g ( α, y ) , y (0) = β exists and remains in Ω for all τ > , and if it tends to φ ( α ) , as τ → + ∞ . Theorem 2.9.
Let Assumptions (A) to (F) be satisfied and let ( α, β ) be apoint in the domain of influence of the root y = φ ( x ) . Then the solution x ( t, ε ) , y ( t, ε ) of the full initial value problem (29), (32) is linked with thesolution ( x ( t ) , y ( t ) = φ ( x ( t ))) of the reduced problem (33), (34) by thelimiting relations lim ε → x ( t, ε ) = x ( t ) , ≤ t ≤ T lim ε → y ( t, ε ) = y ( t ) = φ ( x ( t )) 0 < t ≤ T (35) Here T is any number such that y = φ ( x ( t )) is an isolated stable root of g ( x ( t ) , y ) = 0 for ≤ t ≤ T . The convergence is uniform in ≤ t ≤ T ,for x ( t, ε ) , and in any interval < t ≤ t ≤ T for y ( t, ε ) . Tihonov’s Theorem 2.9 is only the first step in the asymptotic solution ofinitial value problems of the singular perturbation type. The most naturalapproach to this problems is to attempt a solution ( outer solution ) in theform of a series in powers of ε : x = ∞ (cid:88) r =0 x r ( t ) ε r , y = ∞ (cid:88) r =0 y r ( t ) ε r (36)and to determine the coefficients x r ( t ), y r ( t ) by means of formal substitutionand comparison of coefficients.It is clear that we have to relate the series (36) to the behavior of thesolution of (29) in the boundary layer, as shown, for example, in [20, 42,11]. For values of t that are of order O ( ε ) the solution to our perturbationproblem can be found starting from the stretching transformation t = τ ε .Hence, the stretched form of the original problem is dxdτ = εf ( x, y ) , dydτ = g ( x, y ) ,x = α, y = β, for τ = 0 . (37)Also in this case we determine the solution of the transient phase ( innersolution ) in terms of a series in powers of ε .The developments of the passages is beyond the scope of this paper. Forthe other accounts and a more detailed discussion, see [51].13 Main results
Let us consider the enzymatic reaction described in (1) and (3). Clearly(
X, C ) = (0 ,
0) is a fixed point of (3). Following [20], let us first adimension-alize equations (3). Let us observe that we could use different adimension-alization procedures, in particular using the parameter ε SS = E T S T + K M , asproposed in [42]. However, we follow the simpler scheme shown in [20], justin order to test our theoretical results and compare them with the resultsshown in [9]: dudτ = − u + ( u + κ − λ ) v, u (0) = 1 ,ε dvdτ = u − ( u + κ ) v, v (0) = 0 . (38)where τ = k E T t, u = XX T , v = CE T , λ = k k X T , and κ = k + k − k X T = K M X T , ε = E T X T . This is the Heineken-Tsuchiya-Aris system [20]. Carr ([9], p.8, example 3)uses equations dudτ = − u + ( u + c ) v,ε dvdτ = u − ( u + 1) v. (39)To obtain (38) from (39) just impose κ + λ = c and κ = 1. We will startfrom (38) for having a more realistic system. By applying the sQSSA (whichcorresponds to impose ε = 0), we have the reduced system ( outer solution )of (38) dudτ = − λuκ + u ,v = uκ + u . (40)As above remarked, Heineken et al. [20] and Dvoˇr´ak and ˇSiˇska [16] quoteTihonov’s Theorem in order to justify the sQSSA.The aim of this subsection is now to determine the center manifold for(38), using the techniques described in [52, 53] and to show that it is asymp-totically equivalent to the singular points related to Tihonov theory.14o this aim, let us now set τ = εs . Equations (38) can be rewritten inthe equivalent form ( inner solution ) duds = εϕ ( u, v ) ,dvds = u − ( u + κ ) v. (41)where ϕ ( u, v ) = − u + ( u + κ − λ ) v . In order to obtain for (41) a blockform, of the type (13), where the submatrix having eigenvalues with zeroreal parts is separated from the submatrix having eigenvalues with negativereal parts, we operate the substitution w := u − κv , i.e., v = u − wκ . Hence, ϕ ( u, w ) = − u + ( u + κ − λ ) u − wκ and dwds = duds − κ dvds = εϕ ( u, w ) − κu + ( u + κ )( u − w ) = εϕ ( u, w ) + u ( u − w ) − κw. (42)Following [53], the way we will handle parametrized systems consists ofincluding the parameter ε as a new dependent variable as in (43), whichmerely acts to augment the matrix A by adding a new center direction thathas no dynamics. In this way, system (41) becomes duds = εϕ ( u, w ) ,dwds = − κw + u ( u − w ) + εϕ ( u, w ) ,dεds = 0 (43)where the parameter ε is a new variable and the system could have alsoother fixed points.The associated linearized system has a diagonal form, where the eigen-values are given by 0 (multiplicity 2) and − κ .To find a center manifold, all we need to do is to solve (22) for system(43), employing Theorem 2.5, which gives us a method for computing anapproximate solution of (22) to any desired degree of accuracy. Referring to(22) and (13), A = 0, B = − κ , so we search for a function w = h ( u, ε ) suchthat D u h ( u, ε ) (0 + f ( u, h ( u, ε ) , ε )) + κh ( u, ε ) − g ( u, h ( u, ε ) , ε ) = 0 (44)where f ( u, h ( u, ε ) , ε ) = εϕ ( u, h ( u, ε )) ,g ( u, h ( u, ε ) , ε ) = u ( u − h ( u, ε )) + εϕ ( u, h ( u, ε )) . (45)15sing Theorem 2.5 we assume h ( u, ε ) = a u + a uε + a ε + . . . (46)Substituting (46) into (44), one has: ε (2 a u + a ε + . . . ) ϕ ( u, h ( u, ε )) + κ (cid:0) a u + a uε + . . . (cid:1) + − u (cid:0) u − a u − a uε + . . . (cid:1) − εϕ ( u, h ( u, ε )) = 0 (47)where ϕ ( u, h ( u, ε )) = − a u − a uε + · · · − λκ (cid:0) u − a u − a uε + . . . (cid:1) + uκ (cid:0) u − a u − a uε + . . . (cid:1) = − λκ u + (cid:18) − a + 1 κ + λκ a (cid:19) u + (cid:18) − a + λκ a (cid:19) uε + (cid:18) − a + λκ a (cid:19) ε + . . . (48)Accordingly, substituting (48) into (47), one has ε (2 a u + a ε + . . . ) (cid:34) − λκ u + (cid:18) − a + 1 κ + λκ a (cid:19) u + (cid:18) − a + λκ a (cid:19) uε ++ (cid:18) − a + λκ a (cid:19) ε + . . . (cid:35) + κ (cid:0) a u + a uε + a ε + . . . (cid:1) + − u (cid:0) u − a u − a uε − a ε + . . . (cid:1) − ε (cid:34) − λκ u + (cid:18) − a + 1 κ + λκ a (cid:19) u ++ (cid:18) − a + λκ a (cid:19) uε + (cid:18) − a + λκ a (cid:19) ε + . . . (cid:35) = 0 (49)Truncating at second order terms, we obtain:( ka − u + (cid:18) ka + λκ (cid:19) uε + κa ε + · · · = 0Equating terms of the same power to zero gives a = κ , a = − λκ and a = 0. Hence, the center manifold for system (43) is: h ( u, ε ) = 1 κ u − λκ uε + . . . (50)which, for κ = 1, gives the result shown in [9]. Finally, substituting (50) into(43), we obtain the vector field reduced to the center manifold, according16o equation (16) of Theorem 2.3. In fact, if a = κ , a = − λκ and a = 0,formula (48) becomes: ϕ ( u, h ( u, ε )) = − λκ u + λκ u − λκ (cid:18) − λκ (cid:19) uε + . . . Thus: duds = ε (cid:20) − λκ u + λκ u − λκ (cid:18) − λκ (cid:19) uε + . . . (cid:21) ,dεds = 0 (51)or, in terms of the original time scale,˙ u = λκ u (cid:20) − uκ − εκ (cid:18) − λκ (cid:19) + . . . (cid:21) , ˙ ε = 0 (52)Let us now conclude showing that the center manifold obtained followingthis method is asymptotically sufficiently close to (40). We can obtain backthe equation in v . In fact, since v = u − wκ , from (50) and w = h ( u, ε ) = 1 κ u − λκ uε + . . . we have v = uκ (cid:16) − wu (cid:17) = uκ (cid:18) − κ u + λκ ε + . . . (cid:19) (53)Considering ε (cid:28)
1, one has v ∼ uκ (cid:16) − uκ (cid:17) ∼ uκ (cid:18)
11 + uκ (cid:19) = uκ + u , for u → ε (cid:28) ε , respectively. Obviously, the latter gives a better approximation of thenumerical solution of (38), while the former well approximates the sQSSAcurve, which in fact can be considered the zeroth order term of an asymptoticexpansion of the solution in terms of ε .17
10 20 30 40 5000.0050.010.0150.020.0250.030.035 X (Protein) C ( C o m p l e x ) Phase Plane FullsQSSAStandard Center Manifold (Order 0)Standard Center Manifold (Order 1) 0 0.2 0.4 0.6 0.8 100.0050.010.0150.020.0250.030.0350.040.0450.05 X (Protein) C ( C o m p l e x ) Phase Plane FullsQSSAStandard Center Manifold (Order 0)Standard Center Manifold (Order 1)
Figure 2: Comparison in the phase space (
X, C ) of the numerical solu-tion of the system (38) (blue solid line) with its sQSSA (39) (black solidline) and its zeroth order (dashed line) and first order (dashed/dottedline) center manifold (53). The parameter sets are the following. Left: k = 0 . k = 10; k − = 0 . E T = 0 . X T = 50; K M = 100 . K =100; ε HT A = 0 . ε SS = 0 . k = k = 1; k − = 0 . E T = 0 . X T = 1; K M = 1 . K = 1; ε HT A =0 . ε SS = 0 .
05. The set was taken (and modified) from [22]. Since in bothcases the value of ε SS is sufficiently small, the different approximations con-verge to the graph of the numerical solution. In the plot on the right it ispossible to appreciate the different behavior of the zero-th order and thefirst order center manifolds. While the first order manifold approximatesin a better way the numerical solution, the zero-th order converges to thesQSSA, that does not approximate sufficiently well the numerical solution,since it is the zero-th order term of the singular perturbation of the solutionin terms of the parameter ε SS . 18 .2 The total quasi-steady state approximation (tQSSA) Let us now consider system (5). Let us adimensionalize the system, as in[39, 11] dudτ = − v, u (0) = 1 ,ε dvdτ = ησv − ( η + κ m ) v − σuv + u, v (0) = 0 . (55)where ¯ X = uX T , C = (cid:18) E T X T E T + K M + X T (cid:19) v, τ = E T + K M + X T k E T t, and ε = KE T ( E T + K M + X T ) , K = k k with system parameters σ = X T E T + K M + X T , η = E T E T + K M + X T , κ m = k m E T + K M + X T such that σ + η + κ m = 1.By applying the tQSSA (which corresponds to imposing ε = 0) to (55),we have ησv − ( η + κ m ) v − σuv + u = 0and, solving in v : v = η + κ m + σu − (cid:113) ( η + κ m + σu ) − ησu ησ (56)which represents the singular point (or outer solution ) of (55), where η , σ , κ m are viewed as fixed positive constants and ε is the parameter. Its fixedpoint is ( u, v ) = (0 , ε appears already in [32] and was used in [39, 15,11] to determine the asymptotic expansions whose leading term is just thetQSSA. 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(55), using the techniques described in [52, 53] and to show that they areasymptotically equivalent to the singular points related to Tihonov theory.To this aim, let us now set τ = εs , system (55) can be rewritten in theform ( inner solution ) (cid:40) duds = − εv, dvds = u − ( η + κ m ) v + ησv − σuv . (57)19n order to obtain for (57) a block form, of the type (13), we make thesubstitution w := u − ( η + κ m ) v , i.e., v = u − wη + κ m . Hence, dwds = duds − ( η + κ m ) dvds = ε w − uη + κ m + − ( η + κ m ) (cid:32) u − ( η + κ m ) u − wη + κ m + ησ (cid:18) u − wη + κ m (cid:19) − σu u − wη + κ m (cid:33) = − ( η + κ m ) w + σu ( u − w ) + ε w − uη + κ m − ησ ( u − w ) η + κ m Doing so, and introducing the new variable ε , system (57) becomes duds = ε w − uη + κ m dwds = − ( η + κ m ) w + σu ( u − w ) + ε w − uη + κ m − ησ ( u − w ) η + κ m ˙ ε = 0 (58)The associated linearized system has a diagonal form and, in fact, the eigen-values are given by 0 (with multiplicity 2) and − ( η + κ m ).Also in this case the eigenvalue 0 has multiplicity 2.We solve (22) for system (58), employing Theorem 2.5 and determinethe center manifold. Referring to (22) and (13), we have that A = 0, B = − ( η + κ m ). Accordingly, we search a function h ( u, ε ) such that D u h ( u, ε ) (0 + f ( u, h ( u, ε ) , ε )) + ( η + κ m ) h ( u, ε ) − g ( u, h ( u, ε ) , ε ) = 0 f ( u, h ( u, ε ) , ε ) = ε h ( u, ε ) − uη + κ m ,g ( u, h ( u, ε ) , ε ) = σu ( u − h ( u, ε )) + ε h ( u, ε ) − uη + κ m − ησ ( u − h ( u, ε )) η + κ m (59)Using Theorem 2.5 we assume h ( u, ε ) = a u + a uε + a ε + . . . (60)Substituting (60) into (59), one has: ε (2 a u + a ε + . . . ) − u + a u + a uε + a ε + . . .η + κ m ++ ( η + κ m ) (cid:0) a u + a uε + a ε + . . . (cid:1) + σu ( − u + a u + a uε + a ε + . . . )+ − ε − u + a u + a uε + a ε + . . .η + κ m + ησ (cid:0) − u + a u + a uε + a ε + . . . (cid:1) η + κ m = 0(61)20runcating at second order terms, we obtain:( η + κ m ) (cid:0) a u + a uε + a ε (cid:1) − σu + εuη + κ m + ησ u η + κ m = 0so, (cid:20) ( η + κ m ) a − σ + ηση + κ m (cid:21) u ++ (cid:20) ( η + κ m ) a + 1 η + κ m (cid:21) uε + ( η + κ m ) a ε = 0from which: a = ση + κ m − ησ ( η + κ m ) , a = − η + κ m ) , a = 0Hence, the center manifold for system (58) is: h ( u, ε ) = σκ m ( η + κ m ) u − η + κ m ) uε + . . . (62)Finally, substituting (62) into (58) we obtain the vector field reduced to thecenter manifold, according to equation (16) of Theorem 2.3. Then: (cid:40) duds = εη + κ m (cid:104) − u + σκ m ( η + κ m ) u − η + κ m ) uε + . . . (cid:105) , ˙ ε = 0 (63)or, in terms of the original time scale,˙ u = uη + κ m (cid:20) − σκ m ( η + κ m ) u − η + κ m ) ε + . . . (cid:21) , ˙ ε = 0 (64)Let us show that the center manifold obtained in (62) is asymptoticallyclose to the root given by (56), in terms of Tihonov’s Theorem. From (62),and since v = u − wη + κ m , with w = h ( u, ε ) = σκ m ( η + κ m ) u − η + κ m ) uε + . . . we have v = uη + κ m (cid:16) − wu (cid:17) = uη + κ m (cid:18) − σκ m ( η + κ m ) u + 1( η + κ m ) ε + . . . (cid:19) (65)21ince equation (56) is obtained putting ε (cid:28)
1, one has v ∼ uη + κ m (cid:18) − σκ m ( η + κ m ) u (cid:19) ∼ uη + κ m (cid:32)
11 + σκ m ( η + κ m ) u (cid:33) = uη + κ m + σκ m η + κ m u , for u → v = η + κ m + σu ησ (cid:32) − (cid:115) − ησu ( η + κ m + σu ) (cid:33) = η + κ m + σu ησ (cid:2) − √ − ε (cid:3) ∼ η + κ m + σuησ ε = uη + κ m + σu , for u → uη + κ m when u →
0. This means that both the expressions can be intended as two differentapproximations of the center manifold.In Figure (3) we compare the tQSSA of system (55), obtained from (56),with the center manifold (65), at the zeroth order and at the first order in ε , respectively. Obviously, the latter gives a better approximation of thenumerical solution of (55), while the former well approximates the tQSSAcurve, which in fact can be considered the zeroth order term of an asymptoticexpansion of the solution in terms of ε . Let us consider now a more general system of the following form ( outersolution ) (cid:40) dudτ = ϕ ( u, v ) ,ε dvdτ = au + bv + ψ ( u, v ) , a, b ∈ R , b < inner solution (cid:40) duds = εϕ ( u, v ) , dvds = au + bv + ψ ( u, v ) , a, b ∈ R , b < τ = εs ) where ϕ (0 ,
0) = ψ (0 ,
0) = 0 , and ψ u (0 ,
0) = ψ v (0 ,
0) = 0 . (70)The origin is a fixed point for (69). Heineken-Tsuchiya-Aris system (41)and the system obtained by the tQSSA approximation (57), are particularcases of the system (69)-(70). 22igure 3: Comparison in the phase space ( ¯ X, C ) of the numerical solution ofthe system (55) (blue solid line) with its tQSSA (56) (black solid line) andits zeroth order (dashed line) and first order (dashed/dotted line) centermanifold (65). The parameter sets are the following. Left: k = k =1; k − = 3; E T = 1; X T = 1; K M = 4; K = 1; ε HT A = 1; ε SS = 0 . ε = 0 . k = 0 . k = 10; k − = 0 . E T =400; X T = 100; K M = 100 . K = 100; ε HT A = 4; ε SS = 2; ε = 0 .
11. The setwas taken from [14]. In the plot on the left, since the value of ε is sufficientlysmall, the different approximations converge to the graph of the numericalsolution. In the plot on the right it is possible to appreciate the differentbehavior of the zero-th order and the first order center manifolds. While thefirst order manifold approximates in a better way the numerical solution, thezero-th order converges to the tQSSA, that does not approximate sufficientlywell the numerical solution, since it is the zero-th order term of the singularperturbation of the solution in terms of the parameter ε .23e are able to state a more general theorem concerning the center man-ifold, which is the main result of our paper.Let w := au + bv ; hence: dwds = a duds + b dvds = a [ εϕ ( u, v )] + b [ au + bv + ψ ( u, v )]= bw + aεϕ (cid:18) u, w − aub (cid:19) + bψ (cid:18) u, w − aub (cid:19) Doing so, system (69) becomes, for a, b ∈ R and b < duds = εϕ (cid:18) u, w − aub (cid:19) ,dwds = bw + aεϕ (cid:18) u, w − aub (cid:19) + bψ (cid:18) u, w − aub (cid:19) dεds = 0 (71)The associated linearized system has a block form of type (13) and, in fact,the eigenvalues are given by 0 (with multiplicity 2) and b <
0. Thus in everysystem of the form (71) we are in presence of a center manifold.We write equation (22) for system (71), employing Theorem 2.5. Refer-ring to (22) and (13), we have that A = 0, B = b . Accordingly, we searchfor a function w = h ( u, ε ) such that D u h ( u, ε ) (cid:18) εϕ (cid:18) u, w − aub (cid:19)(cid:19) + − bh ( u, ε ) − aεϕ (cid:18) u, w − aub (cid:19) − bψ (cid:18) u, w − aub (cid:19) = 0from which, since D u h ( u, ε ) εϕ (cid:0) u, w − aub (cid:1) is a function at least of third orderin ε and u , while we are interested in a second order expression of function h ( u, ε ), we can neglect this term and focus on bh ( u, ε ) + aεϕ (cid:18) u, w − aub (cid:19) + bψ (cid:18) u, w − aub (cid:19) = 0 (72) 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 h ( u, ε ) = λ u + λ uε + λ ε + . . . (73)24nd it is trivial to prove that h ( u, ε ) satisfies (72) for λ = 0. Moreover,from (70), ψ (cid:18) u, w − aub (cid:19) = 12 (cid:2) Θ( u , uw, w ) (cid:3) + . . . where Θ( u , uw, w ) contains the quadratic terms in u and w .Since the terms in uw and w , with w = h ( u, ε ) = λ u + λ uε + . . . , areat least of third order in ε and u , we consider only term in u . Therefore, ψ (cid:18) u, w − aub (cid:19) = 12 (cid:20) ψ uu (0 , − ab ψ u,v (0 ,
0) + (cid:16) ab (cid:17) ψ v,v (0 , (cid:21) u + . . . while for ϕ it is sufficient to consider the first order expansion in u because,otherwise, in (72) we would have third order terms for εϕ (cid:0) u, w − aub (cid:1) in ε and u . Thus, ϕ (cid:18) u, w − aub (cid:19) = (cid:104) ϕ u (0 , − ab ϕ v (0 , (cid:105) u + . . . where we recall that v = w − aub . Accordingly, equation (72) becomes: b (cid:0) λ u + λ uε + . . . (cid:1) + aεu (cid:104) ϕ u (0 , − ab ϕ v (0 , (cid:105) + b (cid:20) ψ uu (0 , − ab ψ u,v (0 ,
0) + (cid:16) ab (cid:17) ψ v,v (0 , (cid:21) u + · · · = 0 (74)Equating to zero terms of the same power gives λ = − (cid:20) ψ uu (0 , − ab ψ u,v (0 ,
0) + (cid:16) ab (cid:17) ψ v,v (0 , (cid:21) λ = − ab (cid:104) ϕ u (0 , − ab ϕ v (0 , (cid:105) λ = 0 (75)Hence, the center manifold for system (69) is w = h ( u, ε ) = − (cid:20) ψ uu (0 , − ab ψ u,v (0 ,
0) + (cid:16) ab (cid:17) ψ v,v (0 , (cid:21) u − ab (cid:104) ϕ u (0 , − ab ϕ v (0 , (cid:105) uε + . . . (76)Setting in the RHS ε = 0, we obtain the center manifold w = h ( u,
0) of(69).
Step 2. Singular Point Technique [51]
On the other hand, v = w − aub = λ u + λ uε + λ ε + · · · − aub Since, setting ε = 0, we have that v = w − aub = λ u − aub , equation (72)becomes 25 u + ψ ( u, v ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) v = λ u − aub = 0 (77)which gives an identity up to O ( u ), if we substitute λ as above and if weoperate a Taylor expansion around ( u, v ) = (0 , duds = εϕ (cid:18) u, h ( u, ε ) − aub (cid:19) ,dεds = 0 (78)or, in terms of the original time scale,˙ u = ϕ (cid:18) u, h ( u, ε ) − aub (cid:19) ˙ ε = 0 (79)Moreover: ϕ ( u, v ) = ∂ϕ∂u ∂u∂u (cid:124)(cid:123)(cid:122)(cid:125) =1 + ∂ϕ∂v ∂v∂u (cid:124)(cid:123)(cid:122)(cid:125) = − a/b u + ∂ϕ∂u ∂u∂w (cid:124)(cid:123)(cid:122)(cid:125) =0 + ∂ϕ∂v ∂v∂w (cid:124)(cid:123)(cid:122)(cid:125) =1 /b w + . . . where v = w − aub and all the derivatives are calculated in (0 , ϕ ( u, v ) = (cid:104) ϕ u (0 , − ab ϕ v (0 , (cid:105) u + w ϕ v (0 , b + . . . and, since λ = − ab (cid:2) ϕ u (0 , − ab ϕ v (0 , (cid:3) , we have: ϕ ( u, v ) = − ba λ u + w ϕ v (0 , b + . . . for v = w − aub . From (73), the vector field reduced to the center manifold, interms of the original time scale, near the origin, becomes:˙ u = − ba λ u + (cid:0) λ u + λ uε (cid:1) ϕ v (0 , b + o (cid:0) ε + u (cid:1) ˙ ε = 0 (80)for values of λ and λ as above.Summarizing, we have obtained two relations:26 ) From the Center Manifold Theory : considering (72), for w = h ( u, ε ), and setting ε = 0, we have: bw + bψ (cid:18) u, w − aub (cid:19) = 0 (81) b) From Singular Perturbation Techniques : by assumption B ofsection 2.2, and since g ( u, w ) = bw + aεϕ (cid:0) u, w − aub (cid:1) + bψ (cid:0) u, w − aub (cid:1) in (71),we have, putting ε = 0: g ( u, φ ( u )) = 0 ⇒ bφ ( u ) + bψ (cid:18) u, φ ( u ) − aub (cid:19) = 0 (82)Comparing (81) and (82), we observe a relation between h ( u,
0) and φ ( u )but we cannot infer that h ( u,
0) = φ ( u ), due to the non-uniqueness of centermanifold. However, in the above steps we have proven that h ( u, ∼ φ ( u ) , for u → C = 0, as a manifold which is asymptoti-cally equivalent to the center manifold, as confirmed by equation (54) forHeineken-Tsuchiya-Aris system, and by equations (66)-(67) for the tQSSA.This explains why, in order to achieve the center manifold of (38) and(55), it is sufficient to consider - for u → ε = 0). We recallthat in many papers (see, for example, [20, 16, 21], who refer to Tihonov’sTheorem, and [22], according to Singular Perturbation Theory), the centermanifold is obtained equating to zero the right hand side of the equation ofthe form: ε dydt = g ( x, y ) . The quasi-steady state approximation has been a challenge for applied math-ematicians, who had to explain the feasibility of an approximation whichimposes to the complex C both to be constant and to depend on X . SomeBiochemistry texts (see, for example, [24, 54, 37, 19]) mislead the reader,interpreting the QSSA as a true equality, which brings to assert that the ra-tio E ( t ) S ( t ) /C ( t ) is constant during all the quasi-steady state phase. Thisis obviously not true. In [5] the authors solve the apparent incongruence,determining the asymptotic value of E ( t ) S ( t ) /C ( t ), showing that, for everychoice of the kinetic parameters and of the initial conditions,27 as X as C as ( t ) → (cid:18) k − αα (cid:19) E T =: K W ; ( K D < K W < K M ) (84)(where α = k ( K M + E T ) (cid:104) − (cid:113) − k E T k ( K M + E T ) (cid:105) ), differently from whatis wrongly stated.Heineken et al. [20] and successively other authors [42, 39, 11] interpretedsQSSA and tQSSA as leading order expansions of the solutions in terms ofa suitable parameter, which has to be considered small.This interpretation allows us to embed the QSSA theory in a frameworkwhich is related to Tihonov’s Theorem [46, 47, 48, 20, 42, 51, 16, 21], wherethe parameter multiplies the derivative of C and the QSSA can be obtainedas the singular point of the original system, setting ε = 0.In this paper we have shown that, at least in the classical simple scheme(3), the approximation obtained applying Tihonov’s Theorem is asymptot-ically equivalent to the center manifold of the system, which means thatreduced system and center manifold are two sides of the same coin.Once again, the total QSSA has shown to be much more efficient andnatural than the standard one, mainly thanks to the fact that the parameterused for the expansions in the total framework is always less than 14 .In our actual researches we are applying the techniques shown in thispaper to more complex enzyme reactions, as the fully competitive inhibition[4], the phosphorylation-dephosphorylation cycle (or Goldbeter-Koshlandswitch [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 previouspapers [33, 3, 34, 35, 36, 12, 13, 6].The techniques here shown will allow to read the tQSSA as the lead-ing term of an asymptotic expansion in terms of a suitable perturbationparameter, in these more complex cases, too. Acknowledgements
The authors are deeply grateful to Prof. Enzo Orsingher, from SapienzaUniversity (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.