\documentclass[reqno]{amsart}
\usepackage{subfigure,graphicx}
\usepackage{hyperref}

\AtBeginDocument{ {\noindent\small
2006 International Conference in Honor of Jacqueline Fleckinger.
\newline {\em Electronic Journal of Differential Equations},
Conference 16, 2007, pp. 155--184.
\newline ISSN: 1072-6691.
URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu (login: ftp)}
\thanks{\copyright 2007 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document} \setcounter{page}{155}

\title[\hfilneg EJDE/Conf/16 \hfil
       Reduction for Michaelis-Menten-Henri kinetics]
{Reduction for Michaelis-Menten-Henri kinetics
 in the presence of diffusion}

\author[L.~Kalachev, H.~Kaper, T.~Kaper, N.~Popovi\'c, A.~Zagaris
\hfil EJDE/Conf/16 \hfilneg]
{Leonid V.~Kalachev, Hans G.~Kaper, Tasso J.~Kaper,\\
 Nikola Popovi\'c, Antonios Zagaris} % in alphabetical order

\address{Leonid V.~Kalachev \newline
Department of Mathematical Sciences, University of Montana,
Missoula, MT 59812, USA}
\email{kalachevl@mso.umt.edu}

\address{Hans G. Kaper \newline
Mathematics and Computer Science Division,
Argonne National Laboratory, Argonne, IL 60439, USA}
\curraddr{Division of Mathematical Sciences, National Science Foundation,
Arlington, VA 22230, USA}
\email{hkaper@nsf.gov}

\address{Tasso J. Kaper \newline
Department of Mathematics and Statistics and Center for BioDynamics,
Boston University, Boston, MA 02215, USA}
\email{tasso@math.bu.edu}

\address{Nikola Popovi\'c \newline
Department of Mathematics and Statistics and Center for BioDynamics,
Boston University, Boston, MA 02215, USA}
\email{popovic@math.bu.edu}

\address{Antonios Zagaris \newline
Korteweg-de Vries Institute, University of Amsterdam,
1018 TV Amsterdam, The Netherlands}
\curraddr{Modelling, Analysis and Simulation,
          Centrum voor Wiskunde en Informatica (CWI),
1090 GB Amsterdam, The Netherlands}
\email{a.zagaris@cwi.nl}

\thanks{Published May 15, 2007.}
\subjclass[2000]{35K57, 35B40, 92C45, 41A60}
\keywords{Michaelis-Menten-Henri mechanism; diffusion;
        dimension reduction; \hfill\break\indent matched asymptotics}

\dedicatory{Dedicated to Jacqueline Fleckinger on her 65-th birthday}

\begin{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.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{example}[theorem]{Example}
\newtheorem{xca}[theorem]{Exercise}
\newtheorem{remark}[theorem]{Remark}
\newcommand{\abs}[1]{|#1|}

\section{Introduction}\label{section1}

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~\cite{H1903,MM1913}.
The reaction mechanism is represented symbolically by
\begin{equation}\label{1a}
S+E\overset{k_1}{\underset{k_{-1}}{\rightleftharpoons}}C
\stackrel{k_2}{\rightarrow} P+E,
\end{equation}
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 \cite{DW1979,E2006,FB1997,H1982,MZLW2005,P1987,PL1984,SM2003,S1975,SPM1975,WH2004},
and examples
involving nutrient uptake in cells
and heterogeneous
catalytic reactions
are analyzed
in~\cite[Chapter~7.1]{EK2005}
and~\cite{BDM1984}, 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, \emph{e.g.}, \cite{KS1998,LS1974,M1989,OM1991}.

\subsection{MMH Kinetics}
\label{section1.1}

Equations governing the kinetics
of \eqref{1a} may be derived
from the law of mass action,
\begin{subequations}\label{1b}
\begin{align}
\frac{dS}{d\tilde{t}} &=-k_1SE+k_{-1}C, \\
\frac{dE}{d\tilde{t}} &=-k_1SE+(k_{-1}+k_2)C, \\
\frac{dC}{d\tilde{t}} &=k_1SE-(k_{-1}+k_2)C, \\
\frac{dP}{d\tilde{t}} &=k_2C,
\end{align}
\end{subequations}
where $S$, $E$, $C$, and $P$ denote the concentrations of substrate,
enzyme,
complex, and product, respectively.
The initial conditions
are
$S(0)={\bar S},
E(0)={\bar E},
C(0)=0$,
and $P(0)=0$.

Dimension reduction
in \eqref{1b}
is traditionally achieved
in two steps.
The first step uses
the pair
of conservation relations
that exist
for the mechanism
\eqref{1a}.
In particular,
the total concentration
of enzyme (free
and bound in complex) is constant;
that is,
$E(t)+C(t)={\bar E}$ for all $t>0$.
In addition,
$S(t)+C(t)+P(t)={\bar S}$ for all $t>0$.
Therefore,
there is a decrease
from four variables
to two,
and the governing equations \eqref{1b}
reduce to
\begin{subequations}\label{1btwo}
\begin{gather}
\frac{dS}{d\tilde{t}} =-k_1\bar{E}S+(k_1{S}+ k_{-1})C, \\
\frac{dC}{d\tilde{t}} =k_1\bar{E}S-(k_1{S}+k_{-1}+k_2)C.
\end{gather}
\end{subequations}

The second reduction step exploits
the separation of time scales.
In particular,
${\bar E} \ll {\bar S}$.
Hence, there is
a small, positive,
dimensionless parameter,
\begin{equation}\label{def-eps}
\varepsilon=\frac{\bar{E}}{\bar{S}},
\end{equation}
and the nondimensionalized equations
are naturally formulated
as a fast--slow system,
\begin{subequations}\label{full-kinetics}
\begin{gather}
\dot{s} =-s+(s+\kappa -\lambda)c, \label{full-kinetics-1}\\
\varepsilon \dot{c} =s-(s+\kappa)c \label{full-kinetics-2},
\end{gather}
\end{subequations}
where
\begin{equation}
\label{eq-13d}
t=k_1\bar{E}\tilde{t},\quad s=\frac S{\bar{S}},\quad e=\frac E{\bar{E}},
\quad
c=\frac C{\bar{E}}.
\end{equation}
The dimensionless parameters are
\begin{equation}
\label{eq-13e}
\kappa=\frac{k_{-1}+k_2}{k_1\bar{S}},
\quad
\lambda=\frac{k_2}{k_1\bar{S}}.
\end{equation}

During a short $\mathcal{O}(\varepsilon)$
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 \eqref{full-kinetics}
has an asymptotically stable,
invariant,
slow manifold.
During the transient,
the concentrations
relax to the slow manifold,
decaying exponentially
toward it.
Subsequently,
on the $\mathcal{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
\begin{align}\label{critical-man}
c=\frac{s}{s+\kappa},
\end{align}
and the reduced system dynamics on it
are governed
by a single equation,
\begin{align*}
\frac{ds}{dt}=-\frac{\lambda s}{s+\kappa}.
\end{align*}
This leading order slow manifold
is obtained directly
from \eqref{full-kinetics-2}
with $\varepsilon=0$
and is referred to
as the quasi-steady state approximation
\cite{BAO1963,KS1998,LS1974,M1989,PL1984,S1998}
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~\cite{J1994,K1999}.
We emphasize
that the critical manifold is
only approximately invariant
under the dynamics
of \eqref{full-kinetics};
the exact slow manifold
is invariant.
Additional studies
of the accuracy
of the quasi-steady state approximation
are given in
\cite{BAO1963,HTA1967,LS1974,OM1991,P1987,S1998,SS1989}.

\begin{remark} \label{rmk1.1} \rm
For all sufficiently small, positive
$\varepsilon$,
there is a family
of slow manifolds,
all of which
are exponentially
($\mathcal{O}({\rm e}^{-{\tilde k}/\varepsilon}$)
for some $\tilde{k}>0$)
close to each other,
{\it i.e.}, the asymptotic expansions
of these slow manifolds
are the same to all powers of $\varepsilon$,
see~\cite{J1994,K1999},
for example.
For convenience,
we will sometimes
refer to `the'
(rather than to `a')
slow manifold.
\end{remark}

\subsection{MMH Kinetics with Diffusion}
\label{section1.2}

Given the effectiveness
of the two reduction steps
in the kinetics problem \eqref{1b},
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 \eqref{1a} is irreversible.
The governing equations
in one space dimension are
\begin{subequations}\label{1c}
\begin{align}
\frac{\partial S}{\partial\tilde{t}} &=-k_1SE+k_{-1}
C+D_{S}\frac{\partial^2S}{\partial\tilde{x}^2}, \\
\frac{\partial E}{\partial\tilde{t}} &=-k_1SE+
(k_{-1}+k_2)C+D_{E}\frac{\partial^2E}{\partial\tilde{x}^2}, \\
\frac{\partial C}{\partial\tilde{t}} &=k_1SE-
(k_{-1}+k_2)C+D_{C}\frac{\partial^2C}{\partial\tilde{x}^2},
\end{align}
\end{subequations}
with $\tilde x\in [0,\ell]$,
subject to no-flux (Neumann) boundary conditions
\begin{equation}\label{1d}
\frac{\partial S}{\partial\tilde{x}}\Big|_{\tilde{x}=0,\ell}=
\frac{\partial E}{\partial\tilde{x}}\Big|_{\tilde{x}=0,\ell}=
\frac{\partial C}{\partial\tilde{x}}\Big|_{\tilde{x}=0,\ell}=0
\end{equation}
and initial conditions
\begin{equation}\label{1e}
S(0,\tilde{x})=S_i(\tilde{x}),\quad E(0,\tilde{x})=
E_i(\tilde{x}),\quad C(0,\tilde{x})=0.
\end{equation}
Here,
$\ell>0$ is the $\mathcal{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
\eqref{1c}--\eqref{1e}
as follows.
The nondimensional spatial variable is
$x=\tilde{x}/\ell$.
Time and the species' concentrations
are nondimensionalized
as in \eqref{eq-13d},
but now $\bar{S}$ and $\bar{E}$
denote the spatial averages,
\begin{align*}
\bar{S}=\frac1\ell\int_0^\ell S_i(\tilde{x})\,d\tilde{x}, \quad
\quad\bar{E}=\frac1\ell\int_0^\ell E_i(\tilde{x})\,d\tilde{x}.
\end{align*}
The nondimensional parameters
are again given by \eqref{eq-13e},
and the diffusivities
are scaled as
\begin{align}\label{scalings}
\delta=\frac{D_{S}}{k_1\ell^2\bar{E}}, \quad a=\frac{D_{E}}{D_{S}},
\quad b=\frac{D_{C}}{D_{S}}.
\end{align}
Thus, we obtain the equations
\begin{subequations}\label{1f}
\begin{align}
\frac{\partial s}{\partial t} &=-se+(\kappa-\lambda)c+\delta
\frac{\partial^2s}{\partial x^2}, \\
\frac{\partial e}{\partial t} &=-\frac1\varepsilon(se-\kappa c)+a\delta
\frac{\partial^2e}{\partial x^2}, \\
\frac{\partial c}{\partial t} &=\frac1\varepsilon(se-\kappa c)+b\delta
\frac{\partial^2c}{\partial x^2},
\end{align}
\end{subequations}
on the unit interval,
subject to the Neumann boundary conditions
\begin{align}\label{1g}
\frac{\partial s}{\partial x}\Big|_{x=0,1}=\frac{\partial e}{\partial x}
\Big|_{x=0,1}=\frac{\partial c}{\partial x}\Big|_{x=0,1}=0,
\end{align}
and the initial conditions
\begin{align}\label{1i}
s(0,x)=s_i(x),\quad e(0,x)=e_i(x),\quad c(0,x)=0.
\end{align}
Here, $s_i=S_i/\bar{S}$ and $e_i=E_i/\bar{E}$.
We assume $0<\varepsilon\ll 1$ and that
$a$ and $b$ are $\mathcal{O}(1)$.

In vector notation,
equations \eqref{1f} may be written as
\begin{align}\label{1h}
\frac{\partial\mathbf{u}}{\partial t}=\frac 1\varepsilon\mathbf{F}_\varepsilon(\mathbf{u})+\delta D
\frac{\partial^2\mathbf{u}}{\partial x^2},
\end{align}
where
\begin{align}
\label{1FFu}
\mathbf{F}_\varepsilon(\mathbf{u})=\begin{pmatrix}
-\varepsilon se+\varepsilon(\kappa-\lambda)c \\ -(se-\kappa c) \\
se-\kappa c
\end{pmatrix},
\quad
\mathbf{u}=\begin{pmatrix} s \\ e \\ c \end{pmatrix},
\end{align}
and $D=\mathop{\rm diag}(1,a,b)$.
Moreover, we use $[\mathbf{F}_\varepsilon(\mathbf{u})]_k$
to denote the $k$th order terms in the Taylor
expansion of $\mathbf{F}_\varepsilon$ with respect to $\varepsilon$.
Given a formal asymptotic expansion
${\bf u}(\cdot,x,\varepsilon)=\sum_{k=0}^\infty {\bf u}_k (\cdot,x) \varepsilon^k$
of the solution of \eqref{1h},
$[\mathbf{F}_\varepsilon(\mathbf{u})]_k$ will generically be a function of
$\mathbf{u}_0,\dots,\mathbf{u}_k$.


\subsection{Summary of Main Results}
\label{section1.3}

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:

\begin{enumerate}
\item[(i)] Large diffusivities,
$\delta =\mathcal{O}(1/\varepsilon^2)$:
the diffusive time scale is shorter than
both the fast and the slow kinetic
time scales;

\item[(ii)] Moderately large diffusivities,
$\delta = \mathcal{O}(1/\varepsilon)$:
the diffusive time scale is comparable to the time scale of the
fast kinetics; and

\item[(iii)] Small diffusivities,
$\delta =\mathcal{O}(\varepsilon)$:
the diffusive
time scale is longer than both kinetic time scales.


\end{enumerate}
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 ($\tau=t/\varepsilon^2$) time scale.
Then, the dynamics
on the fast ($\eta=t/\varepsilon$)
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 ($\eta=t/\varepsilon$)
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 ($\eta=t/\varepsilon$) 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
($\zeta=\varepsilon t$)
time scales,
the enzyme profile homogenizes.

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

\begin{table}[ht]
\begin{center}
\begin{tabular}{|l|l|l|l|}
\hline
Regime & $\delta$ & Dynamics & Time scale \\
\hline
\hline
(i) & ${1}/{\varepsilon^2}$ & homogenization of $s, e, c$ & super-fast $\tau={t/\varepsilon^2}$ \\
    &                & fast kinetics               & fast $\eta={t/\varepsilon}$ \\
    &                & slow reduced kinetics       & slow $t$ \\
\hline
(ii) & ${1}/{\varepsilon}$ & homogenization \& fast kinetics & fast $\eta={t}/{\varepsilon}$ \\
     &              & slow reduced kinetics       & slow $t$ \\
\hline
(iii) & $\varepsilon$ & pointwise fast kinetics          & fast $\eta={t}/{\varepsilon}$ \\
      &        & pointwise slow reduced kinetics  & slow $t$ \\
      &        & homogenization of $e$            & super-slow $\zeta=\varepsilon t$ \\
\hline
\end{tabular}
\end{center}
\caption{Summary of the observed dynamics and the time scales of \eqref{1f}
in regimes (i)--(iii)}
\end{table}


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~\cite{VBK1995}.)
Moreover,
we present
numerical simulations
in every regime
to further illustrate
our analysis.

\begin{remark} \label{rmk1.2} \rm
The regime of moderately small diffusivities,
$\delta = \mathcal{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 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.
\end{remark}

\begin{remark} \label{rmk1.3} \rm
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.
\end{remark}

\begin{remark} \label{rmk1.4} \rm
The analysis may also be extended to problems
in which the domain length
$\ell$ is not $\mathcal{O}(1)$.
For example,
the analysis of regime (i)
also applies
to problems
in which $\ell$ is small
and the diffusivities
are not large.
In that case,
it is natural to scale
the spatial variable
as $x=\varepsilon\tilde{x}$.
With this scaling,
the diffusion terms
in \eqref{1f}
are of the form
\[
\delta\frac{\partial^2 s}{\partial x^2}
=\frac{\delta}{\varepsilon^2}\frac{\partial^2 s}{\partial\tilde{x}^2}.
\]
Hence,
diffusion dominates again,
even if the actual diffusion coefficients
are $\delta =\mathcal{O}(1)$.
\end{remark}


\begin{remark} \label{rmk1.5} \rm
The influence of diffusion in the MMH mechanism
has also been studied in \cite{YTBMP1995}.
Specifically, the reduced kinetics model
\eqref{full-kinetics}
is considered and a diffusion term
is added for the substrate only,
with $\mathcal{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 concentration-dependent
correction at $\mathcal{O}(\varepsilon)$.
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.
\end{remark}


\begin{remark} \label{rmk1.6} \rm
It has been shown in~\cite{SS1989}
that the effective small parameter
in the MMH mechanism is
$\tilde{\varepsilon}=\bar{E}/(\bar{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.
\end{remark}


This article is organized as follows. The regimes (i)--(iii)
are analyzed in
Sections~\ref{section2}--\ref{section4}, respectively.
In Section~\ref{section5}, 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 \eqref{1f} is linearly stable,
irrespective of the magnitudes of the diffusion coefficients.
Appendix~B contains a technical result
relating to Section~\ref{section3}.


\section{Large Diffusivities}\label{section2}

In this section, we consider the regime in which the diffusivities
of all species are large,
$\delta=\mathcal{O}(1/\varepsilon^2)$; for convenience, we choose
$\delta=1/\varepsilon^2$ in \eqref{1f}.
Here, the time scale
on which diffusion acts is much shorter than that of the fast kinetics.
There is a very short transient period, $\mathcal{O}(\varepsilon^2)$ in duration,
in which the initially heterogeneous species' concentrations,
given by \eqref{1i}, homogenize
and during which essentially no reaction takes place.
After this short transient,
the problem reduces to the well-understood,
classical problem
of pure kinetics for the homogeneous solution,
see, {\it e.g.},~\cite{LS1974}.
We treat this regime in some detail to introduce
the method employed in this article
in an elementary context.

After introduction of $\delta =1/\varepsilon^2$, equations \eqref{1f} become
\begin{subequations}\label{2a}
\begin{gather}
\varepsilon^2\frac{\partial s}{\partial t} =-\varepsilon^2se+\varepsilon^2(\kappa-\lambda)c
+\frac{\partial^2s}{\partial x^2}, \\
\varepsilon^2\frac{\partial e}{\partial t} =-\varepsilon(se-\kappa c)+a
\frac{\partial^2e}{\partial x^2}, \\
\varepsilon^2\frac{\partial c}{\partial t} =\varepsilon(se-\kappa c)+b
\frac{\partial^2c}{\partial x^2}.
\end{gather}
\end{subequations}
Equivalently, in vector form,
\begin{equation}
\label{2b}
\varepsilon^2\frac{\partial\mathbf{u}}{\partial t}=\varepsilon\mathbf{F}_\varepsilon(\mathbf{u})+D\frac{\partial^2\mathbf{u}}
{\partial x^2},
\end{equation}
where $\mathbf{F}$ is defined in \eqref{1FFu}.

\subsection{Homogenization: The Super-Fast (Inner) Time Scale}
\label{section2.1}

To study the initial transient period,
we let $\tau=t/\varepsilon^2$
denote the super-fast (inner) time
and let $\hat{\mathbf{u}}(\tau,x,\varepsilon)=\mathbf{u}(t,x,\varepsilon)$.
The governing equations become
\begin{equation}\label{2c}
\frac{\partial\hat{\mathbf{u}}}{\partial\tau}=\varepsilon\mathbf{F}_\varepsilon(\hat{\mathbf{u}})
 +D\frac{\partial^2\hat{\mathbf{u}}}{\partial x^2},
\end{equation}
with initial and boundary conditions
\begin{equation}\label{2g}
\hat{\mathbf{u}}(0,x,\varepsilon)=\mathbf{u}_i(x)
\quad {\rm and}\quad
\frac{\partial\hat{\mathbf{u}}}{\partial x}\Big|_{x=0,1}=0.
\end{equation}

We consider \eqref{2c} and \eqref{2g}
over an $\mathcal{O}(1)$--interval
of $\tau$ time,
starting at $\tau=0$.
Asymptotically,
as $\varepsilon \to 0^+$,
the solution
can be expressed
using the Ansatz
\[
\hat{\mathbf{u}}(\tau,x,\varepsilon)=\hat{\mathbf{u}}_0(\tau,x)+\varepsilon\hat{\mathbf{u}}_1(\tau,x)+\mathcal{O}(\varepsilon^2).
\]
Hence, we expand both sides of \eqref{2c} in powers of $\varepsilon$ to obtain
a recursive sequence of differential equations for $\hat{\mathbf{u}}_k$, $k=1,2,
\dots$. At $\mathcal{O}(1)$, $\hat{\mathbf{u}}_0$ satisfies
the homogeneous heat equation,
\begin{equation} \label{2d}
\mathcal{O}(1):\quad\mathcal{L}_\tau\hat{\mathbf{u}}_0=0,
\end{equation}
where
$\mathcal{L}_\tau ={\partial}/{\partial\tau} -D{\partial^2}/{\partial x^2}$,
subject to Neumann boundary conditions.
The solution is
\begin{align}\label{2e}
\hat{\mathbf{u}}_0(\tau,x)=\sum_{k=0}^\infty\hat{\mathbf{u}}_{0k}{\rm e}^{-D(k\pi)^2\tau}
\cos{(k\pi x)}.
\end{align}
Here, the coefficients $\hat{\mathbf{u}}_{0k}$ are
the Fourier coefficients of the initial distribution
$\mathbf{u}_i(x)$ with respect to $\{\cos{(k\pi x)}\}_{k\ge 0}$,
and they are constant
during the fast transients on the $\tau$ time scale.


Asymptotically, as $\tau \to \infty$,
we find
\begin{equation} \label{2f}
\hat{\mathbf{u}}_0(\tau,x) \to \hat{\mathbf{u}}_{00}
=\int_0^1\mathbf{u}_i(x)\,dx
=\begin{pmatrix} 1 \\ 1 \\ 0 \end{pmatrix},
\end{equation}
where we used \eqref{1i}.
Therefore, asymptotically,
the effect of diffusion
in (\ref{2a})
is to smear out
the initial distributions
of the reactants
until they are effectively
uniformly distributed
over the entire spatial domain.

\begin{remark} \label{rmk2.1} \rm
It will suffice to consider
the leading-order fast solution (\ref{2e})
to accomplish matching
to lowest order
in the next subsection.
To that end,
it is useful to write (\ref{2e}) as
\begin{align}\label{2h}
\hat{\mathbf{u}}_0(\tau,x)=(1,1,0)^T+\mathcal{O}({\rm e}^{-d\pi^2\tau}),
\end{align}
where $d=\min\{1,a,b\}$
and where we used \eqref{2f}.
\end{remark}


\subsection{Fast Kinetics: The Fast Time Scale}
\label{section2.2}

The fast kinetics take place on the fast
$\eta=t/\varepsilon$
time scale, during which
the system dynamics
are given by
\begin{equation} \label{3a}
\varepsilon\frac{\partial\tilde{\mathbf{u}}}{\partial\eta}=\varepsilon\mathbf{F}_\varepsilon(\tilde{\mathbf{u}})+
D\frac{\partial^2\tilde{\mathbf{u}}}{\partial x^2},
\end{equation}
subject to Neumann boundary conditions.
Here, $\tilde{\mathbf{u}}=\tilde{\mathbf{u}}(\eta,x,\varepsilon)$,
and we assume
$\tilde{\mathbf{u}}(\eta,x,\varepsilon)=\tilde{\mathbf{u}}_0(\eta,x)+\varepsilon\tilde{\mathbf{u}}_1(\eta,x)+
\mathcal{O}(\varepsilon^2)$.
Then, expanding \eqref{3a} in powers of $\varepsilon$ and rearranging the resulting
equations, we find
\begin{subequations}\label{3c}
\begin{align}
\mathcal{O}(1):\quad &-D\frac{\partial^2\tilde{\mathbf{u}}_0}{\partial x^2}=0, \label{3c1} \\
\mathcal{O}(\varepsilon):\quad &-D\frac{\partial^2\tilde{\mathbf{u}}_1}{\partial x^2}=-
\frac{\partial\tilde{\mathbf{u}}_0}{\partial\eta}+[\mathbf{F}_\varepsilon(\tilde{\mathbf{u}})]_0, \label{3c2}
\end{align}
\end{subequations}
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 \eqref{2a}.

Integrating \eqref{3c1} and taking into account the boundary conditions, we
conclude that
\begin{align}\label{31a}
\tilde{\mathbf{u}}_0(\eta,x)=\tilde{\mathbf{u}}_0(\eta),
\end{align}
{\it i.e.}, that $\tilde{\mathbf{u}}_0$ is independent of $x$. Similarly, it
follows from \eqref{3c2} that
\begin{equation} \label{3e}
-D\int_0^1\frac{\partial^2\tilde{\mathbf{u}}_1}{\partial x^2}\,dx=
-D\frac{\partial\tilde{\mathbf{u}}_1}{\partial x}\Big|_{x=0}^1=0=\int_0^1
\Big(-\frac{d\tilde{\mathbf{u}}_0}{d\eta}+[\mathbf{F}_\varepsilon(\tilde{\mathbf{u}})]_0\Big)\,dx.
\end{equation}
Since the integrand in the right member of \eqref{3e} is independent
of $x$, we see that the dynamics of $\tilde{\mathbf{u}}_0$ are governed by the ordinary
differential equation
\[
\frac{d\tilde{\mathbf{u}}_0}{d\eta}=[\mathbf{F}_\varepsilon(\tilde{\mathbf{u}})]_0,
\]
which we write out componentwise,
\begin{subequations}\label{3f}
\begin{align}
\frac{d\tilde{s}_0}{d\eta} &=0, \\
\frac{d\tilde{e}_0}{d\eta} &=-(\tilde{s}_0\tilde{e}_0-\kappa\tilde{c}_0), \\
\frac{d\tilde{c}_0}{d\eta} &=\tilde{s}_0\tilde{e}_0-\kappa\tilde{c}_0.
\end{align}
\end{subequations}
These equations are the same as one finds to leading order
in the classical MMH kinetics problem; hence,
they can be solved explicitly:
$\tilde{s}_0(\eta)\equiv\tilde{s}_0(0)$ and
\begin{equation}
\tilde{e}_0(\eta)+\tilde{c}_0(\eta)=\tilde{e}_0(0)+\tilde{c}_0(0)\quad
\text{for all }\eta\ge 0.
\nonumber
\end{equation}
The initial conditions for \eqref{3f} are determined by matching with
the leading-order equations on the super-fast ($\tau$) scale;
notably, we require
that $\lim_{\tau\to\infty}\hat{\mathbf{u}}_0(\tau,x)=\lim_{\eta\to 0^+}
\tilde{\mathbf{u}}_0(\eta,x)$. Now, by \eqref{2h},
\begin{equation} \label{3h}
\lim_{\tau\to\infty}\hat{\mathbf{u}}_0(\tau,x)=(1,1,0)^T.
\end{equation}
Hence, we have
\begin{equation} \label{31b}
\tilde{s}_0(\eta)\equiv 1\quad\text{and}\quad\tilde{e}_0(\eta)=
1-\tilde{c}_0(\eta).
\end{equation}
In turn, it follows that
\[
\frac{d\tilde{c}_0}{d\eta}=1-(1+\kappa)\tilde{c}_0,\quad\text{with }
\tilde{c}_0(0)=0,
\]
and, hence,
\begin{equation} \label{3g}
\tilde{\mathbf{u}}_0(\eta,x)=\Big(
1,\frac 1{1+\kappa}\big(\kappa +{\rm e}^{-(1+\kappa)\eta}\big),
\dfrac 1{1+\kappa}\big(1-{\rm e}^{-(1+\kappa)\eta}\big)\Big)^T.
\end{equation}
Therefore, on the fast ($\eta$) scale,
the species' concentrations are essentially homogeneous,
and the fast chemistry occurs,
with the binding of enzyme and substrate to form complex.

\subsection{Slow Reduced Kinetics: The Slow (Outer) Time Scale}
\label{section2.3}

The slow, reduced kinetics take place
on the slow ($t$) time scale.
The dynamics are governed
by the original system,
\eqref{2b},
subject to Neumann boundary conditions.
We expand
$\mathbf{u}(t,x,\varepsilon)=\mathbf{u}_0(t,x)+\varepsilon\mathbf{u}_1(t,x)+\varepsilon^2\mathbf{u}_2(t,x)+\mathcal{O}(\varepsilon^3)$
and equate coefficients of equal powers of $\varepsilon$ to obtain
\begin{subequations}\label{4b}
\begin{align}
\mathcal{O}(1):\quad &-D\frac{\partial^2\mathbf{u}_0}{\partial x^2}=0, \label{4b1} \\
\mathcal{O}(\varepsilon):\quad &-D\frac{\partial^2\mathbf{u}_1}{\partial x^2}=[\mathbf{F}_\varepsilon(\mathbf{u})]_0,
\label{4b2} \\
\mathcal{O}(\varepsilon^2):\quad &-D\frac{\partial^2\mathbf{u}_2}{\partial x^2}=
-\frac{\partial\mathbf{u}_0}{\partial t}+[\mathbf{F}_\varepsilon(\mathbf{u})]_1. \label{4b3}
\end{align}
\end{subequations}
Applying the same type of solvability argument used in Section 2.2,
we conclude from \eqref{4b1} that $\mathbf{u}_0(t,x)=\mathbf{u}_0(t)$. Similarly,
the solvability for \eqref{4b2} implies $\mathbf{u}_1(t,x)=\mathbf{u}_1(t)$,
since the right member
of \eqref{4b2}
is independent of $x$
and, hence, $[\mathbf{F}_\varepsilon(\mathbf{u})]_0=0$.
In turn, this yields
\begin{equation} \label{4c}
s_0(t)e_0(t)-\kappa c_0(t)=0.
\end{equation}

Next, by writing out \eqref{4b3} componentwise
and by applying a solvability
argument similar to the one used in \eqref{3e},
we find
\begin{subequations}\label{41a}
\begin{align}
\frac{ds_0}{dt} &=-s_0e_0+(\kappa-\lambda)c_0, \label{41a1} \\
\frac{de_0}{dt} &=-(s_0e_1+s_1e_0-\kappa c_1), \label{41a2} \\
\frac{dc_0}{dt} &=s_0e_1+s_1e_0-\kappa c_1. \label{41a3}
\end{align}
\end{subequations}
In turn, equation \eqref{41a1} may be simplified using \eqref{4c} to obtain
$ds_0/dt=-\lambda c_0. $
In addition, \eqref{41a2} and \eqref{41a3} imply
\begin{equation} \label{41b}
e_0(t)+c_0(t)=e_0(0)+c_0(0)\quad\text{for all }t\ge 0,
\end{equation}
where the constant is to be determined by matching
with the equations
on the fast ($\eta$) scale:
$\lim_{\eta\to\infty}\tilde{\mathbf{u}}_0(\eta)
=\lim_{t\to 0^+}\mathbf{u}_0(t)$.
From (\ref{3g}), we find
\begin{equation} \label{4h}
\lim_{\eta\to\infty}\tilde{\mathbf{u}}_0(\eta)=\Big(1,\frac\kappa{1+\kappa},
\frac1{1+\kappa}\Big)^T.
\end{equation}
Hence,
\begin{equation} \label{4e0c0}
e_0(t)+c_0(t)=1.
\end{equation}

Finally,
we combine \eqref{4c} and \eqref{4e0c0} to obtain
the critical manifold from the classical kinetics problem
with $\varepsilon=0$,
\begin{equation} \label{4f}
c_0(t)=\frac{s_0(t)}{s_0(t)+\kappa}
\quad\text{and}\quad
e_0(t)=\frac{\kappa}{s_0(t)+\kappa}.
\end{equation}
Moreover,
we see that
the reduced equation for $s_0(t)$ on this critical manifold is
\begin{equation} \label{4d}
\frac{ds_0}{dt}=-\lambda\frac{s_0}{s_0+\kappa}\quad\text{with } s_0(0)=1,
\end{equation}
just as is the case for the pure MMH kinetics problem,
see for example~\cite{LS1974}.
The solution of \eqref{4d}
is known implicitly,
\begin{equation} \label{4g}
s_0(t)+\kappa\ln{s_0(t)}=-\lambda t+1.
\end{equation}
Also, the rate of approach toward
the slow manifold
is determined by the dynamics
on the fast
($\eta$) scale, cf. \eqref{3g}.


\subsection{The Uniformly Valid Leading-Order Approximation}
\label{section2.4}

In this regime of large diffusivities, the
leading-order approximation, uniformly valid in time
and space, to the
solution of \eqref{1f} is obtained
by combining the expressions for $\hat{\mathbf{u}}_0$, $\tilde{\mathbf{u}}_0$,
and $\mathbf{u}_0$ (recall \eqref{2e}, \eqref{3g}, \eqref{4f}, and \eqref{4g})
and subtracting their respective common parts (recall \eqref{3h} and
\eqref{4h}). We find
\begin{subequations}\label{13a}
\begin{align}
s(t,x) &=s_0(t)+\mathcal{O}\big({\rm e}^{-\pi^2\delta t}\big)+\mathcal{O}(\varepsilon), \label{13a1} \\
e(t,x) &= \frac{\kappa}{s_0(t)+\kappa}
         +\frac{{\rm e}^{-(1+\kappa)\sqrt{\delta}t}}{1+\kappa}
         +\mathcal{O}\big({\rm e}^{-b\pi^2\delta t}\big)
         +\mathcal{O}(\varepsilon), \label{13a2} \\
c(t,x) &= \frac{s_0(t)}{s_0(t)+\kappa}
         -\frac{{\rm e}^{-(1+\kappa)\sqrt{\delta}t}}{1+\kappa}
         +\mathcal{O}\big({\rm e}^{-b\pi^2\delta t}\big)
         +\mathcal{O}(\varepsilon), \label{13a3}
\end{align}
\end{subequations}
where $s_0$ is defined by \eqref{4g},
and we recall that $\delta=1/\varepsilon^2$,
{\it i.e.,} $\delta$
is asymptotically larger
than the rate constant
corresponding to the faster
of the two kinetics scales.

The physical interpretation
of \eqref{13a}
is that,
during a short
initial time interval
of $\mathcal{O}(\varepsilon^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 $\delta\to\infty$, the expressions in \eqref{13a} agree
with the results for the chemical kinetics problem considered
for example in Lin and Segel~\cite[Equations~(14) and (15)]{LS1974}.
Moreover, the algebraic corrections at $\mathcal{O}(\varepsilon)$ 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 $\mathbf{u}_1$, for example,
in the following subsection.


\subsection{Higher-Order Corrections}
\label{section2.5}

On the slow time scale $t$,
the $\mathcal{O}(\varepsilon)$ corrections to the leading-order solution
are characterized
by the $\mathcal{O}(\varepsilon^3)$--terms in \eqref{2b},
\begin{equation} \label{4e}
\mathcal{O}(\varepsilon^3):\quad -D\frac{\partial^2\mathbf{u}_3}{\partial x^2}=-\frac{\partial\mathbf{u}_1}
{\partial t}+[\mathbf{F}_\varepsilon(\mathbf{u})]_2.
\end{equation}
Equation \eqref{4b3} yields that $\mathbf{u}_2(t,x)=\mathbf{u}_2(t)$;
hence, application
of the solvability condition to \eqref{4e} implies
\begin{align*}
\frac{ds_1}{dt} &=-s_1e_0-s_0e_1+(\kappa -\lambda)c_1, \\
\frac{de_1}{dt} &=-(s_0e_2+s_1e_1+s_2e_0-\kappa c_2), \\
\frac{dc_1}{dt} &= s_0e_2+s_1e_1+s_2e_0-\kappa c_2.
\end{align*}
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+(\kappa -\lambda+s_0)c_1,
\]
which is precisely~\cite[Equation~(25a)]{LS1974}.
One may proceed in
a similar manner to obtain the asymptotics of $\mathbf{u}$ to any order,
as well as the $\mathcal{O}(\varepsilon)$ 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 $\mathcal{O}(\varepsilon)$,
system \eqref{2c} yields
\begin{align*}
\frac{\partial \hat{s}_1}{\partial \tau}
 &= \frac{\partial^2 \hat{s}_1}{\partial x^2}, \\
\frac{\partial \hat{e}_1}{\partial \tau}
 &= -(\hat{s}_0 \hat{e}_0 - \kappa \hat{c}_0)
             + a \frac{\partial^2 \hat{e}_1}{\partial x^2}, \\
\frac{\partial \hat{c}_1}{\partial \tau}
 &= \hat{s}_0 \hat{e}_0 - \kappa \hat{c}_0
             + b \frac{\partial^2 \hat{c}_1}{\partial x^2}.
\end{align*}
Hence, substituting the leading-order
solution \eqref{2e}, one sees that
the constant terms, corresponding to $k=0$,
lead to linearly growing and decaying terms
in $\hat{e}_1$
and $\hat{c}_1$.
These secular terms render the expansion
invalid as an approximation on time scales of $\tau = \mathcal{O}(1/\varepsilon)$.
Therefore, one needs to use the multiple scales method
(or, alternatively, the boundary function method), as follows.
We let
\[
\hat{\mathbf{u}}_0=\hat{\mathbf{u}}_0(\eta, \tau, x)
          =\hat{\mathbf{u}}_{00}(\eta)
            + \Sigma_{k=1}^\infty \hat{\mathbf{u}}_{0k}
               {\rm e}^{-D(k\pi)^2 \tau}\cos(k\pi x),
\]
so that ${\hat {\bf u}}_{00}$
varies on the fast scale,
while
for $k\ge1$
$\hat{\mathbf{u}}_{0k}$
remains constant, as before.
Therefore, the equations at $\mathcal{O}(\varepsilon)$ are now
\begin{align*}
\frac{\partial \hat{s}_1}{\partial \tau}
+\frac{\partial \hat{s}_{00}}{\partial \eta}
 &= \frac{\partial^2 \hat{s}_1}{\partial x^2}, \\
\frac{\partial \hat{e}_1}{\partial \tau}
+\frac{\partial \hat{e}_{00}}{\partial \eta}
 &= -(\hat{s}_0 \hat{e}_0 - \kappa \hat{c}_0)
             + a \frac{\partial^2 \hat{e}_1}{\partial x^2}, \\
\frac{\partial \hat{c}_1}{\partial \tau}
+\frac{\partial \hat{c}_{00}}{\partial \eta}
 &= \hat{s}_0 \hat{e}_0 - \kappa \hat{c}_0
             + b \frac{\partial^2 \hat{c}_1}{\partial x^2}.
\end{align*}
Solvability (or `elimination of the secular terms')
implies that one should choose the $\eta$--dependence
in $\hat{\mathbf{u}}_{00}$ so that
\begin{align*}
\frac{\partial \hat{s}_{00}}{\partial \eta} &= 0, \\
\frac{\partial \hat{e}_{00}}{\partial \eta}
 &= - (\hat{s}_{00} \hat{e}_{00} - \kappa \hat{c}_{00}), \\
\frac{\partial \hat{c}_{00}}{\partial \eta}
 &= \hat{s}_{00} \hat{e}_{00} - \kappa \hat{c}_{00}.
\end{align*}
By this choice, the solution of the system
for $\hat{\mathbf{u}}_1$ is bounded.
Also, we observe that these equations
are exactly the same as equations \eqref{3f}
for $\tilde{\mathbf{u}}_0$,
as expected.
For an exposition of the method of multiple scales,
see~\cite{KC1996}, for example.


\section{Moderately Large Diffusivities}
\label{section3}


In this section, we examine the regime
of moderately large diffusivities,
$\delta=\mathcal{O}(1/\varepsilon)$;
for convenience,
we take $\delta=1/\varepsilon$
in \eqref{1f}.
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~\ref{section2.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
\begin{subequations}\label{42a}
\begin{align}
\varepsilon\frac{\partial s}{\partial t} &=-\varepsilon se+\varepsilon(\kappa-\lambda)c+
\frac{\partial^2s}{\partial x^2}, \label{42a1} \\
\varepsilon\frac{\partial e}{\partial t} &=-(se-\kappa c)+a\frac{\partial^2e}
{\partial x^2}, \label{42a2} \\
\varepsilon\frac{\partial c}{\partial t} &=se-\kappa c+b\frac{\partial^2c}
{\partial x^2}. \label{42a3}
\end{align}
\end{subequations}
Equivalently, in vector form, they are given by
\begin{align*}
\varepsilon\frac{\partial\mathbf{u}}{\partial t}=\mathbf{F}_\varepsilon(\mathbf{u})+D\frac{\partial^2\mathbf{u}}
{\partial x^2}.
\end{align*}

\subsection{Homogenization and Fast Kinetics: The Fast (Inner) Time Scale}
\label{section3.1}


Recall that $\hat{\mathbf{u}}(\eta,x,\varepsilon)
=\mathbf{u}(t,x,\varepsilon)$ on the fast time scale given
by $\eta=t/\varepsilon$.
Then, the governing equations
become
\begin{subequations}\label{42c}
\begin{align}
\frac{\partial\hat s}{\partial\eta} &=-\varepsilon\hat s\hat e+\varepsilon(\kappa-\lambda)
\hat c+\frac{\partial^2\hat s}{\partial x^2}, \\
\frac{\partial\hat e}{\partial\eta} &=-(\hat s\hat e-\kappa\hat c)
+a\frac{\partial^2\hat e}{\partial x^2}, \\
\frac{\partial\hat c}{\partial\eta} &=\hat s\hat e-\kappa\hat c
+b\frac{\partial^2\hat c}{\partial x^2} \label{42c3}
\end{align}
\end{subequations}
subject to the usual initial conditions
and Neumann boundary conditions.

Making the Ansatz
$\hat{\mathbf{u}}(\eta,x,\varepsilon)=\hat{\mathbf{u}}_0(\eta,x)+\varepsilon\hat{\mathbf{u}}_1(\eta,x)+\mathcal{O}(\varepsilon^2), $
we find to lowest order that
\begin{subequations}\label{42d}
\begin{align}
\frac{\partial\hat s_0}{\partial\eta} &=\frac{\partial^2\hat s_0}
{\partial x^2}, \label{42d1} \\
\frac{\partial\hat e_0}{\partial\eta} &=-(\hat s_0\hat e_0-\kappa\hat c_0)
+a\frac{\partial^2\hat e_0}{\partial x^2}, \label{42d2} \\
\frac{\partial\hat c_0}{\partial\eta} &=\hat s_0\hat e_0-\kappa\hat c_0
+b\frac{\partial^2\hat c_0}{\partial x^2}. \label{42d3}
\end{align}
\end{subequations}
Hence,
$\hat s_0$ satisfies
the homogeneous heat equation
with Neumann boundary conditions
and with initial data
given by
$\hat{s}_0(0,x)=s_i(x)$,
which implies that
\begin{equation} \label{55a}
\hat s_0(\eta,x)=\sum_{k=0}^\infty\hat s_{0k}{\rm e}^{-(k\pi)^2\eta}\cos{(k\pi x)},
\end{equation}
where $\hat s_{00}=1$ and
\[
\hat s_{0k}=2\int_0^1\hat s_0(0,x)\cos{(k\pi x)}\,dx=2\int_0^1s_i(x)
\cos{(k\pi x)}\,dx\quad\text{for }k\ge 1.
\]

Next, we find the corresponding expressions
for $\hat e_0$ and $\hat c_0$.
Equations \eqref{42d2} and \eqref{42d3}
are linear in $\hat{e}_0$ and $\hat{c}_0$,
with some of the coefficients depending
on $\eta$ and $x$ through $s_0$.
Hence,
due to the Neumann boundary conditions
on $\hat e_0$ and $\hat c_0$,
one can write
\begin{equation} \label{55b}
\hat e_0(\eta,x)=\sum_{k=0}^\infty\hat e_{0k}(\eta)\cos{(k\pi x)}\quad
\text{and}\quad\hat c_0(\eta,x)=\sum_{k=0}^\infty\hat c_{0k}(\eta)
\cos{(k\pi x)}
\end{equation}
for some sets of coefficients $\{\hat e_{0k}(\eta)\}$ and
$\{\hat c_{0k}(\eta)\}$.
Substituting \eqref{55a} and
\eqref{55b}
into \eqref{42d2} and \eqref{42d3},
making use of the identity
\[
\cos{(m\pi x)}\cos{(n\pi x)}=\frac12\big(\cos{((m+n)\pi x)}+
\cos{((m-n)\pi x)}\big)
\]
as well as of $\hat s_{00}=1$, and collecting coefficients of like cosines,
we find that $\hat e_{0k}$ and $\hat c_{0k}$ must satisfy the following
infinite system of linear ordinary differential equations:
 \begin{subequations}\label{55c}
 \begin{align}
 \frac{d\hat e_{00}}{d\eta}
&=-(\hat e_{00}
 -\kappa\hat c_{00})-\mathcal{F}_0, \\
 \frac{d\hat c_{00}}{d\eta}
 &= \hat e_{00}
 -\kappa\hat c_{00}+\mathcal{F}_0,
 \end{align}
 \end{subequations}
 for the zeroth Fourier mode, and
 \begin{subequations}\label{55c-generalk}
 \begin{align}
 \frac{d\hat e_{0k}}{d\eta} &=
 -\big(1+a(k\pi)^2+\frac{1}{2}\hat s_{0(2k)}{\rm e}^{-4(k\pi)^2t}\big)\hat e_{0k}
 +\kappa\hat c_{0k}-\mathcal{F}_k, \\
 \frac{d\hat c_{0k}}{d\eta} &=
 -\big(\kappa+b(k\pi)^2\big)\hat c_{0k}
 +\big(1+\frac{1}{2}\hat s_{0(2k)}{\rm e}^{-4(k\pi)^2t}\big)\hat e_{0k}+\mathcal{F}_k,
 \end{align}
 \end{subequations}
 for the $k$th Fourier mode, $k\ge1$.
 Here, $\mathcal{F}_k$ is defined via
 \begin{subequations}
 \begin{align}
 \mathcal{F}_0(\eta) &=\frac12\sum_{m\ge1}
 \hat s_{0m}{\rm e}^{-(m\pi)^2\eta}\hat e_{0m}(\eta),\\
 \mathcal{F}_k(\eta) &=\frac12\sum_{\substack{m\ge0 \\ m\ne k}}
 \Big(\hat s_{0(m+k)}{\rm e}^{-(m+k)^2\pi^2\eta}+\hat s_{0\abs{m-k}}{\rm e}^{-(m-k)^2\pi^2\eta}\Big)\hat e_{0m}(\eta),
 \quad k\ge1.
 \end{align}
 \end{subequations}
In particular, \eqref{55c} yields that
$\hat e_{00}(\eta)+\hat c_{00}(\eta)$ is constant
in this regime.
Then, solving \eqref{55c}
and taking into account the identities
\begin{gather*}
\hat e_{00}(0)=\int_0^1\hat e_0(0,x)\,dx=\int_0^1e_i(x)\,dx=1,\\
 \hat c_{00}(0)=\int_0^1\hat c_0(0,x)\,dx=\int_0^1c_i(x)\,dx=0
\end{gather*}
as well as the fact that $\hat e_{0k}$ and $\hat c_{0k}$ must
be bounded in $\eta$, we find
\begin{equation}\label{3.8}
\hat e_{00}(\eta)=\frac\kappa{1+\kappa}+\mathcal{O}({\rm e}^{-\min\{\pi^2,1+\kappa\}\eta}),
\quad\hat c_{00}(\eta)=\frac1{1+\kappa}+
\mathcal{O}({\rm e}^{-\min\{\pi^2,1+\kappa\}\eta}).
\end{equation}
Similar expressions can be derived for $\hat e_{0k}$ and $\hat c_{0k}$ with
$k\ge 1$; in particular, one can show that $\hat e_{0k},\hat c_{0k}
=\mathcal{O}({\rm e}^{-d\pi^2\eta})$, where we recall $d=\min\{1,a,b\}$. Hence, we
conclude that for $\eta$ large,
\begin{equation} \label{42m}
\hat{\mathbf{u}}_0(\eta,x)=\Big(1,\frac\kappa{1+\kappa},\frac1{1+\kappa}\Big)^T
+\mathcal{O}({\rm e}^{-\min\{d\pi^2,1+\kappa\}\eta}).
\end{equation}
See Appendix~B for details.


\subsection{Reduced Kinetics: The Slow (Outer) Time Scale}
\label{section3.2}

On the slow (outer) time scale, the system dynamics are naturally
described by \eqref{42a} in the original time $t$. To
lowest order, we find
\begin{subequations}\label{42g}
\begin{align}
0 &=\frac{\partial^2s_0}{\partial x^2}, \label{42g1} \\
0 &=-(s_0e_0-\kappa c_0)+a\frac{\partial^2e_0}{\partial x^2},
\label{42g2} \\
0 &=s_0e_0-\kappa c_0+b\frac{\partial^2c_0}{\partial x^2}, \label{42g3}
\end{align}
\end{subequations}
where we again make the Ansatz
$\mathbf{u}(t,x,\varepsilon)=\mathbf{u}_0(t,x)+\varepsilon\mathbf{u}_1(t,x)+\mathcal{O}(\varepsilon^2)$.

Integrating \eqref{42g1} with respect to $x$ and making use of
the Neumann boundary conditions, we conclude immediately that
$s_0(t,x)\equiv s_0(t)$. Similarly, it follows from either \eqref{42g2}
or \eqref{42g3} that
\begin{equation} \label{99b}
s_0(t)\int_0^1e_0(t,x)\,dx-\kappa\int_0^1c_0(t,x)\,dx=0.
\end{equation}
To obtain explicit formulae for $e_0$ and $c_0$, we rescale $x$
via $\bar{x}=x/\sqrt{a}$,
and observe that the ratio
\[
\alpha = a/b
\]
is the relevant parameter, rather than $a$ and $b$ separately.
We rewrite \eqref{42g2}
and \eqref{42g3} as a four-dimensional linear system with
$t$--dependent coefficients,
\begin{subequations}\label{99a}
\begin{align}
e_0' &=f_0, \\
f_0' &=s_0(t) e_0-\kappa c_0, \\
c_0' &=d_0, \\
d_0' &=-\alpha(s_0(t) e_0-\kappa c_0),
\end{align}
\end{subequations}
where the prime denotes
(partial) differentiation
with respect to $\bar{x}$.

The general solution
of \eqref{99a}
is given by
\begin{subequations}\label{99c}
\begin{gather}
e_0 (t,\bar{x})=\kappa\gamma_1(t)+\kappa\gamma_2(t)\bar{x}
+\gamma_3(t) {\rm e}^{\sqrt{s_0(t)+\alpha\kappa}\,\bar{x}}
-\gamma_4(t) {\rm e}^{-\sqrt{s_0(t)+\alpha\kappa}\,\bar{x}},
 \label{99c1} \\
c_0 (t,\bar{x})=s_0(t)\gamma_1(t)+s_0(t)\gamma_2(t)\bar{x}
-\alpha\gamma_3(t) {\rm e}^{\sqrt{s_0(t)+\alpha\kappa}\,\bar{x}}
+\alpha\gamma_4(t) {\rm e}^{-\sqrt{s_0(t)+\alpha\kappa}\,\bar{x}},
 \label{99c2}
\end{gather}
\end{subequations}
where $f_0=e_0'$ and $d_0=c_0'$ can be found from \eqref{99c}, and
$\gamma_1,\dots,\gamma_4$ are $t$--dependent constants of integration.
Making use of the Neumann boundary conditions on $e_0$ and $c_0$ in \eqref{99c},
one sees that $\gamma_2$, $\gamma_3$, and $\gamma_4$ must be identically zero.
Hence, the only solution \eqref{99c}
that satisfies the boundary conditions
is given by
\begin{equation} \label{99d}
(e_0,c_0)(t)=\gamma_1(t)(\kappa,s_0(t)),
\end{equation}
where $\gamma_1(t)$ is as yet undetermined, and
we see
that \eqref{99b} reduces to
\begin{equation} \label{99e}
s_0(t)e_0(t)-\kappa c_0(t)=0.
\end{equation}
Summarizing, we see that $e_0$ and $c_0$
are spatially uniform on the slow time scale
and that the constraint \eqref{99e}
is satisfied at all times.

To describe the dynamics of $\mathbf{u}_0$
on the $t$ time scale,
we consider
the $\mathcal{O}(\varepsilon)$--terms in \eqref{42a},
\begin{subequations}\label{42i}
\begin{align}
\frac{ds_0}{dt} &=-s_0e_0+(\kappa-\lambda)c_0+\frac{\partial^2s_1}{
\partial x^2}, \label{42i1} \\
\frac{de_0}{dt} &=-(s_1e_0+s_0e_1-\kappa c_1)+
a\frac{\partial^2e_1}{\partial x^2}, \label{42i2} \\
\frac{dc_0}{dt} &=s_1e_0+s_0e_1-\kappa c_1+
b\frac{\partial^2c_1}{\partial x^2}. \label{42i3}
\end{align}
\end{subequations}
First, an integration of \eqref{42i1} over $x\in[0,1]$ in combination
with \eqref{99e} yields
\begin{equation}\label{55g}
\frac{ds_0}{dt}=-\lambda c_0.
\end{equation}
To find the corresponding dynamics
of $e_0$ and $c_0$,
we add \eqref{42i2}
and \eqref{42i3}
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\ge 0$.
Finally, taking into account that $\lim_{t\to 0^+}\mathbf{u}(t,x)$ must equal
\[
\lim_{\eta\to\infty}\hat{\mathbf{u}}_0(\eta,x)=\Big(1,\frac\kappa{1+\kappa},
\frac1{1+\kappa}\Big)^T,
\]
we obtain
\[
e_0(t)+c_0(t)=1.
\]

Combining this identity with \eqref{99d} and solving the resulting equation
for $\gamma_1$, we find $\gamma_1(t)=(s_0(t)+\kappa)^{-1}$ and therefore
\begin{equation} \label{42l}
e_0(t)=\frac{\kappa}{s_0(t)+\kappa}\quad\text{and}\quad
c_0(t)=\frac{s_0(t)}{s_0(t)+\kappa},
\end{equation}
see \eqref{4f}.
Finally, substitution of \eqref{42l} into \eqref{55g} yields the
governing equation for $s_0$,
\begin{equation}\label{sec4s0}
\frac{ds_0}{dt}=-\lambda\frac{s_0}{s_0+\kappa}\quad\text{with }s_0(0)=1,
\end{equation}
see also \eqref{4d}.


This approximation
is {\it 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 \eqref{42a}.


\subsection{The Uniformly Valid Leading-Order Approximation}
\label{section3.3}


In this regime, $\delta$ is of the same size as the rate constant
corresponding to the faster of the two kinetics time scales,
{\it i.e.}, $\delta$ is moderately large.
Therefore, given the expressions for $\hat{\mathbf{u}}_0$ and $\mathbf{u}_0$ in
\eqref{42m} and \eqref{42l},
respectively,
with $s_0(t)$
given by \eqref{sec4s0},
the uniformly valid
leading-order approximation for $\mathbf{u}$ is
\begin{subequations}\label{13c}
\begin{align}
s(t,x) &=s_0(t)+\mathcal{O}({\rm e}^{-\min\{\pi^2,1+\kappa\}\delta t})+\mathcal{O}(\varepsilon), \label{13c1} \\
e(t,x) &=\frac{\kappa}{s_0(t)+\kappa}
+\mathcal{O}({\rm e}^{-\min\{a\pi^2,1+\kappa\}\delta t})+\mathcal{O}(\varepsilon), \label{13c2} \\
c(t,x) &=\frac{s_0(t)}{s_0(t)+\kappa}
+\mathcal{O}({\rm e}^{-\min\{b\pi^2,1+\kappa\}\delta t})+\mathcal{O}(\varepsilon), \label{13c3}
\end{align}
\end{subequations}
where $\delta=1/\varepsilon$.
Again, we emphasize that
this approximation
is {\it 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.


\subsection{Higher-Order Corrections}
\label{section3.4}

Similar reasoning as used
in Section 3.2
can be applied to determine the asymptotics of $\mathbf{u}_n$,
for any $n \ge 1$.
In particular, one can show that $\mathbf{u}_n$
remains spatially uniform
for all times $t$, $\mathbf{u}_n=\mathbf{u}_n(t)$, and that the state of the system at time
$t$ is determined by the set of equations
\begin{subequations}\label{eq41}
\begin{gather}
\frac{ds_n(t)}{dt} =[\mathbf{F}_{1,\varepsilon}(\mathbf{u}(t))]_n, \\
e_n(t)+c_n(t) =0, \\
s_0(t)e_n(t)-\kappa c_n(t) =\dot{c}_{n-1}+r_{n-1}(t).
\end{gather}
\end{subequations}
Here, the overdot denotes
differentiation with respect to time,
$\mathbf{F}_{1,\varepsilon}(\mathbf{u})=-se+(\kappa-\lambda)c$, and the term
$r_{n-1}=[\mathbf{F}_{2,\varepsilon}(\mathbf{u})]_n+s_0e_n-\kappa c_n$,
with $\mathbf{F}_{2,\varepsilon}(\mathbf{u})=-se+\kappa c$,
is a function of $\mathbf{u}_1,\ldots,\mathbf{u}_{n-1},
s_n$, and $e_0$ exclusively.

The proof is by induction.
First,
at $\mathcal{O}(\varepsilon^n)$,
\eqref{42a1} reads
\[
\frac{ds_{n-1}(t)}{dt} =[\mathbf{F}_{1,\varepsilon}(\mathbf{u})]_{n-1}(t)+
\frac{\partial^2s_n(t,x)}{\partial x^2}.
\]
Now, $ds_{n-1}/dt=[\mathbf{F}_{1,\varepsilon}(\mathbf{u})]_{n-1}$ by the induction hypothesis, and thus
$\partial^2s_n/\partial x^2=0$.
Using the boundary conditions, then, we find
that $s_n(t,x)=s_n(t)$.

Next,
at $\mathcal{O}(\varepsilon^n)$,
\eqref{42a2} and \eqref{42a3} yield
\begin{subequations}
\begin{gather*}
\dot{e}_{n-1}-r_{n-1} =-(s_0e_n-\kappa c_n)
+a\frac{\partial^2e_n}{\partial x^2}, \\
\dot{c}_{n-1}+r_{n-1} =s_0e_n-\kappa c_n
+b\frac{\partial^2c_n}{\partial x^2},
\end{gather*}
\end{subequations}
Here again, we rescale $x$
via $\bar{x}=x/\sqrt{a}$
and rewrite these equations
as a four-dimensional,
inhomogeneous, linear system
with $t$--dependent coefficients,
\begin{subequations}\label{eq42}
\begin{align}
e_n' &=f_n, \\
f_n' &=s_0e_n-\kappa c_n+\dot{e}_{n-1}-r_{n-1}, \\
c_n' &=d_n, \\
d_n' &=-\alpha(s_0e_n-\kappa c_n+\dot{e}_{n-1}-r_{n-1}),
\end{align}
\end{subequations}
where the prime denotes (partial) differentiation with respect to $\bar{x}$
and we have used that $\dot{e}_{n-1}(t)+\dot{c}_{n-1}(t)=0$
by the induction hypothesis.

The general solution of \eqref{eq42} is given by
\begin{subequations}\label{eq43}
\begin{align}
\begin{split}
e_n (t,x) &=\kappa\gamma_1(t)+\kappa\gamma_2(t)\bar{x}
+\gamma_3(t) {\rm e}^{\sqrt{s_0(t)+\alpha\kappa}\,\bar{x}}
-\gamma_4(t) {\rm e}^{-\sqrt{s_0(t)+\alpha\kappa}\,\bar{x}} \\
&\quad +\frac{\dot{e}_{n-1}(t)-r_{n-1}(t)}{s_0(t)+\alpha\kappa}
   \Big(\cosh(\sqrt{s_0(t)+\alpha\kappa}\,\bar{x})-1\Big),
\end{split} \label{eq431} \\
\begin{split}
c_n (t,x) &=s_0(t)\gamma_1(t)+s_0(t)\gamma_2(t)\bar{x}
-\alpha\gamma_3(t) {\rm e}^{\sqrt{s_0(t)+\alpha\kappa}\,\bar{x}}
+\alpha\gamma_4(t) {\rm e}^{-\sqrt{s_0(t)+\alpha\kappa}\,\bar{x}} \\
&\quad -\alpha\frac{\dot{e}_{n-1}(t)-r_{n-1}(t)}{s_0(t)+\alpha\kappa}
\Big(\cosh(\sqrt{s_0(t)+\alpha\kappa}\,\bar{x})-1\Big),
\end{split} \label{eq432}
\end{align}
\end{subequations}
where $f_n=e_n'$ and $d_n=c_n'$ can be found from \eqref{eq43} and
$\gamma_1,\dots,\gamma_4$ are $t$--dependent constants of integration.
Making use of the Neumann boundary conditions on $e_n$ and $c_n$,
one sees that $\gamma_2$, $\gamma_3$, and $\gamma_4$ are such that
\begin{subequations}\label{eq44}
\begin{gather}
e_n (t,x)=e_n(t)
=\kappa\gamma_1(t)-\frac{\dot{e}_{n-1}(t)-r_{n-1}(t)}{s_0(t)+\alpha\kappa},
 \label{eq441} \\
c_n (t,x)=c_n(t)
=s_0(t)\gamma_1(t)+\alpha\frac{\dot{e}_{n-1}(t)-r_{n-1}(t)}{s_0(t)
+\alpha\kappa}.
 \label{eq442}
\end{gather}
\end{subequations}
In summary, these formulas show that $\mathbf{u}_n=\mathbf{u}_n(t)$.

Next, by integrating the $\mathcal{O}(\varepsilon^{n+1})$ terms
in \eqref{42a1}
and using the boundary conditions,
we obtain
$\int_0^1(ds_{n}/dt)\ dx=\int_0^1[\mathbf{F}_{1,\varepsilon}(\mathbf{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, $[\mathbf{F}_{1,\varepsilon}(\mathbf{u})]_{n}$
is a function
of $\mathbf{u}_0,\ldots,\mathbf{u}_{n-1}$
and thus also a function of time only.
Therefore,
$ds_{n}/dt=[\mathbf{F}_{1,\varepsilon}(\mathbf{u})]_{n}$,
as desired.

Finally,
at $\mathcal{O}(\varepsilon^{n+1})$,
\eqref{42a2} and \eqref{42a3} yield
\begin{gather*}
\dot{e}_n =[\mathbf{F}_{2,\varepsilon}(\mathbf{u})]_n
+a\frac{\partial^2e_{n+1}}{\partial x^2}, \\
\dot{c}_n =-[\mathbf{F}_{2,\varepsilon}(\mathbf{u})]_n
+b\frac{\partial^2c_{n+1}}{\partial x^2}.
\end{gather*}
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
$\dot{e}_n(t)+\dot{c}_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 \eqref{eq44}
with the identity
$\dot{e}_n(t)+\dot{c}_n(t)=0$,
we obtain \eqref{eq41},
and the proof is complete.

\begin{remark} \label{tmk3.1} \rm
Equations \eqref{eq41} 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 $\varepsilon$. 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 $\varepsilon$~\cite{K1999}) also apply to this case.
\end{remark}


\section{Small Diffusivities}\label{section4}

In this section, we consider the regime
in which the species' diffusivities are small,
$\delta=\varepsilon$, 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 $\mathcal{O}(\varepsilon^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 \eqref{1f} with $\delta=\varepsilon$ are
\begin{subequations}\label{5a}
\begin{align}
\varepsilon\frac{\partial s}{\partial t} &=-\varepsilon se+\varepsilon(\kappa-\lambda)c
+\varepsilon^2\frac{\partial^2s}{\partial x^2}, \label{5a1} \\
\varepsilon\frac{\partial e}{\partial t} &=-(se-\kappa c)+a\varepsilon^2
\frac{\partial^2e}{\partial x^2}, \\
\varepsilon\frac{\partial c}{\partial t} &=se-\kappa c+b\varepsilon^2
\frac{\partial^2c}{\partial x^2}.
\end{align}
\end{subequations}
In vector notation, \eqref{5a} is given by
\begin{equation}
\label{5b}
\varepsilon\frac{\partial\mathbf{u}}{\partial t}=\mathbf{F}_\varepsilon(\mathbf{u})+\varepsilon^2 D \frac{\partial^2\mathbf{u}}
{\partial x^2}.
\end{equation}
One factor
of $\varepsilon$ in \eqref{5a1}
is redundant;
however, we retain it
for consistency of notation.


\subsection{Pointwise Fast Kinetics: The Fast (Inner) Time Scale}
\label{section4.1}

On the fast ($\eta=t/\varepsilon$)
time scale, the full governing equations
are given by
\begin{equation}\label{51a}
\frac{\partial\hat{\mathbf{u}}}{\partial\eta}=\mathbf{F}_\varepsilon(\hat{\mathbf{u}})+\varepsilon^2D
\frac{\partial^2\hat{\mathbf{u}}}{\partial x^2},
\end{equation}
subject to the usual initial and boundary conditions, see
Section~\ref{section1}. Again, we expand
$\hat{\mathbf{u}}(\eta,x,\varepsilon)=\hat{\mathbf{u}}_0(\eta,x)+\varepsilon\hat{\mathbf{u}}_1(\eta,x)+\mathcal{O}(\varepsilon^2)$.
To lowest order, \eqref{51a} implies
\begin{equation}\label{5c}
\mathcal{O}(1):\quad\frac{\partial\hat{\mathbf{u}}_0}{\partial\eta}=[\mathbf{F}_\varepsilon(\hat{\mathbf{u}})]_0,
\end{equation}
subject to Neumann boundary conditions and the initial condition
\begin{equation}\label{5d}
\hat{\mathbf{u}}_0(0,x)=\mathbf{u}_i(x).
\end{equation}
Solving the system \eqref{5c} and \eqref{5d} componentwise,
we find $\hat{s}_0(\eta,x)=
\hat{s}_0(0,x)$; hence, $\hat{s}_0(\eta,x)\equiv s_i(x)$
for all $\eta \ge 0$. Moreover,
\[
\frac{\partial\hat{e}_0}{\partial\eta}+\frac{\partial\hat{c}_0}
{\partial\eta}=0,
\]
so that
\begin{equation} \label{5ec}
\hat{e}_0(\eta,x)+\hat{c}_0(\eta,x)=\hat{e}_0(0,x)+\hat{c}_0(0,x)\equiv
e_i(x)\quad\text{for all }\eta\ge 0.
\end{equation}
Now,
given that
$\hat{s}_0=s_i$ and $\hat{e}_0=e_i-\hat{c}_0$,
we see
from \eqref{5c} that
\[
\frac{\partial\hat{c}_0}{\partial\eta}=-(s_i+\kappa)\hat{c}_0+s_ie_i\quad
\text{with }\hat{c}_0(0,x)=0,
\]
which has the solution
\[
\hat{c}_0(\eta,x)=\frac{s_i(x)e_i(x)}{s_i(x)+\kappa}\big(1-{\rm e}^{-(s_i(x)
+\kappa)\eta}\big).
\]
In summary, we have
\begin{equation} \label{5e}
\hat{\mathbf{u}}_0(\eta,x)=\bigg(s_i(x),\frac{e_i(x)}{s_i(x)+\kappa}
\big(\kappa+s_i(x){\rm e}^{-(s_i(x)+\kappa)\eta}\big),
\frac{s_i(x)e_i(x)}{s_i(x)+\kappa}\big(1-{\rm e}^{-(s_i(x)+\kappa)\eta}\big)
\bigg)^T.
\end{equation}

The physical interpretation of these results is that,
on the fast time scale,
diffusion has no effect on the concentration amplitudes
to $\mathcal{O}(1)$ and $\mathcal{O}(\varepsilon)$; it only influences the amplitudes
at $\mathcal{O}(\varepsilon^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 $\varepsilon$, $\hat{s}$ remains
unaltered, and $\hat{e}$ and $\hat{c}$ are constrained to satisfy
the linear conservation law \eqref{5ec},
while the concentrations evolve to the appropriate
($x$--dependent) point on the
leading-order slow manifold as $\eta\to\infty$, which is
exactly dictated by the chemical kinetics.


\subsection{Pointwise Slow Reduced Kinetics: The Slow Time Scale}
\label{section4.2}


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~\ref{section2}, we replace $\mathbf{u}$ by
$\tilde{\mathbf{u}}=\tilde{\mathbf{u}}(t,x,\varepsilon)$ in \eqref{5b} and make the Ansatz
$\tilde{\mathbf{u}}(t,x,\varepsilon) = \tilde{\mathbf{u}}_0(t,x) + \varepsilon\tilde{\mathbf{u}}_1(t,x)
+ \mathcal{O}(\varepsilon^2)$.
After expanding with respect to $\varepsilon$ and rearranging, we have
\begin{subequations}\label{5f}
\begin{align}
\mathcal{O}(1):\quad &[\mathbf{F}_\varepsilon(\tilde{\mathbf{u}})]_0=0, \label{5f1} \\
\mathcal{O}(\varepsilon):\quad &\frac{\partial\tilde{\mathbf{u}}_0}{\partial t}=
[\mathbf{F}_\varepsilon(\tilde{\mathbf{u}})]_1, \label{5f2}
\end{align}
\end{subequations}
with Neumann boundary conditions on $\tilde{\mathbf{u}}_0$ and $\tilde{\mathbf{u}}_1$
and with initial conditions to be determined by matching with the equations
on the fast time scale.

Writing \eqref{5f1} componentwise, we see that
$\tilde{s}_0\tilde{e}_0=\kappa\tilde{c}_0$, which we can substitute into
\eqref{5f2} to obtain
\[
\frac{\partial\tilde{s}_0}{\partial t}=-\lambda\tilde{c}_0.
\]
Also, we deduce from \eqref{5f2} that $\tilde{e}_0(t,x)+\tilde{c}_0(t,x)=
\tilde{e}_0(0)+\tilde{c}_0(0)$ must hold. Imposing the matching condition
$\lim_{\eta\to\infty}\hat{\mathbf{u}}_0(\eta,x)=\lim_{t\to 0^+}\tilde{\mathbf{u}}_0(t,x)$
and taking into consideration that
\begin{equation} \label{5i}
\lim_{\eta\to\infty}\hat{\mathbf{u}}_0(\eta,x)=\bigg(s_i(x),\frac{\kappa e_i(x)}
{s_i(x)+\kappa},\frac{s_i(x)e_i(x)}{s_i(x)+\kappa}\bigg)^T,
\end{equation}
we conclude that $\tilde{e}_0(t,x)+\tilde{c}_0(t,x)\equiv e_i(x)$.
Therefore, recalling \eqref{5f1}, we obtain
\begin{equation} \label{5g}
\tilde{e}_0(t,x)=\frac{\kappa e_i(x)}{\tilde{s}_0(t,x)+\kappa}\quad
\text{and}\quad\tilde{c}_0(t,x)=\frac{\tilde{s}_0(t,x)e_i(x)}
{\tilde{s}_0(t,x)+\kappa},
\end{equation}
where $\tilde{s}_0(t,x)$ solves
\begin{equation}\label{5h}
\frac{\partial\tilde{s}_0}{\partial t}=-\lambda\frac{e_i(x)}
{\tilde{s}_0+\kappa}\tilde{s}_0\quad\text{with }\tilde{s}_0(0,x)=s_i(x).
\end{equation}
This reduced equation
for $\tilde{s}_0(t,x)$
is of the same form as \eqref{4d},
the equation for $s_0(t)$
in the large diffusivities regime;
however,
the interpretations differ.
Here,
$\tilde{s}_0(t,x)$ depends on $x$
through $e_i(x)$,
whereas $s_0(t)$ is independent of $x$.
Moreover, we recall that $\tilde{s}_0$
may be found only implicitly,
as was shown
for $s_0(t)$ in \eqref{4g}.

Therefore, we conclude
that, at each point $x$ in the domain,
$\tilde{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 $\tilde{e}_0$ and $\tilde{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 $\tilde{\mathbf{u}}$
takes place outside the unit cube, {\it i.e.}, outside the region
that is feasible for the classical kinetics problem.

\subsection{Homogenization of the Enzyme Concentration:
The Super-Slow \\ (Outer) Time Scale}
\label{section4.3}

Due to the fact
that $\delta=\varepsilon$
in this regime,
diffusion acts
on the super-slow
($\zeta=\varepsilon t$) time scale,
and we express
the governing
equations as
\begin{equation} \label{6b}
\varepsilon^2\frac{\partial\mathbf{u}}{\partial\zeta}=\mathbf{F}_\varepsilon(\mathbf{u})+\varepsilon^2D
\frac{\partial^2\mathbf{u}}{\partial x^2},
\end{equation}
subject to Neumann boundary conditions.
We make the Ansatz
$$
\mathbf{u}(\zeta,x,\varepsilon)=\mathbf{u}_0(\zeta,x)+\varepsilon\mathbf{u}_1(\zeta,x)+\varepsilon^2\mathbf{u}_2(\zeta,x)
+\mathcal{O}(\varepsilon^3),
$$
substitute it into \eqref{6b}, expand, and rearrange
the resulting equations to obtain
\begin{subequations}
\begin{align}
\mathcal{O}(1):\quad &[\mathbf{F}_\varepsilon(\mathbf{u})]_0=0, \label{6c1} \\
\mathcal{O}(\varepsilon):\quad &[\mathbf{F}_\varepsilon(\mathbf{u})]_1=0, \label{6c2} \\
\mathcal{O}(\varepsilon^2):\quad &\mathcal{L}_\zeta\mathbf{u}_0=[\mathbf{F}_\varepsilon(\mathbf{u})]_2, \label{6c3}
\end{align}
\end{subequations}
where
$\mathcal{L}_\zeta =\frac{\partial}{\partial\zeta}-D\frac{\partial^2}{\partial x^2}$,
supplemented by Neumann boundary conditions.

First, observe that \eqref{6c1} gives $s_0e_0=\kappa c_0$, as
on the slow ($t$) scale. Then \eqref{6c2} shows that
$\lambda c_0=0$, which implies $c_0\equiv 0$ (since $\lambda\ne 0$ by
definition). Moreover, we see from that same equation \eqref{6c2} that
\[
s_0e_1+s_1e_0=\kappa c_1,
\]
which we use in \eqref{6c3} to obtain
\[
\frac{\partial s_0}{\partial\zeta}-\frac{\partial^2s_0}{\partial x^2}=
-\lambda c_1.
\]
Similarly, to find $e_0$, we note that \eqref{6c3} and the identity
$c_0\equiv 0$ derived above imply $\mathcal{L}_\zeta c_0=0$; hence,
\[
s_0e_2+s_1e_1+s_2e_0=\kappa c_2.
\]
Therefore, $\mathcal{L}_\zeta e_0=0$ and
\[
e_0(\zeta,x)=\sum_{k=0}^\infty e_{0k}{\rm e}^{-a(k\pi)^2\zeta}\cos{(k\pi x)},
\]
with the constant coefficients $e_{0k}$ to be determined by matching
with the dynamics on the $t$ scale. Since matching requires
$\lim_{t\to\infty}\tilde{\mathbf{u}}_0(t,x)=\lim_{\zeta\to 0^+}\mathbf{u}_0(\zeta,x)$
and 
\begin{equation} \label{6e}
\lim_{t\to\infty}\tilde{\mathbf{u}}_0(t,x)=\big(0,e_i(x),0\big)^T
\end{equation}
by \eqref{5g} and \eqref{5h}, 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}=\int_0^1e_i(x)\,dx =1,
\]
it follows that $e_0(\zeta,x)=1+\mathcal{O}({\rm e}^{-a\pi^2\zeta})$.
All reactions have to leading order been completed
when diffusion sets in, {\it i.e.}, there is only enzyme left to diffuse
until a spatially homogeneous distribution of $e$ has been reached.

Finally, given $e_0(x)\not\equiv 0$ and \eqref{6c1}, we conclude that
$s_0\equiv 0$ must hold. Hence, to summarize, we have
\begin{equation} \label{6f}
\mathbf{u}_0(\zeta,x)=(0,1,0)^T+\mathcal{O}({\rm e}^{-d\pi^2\zeta}),
\end{equation}
where again $d=\min\{1,a,b\}$.

\subsection{The Uniformly Valid Leading-Order Approximation}
\label{section4.4}

We combine the expressions for
$\hat{\mathbf{u}}_0$, $\tilde{\mathbf{u}}_0$, and $\mathbf{u}_0$, which are the leading-order
approximations to the solution $\mathbf{u}$ of \eqref{5b} on the fast, slow,
and super-slow time scales, respectively, into one uniformly valid,
leading-order approximation
for $\mathbf{u}$. Given \eqref{5e}, \eqref{5g} and \eqref{5h}, and \eqref{6f}, as
well as \eqref{5i} and \eqref{6e}, a
straightforward calculation shows
\begin{subequations}\label{13b}
\begin{gather}
s(t,x) = \tilde{s}_0(t,x)+\mathcal{O}({\rm e}^{-\pi^2\delta t})
          +\mathcal{O}(\varepsilon), \label{13b1} \\
e(t,x) = 1-\frac{\tilde{s}_0(t,x)e_i(x)}{\tilde{s}_0(t,x)+\kappa}
         +\frac{s_i(x)e_i(x)}{s_i(x)+\kappa}{\rm e}^{-(s_i(x)+\kappa)t/\delta}
         +\mathcal{O}({\rm e}^{-a\pi^2\delta t})
         +\mathcal{O}(\varepsilon), \label{13b2} \\
c(t,x) = \frac{\tilde{s}_0(t,x)e_i(x)}{\tilde{s}_0(t,x)+\kappa}
         -\frac{s_i(x)e_i(x)}{s_i(x)+\kappa}{\rm e}^{-(s_i(x)+\kappa)t/\delta}
         +\mathcal{O}({\rm e}^{-b\pi^2\delta t})
         +\mathcal{O}(\varepsilon), \label{13b3}
\end{gather}
\end{subequations}
where $\tilde{s}_0$ is defined by \eqref{5h}.
Here, we also recall
that $\delta=\varepsilon$, {\it i.e.,} $\delta$ is asymptotically smaller
than the smaller of the two kinetics rate constants.



\subsection{Higher-Order Corrections}
\label{section4.5}


In this regime, the higher-order algebraic corrections $\mathbf{u}_j$, $j\ge 1$,
may be found explicitly.

The $\mathcal{O}(\varepsilon)$ terms in \eqref{51a} are
\begin{equation}\label{5ho1}
\mathcal{O}(\varepsilon):\quad\frac{\partial\hat{\mathbf{u}}_1}{\partial\eta}=[\mathbf{F}_\varepsilon(\hat{\mathbf{u}})]_1.
\end{equation}
In particular,
the equation for the first component, $\hat{s}_1$, is
\begin{equation}\label{5ho2}
\frac{\partial \hat{s}_1}{\partial \eta}
= -\hat{s}_0 \hat{e}_0
  + (\kappa -\lambda) \hat{c}_0
=\frac{ - s_i(x) e_i(x)}{s_i(x) + \kappa}
\Big( \lambda + e^{-(s_i(x) + \kappa)\eta} (s_i(x)+\kappa -\lambda) \Big),
\end{equation}
where we used \eqref{5e}.
Hence, there will be a secular term in $\hat{s}_1$
that grows linearly in $\eta$,
due to the $\eta$--independent (first) term
in \eqref{5ho2}.

A remedy is at hand using the method of multiple scales.
In particular,
let the terms in the
expansion of $\hat{\mathbf{u}}$ depend not only
on the fast time $\eta$, 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 $\mathcal{O}(1)$ term, $A_0(t,x)\hat{s}_0(\eta,x)$,
in the expansion for $\hat{s}$.
As in Section~\ref{section4.1},
it is again the case that
$\hat{s}_0 = s_i(x)$, independently of $\eta$,
and that \eqref{5ec} holds.
Moreover, the leading-order terms
for the complex and enzyme concentrations
may be found as in Section~\ref{section4.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)+\kappa}
\big(1-{\rm e}^{-(A_0(t,x)s_i(x)+\kappa)\eta}\big).
\]
Substituting these solutions
into \eqref{51a}
and collecting
the $\mathcal{O}(\varepsilon)$ terms,
we find
\begin{equation}\label{5ho3}
s_i(x) \frac{\partial A_0}{\partial t}
+ A_1 \frac{\partial \hat{s}_1}{\partial \eta}
= - \frac{s_i(x) e_i(x) A_0}{A_0 s_i(x) + \kappa}
  \Big( \lambda
         + (A_0 s_i(x)+\kappa-\lambda)e^{-(A_0 s_i(x)+\kappa)\eta}
  \Big).
\end{equation}
Therefore,
if we choose $A_0(t,x)$
such that the $\eta$--independent term
vanishes from the right hand side
of \eqref{5ho3},
then we obtain a regular (bounded) solution
for $\hat{s}_1$.
Specifically,
we choose the dependence of $A_0(t,x)$
on $t$ so that
\begin{equation} \label{5ho4}
\frac{\partial A_0}{\partial t}
= -\lambda \frac{e_i(x) A_0}{A_0 s_i(x) + \kappa}.
\end{equation}
This requirement
is equivalent to applying
a solvability condition
to equation \eqref{5ho3},
which forces the `bad' terms
to drop out.
The resulting equation
for $\hat{s}_1$
may be integrated directly to obtain
a bounded function.
Moreover, as one expects,
the equation \eqref{5ho4}
for $A_0(t,x)$ is precisely
equation \eqref{5h},
which was obtained
from the leading-order analysis
on the slow time scale,
with $\tilde{s}_0(t,x)$
replaced by $A_0(t,x)s_i(x)$.

Higher-order terms $\hat{s}_j$,
$j \ge 1$,
may be found in an analogous manner.
The corresponding amplitude functions
$A_i$ (for $i=1,\ldots,j$)
depend on
$t$, $\varepsilon t$, $\ldots$,
$\varepsilon^{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 $\hat{s}_1$ above.
Of course,
the corrections
to the complex and enzyme concentrations
need to be determined simultaneously
using the multiple scales Ansatz.


\section{Interpretation and Conclusions}\label{section5}

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.

{\bf Regime (i): large diffusivities,
$\delta=1/\varepsilon^2$}.
There is a brief initial
transient period of $\mathcal{O}(\varepsilon^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, $\mathcal{O}(\varepsilon)$--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 \eqref{13a}.

The results of a typical simulation of
\eqref{1f} for this regime are shown in Figure~\ref{figure1}. In particular,
the super-fast initial transient
appears to last for essentially only one
timestep, as is evident from Figures~\ref{figure1a} and \ref{figure1b}.

\begin{figure}[ht!]
\centering
\subfigure[]
{
\includegraphics[scale=0.5]{figures/fig1a} \label{figure1a}
}
\subfigure[]
{
\includegraphics[scale=0.5]{figures/fig1b} \label{figure1b}
}
\subfigure[]
{
\includegraphics[scale=0.5]{figures/fig1c} \label{figure1c}
}
\caption{Regime (i): large diffusivities.
Concentration profiles of (a) $s$, (b) $e$, and
(c) $c$ as functions of $t\ge 0$ and $x\in [0,1]$,
with parameter values $\kappa=3$, $\lambda=5$, $\varepsilon=0.1$,
$\delta=100$, $a=2$, $b=3$, and initial profiles
$(s,e,c)(0,x)=(1+\textstyle{\frac{1}{2}}\cos(2\pi x),
1+\textstyle{\frac{1}{2}}\cos(\pi x),0)$.}
\label{figure1}
\end{figure}


{\bf Regime (ii):
moderately large diffusivities,
$\delta=1/\varepsilon$}.
The substrate concentration
again homogenizes during a brief initial transient,
which is of
$\mathcal{O}(\varepsilon)$ 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 \eqref{13c}.

The results of a typical simulation of
\eqref{1f} are shown in Figure~\ref{figure2}.
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 $\mathcal{O}(\varepsilon)$,
as compared to $\mathcal{O}(\varepsilon^2)$ before.
See especially Figure~\ref{figure2b},
which suggests that the
number of timesteps is approximately five.

\begin{figure}[ht!]
\centering
\subfigure[]
{
\includegraphics[scale=0.5]{figures/fig2a} \label{figure2a}
}
\subfigure[]
{
\includegraphics[scale=0.5]{figures/fig2b} \label{figure2b}
}
\subfigure[]
{
\includegraphics[scale=0.5]{figures/fig2c} \label{figure2c}
}
\caption{Regime (ii): moderately large diffusivities.
Concentration profiles of (a) $s$, (b) $e$, and
(c) $c$ as functions of $t\ge 0$ and $x\in [0,1]$,
with parameter values
$\kappa=3$, $\lambda=5$, $\varepsilon=0.1$, $\delta=10$, $a=2$, $b=0.5$, and
initial profiles $(s,e,c)(0,x)=(1+\textstyle{\frac{1}{2}}\cos(2\pi x),
1+\textstyle{\frac{1}{2}}\cos(\pi x),0)$.}
\label{figure2}
\end{figure}

{\bf Regime (iii): small diffusivities,
$\delta=\varepsilon$}. At every point
in the domain,
the kinetics dominate on fast and slow time intervals on the
orders of $\varepsilon$ 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 converging
to the spatial average of the initial enzyme concentration profile.
The leading-order approximation,
valid uniformly in time and space,
is given in \eqref{13b}.

The results of a typical simulation of
\eqref{1f} for this regime are shown in Figure~\ref{figure3}. It
takes at least ten timesteps before the concentrations are essentially
homogeneous.

\begin{figure}[ht!]
\centering
\subfigure[]
{
\includegraphics[scale=0.5]{figures/fig3a}\label{figure3a}
}
\subfigure[]
{
\includegraphics[scale=0.5]{figures/fig3b}\label{figure3b}
}
\subfigure[]
{
\includegraphics[scale=0.5]{figures/fig3c}\label{figure3c}
}
\caption{Regime (iii): small diffusivities.
Concentration profiles of (a) $s$, (b) $e$, and
(c) $c$ as functions of $t\ge 0$ and $x\in[0,1]$,
with parameter values $\kappa=3$, $\lambda=5$, $\varepsilon=0.1$,
$\delta=0.1$, $a=2$, $b=3$, and initial profiles
$(s,e,c)(0,x)=(1+\textstyle{\frac{1}{2}}\cos(2\pi x),
1+\textstyle{\frac{1}{2}}\cos(\pi x),0)$.}
\label{figure3}
\end{figure}

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~\cite{EK2005}, see also Section~\ref{section1}. 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, {\it 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 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~\cite{BDM1984},
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 

\section{Linear stability of the extinguished state}
\label{appa}


In this appendix, we state and prove
the linear (spectral) stability of the spatially uniform
(homogeneous) equilibria $\mathbf{u}^\ast=(0,e^\ast,0)$ to which solutions
tend in the limit $t\to\infty$, {\it 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~\cite{T1952}; for more general considerations on the
stability of homogeneous solutions in the presence of diffusion, we refer
to \cite{DLAS2004,MS2004}.

\begin{lemma} \label{A.1}
For any $e^\ast\ge 0$ constant, the homogeneous equilibrium solution
$\mathbf{u}^\ast=(0,e^\ast,0)$ of \eqref{1h} is linearly stable.
\end{lemma}

\begin{proof} The argument follows closely the corresponding discussion of
morphogenesis in~\cite[Section~11.4]{EK2005}. First, we linearize
\eqref{1f} about $\mathbf{u}=\mathbf{u}^\ast$. Let $\mathbf{u}'=\mathbf{u}-\mathbf{u}^\ast$; then,
\begin{equation} \label{1j}
\frac{\partial\mathbf{u}'}{\partial t}=\frac 1\varepsilon {\rm d}\mathbf{F}_\varepsilon(\mathbf{u}^\ast)\mathbf{u}'+\delta D
\frac{\partial^2\mathbf{u}'}{\partial x^2},
\end{equation}
where $D={\rm diag}(1,a,b)$ is defined as above and ${\rm d}\mathbf{F}_\varepsilon(\mathbf{u}^\ast)$
is given by
\[
{\rm d}\mathbf{F}_\varepsilon(\mathbf{u}^\ast)=\begin{bmatrix}
-\varepsilon e^\ast & 0 & \varepsilon(\kappa -\lambda) \\
-e^\ast & 0 & \kappa \\
e^\ast & 0 & -\kappa
\end{bmatrix}.
\]
Now, we make the Ansatz $\mathbf{u}'=(s',e',c')^T=(\alpha_1,\alpha_2,\alpha_3)^T
\cos{(k\pi x)}{\rm e}^{\sigma t}$, which we substitute into \eqref{1j}. After
dividing out the common factor $\cos{(k\pi x)}{\rm e}^{\sigma t}$ and
rearranging the resulting equations, we obtain the following
homogeneous system of linear equations for $(\alpha_1,\alpha_2,\alpha_3)^T$,
\begin{align}\label{1k}
(\sigma I_3-A)\cdot (\alpha_1,\alpha_2,\alpha_3)^T=(0,0,0)^T,
\end{align}
where $A$ is defined as $A={\rm d}\mathbf{F}_\varepsilon(\mathbf{u}^\ast)/\varepsilon-\delta (k\pi)^2D$
and $I_n$ denotes the $n\times n$-identity matrix. For \eqref{1k} to have
nontrivial solutions, $\det {(\sigma I_3-A)}$ must vanish. An elementary
calculation shows
\begin{align*}
\det {(\sigma I_3-A)} &=\big(\sigma +\delta a(k\pi)^2\big)
\Big[\sigma^2+\Big(e^\ast+\delta(k\pi)^2+\frac{\kappa}{\varepsilon}
+\delta b(k\pi)^2\Big)\sigma \\
&\quad +\big(e^\ast+\delta(k\pi)^2\big)
\Big(\frac{\kappa}{\varepsilon}+\delta b(k\pi)^2\Big)-\frac{e^\ast}{\varepsilon}
(\kappa -\lambda)\Big],
\end{align*}
which implies that either $\sigma=-\delta a(k\pi)^2<0$, or that the
term in square brackets vanishes. This term, however, is precisely
the determinant of the $2\times 2$-submatrix $\sigma I_2-\bar{A}$ of
$\sigma I_3-A$ obtained by deleting the second row and column, with
\begin{equation} \label{1l}
\bar{A}=\begin{pmatrix}
-\big(e^\ast+\delta (k\pi)^2\big) & \kappa -\lambda \\
\dfrac{e^\ast}{\varepsilon} & -\Big(\dfrac{\kappa}{\varepsilon}+\delta b(k\pi)^2\Big)
\end{pmatrix}.
\end{equation}
Hence, we are within the
framework of~\cite{EK2005}, and it remains to show that $\mathop{\rm tr}
\bar{A}<0$ and $\det \bar{A}>0$ for $\mathop{\rm Re}(\sigma)<0$ to hold, as
clearly
\[
\sigma =\frac{\mathop{\rm tr}\bar{A}}{2}\pm\sqrt{\Big(\frac{\mathop{\rm tr}\bar{A}}
{2}\Big)^2-\det \bar{A}}.
\]
However, this is immediate from \eqref{1l}, since
\[
\mathop{\rm tr}\bar{A}=-\Big( e^\ast+\frac{\kappa}{\varepsilon}
+\delta (k\pi)^2(1+b)\Big)
\]
and
\[
\det {\bar{A}}=\frac1\varepsilon\big(\kappa\delta(k\pi)^2+\lambda\big)
+\delta b(k\pi)^2\big(e^\ast+\delta (k\pi)^2\big),
\]
and since the values of all individual parameters in these
expressions are non-negative.
\end{proof}

\begin{remark} \label{rmkA.2} \rm
The above conclusion also holds for $\kappa=0$,
since $\mathop{\rm tr}\bar{A}$ remains negative and $\det \bar{A}$
positive.
\end{remark}

\section{Estimates on \eqref{55c} and \eqref{55c-generalk} for regime (ii)}
\label{appb}

In this appendix,
we present some of the technical steps
involved in analyzing
the dynamics of $\hat e_{0k}$ and $\hat c_{0k}$,
\eqref{55c} and \eqref{55c-generalk},
and of deriving \eqref{42m},
for regime (ii).

As shown in Section~\ref{section3},
we know that
$\hat e_{00}(\eta)+\hat c_{00}(\eta)=1$
for all $\eta>0$.
Thus, the equation \eqref{55c}
for $\hat e_{00}$
may be rewritten as
\[
  \frac{d}{d\eta}\Big(\hat e_{00}-\frac{\kappa}{\kappa+1}\Big)
=   -(1+\kappa)\Big(\hat e_{00}-\frac{\kappa}{\kappa+1}\Big)-\mathcal{F}_0.
\]
Writing $w=\hat e_{00}-\kappa/(\kappa+1)$ and multiplying both members by $w$,
we obtain
\[
  \frac12\frac{d(w^2)}{d\eta} =
  -(1+\kappa)w^2-\mathcal{F}_0w.
\]
Applying Young's inequality
to the last term
in the right member,
we find
\[
  \frac12\frac{d(w^2)}{d\eta} \le
  -(1+\kappa-\sigma)w^2+\frac{1}{4\sigma}\mathcal{F}_0^2.
\]
(Here, $\sigma>0$
is a suitably chosen parameter.)
Next, application
of Gronwall's inequality yields
\begin{equation} \label{021}
  (w(\eta))^2 \le (w(0))^2{\rm e}^{-2(1+\kappa-\sigma)\eta}
  +\frac{1}{2\sigma}\int_{0}^\eta
       {\rm e}^{-2(1+\kappa-\sigma)(\eta-s)}(\mathcal{F}_0(s))^2\,ds.
\end{equation}
Recalling the definition
of $\mathcal{F}_0$ and
that $\hat s_0 = s_i$,
applying H\"older's inequality twice,
and employing Parseval's identity,
we estimate
\begin{equation} \label{022}
  (\mathcal{F}_0(\eta))^2
  =\frac12{\rm e}^{-2\pi^2\eta}\Big(\sum_{m\ge1}\hat s_{0m}\hat e_{0m}(\eta)
                           \Big)^2
  \le \frac12{\rm e}^{-2\pi^2\eta}\Vert{s_i}\Vert_2^2\Vert{\hat e_0(\eta)}\Vert_2^2.
\end{equation}
Substituting this estimate into the integral in \eqref{021},
we calculate
\[
  (w(\eta))^2 \le \left((w(0))^2-C(\sigma)\right)
{\rm e}^{-2(1+\kappa-\sigma)\eta}+C(\sigma){\rm e}^{-2\pi^2\eta} ,
\]
where
\[
C(\sigma)=\frac{
\Vert{s_i}\Vert_2^2\left(\sup_{\eta\ge0}\Vert{\hat e_0(\eta)}\Vert_2^2\right)
}{8\sigma(1+\kappa-\sigma-\pi^2)} .
\]
Here, $\Vert{\hat e_0(\eta)}\Vert_2$
is guaranteed to remain bounded uniformly in time
by standard results.
This proves that $w(\eta)\to 0$
(or, equivalently, that
$\hat e_{00}(\eta)\to\kappa/(\kappa+1)$)
as $\eta\to\infty$ at an exponential rate.

The same type of estimates can be made to determine the long-term evolution
of the $k$th Fourier modes $\hat e_{0k}(\eta)$ and $\hat c_{0k}(\eta)$.
In particular, one can show that these modes decay to zero exponentially fast.
Hence, we have demonstrated \eqref{42m}.

\subsection*{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\'c was supported by
NSF grant DMS-0109427.

The research of A.~Zagaris was supported by
NWO grant 639.031.617.


\begin{thebibliography}{00}

\bibitem{BAO1963} J. R. Bowen, A. Acrivos, and A. K. Oppenheim;
\emph{Singular perturbation refinement to the quasi-steady state
approximation in chemical kinetics}, Chem. Eng. Sci. \textbf{18} (1963),
177--188.

\bibitem{BDM1984} M. Boudart and G. Dj\'ega-Mariadassou; \emph{Kinetics of
Heterogeneous Catalytic Reactions}, Princeton University Press, Princeton,
1984.

\bibitem{DLAS2004} P. De Leenheer, D. Angeli, and E. D. Sontag;
\emph{Monotone chemical reaction networks}, DIMACS Technical
Report \textbf{16} (2004).

\bibitem{DW1979} M. Dixon and E. C. Webb; \emph{Enzymes}, third edition,
Academic Press, New York, 1979.

\bibitem{EK2005} L. Edelstein-Keshet; \emph{Mathematical models in biology},
Classics in Applied Mathematics \textbf{46}, Society for Industrial and
Applied Mathematics, Philadelphia, 2005. Reprint of the 1988 original.

\bibitem{E2006} 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;
 \emph{Ever-fluctuating
single enzyme molecules: Michaelis-Menten equation revisited}, Nature
Chemical Biology \textbf{2}(2) (2006), 87--94.

\bibitem{FB1997} J. E. Ferrell, Jr. and R. R. Bhatt; \emph{Mechanistic
studies of the dual phosphorylation of mitogenactivated protein kinase},
J. Biol. Chem. \textbf{272}(30) (1997), 19008--19016.

\bibitem{H1982} G. G. Hammes;
\emph{Enzyme Catalysis and Regulation},
Academic Press, New York, 1982.

\bibitem{HTA1967} F. G. Heineken, H. M. Tsuchiya, and R. Aris; \emph{On
the mathematical status of the pseudo-steady state hypothesis of biochemical
kinetics}, Math. Biosci. \textbf{1} (1967), 95--113.

\bibitem{H1903} V. Henri; \emph{Lois g\'en\'erales de l'action des diastases},
Hermann, Paris, 1903.

\bibitem{J1994} C. K. R. T. Jones;
\emph{Geometric Singular Perturbation Theory}
pages 44--118, in \emph{Dynamical Systems, Montecatini Terme, 1994},
R. Johnson, ed.,
LNM \textbf{1609},
Springer, Berlin, 1994.

\bibitem{K1999} T. J. Kaper;
\emph{An introduction to geometric methods
and dynamical systems theory
for singular perturbation problems},
pages 85--132, in \emph{Analyzing Multiscale Phenomena
Using Singular Perturbation Methods},
eds. J. Cronin and R.E. O'Malley, Jr.,
Proc. Symp. Appl. Math, \textbf{56},
American Mathematical Society,
Providence, RI, 1999.

\bibitem{KS1998} J. Keener and J. Sneyd;
\emph{Mathematical Physiology},
Interdisciplinary Applied Mathematics \textbf{8},
Springer-Verlag, New York, 1998.

\bibitem{KC1996} J. Kevorkian and J. D. Cole;
\emph{Multiscale and Singular Perturbation Methods}
Applied Mathematical Sciences \textbf{114},
Springer-Verlag, New York, 1996.

\bibitem{LS1974} C. C. Lin and L. A. Segel; \emph{Mathematics Applied to
Deterministic Problems in the Natural Sciences. With material on elasticity
by G.H. Handelman}, Macmillan Publishing Co., Inc., New York, 1974.

\bibitem{MM1913} L. Michaelis and M. L. Menten; \emph{Die Kinetik der
Invertinwirkung}, Biochem. Zeitsch. \textbf{49} (1913), 333--369.

\bibitem{MZLW2005} P. Miller, A. M. Zhabotinsky, J. E. Lisman, and X.-J. Wang;
\emph{The stability of a stochastic CaMKII Switch: Dependence
on the number of enzyme molecules and protein turnover},
PLoS Biol. \textbf{3}(4) (2005), e107, 705--717.

\bibitem{MS2004} M. Mincheva and D. Siegel; \emph{Stability of mass action
reaction--diffusion systems}, Nonlinear Anal. \textbf{56}(8) (2004),
1105--1131.

\bibitem{M1989} J. D. Murray; \emph{Mathematical Biology}, Biomathematics
\textbf{19}, Springer-Verlag, Berlin, 1989.

\bibitem{OM1991} R. E. O'Malley, Jr.; \emph{Singular Perturbation Methods
for Ordinary Differential Equations}, Applied Mathematical Sciences
\textbf{89}, Springer-Verlag, New York, 1991.

\bibitem{P1987} B. O. Palsson; \emph{On the dynamics of the irreversible
Michaelis-Menten reaction mechanism}, Chem. Eng. Sci. \textbf{42}(3) (1987),
447--458.

\bibitem{PL1984} B. O. Palsson and E. N. Lightfoot; \emph{Mathematical
Modeling of Dynamics and Control in Metabolic Networks. I. On Michaelis-Menten
Kinetics}, J. Theor. Biol. \textbf{111} (1984), 273--302.

\bibitem{SM2003} S. Schnell and P. K. Maini; \emph{A Century of Enzyme
Kinetics: Reliability of the $K_M$ and $v_{\rm max}$ Estimates}, Comm.
Theoret. Biol. \textbf{8} (2003), 169--187.

\bibitem{S1975} I. H. Segel; \emph{Enzyme Kinetics: Analysis
of Rapid Equilibrium and Steady-State Enzyme Systems},
Wiley and Sons, New York, 1975.

\bibitem{S1998} M. Stiefenhofer; \emph{Quasi-steady-state approximation
for chemical reaction networks}, J. Math. Biol. \textbf{36} (1998), 593--609.

\bibitem{SPM1975} S. Strickland, G. Palmer and V. Massey; \emph{Determination
of disassociation constants and specific rate constants of enzyme-substrate
(or protein-ligand) interactions from rapid reaction kinetic data},
J. Biol. Chem. \textbf{250}(11) (1975), 4048--4052.

\bibitem{SS1989} L. A. Segel and M. Slemrod; \emph{The quasi-steady state
assumption: A case study in perturbation}, SIAM Review \textbf{31} (1989),
446--477.

\bibitem{T1952} A. M. Turing; \emph{The chemical basis of morphogenesis},
Phil. Trans. Roy. Soc. Lond. Series B \textbf{237} (1952), 6220--6249.

\bibitem{VBK1995} A. B. Vasil'eva, V. F. Butuzov, and L. V. Kalachev;
\emph{The boundary function method for singular perturbation problems.
With a foreword by Robert E. O'Malley, Jr.}, SIAM Studies in Applied
Mathematics \textbf{14}, Society for Industrial and Applied Mathematics,
Philadelphia, 1995.

\bibitem{WH2004} Y. Wei and M. H. Hecht; \emph{Enzyme-like proteins from
an unselected library of designed amino acid sequences}, Protein Engineering,
Design and Selection \textbf{17}(1), 67--75, 2004.

\bibitem{YTBMP1995}
A. N. Yannacopoulos, A. S. Tomlin, J. Brindley, J. H. Merkin,
and M.J. Pilling;
\emph{The use of algebraic sets in the approximation
of inertial manifolds and lumping in chemical kinetics},
Physica D \textbf{83}, 421--449, 1995.

\end{thebibliography}

\end{document}