Robust numerical method for a singularly perturbed problem arising in the modelling of enzyme kinetics

— A system of two coupled nonlinear initial value equations, arising in the mathematical modelling of enzyme kinetics, is examined. The system is singularly perturbed and one of the components will contain steep gradients. A priori parameter explicit bounds on the two components are established. A numerical method incorporating a specially constructed piecewise-uniform mesh is used to generate numerical approximations, which are shown to converge pointwise to the continuous solution irrespective of the size of the singular perturbation parameter. Numerical results are pre-sented to illustrate the computational performance of the numerical method. The numerical method is also remarkably simple to implement.


I. INTRODUCTION
The Henri-Michaelis-Menten system of nonlinear differential equations arises in the mathematical modelling of enzyme-substrate dynamics, see, for example, [1], [3], [5], [7]. As analytical solutions are not available, it is necessary to solve this system numerically. This may be difficult when the system is singularly perturbed [1]. Asymptotic expansions associated with this singularly perturbed problem are discussed in [9], [7], [8]. Here, to establish a parameter-uniform pointwise error bound on the numerical approximations, we construct a Shishkin decomposition [6] for the solutions. This can be viewed as an alternative to an asymptotic expansion. Based on this decomposition, an efficient finite difference method is constructed, which uses a specially constructed piecewise uniform mesh and an appropriate standard finite difference operator to deal with the steep gradients that appear initially, see for example [6]. This numerical method is shown to be parameter-uniformly convergent, in the sense that its numerical solutions converge to the exact solution, uniformly with respect to the singular perturbation parameter, with essentially first order. A Shishkin mesh was also used in [4] to solve the same system, but it should be noted that the hypotheses [4, (5) and (6)], required by the theory in [4], are not fulfilled by the system considered in the present paper.
Asymptotic expansions yield useful approximations when the singular perturbation parameter ε is sufficiently small for terms of a certain order O(ε n ) to be negligible. However, the accuracy in parameter-uniform approximations is valid irrespective of the size of ε and depends only on the number of mesh points N used in the computations. Moreover, there can be a debate [7], [8] about how to choose the dimensionless variables when using the quasi-steady-state assumption [8], which becomes irrelevant if one has a parameter-uniform numerical method. The Shishkin-mesh method developed below does not require an asymptotic expansion to compute an approximation and, hence, the method is simpler than other approaches which generate individual terms (or approximations) in an asymptotic expansion (e.g., [9]).
The structure of the paper is as follows. In the next section the continuous problem is formulated. In §3 the problem is discretised and the error analysis is performed in §4. The main results are stated in Theorems 5 and 6. The paper concludes with §5 in which some numerical experiments illustrate the form of the solution, and its initial layer, and also support the theoretical error analysis.

II. CONTINUOUS PROBLEM
In the basic model of enzyme reactions (e.g. [7]), a substrate S reacts with an enzyme E to form a complex SE which is converted into a product P and the enzyme. The concentrations of these variables vary with time τ and are denoted here by lower case letters These reactions can be modelled by the following system of four first order equations for four unknowns with given initial conditions where the parameters k −1 , k 1 , k 2 are reaction rate constants. Note that Hence, we have only two unknown variables (s and c) to determine. As in [7], we introduce the following scalings and nondimensionless variables u and v: This leads to the following autonomous system of two coupled nonlinear initial value equations ([7, eq. (6.13)]): We observe the following facts and so there will be an initial layer in v and a weak initial layer in u (as εu (0) = O(1) and u (0) = O(1)). The exact solution of problem (2.2) is unknown.
Lemma 1. For all t > 0, the solutions of (2.2) are bounded as follows: Proof: Given the initial conditions (2.3) and the fact that u , v are continuous, there exists an interval (0, t 1 ), where u (t) < 0 and v (t) > 0 for all t ∈ (0, t 1 ). Hence, there also must exist some 0 < t 0 < t 1 such that 0 < u(t) < 1 and Also, for any positive time Assume that there exists a t * > t 1 such that where v (t * ) = 0 and u (t * ) < 0. If no such t * is reached, then we are done. If this time t * exists, then v will have a maximum at this point and so v(t) < 1 for all time where u(t) + K > 0 . Hence we have that u(t * ) and v(t * ) are both positive. By (2.2a), there does not exist a least time t 2 > t * where u(t 2 ) = 0, u (t 2 ) ≤ 0 and v(t 2 ) > 0. Moreover, by (2.2b), there does not exist a least time t 2 > t * where v(t 2 ) = 0, v (t 2 ) ≤ 0 and u(t 2 ) > 0. Hence, if either u or v were to become negative, then u, v would have to simultaneously become zero at t 2 . Hence, by (2.2a) and (2.2b), we would have that u(t 2 ) = u (t 2 ) = v(t 2 ) = v (t 2 ) = 0. Moreover, by differentiating (2.2a) and (2.2b), we would find that all derivatives of u (and v) were zero at this point t 2 . For a smooth function u ∈ C ∞ (0, T ], it cannot be that all derivatives of a non-trivial function at a certain time t 2 are zero. Hence, this point t 2 does not exist. It follows that u(t) > 0, v(t) > 0, ∀t ≥ 0. Finally, as u(t) + K − α > 0, . From the bounds in this Lemma, we can deduce that, for i = 1, 2, where g := max t∈[0,T ] |g(t)|. However, these bounds are not sufficiently sharp for our purposes, as they do not identify the fact that the large derivatives will only occur initially. To generate sharper bounds, we will construct a Shishkin decomposition [6] of the solution.
(2.5) The reduced solution (u 0 , v 0 ) or outer solution approximates the solution (u, v) outside a neighbourhood of t = 0. Using the stretched variable for t ≤ Cε. This motivates the following Shishkin decomposition of the solution (u, v).
Lemma 2. The solutions of problem (2.2) can be decomposed as follows: and the remainder terms R u , R v are bounded by By inserting this expansion for v into (2.2b), we see that the remainder term R v satisfies the initial value problem Hence, Note that g ≤ C and by the previous lemma, Using the differential equation in (2.7), we conclude that We also have the following decomposition Note that u 0 (t) is implicitly defined in (2.5). By inserting the expansions for u and v into (2.2a), we see that Hence, the remainder term R u satisfies the initial value problem Hence, as (u 0 + K)(u + K) ≥ K 2 > 0, by writing out a closed form representation for the function R u , we deduce that

III. DISCRETE PROBLEM
Consider the following implicit linear finite difference scheme for the continuous problem (2.2): and and the piecewise uniform mesh Ω N is constructed below. Note that for any s > 0, since |u 0 (t)| ≤ 1, Hence, the initial layer function B(t) is bounded by Based on the decomposition (2.6b) and (2.6a), we will use a piecewise uniform Shishkin mesh [6], denoted here by where there the fine and coarse mesh sizes k, K are defined by The transition point σ is taken to be For simplicity of exposition we assume 1 that ε(ln N ) 2 ≤ C then σ 2 ≤ Cε and for N sufficiently large, In other words, outside the initial layer, the solutions (u, v) are computationally close to the reduced solutions (u 0 , v 0 ) when ε is sufficiently small.
We next establish that, within the fine mesh, the discrete solutions U, V are bounded by the same bounds as their continuous counterparts u, v. In the next section, we will extend this result to the mesh points outside the fine mesh. These bounds on the discrete solutions ensure that the linear system (3.8) has a unique solution.
Lemma 4. For the solution of (3.8) and N sufficiently large (independent of ε), we have, for all Proof: By explicitly solving the linear system (3.8) at the first internal mesh point, we see that Define the associated system matrix and write the discrete problem (3.8) in the form We now complete the argument by using induction. Assume the statement is true for all From this, we deduce that We rewrite If (1 + K)V (t j−1 ) < 1 and j ≤ N/2 then k j ≤ CεN −1 ln N and, under these assumptions,

IV. ERROR ANALYSIS
Consider the error (U − u, V − v)(t j ) at each mesh point t j , which satisfies We rewrite these equations in matrix form where the matrix M j is defined in (3.9). The next theorem establishes an error estimate within the fine mesh region [0, σ].
Theorem 5. For the solutions of (2.2), (3.8) and N sufficiently large (independent of ε), we have, for all t j ≤ σ, Proof: For t j ≤ σ, Hence, As in the previous lemma, det(M j ) > 1 and using the explicit form of the inverse M −1 j in (3.10), we deduce that Hence, using an iterative argument, for all j ≤ N/2 It is not as straightforward to establish this error bound in the coarse mesh region. Before we examine this error, we present a discrete analogue to the solution decomposition established in Lemma 2. We first construct the decomposition in the fine mesh region, as it is relatively straightforward to do so. Then, in the next section, we inductively generate this same discrete decomposition and simultaneously establish error bounds at each point t j in the coarse mesh region.
We have the following decomposition of the discrete solution V : The remainder term R N v satisfies R N v (0) = 0 and the finite difference equation For t j ≤ σ and N sufficiently large, Then in the fine mesh U (t j−1 ) + K ε k ≥ 2 ln N N and from [6, Lemma 5.1] Hence, we have the bound where Note that (4.12) Hence, using a discrete comparison principle 2 Since U (t j ) > 0, by Lemma 3, we have that for all t j ≤ σ Theorem 6. For the solution of (3.8) and N sufficiently large (independent of ε), we have, for all t j ≥ σ, 0 < U (t j ) < 1, 0 < V (t j ) < 1 1 + K and the following parameter-uniform error bound Proof: The proof is by induction. By the previous two lemmas the statement is true for t k = σ.
Assume now that it is true for all σ ≤ t k ≤ t j−1 . The argument in Lemma 2 to establish the bounds 0 < U (t j ), 0 < V (t j ) < 1 1 + K is still valid within the coarse mesh. However, the argument to establish the upper bound U (t j ) < 1 requires an alternative argument in the coarse mesh region.
Under the induction assumption, for all σ ≤ t k ≤ t j the decomposition in (4.11a) Returning to the end of the proof of Lemma 2, we have that Using the decomposition (4.11a), we see that for N sufficiently large. Hence, as in the proof of Lemma 2, U (t j ) < 1.
For each t j , using the decompositions (2.6b) and (4.11a) we can write the error in v in the form For t j > σ, using (2.6b) we have that Hence, using (2.6b) again, we get that Also In the coarse mesh The remainder term R u can be further decomposed into the sum Hence, on the coarse mesh, Hence, the error (U − u)(t j ) satisfies By the induction hypothesis and using the expressions in (4.13) and (4.12), we deduce that Since the coefficient of the zero order term in this error equation is strictly bounded below by a positive number, we deduce that The bound on the error in V (t j ) follows from (4.15). The proof is completed by induction.
LetŪ andV denote the piecewise linear interpolants of U and V , respectively. Then, as in [6], we can extend the result to a global error bound.
Remark 8. According to [7, pg. 186], it is of interest to generate an accurate approximation to the rate of reaction u (t). Given the previous results, we have that In the fine mesh as u ≤ C. Hence, we have established the nodal error bound Moreover, as (2.2a) implies that this nodal error bound can be extended to the global error bound. For all 1 ≤ j ≤ N

V. NUMERICAL RESULTS
In this section we present some numerical results for the numerical method (3.8) applied to the system (2.2) where the parameters in (2.1) have been taken to be k 1 = 16847, k −1 = 7, k 2 = 12, This yields the following parameter values for the non-dimensional system (2.2) These values are used in [1] to fit the Henri-Michaelis-Menten system to experimental data derived from acetylcholine hydrolysis by acetylcholinesterase. We examine the performance of the numerical method over an extensive range of the parameter ε, which is equivalent to varying the initial concentration e(0) of the enzyme. In Figure 1 the computed approximations U, V (generated from the finite difference scheme (3.8)) are displayed. The plots confirm that v has an initial layer. Larger values of the parameter ε will result in less steep gradients appearing in the plot of v.
The global orders of convergence of the finite difference scheme (3.8) are estimated using the double-mesh principle [2,Chapter 8,pg. 170]. Note that this principle provides estimates of the orders of convergence despite the fact that the exact solution is unknown.
We denote by U N and U 2N the computed solutions on the Shishkin meshes Ω N and Ω 2N respectively. These solutions are used to compute the maximum two-mesh global differences whereŪ N (Ū 2N ) denotes the linear interpolant of the discrete solutions U N (U 2N ) over the mesh Ω N (Ω 2N ), respectively. The uniform global twomesh differencesD N and their corresponding uni-   Tables I  and II. The uniform two-mesh global differences D N , and their global orders of convergencep N are given in the last two rows of each table. These numerical results are in line with the asymptotic error bound established in Theorems 5 and 6. In the final Table III, the uniform two-mesh global differences and their global orders of convergence for the discrete approximations to the rate of reaction u (t) are displayed, which again show parameter-uniform convergence. In all three Tables we observe that, for N sufficiently large, the global orders of convergencep N are tending towards the rate associated with the bound N −1 (ln N ) 2 .

VI. CONCLUSION
A numerical method is constructed for a singularly perturbed system of two coupled nonlinear initial value equations. Theoretical error bounds are established at all time points, which guarantee that the numerical approximations converge to the continuous solution, irrespective of the size of the singular perturbation parameter. Numerical results support these theoretical error bounds.