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

College of Public Health, Xinjiang Medical University, Urumqi 830011, China School of Science, Chang’an University, Xi’an 710064, China Central Laboratory, Xinjiang Medical University, Urumqi 830011, China College of Medical Engineering and Technology, Xinjiang Medical University, Urumqi 830011, China Department of Infectious Diseases, 0e First Affiliated Hospital of Xinjiang Medical University, Urumqi 830054, China


Introduction
Visceral leishmaniasis (VL) is also called kala-azar, which is a chronic infectious disease caused by Leishmania 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]. e pathogen causing visceral leishmaniasis are L. donovani, L. infantum, and L. chagasi. e 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.
e 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 wildlifeborne 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. e 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 Figures 1(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.
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][7][8][9][10][11][12][13][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 R 0 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. ere 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. e 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. en, sensitivity analysis was performed on the number of leishmaniasis human/dog/ sandfly cases and the basic reproduction number R 0 based on some key parameters. Finally, we explored some effective strategies for the prevention and control of VL in Kashgar. e article is organized as follows. In Section 2, this paper introduces the model of VL and gives the expression of the basic reproduction number and the parameter value. e dynamical behaviors of the model are analyzed to better understand the transmission trends of the disease. Some mathematical analyses are given in Section 3. e 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.

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 by S h , E h , I h , and R h , respectively. e dog population is divided into two groups: the susceptible and the infected, denoted by S r and I r , respectively. We divide the sandfly population into two subclasses: the susceptible and the infected, denoted by S v and I v , respectively. And the total population for humans, dogs, and sandflies is e flowchart of VL transmission is illustrated in Figure 2.
e transmission process of VL is described by the following eight differential equations: We assume that the susceptible individuals S h and dogs S r are bitten by infected sandflies with β 1 and β 4 infection 2 Complexity 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/c indicate the latency and recovery rate of VL, respectively. e natural mortality rates of humans, vectors, and reserves are μ h , μ v , and μ r , respectively. e death rate of humans and reserves caused by VL are δ 1 and δ 2 , respectively. e susceptible human population is governed by the logistic growth with carrying capacity K 1 as well as intrinsic growth rate r 1 . In the absence of disease, the vector population density grows according to a logistic curve carrying capacity K 2 , with an intrinsic growth rate r 2 . In the absence of disease, the dog population density grows according to a logistic curve carrying capacity K 3 , with an intrinsic growth rate r 3 .

Basic Reproduction Number.
e basic reproduction number R 0 refers to the number of people infected by a patient during the average period of illness when all people are susceptible. R 0 was calculated using the next-generation en, us, e reproduction number is given by R 0 � ρ(FV − 1 ), and

Complexity
R 0 can be used as a basic indicator in the study of VL propagation dynamics model. Usually, when R 0 � 1, it can be used as a threshold for the demise of the disease. When R 0 > 1, the disease will not die and eventually turn into endemic disease; when R 0 < 1, the disease will die out naturally.

Dynamic Behaviors for the Model.
e 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 In fact, by using we can deduce that Adding the first four equations of (1) yields . e same conclusion can be drawn for N v and N r . According to the above analysis, Γ is a maximum positive invariant set of (1). ere always exists a disease-free equilibrium P 0 (K 1 , 0, 0, 0, K 2 , 0, K 3 , 0). In the following, we will show the stability for the disease-free equilibrium.

Proof.
e Jacobian matrix of model (1) at P 0 is given by It is easy to see that there are four negative eigenvalues of J(P 0 ): −r 1 , −r 2 , −r 3 , −μ h . e other eigenvalues are determined by the following fourth-order equation: where Complexity e inequality R 0 < 1 implies that From (13), we have a i > 0 for i � 1, 2, 3, 4. On account of the above inequalities, we have Furthermore, erefore, each eigenvalue of equation (10) admits negative real part. When R 0 < 1, the disease-free equilibrium P 0 is locally asymptotically stable. When R 0 > 1, we conclude from a 4 > 0 that equation (10) has at least one positive real root; hence, the disease-free equilibrium P 0 is unstable.
Proof. Let us consider the following Lyapunov function: where By LaSalle's invariance principle, the omega limit set of each solution starting from Γ lies in an invariant set contained in 6 Complexity It can be verified that the only invariant set contained in Ω is the singleton P 0 . en, the disease-free equilibrium P 0 is globally asymptotically stable. is completes the proof. In the following, we will consider the dynamical behavior of model (1) when R 0 > 1.

Theorem 3. For any solution (S
Next, we will show that system (1) is uniformly persistent with respect to (X 0 , zX 0 ). Obviously, X is positively invariant with respect to (1) erefore, X 0 is also a positive invariant set for system (1). Furthermore, by eorem 3, there exists a compact set c in which all solutions of (1) initiated in X will enter and remain forever after. e compactness condition (C 4.2 ) in ieme [17] is easily verified for this set c. Denote We first claim that From the third equation of (1), it may be concluded that for all t > t 0 . is is a contradiction with us, (23) is valid. Denote where ω(S h (0), E h (0), I h (0), R h (0), S v (0), I v (0), S r (0), I r (0)) is the omega limit set of the solution to (1) through (S h (0), E h (0), I h (0), R h (0), S v (0), I v (0), S r (0), I r (0)).
A trivial verification shows that system (27) has a unique equilibrium (K 1 , 0, K 2 , K 3 ). us, P 0 (K 1 , 0, 0, 0, K 2 , 0, K 3 , 0) is the unique equilibrium of (1) in M z . It is easily seen that (K 1 , 0, K 2 , K 3 ) is globally asymptotically stable. erefore, we have Ω � P 0 . And P 0 is a covering of Ω, which is isolated and is acyclic (since there exists no solution in M z which links P 0 to itself ). Finally, the proof will be done if P 0 is a weak repeller for X 0 , i.e., lim sup t⟶∞ dist Ψ(t), P 0 > 0, where is an arbitrary solution with initial value in X 0 . By Leenheer and Smith (Proof of Lemma 3.5, [18]), we only need to prove that W s (P 0 ) ∩ X 0 � ϕ where W s (P 0 ) is the stable manifold of E 0 . Suppose it is not true, then there exists a solution (S h (t), E h (t), I h (t), R h (t), S v (t), I v (t), S r (t), I r (t)) in X 0 , such that as t ⟶ ∞, When R 0 > 1, On account of (30), we may choose ρ > 0, σ > 0, η > 0, and ϵ > 0 such that From (31), we see that For ϵ > 0, by (29), there exists T > 0 such that for all t ≥ T. Let Hence, L(t) ⟶ ∞ as t ⟶ ∞, which contradicts to the boundedness of L(t). is completes the proof.

Remark 1.
eorems 1∼3 state that the basic reproduction number R 0 is a sharp threshold value for model (1). As R 0 < 1, the disease-free equilibrium P 0 is globally stable, i.e., the disease will go to extinction. As R 0 > 1, the disease is uniformly persistent, i.e., the disease will become an endemic in the meaning of persistence.

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. ese parameter values are listed in Table 1 (3) e values of the parameter r 1 are obtained by numerical differentiation according to the population data in Kashgar, and the value of r 2 is calculated according to the life cycle of the sandflies (see [19]); the value of r 3 is assumed because there is no the specific data about dog.

Complexity
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 , β 4 | y, where vector y is derived from (d/dt)Y(t) � ωE h and Y(t) denotes the reported cumulative cases. We conduct sampling via MCMC (Markov Chain Monte Carlo) using the Metropolis-Hastings acceptance rule. e posterior density is e prior density f Θ (θ) is the joint probability of four univariate priors. We consider that β 1 , β 2 , β 3 , and β 4 are distributed according to u(0, 1). e 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. e results are shown in Table 2.
e 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 Figure 3). 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. is model is used to predict the prevalence of VL in the Kashgar Prefecture over the next decade (see Figure 4). e basic reproduction number R 0 � 1.76(95% CI: 1.49-1.93). e result shows that R 0 > 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 number R 0 is a very important parameter in the infectious disease model. In our model, R 0 is determined by the parameters of ω, δ 1 , δ 2 , c, μ 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 and p 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 R 0 , respectively. We assume the significance level α � 0.05. e larger the absolute value of PRCC, the stronger the correlation between the input parameters and R 0 . It can be seen from Table 3 that only δ 2 , μ v , β 3 , and β 4 have significant impact on R 0 . More concretely, parameters δ 2 and μ v have negative impact on R 0 and parameters β 3 and β 4 have positive impact.
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. e 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. e 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. e 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 ). rough 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 R 0 . And β 1 , β 3 , and β 4 are the most sensitive parameters for the number of leishmaniasis human cases. erefore, it is necessary to change the values of the parameters to observe their effects on the number of leishmaniasis human cases and R 0 (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.  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 R 0 (see Figure 8).
It can be seen from the figure that the parameters β 3 and β 4 have a strong influence on R 0 , and the value of the basic reproduction number R 0 increased with the increases of β 3 and β 4 ; when β 3 is less than 6.854878e − 06 or β 4 is less than 2.39872e − 06, R 0 < 1, and the disease can be eliminated.
is suggests that reducing the proportion rate between sandflies and dogs can effectively control the prevalence of VL.    rough 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. erefore, 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.

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. e 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. e 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 number R 0 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 R 0 are β 3 and β 4 , indicating that the infection between sandfly and dog has the greatest impact on R 0 , 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 R 0 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 R 0 will decrease to 1, and VL will not be epidemic but gradually disappear in Kashgar. erefore, 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.
Data Availability e data will be available upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.

Authors' Contributions
Yateng Song and Tailei Zhang contributed equally.