# Mathematical Model Analysis and Simulation of Visceral Leishmaniasis, Kashgar, Xinjiang, 2004–2016

Authors: Yateng Song; Tailei Zhang; Hui Li; Kai Wang; Xiaobo Lu
Journal: Complexity
 (2020)

---

## Abstract

Visceral leishmaniasis (VL), known as kala-azar, is a serious parasitic disease. After malaria, VL is the second largest parasitic killer. This paper focuses on the VL transmission around sandflies, dogs, and people. Kashgar is located on the southwestern edge of Xinjiang, where kala-azar parasite infection occurs every year. According to the cases reported in the Kashgar Prefecture from 2004 to 2016, we proposed a dynamic model based on these three populations. The SEIR model was established for human population, the SI model was established for sandfly population, and the SI model was established for dog population. We fitted the model to cumulative cases from 2004 to 2016 for the epidemic in Kashgar and predicted that the cumulative incidence of kala-azar in Kashgar would continue to increase, but its growth rate would gradually slow down, which means that the number of cases would gradually decrease every year. We also estimated the basic reproduction numberR0 = 1.76 (95% CI: 1.49–1.93). The sensitivity analysis shows that the mutual infection between sandfly and dog contributes the most to the basic reproduction number, while the transmission proportion of sandfly to the susceptible person and the mutual infection between sandfly and dog contribute the most to the number of leishmaniasis human cases. Therefore, according to the sensitivity analysis results, reducing the contact between sandflies and dogs is an effective way to reduce kala-azar.


---

## Body

## 1. Introduction

Visceral leishmaniasis (VL) is also called kala-azar, which is a chronic infectious disease caused byLeishmania infantum. VL has ranked as the second largest parasitic killer after malaria and draws worldwide attention because of its severity; about 12 million people are affected by it around the world [1]. The pathogen causing visceral leishmaniasis are L. donovani, L. infantum, and L. chagasi. The incubation period for VL is generally 3−6 months, at least 10 days, while the longest is 9 years [1]. Humans and dogs are the main infection sources while sandflies are the main carriers.The World Health Organization (WHO) lists leishmaniasis as the most easily neglected disease. It is reported that approximately 30% of new clinical cases and 58,000 deaths of leishmaniasis occur worldwide each year [2]. According to the characteristics of the source of infection, VL is mainly divided into three types: human-borne, canine-borne, and wildlife-borne in China [3]. In the early days of the People’s Republic of China, human-borne kala-azar was the main disease in 16 provinces, such as Shandong and Henan. After prevention and treatment, it was basically eliminated. In recent years, VL has been prevalent in six provinces including Xinjiang, Inner Mongolia, Gansu, Sichuan, Shaanxi, and Shanxi [3]. Among six provinces, the typical VL in Xinjiang is wildlife-borne, and some areas have human-borne VL, while canine-borne is the main transmission way in Gansu.Xinjiang is the largest provincial administrative region in China and is also a high-risk area for VL, especially in the Kashgar Prefecture in southern Xinjiang. The Kashgar Prefecture is located in the southwestern margin of the Xinjiang Uygur Autonomous Region. Due to its unique geographical environment, there are two types of visceral leishmaniasis, human and desert [4]. Since the 1990s, the incidence of VL in Kashgar has been increasing year by year. From 2005 to 2015, Xinjiang’s VL has experienced significant fluctuation. In 2007–2012, the incidence of VL in Kashgar accounted for more than 90% of the total population [5].We utilized the data of human leishmaniasis cases from 2004 to 2016 reported by the Center for Disease Control and Prevention in Xinjiang, and then we plot bar diagrams about total, sex, and age which are presented in Figures1(a)–1(c). From Figure 1(a), we could find that the reported cases of VL in Kashgar Prefecture of Xinjiang fluctuate considerably, with outbreaks occurring in October and November of 2008, 2009, and 2015, respectively. From Figure 1(b), it shows that there are more reported cases of VL in males than in females. Figure 1(c) shows that there are more reported cases of VL in the younger age group.Figure 1
The data of human leishmaniasis cases from 2004 to 2016. (a) The reported data are described at different years. (b) The reported data are described at different genders and (c) different ages from 2004 to 2016.
(a)
(b)
(c)In recent years, a lot of mathematical models have been established around humans, dogs, and sandflies to understand the transmission dynamics of VL (see for instance [6–14]). For example, Muhammad and Ali [6] considered the uniform mixing of the population, established the SEIHR model for the population and the SEIR model for the dog, which used sandfly as the medium, modeled and backwarded the branching of zoonotic visceral leishmaniasis, and calculated the basic reproduction number R0 for optimal control. Elmojtaba et al. [7] developed a mathematical model to study the dynamics of visceral leishmaniasis in Sudan, considering three different groups, analyzing the balance point and its stability and providing a basis for controlling and eliminating diseases. Zamir et al. [8] established SIR mathematical models for humans, hosts, and media, using the Routh–Hurwitz standard and next-generation methods to obtain threshold conditions and give numerical simulation results. There are three main control strategies in the controls of VL model: parameter control strategy, optimal control strategy, and control strategy selection using simulation. However, the most generalized control strategy is the parameter control strategy. When the parameters are adjusted, people can apply it to a real-world control strategy, and there is also a lot of research on optimal control strategy and control strategy selection using simulation. Simulation comparison is the most common method of VL mathematic control modeling in simulation. It can compare the human infected population with control and without control. Meanwhile, it also proves the effectiveness of the combined control strategies [9]. With reference to these established mathematical models, we have established a mathematical model based on logistic process.The purpose of this paper is to know the transmission state about VL among humans, dogs, and sandflies in the Kashgar Prefecture of Xinjiang. Firstly, we proposed a model to simulate the cumulative data of the Kashgar Prefecture and estimate the basic reproduction number and the dynamic behavior of the model. Then, sensitivity analysis was performed on the number of leishmaniasis human/dog/sandfly cases and the basic reproduction numberR0 based on some key parameters. Finally, we explored some effective strategies for the prevention and control of VL in Kashgar.The article is organized as follows. In Section2, this paper introduces the model of VL and gives the expression of the basic reproduction number and the parameter value. The dynamical behaviors of the model are analyzed to better understand the transmission trends of the disease. Some mathematical analyses are given in Section 3. The numerical simulations, prediction of the epidemic trends for the next decades, estimation of the basic reproduction number, and sensitivity analysis of the basic reproduction number and the number of infected humans/sandflies/dogs are presented in Section 4. In Section 5, some brief summaries and discussion are given.


## 2. Mathematical Model

In order to study the transmission of VL in Kashgar, we developed a mathematical model based on humans, dogs, and sandflies, where the human population is divided into four groups: the susceptible, the exposed, the infected, and the recovered denoted bySh, Eh, Ih, and Rh, respectively. The dog population is divided into two groups: the susceptible and the infected, denoted by Sr and Ir, respectively. We divide the sandfly population into two subclasses: the susceptible and the infected, denoted by Sv and Iv, respectively. And the total population for humans, dogs, and sandflies is Nh=Sh+Eh+Ih+Rh; Nv=Sv+Iv; Nr=Sr+Ir. The flowchart of VL transmission is illustrated in Figure 2.Figure 2
Flow diagram of the VL transmission model.The transmission process of VL is described by the following eight differential equations:(1)dShdt=r1Sh1−Sh+Eh+Ih+RhK1−β1ShIv,dEhdt=β1ShIv−μhEh−ωEh,dIhdt=ωEh−μh+δ1Ih−γIh,dRhdt=γIh−μhRh,dSvdt=r2Sv1−Sv+IvK2−β2Ih+β3IrSv,dIvdt=β2Ih+β3IrSv−μvIv,dSrdt=r3Sr1−Sr+IrK3−β4IvSr,dIrdt=β4IvSr−μr+δ2Ir.We assume that the susceptible individualsSh and dogs Sr are bitten by infected sandflies with β1 and β4 infection rates, respectively. And the transmission probability from an infected dog or an infected individual to a susceptible sandfly is β3 and β2, respectively. 1/ω and 1/γ indicate the latency and recovery rate of VL, respectively. The natural mortality rates of humans, vectors, and reserves are μh, μv, and μr, respectively. The death rate of humans and reserves caused by VL are δ1 and δ2, respectively. The susceptible human population is governed by the logistic growth with carrying capacity K1 as well as intrinsic growth rate r1. In the absence of disease, the vector population density grows according to a logistic curve carrying capacity K2, with an intrinsic growth rate r2. In the absence of disease, the dog population density grows according to a logistic curve carrying capacity K3, with an intrinsic growth rate r3.


## 3. Mathematical Analysis

### 3.1. Basic Reproduction Number

The basic reproduction numberR0 refers to the number of people infected by a patient during the average period of illness when all people are susceptible. R0 was calculated using the next-generation matrix [15, 16]. The infected compartments are Eh, Ih, Iv, and Ir, giving m=4, so we have(2)ℱ=β1ShIv0β2Ih+β3IrSvβ4IvSr,ν=μh+ωEhμh+δ1+γIh−ωEhμrIrμr+δ2Ir.Then,(3)F=00β1K1000000β2K20β3K300β4K30,V=μh+ω000−ωμh+δ1+γ0000μv0000μr+δ2.Thus,(4)FV−1=00β1K1μv00000ωβ2K2μh+ωμh+δ1+γβ2K2μh+δ1+γ0β3K2μr+δ200β4K3μv0.The reproduction number is given byR0=ρFV−1, and(5)R0=K2K3β3β4μh2+δ1μh+γμh+δ1ω+γω+μhω+ωK1β1β2δ2+μrμvδ2+μrμh+ωδ1+γ+μh.R0 can be used as a basic indicator in the study of VL propagation dynamics model. Usually, when R0=1, it can be used as a threshold for the demise of the disease. When R0>1, the disease will not die and eventually turn into endemic disease; when R0<1, the disease will die out naturally.


### 3.2. Dynamic Behaviors for the Model

The dynamical behavior of model (1) is analyzed to better understand the transmission trends of the disease. We begin with elementary properties of solutions to the model (1). From biological considerations, we study model (1) in the set(6)Γ=Sh,Eh,Ih,Rh,Sv,Iv,Sr,Ir∈ℝ+8:0<Sh≤K1,0<Sv≤K2,0<Sr≤K3,Nh≤μh+r12K14μhr1,Nv≤μv+r22K24μvr2,Nr≤μr+r32K34μrr3.In fact, by using(7)dShdt≤r1Sh1−ShK1,dSvdt≤r2Sv1−SvK2,dSrdt≤r3Sr1−SrK3,we can deduce that Sht≤K1, Svt≤K2, and Srt≤K3 whenever Sh0≤K1, Sv0≤K2, and Sr0≤K3. Adding the first four equations of (1) yields(8)dNhdt=r1Sh1−NhK1−μhNh−δ1Ih+μhSh≤r1+μhSh−r1Sh2K1−μhNh≤r1+μh2K14r1−μhNh.In view of system (8), it is easy to check that Nh′t≤0 wherever Nh≥r1+μh2K1/4μhr1. The same conclusion can be drawn for Nv and Nr. According to the above analysis, Γ is a maximum positive invariant set of (1). There always exists a disease-free equilibrium P0 (K1, 0, 0, 0, K2, 0, K3, 0). In the following, we will show the stability for the disease-free equilibrium.Theorem 1.
The disease-free equilibriumP0 is locally asymptotically stable if R0<1 and unstable if R0>1.Proof.
The Jacobian matrix of model (1) at P0 is given by(9)JP0=−r1−r1−r1−r10−β1K1000−μh+ω000β1K1000ω−μh+δ1+γ0000000γ−μh000000−β2K20−r2−r20−β3K200β2K200−μv0β3K200000−β4K3−r3−r300000β4K30−μr+δ2.
It is easy to see that there are four negative eigenvalues ofJP0: −r1,−r2,−r3,−μh. The other eigenvalues are determined by the following fourth-order equation:(10)λ4+a1λ3+a2λ2+a3λ+a4=0,where(11)a1=m1+m2+n+μv,a2=m1m2+μvn−β3β4K2K3+m1+m2μv+n,a3=m1m2μv+n+m1+m2μvn−β3β4K2K3−ωβ1β2K1K2,a4=m1m2μvn−β3β4K2K3−nωβ1β2K1K2,(12)m1=μh+ω,m2=μh+δ1+γ,n=μr+δ2.
The inequalityR0<1 implies that(13)β3β4K2K3<μvn,ωβ1β2K1K2<m1m2μv,m1m2β3β4K2K3+nωβ1β2K1K2<m1m2μvn.
From (13), we have ai>0 for i=1,2,3,4. On account of the above inequalities, we have(14)a1a2−a3=m1+m2+n+μvm1m2+μvn−β3β4K2K3+m1+m2μv+n−m1m2μv+n−m1+m2μvn−β3β4K2K3+ωβ1β2K1K2=m1m2m1+m2+μv+nμvn−β3β4K2K3+m1+m2+n+μvm1+m2μv+n+ωβ1β2K1K2>0.
Furthermore,(15)a3a1a2−a3−a12a4=m1m2μv+n+m1+m2μvn−β3β4K2K3−ωβ1β2K1K2m1m2m1+m2+μv+nμvn−β3β4K2K3+m1+m2+n+μvm1+m2μv+n+ωβ1β2K1K2−m1m2μvn−β3β4K2K3m1+m2+n+μv2+nωβ1β2K1K2m1+m2+n+μv2>m1+m2μvn−β3β4K2K3m1m2m1+m2+m1+m2+n+μvm1+m2μv+n+ωβ1β2K1K2−m1m2μvn−β3β4K2K3m1+m2+n+μv2=m12+m22μvn−β3β4K2K3m1+m2μv+n+m12+m22+m1m2μvn−β3β4K2K3μv+n2>0.
Therefore, each eigenvalue of equation (10) admits negative real part. When R0<1, the disease-free equilibrium P0 is locally asymptotically stable. When R0>1, we conclude from a4>0 that equation (10) has at least one positive real root; hence, the disease-free equilibrium P0 is unstable.Theorem 2.
AsR0<1, the disease-free equilibrium is globally asymptotically stable in Γ.Proof.
Let us consider the following Lyapunov function:(16)Vt=ωnβ2K2Eht+m1nβ2K2Iht+m1m2nIvt+m1m2β3K2Irt,where(17)m1=μh+ω,m2=μh+δ1+γ,n=μr+δ2.
The derivative ofVt along solutions of (1) is(18)V˙t=ωnβ1β2K2ShIv−m1m2nβ2K2Ih+m1m2nSvβ2Ih+β3Ir−m1m2μvnIv+m1m2β3β4K2SrIv−m1m2nβ3K2Ir≤m1m2β3β4K2K3+nωβ1β2K1K2−m1m2μvnIv=m1m2μvnR02−1Iv≤0.
By LaSalle’s invariance principle, the omega limit set of each solution starting fromΓ lies in an invariant set contained in(19)Ω=Sh,Eh,Ih,Rh,Sv,Iv,Sr,Ir∈Γ:Iv=0.
It can be verified that the only invariant set contained inΩ is the singleton P0. Then, the disease-free equilibrium P0 is globally asymptotically stable. This completes the proof. In the following, we will consider the dynamical behavior of model (1) when R0>1.Theorem 3.
For any solution (Sht,Eht,Iht,Rht,Svt,Ivt,Srt,Irt) of (1) with initial values Sh0>0, Eh0≥0, Ih0≥0, Rh0≥0, Sv0>0, Iv0≥0, Sr0>0, Ir0≥0, and Ih0+Iv0+Ir0>0, there exists ν>0 such that(20)liminft⟶∞Iht>ν,liminft⟶∞Ivt>ν,liminft⟶∞Irt>ν,as R0>1.Proof.
Let(21)X=Sh,Eh,Ih,Rh,Sv,Iv,Sr,Ir∈ℝ+8:Sh>0,Sv>0,Sr>0,X0=Sh,Eh,Ih,Rh,Sv,Iv,Sr,Ir∈X:Ih+Iv+Ir>0,∂X0=X\X0.
Next, we will show that system (1) is uniformly persistent with respect to (X0,∂X0). Obviously, X is positively invariant with respect to (1). Set a=maxμh+δ1+γ,μv,μr+δ2. The inequality Iht+Ivt+Irt′≥−aIht+Ivt+Irt and Ih0+Iv0+Ir0>0 imply Iht+Ivt+Irt>0 for all t≥0 as Ih0+Iv0+Ir0>0. Therefore, X0 is also a positive invariant set for system (1). Furthermore, by Theorem 3, there exists a compact set c in which all solutions of (1) initiated in X will enter and remain forever after. The compactness condition (C4.2) in Thieme [17] is easily verified for this set c. Denote(22)M∂=Sh0,Eh0,Ih0,Rh0,Sv0,Iv0,Sr0,Ir0:Sht,Eht,Iht,Rht,Svt,Ivt,Srt,Irt∈∂X0,t≥0.
We first claim that(23)M∂=Sh,0,0,Rh,Sv,0,Sr,0:Sh,Rh,Sv,Sr≥0.
Suppose thatSh0,Eh0,Ih0,Rh0,Sv0,Iv0,Sr0,Ir0∈M∂. It suffices to show Eht=0 for all t≥0. If it is not true, then there exists t0>0 such that Eht0>0. Then, Eh′t≥μh+ωEhg implies Eht>0 for all t≥t0. From the third equation of (1), it may be concluded that(24)Iht=e−μh+δ1+γt−t0Iht0+∫t0te−μh+δ1+γt−θωEhθdθ≥∫t0te−μh+δ1+γt−θωEhθdθ>0,for all t>t0. This is a contradiction with(25)Sh0,Eh0,Ih0,Rh0,Sv0,Iv0,Sr0,Ir0∈M∂.
Thus, (23) is valid.
Denote(26)Ω=∪ωSh0,Eh0,Ih0,Rh0,Sv0,Iv0,Sr0,Ir0:Sh0,Eh0,Ih0,Rh0,Sv0,Iv0,Sr0,Ir0∈M∂,where ωSh0,Eh0,Ih0,Rh0,Sv0,Iv0,Sr0,Ir0 is the omega limit set of the solution to (1) through Sh0,Eh0,Ih0,Rh0,Sv0,Iv0,Sr0,Ir0. Restricting (1) on M∂ yields(27)dShdt=r1Sh1−Sh+RhK1,dRhdt=−μhRh,dSvdt=r2Sv1−SvK2,dSrdt=r3Sr1−SrK3.
A trivial verification shows that system (27) has a unique equilibrium K1,0,K2,K3. Thus, P0K1,0,0,0,K2,0,K3,0 is the unique equilibrium of (1) in M∂. It is easily seen that K1,0,K2,K3 is globally asymptotically stable. Therefore, we have Ω=P0. And P0 is a covering of Ω, which is isolated and is acyclic (since there exists no solution in M∂ which links P0 to itself). Finally, the proof will be done if P0 is a weak repeller for X0, i.e.,(28)limsupt⟶∞distΨt,P0>0,where Ψt=Sht,Eht,Iht,Rht,Svt,Ivt,Srt,Irt is an arbitrary solution with initial value in X0. By Leenheer and Smith (Proof of Lemma 3.5, [18]), we only need to prove that WsP0∩X0=ϕ where WsP0 is the stable manifold of E0. Suppose it is not true, then there exists a solution Sht,Eht,Iht,Rht,Svt,Ivt,Srt,Irt in X0, such that as t⟶∞,(29)Sht⟶K1,Eht⟶0,Iht⟶0,Rht⟶0,Svt⟶K2,Ivt⟶0,Srt⟶K3,Irt⟶0.
WhenR0>1,(30)ωβ2K2μr+δ2β1K1+β3K2μh+ωμh+δ1+γβ4K3>μh+ωμh+δ1+γμr+δ2μv.
On account of (30), we may choose ρ>0,σ>0,η>0, and ϵ>0 such that(31)0,<ωβ2K2−ϵμr+δ2−ρ≪1,0,<β3K2−ϵμr+ωμh+δ1+γ−σ≪1,0,<η−μh+ωμh+δ1+γμr+δ2≪1.(32)ρβ1K1−ϵ+σβ4K3−ϵ−ημv>0.
From (31), we see that μh+ω/ωρ<μh+ωμr+δ2β2K2−ϵ<β2K2−ϵ/μh+δ1+γη, ηβ3K2−ϵ>σμr+δ2. Therefore, there exists ξ such that(33)μh+ωωρ<ξ<β2K2−ϵμh+δ1+γη.
Forϵ>0, by (29), there exists T>0 such that(34)K1−ϵ<Sht<K1+ϵ,K2−ϵ<Svt<K2+ϵ,K3−ϵ<Srt<K3+ϵ,for all t≥T. Let(35)Lt=ρEht+ξIht+ηIvt+σIrt.
The derivative ofL along the solution Sht,Eht,Iht,Rht,Svt,Ivt,Srt,Irt is given by(36)L′t=ξω−ρμh+ωEh+ηβ2Sv−ξμh+δ1+γIh+ρβ1Sh+σβ4Sr−ημvIv+ηβ3Sv−σμr+δ2Ir≥ξω−ρμh+ωEh+ηβ2K2−ϵ−ξμh+δ1+γIh+ρβ1K1−ϵ+σβ4K3−ϵ−ημvIv+ηβ3K2−ϵ−σμr+δ2Ir≥ϱLt,for all t≥T, where(37)ϱ=minξω−ρμh+ωρ,ηβ2K2−ϵ−ξμh+δ1+γξ,ρβ1K1−ϵ+σβ4K3−ϵ−ημvη,ηβ3K2−ϵ−σμr+δ2σ>0.
Hence,Lt⟶∞ as t⟶∞, which contradicts to the boundedness of Lt. This completes the proof.Remark 1.
Theorems1∼3 state that the basic reproduction number R0 is a sharp threshold value for model (1). As R0<1, the disease-free equilibrium P0 is globally stable, i.e., the disease will go to extinction. As R0>1, the disease is uniformly persistent, i.e., the disease will become an endemic in the meaning of persistence.


## 3.1. Basic Reproduction Number

The basic reproduction numberR0 refers to the number of people infected by a patient during the average period of illness when all people are susceptible. R0 was calculated using the next-generation matrix [15, 16]. The infected compartments are Eh, Ih, Iv, and Ir, giving m=4, so we have(2)ℱ=β1ShIv0β2Ih+β3IrSvβ4IvSr,ν=μh+ωEhμh+δ1+γIh−ωEhμrIrμr+δ2Ir.Then,(3)F=00β1K1000000β2K20β3K300β4K30,V=μh+ω000−ωμh+δ1+γ0000μv0000μr+δ2.Thus,(4)FV−1=00β1K1μv00000ωβ2K2μh+ωμh+δ1+γβ2K2μh+δ1+γ0β3K2μr+δ200β4K3μv0.The reproduction number is given byR0=ρFV−1, and(5)R0=K2K3β3β4μh2+δ1μh+γμh+δ1ω+γω+μhω+ωK1β1β2δ2+μrμvδ2+μrμh+ωδ1+γ+μh.R0 can be used as a basic indicator in the study of VL propagation dynamics model. Usually, when R0=1, it can be used as a threshold for the demise of the disease. When R0>1, the disease will not die and eventually turn into endemic disease; when R0<1, the disease will die out naturally.


## 3.2. Dynamic Behaviors for the Model

The dynamical behavior of model (1) is analyzed to better understand the transmission trends of the disease. We begin with elementary properties of solutions to the model (1). From biological considerations, we study model (1) in the set(6)Γ=Sh,Eh,Ih,Rh,Sv,Iv,Sr,Ir∈ℝ+8:0<Sh≤K1,0<Sv≤K2,0<Sr≤K3,Nh≤μh+r12K14μhr1,Nv≤μv+r22K24μvr2,Nr≤μr+r32K34μrr3.In fact, by using(7)dShdt≤r1Sh1−ShK1,dSvdt≤r2Sv1−SvK2,dSrdt≤r3Sr1−SrK3,we can deduce that Sht≤K1, Svt≤K2, and Srt≤K3 whenever Sh0≤K1, Sv0≤K2, and Sr0≤K3. Adding the first four equations of (1) yields(8)dNhdt=r1Sh1−NhK1−μhNh−δ1Ih+μhSh≤r1+μhSh−r1Sh2K1−μhNh≤r1+μh2K14r1−μhNh.In view of system (8), it is easy to check that Nh′t≤0 wherever Nh≥r1+μh2K1/4μhr1. The same conclusion can be drawn for Nv and Nr. According to the above analysis, Γ is a maximum positive invariant set of (1). There always exists a disease-free equilibrium P0 (K1, 0, 0, 0, K2, 0, K3, 0). In the following, we will show the stability for the disease-free equilibrium.Theorem 1.
The disease-free equilibriumP0 is locally asymptotically stable if R0<1 and unstable if R0>1.Proof.
The Jacobian matrix of model (1) at P0 is given by(9)JP0=−r1−r1−r1−r10−β1K1000−μh+ω000β1K1000ω−μh+δ1+γ0000000γ−μh000000−β2K20−r2−r20−β3K200β2K200−μv0β3K200000−β4K3−r3−r300000β4K30−μr+δ2.
It is easy to see that there are four negative eigenvalues ofJP0: −r1,−r2,−r3,−μh. The other eigenvalues are determined by the following fourth-order equation:(10)λ4+a1λ3+a2λ2+a3λ+a4=0,where(11)a1=m1+m2+n+μv,a2=m1m2+μvn−β3β4K2K3+m1+m2μv+n,a3=m1m2μv+n+m1+m2μvn−β3β4K2K3−ωβ1β2K1K2,a4=m1m2μvn−β3β4K2K3−nωβ1β2K1K2,(12)m1=μh+ω,m2=μh+δ1+γ,n=μr+δ2.
The inequalityR0<1 implies that(13)β3β4K2K3<μvn,ωβ1β2K1K2<m1m2μv,m1m2β3β4K2K3+nωβ1β2K1K2<m1m2μvn.
From (13), we have ai>0 for i=1,2,3,4. On account of the above inequalities, we have(14)a1a2−a3=m1+m2+n+μvm1m2+μvn−β3β4K2K3+m1+m2μv+n−m1m2μv+n−m1+m2μvn−β3β4K2K3+ωβ1β2K1K2=m1m2m1+m2+μv+nμvn−β3β4K2K3+m1+m2+n+μvm1+m2μv+n+ωβ1β2K1K2>0.
Furthermore,(15)a3a1a2−a3−a12a4=m1m2μv+n+m1+m2μvn−β3β4K2K3−ωβ1β2K1K2m1m2m1+m2+μv+nμvn−β3β4K2K3+m1+m2+n+μvm1+m2μv+n+ωβ1β2K1K2−m1m2μvn−β3β4K2K3m1+m2+n+μv2+nωβ1β2K1K2m1+m2+n+μv2>m1+m2μvn−β3β4K2K3m1m2m1+m2+m1+m2+n+μvm1+m2μv+n+ωβ1β2K1K2−m1m2μvn−β3β4K2K3m1+m2+n+μv2=m12+m22μvn−β3β4K2K3m1+m2μv+n+m12+m22+m1m2μvn−β3β4K2K3μv+n2>0.
Therefore, each eigenvalue of equation (10) admits negative real part. When R0<1, the disease-free equilibrium P0 is locally asymptotically stable. When R0>1, we conclude from a4>0 that equation (10) has at least one positive real root; hence, the disease-free equilibrium P0 is unstable.Theorem 2.
AsR0<1, the disease-free equilibrium is globally asymptotically stable in Γ.Proof.
Let us consider the following Lyapunov function:(16)Vt=ωnβ2K2Eht+m1nβ2K2Iht+m1m2nIvt+m1m2β3K2Irt,where(17)m1=μh+ω,m2=μh+δ1+γ,n=μr+δ2.
The derivative ofVt along solutions of (1) is(18)V˙t=ωnβ1β2K2ShIv−m1m2nβ2K2Ih+m1m2nSvβ2Ih+β3Ir−m1m2μvnIv+m1m2β3β4K2SrIv−m1m2nβ3K2Ir≤m1m2β3β4K2K3+nωβ1β2K1K2−m1m2μvnIv=m1m2μvnR02−1Iv≤0.
By LaSalle’s invariance principle, the omega limit set of each solution starting fromΓ lies in an invariant set contained in(19)Ω=Sh,Eh,Ih,Rh,Sv,Iv,Sr,Ir∈Γ:Iv=0.
It can be verified that the only invariant set contained inΩ is the singleton P0. Then, the disease-free equilibrium P0 is globally asymptotically stable. This completes the proof. In the following, we will consider the dynamical behavior of model (1) when R0>1.Theorem 3.
For any solution (Sht,Eht,Iht,Rht,Svt,Ivt,Srt,Irt) of (1) with initial values Sh0>0, Eh0≥0, Ih0≥0, Rh0≥0, Sv0>0, Iv0≥0, Sr0>0, Ir0≥0, and Ih0+Iv0+Ir0>0, there exists ν>0 such that(20)liminft⟶∞Iht>ν,liminft⟶∞Ivt>ν,liminft⟶∞Irt>ν,as R0>1.Proof.
Let(21)X=Sh,Eh,Ih,Rh,Sv,Iv,Sr,Ir∈ℝ+8:Sh>0,Sv>0,Sr>0,X0=Sh,Eh,Ih,Rh,Sv,Iv,Sr,Ir∈X:Ih+Iv+Ir>0,∂X0=X\X0.
Next, we will show that system (1) is uniformly persistent with respect to (X0,∂X0). Obviously, X is positively invariant with respect to (1). Set a=maxμh+δ1+γ,μv,μr+δ2. The inequality Iht+Ivt+Irt′≥−aIht+Ivt+Irt and Ih0+Iv0+Ir0>0 imply Iht+Ivt+Irt>0 for all t≥0 as Ih0+Iv0+Ir0>0. Therefore, X0 is also a positive invariant set for system (1). Furthermore, by Theorem 3, there exists a compact set c in which all solutions of (1) initiated in X will enter and remain forever after. The compactness condition (C4.2) in Thieme [17] is easily verified for this set c. Denote(22)M∂=Sh0,Eh0,Ih0,Rh0,Sv0,Iv0,Sr0,Ir0:Sht,Eht,Iht,Rht,Svt,Ivt,Srt,Irt∈∂X0,t≥0.
We first claim that(23)M∂=Sh,0,0,Rh,Sv,0,Sr,0:Sh,Rh,Sv,Sr≥0.
Suppose thatSh0,Eh0,Ih0,Rh0,Sv0,Iv0,Sr0,Ir0∈M∂. It suffices to show Eht=0 for all t≥0. If it is not true, then there exists t0>0 such that Eht0>0. Then, Eh′t≥μh+ωEhg implies Eht>0 for all t≥t0. From the third equation of (1), it may be concluded that(24)Iht=e−μh+δ1+γt−t0Iht0+∫t0te−μh+δ1+γt−θωEhθdθ≥∫t0te−μh+δ1+γt−θωEhθdθ>0,for all t>t0. This is a contradiction with(25)Sh0,Eh0,Ih0,Rh0,Sv0,Iv0,Sr0,Ir0∈M∂.
Thus, (23) is valid.
Denote(26)Ω=∪ωSh0,Eh0,Ih0,Rh0,Sv0,Iv0,Sr0,Ir0:Sh0,Eh0,Ih0,Rh0,Sv0,Iv0,Sr0,Ir0∈M∂,where ωSh0,Eh0,Ih0,Rh0,Sv0,Iv0,Sr0,Ir0 is the omega limit set of the solution to (1) through Sh0,Eh0,Ih0,Rh0,Sv0,Iv0,Sr0,Ir0. Restricting (1) on M∂ yields(27)dShdt=r1Sh1−Sh+RhK1,dRhdt=−μhRh,dSvdt=r2Sv1−SvK2,dSrdt=r3Sr1−SrK3.
A trivial verification shows that system (27) has a unique equilibrium K1,0,K2,K3. Thus, P0K1,0,0,0,K2,0,K3,0 is the unique equilibrium of (1) in M∂. It is easily seen that K1,0,K2,K3 is globally asymptotically stable. Therefore, we have Ω=P0. And P0 is a covering of Ω, which is isolated and is acyclic (since there exists no solution in M∂ which links P0 to itself). Finally, the proof will be done if P0 is a weak repeller for X0, i.e.,(28)limsupt⟶∞distΨt,P0>0,where Ψt=Sht,Eht,Iht,Rht,Svt,Ivt,Srt,Irt is an arbitrary solution with initial value in X0. By Leenheer and Smith (Proof of Lemma 3.5, [18]), we only need to prove that WsP0∩X0=ϕ where WsP0 is the stable manifold of E0. Suppose it is not true, then there exists a solution Sht,Eht,Iht,Rht,Svt,Ivt,Srt,Irt in X0, such that as t⟶∞,(29)Sht⟶K1,Eht⟶0,Iht⟶0,Rht⟶0,Svt⟶K2,Ivt⟶0,Srt⟶K3,Irt⟶0.
WhenR0>1,(30)ωβ2K2μr+δ2β1K1+β3K2μh+ωμh+δ1+γβ4K3>μh+ωμh+δ1+γμr+δ2μv.
On account of (30), we may choose ρ>0,σ>0,η>0, and ϵ>0 such that(31)0,<ωβ2K2−ϵμr+δ2−ρ≪1,0,<β3K2−ϵμr+ωμh+δ1+γ−σ≪1,0,<η−μh+ωμh+δ1+γμr+δ2≪1.(32)ρβ1K1−ϵ+σβ4K3−ϵ−ημv>0.
From (31), we see that μh+ω/ωρ<μh+ωμr+δ2β2K2−ϵ<β2K2−ϵ/μh+δ1+γη, ηβ3K2−ϵ>σμr+δ2. Therefore, there exists ξ such that(33)μh+ωωρ<ξ<β2K2−ϵμh+δ1+γη.
Forϵ>0, by (29), there exists T>0 such that(34)K1−ϵ<Sht<K1+ϵ,K2−ϵ<Svt<K2+ϵ,K3−ϵ<Srt<K3+ϵ,for all t≥T. Let(35)Lt=ρEht+ξIht+ηIvt+σIrt.
The derivative ofL along the solution Sht,Eht,Iht,Rht,Svt,Ivt,Srt,Irt is given by(36)L′t=ξω−ρμh+ωEh+ηβ2Sv−ξμh+δ1+γIh+ρβ1Sh+σβ4Sr−ημvIv+ηβ3Sv−σμr+δ2Ir≥ξω−ρμh+ωEh+ηβ2K2−ϵ−ξμh+δ1+γIh+ρβ1K1−ϵ+σβ4K3−ϵ−ημvIv+ηβ3K2−ϵ−σμr+δ2Ir≥ϱLt,for all t≥T, where(37)ϱ=minξω−ρμh+ωρ,ηβ2K2−ϵ−ξμh+δ1+γξ,ρβ1K1−ϵ+σβ4K3−ϵ−ημvη,ηβ3K2−ϵ−σμr+δ2σ>0.
Hence,Lt⟶∞ as t⟶∞, which contradicts to the boundedness of Lt. This completes the proof.Remark 1.
Theorems1∼3 state that the basic reproduction number R0 is a sharp threshold value for model (1). As R0<1, the disease-free equilibrium P0 is globally stable, i.e., the disease will go to extinction. As R0>1, the disease is uniformly persistent, i.e., the disease will become an endemic in the meaning of persistence.


## 4. Results

We applied model (1) to study the infection status of VL in Kashgar, Xinjiang. Most of the parameters were obtained from the literature, and some of them were assumed or simulated to have more realistic results. These parameter values are listed in Table 1. The interpretation of the parameters values are as follows:(1)
According to the population statistics in Kashgar Prefecture in Xinjiang Statistical Yearbook 2004–2016 [20], we have estimated the carrying capacity of population K1 = 5000000.(2)
According to the terrain and climate factors of Kashgar and the survival habits of the sandflies [21–24], the carrying capacity of sandflies K2 = 7000000 has been estimated in Kashgar. In 2015, the number of dogs in Kashgar Prefecture of Xinjiang was close to 80,000 [25], so we have assumed the carrying capacity of dogs in Xinjiang K3 = 100000.(3)
The values of the parameterr1 are obtained by numerical differentiation according to the population data in Kashgar, and the value of r2 is calculated according to the life cycle of the sandflies (see [19]); the value of r3 is assumed because there is no the specific data about dog.Table 1
Parameter values (and their sources).
Parameters
Value
Interpretation∗
Source
r
1
0.019
Intrinsic growth rate of human
Calculated
r
2
0.026
Intrinsic growth rate of vector
[19]
r
3
0.296
Intrinsic growth rate of reserve
Assumption
K
1
5000000
Human carrying capacity
Assumption
K
2
7000000
Vector carrying capacity
Assumption
K
3
100000
Reserve carrying capacity
Assumption
β
1
1.35e − 08
Transmission probability from an infected sandfly to a susceptible human
Fitting
β
2
3.565e − 17
Transmission probability from an infected human to a susceptible sand fly
Fitting
β
3
3.515e − 05
Transmission probability from an infected reserve to a susceptible sandfly
Fitting
β
4
1.23e − 05
Transmission probability from an infected sandfly to a susceptible reserve
Fitting
μ
h
0.0125
Natural mortality rate of human
[10]
μ
v
25.55
Natural mortality rate of vector
[11]
μ
r
0.1
Natural mortality rate of reserve
[12]
1
/
ω
0.33
Human incubation period
[13]
δ
1
4.015
VL-induced death rate of humans
[8]
δ
2
3.65
VL-induced death rate of reserve
[13]
1
/
γ
2.5
Recovery period of VL
[14]We model cumulative cases as a Poisson-distributed random variable because the Poisson distribution describes the number of observed events in an interval of time. We calibrate the model by sampling from the posterior distribution of parameter vectorθy=β1,β2,β3,β4y, where vector y is derived from d/dtYt=ωEh and Yt denotes the reported cumulative cases. We conduct sampling via MCMC (Markov Chain Monte Carlo) using the Metropolis–Hastings acceptance rule. The posterior density is(38)fΘyθy=∏TℒYt|θfΘθ.The prior densityfΘθ is the joint probability of four univariate priors. We consider that β1, β2, β3, and β4 are distributed according to u0,1. The program was implemented in R version 3.6.0. We sampled from 30,000 MCMC iterations and discarded the first 10,000 samples as a burn-in period. On the basis of these 20,000 samples, the point estimates and 95% confidence intervals for the transmission coefficients were calculated. The results are shown in Table 2.Table 2
Parameter values for point estimation and 95% interval estimation.
Parameters
Point estimation
95% confidence interval
β
1
1.35e − 08
[1.06e − 08, 1.94e − 08]
β
2
3.565e − 17
[2.87e − 19, 2.44e − 15]
β
3
3.515e − 05
[2.067e − 05, 1.002e − 04]
β
4
1.23e − 05
[3.18e − 06, 2.36e − 05]The MCMC method was used to fit the cumulative incidence data of VL in Kashgar from 2004 to 2016, and the 95% confidence interval of the fitted curve was obtained (see Figure3). It can be seen from the figure that the fitted values of the model matches are in accordance with the accumulated data values in the Kashgar Prefecture. Only some data fluctuate, but in any case, the cumulative data of VL in Kashgar are increasing year by year. This model is used to predict the prevalence of VL in the Kashgar Prefecture over the next decade (see Figure 4).Figure 3
Cumulative incidence simulation of VL in Kashgar from 2004 to 2016.Figure 4
Prediction of cumulative incidence in 10 years from 2016 to 2026.The basic reproduction numberR0 = 1.76(95% CI: 1.49–1.93). The result shows that R0>1; according to the threshold theory, the disease will not disappear in the Kashgar Prefecture in a short period of time, and an endemic disease will be formed.It is well known that the basic reproduction numberR0 is a very important parameter in the infectious disease model. In our model, R0 is determined by the parameters of ω, δ1, δ2, γ, μh, μr, μv, β1, β2, β3, and β4. We use the Latin hypercube sampling (LHS) and partial rank correlation coefficients (PRCCs) to examine parameters which have a significant influence on the transmission of VL [26]. Using model (1), 2000 and 3000 samples are randomly generated by assuming a uniform distribution for each parameter based on values from Table 1.We select eleven parameters as input variables and calculate the corresponding PRCC andp values for the eleven parameters. Table 3 and Figure 5 show the exact PRCC and p values of each input parameter and the effect on the basic reproduction number R0, respectively. We assume the significance level α=0.05. The larger the absolute value of PRCC, the stronger the correlation between the input parameters and R0. It can be seen from Table 3 that only δ2, μv, β3, and β4 have significant impact on R0. More concretely, parameters δ2 and μv have negative impact on R0 and parameters β3 and β4 have positive impact.Table 3
Partial rank correlation coefficients (PRCCs) forR0 on each input parameter variable.
Parameters
Distribution
2000 samples
3000 samples
PRCC
p value
PRCC
p value
ω
N
3.03
,
0.15
2
9.295e − 04
0.504
4.379e − 04
0.499
δ
1
N
4.015
,
0.2
2
−1.995e − 04
0.493
−1.619e − 04
0.496
δ
2
N
3.65
,
0.2
2
−0.247
<
0.001
−0.2474
<
0.001
γ
N
0.4
,
0.02
2
−1.281e − 03
0.493
3.011e − 04
0.502
μ
h
N
0.0125
,
0.0006
2
8.376e − 04
0.511
−1.041e − 03
0.505
μ
r
N
0.1
,
0.005
2
−6.862e − 03
0.481
−6.146e − 03
0.495
μ
v
N
25.55
,
3.28
2
−0.2327
<
0.001
−0.2334
<
0.001
β
1
N (1.35e − 08, 1e − 08)
3.792e − 04
0.506
−5.341e − 05
0.502
β
2
N (3.565e − 08, 2e − 172)
−3.859e − 04
0.493
−3.456e − 04
0.509
β
3
N (3.515e − 08, 1.5e − 052)
0.914
<
0.001
0.914
<
0.001
β
4
N (1.23e − 08, 5e − 062)
0.906
<
0.001
0.906
<
0.001Figure 5
Partial rank correlation coefficient (PRCC) results for the dependence ofR0 on each parameter.We perform sensitivity analysis of infected humans, infected sandflies, and infected dogs through evaluating the PRCCs with the parameters of interest of model (1) over time by choosing a normal distribution with mean value and standard deviation shown in Figure 6. In Figures 6(a)–6(c), we plot the PRCCs over time with respect to the infected humans, sandflies, and infected dogs, respectively. Figure 6(a) indicates that there are three PRCC values that are significantly different from zero. The first three parameters with most impact on the outcome (the number of infected humans) are the sandfly-to-human transmission (β1), transmission from an infected dog to a susceptible sandfly (β3), and transmission from an infected sandfly to a susceptible dog (β4). Figure 6(c) indicates that there are two PRCC values that are significantly different from zero. The first two parameters with most impact on the outcome (the number of infected humans) are the transmission from an infected dog to a susceptible sandfly (β3) and transmission from an infected sandfly to a susceptible dog (β4). Figure 6(c) indicates that there are two PRCC values that are significantly different from zero. The first two parameters with most impact on the outcome (the number of infected humans) are the transmission from an infected sandfly to a susceptible dog (β4) and transmission from an infected dog to a susceptible sandfly (β3).Figure 6
Sensitivity analysis: PRCC of the eleven parameters for (a) the number of infected humans, (b) the number of infected sandflies, and (c) the number of infected dogs.
(a)
(b)
(c)Through sensitivity analysis, we demonstrate that the infection rate from infected dogs to susceptible sandfliesβ3 and infected sandflies to susceptible dogs β4 is the most sensitive parameter of R0. And β1, β3, and β4 are the most sensitive parameters for the number of leishmaniasis human cases. Therefore, it is necessary to change the values of the parameters to observe their effects on the number of leishmaniasis human cases and R0 (see Figure 7). From Figure 7(a), we observe that the number of leishmaniasis human cases decrease with the decrease of β1. We change the values of β3 and β4 to 4, 2, 1/2, and 1/4 times the original values. We observe the effects of different values on the number of leishmaniasis human cases, and we find that the change of parameter β3, β4 can influence not only the number of leishmaniasis human cases, but also the peak time. As Figures 7(b) and 7(c) illustrate, when fixing other parameters at constant, the number of leishmaniasis human cases fall with the decrease of β3 and β4, respectively. And the peak of outbreak will be postponed.Figure 7
The influence of different values of parameters on the number of diseases of VL. (a) Versusβ1. (b) Versus β3. (c) Versus β4.
(a)
(b)
(c)Finally, in order to find better control strategies for VL transmission, we focus on changing the values of the parametersβ3 and β4 to confirm the influence on the basic reproduction number R0 (see Figure 8).Figure 8
The influence of parameters onR0 (red horizontal line is contour of R0=1). (a) Versus β3. (b) Versus β4.
(a)
(b)It can be seen from the figure that the parametersβ3 and β4 have a strong influence on R0, and the value of the basic reproduction number R0 increased with the increases of β3 and β4; when β3 is less than 6.854878e − 06 or β4 is less than 2.39872e − 06, R0<1, and the disease can be eliminated. This suggests that reducing the proportion rate between sandflies and dogs can effectively control the prevalence of VL.Through the above sensitivity analysis results, we find that the main cause of the outbreak of VL is the mutual contact infection between dogs and sandflies. Therefore, we can get some effective strategies to reduce the prevalence of VL: the effective way is to reduce the contact between dogs and sandflies, and thus we can give the dog a collar with impregnated insecticide, increase the control of dog and reduce its active area, spray insecticides vigorously, etc.


## 5. Conclusion and Discussion

VL is a serious parasitic disease. It has been endangered for several decades in Xinjiang and has become a major parasitic disease affecting the local social and economic development. Kashgar is a high-risk area for the occurrence of VL. In order to reveal the spread of VL in Kashgar and predict the prevalence of VL, this paper proposes a dynamic model of VL propagation with logistic growth. The model describes the transmission of VL among humans, dogs, and sandflies.We use model (1) to fit the cumulative data of VL in Kashgar. As it can be seen from the simulation results in Figure 2, our model is consistent with the actual data of the cumulative cases of VL in Kashgar. The results show that there are certain reliability and rationality to study the prevalence of VL in Kashgar using logistic growth of the VL model. Using model (1), the basic reproduction number is estimated to be 1.76 (95% CI: 1.49–1.93) in the Kashgar Prefecture of Xinjiang. According to the threshold theory, it shows that VL will not disappear in the Kashgar area in a short time, and it may gradually become an endemic disease. According to the predictions for the next decade (see Figure 4), the cumulative incidence of VL in Kashgar is growing slowly, which means that the number of cases will gradually decrease over the next decade.By selecting the sensitivity analysis of the parameters of interest for the basic reproduction numberR0 and the number of leishmaniasis human cases (see Figures 5 and 6), we can find that the most important factors affecting the basic reproduction number R0 are β3 and β4, indicating that the infection between sandfly and dog has the greatest impact on R0, and the most sensitive parameters affecting the number of leishmaniasis human cases are β1, β3, and β4, which explains that the infection among the infected person, the dog, and the sandfly has the strongest impact on the number of diseases of VL. From Figure 7, we find that reducing the values of β1, β3, and β4 can effectively reduce the value of the basic reproduction number R0 and the number of leishmaniasis human cases. When the values of β3 and β4 drop to 6.854878e − 06 and 2.39872e − 06, respectively, the value of the basic reproduction number R0 will decrease to 1, and VL will not be epidemic but gradually disappear in Kashgar. Therefore, in order to control the spread of the disease, effective strategies should be taken to prevent and control leishmaniasis in Kashgar; we can decrease the incidence of leishmaniasis in humans by reducing the contact between sandflies and dogs. Without considering costs, we can vector controls (e.g., environmental clean-up and insecticide sprayed around buildings) and dog controls (e.g., insecticide releasing dog collars and dog vaccinations). In order to reduce the bite of the sandflies, we can use insecticidal bed nets. Meanwhile, we should conduct public education on the dangers and prevention of leishmaniasis for people. In short, the model we have established now can reflect the dynamics of VL in Kashgar Prefecture.



---

*Source: 5049825-2020-04-24.xml*