2006 International Conference in Honor of Jacqueline Fleckinger. Electronic Journal of Differential Equations, Conference 16, 2007, pp. 155–184. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu ftp ejde.math.txstate.edu (login: ftp)

REDUCTION FOR MICHAELIS-MENTEN-HENRI KINETICS IN THE PRESENCE OF DIFFUSION

LEONID V. KALACHEV, HANS G. KAPER, TASSO J. KAPER, NIKOLA POPOVIĆ, ANTONIOS ZAGARIS

Dedicated to Jacqueline Fleckinger on her 65-th birthday

Abstract. The Michaelis-Menten-Henri (MMH) mechanism is one of the paradigm reaction mechanisms in biology and chemistry. In its simplest form, it involves a substrate that reacts (reversibly) with an enzyme, forming a complex which is transformed (irreversibly) into a product and the enzyme. Given these basic kinetics, a dimension reduction has traditionally been achieved in two steps, by using conservation relations to reduce the number of species and by exploiting the inherent fast–slow structure of the resulting equations. In the present article, we investigate how the dynamics change if the species are additionally allowed to diffuse. We study the two extreme regimes of large diffusivities and of small diffusivities, as well as an intermediate regime in which the time scale of diffusion is comparable to that of the fast reaction kinetics. We show that reduction is possible in each of these regimes, with the nature of the reduction being regime dependent. Our analysis relies on the classical method of matched asymptotic expansions to derive approximations for the solutions that are uniformly valid in space and time.

1. Introduction

One of the paradigm reaction mechanisms in biology and chemistry—often referred to as the Michaelis-Menten-Henri (MMH) mechanism—involves a substrate (S) that reacts (reversibly) with an enzyme (E) to form a complex (C) which, in turn, is transformed (irreversibly) into a product (P) and the enzyme [10, 16]. The reaction mechanism is represented symbolically by

S + E ⇌_{k_{−1}}^{k_{1}} C →^{k_{2}} P + E,  (1.1)

where k_{1}, k_{−1}, and k_{2} are rate constants.

The MMH mechanism models the kinetics of many fundamental reactions. Examples from biochemistry include those discussed in [4, 6, 7, 8, 17, 21, 22, 23, 24, 26,

2000 Mathematics Subject Classification. 35K57, 35B40, 92C45, 41A60. Key words and phrases. Michaelis-Menten-Henri mechanism; diffusion; dimension reduction; matched asymptotics. \begin{matrix}c \\ ◯ \end{matrix}2007 Texas State University - San Marcos. Published May 15, 2007.

156 L. KALACHEV, H. KAPER, T. KAPER, N. POPOVIĆ, A. ZAGARIS EJDE/CONF/16

30], and examples involving nutrient uptake in cells and heterogeneous catalytic reactions are analyzed in [5, Chapter 7.1] and [2], respectively. The MMH mechanism is also presented as a prototypical mechanism exhibiting fast and slow dynamics— and, hence, the potential for dimension reduction—in numerous textbooks, see, e.g., [13, 15, 19, 20].

1.1. MMH Kinetics. Equations governing the kinetics of (1.1) may be derived from the law of mass action,

\frac{dS}{dt̃} = −k_{1}SE + k_{−1}C,  (1.2a)
\frac{dE}{dt̃} = −k_{1}SE + (k_{−1} + k_{2})C,  (1.2b)
\frac{dC}{dt̃} = k_{1}SE − (k_{−1} + k_{2})C,  (1.2c)
\frac{dP}{dt̃} = k_{2}C,  (1.2d)

where S, E, C, and P denote the concentrations of substrate, enzyme, complex, and product, respectively. The initial conditions are S(0) = S̄, E(0) = Ē, C(0) = 0, and P (0) = 0.

Dimension reduction in (1.2) is traditionally achieved in two steps. The first step uses the pair of conservation relations that exist for the mechanism (1.1). In particular, the total concentration of enzyme (free and bound in complex) is constant; that is, E(t) + C(t) = Ē for all t > 0. In addition, S(t) + C(t) + P (t) = S̄ for all t > 0. Therefore, there is a decrease from four variables to two, and the governing equations (1.2) reduce to

\frac{dS}{dt̃} = −k_{1}ĒS + (k_{1}S + k_{−1})C,  (1.3a)
\frac{dC}{dt̃} = k_{1}ĒS − (k_{1}S + k_{−1} + k_{2})C.  (1.3b)

The second reduction step exploits the separation of time scales. In particular, Ē ≪ S̄. Hence, there is a small, positive, dimensionless parameter,

ε = \frac{Ē}{S̄} ,  (1.4)

and the nondimensionalized equations are naturally formulated as a fast–slow system,

ṡ = −s + (s + κ − λ)c,  (1.5a)
εċ = s − (s + κ)c,  (1.5b)

where

t = k_{1}Ēt̃, s = \frac{S}{S̄}, e = \frac{E}{Ē} , c = \frac{C}{Ē} .  (1.6)

The dimensionless parameters are

κ = \frac{k_{−1} + k_{2}}{k_{1}S̄} , λ = \frac{k_{2}}{k_{1}S̄}.  (1.7)

During a short O(ε) initial transient period, the variable c is fast and rises rapidly to its maximum value, while the variable s is slow and remains essentially constant. Subsequently, the evolution of c is slaved to that of s, and both c and s evolve slowly toward their equilibrium value, zero. This slaving is often referred to as reduced kinetics. From the point of view of dynamical systems, the system (1.5) has an asymptotically stable, invariant, slow manifold. During the transient, the concentrations relax to the slow manifold, decaying exponentially toward it. Subsequently, on the O(1) time scale, the reaction kinetics play out near this manifold. Since the slow manifold is only one-dimensional, a reduction is achieved.

To leading order, the slow manifold is given by

c = \frac{s}{s + κ},  (1.8)

and the reduced system dynamics on it are governed by a single equation,

\frac{ds}{dt} = −\frac{λs}{s + κ}.

This leading order slow manifold is obtained directly from (1.5b) with ε = 0 and is referred to as the quasi-steady state approximation [1, 13, 15, 19, 22, 25] in chemistry and as the critical manifold in mathematics. Higher-order corrections to the critical manifold, which is sufficiently accurate for many applications, may be calculated using geometric singular perturbation theory, see for example [11, 12]. We emphasize that the critical manifold is only approximately invariant under the dynamics of (1.5); the exact slow manifold is invariant. Additional studies of the accuracy of the quasi-steady state approximation are given in [1, 9, 15, 20, 21, 25, 27].

Remark 1.1. For all sufficiently small, positive ε, there is a family of slow manifolds, all of which are exponentially (O(e^{−k̃/ε}) for some k̃ > 0) close to each other, i.e., the asymptotic expansions of these slow manifolds are the same to all powers of ε, see [11, 12], for example. For convenience, we will sometimes refer to ‘the’ (rather than to ‘a’) slow manifold.

1.2. MMH Kinetics with Diffusion. Given the effectiveness of the two reduction steps in the kinetics problem (1.2), one is naturally led to ask what happens when the species are simultaneously permitted to diffuse, and whether any similar reduction can be achieved. The conservation relations used in the first reduction step of the kinetics analysis do not generalize. However, there is still a separation of time scales in the reaction kinetics, and the process of diffusion introduces one (or more) additional time scale(s). Therefore, one expects that dimension reduction may still be achieved by exploiting the separation of time scales, and the purpose of this article is to investigate this possibility.

The problem with diffusion is governed by the evolution of the concentrations of substrate, enzyme, and complex in time and space. The concentration of product can be found by quadrature as a function of these other concentrations, since the second reaction in (1.1) is irreversible. The governing equations in one space dimension are

\frac{∂S}{∂t̃} = −k_{1}SE + k_{−1}C + D_{S} \frac{∂^{2}S}{∂x̃^{2}} ,  (1.9a)
\frac{∂E}{∂t̃} = −k_{1}SE + (k_{−1} + k_{2})C + D_{E} \frac{∂^{2}E}{∂x̃^{2}} ,  (1.9b)
\frac{∂C}{∂t̃} = k_{1}SE − (k_{−1} + k_{2})C + D_{C} \frac{∂^{2}C}{∂x̃^{2}} ,  (1.9c)

with x̃ ∈ [0, ℓ], subject to no-flux (Neumann) boundary conditions

\frac{∂S}{∂x̃} |_{x̃=0,ℓ} = \frac{∂E}{∂x̃} |_{x̃=0,ℓ} = \frac{∂C}{∂x̃} |_{x̃=0,ℓ} = 0  (1.10)

and initial conditions

S(0, x̃) = S_{i}(x̃), E(0, x̃) = E_{i}(x̃), C(0, x̃) = 0.  (1.11)

Here, ℓ > 0 is the O(1) size (length) of the reactor; D_{S}, D_{E}, and D_{C} denote the diffusivities of S, E, and C, respectively; and S_{i} and E_{i} are given, smooth functions describing the initial spatial profiles of substrate and enzyme, respectively.

We nondimensionalize (1.9)–(1.11) as follows. The nondimensional spatial variable is x = x̃/ℓ. Time and the species’ concentrations are nondimensionalized as in (1.6), but now S̄ and Ē denote the spatial averages,

S̄ = \frac{1}{ℓ} ∫_{0}^{ℓ} S_{i}(x̃) dx̃,  Ē = \frac{1}{ℓ} ∫_{0}^{ℓ} E_{i}(x̃) dx̃.

The nondimensional parameters are again given by (1.7), and the diffusivities are scaled as

δ = \frac{D_{S}}{k_{1}ℓ^{2}Ē} , a = \frac{D_{E}}{D_{S}} , b = \frac{D_{C}}{D_{S}} .  (1.12)

Thus, we obtain the equations

\frac{∂s}{∂t} = −se + (κ − λ)c + δ \frac{∂^{2}s}{∂x^{2}} ,  (1.13a)
\frac{∂e}{∂t} = −\frac{1}{ε}(se − κc) + aδ \frac{∂^{2}e}{∂x^{2}} ,  (1.13b)
\frac{∂c}{∂t} = \frac{1}{ε}(se − κc) + bδ \frac{∂^{2}c}{∂x^{2}} ,  (1.13c)

on the unit interval, subject to the Neumann boundary conditions

\frac{∂s}{∂x}|_{x=0,1} = \frac{∂e}{∂x}|_{x=0,1} = \frac{∂c}{∂x}|_{x=0,1} = 0,  (1.14)

and the initial conditions

s(0, x) = s_{i}(x), e(0, x) = e_{i}(x), c(0, x) = 0.  (1.15)

Here, s_{i} = S_{i}/S̄ and e_{i} = E_{i}/Ē. We assume 0 < ε ≪ 1 and that a and b are O(1). In vector notation, equations (1.13) may be written as

\frac{∂u}{∂t} = \frac{1}{ε}F_{ε}(u) + δD \frac{∂^{2}u}{∂x^{2}},  (1.16)

where

F_{ε}(u) = (\begin{matrix}−εse + ε(κ − λ)c \\ −(se − κc) \\ se − κc \end{matrix}) , u = (\begin{matrix}s \\ e \\ c \end{matrix}) ,  (1.17)

and D = diag(1, a, b). Moreover, we use [F_{ε}(u)]_{k} to denote the kth order terms in the Taylor expansion of F_{ε} with respect to ε. Given a formal asymptotic expansion u(·, x, ε) = ∑_{k=0}^{∞} u_{k}(·, x)ε^{k} of the solution of (1.16), [F_{ε}(u)]_{k} will generically be a function of u_{0}, …, u_{k}.

1.3. Summary of Main Results. The impact of diffusion depends on the time scales associated with the species’ diffusivities relative to those of the reaction kinetics. We examine a spectrum of species’ diffusivities here:

(i) Large diffusivities, δ = O(1/ε^{2}): the diffusive time scale is shorter than both the fast and the slow kinetic time scales;

(ii) Moderately large diffusivities, δ = O(1/ε): the diffusive time scale is com-

parable to the time scale of the fast kinetics; and (iii) Small diffusivities, δ = O(ε): the diffusive time scale is longer than both

kinetic time scales. Our principal findings are that reduction is possible in all regimes under consideration.

In regime (i), diffusion effectively homogenizes the concentrations of all three species on the super-fast (τ = t/ε^{2}) time scale. Then, the dynamics on the fast (η = t/ε) and slow (t) time scales are given by the classical MMH kinetics mechanism, with the fast reactions occurring on the fast scale and the reduced kinetics taking place on the slow time scale. We treat this regime primarily to introduce the method we use throughout.

In regime (ii), the species undergo both diffusion and the fast reaction on the fast (η = t/ε) time scale. In particular, the substrate concentration satisfies the homogeneous heat equation to leading order; hence, it homogenizes exponentially in time. The enzyme and complex concentrations satisfy nonautonomous, linear reaction–diffusion equations to leading order, and they also homogenize exponentially in time, approaching points on the classical critical manifold. Then, on the slow (t) time scale, the solution is essentially spatially homogeneous. The concentration of substrate evolves according to the classical reduced equation, while the enzyme and complex concentrations are constrained to lie on the critical manifold, to leading order. Most significantly, these leading-order results are independent of the diffusivities of the enzyme and complex, even when these diffusivities are unequal.

In regime (iii), the MMH reaction kinetics take place at every point in the domain effectively decoupled from the kinetics at any other point. On the fast (η = t/ε) time scale, enzyme and substrate bind to form complex with the amount of complex at each point x depending on the initial enzyme concentration e_{i}(x), while on the slow (t) time scale the substrate and complex concentrations slowly approach equilibrium in an x-dependent manner. We label these dynamics as pointwise fast kinetics and pointwise slow, reduced kinetics, respectively. Also, on the slow time scale, the enzyme concentration returns essentially to the initial enzyme profile. Then, on asymptotically large or super-slow (ζ = εt) time scales, the enzyme profile homogenizes.

The observed dynamics and the time scales in these regimes are summarized in Table 1.

We use matched asymptotic expansions in time in a straightforward manner in each of the regimes identified above. (Equivalent results could, for example, be obtained via the so-called boundary function method [29].) Moreover, we present numerical simulations in every regime to further illustrate our analysis.

Remark 1.2. The regime of moderately small diffusivities, δ = O(1), will be analyzed in a separate article. In this regime, the diffusive time scale is comparable to that of the slow kinetics. Preliminary results suggest that the fast dynamics are

Regime δ  Dynamics  Time scale
(i)  1/ε^{2} homogenization of s, e, c  super-fast τ = t/ε^{2}
fast kinetics  fast η = t/ε
slow reduced kinetics  slow t
(ii)  1/ε homogenization & fast kinetics fast η = t/ε
slow reduced kinetics  slow t
(iii) ε pointwise fast kinetics fast η = t/ε
pointwise slow reduced kinetics slow t
homogenization of e  super-slow ζ = εt

Table 1. Summary of the observed dynamics and the time scales of (1.13) in regimes (i)–(iii)

similar to those in regime (iii), but without the concentrations becoming homogenized, while the slow dynamics are governed by a reaction-diffusion equation for the concentration of substrate, with the concentrations of enzyme and complex slaved to it.

Remark 1.3. The MMH mechanism in the presence of diffusion is analyzed here as a prototype problem. The method we employ here may be used for other mechanisms with one or more kinetics time scales.

Remark 1.4. The analysis may also be extended to problems in which the domain length ℓ is not O(1). For example, the analysis of regime (i) also applies to problems in which ℓ is small and the diffusivities are not large. In that case, it is natural to scale the spatial variable as x = εx̃. With this scaling, the diffusion terms in (1.13) are of the form

δ \frac{∂^{2}s}{∂x^{2}} = \frac{δ}{ε^{2}} \frac{∂^{2}s}{∂x̃^{2}} .

Hence, diffusion dominates again, even if the actual diffusion coefficients are δ = O(1).

Remark 1.5. The influence of diffusion in the MMH mechanism has also been studied in [31]. Specifically, the reduced kinetics model (1.5) is considered and a diffusion term is added for the substrate only, with O(1) diffusion coefficient. Via an inertial manifold approach, it is shown that this system may be reduced to a single reaction-diffusion equation for s, in which the diffusivity has a concentrationdependent correction at O(ε). The fast transients for this model are also calculated, and extensions are given for general systems with fast–slow kinetics in which the slow species also diffuse.

Remark 1.6. It has been shown in [27] that the effective small parameter in the MMH mechanism is ε̃ = Ē/(S̄+K_{M}), where K_{M} = (k_{−1} +k_{2})/k_{1} is the Michaelis– Menten constant. Hence, there is a wider range of physical parameters for which one has a separation of time scales. Our method can be applied to the equations with this small parameter as well; however, here we use the traditional scaling, since it is still the one that is most commonly used.

This article is organized as follows. The regimes (i)–(iii) are analyzed in Sections 2–4, respectively. In Section 5, the results of the preceding sections are discussed, and the theoretical results are further illustrated using numerical simulations. In Appendix A, it is shown via a Turing analysis that the homogeneous attractor of (1.13) is linearly stable, irrespective of the magnitudes of the diffusion coefficients. Appendix B contains a technical result relating to Section 3.

2. Large Diffusivities

In this section, we consider the regime in which the diffusivities of all species are large, δ = O(1/ε^{2}); for convenience, we choose δ = 1/ε^{2} in (1.13). Here, the time scale on which diffusion acts is much shorter than that of the fast kinetics. There is a very short transient period, O(ε^{2}) in duration, in which the initially heterogeneous species’ concentrations, given by (1.15), homogenize and during which essentially no reaction takes place. After this short transient, the problem reduces to the wellunderstood, classical problem of pure kinetics for the homogeneous solution, see, e.g., [15]. We treat this regime in some detail to introduce the method employed in this article in an elementary context.

After introduction of δ = 1/ε^{2}, equations (1.13) become

ε^{2} \frac{∂s}{∂t} = −ε^{2}se + ε^{2}(κ − λ)c + \frac{∂^{2}s}{∂x^{2}} ,  (2.1a)
ε^{2} \frac{∂e}{∂t} = −ε(se − κc) + a\frac{∂^{2}e}{∂x^{2}} ,  (2.1b)
ε^{2} \frac{∂c}{∂t} = ε(se − κc) + b\frac{∂^{2}c}{∂x^{2}} .  (2.1c)

Equivalently, in vector form,

ε^{2} \frac{∂u}{∂t} = εF_{ε}(u) + D \frac{∂^{2}u}{∂x^{2}},  (2.2)

where F is defined in (1.17).

2.1. Homogenization: The Super-Fast (Inner) Time Scale. To study the initial transient period, we let τ = t/ε^{2} denote the super-fast (inner) time and let û(τ, x, ε) = u(t, x, ε). The governing equations become

\frac{∂û}{∂τ} = εF_{ε}(û) + D \frac{∂^{2}û}{∂x^{2}},  (2.3)

with initial and boundary conditions

û(0, x, ε) = u_{i}(x) and \frac{∂û}{∂x}|_{x=0,1} = 0.  (2.4)

We consider (2.3) and (2.4) over an O(1)–interval of τ time, starting at τ = 0. Asymptotically, as ε → 0^{+}, the solution can be expressed using the Ansatz

û(τ, x, ε) = û_{0}(τ, x) + εû_{1}(τ, x) + O(ε^{2}).

Hence, we expand both sides of (2.3) in powers of ε to obtain a recursive sequence of differential equations for û_{k}, k = 1, 2, . . .. At O(1), û_{0} satisfies the homogeneous heat equation,

O(1) : L_{τ} û_{0} = 0,  (2.5)

where L_{τ} = ∂/∂τ − D∂^{2}/∂x^{2}, subject to Neumann boundary conditions. The solution is

û_{0}(τ, x) = ∑_{k=0}^{∞} û_{0k}e^{−D(kπ)^{2}τ} cos (kπx).  (2.6)

Here, the coefficients û_{0k} are the Fourier coefficients of the initial distribution u_{i}(x) with respect to {cos (kπx)}_{k≥0}, and they are constant during the fast transients on the τ time scale.

Asymptotically, as τ → ∞, we find

û_{0}(τ, x) → û_{00} = ∫_{0}^{1} u_{i}(x) dx = (\begin{matrix}1 \\ 1 \\ 0 \end{matrix}) ,  (2.7)

where we used (1.15). Therefore, asymptotically, the effect of diffusion in (2.1) is to smear out the initial distributions of the reactants until they are effectively uniformly distributed over the entire spatial domain.

Remark 2.1. It will suffice to consider the leading-order fast solution (2.6) to accomplish matching to lowest order in the next subsection. To that end, it is useful to write (2.6) as

û_{0}(τ, x) = (1, 1, 0)^{T} + O(e^{−dπ^{2}τ}),  (2.8)

where d = min{1, a, b} and where we used (2.7).

2.2. Fast Kinetics: The Fast Time Scale. The fast kinetics take place on the fast η = t/ε time scale, during which the system dynamics are given by

ε\frac{∂ũ}{∂η} = εF_{ε}(ũ) + D \frac{∂^{2}ũ}{∂x^{2}},  (2.9)

subject to Neumann boundary conditions. Here, ũ = ũ(η, x, ε), and we assume ũ(η, x, ε) = ũ_{0}(η, x) + εũ_{1}(η, x) + O(ε^{2}). Then, expanding (2.9) in powers of ε and rearranging the resulting equations, we find

O(1) : − D \frac{∂^{2}ũ_{0}}{∂x^{2}} = 0,  (2.10a)
O(ε) : − D \frac{∂^{2}ũ_{1}}{∂x^{2}} = −\frac{∂ũ_{0}}{∂η} + [F_{ε}(ũ)]_{0},  (2.10b)

subject to Neumann boundary conditions. It will suffice to consider the dynamics to this order to obtain a uniform leading-order approximation to the solution of the original system (2.1).

Integrating (2.10a) and taking into account the boundary conditions, we conclude that

ũ_{0}(η, x) = ũ_{0}(η),  (2.11)

i.e., that ũ_{0} is independent of x. Similarly, it follows from (2.10b) that

−D ∫_{0}^{1} \frac{∂^{2}ũ_{1}}{∂x^{2}} dx = −D \frac{∂ũ_{1}}{∂x} |_{x=0}^{1} = 0 = ∫_{0}^{1} ( − \frac{dũ_{0}}{dη} + [F_{ε}(ũ)]_{0}) dx.  (2.12)

Since the integrand in the right member of (2.12) is independent of x, we see that the dynamics of ũ_{0} are governed by the ordinary differential equation

\frac{dũ_{0}}{dη} = [F_{ε}(ũ)]_{0},

which we write out componentwise,

\frac{ds̃_{0}}{dη} = 0,  (2.13a)
\frac{dẽ_{0}}{dη} = −(s̃_{0}ẽ_{0} − κc̃_{0}),  (2.13b)
\frac{dc̃_{0}}{dη} = s̃_{0}ẽ_{0} − κc̃_{0}.  (2.13c)

These equations are the same as one finds to leading order in the classical MMH kinetics problem; hence, they can be solved explicitly: s̃_{0}(η) ≡ s̃_{0}(0) and

ẽ_{0}(η) + c̃_{0}(η) = ẽ_{0}(0) + c̃_{0}(0) for all η ≥ 0.

The initial conditions for (2.13) are determined by matching with the leading-order equations on the super-fast (τ) scale; notably, we require that lim_{τ→∞} û_{0}(τ, x) = lim_{η→0^{+}} ũ_{0}(η, x). Now, by (2.8),

lim_{τ→∞} û_{0}(τ, x) = (1, 1, 0)^{T} .  (2.14)

Hence, we have

s̃_{0}(η) ≡ 1 and ẽ_{0}(η) = 1 − c̃_{0}(η).  (2.15)

In turn, it follows that

\frac{dc̃_{0}}{dη} = 1 − (1 + κ)c̃_{0}, with c̃_{0}(0) = 0,

and, hence,

ũ_{0}(η, x) = (1, \frac{1}{1 + κ}(κ + e^{−(1+κ)η}), \frac{1}{1 + κ}(1 − e^{−(1+κ)η}))^{T} .  (2.16)

Therefore, on the fast (η) scale, the species’ concentrations are essentially homogeneous, and the fast chemistry occurs, with the binding of enzyme and substrate to form complex.

2.3. Slow Reduced Kinetics: The Slow (Outer) Time Scale. The slow, reduced kinetics take place on the slow (t) time scale. The dynamics are governed by the original system, (2.2), subject to Neumann boundary conditions. We expand u(t, x, ε) = u_{0}(t, x) + εu_{1}(t, x) + ε^{2}u_{2}(t, x) + O(ε^{3}) and equate coefficients of equal powers of ε to obtain

O(1) : − D \frac{∂^{2}u_{0}}{∂x^{2}} = 0,  (2.17a)
O(ε) : − D \frac{∂^{2}u_{1}}{∂x^{2}} = [F_{ε}(u)]_{0},  (2.17b)
O(ε^{2}) : − D \frac{∂^{2}u_{2}}{∂x^{2}} = −\frac{∂u_{0}}{∂t} + [F_{ε}(u)]_{1}.  (2.17c)

Applying the same type of solvability argument used in Section 2.2, we conclude from (2.17a) that u_{0}(t, x) = u_{0}(t). Similarly, the solvability for (2.17b) implies u_{1}(t, x) = u_{1}(t), since the right member of (2.17b) is independent of x and, hence, [F_{ε}(u)]_{0} = 0. In turn, this yields

s_{0}(t)e_{0}(t) − κc_{0}(t) = 0.  (2.18)

Next, by writing out (2.17c) componentwise and by applying a solvability argument similar to the one used in (2.12), we find

\frac{ds_{0}}{dt} = −s_{0}e_{0} + (κ − λ)c_{0},  (2.19a)
\frac{de_{0}}{dt} = −(s_{0}e_{1} + s_{1}e_{0} − κc_{1}),  (2.19b)
\frac{dc_{0}}{dt} = s_{0}e_{1} + s_{1}e_{0} − κc_{1}.  (2.19c)

In turn, equation (2.19a) may be simplified using (2.18) to obtain ds_{0}/dt = −λc_{0}. In addition, (2.19b) and (2.19c) imply

e_{0}(t) + c_{0}(t) = e_{0}(0) + c_{0}(0) for all t ≥ 0,  (2.20)

where the constant is to be determined by matching with the equations on the fast (η) scale: lim_{η→∞} ũ_{0}(η) = lim_{t→0^{+}} u_{0}(t). From (2.16), we find

lim_{η→∞} ũ_{0}(η) = (1, \frac{κ}{1 + κ}, \frac{1}{1 + κ})^{T} .  (2.21)

Hence,

e_{0}(t) + c_{0}(t) = 1.  (2.22)

Finally, we combine (2.18) and (2.22) to obtain the critical manifold from the classical kinetics problem with ε = 0,

c_{0}(t) = \frac{s_{0}(t)}{s_{0}(t) + κ} and e_{0}(t) = \frac{κ}{s_{0}(t) + κ}.  (2.23)

Moreover, we see that the reduced equation for s_{0}(t) on this critical manifold is

\frac{ds_{0}}{dt} = −λ\frac{s_{0}}{s_{0} + κ} with s_{0}(0) = 1,  (2.24)

just as is the case for the pure MMH kinetics problem, see for example [15]. The solution of (2.24) is known implicitly,

s_{0}(t) + κ ln s_{0}(t) = −λt + 1.  (2.25)

Also, the rate of approach toward the slow manifold is determined by the dynamics on the fast (η) scale, cf. (2.16).

2.4. The Uniformly Valid Leading-Order Approximation. In this regime of large diffusivities, the leading-order approximation, uniformly valid in time and space, to the solution of (1.13) is obtained by combining the expressions for û_{0}, ũ_{0}, and u_{0} (recall (2.6), (2.16), (2.23), and (2.25)) and subtracting their respective common parts (recall (2.14) and (2.21)). We find

s(t, x) = s_{0}(t) + O(e^{−π^{2}δt}) + O(ε),  (2.26a)
e(t, x) = \frac{κ}{s_{0}(t) + κ} + \frac{e^{−(1+κ)\sqrt{δ}t}}{1 + κ} + O(e^{−bπ^{2}δt}) + O(ε),  (2.26b)

c(t, x) = \frac{s_{0}(t)}{s_{0}(t) + κ} − \frac{e^{−(1+κ)\sqrt{δ}t}}{1 + κ} + O(e^{−bπ^{2}δt}) + O(ε),  (2.26c)

where s_{0} is defined by (2.25), and we recall that δ = 1/ε^{2}, i.e., δ is asymptotically larger than the rate constant corresponding to the faster of the two kinetics scales.

The physical interpretation of (2.26) is that, during a short initial time interval of O(ε^{2}), diffusion effectively homogenizes the species’ concentrations. Thereafter, the concentrations of all three species are essentially uniform and independent of the fine structure of the initial distributions, and they evolve as in the classical MMH kinetics problem. On the fast time scale, enzyme rapidly binds to form complex, while in the phase space the concentrations quickly approach the slow manifold. Then, on the slow (outer) time scale, one observes the reduced kinetics; the concentrations evolve toward equilibrium along the slow manifold, with the concentrations of enzyme and complex being slaved to that of the substrate.

In the limit as δ → ∞, the expressions in (2.26) agree with the results for the chemical kinetics problem considered for example in Lin and Segel [15, Equations (14) and (15)]. Moreover, the algebraic corrections at O(ε) and upwards are independent of x, and they are governed by ordinary differential equations in t. These are obtained from solvability conditions, as is shown for u_{1}, for example, in the following subsection.

2.5. Higher-Order Corrections. On the slow time scale t, the O(ε) corrections to the leading-order solution are characterized by the O(ε^{3})–terms in (2.2),

O(ε^{3}) : −D \frac{∂^{2}u_{3}}{∂x^{2}} = −\frac{∂u_{1}}{∂t} + [F_{ε}(u)]_{2}.  (2.27)

Equation (2.17c) yields that u_{2}(t, x) = u_{2}(t); hence, application of the solvability condition to (2.27) implies

\frac{ds_{1}}{dt} = −s_{1}e_{0} − s_{0}e_{1} + (κ − λ)c_{1},
\frac{de_{1}}{dt} = −(s_{0}e_{2} + s_{1}e_{1} + s_{2}e_{0} − κc_{2}),
\frac{dc_{1}}{dt} = s_{0}e_{2} + s_{1}e_{1} + s_{2}e_{0} − κc_{2}.

Therefore, e_{1}(t) + c_{1}(t) = e_{1}(0) + c_{1}(0) = 0, since matching with the next-order approximation on the super-fast and fast scales shows that this constant is zero. Hence, e_{1} = −c_{1}; and, recalling that e_{0} = 1 − c_{0}, we see that

\frac{ds_{1}}{dt} = −(1 − c_{0})s_{1} + (κ − λ + s_{0})c_{1},

which is precisely [15, Equation (25a)]. One may proceed in a similar manner to obtain the asymptotics of u to any order, as well as the O(ε) and higher-order corrections to the slow manifold, as in the classical MMH kinetics problem.

We observe that one may also calculate the higher-order corrections to the leading-order solution on the super-fast time scale. At O(ε), system (2.3) yields

\frac{∂ŝ_{1}}{∂τ} = \frac{∂^{2}ŝ_{1}}{∂x^{2}} ,
\frac{∂ê_{1}}{∂τ} = −(ŝ_{0}ê_{0} − κĉ_{0}) + a\frac{∂^{2}ê_{1}}{∂x^{2}} ,
\frac{∂ĉ_{1}}{∂τ} = ŝ_{0}ê_{0} − κĉ_{0} + b\frac{∂^{2}ĉ_{1}}{∂x^{2}} .

Hence, substituting the leading-order solution (2.6), one sees that the constant terms, corresponding to k = 0, lead to linearly growing and decaying terms in ê_{1} and ĉ_{1}. These secular terms render the expansion invalid as an approximation on time scales of τ = O(1/ε). Therefore, one needs to use the multiple scales method (or, alternatively, the boundary function method), as follows. We let

û_{0} = û_{0}(η, τ, x) = û_{00}(η) + Σ_{k=1}^{∞}û_{0k}e^{−D(kπ)^{2}τ} cos(kπx),

so that û_{00} varies on the fast scale, while for k ≥ 1 û_{0k} remains constant, as before. Therefore, the equations at O(ε) are now

\frac{∂ŝ_{1}}{∂τ} + \frac{∂ŝ_{00}}{∂η} = \frac{∂^{2}ŝ_{1}}{∂x^{2}} ,
\frac{∂ê_{1}}{∂τ} + \frac{∂ê_{00}}{∂η} = −(ŝ_{0}ê_{0} − κĉ_{0}) + a\frac{∂^{2}ê_{1}}{∂x^{2}} ,
\frac{∂ĉ_{1}}{∂τ} + \frac{∂ĉ_{00}}{∂η} = ŝ_{0}ê_{0} − κĉ_{0} + b\frac{∂^{2}ĉ_{1}}{∂x^{2}} .

Solvability (or ‘elimination of the secular terms’) implies that one should choose the η–dependence in û_{00} so that

\frac{∂ŝ_{00}}{∂η} = 0,
\frac{∂ê_{00}}{∂η} = −(ŝ_{00}ê_{00} − κĉ_{00}),
\frac{∂ĉ_{00}}{∂η} = ŝ_{00}ê_{00} − κĉ_{00}.

By this choice, the solution of the system for û_{1} is bounded. Also, we observe that these equations are exactly the same as equations (2.13) for ũ_{0}, as expected. For an exposition of the method of multiple scales, see [14], for example.

3. Moderately Large Diffusivities

In this section, we examine the regime of moderately large diffusivities, δ = O(1/ε); for convenience, we take δ = 1/ε in (1.13). Here, the diffusive time scale is of the same order of magnitude as that of the fast kinetics.

We show that, on the fast time scale, s satisfies the homogeneous heat equation to leading order and hence approaches one exponentially in time. At the same time, to leading order, e and c satisfy inhomogeneous, linear reaction–diffusion equations, and they approach constant values exponentially in time. Moreover, the homogeneous values of e and s are, to leading order, precisely those values corresponding to the point on the slow kinetics manifold, which is to be expected.

The rate constant for the exponential convergence in time to this homogeneous state is one order of magnitude smaller than in regime (i), see Section 2.1.

On the long time scale, s is the reaction progress variable. It satisfies the classical slow reduced MMH kinetics equation. At the same time, the homogeneous enzyme and complex concentrations are slaved to that of s and lie on the critical manifold to leading order. Most significantly, the leading order dynamics in this regime turn out to be independent of the diffusivities a and b.

The equations in this regime are

ε\frac{∂s}{∂t} = −εse + ε(κ − λ)c + \frac{∂^{2}s}{∂x^{2}} ,  (3.1a)
ε\frac{∂e}{∂t} = −(se − κc) + a\frac{∂^{2}e}{∂x^{2}} ,  (3.1b)
ε\frac{∂c}{∂t} = se − κc + b\frac{∂^{2}c}{∂x^{2}} .  (3.1c)

Equivalently, in vector form, they are given by

ε\frac{∂u}{∂t} = F_{ε}(u) + D \frac{∂^{2}u}{∂x^{2}}.

3.1. Homogenization and Fast Kinetics: The Fast (Inner) Time Scale. Recall that û(η, x, ε) = u(t, x, ε) on the fast time scale given by η = t/ε. Then, the governing equations become

\frac{∂ŝ}{∂η} = −εŝê + ε(κ − λ)ĉ + \frac{∂^{2}ŝ}{∂x^{2}},  (3.2a)
\frac{∂ê}{∂η} = −(ŝê − κĉ) + a\frac{∂^{2}ê}{∂x^{2}},  (3.2b)
\frac{∂ĉ}{∂η} = ŝê − κĉ + b\frac{∂^{2}ĉ}{∂x^{2}}  (3.2c)

subject to the usual initial conditions and Neumann boundary conditions.

Making the Ansatz û(η, x, ε) = û_{0}(η, x) + εû_{1}(η, x) + O(ε^{2}), we find to lowest order that

\frac{∂ŝ_{0}}{∂η} = \frac{∂^{2}ŝ_{0}}{∂x^{2}} ,  (3.3a)
\frac{∂ê_{0}}{∂η} = −(ŝ_{0}ê_{0} − κĉ_{0}) + a\frac{∂^{2}ê_{0}}{∂x^{2}} ,  (3.3b)
\frac{∂ĉ_{0}}{∂η} = ŝ_{0}ê_{0} − κĉ_{0} + b\frac{∂^{2}ĉ_{0}}{∂x^{2}} .  (3.3c)

Hence, ŝ_{0} satisfies the homogeneous heat equation with Neumann boundary conditions and with initial data given by ŝ_{0}(0, x) = s_{i}(x), which implies that

ŝ_{0}(η, x) = ∑_{k=0}^{∞} ŝ_{0k}e^{−(kπ)^{2}η} cos (kπx),  (3.4)

where ŝ_{00} = 1 and

ŝ_{0k} = 2 ∫_{0}^{1} ŝ_{0}(0, x) cos (kπx) dx = 2 ∫_{0}^{1} s_{i}(x) cos (kπx) dx for k ≥ 1.

Next, we find the corresponding expressions for ê_{0} and ĉ_{0}. Equations (3.3b) and (3.3c) are linear in ê_{0} and ĉ_{0}, with some of the coefficients depending on η and x through s_{0}. Hence, due to the Neumann boundary conditions on ê_{0} and ĉ_{0}, one can write

ê_{0}(η, x) = ∑_{k=0}^{∞} ê_{0k}(η) cos (kπx) and ĉ_{0}(η, x) = ∑_{k=0}^{∞} ĉ_{0k}(η) cos (kπx)  (3.5)

for some sets of coefficients {ê_{0k}(η)} and {ĉ_{0k}(η)}. Substituting (3.4) and (3.5) into (3.3b) and (3.3c), making use of the identity

cos (mπx) cos (nπx) = \frac{1}{2}( cos ((m + n)πx) + cos ((m − n)πx))

as well as of ŝ_{00} = 1, and collecting coefficients of like cosines, we find that ê_{0k} and ĉ_{0k} must satisfy the following infinite system of linear ordinary differential equations:

\frac{dê_{00}}{dη} = −(ê_{00} − κĉ_{00}) − F_{0},  (3.6a)
\frac{dĉ_{00}}{dη} = ê_{00} − κĉ_{00} + F_{0},  (3.6b)

for the zeroth Fourier mode, and

\frac{dê_{0k}}{dη} = −(1 + a(kπ)^{2} + \frac{1}{2}ŝ_{0(2k)}e^{−4(kπ)^{2}t})ê_{0k} + κĉ_{0k} − F_{k},  (3.7a)
\frac{dĉ_{0k}}{dη} = −(κ + b(kπ)^{2})ĉ_{0k} + (1 + \frac{1}{2}ŝ_{0(2k)}e^{−4(kπ)^{2}t})ê_{0k} + F_{k},  (3.7b)

for the kth Fourier mode, k ≥ 1. Here, F_{k} is defined via

F_{0}(η) = \frac{1}{2} ∑_{m≥1} ŝ_{0m}e^{−(mπ)^{2}η}ê_{0m}(η),  (3.8a)

F_{k}(η) = \frac{1}{2} ∑_{\begin{matrix}m≥0 \\ m≠k \end{matrix}} (ŝ_{0(m+k)}e^{−(m+k)^{2}π^{2}η} + ŝ_{0|m−k|}e^{−(m−k)^{2}π^{2}η})ê_{0m}(η), k ≥ 1. (3.8b)

In particular, (3.6) yields that ê_{00}(η) + ĉ_{00}(η) is constant in this regime. Then, solving (3.6) and taking into account the identities

ê_{00}(0) = ∫_{0}^{1} ê_{0}(0, x) dx = ∫_{0}^{1} e_{i}(x) dx = 1,
ĉ_{00}(0) = ∫_{0}^{1} ĉ_{0}(0, x) dx = ∫_{0}^{1} c_{i}(x) dx = 0

as well as the fact that ê_{0k} and ĉ_{0k} must be bounded in η, we find

ê_{00}(η) = \frac{κ}{1 + κ} + O(e^{− min\{π^{2},1+κ\}η}), ĉ_{00}(η) = \frac{1}{1 + κ} + O(e^{− min\{π^{2},1+κ\}η}). (3.9)

Similar expressions can be derived for ê_{0k} and ĉ_{0k} with k ≥ 1; in particular, one can show that ê_{0k}, ĉ_{0k} = O(e^{−dπ^{2}η}), where we recall d = min{1, a, b}. Hence, we conclude that for η large,

û_{0}(η, x) = (1, \frac{κ}{1 + κ}, \frac{1}{1 + κ})^{T} + O(e^{− min\{dπ^{2},1+κ\}η}).  (3.10)

See Appendix B for details.

3.2. Reduced Kinetics: The Slow (Outer) Time Scale. On the slow (outer) time scale, the system dynamics are naturally described by (3.1) in the original time t. To lowest order, we find

0 = \frac{∂^{2}s_{0}}{∂x^{2}} ,  (3.11a)
0 = −(s_{0}e_{0} − κc_{0}) + a\frac{∂^{2}e_{0}}{∂x^{2}} ,  (3.11b)
0 = s_{0}e_{0} − κc_{0} + b\frac{∂^{2}c_{0}}{∂x^{2}} ,  (3.11c)

where we again make the Ansatz u(t, x, ε) = u_{0}(t, x) + εu_{1}(t, x) + O(ε^{2}).

Integrating (3.11a) with respect to x and making use of the Neumann boundary conditions, we conclude immediately that s_{0}(t, x) ≡ s_{0}(t). Similarly, it follows from either (3.11b) or (3.11c) that

s_{0}(t) ∫_{0}^{1} e_{0}(t, x) dx − κ ∫_{0}^{1} c_{0}(t, x) dx = 0.  (3.12)

To obtain explicit formulae for e_{0} and c_{0}, we rescale x via x̄ = x/\sqrt{a}, and observe that the ratio

α = a/b

is the relevant parameter, rather than a and b separately. We rewrite (3.11b) and (3.11c) as a four-dimensional linear system with t–dependent coefficients,

e′_{0} = f_{0},  (3.13a)
f′_{0} = s_{0}(t)e_{0} − κc_{0},  (3.13b)
c′_{0} = d_{0},  (3.13c)
d′_{0} = −α(s_{0}(t)e_{0} − κc_{0}),  (3.13d)

where the prime denotes (partial) differentiation with respect to x̄.

The general solution of (3.13) is given by

e_{0}(t, x̄) = κγ_{1}(t) + κγ_{2}(t)x̄ + γ_{3}(t)e^{\sqrt{s_{0}(t)+ακ}} x̄ − γ_{4}(t)e^{−\sqrt{s_{0}(t)+ακ} x̄}, (3.14a)

c_{0}(t, x̄) = s_{0}(t)γ_{1}(t) + s_{0}(t)γ_{2}(t)x̄ − αγ_{3}(t)e^{\sqrt{s_{0}(t)+ακ}} x̄ + αγ_{4}(t)e^{−\sqrt{s_{0}(t)+ακ} x̄}, (3.14b)

where f_{0} = e′_{0} and d_{0} = c′_{0} can be found from (3.14), and γ_{1}, …, γ_{4} are t–dependent constants of integration. Making use of the Neumann boundary conditions on e_{0} and c_{0} in (3.14), one sees that γ_{2}, γ_{3}, and γ_{4} must be identically zero. Hence, the only solution (3.14) that satisfies the boundary conditions is given by

(e_{0}, c_{0})(t) = γ_{1}(t)(κ, s_{0}(t)),  (3.15)

where γ_{1}(t) is as yet undetermined, and we see that (3.12) reduces to

s_{0}(t)e_{0}(t) − κc_{0}(t) = 0.  (3.16)

Summarizing, we see that e_{0} and c_{0} are spatially uniform on the slow time scale and that the constraint (3.16) is satisfied at all times.

To describe the dynamics of u_{0} on the t time scale, we consider the O(ε)–terms in (3.1),

\frac{ds_{0}}{dt} = −s_{0}e_{0} + (κ − λ)c_{0} + \frac{∂^{2}s_{1}}{∂x^{2}} ,  (3.17a)
\frac{de_{0}}{dt} = −(s_{1}e_{0} + s_{0}e_{1} − κc_{1}) + a\frac{∂^{2}e_{1}}{∂x^{2}} ,  (3.17b)
\frac{dc_{0}}{dt} = s_{1}e_{0} + s_{0}e_{1} − κc_{1} + b\frac{∂^{2}c_{1}}{∂x^{2}} .  (3.17c)

First, an integration of (3.17a) over x ∈ [0, 1] in combination with (3.16) yields

\frac{ds_{0}}{dt} = −λc_{0}.  (3.18)

To find the corresponding dynamics of e_{0} and c_{0}, we add (3.17b) and (3.17c) and integrate the resulting expression with respect to x, which gives

\frac{de_{0}}{dt} + \frac{dc_{0}}{dt} = 0.

Hence, e_{0}(t) + c_{0}(t) = e_{0}(0) + c_{0}(0) for all t ≥ 0. Finally, taking into account that lim_{t→0^{+}} u(t, x) must equal

lim_{η→∞} û_{0}(η, x) = (1, \frac{κ}{1 + κ}, \frac{1}{1 + κ})^{T} ,

we obtain

e_{0}(t) + c_{0}(t) = 1.

Combining this identity with (3.15) and solving the resulting equation for γ_{1}, we find γ_{1}(t) = (s_{0}(t) + κ)^{−1} and therefore

e_{0}(t) = \frac{κ}{s_{0}(t) + κ} and c_{0}(t) = \frac{s_{0}(t)}{s_{0}(t) + κ},  (3.19)

see (2.23). Finally, substitution of (3.19) into (3.18) yields the governing equation for s_{0},

\frac{ds_{0}}{dt} = −λ\frac{s_{0}}{s_{0} + κ} with s_{0}(0) = 1,  (3.20)

see also (2.24).

This approximation is independent of the values of a and b, and hence it coincides with the one in Section 2. The unequal diffusivities of e and c do not influence the leading-order asymptotics of (3.1).

3.3. The Uniformly Valid Leading-Order Approximation. In this regime, δ is of the same size as the rate constant corresponding to the faster of the two kinetics time scales, i.e., δ is moderately large. Therefore, given the expressions for û_{0} and u_{0} in (3.10) and (3.19), respectively, with s_{0}(t) given by (3.20), the uniformly valid leading-order approximation for u is

s(t, x) = s_{0}(t) + O(e^{− min{π^{2},1+κ}δt}) + O(ε),  (3.21a)
e(t, x) = \frac{κ}{s_{0}(t) + κ} + O(e^{− min\{aπ^{2},1+κ\}δt}) + O(ε),  (3.21b)
c(t, x) = \frac{s_{0}(t)}{s_{0}(t) + κ} + O(e^{− min\{bπ^{2},1+κ\}δt}) + O(ε),  (3.21c)

where δ = 1/ε. Again, we emphasize that this approximation is independent of the values of a and b. Correspondingly, the leading-order dynamics may be interpreted physically in the same manner as that found in the regime of large diffusivities (Section 2.4), although of course the error terms here are more significant.

3.4. Higher-Order Corrections. Similar reasoning as used in Section 3.2 can be applied to determine the asymptotics of u_{n}, for any n ≥ 1. In particular, one can show that u_{n} remains spatially uniform for all times t, u_{n} = u_{n}(t), and that the state of the system at time t is determined by the set of equations

\frac{ds_{n}(t)}{dt} = [F_{1,ε}(u(t))]_{n},  (3.22a)
e_{n}(t) + c_{n}(t) = 0,  (3.22b)
s_{0}(t)e_{n}(t) − κc_{n}(t) = ċ_{n−1} + r_{n−1}(t).  (3.22c)

Here, the overdot denotes differentiation with respect to time, F_{1,ε}(u) = −se + (κ − λ)c, and the term r_{n−1} = [F_{2,ε}(u)]_{n} + s_{0}e_{n} − κc_{n}, with F_{2,ε}(u) = −se + κc, is a function of u_{1}, …, u_{n−1}, s_{n}, and e_{0} exclusively.

The proof is by induction. First, at O(ε^{n}), (3.1a) reads

\frac{ds_{n−1}(t)}{dt} = [F_{1,ε}(u)]_{n−1}(t) + \frac{∂^{2}s_{n}(t, x)}{∂x^{2}}.

Now, ds_{n−1}/dt = [F_{1,ε}(u)]_{n−1} by the induction hypothesis, and thus ∂^{2}s_{n}/∂x^{2} = 0. Using the boundary conditions, then, we find that s_{n}(t, x) = s_{n}(t).

Next, at O(ε^{n}), (3.1b) and (3.1c) yield

ė_{n−1} − r_{n−1} = −(s_{0}e_{n} − κc_{n}) + a\frac{∂^{2}e_{n}}{∂x^{2}} ,

ċ_{n−1} + r_{n−1} = s_{0}e_{n} − κc_{n} + b\frac{∂^{2}c_{n}}{∂x^{2}} ,

Here again, we rescale x via x̄ = x/\sqrt{a} and rewrite these equations as a fourdimensional, inhomogeneous, linear system with t–dependent coefficients,

e′_{n} = f_{n},  (3.24a)
f′_{n} = s_{0}e_{n} − κc_{n} + ė_{n−1} − r_{n−1},  (3.24b)
c′_{n} = d_{n},  (3.24c)
d′_{n} = −α(s_{0}e_{n} − κc_{n} + ė_{n−1} − r_{n−1}),  (3.24d)

where the prime denotes (partial) differentiation with respect to x̄ and we have used that ė_{n−1}(t) + ċ_{n−1}(t) = 0 by the induction hypothesis.

The general solution of (3.24) is given by

\begin{aligned}e_{n}(t, x) & = κγ_{1}(t) + κγ_{2}(t)x̄ + γ_{3}(t)e^{\sqrt{s_{0}(t)+ακ}} x̄ − γ_{4}(t)e^{−\sqrt{s_{0}(t)+ακ} x̄} \\ & + \frac{ė_{n−1}(t) − r_{n−1}(t)}{s_{0}(t) + ακ}( cosh(\sqrt{s_{0}(t) + ακ} x̄) − 1), \end{aligned}  (3.25a)

\begin{aligned}c_{n}(t, x) & = s_{0}(t)γ_{1}(t) + s_{0}(t)γ_{2}(t)x̄ − αγ_{3}(t)e^{\sqrt{s_{0}(t)+ακ}} x̄ + αγ_{4}(t)e^{−\sqrt{s_{0}(t)+ακ} x̄} \\ & − α\frac{ė_{n−1}(t) − r_{n−1}(t)}{s_{0}(t) + ακ}( cosh(\sqrt{s_{0}(t) + ακ} x̄) − 1), \end{aligned} (3.25b)

where f_{n} = e′_{n} and d_{n} = c′_{n} can be found from (3.25) and γ_{1}, …, γ_{4} are t–dependent constants of integration. Making use of the Neumann boundary conditions on e_{n} and c_{n}, one sees that γ_{2}, γ_{3}, and γ_{4} are such that

e_{n}(t, x) = e_{n}(t) = κγ_{1}(t) − \frac{ė_{n−1}(t) − r_{n−1}(t)}{s_{0}(t) + ακ},  (3.26a)
c_{n}(t, x) = c_{n}(t) = s_{0}(t)γ_{1}(t) + α\frac{ė_{n−1}(t) − r_{n−1}(t)}{s_{0}(t) + ακ}.  (3.26b)

In summary, these formulas show that u_{n} = u_{n}(t).

Next, by integrating the O(ε^{n+1}) terms in (3.1a) and using the boundary conditions, we obtain ∫_{0}^{1}(ds_{n}/dt) dx = ∫_{0}^{1}[F_{1,ε}(u)]_{n} dx. To evaluate these integrals, we observe that s_{n}, and thus also ds_{n}/dt, is a function of time only. Moreover, [F_{1,ε}(u)]_{n} is a function of u_{0}, …, u_{n−1} and thus also a function of time only. Therefore, ds_{n}/dt = [F_{1,ε}(u)]_{n}, as desired.

Finally, at O(ε^{n+1}), (3.1b) and (3.1c) yield

ė_{n} = [F_{2,ε}(u)]_{n} + a\frac{∂^{2}e_{n+1}}{∂x^{2}} ,

ċ_{n} = −[F_{2,ε}(u)]_{n} + b\frac{∂^{2}c_{n+1}}{∂x^{2}} .

Integrating both members of these equations over the spatial domain [0, 1], using the boundary conditions and recalling that e_{n} and c_{n} are functions of time only, we obtain the identity ė_{n}(t) + ċ_{n}(t) = 0. Recalling also that e(t, x) + c(t, x) = 1 and that e_{0}(t) + c_{0}(t) = 1, we conclude that e_{n}(t) + c_{n}(t) = 0 as desired. Therefore, combining (3.26) with the identity ė_{n}(t) + ċ_{n}(t) = 0, we obtain (3.22), and the proof is complete.

Remark 3.1. Equations (3.22) show that the system dynamics are governed solely by the chemical kinetics together with the conservation law e(t) + c(t) = 1, to within all algebraic orders in ε. Thus, all of the results that are valid for the usual, nondiffusive MMH kinetics (such as the existence of a slow invariant manifold with an asymptotic expansion in powers of ε [12]) also apply to this case.

4. Small Diffusivities

In this section, we consider the regime in which the species’ diffusivities are small, δ = ε, so that the time scale on which diffusion acts is much longer than that of the slow kinetics.

The solution in this regime depends on the fine structure of the initial distributions, as well as on the local distributions of the species. We will show that, on the fast and slow time scales, the effects of diffusion may be neglected to a fairly good approximation (to within O(ε^{2})) and that the kinetics are largely decoupled at each point x in space. In particular, for each fixed x, the kinetics follow the classical MMH kinetics. Complex is formed on the fast time scale, with the species’ concentrations rapidly approaching an x–dependent point on the slow manifold. We label the dynamics on the fast time scale as pointwise fast kinetics. Then, the reaction proceeds on the slow time scale, with the concentrations of substrate and complex evolving along the slow manifold to zero, pointwise in x. Also, the enzyme concentration evolves to leading order toward the initial profile e_{i}(x), which still depends on x. We label the dynamics on the slow time scale as pointwise slow reduced kinetics. Finally, on the super-slow time scale, diffusion homogenizes the enzyme concentration profile.

Equations (1.13) with δ = ε are

ε\frac{∂s}{∂t} = −εse + ε(κ − λ)c + ε^{2} \frac{∂^{2}s}{∂x^{2}} ,  (4.1a)
ε\frac{∂e}{∂t} = −(se − κc) + aε^{2} \frac{∂^{2}e}{∂x^{2}} ,  (4.1b)
ε\frac{∂c}{∂t} = se − κc + bε^{2} \frac{∂^{2}c}{∂x^{2}} .  (4.1c)

In vector notation, (4.1) is given by

ε\frac{∂u}{∂t} = F_{ε}(u) + ε^{2}D \frac{∂^{2}u}{∂x^{2}}.  (4.2)

One factor of ε in (4.1a) is redundant; however, we retain it for consistency of notation.

4.1. Pointwise Fast Kinetics: The Fast (Inner) Time Scale. On the fast (η = t/ε) time scale, the full governing equations are given by

\frac{∂û}{∂η} = F_{ε}(û) + ε^{2}D \frac{∂^{2}û}{∂x^{2}},  (4.3)

subject to the usual initial and boundary conditions, see Section 1. Again, we expand û(η, x, ε) = û_{0}(η, x) + εû_{1}(η, x) + O(ε^{2}). To lowest order, (4.3) implies

O(1) : \frac{∂û_{0}}{∂η} = [F_{ε}(û)]_{0},  (4.4)

subject to Neumann boundary conditions and the initial condition

û_{0}(0, x) = u_{i}(x).  (4.5)

Solving the system (4.4) and (4.5) componentwise, we find ŝ_{0}(η, x) = ŝ_{0}(0, x); hence, ŝ_{0}(η, x) ≡ s_{i}(x) for all η ≥ 0. Moreover,

\frac{∂ê_{0}}{∂η} + \frac{∂ĉ_{0}}{∂η} = 0,

so that

ê_{0}(η, x) + ĉ_{0}(η, x) = ê_{0}(0, x) + ĉ_{0}(0, x) ≡ e_{i}(x) for all η ≥ 0.  (4.6)

Now, given that ŝ_{0} = s_{i} and ê_{0} = e_{i} − ĉ_{0}, we see from (4.4) that

\frac{∂ĉ_{0}}{∂η} = −(s_{i} + κ)ĉ_{0} + s_{i}e_{i} with ĉ_{0}(0, x) = 0,

which has the solution

ĉ_{0}(η, x) = \frac{s_{i}(x)e_{i}(x)}{s_{i}(x) + κ}(1 − e^{−(s_{i}(x)+κ)η}).

In summary, we have

û_{0}(η, x) = (s_{i}(x), \frac{e_{i}(x)}{s_{i}(x) + κ}(κ+s_{i}(x)e^{−(s_{i}(x)+κ)η}), \frac{s_{i}(x)e_{i}(x)}{s_{i}(x) + κ}(1−e^{−(s_{i}(x)+κ)η}))^{T} .

(4.7)

The physical interpretation of these results is that, on the fast time scale, diffusion has no effect on the concentration amplitudes to O(1) and O(ε); it only influences the amplitudes at O(ε^{2}). Instead, on the fast time scale, the kinetics have center stage and are essentially independent at each point x. In particular, pointwise in x and to leading order in ε, ŝ remains unaltered, and ê and ĉ are constrained to satisfy the linear conservation law (4.6), while the concentrations evolve to the appropriate (x–dependent) point on the leading-order slow manifold as η → ∞, which is exactly dictated by the chemical kinetics.

4.2. Pointwise Slow Reduced Kinetics: The Slow Time Scale. The slow time scale is given by the original time scale t itself; hence, with a slight abuse of notation to ensure consistency with Section 2, we replace u by ũ = ũ(t, x, ε) in (4.2) and make the Ansatz ũ(t, x, ε) = ũ_{0}(t, x)+εũ_{1}(t, x)+O(ε^{2}). After expanding with respect to ε and rearranging, we have

O(1) : [F_{ε}(ũ)]_{0} = 0,  (4.8a)
O(ε) : \frac{∂ũ_{0}}{∂t} = [F_{ε}(ũ)]_{1},  (4.8b)

with Neumann boundary conditions on ũ_{0} and ũ_{1} and with initial conditions to be determined by matching with the equations on the fast time scale.

Writing (4.8a) componentwise, we see that s̃_{0}ẽ_{0} = κc̃_{0}, which we can substitute into (4.8b) to obtain

\frac{∂s̃_{0}}{∂t} = −λc̃_{0}.

Also, we deduce from (4.8b) that ẽ_{0}(t, x) + c̃_{0}(t, x) = ẽ_{0}(0) + c̃_{0}(0) must hold. Imposing the matching condition lim_{η→∞} û_{0}(η, x) = lim_{t→0^{+}} ũ_{0}(t, x) and taking into consideration that

lim_{η→∞} û_{0}(η, x) = (s_{i}(x), \frac{κe_{i}(x)}{s_{i}(x) + κ}, \frac{s_{i}(x)e_{i}(x)}{s_{i}(x) + κ})^{T} ,  (4.9)

we conclude that ẽ_{0}(t, x) + c̃_{0}(t, x) ≡ e_{i}(x). Therefore, recalling (4.8a), we obtain

ẽ_{0}(t, x) = \frac{κe_{i}(x)}{s̃_{0}(t, x) + κ} and c̃_{0}(t, x) = \frac{s̃_{0}(t, x)e_{i}(x)}{s̃_{0}(t, x) + κ},  (4.10)

where s̃_{0}(t, x) solves

\frac{∂s̃_{0}}{∂t} = −λ\frac{e_{i}(x)}{s̃_{0} + κ}s̃_{0} with s̃_{0}(0, x) = s_{i}(x).  (4.11)

This reduced equation for s̃_{0}(t, x) is of the same form as (2.24), the equation for s_{0}(t) in the large diffusivities regime; however, the interpretations differ. Here, s̃_{0}(t, x) depends on x through e_{i}(x), whereas s_{0}(t) is independent of x. Moreover, we recall that s̃_{0} may be found only implicitly, as was shown for s_{0}(t) in (2.25).

Therefore, we conclude that, at each point x in the domain, s̃_{0} is the natural reaction progress variable and that the enzyme and complex concentrations are given as functions of it, as in the chemical kinetics problem. However, we emphasize that, in this small-diffusivity regime, the sum of ẽ_{0} and c̃_{0} depends on x and is equal to the local initial enzyme profile e_{i}(x), not uniformly equal to one, as in the classical kinetics problem. This means, among other things, that, for certain spatial profiles e_{i}(x), the dynamic evolution of ũ takes place outside the unit cube, i.e., outside the region that is feasible for the classical kinetics problem.

4.3. Homogenization of the Enzyme Concentration: The Super-Slow (Outer) Time Scale. Due to the fact that δ = ε in this regime, diffusion acts on the super-slow (ζ = εt) time scale, and we express the governing equations as

ε^{2} \frac{∂u}{∂ζ} = F_{ε}(u) + ε^{2}D \frac{∂^{2}u}{∂x^{2}},  (4.12)

subject to Neumann boundary conditions. We make the Ansatz

u(ζ, x, ε) = u_{0}(ζ, x) + εu_{1}(ζ, x) + ε^{2}u_{2}(ζ, x) + O(ε^{3}),

substitute it into (4.12), expand, and rearrange the resulting equations to obtain

O(1) : [F_{ε}(u)]_{0} = 0,  (4.13a)
O(ε) : [F_{ε}(u)]_{1} = 0,  (4.13b)
O(ε^{2}) : L_{ζ}u_{0} = [F_{ε}(u)]_{2},  (4.13c)

where L_{ζ} = \frac{∂}{∂ζ} − D \frac{∂^{2}}{∂x^{2}} , supplemented by Neumann boundary conditions.

First, observe that (4.13a) gives s_{0}e_{0} = κc_{0}, as on the slow (t) scale. Then (4.13b) shows that λc_{0} = 0, which implies c_{0} ≡ 0 (since λ ≠ 0 by definition). Moreover, we see from that same equation (4.13b) that

s_{0}e_{1} + s_{1}e_{0} = κc_{1},

which we use in (4.13c) to obtain

\frac{∂s_{0}}{∂ζ} − \frac{∂^{2}s_{0}}{∂x^{2}} = −λc_{1}.

Similarly, to find e_{0}, we note that (4.13c) and the identity c_{0} ≡ 0 derived above imply L_{ζ}c_{0} = 0; hence,

s_{0}e_{2} + s_{1}e_{1} + s_{2}e_{0} = κc_{2}.

Therefore, L_{ζ}e_{0} = 0 and

e_{0}(ζ, x) = ∑_{k=0}^{∞} e_{0k}e^{−a(kπ)^{2}ζ} cos (kπx),

with the constant coefficients e_{0k} to be determined by matching with the dynamics on the t scale. Since matching requires lim_{t→∞} ũ_{0}(t, x) = lim_{ζ→0^{+}} u_{0}(ζ, x) and

lim_{t→∞} ũ_{0}(t, x) = (0, e_{i}(x), 0)^{T}  (4.14)

by (4.10) and (4.11), it follows that e_{0} solves the homogeneous heat equation with nonzero initial conditions and with Neumann boundary conditions. Hence, the enzyme concentration will generically be nonzero. More precisely, since

e_{00} = ∫_{0}^{1} e_{i}(x) dx = 1,

it follows that e_{0}(ζ, x) = 1 + O(e^{−aπ^{2}ζ}). All reactions have to leading order been completed when diffusion sets in, i.e., there is only enzyme left to diffuse until a spatially homogeneous distribution of e has been reached.

Finally, given e_{0}(x) ≢ 0 and (4.13a), we conclude that s_{0} ≡ 0 must hold. Hence, to summarize, we have

u_{0}(ζ, x) = (0, 1, 0)^{T} + O(e^{−dπ^{2}ζ}),  (4.15)

where again d = min{1, a, b}.

4.4. The Uniformly Valid Leading-Order Approximation. We combine the expressions for û_{0}, ũ_{0}, and u_{0}, which are the leading-order approximations to the solution u of (4.2) on the fast, slow, and super-slow time scales, respectively, into one uniformly valid, leading-order approximation for u. Given (4.7), (4.10) and (4.11), and (4.15), as well as (4.9) and (4.14), a straightforward calculation shows

s(t, x) = s̃_{0}(t, x) + O(e^{−π^{2}δt}) + O(ε),  (4.16a)

e(t, x) = 1 − \frac{s̃_{0}(t, x)e_{i}(x)}{s̃_{0}(t, x) + κ} + \frac{s_{i}(x)e_{i}(x)}{s_{i}(x) + κ}e^{−(s_{i}(x)+κ)t/δ} + O(e^{−aπ^{2}δt}) + O(ε), (4.16b)

c(t, x) = \frac{s̃_{0}(t, x)e_{i}(x)}{s̃_{0}(t, x) + κ} − \frac{s_{i}(x)e_{i}(x)}{s_{i}(x) + κ}e^{−(s_{i}(x)+κ)t/δ} + O(e^{−bπ^{2}δt}) + O(ε), (4.16c)

where s̃_{0} is defined by (4.11). Here, we also recall that δ = ε, i.e., δ is asymptotically smaller than the smaller of the two kinetics rate constants.

4.5. Higher-Order Corrections. In this regime, the higher-order algebraic corrections u_{j}, j ≥ 1, may be found explicitly.

The O(ε) terms in (4.3) are

O(ε) : \frac{∂û_{1}}{∂η} = [F_{ε}(û)]_{1}.  (4.17)

In particular, the equation for the first component, ŝ_{1}, is

\frac{∂ŝ_{1}}{∂η} = −ŝ_{0}ê_{0} + (κ − λ)ĉ_{0} = \frac{−s_{i}(x)e_{i}(x)}{s_{i}(x) + κ}(λ + e^{−(s_{i}(x)+κ)η}(s_{i}(x) + κ − λ)), (4.18)

where we used (4.7). Hence, there will be a secular term in ŝ_{1} that grows linearly in η, due to the η–independent (first) term in (4.18).

A remedy is at hand using the method of multiple scales. In particular, let the terms in the expansion of û depend not only on the fast time η, but also on the slow time t and on a sequence of successively longer time scales as needed, depending on the order of the truncation of the asymptotic expansion. We start with the O(1) term, A_{0}(t, x)ŝ_{0}(η, x), in the expansion for ŝ. As in Section 4.1, it is again the case that ŝ_{0} = s_{i}(x), independently of η, and that (4.6) holds. Moreover, the leading-order terms for the complex and enzyme concentrations may be found as in Section 4.1. For example, to leading order, the complex concentration is

\frac{A_{0}(t, x)s_{i}(x)e_{i}(x)}{A_{0}(t, x)s_{i}(x) + κ}(1 − e^{−(A_{0}(t,x)s_{i}(x)+κ)η}).

Substituting these solutions into (4.3) and collecting the O(ε) terms, we find

s_{i}(x)\frac{∂A_{0}}{∂t} +A_{1} \frac{∂ŝ_{1}}{∂η} = −\frac{s_{i}(x)e_{i}(x)A_{0}}{A_{0}s_{i}(x) + κ} (λ+(A_{0}s_{i}(x)+κ−λ)e^{−(A_{0}s_{i}(x)+κ)η}). (4.19)

Therefore, if we choose A_{0}(t, x) such that the η–independent term vanishes from the right hand side of (4.19), then we obtain a regular (bounded) solution for ŝ_{1}. Specifically, we choose the dependence of A_{0}(t, x) on t so that

\frac{∂A_{0}}{∂t} = −λ\frac{e_{i}(x)A_{0}}{A_{0}s_{i}(x) + κ}.  (4.20)

This requirement is equivalent to applying a solvability condition to equation (4.19), which forces the ‘bad’ terms to drop out. The resulting equation for ŝ_{1} may be integrated directly to obtain a bounded function. Moreover, as one expects, the equation (4.20) for A_{0}(t, x) is precisely equation (4.11), which was obtained from the leading-order analysis on the slow time scale, with s̃_{0}(t, x) replaced by A_{0}(t, x)s_{i}(x).

Higher-order terms ŝ_{j}, j ≥ 1, may be found in an analogous manner. The corresponding amplitude functions A_{i} (for i = 1, …, j) depend on t, εt, …, ε^{j−1}t, as well as on x, and are determined by solvability conditions, just as A_{0}(t, x) was determined by the solvability condition on the equation for ŝ_{1} above. Of course, the corrections to the complex and enzyme concentrations need to be determined simultaneously using the multiple scales Ansatz.

5. Interpretation and Conclusions

In the present article, we have studied the effects of diffusion on the classical MMH reaction mechanism. We have investigated the two extreme regimes of large diffusivities and of small diffusivities, as well as an intermediate regime of moderately large diffusivities. A brief summary of the observed dynamics and time scales in each regime was given in Table 1. Here, we summarize our findings in more detail and illustrate the analytical results with direct numerical simulations.

Regime (i): large diffusivities, δ = 1/ε^{2}. There is a brief initial transient period of O(ε^{2}) during which the concentrations become spatially homogeneous, with values equal to the spatial averages of the initial concentration profiles. This period is followed by another short, O(ε)–period during which the fast kinetics act. During this second period, the homogeneous concentrations everywhere in the domain are brought into the regime where the MMH reduction holds. Then, for large times, the dynamics are exactly the same as those found in the pure kinetics problem, as described above. The leading-order asymptotics, uniformly valid in time and space, are given in (2.26).

The results of a typical simulation of (1.13) for this regime are shown in Figure 1. In particular, the super-fast initial transient appears to last for essentially only one timestep, as is evident from Figures 1(a) and 1(b).

Regime (ii): moderately large diffusivities, δ = 1/ε. The substrate concentration again homogenizes during a brief initial transient, which is of O(ε) here. However, since the fast reaction kinetics act on the same scale as diffusion, there is an interplay between the two on this fast time scale. Then, as the species’ concentrations homogenize, the kinetics start to take over until an MMH reduction is again applicable. The uniformly valid approximation is given by (3.21).

The results of a typical simulation of (1.13) are shown in Figure 2. The numerical outcome is qualitatively similar to that for regime (i); however, homogenization is one order of magnitude slower, since diffusion acts on a time scale of O(ε), as compared to O(ε^{2}) before. See especially Figure 2(b), which suggests that the number of timesteps is approximately five.

Regime (iii): small diffusivities, δ = ε. At every point in the domain, the kinetics dominate on fast and slow time intervals on the orders of ε and 1, respectively. In other words, at every point in space, the reaction proceeds exactly as in the pure kinetics problem, with no communication to leading order between neighboring points. Then, after the substrate has been used up and all of the complex has decayed into product and enzyme at every point, on the super-slow time scale, the enzyme profile becomes homogenized due to diffusion, ultimately

(a)  (b)

(c)

Figure 1. Regime (i): large diffusivities. Concentration profiles of (a) s, (b) e, and (c) c as functions of t ≥ 0 and x ∈ [0, 1], with parameter values κ = 3, λ = 5, ε = 0.1, δ = 100, a = 2, b = 3, and initial profiles (s, e, c)(0, x) = (1 + \frac{1}{2} cos(2πx), 1 + \frac{1}{2} cos(πx), 0).

converging to the spatial average of the initial enzyme concentration profile. The leading-order approximation, valid uniformly in time and space, is given in (4.16).

The results of a typical simulation of (1.13) for this regime are shown in Figure 3. It takes at least ten timesteps before the concentrations are essentially homogeneous.

The above case studies are of interest in a variety of applications modeled by the MMH mechanism. One interesting example is the problem of nutrient uptake in cells [5], see also Section 1. There, the role of the enzyme is played by the unoccupied receptor sites on the outer boundary of the cell membrane, and the role of the substrate is played by nutrients outside the cell. When a nutrient molecule binds to a receptor site, the receptor site is said to be occupied; this occupied site plays the role of the complex. Moreover, this binding is a reversible reaction; after binding, nutrient molecules may be transported into the cell, which renders the binding site again unoccupied, i.e., the “enzyme” is recycled. Naturally, nutrient inside the cell corresponds to the MMH product. In sum, the classical MMH mechanism may be used to model this type of nutrient uptake process.

Moreover, this class of problems constitutes an interesting case study in which the diffusivities of “enzyme” and “complex” are identical, since the receptor sites (occupied and unoccupied) are attached to the cell membrane and hence move with

(a)  (b)

(c)

Figure 2. Regime (ii): moderately large diffusivities. Concentration profiles of (a) s, (b) e, and (c) c as functions of t ≥ 0 and x ∈ [0, 1], with parameter values κ = 3, λ = 5, ε = 0.1, δ = 10, a = 2, b = 0.5, and initial profiles (s, e, c)(0, x) = (1 + \frac{1}{2} cos(2πx), 1 + \frac{1}{2} cos(πx), 0).

the cell. Preliminary results suggest that at least some of the more complex regimes of mixed diffusivities can be analyzed in a manner similar to that outlined above.

The investigation of various mixed regimes and of other related reaction mechanisms, including catalysis [2], as well as the study of the effects of Dirichlet and mixed boundary conditions (which give rise to time-dependent boundary layers) constitute possible topics for future research.

Appendix A. Linear stability of the extinguished state

In this appendix, we state and prove the linear (spectral) stability of the spatially uniform (homogeneous) equilibria u^{∗} = (0, e^{∗}, 0) to which solutions tend in the limit t → ∞, i.e., after all reactions have been completed and the effect of diffusion has been accounted for. We show that these equilibria are (linearly) stable with respect to small inhomogeneous perturbations, irrespective of the relative time scales of reaction and diffusion. Our approach relies on a classical Turing stability analysis [28]; for more general considerations on the stability of homogeneous solutions in the presence of diffusion, we refer to [3, 18].

Lemma A.1. For any e^{∗} ≥ 0 constant, the homogeneous equilibrium solution u^{∗} = (0, e^{∗}, 0) of (1.16) is linearly stable.

(a)  (b)

(c)

Figure 3. Regime (iii): small diffusivities. Concentration profiles of (a) s, (b) e, and (c) c as functions of t ≥ 0 and x ∈ [0, 1], with parameter values κ = 3, λ = 5, ε = 0.1, δ = 0.1, a = 2, b = 3, and initial profiles (s, e, c)(0, x) = (1 + \frac{1}{2} cos(2πx), 1 + \frac{1}{2} cos(πx), 0).

Proof. The argument follows closely the corresponding discussion of morphogenesis in [5, Section 11.4]. First, we linearize (1.13) about u = u^{∗}. Let u′ = u − u^{∗}; then,

\frac{∂u′}{∂t} = \frac{1}{ε}dF_{ε}(u^{∗})u′ + δD \frac{∂^{2}u′}{∂x^{2}} ,  (A.1)

where D = diag(1, a, b) is defined as above and dF_{ε}(u^{∗}) is given by

dF_{ε}(u^{∗}) = [\begin{matrix}−εe^{∗} & 0 & ε(κ − λ) \\ −e^{∗} & 0 & κ \\ e^{∗} & 0 & −κ \end{matrix}] .

Now, we make the Ansatz u′ = (s′, e′, c′)^{T} = (α_{1}, α_{2}, α_{3})^{T} cos (kπx)e^{σt}, which we substitute into (A.1). After dividing out the common factor cos (kπx)e^{σt} and rearranging the resulting equations, we obtain the following homogeneous system of linear equations for (α_{1}, α_{2}, α_{3})^{T},

(σI_{3} − A) · (α_{1}, α_{2}, α_{3})^{T} = (0, 0, 0)^{T} ,  (A.2)

where A is defined as A = dF_{ε}(u^{∗})/ε − δ(kπ)^{2}D and I_{n} denotes the n × n-identity matrix. For (A.2) to have nontrivial solutions, det (σI_{3} − A) must vanish. An elementary calculation shows

det (σI_{3} − A) = (σ + δa(kπ)^{2})[σ^{2} + (e^{∗} + δ(kπ)^{2} + \frac{κ}{ε} + δb(kπ)^{2})σ

+ (e^{∗} + δ(kπ)^{2})(\frac{κ}{ε} + δb(kπ)^{2}) − \frac{e^{∗}}{ε} (κ − λ)],

which implies that either σ = −δa(kπ)^{2} < 0, or that the term in square brackets vanishes. This term, however, is precisely the determinant of the 2 × 2-submatrix σI_{2} − Ā of σI_{3} − A obtained by deleting the second row and column, with

Ā = (\begin{matrix}−(e^{∗} + δ(kπ)^{2}) & κ − λ \\ \frac{e^{∗}}{ε} & −(\frac{κ}{ε} + δb(kπ)^{2}) \end{matrix}) .  (A.3)

Hence, we are within the framework of [5], and it remains to show that tr Ā < 0 and det Ā > 0 for Re(σ) < 0 to hold, as clearly

σ = \frac{tr Ā}{2} ± \sqrt{(\frac{tr Ā}{2})^{2} − det Ā}.

However, this is immediate from (A.3), since

tr Ā = −(e^{∗} + \frac{κ}{ε} + δ(kπ)^{2}(1 + b))

and

det Ā = \frac{1}{ε}(κδ(kπ)^{2} + λ) + δb(kπ)^{2}(e^{∗} + δ(kπ)^{2}),

and since the values of all individual parameters in these expressions are nonnegative. □

Remark A.2. The above conclusion also holds for κ = 0, since tr Ā remains negative and det Ā positive.

Appendix B. Estimates on (3.6) and (3.7) for regime (ii)

In this appendix, we present some of the technical steps involved in analyzing the dynamics of ê_{0k} and ĉ_{0k}, (3.6) and (3.7), and of deriving (3.10), for regime (ii).

As shown in Section 3, we know that ê_{00}(η) + ĉ_{00}(η) = 1 for all η > 0. Thus, the equation (3.6) for ê_{00} may be rewritten as

\frac{d}{dη} (ê_{00} − \frac{κ}{κ + 1}) = −(1 + κ)(ê_{00} − \frac{κ}{κ + 1}) − F_{0}.

Writing w = ê_{00} − κ/(κ + 1) and multiplying both members by w, we obtain

\frac{1}{2} \frac{d(w^{2})}{dη} = −(1 + κ)w^{2} − F_{0}w.

Applying Young’s inequality to the last term in the right member, we find

\frac{1}{2} \frac{d(w^{2})}{dη} ≤ −(1 + κ − σ)w^{2} + \frac{1}{4σ} F_{0}^{2}.

(Here, σ > 0 is a suitably chosen parameter.) Next, application of Gronwall’s inequality yields

(w(η))^{2} ≤ (w(0))^{2}e^{−2(1+κ−σ)η} + \frac{1}{2σ} ∫_{0}^{η} e^{−2(1+κ−σ)(η−s)}(F_{0}(s))^{2} ds.  (B.1)

Recalling the definition of F_{0} and that ŝ_{0} = s_{i}, applying Hölder’s inequality twice, and employing Parseval’s identity, we estimate

(F_{0}(η))^{2} = \frac{1}{2}e^{−2π^{2}η}( ∑_{m≥1} ŝ_{0m}ê_{0m}(η))^{2} ≤ \frac{1}{2}e^{−2π^{2}η}∥s_{i}∥_{2}^{2}∥ê_{0}(η)∥_{2}^{2}.  (B.2)

Substituting this estimate into the integral in (B.1), we calculate

(w(η))^{2} ≤ ((w(0))^{2} − C(σ)) e^{−2(1+κ−σ)η} + C(σ)e^{−2π^{2}η},

where

C(σ) = \frac{∥s_{i}∥_{2}^{2} (sup_{η≥0} ∥ê_{0}(η)∥_{2}^{2})}{8σ(1 + κ − σ − π^{2})}.

Here, ∥ê_{0}(η)∥_{2} is guaranteed to remain bounded uniformly in time by standard results. This proves that w(η) → 0 (or, equivalently, that ê_{00}(η) → κ/(κ + 1)) as η → ∞ at an exponential rate.

The same type of estimates can be made to determine the long-term evolution of the kth Fourier modes ê_{0k}(η) and ĉ_{0k}(η). In particular, one can show that these modes decay to zero exponentially fast. Hence, we have demonstrated (3.10).

Acknowledgments. The authors of this manuscript thank Mike Davis (Argonne National Laboratory) and Gene Wayne (Boston University) for stimulating conversations.

The research of H. Kaper was funded in part by Award No. DMS-0549430-001 from the National Science Foundation, and by Contract DE-AC02-06CH11357 from the US Department of Energy.

The research of T. Kaper was supported by NSF grant DMS-0606343. T. Kaper thanks the Department of Mathematical Sciences at the University of Montana and the CWI for their hospitality.

The research of N. Popović was supported by NSF grant DMS-0109427. The research of A. Zagaris was supported by NWO grant 639.031.617.

References

[1] J. R. Bowen, A. Acrivos, and A. K. Oppenheim; Singular perturbation refinement to the quasi-steady state approximation in chemical kinetics, Chem. Eng. Sci. 18 (1963), 177–188.

[2] M. Boudart and G. Djéga-Mariadassou; Kinetics of Heterogeneous Catalytic Reactions, Princeton University Press, Princeton, 1984.

[3] P. De Leenheer, D. Angeli, and E. D. Sontag; Monotone chemical reaction networks, DIMACS Technical Report 16 (2004).

[4] M. Dixon and E. C. Webb; Enzymes, third edition, Academic Press, New York, 1979.

[5] L. Edelstein-Keshet; Mathematical models in biology, Classics in Applied Mathematics 46, Society for Industrial and Applied Mathematics, Philadelphia, 2005. Reprint of the 1988 original.

[6] B. P. English, W. Min, A. M. van Oijen, K. T. Lee, G. Luo, H. Sun, B.J. Cherayil, S. C. Kou, and X. S. Xie; Ever-fluctuating single enzyme molecules: Michaelis-Menten equation revisited, Nature Chemical Biology 2(2) (2006), 87–94.

[7] J. E. Ferrell, Jr. and R. R. Bhatt; Mechanistic studies of the dual phosphorylation of mitogenactivated protein kinase, J. Biol. Chem. 272(30) (1997), 19008–19016.

[8] G. G. Hammes; Enzyme Catalysis and Regulation, Academic Press, New York, 1982.

[9] F. G. Heineken, H. M. Tsuchiya, and R. Aris; On the mathematical status of the pseudosteady state hypothesis of biochemical kinetics, Math. Biosci. 1 (1967), 95–113.

[10] V. Henri; Lois générales de l’action des diastases, Hermann, Paris, 1903.

[11] C. K. R. T. Jones; Geometric Singular Perturbation Theory pages 44–118, in Dynamical Systems, Montecatini Terme, 1994, R. Johnson, ed., LNM 1609, Springer, Berlin, 1994.

[12] T. J. Kaper; An introduction to geometric methods and dynamical systems theory for singular perturbation problems, pages 85–132, in Analyzing Multiscale Phenomena Using Singular Perturbation Methods, eds. J. Cronin and R.E. O’Malley, Jr., Proc. Symp. Appl. Math, 56, American Mathematical Society, Providence, RI, 1999.

[13] J. Keener and J. Sneyd; Mathematical Physiology, Interdisciplinary Applied Mathematics 8, Springer-Verlag, New York, 1998.

[14] J. Kevorkian and J. D. Cole; Multiscale and Singular Perturbation Methods Applied Mathematical Sciences 114, Springer-Verlag, New York, 1996.

[15] C. C. Lin and L. A. Segel; Mathematics Applied to Deterministic Problems in the Natural Sciences. With material on elasticity by G.H. Handelman, Macmillan Publishing Co., Inc., New York, 1974.

[16] L. Michaelis and M. L. Menten; Die Kinetik der Invertinwirkung, Biochem. Zeitsch. 49 (1913), 333–369.

[17] P. Miller, A. M. Zhabotinsky, J. E. Lisman, and X.-J. Wang; The stability of a stochastic CaMKII Switch: Dependence on the number of enzyme molecules and protein turnover, PLoS Biol. 3(4) (2005), e107, 705–717.

[18] M. Mincheva and D. Siegel; Stability of mass action reaction–diffusion systems, Nonlinear Anal. 56(8) (2004), 1105–1131.

[19] J. D. Murray; Mathematical Biology, Biomathematics 19, Springer-Verlag, Berlin, 1989.

[20] R. E. O’Malley, Jr.; Singular Perturbation Methods for Ordinary Differential Equations, Applied Mathematical Sciences 89, Springer-Verlag, New York, 1991.

[21] B. O. Palsson; On the dynamics of the irreversible Michaelis-Menten reaction mechanism, Chem. Eng. Sci. 42(3) (1987), 447–458.

[22] B. O. Palsson and E. N. Lightfoot; Mathematical Modeling of Dynamics and Control in Metabolic Networks. I. On Michaelis-Menten Kinetics, J. Theor. Biol. 111 (1984), 273–302.

[23] S. Schnell and P. K. Maini; A Century of Enzyme Kinetics: Reliability of the K_{M} and v_{max} Estimates, Comm. Theoret. Biol. 8 (2003), 169–187.

[24] I. H. Segel; Enzyme Kinetics: Analysis of Rapid Equilibrium and Steady-State Enzyme Systems, Wiley and Sons, New York, 1975.

[25] M. Stiefenhofer; Quasi-steady-state approximation for chemical reaction networks, J. Math. Biol. 36 (1998), 593–609.

[26] S. Strickland, G. Palmer and V. Massey; Determination of disassociation constants and specific rate constants of enzyme-substrate (or protein-ligand) interactions from rapid reaction kinetic data, J. Biol. Chem. 250(11) (1975), 4048–4052.

[27] L. A. Segel and M. Slemrod; The quasi-steady state assumption: A case study in perturbation, SIAM Review 31 (1989), 446–477.

[28] A. M. Turing; The chemical basis of morphogenesis, Phil. Trans. Roy. Soc. Lond. Series B 237 (1952), 6220–6249.

[29] A. B. Vasil’eva, V. F. Butuzov, and L. V. Kalachev; The boundary function method for singular perturbation problems. With a foreword by Robert E. O’Malley, Jr., SIAM Studies in Applied Mathematics 14, Society for Industrial and Applied Mathematics, Philadelphia, 1995.

[30] Y. Wei and M. H. Hecht; Enzyme-like proteins from an unselected library of designed amino acid sequences, Protein Engineering, Design and Selection 17(1), 67–75, 2004.

[31] A. N. Yannacopoulos, A. S. Tomlin, J. Brindley, J. H. Merkin, and M.J. Pilling; The use of algebraic sets in the approximation of inertial manifolds and lumping in chemical kinetics, Physica D 83, 421–449, 1995.

Leonid V. Kalachev Department of Mathematical Sciences, University of Montana, Missoula, MT 59812, USA

E-mail address: kalachevl@mso.umt.edu

Hans G. Kaper Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, IL 60439, USA

Current address: Division of Mathematical Sciences, National Science Foundation, Arlington, VA 22230, USA

E-mail address: hkaper@nsf.gov

Tasso J. Kaper Department of Mathematics and Statistics and Center for BioDynamics, Boston University, Boston, MA 02215, USA

E-mail address: tasso@math.bu.edu

Nikola Popović Department of Mathematics and Statistics and Center for BioDynamics, Boston University, Boston, MA 02215, USA

E-mail address: popovic@math.bu.edu

Antonios Zagaris Korteweg-de Vries Institute, University of Amsterdam, 1018 TV Amsterdam, The Netherlands

Current address: Modelling, Analysis and Simulation, Centrum voor Wiskunde en Informatica (CWI), 1090 GB Amsterdam, The Netherlands

E-mail address: a.zagaris@cwi.nl