Mathematical analysis of a model for zoonotic visceral leishmaniasis.

Zoonotic visceral leishmaniasis (ZVL), caused by the protozoan parasite Leishmania infantum and transmitted to humans and reservoir hosts by female sandflies, is endemic in many parts of the world (notably in Africa, Asia and the Mediterranean). This study presents a new mathematical model for assessing the transmission dynamics of ZVL in human and non-human animal reservoir populations. The model undergoes the usual phenomenon of backward bifurcation exhibited by similar vector-borne disease transmission models. In the absence of such phenomenon (which is shown to arise due to the disease-induced mortality in the host populations), the nontrivial disease-free equilibrium of the model is shown to be globally-asymptotically stable when the associated reproduction number of the model is less than unity. Using case and demographic data relevant to ZVL dynamics in Arac̣atuba municipality of Brazil, it is shown, for the default case when systemic insecticide-based drugs are not used to treat infected reservoir hosts, that the associated reproduction number of the model (ℛ0) ranges from 0.3 to 1.4, with a mean of ℛ0=0.85 . Furthermore, when the effect of such drug treatment is explicitly incorporated in the model (i.e., accounting for the additional larval and sandfly mortality, following feeding on the treated reservoirs), the range of ℛ0 decreases to ℛ0∈[0.1,0.6] , with a mean of ℛ0=0.35 (this significantly increases the prospect of the effective control or elimination of the disease). Thus, ZVL transmission models (in communities where such treatment strategy is implemented) that do not explicitly incorporate the effect of such treatment may be over-estimating the disease burden (as measured in terms of ℛ0 ) in the community. It is shown that ℛ0 is more sensitive to increases in sandfly lifespan than that of the animal reservoir (so, a strategy that focuses on reducing sandflies, rather than the animal reservoir (e.g., via culling), may be more effective in reducing ZVL burden in the community). Further sensitivity analysis of the model ranks the sandfly removal rate (by natural death or by feeding from insecticide-treated reservoir hosts), the biting rate of sandflies on the reservoir hosts and the progression rate of exposed reservoirs to active ZVL as the three parameters with the most effect on the disease dynamics or burden (as measured in terms of the reproduction number ℛ0 ). Hence, this study identifies the key parameters that play a key role on the disease dynamics, and thereby contributing in the design of effective control strategies (that target the identified parameters).

Introduction

The protozoan Leishmania infantum (syn., L. chagasi) is the causative agent of zoonotic visceral leishmaniasis (ZVL) in humans and canine leishmaniasis (CanL) in dogs (Hartemink et al., 2011, Podaliri Vulpiani et al., 2011, Ribas et al., 2013). The protozoan parasite is transmitted from infected animal hosts (domestic dogs serve as principal reservoirs) to susceptible female sandflies (Diptera: Phlebotomine) and then to susceptible humans (who are regarded as dead-end hosts of the disease) (Elnaiem et al., 2001, Hoogstraal and Heyneman, 1969, Hussaini et al., 2016, Kirk, 1939, Podaliri Vulpiani et al., 2011, Ribas et al., 2013). ZVL, which is endemic in Africa, Europe (particularly the Mediterranean region) and Asia (particularly the Indian subcontinent) (European Centre for Disease Prevention and Control,, Podaliri Vulpiani et al., 2011), is an acute and life-threatening emerging disease with estimated yearly incidence in the range 200 000 to 400 000 (Leta et al., 2014, World Health Organization,). Furthermore, increase in risk factors associated with climate change and other environmental challenges makes ZVL to be a growing major public health concern (Hartemink et al., 2011).

An adult female sandfly lays about eggs during a single gonotrophic cycle (these eggs are typically laid in damp dark places in the cattle sheds, animal burrows, tree roots and in soil rich in organic matter) (European Centre for Disease Prevention and Control,, Sand fly life cycle, ). The eggs laid in these micro-habitats hatch into larvae in days (European Centre for Disease Prevention and Control, ). Larvae develop into four instar stages (each one larger than the one before; the newly hatched first instar larvae have two rear bristles, while all later larval developments have four rear bristles) (European Centre for Disease Prevention and Control, ). Larvae are mainly scavengers found in moist areas, such as animal burrows, feeding on organic matter (e.g., fungi, decaying leaves and animal faeces) (European Centre for Disease Prevention and Control,, Sand fly life cycle, ). During the fourth molt, the larva matures into a pupa (the whole process of maturation from larvae to pupae takes about days depending on species, temperature and nutrient availability) (European Centre for Disease Prevention and Control, ). Pupae then develop into adult sandflies in about days (European Centre for Disease Prevention and Control,, Sand fly life cycle, ). Thus, the duration of the whole cycle, from egg laying to an adult sandfly, varies between 30 and 63 days depending on species, temperature and nutrient availability (European Centre for Disease Prevention and Control, ). Adult sandflies usually mate within a few days after emerging from the pupal stage, after which the female sandfly moves to quest for blood meal required to produce eggs (European Centre for Disease Prevention and Control, ). The feeding activity of the female adult sandfly is influenced by temperature, humidity and air movement (European Centre for Disease Prevention and Control,, Sand fly life cycle, ). Sandflies, which are active and feed during the early morning and evening hours when temperature falls and humidity rises, have an average lifespan of about 14 days (European Centre for Disease Prevention and Control,, Sand fly life cycle, ). A schematic description of the life-cycle of the sandfly is depicted in Fig. 1. Although there is a vaccine against ZVL in animal populations (CaniLeish) (CaniLeish, 2017, Vetlife,), no such vaccine currently exists for use in humans (although a number of candidate vaccines are at various stages of development and clinical trials) (Gillespie et al., 2016, Kumar and Engwerda, 2014, Mcallister, 2014) (it is however, known that an effective vaccine against leishmaniasis will prompt long-lasting immunity in humans (Bertholet et al., 2009, Gillespie et al., 2016, Mcallister, 2014, Nagill and Kaur, 2011)). Furthermore, although ZVL is curable using drugs such as miltefosine, paromomycin and liposomal amphotericin B (Chappuis et al., 2007), basic anti-ZVL preventive measures (such as personal protection against sandfly bites and sandfly-reduction strategies focused on spraying anti-sandfly insecticides in human and animal reservoir habitats) remain perhaps the most effective method for combating ZVL spread in humans (World Health Organization). Treatment of animal reservoir (with systemic insecticide-based drugs, such as fipronil) are implemented in places like Bihar, India (Poché, Grant, & Wang, 2016). An additional benefit of the treatment strategy is that it reduces the number of larvae and adult sandflies who feed on the faeces of (insecticide-based) treated infected reservoirs (Poché et al., 2016).

Fig. 1

Schematic diagram of the life-cycle of the sandfly (Sharma, ).

A number of modeling studies have been carried out to gain insight into ZVL transmission dynamics in human and/or reservoir populations (see, for instance (Burattini et al., 1998, Carvalho et al.,, Hartemink et al., 2011, Hussaini et al., 2016, Ribas et al., 2013, Shimozako et al., 2017, Zhao et al., 2016), and some of the references therein). Burattini et al. (Burattini et al., 1998). proposed mechanistic model for ZVL transmission within the human and animal reservoir populations, and used the model to evaluate control strategies. Ribas et al. (Ribas et al., 2013). added control terms to the model in (Burattini et al., 1998) to estimate the optimal control strategies for ZVL. Zhao et al. (Zhao et al., 2016). developed a model to describe the ZVL transmission dynamic using a modified SEIR model and the model exhibited backward bifurcation phenomenon. Shimozako et al (Shimozako et al., 2017). updated most of parameters in (Burattini et al., 1998, Ribas et al., 2013) and calculated new value of . The current study focuses on the design and analysis of a novel model, which extends some of the aforementioned modeling studies, for assessing the transmission dynamics of ZVL in human and non-human animal reservoir populations. The paper is organized as follows. The model is formulated in Section 2 and rigorously analyzed in Section 3. Sensitivity uncertainty analysis and numerical simulations are reported in Section 4.

Model formulation

The model to be developed monitors the transmission dynamics of zoonotic visceral leishmaniasis (ZVL) within the human and animal (reservoir) host populations. Unlike anthroponotic visceral leishmaniasis (which is transmitted from human to vector to human), ZVL is transmitted from infected animals (reservoir) to susceptible vectors (sandflies) and then back to humans. The total human population at time t, denoted by , is sub-divided into four compartments of susceptible , infected but not infectious (i.e., asymptomatically-infected humans) , symptomatically-infected and recovered humans, so that:

Furthermore, the total sandfly population at time t, denoted by , is split into two main classes of immature and mature adult female phlebotomine sandfly classes. The total immature sandfly population at time t, denoted by , consists of the first three stages of sandfly life-cycle (i.e., eggs , larvae and pupae ). Further, the total mature female sandfly population at time t, denoted by , is split into compartments for susceptible female sandflies and infected female sandflies , so that

Finally, the total animal reservoir population at time t, denoted by , is sub-divided into compartments for susceptible exposed , infected and treated reservoirs, so that

The model for ZVL transmission in human and reservoir animal populations is given by the following deterministic system of non-linear differential equations (a flow diagram of the model is depicted in Fig. 2; the state variables and parameters of the model are described in Table 1, Table 2, respectively):

Fig. 2

Flow chart of model (2.1), where , , , .

Table 1

Description of the variables of the model (2.1).

Variable	Interpretation
S_H(t)	Population of susceptible humans
E_H(t)	Population of humans exposed to ZVL
I_H(t)	Population of humans with clinical symptoms of ZVL
R_H(t)	Population of humans who recovered from ZVL
E_V(t)	Population of sandfly eggs
L_V(t)	Population of sandfly larvae
P_V(t)	Population of sandfly pupae
S_V(t)	Population of susceptible adult female sandflies
I_V(t)	Population of ZVL-infected adult female sandflies
S_R(t)	Population of susceptible ZVL reservoirs
E_R(t)	Population of reservoirs exposed to ZVL
I_R(t)	Population of infected reservoirs with clinical symptoms of ZVL
T_R(t)	Population of ZVL-treated reservoirs

Table 2

Description of parameters of the model (2.1).

Parameter	Interpretation
Π_H(Π_R)	Recruitment rate of humans (reservoirs)
μ_H(μ_R)	Natural death rate of humans (reservoirs)
ψ_V	Oviposition rate
μ_E,μ_L,μ_P,μ_M	Natural death rate of eggs, larvae, pupae and adult sandflies, respectively
β_H(β_R)	Transmission probability from infected sandflies to susceptible human (reservoir) hosts
β_V	Transmission probability from infected reservoirs to susceptible sandflies
b_H(b_R)	Per capita biting rate of sandflies on the human (reservoir) hosts
γ_H(γ_R)	Progression rate of exposed human (reservoir) hosts to active ZVL class
τ_H(τ_R)	Treatment rates of human (reservoir) hosts
σ_E	Average maturation rate from eggs to larvae
σ_L	Average maturation rate from larvae to pupae
σ_P	Average maturation rate from pupae to adult sandflies
ρ_R	Rate of relapse of treated reservoirs
K_M	Carrying capacity of adult sandflies
η_R	Modification parameter for relative of infectiousness of reservoirs
f	Fraction of newly-emerged sandflies that are females
δ_H(δ_R)	Disease-induced death rates of human (reservoir) hosts
ξ_L(ξ_M)	Additional death rate of larvae (adult sandflies) due to feeding on faeces of treated reservoir

In the model (2.1), is the recruitment rate for human (reservoir), is the biting rate of adult female sandflies on the human (reservoir) host, is the probability of infection per bite from an infected adult female sandfly (human) to a susceptible human (sandfly), is the probability of infection per bite from an infected reservoir to a susceptible adult female sandfly, is the natural death rate in humans (reservoir hosts) and accounts for the reduction of infectiousness of exposed reservoirs. Laboratory experiments by Laurenti et al. (Laurenti et al., 2013), show that asymptomatic reservoir transmits ZVL to susceptible sandflies at a rate greater than that of symptomatic reservoir hosts (i.e., ). The parameter measures the rate at which humans (reservoir hosts) in the class develop clinical symptoms of ZVL, while the parameter measures the treatment rate of symptomatic humans (reservoir hosts). The parameters and account, respectively, for the disease-induced death rate and failure rate of treatment received by infected reservoir hosts. It is assumed that recovery confers permanent immunity against ZVL re-infection in humans (Roberts, 2005).

Eggs are laid by adult female sandflies (usually on the surface an organic matter), assumed to be at a logistic rate (where for all t is the carrying capacity of female adult sandflies and is the egg deposition rate). Eggs hatch into the larvae (at a rate ) which, in turn, mature into pupae (at a rate ) and, finally, pupae mature into adult female sandflies (at a rate ). Susceptible adult female sandflies acquire ZVL infection (at the rate , as defined in the caption of Fig. 2) and suffer natural death (at a rate ). Furthermore, adult female sandflies die due to feeding on infected reservoir hosts that have been treated with systemic insecticide-based drugs (Poché et al., 2016) (at a rate ).

The parameters and represent, respectively, the natural death rate for eggs, larvae and pupae, while is the density-dependent mortality rate for larvae (accounting for the cannibalism that occurs during larval competition for resources (nutrients) and space) (Poché et al., 2016, Srinivasan and Panicker, 1992). Finally, as in the case of adult sandflies, larvae also suffer additional mortality by feeding on organic material from (insecticide-based) treated infected reservoir hosts, at a rate (Poché et al., 2016). Following (Poché et al., 2016), the parameters and are defined, respectively, as:where, is the number of days of post-treatment of infected reservoir and is the number of post-defecation days of reservoir. The model (2.1) accounts for the conservation law of sandfly bites on human and reservoir hosts (the consequence of which is that the human (and reservoir) hosts are always sufficient in abundance and the total number of bites made by sandflies balances the total number of bites received by the human (and reservoir) hosts; see also (Bergsman et al., 2016, Bowman et al., 2005, Cruz-Pacheco et al., 2012, Marini et al., 2017, Subramanian et al., 2015) for models of similar diseases with one vector and multiple hosts).

Some of the main assumptions made in the formulation of the model (2.1) are:

Humans are dead-end hosts (i.e., they acquire, but do not transmit, ZVL infection) (Hartemink et al., 2011).

Humans who recovered from ZVL infection acquire permanent immunity against re-infection (i.e., ) (Roberts, 2005).

Treated infected reservoir hosts do not usually get cured but develop an immune response that prevents them from becoming infectious (Baneth and Shaw, 2002, Espejo et al., 2015).

Treated reservoir hosts can relapse to active ZVL class due to treatment failure (i.e., ) (Petersen and Barr, 2009, Quinnell and Courtenay, 2009).

Recovery confers permanent immunity against ZVL re-infection in humans (Hussaini et al., 2016).

No direct transmission between reservoirs or between sandflies is assumed (Hartemink et al., 2011).

The model (2.1) extends the deterministic models for ZVL transmission developed in (Burattini et al., 1998, Ribas et al., 2013, Shimozako et al., 2017, Zhao et al., 2016) by, inter alia,

adding the compartments of immature sandflies (i.e., the compartments , and ).

allowing for the relapse of treated reservoir hosts to active ZVL class due to treatment failure (i.e., ).

adding density-dependent larval mortality (i.e., ).

allowing for additional mortality of sandfly larvae (i.e., ) and adult female sandfly due to feeding on the faeces of treated infected reservoir hosts.

using varying total populations of the human and reservoir hosts (constant population was used in (Burattini et al., 1998, Ribas et al., 2013, Shimozako et al., 2017)).

The model (2.1) is, first of all, fitted using the ZVL case and demographic data from Arac̣atuba municipality, Brazil for the period (tabulated in Table 3) (Centre of Epidemiological Surveillance of Sao Paulo State (CES-SP) and Brazil, 2016, Shimozako et al., 2017). The results obtained, depicted in Fig. 3, show a reasonably good fit to the data (expressed in terms of cumulative number of yearly cases). It is worth mentioning that, for the model fitting, the human demographic parameters ( and ) are parameterize as follows. Since the average total population of Arac̣atuba municipality is (see Table 3), and the average lifespan in Brazil is 75 years (World Bank ata, 2015) (i.e.,  = 75 years, so that per day), it follows from the relation that per day. Furthermore, since systemic insecticide-based drugs were not used to treat infected reservoir hosts in Arac̣atuba municipality, Brazil during the period 1999 to 2015, the associated parameters and (for the treatment of infected reservoir hosts) were set to zero, while all other parameters of the model are set at their baseline values in Table 4.

Table 3

Human reported ZVL cases in Arac̣atuba municipality, Brazil (Shimozako et al., 2017).

year	Number of cases	Cumulative cases
1999	15	15
2000	12	27
2001	29	58
2002	52	110
2003	40	150
2004	41	191
2005	16	207
2006	20	227
2007	42	269
2008	27	296
2009	15	311
2010	4	315
2011	5	321
2012	6	327
2013	3	330
2014	12	342
2015	4	346

Fig. 3

Comparison of observed ZVL cumulative data from Arac̣atuba municipality, Brazil (dotted lines) and model prediction (solid curve). Parameter values used are as given in Table 4, with and the following initial conditions: (0) = 176000; (0) = 4000; (0) = 15; (0) = 9; (0) = 1000; (0) = 100; (0) = 50; (0) = 10; (0) = 1000; (0) = 2000; (0) = 300; (0) = 100; (0) = 10.

Table 4

Values and ranges of the parameters of the model (2.1).

Parameter	Range	Baseline	Reference
Π_H	4−7 day−1	6 day−1	(World Bank ata, 2015)
μ_H	3.67×10^−5−5.07×10^−5 day−1	3.67×10^−5 day−1	(World Bank ata, 2015)
b_R	0.03−0.2	0.16	(Hartemink et al., 2011, Shimozako et al., 2017)
β_H	0.2−0.8	0.56	(Hartemink et al., 2011, Zhao et al., 2016)
δ_H	2.37×10^−4−5.03×10^−4 day−1	0.0003 day−1	(Stauch et al., 2011, Zhao et al., 2016)
τ_H	0.12−0.95 day−1	0.5294 day−1	(Stauch et al., 2011, Zhao et al., 2016)
γ_H	0.00556−0.01667 day−1	0.0111 day−1	(Shimozako et al., 2017)
ψ_V	30−70 egg oviposition−1	50 egg oviposition−1	(European Centre for Disease Prevention and Control, )
K_M	9000−1.1×10^9	5.5×10^6	(Hussaini et al., 2016)
μ_E	0.05−0.25 day−1	0.143 day−1	(European Centre for Disease Prevention and Control, )
μ_L	0.0333−0.05 day−1	0.0455 day−1	(European Centre for Disease Prevention and Control, )
μ_P	0.0769−0.167 day−1	0.143 day−1	(European Centre for Disease Prevention and Control, )
μ_M	0.0416−0.083 day−1	0.0714 day−1	(European Centre for Disease Prevention and Control, )
σ_E	0.05−0.25 day−1	0.0833 day−1	(European Centre for Disease Prevention and Control,)
σ_L	0.0333−0.05 day−1	0.04 day−1	(European Centre for Disease Prevention and Control, 2017)
r_L	0.0009−0.011 day−1	0.00893 day−1	Fitted
σ_P	0.07−0.1667 day−1	0.0833 day−1	(European Centre for Disease Prevention and Control,)
f	0.413−0.9	0.5	Assumed
ξ_L	0.0456−0.564 day−1	0.1 day−1	(Poché et al., 2016)
ξ_M	0.0192−0.469 day−1	0.0923 day−1	(Poché et al., 2016)
β_V	0.2−0.8	0.7	(Hartemink et al., 2011, Zhao et al., 2016)
Π_R	7.49−11.4 day−1	8.33 day−1	(The Africa Report, 2013)
μ_R	1.522×10^−4−5.48×10^−4 day−1	2.28×10^−4 day−1	(PetCareRx, 2013)
b_R	0.03−0.2	0.16	(Hartemink et al., 2011, Shimozako et al., 2017)
δ_R	0.0099−0.0121 day−1	0.011 day−1	(Hartemink et al., 2011)
η_R	1.0−1.75	1.39	Fitted
ρ_R	0.00137−0.011 day−1	7.083×10^−3 day−1	Fitted
τ_R	0.01−0.04 day−1	0.0233 day−1	(Hartemink et al., 2011)
γ_R	3.9×10^−4−0.0167 day−1	0.011 day−1	(Parasites - Leishmaniasis, 2017, Parnell et al., 2008)

Basic properties

The basic properties of the model (2.1) will now be explored. It should be noted, first of all, that all parameters of the model are non-negative (with the death rates (, , , , ), recruitment rates , transmission probabilities and the biting rates assumed to be strictly positive). It is convenient to let . Consider the following equations for the rate of change of the total human, vector and reservoir host populations:

Furthermore, consider the region:

It can be shown (by solving for and in (2.2), (2.3), (2.4)) that all solutions of the system starting in the region remain in for all . Thus, the region is positively-invariant, and it is sufficient to consider solutions in . In this region, the usual existence, uniqueness and continuation results hold for the system (Forouzannia & Gumel, 2014).

Mathematical analysis

Disease-free equilibria

The model (2.1) has two disease-free equilibria, namely the trivial disease-free equilibrium (TDFE, denoted by ) and a non-trivial disease-free equilibrium (NDFE, denoted by ), as described below.where,with

TDFE (where no sandflies exist):

NDFE:

It follows that the NDFE exists if and only if . Furthermore, the NDFE reduces to the TDFE when . The threshold is similar to the vectorial reproduction number described in (Okuneye, Abdelrazec, & Gumel, 2018). It measures the average number of new adult female sandflies produced by one reproductive sandfly during its entire reproductive period. It is the product of the eggs oviposition rate , the fraction of eggs that survives and develops into larvae , the fraction of these larvae that survives and develops into pupae , the fraction of pupae that survives and develops into female adult sandflies and the average lifespan of adult female sandfly .

Asymptotic stability of TDFE

The TDFE of the model (2.1), , is globally-asymptotically stable (GAS) in whenever .

Let . Consider, first of all, the sandfly-only system of the model (2.1):

The system (3.2) has a unique trivial equilibrium (whenever ), given byin the invariant region

Furthermore, consider the following Lyapunov function for the system (3.2):where, and , with Lyapunov derivative given by (where a dot represents differentiation with respect to time t):

Thus, it follows, for in , that the Lyapunov derivative . Furthermore, it follows from the LaSalle’s Invariance Principle (Theorem 6.4 of (LaSalle, 1976)) that the maximal invariant set contained in , , , , is the singleton is the singleton . Hence, the unique trivial equilibrium of the system (3.2) is GAS in whenever . Thus, for ,

Since the model (2.1) is Type K (Smith, 1986), it follows, by substituting (3.3) into (2.1), that

Thus, by combining Equations (3.4), (3.3), it follows that the TDFE of the model (2.1) is GAS in whenever .

It is worth stating that the trivial equilibrium is ecologically unrealistic, since it is associated with the (unrealistic) scenario where sandflies do not exist.

Asymptotic stability of NDFE

Let (so that the NDFE, , of the model (2.1) exists). It can be shown, using the next generation operator method (Diekmann et al., 1990, Van den Driessche and Watmough, 2002), that the associated reproduction number of the model (2.1) (denoted by ) is given by:where,and,with , , , , , and . The result below follows from Theorem 2 of (Van den Driessche & Watmough, 2002).

The NDFE, , of the model (2.1), with , is locally-asymptotically stable (LAS) in if , and unstable if .

The epidemiological implication of Theorem 3.2 is that ZVL can be effectively controlled in the two hosts populations (humans and non-humans animal reservoir hosts) if the initial number infected hosts and vector are small enough (i.e., in the basin of attraction of the non-trivial disease-free equilibrium, ).

Interpretation of

The threshold quantity is ecologically and epidemiologically interpreted as follows.where,where is the total number of the cycles at which infectious reservoir received and failed treatment (and returned to the symptomatically-infectious class), accounting for the average number of new infectious sandflies generated by exposed (asymptomatically-infectious) reservoirs , measures the average number of new infectious sandflies generated by symptomatically-infectious reservoirs that have not undergone any treatment and accounts for the average number of new infectious sandflies generated by symptomatically-infectious reservoirs that have undergone (and failed) treatment (and return to the symptomatically-infectious class) at least once. In particular,

Interpretation of : The quantity , given in (3.6), is associated with the infection of susceptible reservoirs by infectious sandflies. It is the product of the infection rate of susceptible reservoirs by infectious sandflies and the average duration of infectious sandflies in the class, .

Interpretation of : The quantity , given in (3.7), is associated with the infection of susceptible sandflies by exposed (asymptotically infectious) and symptomatically infectious reservoirs. It can further be expressed as

is the product of the infection rate of susceptible sandflies by exposed (asymptomatically-infected) reservoirs , and the average duration in the class, .

is the product of the infection rate of susceptible sandflies by symptomatically-infected reservoirs , the probability that an exposed reservoir becomes symptomatic (i.e., survived the class and move to the class) , and the average duration in the class, .

is the product of the infection rate of susceptible sandflies by symptomatically infected reservoirs (described above in ii.), and the probability that such infectious reservoirs have received and failed treatment(s) at least once given by (where is the fraction of symptomatic reservoir hosts who received treatment (and progressed to the class), and is the fraction of reservoir hosts who failed treatment and reverts to the class). It is worth mentioning that the total number of the cycle at which infectious reservoir hosts received and failed treatments (and return to the symptomatically infectious class) is finite (i.e., ). Although ZVL is not completely curable (as relapses are common when treatment ceases), euthanasia1 is considered in some cases where the animal is chronically infected (and cannot be cured) (Petersen, 2009).

Backward bifurcation analysis

Backward bifurcation, which has been observed in numerous models for vector-borne diseases (see, for instance (Forouzannia and Gumel, 2014, Garba et al., 2008),), typically occurs when the asymptotically-stable disease-free equilibrium of the model co-exists with a stable endemic equilibrium when the associated reproduction number of the model is less than unity (Castillo-Chavez & Song, 2004). The epidemiological consequence of backward bifurcation is that having the associated basic reproduction number of the model to be less than unity, while necessary, is no longer sufficient for ZVL control (or elimination). In a backward bifurcation situation, effective community-wide control of ZVL (when ) is dependent on the initial sizes of the subpopulations of the model. In other words, backward bifurcation makes effective ZVL control in the community difficult. It is instructive, therefore, to explore the possibility of backward bifurcation in the model (2.1).

Let represents an arbitrary non-trivial equilibrium point (EEP) of the model (2.1) and,

Solving the equations of the model (2.1) at gives:where , and is as given in Equation (3.1).

For mathematical tractability, the computations will be carried out for the special case of the model (2.1) in the absence of disease-induced mortality in humans (i.e., ) and larval density-dependence (i.e., ). Let . It can be shown, by solving for the variables of the resulting reduced version of the model (2.1) at steady-state and (and simplifying), that the solutions of the resulting model (at steady-state, ) satisfy the following quadratic (in terms of ):where,with,

The results below follows from Equation (3.10).

Let . The model (2.1) with has:

a unique endemic equilibrium if ;

(ii)a unique endemic equilibrium if , and or ;

two endemic equilibria if , and ;

no endemic equilibrium otherwise.

Item (iii) of Theorem 3.3 suggests the possibility of a backward bifurcation in the model (2.1) (since the model could have two endemic equilibria when ). This is explored below.

The special case of the model (2.1) with undergoes a backward bifurcation at whenever the Inequality (A-3), given in Appendix A, holds.

The proof of Theorem 3.4, based on using Centre Manifold theory (Castillo-Chavez and Song, 2004, Forouzannia and Gumel, 2014), is given in Appendix A. Fig. 4 depicts the backward bifurcation diagram of model (2.1) for the cases with (Fig. 4(a)), (Fig. 4(b)) and and (Fig. 4(c)).

Fig. 4

Backward bifurcation diagrams of the model (2.1) in the absence of disease-induced death in humans (i.e., ). Parameter values used are as given by their baseline values in Table 4 with .

It is worth mentioning that, for a special case of the model with negligible disease-induced mortality in the host populations (such as ), the expressions for the backward bifurcation coefficients a and b given by Equations (A-2) and (A-3), respectively, in Appendix A, reduce to (it should be noted from Appendix A that eigenvectors , and are all positive, while and are negative):

Hence, it follows from Theorem 4.1 in (Castillo-Chavez & Song, 2004), that the special case of the model (2.1) with will not undergo a backward bifurcation at . This result is consistent with those reported for the dynamics of vector borne diseases, such as those in (Bowman et al., 2005, Forouzannia and Gumel, 2014, Garba et al., 2008, Hussaini et al., 2016). The global asymptotic stability of the non-trivial equilibrium of the model (2.1) is proved below for the aforementioned special case.

The NDFE, , of the special case of the model (2.1) with is GAS in whenever and .

The proof, based on the approach in (Dumont and Chiroleu, 2010, Okuneye and Gumel, 2017), is given in Appendix B. The epidemiological implication of Theorem 3.5 is that, for the special case of the model (2.1) with negligible disease-induced mortality in the host populations (i.e., ), bringing (and maintaining) the threshold quantity to a value less than unity is necessary and sufficient for the effective control (or elimination) of ZVL in the human and animal reservoir populations.

Sensitivity analysis and numerical simulations

The model (2.1) contains 30 parameters, and uncertainties in the estimates of these parameters are expected to arise. The effect of such uncertainties is assessed using Latin Hypercube Sampling (LHS) (Blower and Dowlatabadi, 1994, Marino et al., 2008, McLeod et al., 2006, Mckay et al., 1979). Furthermore, sensitivity analysis (using Partial Rank Correlation Coefficients (PRCC)) is carried out to determine the parameters that have the greatest influence on the dynamics of the disease (using the basic reproduction number as response function (as shown in Table 5)). The ranges and baseline values of the parameters tabulated in Table 4 will be used in this analysis.

Table 5

PRCC values for the parameters of the model (2.1) using the basic reproduction number as response function (the top three (most dominant) parameters that affect the dynamics of the model with respect to are highlighted in bold font). Parameter values and ranges used are as given in Table 4.

Parameters	PRCC_ℛ_0	Parameters	PRCC_ℛ_0	Parameters	PRCC_ℛ_0
Π_H	−0.15	ρ_R	+0.016	σ_P	+0.36
μ_H	+0.14	δ_R	−0.049	μ_P	−0.31
Π_R	+0.053	ψ_V	+0.42	f	+0.39
μ_R	−0.23	σ_E	+0.35	β_V	+0.37
η_R	+0.14	μ_E	−0.36	μ_MT	−0.91
b_R	+0.74	σ_L	+0.33	r_L	−0.061
β_R	+0.38	μ_LT	−0.28	K_M	−0.0016
γ_R	−0.51	τ_R	−0.0047

The top three PRCC-ranked parameters are the sandfly removal rate (given by the aggregated parameter , defined as ), the biting rate of sandflies on reservoir hosts , the progression rate of exposed reservoirs to active ZVL class . Furthermore, parameters such as the sandfly oviposition rate , fraction of female sandfly reaching adult stage , probabilities of infection per bite , progression rate of immature sandfly , death rates of immature sandfly are also influential (but not as dominant as the aforementioned top three PRCC-ranked parameters). Thus, this study shows that effective disease control entails a multi-faceted approach based on minimizing the contact reservoirs have with sandflies (i.e., minimizing and by clearing sandfly breeding sites around the reservoirs and spraying of sandfly repellents), reducing sandfly population (i.e., increasing and reducing by clearing sandfly breeding sites around the reservoirs) and early diagnosis of ZVL cases in reservoirs (i.e., increasing by ZVL screening to high-risk individuals).

The effect of the average lifespan (survival) of sandflies and animal reservoir hosts is monitored by simulating the model (2.1) using the baseline parameter values in Table 4 (relevant to ZVL dynamics in Arac̣atuba municipality, Brazil (Shimozako et al., 2017)). A contour plot of , as a function of and , shows that (i.e., disease burden) increases with increasing survival of both the vector and animal reservoir hosts (), as expected (Fig. 5). In particular, the range of values now increases to , with a mean of . It should be noted that these simulations were generated for the case when (since insecticide-based treatment strategy of the reservoir hosts was not implemented in the Arac̣atuba municipality during the 1999–2015 study period). However, in the hypothetical scenario where such treatment is used (and at the baseline rates given in Table 4), the range of significantly decreases to , with a mean of (Fig. 6). This represents about 60% reduction in the mean value of . Although the default scenario also suggests the feasibility of effective disease control (since the mean value of is ; and Theorem 3.5 shows that disease elimination is feasible, if and are small enough and ), the latter (hypothetical) scenario, where infected reservoirs are treated, significantly enhances the prospect of disease elimination in the municipality (since the mean value of is 0.35). This is quite intuitive, since the population of sandflies obviously decreases if the larvae and adult sandflies continue to feed from the faeces of infected reservoirs. But this poses an ecological dilemma, since treatment of reservoirs can lead to perhaps the removal of sandflies from the local ecosystem (albeit it serves a major epidemiological function of minimizing, or even eliminating, ZVL burden in the community). These simulations show that ZVL modeling studies in communities where such insecticide-based treatment strategy of infected reservoirs is implemented may well be over-estimating the disease burden if they failed to explicitly incorporate the effect of such treatment (i.e., additional larval and adult sandfly mortality due to their feeding on the faeces of the treated infected reservoirs) in the model. Fig. 5 further shows that is more sensitive to increases in sandfly lifespan than that of the animal reservoir (so, a strategy that focuses on reducing sandflies, rather than the animal reservoir (e.g., via culling), may be more effective in reducing ZVL burden in the community).

Fig. 5

Contour plot of , as a function of the average life expectancy of sandflies and animal reservoir hosts . Parameter values used are as in Table 4.

Fig. 6

Contour plot of , as a function of the average life expectancy of sandflies and animal reservoirs . Parameter values used are as in Table 4.

Conclusions

This study is based on the design, analysis and numerical simulations of a new deterministic model for assessing the transmission dynamics of zoonotic visceral leishmaniasis (ZVL) in a community. The model is fitted using case and demographic data relevant to ZVL dynamics in Arac̣atuba municipality in Brazil. The main theoretical and epidemiological findings of the study are summarized below.

The model has a trivial disease-free equilibrium (TDFE) which is globally-asymptotically stable if a certain vectorial threshold quantity is less than unity. It also has a non-trivial disease-free equilibrium (NDFE; whenever is greater than unity) which undergoes a backward bifurcation under certain conditions. In the absence of backward bifurcation, the NDFE is globally-asymptotically stable, for a special case whenever the associated reproduction number is less than unity.

Sensitivity analysis of the model (using the basic reproduction number as the response function) show that the top three PRCC-ranked parameters are the sandfly removal rate ( and ), the biting rate of sandflies on reservoir hosts , the progression rate of exposed reservoirs to active ZVL class . Hence, this study identifies the parameters that should be targeted for effective anti-ZVL control in the community. Other parameters with high PRCC ranking (but not as high as the aforementioned three) are sandfly oviposition rate , the fraction of pupae that became adult female sandflies , and the infection probabilities and ).

Numerical simulations, using the data for ZVL dynamics in Arac̣atuba municipality during the 1999–2015 study period, show that the associated reproduction number ranges from 0.3 to 1.4, with a mean of 0.85. This range dramatically decreases, to (with a mean of 0.35), when insecticide-based treatment of the animal reservoir hosts is implemented. Thus, the prospect of the effective control of ZVL in the community is greatly enhanced if a control strategy based on using insecticide-based treatment of the animal reservoir is implemented. Furthermore, ZVL modeling studies in communities where such treatment is used may be over-estimating the disease burden if they fail to explicitly incorporate the effect of such treatment (i.e., resulting in additional larval and adult sandfly mortality) in the model formulation.

The reproduction number is more sensitive to increases in sandfly lifespan than that of the animal reservoir (so, a strategy that focuses on reducing sandflies, rather than the animal reservoir (e.g., via culling), may be more effective in reducing ZVL burden in the community).