J. Nig. Soc. Phys. Sci. 1 (2019) 20–29

Journal of the
Nigerian Society

of Physical
Sciences

Original Research

An Epidemic Model of Zoonotic Visceral Leishmaniasis with
Time Delay

L. Adamu, N. Hussaini∗

Department of Mathematical Sciences, Bayero University, Kano, Nigeria.

Abstract

This paper presents a mathematical model with time delay for the transmission dynamics of zoonotic visceral leishmaniasis (ZVL which is caused
by protozoan parasite leishmania infantum and transmitted by female sandflies). Qualitative analysis of the ODE version of the model reveals that
the disease-free equilibrium of the model is globally asymptotically stable when the basic reproduction number, R0, is less than unity. Using time
delay as a bifurcation parameter in the delay-differential version of the model, it has been shown that the incubation period has significant effect
on the stability of the equilibria and we derived the condition for Hopf bifurcation to occur.

Keywords: ZVL, Stability, Hopf bifurcation, Time delay

Article History :
Received: 01 April 2019
Received in revised form: 05 May2019
Accepted for publication: 06 May 2019
Published: 09 May 2019

c©2019 Journal of the Nigerian Society of Physical Sciences. All rights reserved.
Communicated by: B. J. Falaye

1. Introduction

Leishmaniasis is a disease caused by a group of protozoan
parasite called leishmania. It is transmitted to humans by the
bite of infected female phlebotomine sandfly that fed previously
on an infected reservoir or infected human. More than 20 leish-
mania species are known to be transmitted to humans and over
90 sandfly species are known to transmit leishmania parasites
[1].

There are basically three forms of the diseases namely: Vis-
ceral leishmaniasis (VL), Cutaneous leishmaniasis (CL), and
Mucocutaneous leishmaniasis. VL also known as kala azar
has two major forms which differ in their characteristics trans-
mission, namely (i) zoonotic visceral leishmaniasis (ZVL) and
(ii) Anthroponotic visceral leishmaniasis (AVL). ZVL, infects

∗Corresponding author tel. no: +2348037594137
Email address: nhussaini.mth@buk.edu.ng (N. Hussaini )

mostly children and immunosuppressed individuals, is trans-
mitted from animal reservoir to female phlebotomine sandflies
and then to humans [2]. ZVL posed a serious public health
threat due to the fact that 200,000–400,000 incidences are esti-
mated yearly [3, 4]. As noted in [5], ZVL is endemic in Africa,
Europe (particularly the Mediterranean region) and Asia (par-
ticularly the Indian subcontinent).

A complication of visceral leishmaniasis is known as post
kala-azar dermal leishmaniasis (PKDL) usually appear 6 months
to 1 or more years after treatment. PKDL is self-healing [6, 7,
8]. There is no licensed vaccine against any form of leishma-
niasis including the ZVL (although a number of candidate vac-
cines are at various stages of development and clinical trials) [9-
12]. On the other hand, vaccine against ZVL for reservoir ex-
ists [13, 14]. Furthermore, drugs like miltefosine, paromomycin
and liposomal amphotericin B can be used to cure ZVL patients
[15, 16].

Evidence has shown that when a susceptible human, sand-

20



L. Adamu & N. Hussaini / J. Nig. Soc. Phys. Sci. 1 (2019) 20–29 21

fly or reservoir is infected, there is an incubation period during
which the infectious agent develops. The incubation periods for
humans, sandflies and reservoir range from 2 to 6 months, 8 to 6
days, and 3 to 7 years, respectively [6, 17, 18]. Furthermore, to
capture the effect of the incubation period on the transmission
dynamics of ZVL it is important to incorporate delay caused by
the incubation period in the ZVL model.

Mathematical models have been developed to study the dy-
namics of zoonotic visceral leishmaniasis (for instance, see Refs.
[12, 19-23] and some reference therein). The aforementioned
models do not, however, incorporate an important aspect of the
disease, that is the incubation period. A few mathematical mod-
eling studies, such as those by [11, 24, 25] incorporated time
delay in the transmission dynamics of the vector. Since the in-
cubation period in the reservoir is the longest in comparison to
the time the parasite takes to develop in humans or sandflies.
Thus, it is instructive to study the effect of ZVL latency period
in the reservoir population on the transmission dynamics of the
disease.

The paper is organized as follows. The model is formulated
in Section 2. The ODE version of the model is introduced and
analyzed in Section 3. The analysis of the model with time
delay is carried out in Section 4. The paper is concluded in
Section 5.

2. Model Formulation

Three different populations, at time t, are considered, namely:
human host population denoted by NH (t), vector population de-
noted by NV (t), and reservoir host population denoted by NR(t).
The total human host population is compartmentalize into four
subpopulation namely: susceptible humans (S H (t)), ZVL in-
fected individuals (IH (t)), individuals who developed PKDL
(PH (t)), and individuals who recovered from ZVL (RH (t)), so
that:

NH (t) = S H (t) + IH (t) + PH (t) + RH (t).

Similarly, the total population of vector (NV (t)) is sub-divided
into susceptible vector (S V (t)), and ZVL-infected vector (IH (t)),
so that:

NV (t) = S V (t) + IV (t).

Furthermore, the reservoir population is categorized into sus-
ceptible reservoir (S R(t)), and ZVL-infected reservoir (IR(t)),
such that:

NR(t) = S R(t) + IR(t).

The delay-differential model for the transmission dynamics of
zoonotic visceral leishmaniasis in human, sandfly and reservoir
populations is given by the following deterministic autonomous
system of non-linear delay differential equations (a schematic
diagram of the model is depicted in fig 1, the state variables
and the parameters of the model are described in Table 1 and 2,
respectively):

S ′H = ΠH −λH S H IV −µH S H,

I′H = λH S H IV − (ζH + µH + δH )IH,

P′H = (1 − f1)ζH IH − (β + γ + µH )PH,
R′H = f1ζH IH + (γ + β)PH −µH RH,
S ′V = ΠV −λV S V IR −µV S V,

I′V = λV S V IR −µV IV,

S ′R = ΠR −λRS R(t −τ)IV (t −τ) −µRS R,
I′R = λRS R(t −τ)IV (t −τ) −µR IR, (1)

with

N′H = ΠH −µH NH,

N′V = ΠV −µH NV,

N′R = ΠR −µR NR,

where τ ≥ 0 represents the incubation time that the reservoir
needed to become infectious.

The susceptible human population (S H (t)) is recruited at a
constant rate ΠH and diminished as a result of infection from
an infected sandfly (IV ) at per capita rate (λH ) so that the in-
cidence of new infection is given by λH S H IV . The population
is decreased by natural death µH (assumed to be the same in
all compartments of humans). Here we neglect the incubation
period in humans, since in ZVL humans are not contributing to
the disease transmission unlike in AVL [26] ZVL symptomatic
individuals IH are generated at the rate λH and decreased by re-
covery at a rate ζH ;. A fraction f1 of individuals recovered and
the remaining fraction (1 − f1) developed PKDL. The popula-
tion of individuals with PKDL PH is decreased by treatment at
a rate γ or natural recovery at rate β. The population of individ-
uals recovered from ZVL is generated at the rates f1ζH , β and
γ.

The population of susceptible female sandflies is generated
at a constant rate ΠV and diminished following contact with
infected reservoir IR at per capita infection rate λV so that the
incidence of new infection is given by the mass action term
λV S V IR. We neglect the incubation period in sandflies since
it is very short. Furthermore, the population is decreased by
natural death at a rate µV (assumed to be the same in both the
compartments of sandflies).

The latency period in the reservoir is very long, as such we,
therefore introduced time delay τ to capture this period. At time
t, susceptible reservoirs (that are recruited into the population
at constant rate ΠR) get ZVL infection as a result of contact
with sandflies that have been infected at time t − τ (assumed
reservoirs survived the incubation period τ) when bitten by an
infectious sandfly at a rate λR.

Thus, the population of infected reservoirs is generated at
the rate λRS R(t − τ)IV (t − τ) and diminished by natural death
(at a rate µR). We assume that all the parameters of the model
to be non-negative. For the human, sandfly, and reservoir the
corresponding population sizes are asymptotically constant i.e,
NH (t) →

ΠH
µH

as t → ∞, NV (t) →
ΠV
µV

, NR(t) →
ΠR
µR

as t → ∞.
Without loss of generality we assume that NH (t) =

ΠH
µH

, NV (t) =
ΠV
µV

, NR(t) =
ΠR
µR
∀t ≥ 0.

The dynamics of system (1) are qualitatively equivalent to

21



L. Adamu & N. Hussaini / J. Nig. Soc. Phys. Sci. 1 (2019) 20–29 22

the dynamics of the following system given by

S ′H = ΠH −λH S H IV −µH S H,

I′H = λH S H IV − (ζH + µH )IH,
P′H = (1 − f1)ζH IH − (β + γ + µH )PH,

I′V = λV

(
ΠV

µV
− IV

)
IR −µV IV,

I′R = λR

(
ΠR

µR
− IR(t −τ)

)
IV (t −τ) −µR IR. (2)

The initial conditions for the above system take the form: S H (θ) =
φ1(θ), IH (θ) = φ2(θ), PH (θ) = φ3(θ), IV (θ) = φ4(θ), IR(θ) =
φ5(θ), φi ≥ 0, i = 1, 2, 3, 4, 5, θ ∈ [−τ, 0], φi(0) > 0 where
φ = (φi(θ)) ∈ C+ × C+. Here, C denote the Banach space
C = C([−τ, 0],R) of continuous function mapping the inter-
val [−τ, 0] into R. The non-negativity of the cone is defined as
C+ = C([−τ, 0],R+).

Some of the main assumptions made in the formulation of
the model are as follows;

i. Homogeneous mixing among the vector-hosts populations
(that is, a sandfly has an equal chance of biting any hu-
man or reservoir host);

ii. PKDL infected individuals can recover naturally or due
to treatment [7];

iii. Rate of PKDL relapse to ZVL is assumed negligible (hence
not considered)[27];

iv. PKDL is not life-threatening (no PKDL-induced death is
assumed)[12];

v. ZVL-infected individuals recover successfully or devel-
ops PKDL [7];

vi. Treated infected reservoirs do not usually get cured but
develop an immune response that prevents them from be-
coming infectious [28, 29];

vii. No direct transmission within reservoirs and also within
sandflies [2].

viii. Death rate for humans due to ZVL-infection is assumed
to be negligible (hence not accounted for)[18]

3. Analysis of the model without delay

We consider a special case of the delay-differential model
where there is no incubation period (i.e τ = 0) . The ODE
version of the model is given by

S ′H = ΠH −λH S H IV −µH S H,

I′H = λH S H IV − (ζH + µH )IH,
P′H = (1 − f1)ζH IH − (β + γ + µH )PH,

I′V = λV

(
ΠV

µV
− IV

)
IR −µV IV,

I′R = λR

(
ΠR

µR
− IR

)
IV −µR IR. (3)

3.1. Invariant region:
It can be shown that the biologically-feasible region:

Ω =

{
(S H, IH, PH, IV, IR) ∈ R5+ : S H, IH, PH, IV, IR

≥ 0, NH ≤
ΠH

µH
, NV ≤

ΠV

µV
, NR ≤

ΠR

µR

}
,

is positively-invariant (see for instance [5, 30]) i.e all solutions
starting in Ω remain in Ω ∀t ≥ 0. Thus it is sufficient to consider
the dynamics of system (3) in Ω.

3.2. Stability of disease-free equilibrium (DFE)
The model (3) has a DFE obtained by setting the right-hand

sides of the equations in the model to zero, given by

E0 = (S
o
H, I

o
H, P

o
H I

o
V, I

o
R) =

(
ΠH

µH
, 0, 0, 0, 0

)
. (4)

It follows that the associated basic reproduction number of the
model (2) denoted by R0 is defined as:

R0 =

√
λRΠRλV ΠV

µR2µV 2
(5)

3.3. Epidemiological interpretation of R0
We can epidemiologically interpret the terms in the expres-

sion for the threshold quantity R0 as follows. The term gives
the number of secondary infections that one infectious host can
generate only through a vector and reservoir transmission. Sus-
ceptible sandflies get ZVL-infection following effective contact
with an infected reservoir (IR ), the number of susceptible sand-
fly ΠV

µV
generated by infected reservoir near DFE is the prod-

uct of infection rate of infected reservoir λR
ΠV
µV

and the average
life expectancy of the infected reservoir ( 1

µR
). Thus the aver-

age number of new sandfly infections generated by an infected
reservoir is given by λR

µR

ΠV
µV

.

Similarly, the expression λV
µV

ΠR
µR

account for the number of
new reservoir infections caused by an infected sandfly over the
expected infection period. The humans contribute nothing to
the ZVL transmission.

Theorem 3.1. If R0 ≤ 1, then the disease-free equilibrium E0
of (3) is globally asymptotically stable in Ω

Proof. Consider the following Liapunov function

L = λRΠR IH + λRΠR PH + µVµR IV + λV ΠV IR
L′ = λRΠR

[
λH S H IV − (ζH + µH +

]
+λRΠR

[
(1 − f1)ζH IH − (β + γ + µH )PH

]
+µVµR

[
λV

(
ΠV

µV
− IV

)
IR −µV IV

]
+λV ΠV

[
λR

(
ΠR

µR
− IR

)
IV −µR IR

]
L′ = λRλRΠRS H IV − (µHλRΠR + λRΠR f1ζH )IH

22



L. Adamu & N. Hussaini / J. Nig. Soc. Phys. Sci. 1 (2019) 20–29 23

μH μH 

SH 

IH SV 

IV 

IR 

SR 

RH PH 

ΠH 

μH 

μH 

μV 

μV 

μR 

μR 

ΠR 

ΠV 

δHf1IH (1-f1)δHIH 

λH λV λR 

(β+γ)PH 

Figure 1: Schematic diagram of the model (1). Solid arrows indicate transitions and dashed arrow indicates interaction. Expressions next to arrows
show the per capita flow rate between compartments.

Table 1: Description of the state variables of the model.

Variable Interpretation
NH Total population of humans
S H Population of susceptible humans with risk of ZVL infection
IH Population of ZVL-infected humans with symptoms of ZVL
PH ZVL-infected humans with PKDL
RH ZVL-recovered individuals
NV Total population of Sandflies

S V Population of susceptible Sandflies
IV Population of ZVL-infected Sandflies
NR Total population of reservoir

S R Population of susceptible reservoir with risk of ZVL-infection
IR ZVL-infected reservoir with symptoms of ZVL

+

(
µVµRλVλRΠV

µV
− ΠVλVµR

)
IR

−(µVµRλVλR)IR IV

+µ2VµR

 ΠVλVλRΠR
µ2Vµ

2
R

− 1
 IV

≤ µ2µR

 ΠVλVλRΠR
µ2vµ

2
R

− 1
 IV

L′ = µ2µR(R
2
0 − 1)IV ≤ 0

if R0 ≤ 1. Thus if R0 ≤ 1 the derivative L′ = 0 if and only if and
IV = 0. Furthermore the case R0 = 1 the derivative L′ = 0, con-
sequently the largest compact invariant set in {(S H, IH, PH, IV, IR) ∈

23



L. Adamu & N. Hussaini / J. Nig. Soc. Phys. Sci. 1 (2019) 20–29 24

Table 2: Description of the state variables of the model.

Parameter Interpretation
ΠH ,ΠV , ΠR Recruitment rate for humans,vector and reservoir
λH Rates of direct transmission in humans
λV Rates of direct transmission in sandflies
λR Rates of direct transmission in reservoir
µH,µV,µR Natural death rates of humans, sandflies and reservoir ,

respectively
ζH Recovery rate of a human from ZVL-infection
f1 Fraction of those who recovered from ZVL
τ Incubation period of the disease
β PKDL natural recovery rate
γ PKDL treatment rate

Ω : L′ = 0} when R0 ≤ 1, is the singleton set {E0}, therefore by
Lasalle-invariance principle [31], every solution that starts in Ω
approaches E0 as t →∞.

3.4. Existence of an endemic equilibrium and its stability

In order to obtain an endemic equilibrium where at least one
of the infected component is non zero, we let E∗ = (S ∗H, I

∗
H, P

∗
H, I

∗
V, I

∗
R)

be an arbitrary endemic equilibrium of the model (3), solving
(3) at equilibrium, the non-trivial equilibrium is given by:

S ∗H =
λRµV ΠH (λV ΠR + µVµR)

g1
,

I∗H =
ΠHλHµ

2
Vµ

2
R(R0 − 1)

(ζH + µH + δH )g1
,

P∗H =
λH ΠHµ

2
Vµ

2
R(R0 − 1)ζH (1 − f1)

g2g1
,

I∗V =
µ2Vµ

2
R(R0 − 1)

(µVµR + ΠRλV )λRµV
,

I∗R =
µ2Vµ

2
R(R0 − 1)

(ΠVλR + µVµR)λVµR
,

where,

g1 = µ
2
Vµ

2
RλH (R0 − 1) + λRµVµH (λV ΠR + µVµR),

g2 = (ζH + µH )(β + γ + µH ).

Lemma 3.2. The endemic equilibrium E∗ of model (3) is locally
asymptotically stable if R0 ≥ 1.

Proof. We linearize the system (3) at an endemic equilibrium
to obtain the corresponding characteristic equation as follows:

det



K1 −λ 0 0 −λH S ∗H 0

−λH I∗V K2 −λ 0 λH S
∗
H 0

0 (1 − f1)ζH K3 −λ 0 0

0 0 0 K4 −λ λV (
ΠV
µV
− I∗V )

0 0 0 λR(
ΠR
µR
− I∗R) K5 −λ


= 0,

(6)

where,
K1 = −λH I∗V − µH, K2 = −(ζH + µH ), K3 = −(β + γ + µH ),
K4 = −(λV I∗R + µV ), K5 = −(λR I

∗
R + µR).

Then we have the following characteristics equation

(K1 −λ) (K2 −λ) (K3 −λ) (λ
2

+ (λV I
∗
R

+ λR I
∗
V + µV + µR)λ + λVµR I

∗
R + µVλR I

∗
V

+ λVλR I
∗
V I
∗
R) = 0,

(7)

which has negative roots λ1 = −(ζH + µH ) < 0, and λ2 = −(β +
γ + µH ) < 0. Furthermore, λ3 = −(λH I∗V + µH ) < 0 if R0 > 1
(i.e., I∗V > 0).
Other roots are given by the following equation

λ2 + (λV I
∗
R + λR I

∗
V + µV + µR)λ

+ λVµR I
∗
R + µVλR I

∗
V + λVλR I

∗
V I
∗
R = 0.

(8)

According to Descarte’s Rule of sign changes, if R0 ≥ 1 (i.e.,
I∗V ≥ 0, and I

∗
R ≥ 0), then the above equation has no sign

changes in its coefficients and therefore there will be no roots
whose real part is positive. This shows that the the endemic
equilibrium is stable. Furthermore, if R0 ≤ 1, then the above
equation has sign changes in its coefficients. Therefore one of
the roots has its real part positive and the endemic equilibrium
is unstable.

4. Analysis of the model with delay

Next we analyze the dynamics of system (2) where the time
delay is non-zero (τ , 0).

Theorem 4.1. The disease-free equilibrium E0 of model (2) is
locally asymptotically stable if R0 < 1 and unstable if R0 > 1.

Proof. We linearize the system (2) with delay around the disease-
free equilibrium, then the characteristics equation correspond-

24



L. Adamu & N. Hussaini / J. Nig. Soc. Phys. Sci. 1 (2019) 20–29 25

ing to the jacobian matrix is given by:

det



λ + µH 0 0 λH
ΠH
µH

0

0 λ + P1 0 −λH
ΠH
µH

0

0 − (1 − f ) ζH λ + P2 0 0

0 0 0 λ + µV −
λV ΠV
µV

0 0 0 −λR ΠR e
−λτ

µR
λ + µR


= 0,

(9)
where,

P1 = ζH + µH, P2 = β + γ + µH. Thus,

(λ + µH )(λ + ζHµH )(λ + β + γ + µH )(−µV
µRλ

2
− hλ−µV

2µR
2

+ λV ΠVλRΠRe
−λτ) = 0,

(10)

where h = µVµR2 + µV 2µR.
The equation (10) has the following roots with negative real
part, namely: λ = −µH , λ = −(µH + ζH ), λ = −(β + γ + µH ) and
the remaining roots of

F(λ,τ) = −µVµRλ
2
− (µVµR

2
+ µV

2µR)λ

−µV
2µR

2
+ λV ΠVλRΠRe

−λτ

= λ2 + (µV + µR)λ + µVµR(1 − R0e
−λτ) = 0.

(11)

For τ = 0 we obtain the same characteristics polynomial as in
ODE case (without delay) i.e,

F(λ, 0) = λ2 + (µV + µR)λ + µVµR(1 − R0) = 0. (12)

Again, it follows from Descarte’s rule of signs changes that
equation (12) has roots with negative real part whenever R0 < 1
and thus, the disease free equilibrium is stable. If R0 > 1 then
equation (12) has at most one of its root whose real part is pos-
itive and the disease free equilibrium is unstable.

For τ , 0, suppose R0 > 1, we show that equation (11) has
a positive roots and the disease free equilibrium is unstable. To
see this, we rearrange the equation (11) in form:

λ2 + (µV + µR)λ = µVµR(R0e
−λτ
− 1). (13)

Suppose λ ∈ R, let H(λ) be the left-hand side of (13) and G(λ)
be the right-hand side. Then H(0) = 0 and limλ→∞ H(λ) =
∞. Moreso, the function G(λ) is decreasing function of λ and
G(0) = µVµR(R0 − 1) > 0 as such the two functions must inter-
sect for some λ∗ > 0. As such equation (13) has a positive real
solution and the disease-free equilibrium is unstable.

Now for R0 < 1 we claim equation (13) does not have a
non-negative real roots. In this case for λ ≥ 0, H(λ) is in-
creasing and G(λ) is still decreasing function of λ but G(0) =
µVµR(R0 − 1) < 0. Thus if equation (13) has roots whose real
parts is non-negative, there must be complex conjugate roots
which cross the imaginary axis. Also by Rouché’s theorem, if
instability occur for a particular value of delay τ, a character-
istic root of equation (13) intersect the imaginary axis. Con-
sequently, equation (13) must have a pair of purely imaginary

solutions for some τ > 0. Assume that λ = iω, and without
loss of generality assume that ω > 0 is a root of (13) that is ω
satisfies the following equation:

−ω2 + (µV + µR)iω + µVµR −
λV ΠVλRΠR
µVµR

(cos(ωτ) − i sin(ωτ)) = 0.

Separate the real and the imaginary part and obtain the follow-
ing:

−ω2 + µVµR =
λV ΠV λR ΠR

µV µR
cos(ωτ), (14)

(µV + µR)ω = −
λV ΠV λR ΠR

µV µR
sin(ωτ). (15)

Square both sides of each equation above and add the squared
equations to obtain the following forth order equations in ω :

ω4 + (µV + µR)ω
2

+ µ2Vµ
2
R −

λ2V Π
2
Vλ

2
RΠ

2
R

µ2Vµ
2
R

= 0. (16)

Let σ = ω2 so that equation (16) can be reduced to quadratic
form as follows:

σ2 + (µV + µR)σ + µ
2
Vµ

2
R −

λ2V Π
2
Vλ

2
RΠ

2
R

µ2Vµ
2
R

= 0. (17)

We denote the coefficient as:

a11 = (µ
2
V + µ

2
R) > 0,

a12 = µ
2
Vµ

2
R −

λ2V Π
2
Vλ

2
RΠ

2
R

µ2Vµ
2
R

=

(
µVµR −

λV ΠVλRΠR
µVµR

) (
µVµR +

λV ΠVλRΠR
µVµR

)
= µVµR

1 − λV ΠVλRΠR
µ2Vµ

2
R

 (µVµR + λV ΠVλRΠR
µVµR

)
= µVµR(1 − R0)

(
µVµR +

λV ΠVλRΠR
µVµR

)
> 0,

i f R0 < 1.

Then equation (17) can be written as:

σ2 + a11σ + a12 = 0. (18)

It follows that a12 is positive whenever R0 < 1. Thus, the two
roots of equation (18) have positive product which means that
they are complex or they are real but they have the same sign.
Moreso, they have negative sum which implies they are either
real and negative or complex conjugate with negative real parts.
Hence, equation (18) does not have positive real roots which
leads to the conclusion that there is no ω such that iω is a solu-
tion of equation (18). Therefore, it follows from Kuang’s the-
orem [Ref. [32], pp 83] that all the eigenvalues of the char-
acteristics equation (13) have negative real parts for all delay
values τ ≥ 0. This implies that the disease free equilibrium E′
is locally asymptotically stable if R0 < 1.

25



L. Adamu & N. Hussaini / J. Nig. Soc. Phys. Sci. 1 (2019) 20–29 26

4.1. Hopf bifurcation analysis

Let R0 > 1 (i.e, the endemic equilibrium E∗ exists) and τ be
bifurcation parameter. The characteristics equation is obtained
in form of transcendental equation after linearizing system (2)
at an endemic equilibrium E∗ and obtained the roots (−λH I∗V −
µH ), −(ζH + µH ), −(β + γ + µH ) and the roots of equation (19)
can be written in the form:

λ2 + b2λ + b3 = e
−λτ (G1λ + G2) (19)

where,

b2 = µR + µV + λV I
∗
R,

b3 = µRλV I
∗
V + µVµR,

G1 = −λR I
∗
V,

G2 = −µVλR I
∗
V +

λVλRΠV ΠR
µVµR

−
λVλRΠV I∗R

µV
−
λVλRΠR I∗V

µR

Now from system (3) we have

λV ΠV
µV

= λV I
∗
V +

µV I∗V
I∗R

and

λRΠR
µR

= λR I
∗
R +

µR I∗R
I∗V

(20)

Now,

G2 = −µVλR I
∗
V +

λV ΠV
µV

λRΠR
µR

−
λV ΠV
µV

λR I
∗
R −

λRΠR
µR

λV I
∗
V

Substituting (20) in G2 we obtained:

G2 = −µVλR I
∗
V + (λV I

∗
V +

µV I∗V
I∗R

)(λR I
∗
R +

µR I∗R
I∗V

)

− (λV I
∗
V +

µV I∗V
I∗R

)λR I
∗
R − (λR I

∗
R +

µR I∗R
I∗V

)λV I
∗
V

= −µVλR I
∗
V + λVλR I

∗
V I
∗
R + λVµR I

∗
R + λRµV I

∗
V + µVµR

−λVλR I
∗
V I
∗
R −λRµV I

∗
V −λVλR I

∗
V I
∗
R −λVµR I

∗
R

= −µVλR I
∗
V + µVµR −λVλR I

∗
V I
∗
R.

Therefore,
G2 = µVµR −λR I

∗
V (λV I

∗
V + µV ).

For τ = 0, it follow from Lemma (3.2) above that the en-
demic equilibrium is locally asymptotically stable. Moreover
we claim that for any τ > 0, equation (19) does not have pos-
itive real solution. To see this, we have b2 > 0, b3 > 0. We
move the positive terms from the right-hand side to left-hand
side. The rewritten equation (19) takes the form:

λ2 + b2λ + b̃3 = e
−λτ

(
G1λ + G̃2

)
, (21)

where b̃3 = b3 − e−λτ (µVµR) > 0, ∀λ ≥ 0, ∀τ > 0. On the
other hand G1 < 0, and G̃2 = −λR I∗V (λV I

∗
V + µV ) < 0. Con-

sequently, for all λ ≥ 0 the left-hand side in equation (21) is
positive while the right hand side is negative and the two can-

not be equal for all λ ≥ 0. To obtain the complex conjugate
solutions with positive real parts. Let λ = iω with ω > 0 be the
root of equation (19). Substituting λ = iω into equation (19),
the following equation is obtained:

−ω2 + b2iω + b3 = (cos(ωτ) − i sin(ωτ))(G1iω + G2). (22)

Separating the real and imaginary parts, we obtain that

b3 −ω
2

= G2 cos(ωτ) + G1ω sin(ωτ), (23)
b2ω = G1ω cos(ωτ) − G2 sin(ωτ), (24)

We obtained a polynomial equation by eliminating the trigono-
metric term, this is done by squaring both side of each equation
above and adding the resulting equation we have:

ω4 + (b22 − 2b3 − G
2
1)ω

2
+ b23 − G2

2
= 0. (25)

Notice that this a forth degree polynomial and the delay, τ, has
been eliminated.
Let ω2 = σ ∈ R the equation (25) become a quadratic in σ:

σ2 + ασ + κ = 0, (26)

where,

α = b22 − 2b3 − G
2
1,

κ = b23 − G
2
2.

(27)

We establish the conditions for endemic equilibrium E∗ to be
locally stable, that is equation 25 cannot have a purely imagi-
nary solutions.

Lemma 4.2. The following results follows from polynomial equa-
tion (26)
(i). If κ < 0 or κ ≥ 0 and α < 0, then equation (26) has only
one positive root.
(ii). If κ > 0 and α > 0 or ∆ = α2 − 4κ < 0, then equation (26)
has no positive root.
(iii). If κ > 0, α < 0 and ∆ = α2 − 4κ ≥ 0, then equation (26)
has at least one positive root.

Suppose that condition (ii) of Lemma (4.2) is satisfied then
endemic equilibrium is stable, that is equation (26) does not
have a positive real solution.

Conditon (ii) of lemma 4.2 shows that for τ > 0 there is no
positive σ such that iσ is an eigenvalue of the characteristics
equation (19). Therefore, it follows from Kuang [32] any root
λ of (19) satisfies the relation Reλ < 0.

Theorem 4.3. Assume that condition (ii) of lemma(4.2) is sat-
isfied and R0 > 1, then the endemic equilibrium E∗ of system
(2) is absolutely stable, that is E∗ is asymptotically stable for
all values of delay τ ≥ 0 .

Remarks. Theorem 4.3 indicates that for all values of de-
lay, the endemic equilibrium E∗ of system (2) is asymptotically
stable if the parameters satisfy conditions (ii) of Lemma(4.2),
The stability of the endemic equilibrium will depend on the

26



L. Adamu & N. Hussaini / J. Nig. Soc. Phys. Sci. 1 (2019) 20–29 27

value of the delay if the conditions in Theorem 4.3 are not sat-
isfied. As the delay varies the endemic equilibrium can lose
stability which can lead to oscillations.

For instance equation (26) has at least one positive root say
σ0 if κ < 0 denoted by ω =

√
σ0.

Using time delay (i.e incubation time) τ as the bifurcation
parameter, we drive the condition for Hopf bifurcation to occur.
We can think of the roots of (21), λ as a continuous function in
terms of the delay parameter λ(τ).

Theorem 4.4. Suppose that R0 > 1 and if condition (i) of
Lemma (4.2) is satisfied, then the endemic equilibrium E∗ of
the delay model (2) is asymptotically stable when τ ∈ [0,τ0),
and unstable when τ > τ0 provided that ω0 is the largest pos-
itive simple root of equation (25), furthermore the delay model
(2) undergoes a Hopf bifurcation at an endemic equilibrium E∗

when τ = τ0 where,

τ0 =
1
ω0

arcot
 G2(b3 −ω20) + G1b2ω20
G1(ω30 − b3ω0) + G2b2ω0


Proof. Let λ(τ) = η(τ) + iσ(τ) be the eigenvalue of equation
(21) such that for some initial value of the bifurcation τ0 we
have η(τ0) = 0, and ω(τ0) = ω0 (without loss of generality we
may assume ω0 > 0 ) from equation (23) and (24) we obtained a
sequence of positive values of τn corresponding to any positive
root ω given by:

τn =
1
ω0

arcot
 G2(b3 −ω20) + G1b2ω20
G1(ω30 − b3ω0) + G2b2ω0

 + 2nπ
ω0

,

n = 0, 1, 2, ...

It follows from lemma (3.2) that for τ = 0, E∗ is asymptoti-
cally stable, by Kuang’s theorem [Ref[32], P.83], the endemic
equilibrium E∗ is asymptotically stable for τ ∈ [0,τ0) and un-
stable for all τ > τ0. Now we show the transversal condition

dReλ(τ)
dτ

|τ=τo> 0

holds. When τ passes the critical value τ0 (i.e τ > τ0) by conti-
nuity, the real part of λ(τ) becomes positive and the steady state
become unstable. Moreover, a Hopf bifurcation occur when τ
passes through a critical value τ0 (see Ref. [33])
To establish Hopf bifurcation at τ = τ0 we need to show that
dRe(λ(τ))

dτ > 0. Differentiating equation (21) with respect to τ we
have

(2λ + b2)
dλ
dτ

=[−τe−λτ(G1λ + G2) + G1e
−λτ]

dλ
dτ

−λe−λτ(G1λ + G2).

This gives(
dλ
dτ

)−1
=

2λ + b2 + τe−λτ(G1λ + G2) − G1e−λτ

−λe−λτ(G1λ + G2)

=
2λ + b2

−λe−λτ(G1λ + G2)
+

G1
λ(G1λ + G2)

−
τ

λ

=
λ2 − b3 + e−λτ(G1λ + G2)
−λ2(e−λτ(G1λ + G2))

+
G1

λ(G1λ + G2)
−
τ

λ

=
λ2 − b3

λ2(λ2 + b2λ + b3)
−

G2
λ2(G1λ + G2)

−
τ

λ

Thus,

sign
{( d(Reλ)

dτ

)}
= sign

{
Re

(
dλ
dτ

)−1}
λ=iω0

= sign
{
Re

[
λ2 − b3
−λ2(Z(λ))

]
λ=iω0

+ Re
[

−G2
λ2(G1λ + G2)

]
λ=iω0

}
= sign

{
Re

 −ω20 − b3
ω20(b3 −ω

2
0 + b2ω0i)


+ Re

 G2
ω20(G2 + G1ω0i)

}
= sign

{ ω40 + G22 − b23
ω20((b3 −ω

2
0)

2 + b22ω
2
0)

}
= sign

{ 2ω20 + (b22 − 2b3 − G21)
(b2 −ω20)

2 + b22ω
2
0

}
,

where Z(λ) = λ2 + b2λ + b3. Since f (σ) = σ2 + ασ + κ, thus,

d f (σ)
dσ

= 2σ + α = 2σ + (b22 − 2b3 − G
2
1)

Since ω0 is the simple largest positive root, we have

d f (σ)
dσ

∣∣∣∣∣
σ=ω20

> 0.

Hence

dReλ
dτ

∣∣∣∣∣
ω=ω0,τ=τ0

=

d f (ω20 )
dσ

(b2 −ω20)
2 + b22ω

2
0

> 0,

Remarks If an endemic equilibrium exists and the param-
eters κ < 0 or κ ≥ 0 and α < 0 is satisfied, an increase in the
length of the time delay of ZVL disease transmission below the
critical value of the time delay τ0 will cause an endemic equilib-
rium be stable. As the time delay progresses beyond the critical
delay τ0, a locally asymptotically stable endemic equilibrium
will lose its stability.

27



L. Adamu & N. Hussaini / J. Nig. Soc. Phys. Sci. 1 (2019) 20–29 28

5. Numerical simulation

In this section, numerical simulations were performed for
the model (2). The following parameter values λH = 0.5,β =
0.016,λH = 0.5,µH = 0.0067,µV = 0.087,µR = 0.1,γ =
0.033,δH = 0.011, f1 = 0.64, ΠH = 0.05,ζH = 0.036 are
used. We observed that the number of infected sandflies, in-
fected reservoir and infected humans decay faster for higher
values of the delay (incubation period) as shown in Figures 2,
3, and 4. Therefore, the delay has impact on the transmission
dynamics of the disease.

0 10 20 30 40 50
0

0.5

1

1.5

Time

N
u

m
b

e
r 

o
f 

In
fe

c
te

d
 H

u
m

a
n

s

 

 
tau=1
tau=2
tau=3
tau=4
tau=6

Figure 2: Numerical simulation of model (2) with parameter values
λV = 0.250,λR = 0.20, ΠV = 0.0026, ΠR = 0.022,.

0 10 20 30 40 50
0

0.2

0.4

0.6

0.8

1

Time

N
u

m
b

e
r 

In
fe

c
te

d
 s

a
n

d
fl

ie
s

 

 
tau=1
tau=2
tau=3
tau=4
tau=5

Figure 3: Numerical simulation of model (2) with parameter values
λV = 0.50,λR = 0.1, ΠV = 0.0016, ΠR = 0.052.

0 10 20 30 40 50
0

0.2

0.4

0.6

0.8

1

Time

N
u

m
b

e
r 

o
f 

In
fe

c
te

d
 r

e
se

rv
o

ir

 

 
tau=0.01
tau=1
tau=3
tau=5
tau=8

Figure 4: Numerical simulation of model (2) with parameter values
λV = 0.250,λR = 0.20, ΠV = 0.005, ΠR = 0.000523,.

6. Conclusions

This paper presents a new deterministic delay-deferential
model for assessing the effect of time delay on the dynamics
of ZVL by looking at the stability of the equilibria. The time
delay accounts for the incubation period of reservoirs needed to
become infectious. The main theoretical and epidemiological
findings of the study are summarized as follows:

The ODE version of the model has a DFE (ZVL-free equi-
librium) which is locally-asymptotically stable whenever a cer-
tain threshold quantity (R0) is less than unity. If R0 < 1 the
disease-free equilibrium is globally asymptotically stable. Also
a unique endemic equilibrium exists and is locally asymptoti-
cally stable in the interior of the feasible region.

Base on the delay-deferential model we determined the cri-
teria for Hopf-bifurcation to occur using time delay as the bifur-
cation parameter, when time delay is small the positive equilib-
rium is locally asymptotically stable, while instability can occur
by Hopf-bifurcation as the delay increases. Hopf-bifurcation
has helped us in finding the region of instability in a neighbor-
hood of endemic equilibrium where the population will undergo
regular oscillations.

Acknowledgments

We thank the referees for the positive enlightening com-
ments and suggestions, which have greatly helped us in making
improvements to this paper.

References

[1] World Health Organization (2016), fact sheet http://www.who.int/
mediacentre/factsheets/fs375/en/. (Accessed September 2016).

[2] N. Hartemink, S. O. Vanwambeke, H. Heesterbeek, D. Rogers, D. Mor-
ley, B. Pesson, C. Davies, S. Mahamdallie & P. Ready, “Integrated Map-
ping of Establishment Risk for Emerging Vector-Borne Infections: A

28



L. Adamu & N. Hussaini / J. Nig. Soc. Phys. Sci. 1 (2019) 20–29 29

Case Study of Canine Leishmaniasis in Southwest France”, PLoS ONE 6
(2011) 20817.

[3] H. Seifu, L. Torleif, H.M. Damen, W. Tassew & L. Bernt, “Climate
change, crop production and child under nutrition in Ethiopia, a longi-
tudinal panel study”, BMC public Health 14 (2014) 884.

[4] H.R. Thieme, “Convergence results and Poincare-Bendixson trichtomy
for asymptotically autonomous differential equation”, Journal of Mathe-
matical Biology 30 (1992) 755.

[5] N. Hussaini, J. M-S Lubuma, K. Barley & A. B. Gumel, “Mathematical
analysis of a model for AVL–HIV co-endemicity”, Mathematical Bio-
sciences 271 (2016) 80.

[6] Spickler, Anna Rovid, Technical factsheet (2009) @ http:
//www.cfsph.iastate.edu/DiseaseInfo/disease.php?name=

leishmaniasis&lang=en (Accessed June 2017).
[7] World Health Organization Post-kala-azar dermal leishmaniasis: a man-

ual for case management and control: report of a WHO consultative
meeting, Kolkata, India, 2–3 July 2012.

[8] B. M. Younis , H. A. A. Mohammed, M. M. M. Dafalla, A. O. A. Adam,
M. Y. Elamin, A. M. Musa, A. M. El-Hassan & E. A. G. Khalil, “Cure of
post Kala-azar dermal leishmaniasis with paromomycin/sodium stiboglu-
conate combination: a proof of concept”, Int. J. Res. Med. Sci. 3 (2015)
16.

[9] J. K.Evans and L. Kedzierski, “Development of vaccines against visceral
leishmaniasis”, Journal of tropical medicine, Article ID 892817 (2012)
14.

[10] E. Handman, “Development Leishmaniasis: Current Status of Vaccine”,
Clin. Microbiol. Rev. 14 (2001) 229.

[11] L. Ribas, V.L. Zaher, H.J. Shimozako & E. Massad, Estimating the Op-
timal Control of Zoonotic Visceral Leishmaniasis by the Use of a Mathe-
matical Model, The Scientific World Journal , 2013.

[12] A. Stauch, R. R. Sarkar, A. Picado, B. Ostyn, S. Sundar, S. Rijal, M.
Boelaert, J. C. Dujardin & H. P. Duerr, “Visceral leishmaniasis in the
indian subcontinent: modelling epidemiology and control”, PLoS Negl.
Trop. Dis.5 (2011) 1405.

[13] C.B. Palatnik-de-sousa, “Vaccines for canine leishmaniasis, frontiers in
immunology”, 3 (2012) 69

[14] A. P. Seva, F. G Ovallos, M. Amaku, E. Carrillo, J. Moreno, E. A. Galati,
E. G. Lopes, R. M. Soares & F. Ferreira “Canine-Based Strategies for
Prevention and Control of Visceral Leishmaniasis in Brazil”, PLoS ONE
11 (2016) 0160058.

[15] F. Chappuis, S. Sundar, A. Hailu, H. Ghalib, S. Rijal, W. R. Peeling, J.
Alvar & M. Boelaert, “Visceral leishmaniasis: what are the needs for di-
agnosis,treatment and control?”, Nature Reviews Microbiology 5 (2007)
7

[16] E. M. Moore & D. N. Lockwood, Treatment of Visceral Leishmania-
sis, Hospital for Tropical Diseases, University College London Hospital,
London School of Hygiene and Tropical Medicine 2 (2010).

[17] S. S. Menon, R. Rossi, L. Nshimyumukiza & K. Zinszer, “Decentralized
control of human visceral leishmaniasis in endemic urban areas of Brazil:
a literature review”, Tropical Medicine and Health 44 (2016) 9.

[18] P. D. Ready, “Epidemiology of visceral leishmaniasis”, J. Clinical Epi-
demiology 6 (2014) 147.

[19] I. M. ELmojtaba, J. Y. T. Mugisha & M. H. A. Hashim, “Mathematical
analysis of the dynamics of visceral leishmaniasis in the Sudan”, Appl.
Math. Comput. 217(2010) 2567.

[20] I. M. ELmojtaba, J. Y. T. Mugisha & M. H. A. Hashim, “Vaccination
model for visceral leishmaniasis with infective immigrants”, Math. Meth.
Appl. Sci. 36 (2013) 216.

[21] Z. Muhammad, R. Ali, “Zoonotic Visceral Leishmania: Modeling and
Control”, J. Appl. Comput. Math. 4 (2015) 4.

[22] A. Subramanian, V. Singh, R.R. Sarkar, “Understanding Visceral Leish-
maniasis Disease Transmission and its Control–A Study Based on Math-
ematical Modeling”, J. Mathematics 3 (2015) 913.

[23] S. Zhao, Y. Kuang, C. Wu, D. Ben-Arieh, M. Ramalho-Ortigao & K.
Bi, “Zoonotic visceral leishmaniasis transmission: modeling, backward
bifurcation, and optimal control”, Journal of Mathematical Biology 73
(2016) 1525.

[24] M. N. Burattini, F. A. B. Coutinho, L. F. Lopez & E. Massad, “Modelling
the dynamics of leishmaniasis considering human, animal host and vector
populations”, J. Biol. Syst. 6 (1998) 337.

[25] H. J. Shimozako, J. Wu & E. Massad, “Mathematical modelling for
Zoonotic Visceral Leishmaniasis.A new analysis considering updated pa-
rameters and notified human Brazilian data”, Infectious Disease Mod-
elling 2 (2017) 143.

[26] R. J. Quinnell, O. Couetenay, “Transmission, reservoir hosts and control
of zoonotic visceral leishmaniasis,Institute of Integrative and Compara-
tive Biology”, University of Leeds, Leeds LS2 9JT, UK.

[27] E. E. Zijlstra, A. M. Musa, E. A. G. Khalil, I. M. El-Hassan & A. M.
El-Hassan, “Post-kala- azar dermal leishmaniasis”, Lancet Infect. Dis. 3
(2003) 87.

[28] G. Baneth & S. E. Shaw, “Chemotherapy of canine leishmaniosis”, Vet-
erinary Parasitology 106 (2002) 315.

[29] L. A. Espejo, S. Costard & F. J. Zagmutt, “Modeling canine leishmaniasis
spread to non-endemic areas of Europe”, Epidemiol. Infect. 143 (2015)
1936.

[30] N. Hussaini, M. Winter & A. B. Gumel, “Qualitative assessment of the
role of public health education program on HIV transmission dynamics”,
Math. Med. Biol.: J. IMA 28 (2011) 245.

[31] J. K. Hale, Ordinary Differential Equations, Wiley, NewYork, 1969.
[32] Y. Kuang, Delay Differential Equations with Applications in Population

Dynamics, Academic Press, New York, 1993.
[33] B. D. Hassard, N. D. K. Azarinoff & Y. H. Wan, Theory and Applications

of Hopf Bifurcation, Cambridge University, Cambridge, (1981).

29