key: cord-0007298-rdoi3psr
authors: Anguelov, R.; Dumont, Y.; Lubuma, J. M.-S.
title: On nonstandard finite difference schemes in biosciences
date: 2012-10-02
journal: AIP Conf Proc
DOI: 10.1063/1.4758961
sha: f9021262d6cab003df28424527f9e2f228de838c
doc_id: 7298
cord_uid: rdoi3psr

We design, analyze and implement nonstandard finite difference (NSFD) schemes for some differential models in biosciences. The NSFD schemes are reliable in three directions. They are topologically dynamically consistent for onedimensional models. They can replicate the global asymptotic stability of the disease-free equilibrium of the MSEIR model in epidemiology whenever the basic reproduction number is less than 1. They preserve the positivity and boundedness property of solutions of advection-reaction and reaction-diffusion equations.

The ever growing complexity of natural processes leads to challenging mathematics for the analysis of models in biosciences. Even in those cases when the mathematical theory is well-established, the differential models cannot be completely solved by analytical techniques. Consequently, numerical simulations are of fundamental importance in gaining some useful insights on the solutions of the differential equations.

In recent years, the authors made important contributions on the design, analysis and implementation of nonstandard finite difference (NSFD) schemes in biosciences. However, the results were written in diverse contexts and are to some extent scattered in the literature. The purpose of this paper is to provide a unified and focused presentation that could be user friendly. (For works devoted to the applications of NSFD schemes, we can mention the edited volumes [1, 2] and the references therein).

We report on reliable numerical schemes for bio-models in the following three directions/sections: One-dimensional models, whereby NSFD schemes that are topologically dynamically consistent are considered; Epidemiological models, whereby the focus is on NSFD schemes that replicate the global asymptotic stability of the disease-free equilibrium; Advection-reaction-diffusion models, whereby the interest is in NSFD schemes that preserve the positivity and the boundedness property of solutions. The last section is devoted to concluding remarks where potential research problems for future work are mentioned.

In biosciences, simple models, specifically one-dimensional ones, constitute a convenient setting for developing, demonstrating and testing new methods and approaches, which then can be extended to multidimensional systems. Here, we restrict ourselves to three important examples. We start with the decay equation

which arises in many applied settings such as population dynamics and compartmental models in epidemiology. The solution of the decay equation is u(t) = u 0 e −at (2) and, thus, the corresponding exact scheme is [3] :

Here and after, for k = 0, 1, 2, ···, u k denotes an approximation to the solution y(t) at the discrete time t = t k = kh where h = Δt is the step size. The second example is the logistic equation

which is abundantly used in biosciences. Its solution is

which yields the exact scheme [3] :

As a third example, we consider a less trivial one, namely the Michaelis-Menten (M-M) model:

Its exact scheme, constructed recently in terms of the Lambert W function, is [4] :

More generally, we consider a dynamical system on an interval D ⊂ R defined by a scalar differential equation

as well as a numerical scheme of the form

which is a discrete dynamical system on D. Obviously, the function G = G(h) = G(h; y) must satisfy the usual consistency conditions for the argument u:

It is also assumed that involved functions have the needed smoothness properties. The exact schemes of our three examples and in fact the exact schemes of the many models investigated in [3, 5] involve two characteristics that led two of the authors [6] to the following definitions: (10) that determines approximate solutions u k to the solution u(t) of a differential equation (9) is called a nonstandard finite difference (NSFD) scheme if at least one of the following two conditions is met: • The classical denominator h of the discrete derivative is replaced by a nonnegative function φ (h) such that

• Nonlinear terms that occur in the right-hand side f (u) of the differential equation (9) are approximated in a nonlocal way, i.e., by a suitable function of several points of the mesh.

To go further with the qualitative feature of NSFD schemes, we assume that the following condition holds for all arguments h and u:

We have the following result [7] . The first part is an important "practical" characterization of Condition (13) , while the second part is a substitute for elementary stability in the case of non-hyperbolic fixed points.

The condition (13) is necessary and sufficient for the scheme (10) to be qualitatively stable with respect to monotone dependence on the initial value, i.e.,

Furthermore, under the condition (13) , provided that the scheme (10) has no spurious fixed-points compared with the equilibrium points of the ODE (9), this scheme is qualitatively stable with respect to the property of monotonicity of solutions in the sense that the sequence (u k ) is increasing or decreasing when the exact solution u(t) is increasing or decreasing, respectively.

Theorem 2 makes the relevance of the condition (13) clear. This is further strengthen in the following precise manner [8] : (11) and (13) . Assume that for every h > 0 the fixed points of the map G(h) : D → D are the roots of f (u) = 0. Then the numerical scheme (10) is topologically dynamically consistent with the continuous dynamical systems (9) in the following sense: every map G(h) is topologically equivalent to each evolution operator S(t) of the dynamical system (9) .

In the light of Theorem 3, Condition (13) is the key requirement if numerical schemes are to replicate topological properties of the continuous dynamical systems. Coming back to our examples, we add to their exact schemes some NSFD schemes to which Theorem 3, including Condition (13), applies. For the logistic equation, we have the NSFD scheme [3] :

For the M-M model, we mention the following NSFD schemes [9, 10] :

Despite its apparent implicit nature, the NSFD scheme (17) is a discrete dynamical system on the interval D = [0, ∞) that admits the representation (10) and meets the requirement (13) , since this scheme is equivalent to a quadratic equation in u k+1 that admits only one nonnegative root whenever u k ≥ 0. In contrast, the forward Euler method fails to meet the requirement (13) when applied to the decay, the logistic and the M-M equations. The following figures speak for themselves as far as the performance of the NSFD and standard schemes is concerned: Figure 1 for the decay equation and Figure 2 for the logistic equation with a = 1 and h = 1.8; Figure 3 for the M-M equation where on each picture, the three NSFD schemes (16)-(18) for b = 1 and b = 0 (decay equation) are compared to the exact scheme (8). 

We consider a class of dynamical systems on Ω ⊂ R n + of the form where x = x(t) : [0, ∞) → R n represents the number/density of the populations, x 0 ∈ R n represents the initial populations, some coordinates of the vector f could be the recruitment/birth rates and A (x) is a nonlinear real n × n Metzler matrix.

We make two typical assumptions. These are: the population x is decomposable into non-infected y and infected z individuals, i.e., x = (y, z) ∈ R n 1 × R n 2 , with n 1 + n 2 = n; there exists a unique disease-free equilibrium (DFE) of the form x * = (y * , 0), y * = 0.

We are interested in the global asymptotical stability (GAS) of the DFE, a property that is highly desirable in order to design an intervention to reduce the spread of disease or to eradicate it. More precisely, we study numerical schemes for (19) that are reliable in that they replicate the GAS of the DFE. The schemes under consideration are discrete dynamical systems on Ω of the following form having the disease-free equilibrium x * as fixed-point:

Checking the GAS of an equilibrium/fixed-point is usually based on a Lyapunov function, which is unfortunately not easy to construct. Thus, we use an alternative approach that is based on the availability of the following representations of the systems (19) and (20):

Here A 1 (x), A 2 (x) and A 12 (x) are n 1 × n 1 , n 2 × n 2 and n 1 × n 2 matrices that need not be the same for the discrete and the continuous systems.

For the continuous system, the result due to [11] reads as follows:

Theorem 4 Assume that H 1 . The dynamical system (21) or (19) is dissipative on Ω; H 2 . The equilibrium point y * of the sub-systemẏ

is GAS on the canonical projection of Ω on R n 1 + ; H 3 . The matrix A 2 (x) is a Metzler matrix and is irreducible for each x ∈ Ω; H 4 . There exists an upper bound matrixĀ 2 for the set M = {A 2 (x) : x ∈ Ω} with the following property: either A 2 / ∈ M , orĀ 2 ∈ M and for eachx ∈ Ω satisfyingĀ 2 = A 2 (x) , we necessarily havex ∈ R n 1 + × {0}; H 5 . α Ā 2 ≤ 0, where α Ā 2 is the maximum of the real parts of eigenvalues of the matrixĀ 2 .

Then, the disease-free equilibrium x * is GAS in Ω.

The discrete counterpart of Theorem 4 is established in [12] :

Theorem 5 Assume that D 1 . The discrete dynamical system (22) or (20) is dissipative on Ω; D 2 . The fixed-point y * of the sub-system

is GAS on the set Ω 1 = {y ∈ R n 1 : (y, 0) ∈ Ω}; D 3 . The matrix A 2 (y, z) is nonnegative for all (y, z) ∈ Ω; D 4 . There exists an irreducible matrix A 2 , which is an upper bound for the set M = {A 2 (y, z) : (y, z) ∈ Ω}. D 5 . The matrix A 2 is such that its spectral radius satisfies ρ(A 2 ) ≤ 1. In the case when ρ(A 2 ) = 1, we assume that A 2 (y, z) >> 0, ∀(y, z) ∈ Ω and, in addition, z reduces to the zero vector in R n 2 if there exists a vector (y, z) ∈ Ω such that A 2 (y, z) = A 2 .

Then, the fixed-point (y * , 0) is GAS for discrete system on Ω.

Theorems 4 and 5 apply to a variety of epidemiological models (see for example [11, 13, 14, 15] ). To be more specific, we consider the MSEIR epidemiological model that involves a total population of size N = N(t) divided in five compartments: M = mN -infants with passive immunity, S = sN -susceptible, E = eN -exposed, I = iNinfective, and R = rN -recovered individuals. Figure 4 shows the flow and the corresponding constant rates between compartments, with b and d being the recruitment and death rates, respectively. 

where x = (m, s, e, i, r) T and the compartmental matrix is given by

The disease-free equilibrium is DFE ≡ x * = (0, 1, 0, 0, 0) T . On setting

x = (y, z) T , y = (m, s, r) T , z = (e, i) T , y * = (0, 1, 0) T ,

and

it follows that the model (24) admits the matrix representation (21) . More so, the requirements (H1)-(H4) in Theorem 4 are met, the requirement (H5) being equivalent to (26) below. Thus, we rediscover the following known result [16] :

Theorem 6 For MSEIR model, a dissipative dynamical system on the biologically feasible region Ω in (25) , the DFE is GAS whenever the basic reproduction number is less than 1:

Given θ ∈ [0, 1] andθ ∈ [0, 1], we consider the NSFD scheme:

The denominator function φ (h), chosen in accordance with (12) , must satisfy some additional conditions in order to guarantee the following important facts :

• The NSFD scheme amounts to solving at each step a linear system

whereC(i k , s k ) is an M-matrix,D ≥ 0 and thus the scheme is a dissipative discrete dynamical system on Ω;

• The NSFD scheme is second-order convergent when θ =θ = 1/2;

• The NSFD scheme satisfies the conditions in Theorem 5, specifically (D3), (D4) and (D5).

Typically, as illustrated in the numerical examples below, we take

where the number Q captures the feature of the model through the condition

In summary, for such a suitable φ (h), we obtain the following result [12] :

For θ andθ such that either θ = 0 or θ +θ ≥ 1, the fixed-point x * = (0, 1, 0, 0, 0) of the NSFD scheme, on the set Ω, is GAS whenever R 0 < 1. Furthermore, the NSFD scheme is second-order convergent if θ =θ =1/2 and first-order convergent otherwise.

Remark 8 As mentioned earlier, the authors have applied Theorems 4 and 5 to other compartmental models (see for instance [13, 14, 15] ). The more critical part is to check the conditions (H5) and (D5). In general, this amounts to checking that the basic reproduction number satisfies the condition R 0 ≤ 1, as stated in (26) . However, an additional threshold number ξ < 1 may arise so that the GAS of the DFE is rather subjected to a more restrictive condition R 0 ≤ ξ . For ξ < R 0 ≤ 1, we can have the backward bifurcation phenomenon whereby there is co-existence of the stable DFE with a stable endemic equilibrium. We investigated the backward bifurcation in [13] for a malaria model and in [17] for a disease, like the HIV/AIDS, with multiple strains.

In the numerical simulations below, where h = 2 or h = 10, we use the following values: β = 0.14, δ = 1/180, ε = 1/14, γ = 1/7, d = 1/(40 * 365), q = 0 and thus R 0 = 0.97857. While the theory on the GAS of the DFE is perfectly supported by Figure 5 on the NSFD scheme, Figure 6 displays so well the poor performance of classical methods to the extent that solutions can be negative. For the second-order convergence, the result on Figure 5 (d) are further tabulated in [12] . Note that in each figure, we plot the fractions of populations against days as time unit.

The advection, reaction and diffusion processes are naturally involved in differential models in biosciences. We start with the two-dimensional advection-reaction equation

where the reaction term r(u), which is assumed to have a nonzero finite number of roots, and the initial function f (x, y) are as smooth as needed, a and b being positive constants. Here we assume the following regarding the solution u and the initial condition f :

We are then interested in constructing NSFD schemes that preserve the property (31). By the method of characteristics, the solution of (30) is:

This leads to the exact scheme [3] 

Here, u k m,n denotes an approximation of u(x, y,t) at the grid point t = kΔt, x = mΔx and y = nΔy. Thus for a general reaction term, we can consider the NSFD scheme 

Similarly to (28) and (29), the feature of the model is captured by a number Q ≥ max{|r ( u)|; r( u) = 0} as well as by a function φ like in (12) that also satisfies 0 < φ (z) < 1 for z > 0.

It is proved in [18] , that apart from being elementary stable, the NSFD scheme (35) preserves the property (31) whenever (34) holds and −s ≤ r(s)/Q ≤ 1 − s for 0 ≤ s ≤ 1. This result is illustrated in Figure 7 c. Next, we consider the one-dimensional reaction diffusion equation

for which we make the assumption (31). NSFD schemes that preserve the property (31) for the general equation (37) are discussed in [19] . In the particular case of the Fisher equation

which satisfies the property (31) (see [20] ), a NSFD scheme of the form

is proposed in [19] with various expressions of the denominator function φ (Δt). Further extension to the SIS-diffusion model is done in [21] . It is proved in [19] that the NSFD scheme (39) replicates the property (31) whenever we have the following functional relation between the step sizes:

The result is illustrated on Figure 8 for the initial function f (x) = 0.5 + 0.5 sin 2x, λ = 25 and Δx = .25 for various φ (Δt). The pictures confirm the superiority of the NSFD schemes over the standard explicit finite difference scheme, which further displays unwanted oscillations. NSFD schemes and other methods for the full advection-reaction-diffusion equation as well as their impact in biosciences are investigated in [9, 10, 22, 23] . 

We gave an overview of our recent work on NSFD schemes in biosciences, with the emphasis on three aspects of the reliability and excellence performance of these schemes. Each aspect of the study leads on its own to research problems on which we are currently working. These include: the extension of the concept of topologically dynamically consistent NSFD schemes to multi-dimensional models [8] , the design of positivity-preserving NSFD schemes for cross diffusion equations (or other models in ecology) [24, 25, 26] and the construction of NSFD schemes that replicate the backward bifurcation phenomenon for Volterra integral equation epidemiological models [21, 27] .

Applications of Nonstandard Finite Difference Schemes

Advances in the Applications of Nonstandard Finite Difference Schemes

Nonstandard Finite Difference Models of Differential Equations

Contributions to the theory of non-standard finite difference methods and applications to singular perturbation problems

Nonstandard finite difference schemes for Michaelis-Menten type reaction-diffusion equations

Dynamically consistent nonstandard finite difference schemes for epidemiological models

Stability analysis and dynamics preserving nonstandard finite difference schemes for a malaria model

An Introduction to Nonlinear Partial Differential Equations

Analysis and dynamically consistent numerical schemes for the SIS model and related reaction diffusion equation

Spatio-temporal modeling of mosquito distribution

Mathematical Biology II: Spatial Models and Biomedical Applications

Mathematical studies on the Sterile Insect Technique for the Chikungunya disease and Aedes albopictus

The authors would like to thank the late Prof. C.I. Christov as well as Prof. M.D. Todorov for the invitation to present this work at the AMiTaNS'12 conference in Varna, Bulgaria. We thank the anonymous referee for comments, which have improved the readability of the paper and the consistency in notation.The first and third authors were supported by the National Research Foundation of South Africa. The second author acknowledges the support of the TIS Project, funded by the French Ministry of Health and the European Union.