An Epidemic Model of Zoonotic Visceral Leishmaniasis with Time Delay

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.


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 leishmania 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: Visceral leishmaniasis (VL), Cutaneous leishmaniasis (CL), and Mucocutaneous leishmaniasis. VL also known as kala azar has two major forms which differ in their characteristics transmission, namely (i) zoonotic visceral leishmaniasis (ZVL) and (ii) Anthroponotic visceral leishmaniasis (AVL). ZVL, infects 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 dynamics of zoonotic visceral leishmaniasis (for instance, see Refs. [12,[19][20][21][22][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 modeling studies, such as those by [11,24,25] incorporated time delay in the transmission dynamics of the vector. Since the incubation 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.

Model Formulation
Three different populations, at time t, are considered, namely: human host population denoted by N H (t), vector population denoted by N V (t), and reservoir host population denoted by N R (t). The total human host population is compartmentalize into four subpopulation namely: susceptible humans (S H (t)), ZVL infected individuals (I H (t)), individuals who developed PKDL (P H (t)), and individuals who recovered from ZVL (R H (t)), so that: Similarly, the total population of vector (N V (t)) is sub-divided into susceptible vector (S V (t)), and ZVL-infected vector (I H (t)), so that: Furthermore, the reservoir population is categorized into susceptible reservoir (S R (t)), and ZVL-infected reservoir (I R (t)), such that: 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): with 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 (I V ) at per capita rate (λ H ) so that the incidence of new infection is given by λ H S H I V . 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 I H are generated at the rate λ H and decreased by recovery at a rate ζ H ;. A fraction f 1 of individuals recovered and the remaining fraction (1 − f 1 ) developed PKDL. The population of individuals with PKDL P H is decreased by treatment at a rate γ or natural recovery at rate β. The population of individuals recovered from ZVL is generated at the rates f 1 ζ H , β and γ.
The population of susceptible female sandflies is generated at a constant rate Π V and diminished following contact with infected reservoir I R at per capita infection rate λ V so that the incidence of new infection is given by the mass action term λ V S V I R . 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 λ R S R (t − τ)I V (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, (1) are qualitatively equivalent to the dynamics of the following system given by The initial conditions for the above system take the form: The non-negativity of the cone is defined as 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 human 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 develops PKDL [7]; vi. Treated infected reservoirs do not usually get cured but develop an immune response that prevents them from becoming 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

Invariant region:
It can be shown that the biologically-feasible region: 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 Ω.

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 It follows that the associated basic reproduction number of the model (2) denoted by R 0 is defined as:

Epidemiological interpretation of R 0
We can epidemiologically interpret the terms in the expression for the threshold quantity R 0 as follows. The term gives the number of secondary infections that one infectious host can generate only through a vector and reservoir transmission. Susceptible sandflies get ZVL-infection following effective contact with an infected reservoir (I R ), the number of susceptible sandfly Π V µ V generated by infected reservoir near DFE is the product of infection rate of infected reservoir λ R Π V µ V and the average life expectancy of the infected reservoir ( 1 µ R ). Thus the average number of new sandfly infections generated by an infected reservoir is given by Lasalle-invariance principle [31], every solution that starts in Ω approaches E 0 as t → ∞.

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: where, Lemma 3.2. The endemic equilibrium E * of model (3) is locally asymptotically stable if R 0 ≥ 1.
Proof. We linearize the system (3) at an endemic equilibrium to obtain the corresponding characteristic equation as follows: where, . Then we have the following characteristics equation which has negative roots . Other roots are given by the following equation According to Descarte's Rule of sign changes, if R 0 ≥ 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 R 0 ≤ 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.

Analysis of the model with delay
Next we analyze the dynamics of system (2) where the time delay is non-zero (τ 0). ing to the jacobian matrix is given by: where, where h = µ V µ R 2 + µ V 2 µ R . The equation (10) has the following roots with negative real part, namely: λ = −µ H , λ = −(µ H + ζ H ), λ = −(β + γ + µ H ) and the remaining roots of For τ = 0 we obtain the same characteristics polynomial as in ODE case (without delay) i.e, Again, it follows from Descarte's rule of signs changes that equation (12) has roots with negative real part whenever R 0 < 1 and thus, the disease free equilibrium is stable. If R 0 > 1 then equation (12) has at most one of its root whose real part is positive and the disease free equilibrium is unstable. For τ 0, suppose R 0 > 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: 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 (R 0 − 1) > 0 as such the two functions must intersect for some λ * > 0. As such equation (13) has a positive real solution and the disease-free equilibrium is unstable. Now for R 0 < 1 we claim equation (13) does not have a non-negative real roots. In this case for λ ≥ 0, H(λ) is increasing and G(λ) is still decreasing function of λ but G(0) = µ V µ R (R 0 − 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 characteristic root of equation (13) intersect the imaginary axis. Consequently, 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: Separate the real and the imaginary part and obtain the following: Square both sides of each equation above and add the squared equations to obtain the following forth order equations in ω : Let σ = ω 2 so that equation (16) can be reduced to quadratic form as follows: We denote the coefficient as: Then equation (17) can be written as: It follows that a 12 is positive whenever R 0 < 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 solution of equation (18). Therefore, it follows from Kuang's theorem [Ref. [32], pp 83] that all the eigenvalues of the characteristics equation (13) have negative real parts for all delay values τ ≥ 0. This implies that the disease free equilibrium E is locally asymptotically stable if R 0 < 1.

Hopf bifurcation analysis
Let R 0 > 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: where, Now, Substituting (20) in G 2 we obtained: . For τ = 0, it follow from Lemma (3.2) above that the endemic equilibrium is locally asymptotically stable. Moreover we claim that for any τ > 0, equation (19) does not have positive real solution. To see this, we have b 2 > 0, b 3 > 0. We move the positive terms from the right-hand side to left-hand side. The rewritten equation (19) takes the form: Consequently, 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 + b 2 iω + b 3 = (cos(ωτ) − i sin(ωτ))(G 1 iω + G 2 ). (22) Separating the real and imaginary parts, we obtain that We obtained a polynomial equation by eliminating the trigonometric term, this is done by squaring both side of each equation above and adding the resulting equation we have: Notice that this a forth degree polynomial and the delay, τ, has been eliminated.
Let ω 2 = σ ∈ R the equation (25) become a quadratic in σ: where, We establish the conditions for endemic equilibrium E * to be locally stable, that is equation 25 cannot have a purely imaginary solutions.
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 satisfied and R 0 > 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 delay, 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 value of the delay if the conditions in Theorem 4.3 are not satisfied. As the delay varies the endemic equilibrium can lose stability which can lead to oscillations.
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 R 0 > 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 positive simple root of equation (25), furthermore the delay model (2) undergoes a Hopf bifurcation at an endemic equilibrium E * when τ = τ 0 where, 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: It follows from lemma (3.2) that for τ = 0, E * is asymptotically stable, by Kuang's theorem [Ref [32], P.83], the endemic equilibrium E * is asymptotically stable for τ ∈ [0, τ 0 ) and unstable for all τ > τ 0 . Now we show the transversal condition holds. When τ passes the critical value τ 0 (i.e τ > τ 0 ) by continuity, 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 Since ω 0 is the simple largest positive root, we have Remarks If an endemic equilibrium exists and the parameters κ < 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 equilibrium be stable. As the time delay progresses beyond the critical delay τ 0 , a locally asymptotically stable endemic equilibrium will lose its stability.

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 equilibrium) which is locally-asymptotically stable whenever a certain threshold quantity (R 0 ) is less than unity. If R 0 < 1 the disease-free equilibrium is globally asymptotically stable. Also a unique endemic equilibrium exists and is locally asymptotically stable in the interior of the feasible region.
Base on the delay-deferential model we determined the criteria for Hopf-bifurcation to occur using time delay as the bifurcation parameter, when time delay is small the positive equilibrium 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 neighborhood of endemic equilibrium where the population will undergo regular oscillations.