An exact discretization of Michaelis–Menten type population equations

Abstract

Many single population models involve the inclusion of Michaelis–Menten type terms. We derive an exact discretization for such rational functions and demonstrate how these relations may be incorporated into more complex models representing the dynamics of single populations. The Lambert W-function plays an important role in this work.

Keywords:

• population models
• Michaelis–Menten dynamics
• numerical integration

AMS Subject Classification :

• 34A05
• 34A50
• 65L12
• 92D25

1. Introduction

The first-order differential equation

provides a good mathematical model for a broad range of phenomena in population dynamics, where the ‘populations’ may be humans, fish, chemicals, cells, etc. Terms having the form on the right-hand side of Equation (1), in particular appear in studies of chemostats, chemical reaction kinetics, morphogenesis, and continuous ventilation-volume models. References 1 3 5 9 13 14 give the relevant background materials on these topics as they relate to Equation (1) and its generalization to systems of equations describing interacting populations. For our purposes, we take the parameters (a, b) in Equation (1) to be non-negative, and select n to be either one or two. Generally, the n=1 case is called the Michaelis–Menten (M–M) Equation 9, i.e.

The goal of this study is to derive an exact discretization for the M–M equation and show how this result can be applied to the formulation of discretizations for other differential equations containing a term similar to that appearing on the right-hand side of Equation (2). Our discretization is constructed within the framework of the finite difference methodology 6 as generalized in the work of Mickens 10.

The importance of the results presented here is that the differential equations occurring in population dynamics models containing M–M type terms generally do not possess exact analytical solutions expressible in terms of a finite number of elementary functions 17. Therefore, numerical integration methods must be applied to obtain precise details on the solutions for given a priori values of the parameters 6 10. (Of course, certain properties of the solutions may be determined analytically by use of qualitative methods; these include: location of fixed-points, linear stability properties of fixed-points and special solutions, boundedness, etc. 7 17.) Of equal importance is the fact that knowledge of the exact discretization for the M–M equation can aid in the construction of improved finite difference schemes for other differential equations containing the same or similar terms modelling their dynamics 6 10.

At this point, it is worthwhile to allow a brief discussion on the use of standard (black-box) numerical integration packages (SNIP) or algorithms. SNIPs are based on the assumption that general, valid algorithms may be constructed for particular classes of ODEs. The user needs essentially no or little knowledge of the general properties of the equation's possible solution behaviours. However, there is clear evidence that the use of this approach can lead to incorrect results, i.e. numerical solutions, that are not consistent with the actual solutions of the original ODEs. Several examples of this situation may be found in the work of Dimitrov and Kojouharov 4 who investigate the use of finite different schemes to determine numerical solutions to general predator–prey population models. Their numerical simulations demonstrate that standard Rouge–Kutta numerical methods 6 can give erroneous results in comparison with a priori known mathematical results. However, when they apply ‘positive and elementary stable non-standard finite-difference (NSFD) methods’, the correct numerical behaviour is found.

Since population models involving M–M-type terms appear in a broad range of bioscience fields, it is of value to have a set of model equations for which the exact finite-difference schemes are known. In general, the use of these schemes will not allow the construction of exact schemes for an arbitrary set of differential equations containing M–M-type expressions; however, the use of the knowledge obtained from the previous model schemes may allow us to have increased confidence in our new discretizations as compared with the application of standard methods 6 10 17.

In the next section, we derive the exact solution for the M–M equation (2). This result is then used, in Section 3, to construct the corresponding exact finite difference representation. Finally, in the last section, we discuss our results and some related issues.

2. Exact solution

Examination of Equation (2) shows that it is a separable first-order differential equation. Its implicit solution is, for x(t ₀)=x ₀, given by the expression

However, an explicit solution may be written if use is made of the so-called Lambert W-function which is defined as the solution to the following transcendental equation 2 16:

Further, given

the solution in terms of the W-function is 2 16

Comparing Equations (3) and (5) gives

Note that D can be written as

Therefore, the solution to the M–M differential equation, expressible in terms of the Lambert W-function, is

Using the result, |ζ|≪1,

we can show

which is the correct limiting solution for Equation (2) when b=0 and the initial condition is x(t ₀)=x ₀. This answer is a check on the correctness of the general result given in Equation (9).

3. Exact finite difference scheme

It has been shown by Mickens 10 that if the solution to

is

then the exact finite difference scheme for Equation (11) is

where t _k =hk, with h=Δ t; and x _k =x(t _k ), and k=0, 1, 2, … (see Appendix 2). Therefore, the exact finite difference discretization to the M–M equation is, from Equation (9),

This relation can be rewritten to the form

where φ is the ‘denominator function’ and has the property 6

Note that at this stage of the derivation, any φ satisfying the condition in Equation (17) is valid. However, we select 11 12

since this choice is consistent with Equation (11).

In summary, the exact finite difference scheme or discretization, for the M–M equation (2), is given by Equation (16), where the denominator function φ (a, h) is taken to be the result expressed by Equation (18).

4. Discussion

It should be noted that a standard finite difference scheme for the M–M equation might take a form such as 6

On the other hand, simple NSFD representations include 10 11

where

and

where

An easy and direct calculation shows that if x ₀>0, then solutions to Equation (19) may become negative, while if x ₀>0, solutions to both Equations (20) and (22) remain non-negative, the same condition as satisfied by the solutions to the differential Equation (2).

A reason for rewriting Equation (15), to the expression in Equation (16), is to have it in a form such that a discrete derivative appears on its left side. Comparing Equations (2) and (16), we may identify the following functional discretizations:

While these structures are not readily amenable to guesswork, our calculations clearly demonstrate that they can be directly derived from the original M–M differential equation using the NSFD methodology 11 12. One lesson following from this work is that exact finite difference schemes, for even rather elementary non-linear differential equations, may have a very complex structure. However, this is precisely where the NSFD methodology comes to the forefront 10 11 12.

An application of our result for the M–M equation is to consider

where λ>0. This equation might represent the constant production of a hormone into the blood, followed by its removal by the body from the blood. Based on our results, presented in the previous section, the following finite difference discretization can be generated

where φ is given by Equation (18) and W(···) is the function on the right-hand side of Equation (15)

A related differential equation is

Note that all solutions decrease monotonically to zero. The exact solution for this initial condition is 10

and the corresponding exact finite difference discretization is

This is a remarkable result and is of some interest since there exist real world applications in which terms such as that on the right-hand side of Equation (28) appear. Particular examples include a model for spruce budworm outbreak 8

and a model for the calcium-stimulated-calcium-release mechanism 15

All parameters in Equations (31) and (32) are non-negative.

Our future work in this area will consist of investigating the numerics of discretizations of Equations (26), (31), and (32) based on the derivations presented in this paper. Further, we will study the influences on the discretizations coming from the addition of diffusion terms to these equations 8 14 15.

Finally, it should be realized that the general NSFD methodology can be easily understood and mastered with little effort if the potential user is willing to put in a small investment of their time. The major payoff will be the ability to generate numerical schemes that are dynamically consistent with the original differential equations 12. In particular, important features such as a condition of positivity for the populations and the satisfaction of existing conservation laws are easily incorporated into the discretizations. The standard schemes may not possess these properties (for arbitrary values of the step-size).

Acknowledgements

This research was supported, in part, by funds from the Clark Atlanta University, School of Arts and Sciences, Professional Development Funds. We also thank Professors Kale Oyedeji and Sandra Rucker for several valuable discussions on the topic of this paper.

Appendix 1. The Lambert W-function

The Lambert W-function is the inverse function of . Thus, it is defined as the solution to the transcendental equation 7

It satisfies the differential equation

and has the following Taylor series around z=0,

and this series has a radius of convergence of 1/e. The function defined by this series is an analytic function of z defined in the complex plane with a branch cut along the interval , and this analytic function defines the principal branch of the Lambert W-function.

Further details on this function, including its precise definition, mathematical properties, and various applications are given in Liu 7 and Corless et al. 2.

Appendix 2. Exact schemes

Consider the scalar ODE

where p denotes the parameters appearing in the equation, and f(x, t, p) is such that a unique solution exists over the interval, t ₀<t<T, for some T>t ₀, and p in the interval, . (For population models, in general, t ₀=0 and T=∞, i.e. the solution exists for all times.) The solution to Equation (A4) can be written as

with

Let a particular finite difference discretization of Equation (A4) be

where , , where h=Δ t, and k=0, 1, 2, 3˙s. Its solution can be expressed as

with

Definition A1

Equations (A4) and (A7) are said to have the same general solution, if and only if

where x(t _k ) and x _k are solutions, respectively, of Equations (A4) and (A7), i.e. Equations (A5) and (A8), for arbitrary values of h, h>0.

Definition A2

An exact difference scheme is one for which the solution to the difference equation has the same general solution as the associated differential equation.

Using these definitions, the following theorem holds 10 12:

Theorem A1

The ODE, given by Equation (A4), has the exact finite-difference scheme

where φ is that in Equation (A5).