Mathematical modelling for Zoonotic Visceral Leishmaniasis dynamics: A new analysis considering updated parameters and notified human Brazilian data.

Brazil is one of the highest endemic countries for Zoonotic Visceral Leishmaniasis: according to the Brazilian Ministry of Health, the annual number of new human cases and deaths due to this disease has been increasing for the last 20 years. In addition, regarding the Americas, the specific relationship between canine and human for Visceral Leishmaniasis dynamics is still not well understood. In this work we propose a new model for Zoonotic Visceral Leishmaniasis, based on the models previously published by Burattini et al. (1998) and Ribas et al. (2013). Herein, we modeled the disease dynamics using a modified set of differential equations from those two authors, considering the same assumptions (inclusion of human, dog and sandfly populations, all constants over time). From this set of equations we were able to calculate the basic reproduction number R0 and to analyze the stability and sensitivity of the system to the parameters variability. As main result, when the stability of the system is reached, the normalized reporting human cases rate is estimated in 9.12E-08/day. This estimation is very close to the 2015 report from Araçatuba city, 5.69E-08/day. We also observed from stability and sensitivity analysis that the activity of sandfly population is critical to introduction and maintenance of Zoonotic Visceral Leishmaniasis in the population. In addition, the importance of dog as source of infection concentrates on latent dog, since it does not show clinical symptoms and signs and, therefore, has a great contribution to disease dissemination. As conclusion, considering the presently ethical issues regarding to elimination of positive dog in Brazil and the highly sensitivity of disease dynamics on sandfly population, we recommend that the sandfly population control should be prioritized.

Introduction

Zoonotic Visceral Leishmaniasis is one of the world deadliest and neglected infectious diseases, according to World Health Organization. This disease is endemic in 80 countries worldwide, in which 90% of all cases occur in Bangladesh, Brazil, India, Nepal and Sudan. Thus, about 360 million of people are exposed to risk of infection in the world (Duthie et al., 2012, Killick-Kendrick, 2010, Pan American Health Organization, 2001, World Health Organization, 2017). The Zoonotic Visceral Leishmaniasis is a disease of major human and veterinary medical significance that involves a complex interplay between trypanosomatids protozoan from Leishmania complex, arthropod vectors (in Brazil, we find the female sandfly Lutzomyia longipalpis and Lutzomyia cruzi), environmental influence on vector distribution, small companion animal (dog) reservoir of infection and susceptible human populations. In American continent, Leishmania infantum chagasi is the most important specie from Leishmania complex.

From the last few years, Zoonotic Visceral Leishmaniasis has been emerging within non-endemic areas, mostly because of transportation of dogs from endemic areas and climatic changes with the expansion of the geographical range of the sandfly vector. Thus, the effective control will essentially involve interdisciplinary teams of microbiologists, parasitologists, entomologists, ecologists, epidemiologists, immunologists, veterinarians, public health officers and human physicians (Palatnik-de-Souza & Day, 2011).

Besides the publication of guidelines of Zoonotic Visceral Leishmaniasis control and the investments made in general surveillance activities, the sandfly and the reservoir in urban areas remains among the major challenges for the control activities. These challenges are due to (1) the necessity to better understand the vector behavior in urban environment; (2) the operational and logistic difficulties to carry out activities in sufficient time to obtain good results; and (3) the high costs involved in these activities (Killick-Kendrick, 2010, Maia-Elkhoury et al., 2008). In addition, regarding the Americas, the specific relationship between canine and human for Visceral Leishmaniasis dynamics is still not well understood. Thus, the control of the animal reservoir is complex and often needs to combine different ways of interventions. In particular, the Brazilian Control Program recommends a strategy based on canine culling and vector control with insecticide spraying (Ministry of Health, 2006, Nunes et al., 2008). Therefore, dog treatment is not recommended, since it is difficult to eliminate the parasitemia from infected dogs (Athanasiou et al., 2013, Ministry of Health, 2006). Furthermore, insecticide-impregnated collars for dogs and canine vaccination are not currently recommended as public health control measures (Palatnik-de-Souza and Day, 2011, Romero and Boelaert, 2010).

In this work we propose a new model for Zoonotic Visceral Leishmaniasis, based on the models previously published by Burattini, Coutinho, Lopez, and Massad (1998) and Ribas, Zaher, Shimozako, and Massad (2013). In this new model we updated most of parameters, calculated the new value and analyzed the stability and sensitivity of the system. Then, we discussed the disease dynamics based on those mathematical analyses and addressed the critical points that benefit the introduction and maintenance of this disease in the population.

The model

We used a mathematical model that is an adaptation of the one proposed by Burattini et al. (1998). In our model, we assume:

A human and a dog population, with the biological vector transmitting the infection within and between the two populations;

Those three populations (humans, dogs, and vectors) are constants;

Both human (indexed as h) and dog (indexed as d) populations are divided into four categories: susceptible (x and x), infected but without noticeable disease (l and l) (i.e., “latent”), clinically ill (y and y), and recovered immunes (z and z). On the other hand, the vector population is divided into three categories: noninfected, infected but not infective, and infective individuals, denoted as s, s, and s, respectively.

The flowchart and compartment model (Fig. 1) and the set of differential equations describing the model's dynamics (System 1) are presented as following.

Fig. 1

The compartment model and the flowchart. Note that only dogs are source of infection and sandflies transmits the Leishmania sp. to both, dogs and humans.

The definition, biological meaning, and values of each of parameter are described in Table 1.

Table 1

Parameters adopted in our model. The indexes h, d and s stand for humans, dogs and sandflies, respectively.

Parameter	Meaning	Value	Dimension	Source
μh	Natural mortality rate	3.67 × 10−5	1/day	Brazilian Institute of Geography and Statistics, Brazil (2013)
αh	Kalazar specific lethality	6.31 × 10−3	1/day	World Health Organization (2017)
ah	Average daily biting rate	2.00 × 10−1	human/(sandfly × day)	Epidemiological Surveillance Direction, Santa Catarina State, Brazil (2008)
mh(t)	Vector density per host (time-dependent)	Variable	sandfly/human	Fitted
whc	Ratio human:house	3	human/house	Brazilian Institute of Geography and Statistics, Brazil (2013)
bh	Proportion of infective bites	1.00 × 10−2	dimensionless	Molineaux and Gramiccia (1980)
rh	Spontaneous recovery rate	5.48 × 10−4	1/day	Badaró et al. (1986)
γh	Loss of immunity rate	5.48 × 10−4	1/day	Kault and March (1991)
δh	Latent recovery rate	1.10 × 10−2	1/day	Bardaró et al. (1986)
φh	Inverse of incubation period	4.00 × 10−4	1/day	Pearson and Souza (1990)
σh	Recovery rate to immunes	2.50 × 10−3	1/day	Ministry of Health, Brazil (2006)
ηh	Proportion of unreported cases	0.705	dimensionless	Maia-Elkhoury, Carmo, Sousa-Gomes, and Mota (2007)
μd	Natural mortality rate	2.28 × 10−4	1/day	http://www.sciencedirect.com/science/article/pii/S0960982213004132Selman, Nussey, and Monaghan (2013)
αd	Kalazar specific lethality	1.81 × 10−3	1/day	Lanotte, Rioux, Perieres, and Vollhardt (1979)
ad	Average daily biting rate	2.00 × 10−1	dog/(sandfly × day)	Epidemiological Surveillance Direction, Santa Catarina State, Brazil (2008)
wdh	Ratio human:dog for Araçatuba/SP city	10/1.8	human/dog	Andrade, Queiroz, Perri, and Nunes (2008)
md(t)	Vector density per host	wdh× mh(t)	sandfly/dog	–
φd	Inverse of incubation period	3.78 × 10−4	1/day	Greene (2011)
bd	Proportion of infective bites	1.00 × 10−2	dimensionless	Molineaux and Gramiccia (1980)
rd	Spontaneous recovery rate	2.74 × 10−4	1/day	Lanotte et al. (1979)
γd	Loss of immunity rate (recovery to susceptible)	2.74 × 10−3	1/day	Kault and Marsh (1991)
σd	Recovery rate from clinically ill to immunes	9.04 × 10−4	1/day	Lanotte et al. (1979)
δd	Latent recovery rate	8.22 × 10−3	1/day	Lanotte et al. (1979)
ξd	Dog elimination rate	3.36 × 10−4	1/day	Camargo-Neves (2004)
μs	Natural mortality rate	5.00 × 10−2	1/day	Ministry of Health, Brazil (2006)
τ	Extrinsic incubation period	7	day	Neva and Sacks (1990)
aS	Average daily biting rate (on dogs)	2.00 × 10−1	1/day	Estimated as Epidemiological Surveillance Direction, Santa Catarina State, Brazil (2008)
cl	Probability of latent dog to infect the sandfly	0.385	dimensionless	Laurenti et al. (2013)
cy	Probability of clinically ill dog to infect the sandfly	0.247	dimensionless	Laurenti et al. (2013)

A brief description of system (1) should clarify their meaning.

Let S be the total number of sandflies. The number of bites inflicted in the human host population in an infinitesimal time interval is , where a is the biting rate on humans. The number of bites inflicted by infected flies is , where is the number of infected flies.

Let now be the total number of susceptible individuals in the human population. In an infinitesimal time interval , varies as follows:

The infected flies are able to bite on any category of human population. Thus, only a fraction of the infected bites are on uninfected individuals: , where is the fraction of uninfected humans. But, a fraction b of becomes latent, so X diminishes by ;

Simultaneously, individuals, latent and immune, revert to the susceptible condition, and die by natural causes other than the disease.

We must add an entrance term, due to natality, which we choose to be , where is the disease-induced mortality rate, is the number of infected humans (clinically ill humans), and is the total number of humans needed to maintain a constant population (where , with as the number of latent humans and as the number of recovered humans).

Thus we have:

Dividing this equation by and calling , we get the first equation of System (1).

Observe that is a function time-dependent: . This expression is the simplest way to simulate the changings on sandfly population size dynamics between 1999 and 2015.

We can apply the same process in order to obtain the equation for the dynamic of susceptible dogs (). However, observe from Table 1 that the ratio sandfly:dog depends on the ratio sandfly:human and on the ratio human:dog: . Although all the populations are constant, if we consider the real number of individuals, we expect more humans than dogs. Thus, if the sandfly population is constant, we have different values for and .

The last three equations of system 1 refer to the flies. When infected, a fly remains in a latent stage for a period of time . This time corresponds to the extrinsic incubation period of the parasite inside the vector fly. Numerically it lasts for about half the life expectancy of the flies.

Let S be the number of susceptible flies. In an infinitesimal period of time , bites due to uninfected flies occur on latent and infected dogs (humans are not considered to be infective for flies; see Tesh (1995)). A fraction and of the flies (who bites latents and clinically ill dogs, respectively) becomes latently infected as a result. Therefore, we have:

Dividing by and by , we get the equation for non-infected sandflies ().

Basically, we adopted different mathematical techniques for dogs and humans in latency stage and for sandflies in latency stage because there are differences between their biological characteristics. In the case of human and dogs, once they are in the incubation period, this latency time usually presents an exponential distribution (in average approach). That's why we chose modelling this incubation process using a latent compartment instead of a delay term. On the other hand, the incubation process regarding to sandflies is biologically different from humans and dogs. When a sandfly is infected with leishmania parasite, this sandfly becomes infective only after a constant latent period (Bocharov & Rihan, 2000). Therefore, in this case, it is much more feasible using delay term to model sandfly population dynamics instead of latent compartments. In addition, the infected dogs get infective in an exponential fasion whereas the sandflies get infective immediately after the extrinsic incubation period elapses.

Although this is a brief but detailed description about the non-infected categories equations (that is, , and ), we can note that each term of our system equation has a biological meaning. The meaning of each term is in agreement with the parameter in which is together with.

Although Zoonotic Visceral Leishmaniais (ZVL) is a disease that causes immunological damages (and, therefore potencialize the probability of co-infections), those co-infections were possible and lethal only because of primary leishmaniasis infection. Thus, we considered that ZVL - infected individuals die because of the ZVL infection.

The positivity and boundedness of the solutions

We argue about the positivity and boundeness of the solutions considering the proof provided by Burattini et al. (1998), which our model was based on.

An important point about system (1) is that positivity is preserved. Thus, given positive initial conditions, the variables remain non-negative.

Firstly, let us consider the equations for sandfly population dynamics (that is, the equations for , and ). Observe that the term is always positive in those three equations. However, it may appear at first sight that positivity is not preserved by system (1), because if in the equation for , the delayed term could be non-zero. But, we can demonstrate that if becomes zero, then this delayed term (which is a term of removal from the compartment ) must necessarily be zero. Also, observe that the equations for susceptible category (in both human and dog populations) can be rewrite as (4):

As , those susceptible categories can never become negative. As a consequence, the other categories , and can never become negative. In addition if for some it will never become negative afterwards. Note also that if for some it Will also never become negative since .

Let now be the proportion of flies which become infected between and and have survived for the period . This expression is positive and is a fraction of the total latent flies.

Note that the proportion of latent flies is the sum of in , from to , that is:

Therefore, if then for all . Now, replacing by in , we have that if , then .

Once the positivity of the solutions was proven, we can discuss about the boundedness of the solutions. Here in, we will consider a similar approach as that one presented by Cai, Lashari, Jung, Okosun, and Seo (2013). Thus, it can be shown that the region given by (6):is positively invariant with respect to system (1). Thus, every solution of (1), with initial conditions in remains there for Therefore, it is sufficient to consider the dynamics of the flow generated by (1) in . In this region the model can be considered as been epidemiologically and mathematically well posed.

The reported cases

In Brazil, Zoonotic Visceral Leishmaniasis is a notifiable disease (Ministry of Health, 2006, Day et al., 2012). Thus, we can assume:

A infected human should look for medical treatment when he/she will become clinically ill ();

Only a fraction of those humans that are clinically ill will be reported to sanitary authorities. The remaining fraction (I) will not look for medical help, even if the clinical symptoms and signs appear; or (II) will not be correctly reported in the hospitals.

Now, let's see again the equation for in system (1):

The term from (7) means the rate of latent humans who become clinically ill per day. Thus, per day, those amounts of humans are eligible to look for medical help. However, only a fraction of those clinically ill humans will be correctly notified to sanitary authorities. Therefore, the daily rate of reported human cases is defined by equation (8):

The Centre of Epidemiological Surveillance of São Paulo State (CES-SP) (2016) is the institution who administrates the data about ZoonoticVisceral Leishmaniasis in São Paulo State. In order to validate our model, we decided to use the data of human reported cases from the municipality of Araçatuba (São Paulo State – Brazil) as reference, because it is an endemic city for this disease. Those data are presented in Table 2 and are available on Centre of Epidemiological Surveillance of São Paulo State website (Centre of Epidemiological Surveillance of São Paulo State, Brazil, 2016).

Table 2

Human reported cases in Araçatuba municipality: average of the normalized rate per day for each year.

Year	Human reported cases per year (CES-SP)	Araçatuba's Human population size (BIGS)	Average of normalized human reported cases per day
1999	15	169303	2.43E-07
2000	12	170296	1.93E-07
2001	29	171289	4.64E-07
2002	52	172768	8.25E-07
2003	40	174399	6.28E-07
2004	41	177823	6.32E-07
2005	16	179717	2.44E-07
2006	20	181598	3.02E-07
2007	42	181371	6.34E-07
2008	27	181143	4.08E-07
2009	15	182204	2.26E-07
2010	4	182365	6.01E-08
2011	5	182526	7.51E-08
2012	6	183441	8.96E-08
2013	3	190536	4.31E-08
2014	12	191662	1.72E-07
2015	4	192757	5.69E-08

Note that we have the total of reported cases per year. Thus, since our time scale is day, we estimated an average of human reported cases per day for each year (dividing the total from each year by 365). Finally, we also have to consider that our model works with normalized population (all three populations are constant). Thus, as a last step, we have to divide each rate of human reported cases per day by the official population size of Araçatuba municipality. The population size of Araçatuba municipality is available on Brazilian Institute of Geography and Statistics (BIGS) website (Brazilian Institute of Geography and Statistics, Brazil, 2016).

In order to fit and compare our results to real data, we also calculated a normalized average of reported cases per day from each 365 days of simulation. This simulation was run considering 60 years and the obtained curve was fitted by simple handling along the time-axis (for instance, we could assume the initial day as the first day of 1970 or 1980, depending on how best the simulated curve fits on the real data). Thus, we could obtain the yearly average of reported human cases per day and compare it to the real yearly average provided by CES-SP (Table 2).

The basic reproduction number ()

According Anderson and May (2010), the Disease Free Equilibrium (DFE) state means the state in which there is no disease in the population. Also, there are no infected individuals, even latent ones.

The basic reproduction number, denoted , is ‘the expected number of secondary cases produced, in a completely susceptible population, by a typical infective individual’. If , then on average an infected individual produces less than one new infected individual over the course of its infectious period, and the infection cannot grow. Conversely, if , then each infected individual produces, on average, more than one new infection, and the disease can invade the population (van den Driessche & Watmough, 2002).

In order to estimate the of our model, we applied the method presented by van den Driessche and Watmough (2002). Here, we will avoid any mathematical demonstration regarding to Driessche & Watmough method, since this is not the focus of this work. However, the idea is to calculate the relationship between the parameters that causes instability to the trivial solution of system (1), which represents the absence of disease in the populations considered. That is, we studied the stability of the solution , , , (where ) and , , . If this solution turns out to be stable, that is, when , the disease cannot invade the population (Burattini et al., 1998). We strongly suggest the reader to see all the details in their publication (see van den Driessche and Watmough (2002) in the references list).

Herein, we remark that the Driessche & Watmough method was not developed for delay differential equations. Thus, we refer to Burattini et al. (1998) and Wei (2004) to explain about the delayed terms. Since this system contains a time delay in the population of flies, the linearization around the trivial solution is not straightforward. Basically, the analysis is a particular case of stability analysis on free-disease state condition. In the case of delayed term in our system, we need to consider two Jacobian Matrices separately: one is the usual Jacobian (for term without delay) and the other one is the Jacobian for delayed terms only. We suggest the reader to see Burattini et al. (1998) and Wei (2004), in order to understand how to handle delayed differential equations for analysis.

We will suppress all steps of the calculation and present the final equation for . Thus, we get (Equation (9)):

It is possible to observe that equation (9) splits naturally into two terms: the contributions to from latent dog population and clinically ill dog population. Thus, we have system (10):

Let's explain biologically the meaning of in a similar approach as that one presented by Burattini et al. (1998). Given a population in a DFE state, someone could ask if the introduction of a small amount of infective individuals would start an epidemic outbreak. Once this epidemic appears, it eventually would converge to an endemic equilibrium state (that is, equilibrium with the disease). In this case, if , even a small amount of infective individuals would start an epidemic, which it would be in a endemic level different from zero. However, observe that the expression from equation (9) includes 2 components, e (system (10)). As an example, once a small amount of latent dogs () is introduced, if the epidemic would be installed due to the introduction of those small amount of latent dogs. The same idea occurs if we consider clinically ill dogs (), since . And, if the introduced individuals are from different subpopulations, the epidemic would be installed if the sum of the all contributions is more than 1. In other words, we should have .

Fitting the ratio Human:Sandflies (m(t)) and model dynamics

Previously, we demonstrated by Equation (9) that depends on the parameter (dog/sandfly ratio). But, we can assume that dog population size is related to human population's habits and culture (Beck, 1973; Molineaux & Gramiccia, 1980). Therefore, we can estimate the human:sandfly ratio if we consider the human:dog ratio () for the municipality of Araçatuba. According to Andrade et al. (2008) this ratio was estimated as human/dog. As consequence, we have .

Once we need to fit the value for our model and considering we are interesting to understand the disease dynamic, we can consider the relation between and according to Equation (9). In this case we have to suppose . Thus, we obtain:where numerically we need sandfly/human.

The real data provided in Table 2 suggests that the incidence was not constant along those years in which the data was collected (1999–2015). One reasonable hypothesis is the climate changings that have been occurring for the last years (Massad, Coutinho, Lopez, & da Silva, 2011). Thus, since the sandfly population dynamics depends on climate and geographical conditions, we can include this idea in our model by fitting as time-function. It is not the scope of this paper to model the sandfly population dynamics according to climatic and geographic variations. Therefore, we will assume that a simple function for that can fit the simulation data to the real data, should include those climatic and geographic variabilities.

Let's consider the following function for :

The parameter values for (12) are in Table 3. Biologically, we can suppose that sandfly population reaches stability and oscillations decrease overtime. Thus, note that for t → + ∞ we have trending to .

Table 3

Parameter values for (12) and their biological meaning.

Parameter	Meaning	Value	Dimension	Source
mh0	Vector density per host (baseline value)	0.75	sandfly/human	Fitted
A	Vector density per host	3.4	sandfly/human	Fitted
B	Vector density per host	8.3	sandfly/human	Fitted
L	Linear constant	3.0	dimensionless	Fitted
K1	Proportionality constant	3.5 × 365	day	Fitted
T	Sandfly population dynamics period	5.5 × 365	day	Fitted

Although the ratio shows an oscillating behavior over time, we remind that is the ratio sandfly/human. In other words, in our model, the absolute size of sandfly population is time dependent. However, we normalized our system. Once the population is normalized, we can assume (that is, the sum of proportions of each category in the sandfly population is always equals to 1). Thus, in a proportional approach, all three populations (humans, dogs, and vectors) can be considered constants.

Let's consider the equilibrium condition, that is, . If we substitute the parameters that compose the expressions in (9) and (10), we obtain and . Therefore, the sum of those two values provides us the total contribution from those two classes of dogs, . The difference between and values could be explained by the skin integrity of the infected dogs. In other words, the clinically ill dogs () usually present skin lesions and the skin of a dog from this category is more damaged than that one from a latent dog (). Because of this, we can suppose that the sandflies are less probable to acquire available parasites from dogs of category. On the other hand, the opposite occurs with the latent dog, since their skins are heathier than the clinically ill dog's skin (Laurenti et al., 2013).

Stability analysis

Mathematically, our model is a nonlinear delay differential equation system. It is very usual to model the dynamics of natural phenomena using nonlinear systems, because most of them are ruled by nonlinear behavior. However, in opposition to linear systems, the dynamics of nonlinear systems commonly are not simple and they may appear chaotic. Thus, because of the behavior of nonlinear differential systems, it is useful to study the stability of this system. Herein, we follow the method describe by Wei (2004).

A nonlinear system is considered stable when the variable's derivatives are zero (, is the vector of variables). When the stability is reached, the variables assume constant values, and they are the fixed points of the system. In order to determine how stable the system is when it reaches the fixed points, we need to obtain the Jacobian Matrix (J, as the usual Jacobian Matrix, and J, for the time-delayed terms) of the system and calculate its values on the fixed points. Following, we calculate the eigenvalues λ of the determinant below (Equation (13)).where I is the identity matrix. If all eigenvalues λ has negative real part, the equilibrium of the system at the fixed point is stable. On the other hand, if there is at least one eigenvalue λ with positive real part, the system at the fixed point is unstable.

Basically, our model is composed by three populations, but the human population dynamics is directly dependent on the sandfly population dynamics. On the other hand, the sandfly and dog populations are mutually dependent. Since the human population does not interfere on dog or sandfly dynamics, we do not need to include the humans' equations in the fixed point calculation (the humans’ population fixed points will naturally be solved once we obtain fixed point). Thus, we would work with seven equations only.

We also have to consider that the populations are constants: and . If we use those two conditions, we are allowed to eliminate one differential equation of each population by substituting one category of each population by the respective condition. For convenience, we adopted and . Thus, we obtain the following system (14), reduced to five equations.where the terms with the index τ are the time-delayed terms. The meaning of each parameter is in Table 4.

Table 4

Parameter meanings for (14). All parameters are real positive values.

Parameter	Meaning
Ξ	μd + ξd + γd
β	bdadmd
C1	ascl
C2	ascl exp(-μsτ)
D1	as cy
D2	ascy exp(-μsτ)
E	γd- rd
F	γd- αd
G	μd + ξd + rd + δd + φd
H	μd + ξd + αd + σd

Before obtain the Jacobian Matrices J and J, we need to linearize system (14) around the fixed points, applying Taylor series for differential equation systems (Fiedler-Ferrara & Prado, 1995). Thus, considering the expansion until the first order terms, we have system (15).where the terms with index refer to time-delay terms, the star ‘*’ refers to fixed points and the tilde ‘∼’ indicates the first order approximation for the distances between the variable's value and the fixed points, for instance . Therefore, the tilde marked variables describe the local behavior of solutions close to fixed point and it help us to understand how the system progresses when the initial conditions (in this case, we suppose the trivial solution as initial conditions) are lightly disturbed from the equilibrium state (Fiedler-Ferrara & Prado, 1995).

Numerically, once reached the equilibrium state, the time-delayed terms have the same value as the usual terms. That is, .

The following steps depend on which fixed point we are evaluating. Let us start from the Disease Free Equilibrium (DFE) state. Following, we analyze the Endemic Equilibrium (EE) state.

Stability of Disease Free Equilibrium (DFE)

Our model is considered in DFE state if we consider the trivial solution for fixed points: x = s = 1 and l = y = s = 0. If we substitute those fixed points on system (15), we can uncoupled the equation for and from the remained system, since those two variables disappear on the other three equations. Therefore, the system we need to analyze is (16).

And, applying (13) on system (16), we obtain (17).

equation (17) is the characteristic equation of fixed points of the system. This equation is very similar to usual polynomial equations, exception to exponential terms . This kind of equation is classified as quasi-polynomials and, in opposition to polynomial equations, they usually have infinite solutions in the complex plane. Because of this natural difficult to handle quasi-polynomial equations, we adopted the approximation . Once substituting the exponential terms by this approximation, we obtain a third order polynomial equation. Finally, using the numerical values from Table 1 on Table 4, we were able to calculate the eigenvalues λ from (17): . Since we had at least one eigenvalue greater than 0, we conclude that when the system is on DFE state, the equilibrium is unstable.

Stability of endemic equilibria (EE) and backward bifurcation

The stability analysis of the Endemic Equilibrium (EE) is exactly the same process. However, in this case, all fixed points are non-zero. Therefore, we are not able to simplify the system as we did before, because the all five equations will be coupled among them. For this analysis, we have to consider the non-trivial fixed points: and . Thus, we have to analyze system (10) with their five equations.

Applying (13) on (15), we obtain (18).

Considering the approximation , we found the following values for eigenvalue λ: We observed all eigenvalues have negative values. Therefore, the Endemic Equilibrium state of this system is stable.

The mathematical condition to observe the backward bifurcation in our model (when ) is the existence of two positive real solutions. Thus, when the equilibrium points were calculated considering System 14 and Table 4, we obtain the following s solutions (we suppressed the full calculations):

From Table 1, we have all parameter values and all of them are real positive values. Therefore, for we will always have a negative solution, since there is a minus signal in front of the expression. On the other hand, for we found that the positivity depends on numerical values for numerator, since the denominator is already positive (from we observed in denominator). Therefore, the existence of should obey the following condition:

In our model, we have only one physical solution for s3. Thus, there is no occurrence of backward bifurcation in our system.

Sensitivity analysis

The precise of the results of mathematical and computational models of biological systems is directly dependent of how certain the parameters are. Such parameters are usually estimated from experimental approaches. In some situations, we have some parameters subject to uncertainty due to the lack of complete information about their source (Adhikari & Supakankunti, 2010). Thus, the presence of uncertainty in the experimental data may lead to uncertainties on the estimated parameters. Consequently, the uncertain parameters can propagate their uncertainties onto mathematical models’ results (Vuolo, 1996). Even when a parsimonious approach is followed during model building, available knowledge of phenomena is often incomplete, and experimental measures are lacking, ambiguous, or contradictory. So the question of how to address uncertainties naturally arises as part of the process. Uncertainty and sensitivity analysis techniques help to assess and control these uncertainties. Uncertainty analysis is performed to investigate the uncertainty in the model output that is generated from uncertainty in parameter inputs. Sensitivity analysis naturally follows uncertainty analysis as it assesses how variations in model outputs can be apportioned, qualitatively or quantitatively, to different input sources (Marino, Hogue, Ray, & Kirschner, 2008).

In this study, we used the Latin Hypercube Sampling (LHS) as uncertainty analysis and Partial Rank Correlation Coefficient (PRCC) as index for sensitivity analysis. It is not the focus of this paper to show the mathematical demonstration of those techniques. So, we strongly recommend to the reader to see Marino et al. (2008) for more information.

We analyzed the parameters sensitivity in the nontrivial equilibrium (Endemic Equilibrium) state. As we presented in (7), when we have trending to . Thus, we are able to calculate the equilibrium points. Once we obtain the equilibrium points expressions, we can use them to evaluate their sensitivity to the parameters. We evaluated the sensitivity considering a range of ±1% of each parameter value.

In Table 5 we present a summary of the sensitivity analysis. Fig. 2, Fig. 3, Fig. 4, Fig. 5 we graphically illustrate the relation between variable sensitivity and parameters.

Table 5

Equilibrium state values and the set of parameters that are most sensitive for each variable in the non-trivial equilibrium state.

	Disease Free Equilibrium State (Trivial Equilibrium Point)	Disease Equilibrium State (Nontrivial Equilibrium Point)	Parameter set related to human population, in which the variable is sensitive (on EE state)	Parameter set related to dog population, in which the variable is sensitive (on EE state)	Parameter set related to sandfly population, in which the variable is sensitive (on EE state)
xh	1.0	9.84E-01	-	wdh, δd, φd, σd, ξd	bd, ad, as,mh0, μs, τ, cl
lh	0.0	7.73E-04	φh	wdh, δd, φd, σd, ξd	bd, ad, as,mh0, μs, τ, cl
yh	0.0	3.50E-05	-	wdh, δd, φd, σd, ξd	bd, ad, as,mh0, μs, τ, cl
zh	0.0	1.47E-02	φh	wdh, δd, φd, σd, ξd	bd, ad, as,mh0, μs, τ, cl
xd	1.0	9.80E-01	φh	wdh, δd, φd, σd, ξd	bd, ad, as,mh0, μs, τ, cl, cy
ld	0.0	5.43E-03	φh	wdh, δd, φd, σd, ξd	bd, ad, as,mh0, μs, τ, cl, cy
yd	0.0	6.26E-04	φh	wdh, δd, φd, σd, ξd	bd, ad, as,mh0, μs, τ, cl, cy
zd	0.0	1.37E-02	φh	wdh, δd, φd, σd, ξd	bd, ad, as,mh0, μs, τ, cl, cy
s1	1.0	9.91E-01	φh	wdh, δd, φd, σd, ξd	bd, ad, as,mh0, μs, τ, cl, cy
s2	0.0	2.63E-03	φh	wdh, δd, φd, σd, ξd	bd, ad, as,mh0, μs, τ, cl, cy
s3	0.0	6.27E-03	φh	wdh, δd, φd, σd, ξd	bd, ad, as,mh0, μs, τ, cl, cy
Rep [1/day]	0.0	9.12E-08	μh	wdh, δd, φd, σd, ξd	bd, ad, as,mh0, μs, τ, cl
R0dl	-	9.58E-01	φh	wdh, δd, ξd	bd, ad, as,mh0, μs, τ, cl
R0dy	-	7.09E-02	-	μd, αd, wdh, δd, φd, σd, ξd	bh, bd, as,ad, mh0, μs, τ, cy
R_0	-	1.03E+00	φh	wdh, δd, φd, σd, ξd	bd, ad, as,mh0, μs, τ, cl, cy

Fig. 2

PRCC values in respect to human population categories. Parameters that are significant, (p < 0.05) are marked with a star.

Fig. 3

PRCC values in respect to dog population categories. Parameters that are significant, (p < 0.05) are marked with a star.

Fig. 4

PRCC values in respect to sandfly population categories. Parameters that are significant, (p < 0.05) are marked with a star.

Fig. 5

PRCC values in respect to reported cases rate and . Parameters that are significant, (p < 0.05) are marked with a star.

From this sensitivity analysis method, we obtained some interesting results. This method allows us to check if some variable is sensitive for any parameter of the system, even if this parameter is not directly related to a specific population. This characteristic is different from that one presented by Burattini et al. (1998), in which was found a relationship between the parameter and variable are from the same population. However, we have to stress that Burattini et al. (1998) did not used the same sensitivity analysis method as here.

From Table 5, we observed that all variables listed in this table are sensitive for most of parameters related to sandfly population dynamics. In particular, the parameters that compose the force of infection - , and - resulted in high correlation with the variables (Fig. 2, Fig. 3, Fig. 4, Fig. 5). This dominance of parameters related to sandfly population in the Zoonotic Visceral Leishmaniasis dynamics sensitivity suggests how dependent from the contact between sandfly and dog is. This fact may be very useful for planning activities regarding to sandfly population control.

We also observed that there are some parameters related to dog population in which the model is sensitive - , , , , and . Those, exception of , are all related to category. According to Laurenti et al. (2013), the latent dog has greater probability to infect sandflies than the clinically ill dogs. Thus, we are able to better understand why , , and are sensitive for our model, since those parameters model category dynamics. When any variation on those parameter values occurs, we are directly handling on the main category that composes the real source of infection. On the other hand, the parameter and are related to category, but the model is not so sensitive as that one related to category, since the contribution of as source of infection is lower than .

Numerical simulation

Finally, in order to analyze the dynamics of our model, we simulated our set of equations from (1) considering the parameters on Table 1. We focused on human reported rate, since we can use the real data as reference. Fig. 6, Fig. 7 illustrate the reported human cases rate and the sandflies per human ratio dynamics (yearly average).

Fig. 6

Dynamics of reported human cases rate. The available real data are from 1999 to 2015 (bars) and our model was fitted for the same period (line). Observe that the real data shows three peaks that decrease over time: 2002, 2007 and 2014. Source: CES-SP and BIGE.

Fig. 7

Dynamics of sandflies per human ratio. This curve was obtained from simulation of equation (5). Observe that there is a cycle and the peaks decrease over time. This curve becomes stable according to time progress.

First of all, we need to make some comments about Fig. 5, Fig. 6. First of all, although we based our modelling on the previous studies published by Burattini et al. (1998) and Ribas et al. (2013), our results are clearly different from them, in special, because we not only made changes on model, but also we used different parameter values. In addition, our approach is different, since we focus on incidence rate (human reported cases per day), instead of prevalence.

In order to compare our results, we were able to find some descriptions about some demographic characterization of those populations regarding to Zoonotic Visceral Leishmaniasis (see the endemic equilibrium states from Table 5). For instance, our sum of latent () and recovery () dog was around and the sum of latent () and clinically ill () dogs was around . On the other hand, Cabral et al. (1998) found that the sum of density of latent and recovered dogs was around 0.50 and Quinnell et al. (2001) found the density of 0.579 for the sum of latent and clinically ill dogs. Observe that our densities values are very different from those other authors. However, we remark that those works were not developed at Araçatuba/SP city. In addition, in Brazil the distribution of zoonotic visceral leishmaniasis is not homogeneus. Thus, we can find states with no human reported cases for 2012 year (as Rio Grande do Sul State) (Brasil, 2015). This non-homogeneus distribution is related to geographic featuring and climate changes, in which probably has influenced on the sandfly population dynamics. Therefore, this is one of the arguments to explain the discrepancy regarding to our result and the real data.

Discussion

The model presented in this article provided a new and interesting view about Zoonotic Visceral Leishmaniasis transmission dynamics among humans, dogs and sandflies population. Although the first work developed by our research group about this disease dynamics addresses from 1998 (Burattini et al., 1998), in this paper we included some new approaches. Herein, we not only updated most of parameters, but also analyzed the reported human cases rate dynamics and evaluated the model's sensitivity for parameters using LHS and PRCC as uncertainty and sensitivity analysis, respectively.

In our model we considered only dog population as source of infection, whereas Burattini et al. (1998) considered that the sandfly can acquire the protozoan from both human and dog populations. This explains why only dog population (in particular, latent and clinically ill dogs) composes the expressions (9) and (10). Once splitting this expression, we have the numerical result that latent dog contribution () is greater than clinically ill contribution (). We obtained this result due to the assumptions we adopted (the fraction of sandflies that become infected when bite a latent dog was and a clinically ill dog was ). This result addresses how important the latent dogs are for maintenance and introduction of Zoonotic Visceral Leishmaniasis into population. Also, if we consider that latent dogs are visually healthy, it is very difficult to detect them. Thus, latent dogs stay free to act as source of infection (Ministry of Health, Brazil, 2006).

In our model we could estimate the proportion of the individuals in each stage of disease dynamic. This evaluation is very important, since latent individuals (humans and dogs) are difficult to detect. In particular, the contribution of latent dog to disease maintenance is greater than the clinically ill one. Furthermore, here we estimated the reported human cases rate for our model in the equilibrium state /day. This value is at the same order of magnitude of the 2015 normalized yearly average reported rate (Table 2).

The real importance of reporting human cases of Zoonotic Visceral Leishmaniasis was well illustrated by Killick-Kendrick (2010). They highlighted that the key to reach the Visceral Leishmaniasis control is education. Although education was not considered in our model but we see that public health education and epidemiological surveillance system are very close and work together. As example, infected people that neglected the visceral Leishmaniasis provide unreported cases of this disease. As consequence, few cases are reported and the control programs are undervalued. Thus, if we know the proportion of unreported cases, we could evaluate the efficacy of the surveillance service. Finally, this can indirectly influence the control programs, since the strengthening of the surveillance system capacity is essential to avoid the underreporting of human cases and to follow-up the infection behavior in canine and human population. Strong surveillance will certainly contribute to improve data quality for decision-makers in this complex scenario (Romero and Boelaert, 2010, Maia-Elkhoury et al., 2007).

Classically, Zoonotic Visceral Leishmaniasis transmission is intensified when the prevalence on dog population is high and there is sandly population available. However, all mathematical analysis on our model provided us a better understanding about how dog and sandfly populations influence on disease dynamics. In particular, we observed that the model is highly sensitive to sandfly population parameters and this is related to our findings on stability analysis, in which the Disease Free Equilibrium state is broken when some infective sandfly is introduced. We can also observe the importance of sandfly population on this dynamics when we compare Fig. 6, Fig. 7, where dynamics over time clearly influence on human reported cases curve. At the same time, our calculation showed that it depends on dog (latent and clinically ill categories) and sandfly populations. The sensitivity analysis also indicated that parameters related to latent and clinically ill dog dynamics have some influence on disease dynamics. But, sandfly population is more important than dog population regarding to disease dissemination.

The recent work published by Zhao et al. (2016), although it showed some similar conclusions to that one provided by our model, it also presented some important differences. Firstly, there are some differences on adopted assumptions. For instance, Zhao et al. (2016) assumed that all recovered dogs are always under treatment. In our model, we considered that a dog can become naturally recovered (Table 1). In addition, Zhao et al. (2016) considered that there is a migration rate regarding to sandfly population. In our model we did not include this assumption (Fig. 1).

Although both models were elaborated from a classical compartmental model approach (SEIR model), we have some differences. Mathematically, Zhao et al. (2016) modeled the disease dynamics in a simplest way, since they considered fewer parameters and did not applied delay terms for sandfly population dynamics. This structural difference explains the differences regarding to results. As example, in our model it was not found the coexistence among DFE and EE state. Therefore, according to our model, the phenomenon of backward bifurcation (which Zhao et al. (2016) have proven in their work) was not presented. On the other hand, Zhao et al. (2016) have proven the existence of backward bifurcation in their model.

Finally, Zhao et al. (2016) demonstrated mathematically the optimal control based on their model. Although we did not develop a deep optimal control analysis, our sensitivity analysis also allowed us to conclude that the control strategy should focus on sandfly population, since our system was very sensitive to parameters related sandfly population dynamics (Table 5 and Fig. 2, Fig. 3, Fig. 4, Fig. 5). On the other hand, our model was able to represent the trend pattern of ZVL in a Brazilian city (Araçatuba, SP), since we obtained the real data and, therefore, it allowed us to fit our simulated results to real data (Fig. 6).

In this work, we presented a new model for Zoonotic Visceral Leishmaniasis, considering only dogs as source of infection, different probabilities of infecting sandflies for latent and clinically ill dogs and updated parameters. Since our analysis pointed that the introduction and maintenance of this disease is related to sandfly population and latent and clinically ill dogs, the preventive control activities should be focused on them. In special, considering the presently ethical issues regarding to elimination of positive dog in Brazil and the highly sensitivity of disease dynamics on sandfly population, we recommend that the sandfly population control should be prioritized. The evaluation of preventive activities on Zoonotic Visceral Leishmaniasis control is in our upcoming works.

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

HJS and EM were responsible for study design and planning. HJS and JW conducted the mathematical analysis. HJS conducted model simulations. HJS, JW and EM contributed to results interpretation and discussion. HJS contributed to writing the manuscript. All authors read and approved the final version of the manuscript.