Nonstandard finite difference models of some differential equations - - PowerPoint PPT Presentation

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Nonstandard finite difference models of some differential equations in Life Sciences

Jean M-S Lubuma

Department of Mathematics and Applied Mathematics Faculty of Natural and Agricultural Sciences University of Pretoria (South Africa)

10th Blackwell-Tapia Conference Institute for Computational and Experimental Research in Mathematics (ICERM), Brown University, Providence, USA, 9-10 November 2018,

Lubuma (UP) 2018 1 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

My first contact with the NSFD method

Lubuma (UP) 2018 2 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Motivation and outline

1

Complexity of natural processes leading to challenging mathematical models;

2

Differential Equations cannot be completely solved by analytical techniques;

3

Fundamental importance of reliable numerical methods to gain insight

4

Aim: Overview of some NSFD schemes that are dynamically consistent :

One-dimensional first-order ODEs models: NSFD schemes are topologically dynamically consistent; One-dimensional second-order ODEs models: NSFD schemes preserve the principle of conservation of energy Ebola Virus Disease models: NSFD schemes replicate their dynamics; Reaction-diffusion & Cross-diffusion models: NSFD schemes preserve boundedness and positivity. Colony Collapse Desorder: NSFD schemes preserve the (fast) declines of honeybee colonies. More on NSFD schemes

5

NSFD method in Africa

Lubuma (UP) 2018 3 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

One-dimensional 1st-order examples: Decay and Logistic

Exact scheme (Mickens, 1994) of the decay equation dy dt = −ay : (1) yk+1 − yk (1 − e−ah)/a = −ayk or yk+1 − yk (eah − 1)/a = −ayk+1 (2) Exact scheme (Mickens, 1994) of the logistic equation dy dt = y(1 − y) : (3) yk+1 − yk eh − 1 = yk(1 − yk+1) or yk+1 − yk 1 − e−h = yk+1(1 − yk) (4) Nonstandard scheme: yk+1 − yk h = yk(1 − yk+1). (5)

Lubuma (UP) 2018 4 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

One-dimensional 1st-order examples: M-M equation

Exact scheme (Mickens, 2011) of the Michaelis-Menten (1913) (M-M) model du dt = − au 1 + bu : (6) uk+1 − uk (1 − e−ah)/a = W [ (b exp(−ah)uk exp(buk) ] − buk b(1 − e−ah)/a . (7) NSFD schemes include (Chapwanya, JL and Mickens, 2012): uk+1 − uk (eah − 1)/a = − auk+1 1 + buk ; (8) uk+1 − uk (eah − 1)/a = − auk+1 1 + buk+1 ; (9) uk+1 − uk (1 − e−ah)/a = − auk 1 + buk . (10) Exact schemes of various equations: see, for instance, Roeger (2014)

Lubuma (UP) 2018 5 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

One-dimensional 1st-order examples: dynamic consistency

Each one of the previous NSFD scheme can be written as uk+1 = G(h; uk) where the function G satisfies the conditions G(0; u) = u and dG(h; u) du > 0. (11) This guarantees that the schemes: Are monotone dependent on initial values (Anguelov, JL, 2003): u0 ≤ v 0 = ⇒ uk ≤ v k Preserve monotonicity: u(t) increasing/decreasing = ⇒ uk increasing/decreasing Are topologically dynamically consistent, thus elementary stable (Anguelov, JL and Shillor, 2011): Every map G(h) is topologically equivalent to each evolution operator S(t) i.e. µ ◦ S(t) = G(h) ◦ µ for some homeomorphism µ.

Lubuma (UP) 2018 6 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

One-dimensional 1st-order examples: Simulations D

(a)

1 2 3 4 5 6 7 8 −1 −0.5 0.5 1 t y

(b)

1 2 3 4 5 6 7 8 −1 −0.5 0.5 1 t y

(c)

1 2 3 4 5 6 7 8 −1 −0.5 0.5 1 t y

(d)

2 4 6 8 10 −2.5 −2 −1.5 −1 −0.5 0.5 1 1.5 2 2.5 t y

Figure: Decay equation : (a) Exact solution; (b) Forward Euler method; (c) NSFD with complex denominator; (d) Runge-Kutta method

Lubuma (UP) 2018 7 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

One-dimensional 1st-order examples: Simulations L

(a)

1 2 3 4 5 6 7 8 −0.5 0.5 1 1.5

t y

(b)

1 2 3 4 5 6 7 8 −0.5 0.5 1 1.5

t y

(c)

1 2 3 4 5 6 7 8 −0.5 0.5 1 1.5

t y

Figure: Logistic equation: (a) Exact solution; (b) Forward Euler method; (c) NSFD

Lubuma (UP) 2018 8 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

One-dimensional 1st-order examples: Simulations M-M

10 20 30 40 50 60 0.2 0.4 0.6 0.8 1 t u NSFD scheme, b = 1 NSFD scheme, b = 0 Exact scheme, b = 1 10 20 30 40 50 60 0.2 0.4 0.6 0.8 1 t u NSFD scheme, b = 1 NSFD scheme, b = 0 Exact scheme, b = 1 10 20 30 40 50 60 0.2 0.4 0.6 0.8 1 t u NSFD scheme, b = 1 NSFD scheme, b = 0 Exact scheme, b = 1

Figure: M-M model: NSFD schemes (8)-(10) versus exact scheme (7). Thus M-M equation cannot be approximated by the decay equation: Michaelis-Menten (1913)

Lubuma (UP) 2018 9 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

NSFD Scheme: Quid?

Mickens’ rules: Rule 1. The orders of the discrete derivatives should be equal to the orders

• f the corresponding derivatives of the differential equation.

Rule 2. Denominator functions for the discrete derivatives must, in gen- eral, be expressed in terms of more complicated functions of the step-sizes than those conventionally used. Rule 3. Nonlinear terms should, in general, be replaced by nonlocal dis- crete representations. Rule 4. Special conditions that hold for the solutions of the differential equations should also hold for the solutions of the difference scheme. Rule 5. The scheme should not introduce ghost or spurious solutions. (e.g. the second order Runge-Kutta method, 2(yk+1 − yk)/h = f (yk) + f [yk + hf (yk)], for the decay equation, has spurious fixed-points for h = 2/a)

Lubuma (UP) 2018 10 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

NSFD Scheme: Quid?

Definition (Anguelov and JL, 2001)

1

A difference equation to determine approximate solutions yk to the solution y(t) of a differential equation is called a nonstandard finite difference 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 φ(h) = h + O(h2) as h → 0+. Nonlinear terms that occur in the right-hand side of the ODE are approximated in a nonlocal way, i.e., by a suitable function of several points of the mesh.

2

A NSFD scheme is called qualitatively stable or dynamically consistent with respect to some property P of the differential equation, whenever the discrete equation replicates the property P for every value of h. Remark 1: Our formal definition retains only Mickens’ Rules 2 and 3; Remark 2: Sufficent conditions for topological dynamic consistency in scalar case: uk+1 = G(h; uk), G(0; u) = u and dG(h; u) du > 0. (12)

Lubuma (UP) 2018 11 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

One-dimensional 2nd-order models

Exact scheme of the harmonic oscillator (Mickens, 1994) d2y dt2 + ω2y = 0 : yk+1 − 2yk + yk−1 (4/ω2 sin2(hω/2) + ω2yk = 0,

• r (conservation of energy: Anguelov, Kama, JL; 2005):

( yk+1 − yk 2/ω sin(hω/2) )2 + ω2yk+1yk = ( yk − yk−1 2/ω sin(hω/2) )2 + ω2ykyk−1. NSFD scheme for the Karman-Guderley-type equation (Buckmire, 1994-2004) 1 r d dr ( r dy dr ) − ω2y = 0 : 1 rk ( yk+1 − yk log(rk+1/rk) − yk − yk−1 log(rk/rk−1) ) − ω2yk = 0.

Lubuma (UP) 2018 12 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Ebola Virus Disease: fast and slow transmissions

JL, Tsanou et al. (2016): Can the consumption of contaminated bush meat, the funeral practices and the environmental contamination explain the recurrence and persistence of EVD outbreaks in Africa? S" I" R" P" D"

π"

σ" µ" b" infec*on" γ" µ+δ" """"""α" (shedding)" ξ""(shedding)" µ" η"

Figure: (1) Transmission of deceased individuals and demographic process (Agusto, Gumel et al. (2015), Fasina et al. (2014), Ivorra et al. (2015)); (2) Infection through contaminated environment; (3) Provision of Ebola virus from consumption of bats, hunted meat, etc.; (4) Not included: compartments of vaccinated and trained individuals

Lubuma (UP) 2018 13 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

EVD : Dynamical system

Biologically meaningful system in the region Kσ ⊂ R5

+ of points

(S(t), I(t), R(t), D(t), P(t)) such that H(t) ≤ Π µ , D(t) ≤ (µ + δ)Π bµ , P(t) ≤ b(σµ + ξΠ) + α(µ + δ)Π bηµ :                                          dS(t) dt = Π − (β1I + β2D + β3P) S − µS dI(t) dt = (β1I + β2D + β3P) S − (µ + δ + γ)I dR(t) dt = γI − µR dD(t) dt = (µ + δ)I − bD dP(t) dt = σ + ξI + αD − ηP. (13)

Lubuma (UP) 2018 14 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Symbols Biological descriptions S Susceptible human individuals I Infectious human individuals D Ebola infected and deceased human individuals R Recovered human individuals P Ebola virus pathogens in the environment Π Recruitment rate of susceptible human individuals η Decay rate of Ebola virus in the environment ξ Shedding rate of infectious human individuals α Shedding rate of deceased human individuals δ Disease-induced death rate of human individuals β1 Effective contact rate of infectious human individuals β2 Effective contact rate of deceased human individuals β3 Effective contact rate of Ebola virus γ Recovered rate of human individuals µ Natural death rate of human individuals 1/b Mean caring duration of deceased human individuals σ Recruitment rate of Ebola virus in the environment

Table: Variables and parameters for model (13)

Lubuma (UP) 2018 15 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

EVD: Virus-free versus virus environment σ ≥ 0

Virus-free environment σ = 0: DFE and basic reproduction number R0 : E0 = (S, I, R, D, P) = (Π/µ, 0, 0, 0, 0) and R0 = Πβ1 µ(µ + δ + γ) + (µ + δ)Πβ2 bµ(µ + δ + γ) + Π(bξ + αδ + αµ)β3 bηµ(µ + δ + γ) . (14) Virus environment σ > 0: No DFE. There is a unique EE denoted by E # = ( S#, I #, R#, D#, P#) , with I # ≡ I #(σ) the unique positive root of A2(I #)2 − A1I # − A0 = 0, (15) where A2 = (µ + δ + γ) [η (bβ1 + β2(µ + δ)) + (bξ + αδ + αµ) β3] , A1 = A1(σ) = bηµ(µ + δ + γ) [ R0 − 1 − σβ3 ηµ ] , A0 = A0(σ) = bπσβ3. (16)

Lubuma (UP) 2018 16 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

EVD: Main Result

Theorem Ebola virus-free environment (σ = 0):

1

The DFE E0 is GAS when R0 ≤ 1 and unstable otherwise.

2

If R0 > 1, there is a LAS EE E ∗ = (S∗, I ∗, R∗, D∗, P∗). It is GAS in the absence of shedding (α = 0) or manipulation of deceased human individuals before burial (ξ = 0). Ebola virus-environment (σ > 0):

1

The EE E # is LAS. It is GAS when α = 0 or when ξ = 0.

2

The infectious component I # = I #(σ) of the EE is an increasing function on the interval 0 ≤ σ < ∞, with I #(0) = I ∗ denoting the infectious component

• f the unique EE of the model when σ = 0.

N.B. The severity of the disease decreases with the reduction of the consumption

• f contaminated bush meat: I ∗ < I #(σ).

Lubuma (UP) 2018 17 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

EVD: NSFD Scheme

                                           Sn+1 − Sn ϕ = Π − (β1In + β2Dn + β3Pn) Sn+1 − µSn+1, In+1 − In ϕ = (β1In + β2Dn + β3Pn) Sn+1 − (µ + δ + γ)In+1, Rn+1 − Rn ϕ = γIn+1 − µRn+1, Dn+1 − Dn ϕ = (µ + δ)In+1 − bDn+1, Pn+1 − Pn ϕ = σ + ξIn + αDn − ηPn+1, (17) ϕ = ϕ(h) = ( 1 − e−(µ+δ+γ)h) /(µ + δ + γ) = h + O(h2). (18)

Lubuma (UP) 2018 18 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Dynamic consitency of the NSFD scheme from its explicit form:                                                              Sn+1 = [Πϕ + Sn]/[1 + ϕ(µ + Bn)], In+1 = [1 + ϕ(µ + Bn)]In + ϕ(πϕ + Sn)Bn [1 + ϕ(µ + δ + γ)] [1 + ϕ(µ + Bn)] Rn+1 = [1 + ϕ(µ + δ + γ)] [1 + ϕ(µ + Bn)] Rn (1 + µϕ) [1 + ϕ(µ + δ + γ)] [1 + ϕ(µ + Bn)] + γϕ [(1 + ϕ(µ + Bn))In + ϕ(πϕ + Sn)Bn] (1 + µϕ) [1 + ϕ(µ + δ + γ)] [1 + ϕ(µ + Bn)] Dn+1 = [1 + ϕ(µ + δ + γ)] [1 + ϕ(µ + Bn)] Dn (1 + bϕ) [1 + ϕ(µ + δ + γ)] [1 + ϕ(µ + Bn)] +(µ + δ)ϕ [(1 + ϕ(µ + Bn))In + ϕ(πϕ + Sn)Bn] (1 + bϕ) [1 + ϕ(µ + δ + γ)] [1 + ϕ(µ + Bn)] Pn+1 = [ϕ (σ + ξIn + αDn) + Pn]/(1 + ηϕ), (19) where Bn := β1In + β2Dn + β3Pn.

Lubuma (UP) 2018 19 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

EVD: Simulations - Ebola virus-free environment

200 400 600 −200 200 400 600 800 Temps(days) Populations

Suscept Infected Dead Virus

(a)

200 400 600 5 10 15 20 25 Time(days) Populations

Suscept Infected Dead Virus

(b)

Figure: (a) - Example of spurious negative solutions obtained with the ode45 implemented in MatLab. (b) - Positive solutions obtained with the NSFD scheme (17) implemented in MatLab. Parameters: µ = 0.5; δ = 0.05; γ = 0.06; Π = 10; b = 0.8; ξ = 0.04; α = 0.04; η = 0.03; β1 = 0.006; β2 = 0.012; β3 = 0.01 ; σ = 0 : R0 = 0.0602.

Lubuma (UP) 2018 20 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

EVD: Simulations - Ebola virus-free environment

100 200 300 400 500 5 10 15 Time(days) Populations

Suscept Infected Dead Virus

(a)

200 400 600 10 20 30 40 50 60 Time(days) Populations

Suscept Infected Dead Virus

(b)

Figure: (a) - GAS of the disease free E0. (b) - GAS of the endemic equilibrium E ∗, when σ = 0 ( R0 = 24.7 > 1) . For (a), the parameters are: µ = 0.03; δ = 0.5; γ = 0.006; Π = 10; b = 0.8; ξ = 0.0004; α = 0.004; η = 0.03; β1 = 0.006; β2 = 0.012; β3 = 0.0001. For (b), the parameters are: µ = 0.02; δ = 0.9; γ = 0.06; π = 10; b = 0.8; ξ = 0.04; α = 0.04; η = 0.03; β1 = 0.006; β2 = 0.012; β3 = 0.01;.

Lubuma (UP) 2018 21 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

EVD: Ebola virus-free versus virus environment

200 400 600 10 20 30 40 50 Time(days) Populations

Suscept Infected Dead Virus

(a)

200 400 600 2 4 6 8 Time(days) Infected σ = 0 σ = 0.6

(b)

Figure: (a) - Stability of E #, when σ = 0.6 ; (b) - Comparison of I ∗ and I #. µ = 0.02; δ = 0.9; γ = 0.06; Π = 10; b = 0.8; ξ = 0.04; α = 0.04; η = 0.03; β1 = 0.006; β2 = 0.012; β3 = 0.01. Moreover, for the blue dotted curve, σ = 0 and R0 = 24.7 > 1, whereas for the red curve, σ = 0.6.

Lubuma (UP) 2018 22 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Cross-diffusion model: malignant invasion, tumor growth

Avascular and Vascular Phases (Pepper & Lolas 2008):

1

Early avascular phase: solid tumors confined to tissues where they arise;

2

Avascular tumor needs to grow: supply of nutrients from the blood;

3

Tumor secretes Vascular Endothelial Growth Factors (VEGF), which induce

4

Angiogenesis: Formation of new blood vessels from preexisting vasculature;

5

Cells of a vascularized tumor can:

Invade and destroy the surrounding tissues, via MDEs, eg. Metallo proteases; Establish new colonies (metastasis). Figure: Angiogenesis

Lubuma (UP) 2018 23 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Cross-diffusion model: malignant invasion cells

Evolution in time t and space x of concentration of three physical variables: u, invasive cells; c, connective tissues; p, proteases. Model, with cross-diffusion in first equation (Marchant, Norbury and Perumpanani 2000): ut = u(1 − u) − (ucx)x (20) ct = −pc (21) pt = ϵ−1(uc − p) : (22) Some facts:

1

Positivity of solutions: for u0(x) ≥ 0, c0(x) ≥ 0 and p0(x) ≥ 0, we have u(x, t) ≥ 0, c(x, t) ≥ 0, p(x, t) ≥ 0 with c decreasing in time.

2

The outstanding problem solved recently (Chapwanya, JL, Mickens (2014)) is to obtain positive discrete solutions.

Lubuma (UP) 2018 24 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Malignant invasion model: classical NSFD Scheme

The space-independent case of Eq. (20) of the invasive cells is the logistic equation with exact scheme (Mickens 1994) uk+1 − uk e∆t − 1 = uk(1 − uk+1). Notation: ϕ(∆t) = e∆t − 1 and let 0 ≤ ψ(∆x) = ∆x + O[(∆x)2]. NSFD schemes by Mickens’ rules for Eq (21) of connective tissue and Eq (22) of proteases : ck+1

m

− ck

m

ϕ(∆t) = −pk

mck+1 m

• r ck+1

m

= ck

m

1 + ϕ(∆t)pk

m

; (23) pk+1

m

− pk

m

ϵϕ(ϵ−1∆t) = ϵ−1(uk

mck+1 m

− pk+1

m

) or pk+1

m

= pk

m + ϕ(ϵ−1∆t)uk mck+1 m

1 + ϕ(ϵ−1∆t) . (24)

Lubuma (UP) 2018 25 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Malignant invasion model: NSFD Scheme

Use 2uk+1

m

/(uk+1

m

+ uk

m) ≈ 1. Thus the scheme

uk+1

m

− uk

m

ϕ(∆t) = uk

m(1−uk+1 m

)+ (uk

m + uk m−1)ck m

ψ(∆x)2 − uk

mck m+1 + uk m−1ck m−1

ψ(∆x)2 × 2uk+1

m

uk+1

m

+ uk

m

, (25) which corresponds to the following quadratic equation in uk+1

m

: Ak

m(uk+1 m

)2 + Bk

muk+1 m

+ Dk

m = 0

where Ak

m > 0, Dk m ≤ 0, R = ϕ(∆t)/ψ(∆x)2 and

Bk

m

= −ϕ(∆t)uk

m − R(uk m + uk m−1)ck m + 2R(uk mck m+1 + uk m−1ck m−1). (26)

The only nonnegative root of the quadratic equation is uk+1

m

= −Bk

m +

√ (Bk

m)2 − 4Ak mDk m

2Ak

m

. (27) Theorem The NSFD scheme (27) preserves unconditionally the positivity of solutions.

Lubuma (UP) 2018 26 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Malignant invasion model: NSFD simulations

10 20 30 40 50 0.2 0.4 0.6 0.8 1 t Concentration u(x,t) c(x,t) p(x,t)

10 20 30 40 50 5 10 15 0.5 1 t x u(x,t) 10 20 30 40 50 5 10 15 0.5 1 t x c(x,t) 10 20 30 40 50 5 10 15 0.5 1 t x p(x,t)

Figure: NSFD scheme with ψ(∆x) = ∆x = 0.5; Top: all profiles at x = 10 (left) and profile u(x, t) (right); Bottom: profile c(x, t) (left) and profile p(x, t) (right)

Lubuma (UP) 2018 27 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Malignant invasion model: Standard FD scheme

10 20 30 40 50 0.2 0.4 0.6 0.8 1 t Concentration u(x,t) c(x,t) p(x,t) 10 20 30 40 50 10 20 0.5 1 t x u(x,t) 10 20 30 40 50 10 20 0.5 1 t x c(x,t) 10 20 30 40 50 10 20 0.5 1 t x p(x,t)

Figure: Standard finite difference analogue with ∆x = .5. Top: all profiles at x = 10 (left) and profile u(x, t) (right); Bottom: profile c(x, t) (left) and profile p(x, t) (right)

Lubuma (UP) 2018 28 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Declines of honeybee colonies: CCD

Figure: Khoury, Myerscough, Barron (2011)

dH dt = L(H + F) (H + F + ω) − ( α − σ F H + F ) H; dF dt = ( α − σ F H + F ) H − mF (28) where the eclosion and recruitment functions are: E(H, F) = L(H + F) (H + F + ω) and R(H, F) = ( α − σ F H + F ) .

Lubuma (UP) 2018 29 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Social parasitic model: CC

Figure: Social parasites: Lerata, JL, Yusuf (2018)

dH dt = L(H + F) (H + F + ω) − (α + ˜ m)H; dF dt = αH − mF. (29) N.B.: The number of clone workers who leave the H class is larger than the Savannah workers transiting to F class.

Lubuma (UP) 2018 30 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

NSFD Schemes

KMB NSFD Scheme:      Hn+1 − Hn ϕ = L(Hn + Fn) Hn + Fn + ω − αHn+1 + σ FnHn+1 Fn + Hn+1 Fn+1 − Fn ϕ = αHn+1 − σ FnHn+1 Fn + Hn+1 − mFn+1, (30) where ϕ(∆t) = 1 − exp(−m − σ − α)∆t m + σ + α . SP- NSFD Scheme:      Hn+1 − Hn ϕ = L(Hn + Fn) Hn + Fn + ω − (α + ˜ m)Hn+1 Fn+1 − Fn ϕ = αHn+1 − mFn+1, (31) where ϕ = [1 − exp(−(α + m + ˜ m)∆t)]/(α + m + ˜ m).

Lubuma (UP) 2018 31 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

NSFD Scheme: GAS of trivial/interior equilibrium

100 200 300 400 2000 4000 6000 8000 10000 12000 14000 time (days) Number of bees in colony 100 200 300 400 2000 4000 6000 8000 10000 12000 14000 time (days) Number of bees in colony 100 200 300 400 5000 6000 7000 8000 9000 10000 11000 12000 13000 14000 time (days) Number of bees in colony

Figure: ∆t = 2; KMB: No CCD (m = 0.24, bottom left ); CCD (m = 0.4, top left); SP: CC (faster CCD)(m = 0.24 = ˜ m, top right)

Lubuma (UP) 2018 32 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

More on NSFD Schemes

1

How to construct NSFD schemes?

Mickens and Washington: Denominator function from conservation laws; Anguelov and JL: Nonlocal approximation by perturbation techniques; K Patidar (2016): Review paper, J Difference Eq Appl

2

Compartmental models:

R Mickens (Ed; 2000; 2005), Applications of NSFD schemes; A Gumel (Ed; 2003; 2014), J Difference Eq Appl; Contemporary Maths S Bowong, Y. Dumont, etc. Water and vector borne diseases

3

Models exhibiting backward bifurcation:

Malaria: Anguelov, Dumont, JL & Mureithi 2013; Dukuza 2018 Model for disease with multiple strains: Garba, Gumel and JL 2011.

4

Advection-Reaction-diffusion: Anguelov, Kama, Kojouharov, JL, etc.

5

Higher-order NSFD schemes: Kama, JL, Mureithi, Roux, Terefe; 2003, 2015.

6

Volterra integral eqs, Delay eqs, impact problems, singularly perturbed problems, etc.: Gumel, Hassan, Dumont, JL, Patidar,Terefe, etc.

7

Challenge: NSFD schemes for Hopf bifurcation, limit cycle, etc.

Lubuma (UP) 2018 33 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

NSFD Schemes in Africa

Lubuma (UP) 2018 34 / 34