Mathematical Modeling of Echinococcosis in Humans, Dogs, and Sheep

In this paper, a model for the transmission dynamics of cystic echinococcosis in the dog, sheep, and human populations is developed and analyzed. We first model and analyze the predator-prey interaction model in these populations; then, we propose a mathematical model of the transmission dynamics of cystic echinococcosis. We calculate the basic reproduction number R0 and prove that the disease-free equilibrium is globally asymptotically stable, and hence, the disease dies out if R0 < 1. We further show that the endemic equilibrium is globally asymptotically stable, and hence, the disease persists if R0 > 1. Numerical simulations are performed to illustrate our analytic results. We give sensitivity analysis of the key parameters and give strategies that are helpful to control the transmission of cystic echinococcosis, from which the most sensitive parameter is the transmission rate of Echinococcus’ eggs from the environment to sheep (βes). Thus, the effective controlling strategies are associated with this parameter.


Introduction
Echinococcosis is a parasitic disease caused by ingesting the eggs of the tapeworm genus Echinococcus through contaminated food and water in the environment. The parasite's life cycle is maintained by carnivores, which act as definitive hosts (a dog, fox, canine, felid, or hyenid), and by intermediate hosts, which are usually herbivores (e.g., sheep, goats, cattle, camels, and cervids) [1]. Different species of Echinococcus cause different diseases in humans. Of the different forms of echinococcosis occurring in humans, cystic echinococcosis (CE) and alveolar echinococcosis (AE) are of special importance due to their wide geographic distribution and their medical and economic impact [2,3]. Echinococcus granulosus (E. granulosus) causes cystic echinococcosis while E. multilocularis causes a type of echinococcosis known as alveolar echinococcosis [2][3][4].
Cystic echinococcosis (CE) is a parasitic disease, also called "cystic hydatid disease" caused by the larval stage of small tapeworms (dog tapeworms) known as Echinococcus granulosus [2]. In its transmission dynamics, the domestic dog is the principal definitive host. The parasite is transmit-ted to dogs when they ingest the organs of infected animals (such as sheep) that contain hydatid cysts. The cysts develop into adult tapeworms in the dog. Infected dogs shed tapeworm eggs in their feces which contaminate the environment. Sheep, cattle, goats, and pigs ingest tapeworm eggs in the contaminated ground; once ingested, the eggs hatch and develop into cysts in the internal organs. The most common mode of transmission to humans is by the accidental consumption of water or food that has been contaminated by the fecal matter of an infected dog. Echinococcus eggs that have been deposited in the soil stay viable for up to a year [5][6][7].
Geographically, the greatest prevalence of cystic echinococcosis in human and animal hosts is found in Australia, some parts of America (especially South America), central Asia, northern and eastern Africa, and the Mediterranean Basin [4,5,7]. Cystic echinococcosis typically occurs in poor pastoral regions in which sheep or other livestock are raised and in which dogs are kept, for herding or property guarding, in close proximity to households. Dogs in such regions are frequently fed offal, and, for religious and other reasons, their populations might not be curtailed [8]. The prevalence of the parasite within the offal and the contact rate between sheep and dogs will play an important role in the rate of infection within these populations [9][10][11].
In an ecosystem, predator-prey interaction between species has been understood as a determining factor in the population dynamics. Predation is a biological interaction where a predator feeds on its prey. Predators reduce prey densities while consuming individual prey. In the absence of predators, the population of prey can grow until it approaches the carrying capacity of the environment. In turn, prey densities affect the number of predators because, when prey becomes scarce, predators die of starvation or fail to reproduce [12]. There are considerable and significant efforts to study the population dynamics in predator-prey relationships by using mathematical modeling. Lotka and Volterra initially proposed the prey-predator model [13,14]. Afterwards, prey-predator models became an important research area in applied mathematics.
Mathematical modeling has played a crucial role in understanding the dynamics of an epidemic outbreak and in developing different control measures [15][16][17]. In the past, various models have been developed to study the transmission dynamics of this parasitic disease, principally in the sheep-dog interaction, and much progress has been made over the last 30 years in modeling [6,[18][19][20][21][22][23]. To better understand the dynamics of transmission and to propose effective prevention and intervention strategies against cystic echinococcosis outbreak, an advanced study by inclusion of a human transmission component to Echinococcus granulosus is essential. In [24], a mathematical model of the transmission dynamics of alveolar echinococcosis (AE) that took into account the predator-prey interaction between the populations of fox and voles with analysis was presented. Based on the framework of this paper, we have developed a model that helps us to understand the transmission dynamics of cyst echinococcosis. As life cycle of Echinococcus granulosus' eggs involves dog, sheep, and human populations, the effect of predator-prey interaction between these populations in the transmission of cystic echinococcosis needs to be considered. We formulate and analyze a mathematical model for the transmission dynamics of cystic echinococcosis, which takes into account predator-prey interactions of dog, sheep, and human populations. We consider the populations of dog to be definite host, the populations of sheep and human to be intermediate hosts, and the concentration of parasites in the environment to be the source of infection for intermediate hosts. This paper is organized as follows. In Section 2, mathematical models of predator-prey interaction between the populations of sheep, dog, and human with analysis are presented, and cystic echinococcosis with assumptions is given in detail. In Section 3, well-posedness of the system is discussed, and in Section 4, both the disease-free and endemic equilibrium points are determined; local and global stabilities of these equilibrium points with calculation of the basic reproduction number are presented. Moreover, local and global sensitivity analyses, numerical simulations, and controlling strategies are presented in Section 5, and finally, conclusions are drawn in Section 6.

Predator-Prey
Model. The dynamics of sheep, dog, and human populations can be represented in a predator-prey relationship. We consider the case where dogs feed only on organs of sheep. Without sheep as the source of food, the dog population (N d ) would exponentially decrease and becomes extinct. However, N d is increased by a rate proportional to the number of encounters between the sheep population (N s ) and dog population (N d ). The growth rate of the dog population due to consumption of sheep is ωN s N d , where ω denotes the conversion efficiency of consumed sheep into the dog reproduction rate. Sheep are the main food source for human. In the absence of sheep, it is assumed that there exists some alternative food source for growth of the human population. With sheep as the source of food, the human population (N h ) is increased by a rate proportional to the number of encounters between the sheep population (N s ) and human population (N h ). The growth rate of the human population due to consumption of sheep is θN s N h , where θ denotes the conversion efficiency of consumed sheep into the human reproduction rate. Without sheep as the source of food, the human population N h grows exponentially at a per-capita growth rate r h but eventually increases up to the carrying capacity of the environment K h . As the number of dog and human populations increases, the population decreases in response to increased rates of consumption. The rate of consumption of sheep by dogs and humans is represented by a and c, respectively. However, in the absence of dog and human populations, the sheep population N s increases logistically. It initially increases exponentially at a per-capita growth rate r s , and it is self-limited by the carrying capacity of the environment K s . Furthermore, the sheep, dog, and human populations die naturally at rates represented by μ s , μ d , and μ h , respectively. Hence, the predator-prey interaction between the sheep, dog, and human populations is represented by a system of nonlinear ordinary differential equations: with initial conditions N s ð0Þ ≥ 0, N d ð0Þ ≥ 0, and N h ð0Þ ≥ 0, where r s > 0, a > 0, ω > 0, c > 0, μ s > 0, μ d > 0, μ h > 0, θ > 0, and r h > 0.
Equilibrium points are as follows: (1) The trivial fixed point E 0 = ð0, 0, 0Þ, which represents the extinction of the three populations (2) E 1 = ðK s ð1 − ðμ s /r s ÞÞ, 0, 0Þ, which represent the coexistence of the sheep which represent the extinction of humans and the extinction of dogs, respectively (5) Equilibrium point of the coexistence: which represents the coexistence of all the three populations 2.1.2. Stability of Equilibrium Points. In the next sections, the transmission dynamics of cystic echinococcosis is studied. Since the dog, sheep, and human populations are important for the transmission dynamics of the disease, we are mainly concerned with the equilibrium point of the coexistence. Thus, the detail analysis of the local stability conditions of the equilibrium points is presented in Appendix B, and the stability conditions of E 5 are presented in the following theorem.

Theorem 2.
(1) E 5 is locally asymptotically stable provided that a 1 > 0, a 3 > 0, and a 1 a 2 − a 3 > 0, where (2) The equilibrium point E 5 is globally asymptotically stable in the Int:ℝ 3 + Proof. The proof of this theorem is presented in Appendix B.

Modeling the Transmission Dynamics of Cystic
Echinococcosis. We model the dynamics of echinococcosis as a compartmental framework. In the human population, there are four classes: the susceptible (S h ), exposed (E h ), infectious (I h ), and removed (R h ) classes. The dog population has three classes: the susceptible (S d ), exposed (E d ), and infectious (I d ) classes. The sheep population has also three classes: the susceptible (S s ), exposed (E s ), and infectious (I s ) classes.
In the dynamics of the disease transmission, susceptible sheep are infected by ingesting parasite eggs in the feces of infected definitive hosts (dogs), while humans are infected by accidentally ingesting Echinococcus granulosus' eggs from the environment. Echinococcus granulosus' egg contamination rate in the environment by infected dogs is represented by δ. The incorporation of the saturation is relevant since for the disease to be spread from the environment to sheep, there must be a sufficient number (saturation) from the environment that can cause infection. It measures the crowding effect of the parasite in the environment, the inhibition effect of susceptible human and sheep when their number increases, and the probability of contracting the disease by sheep and human populations from the environment. Susceptible humans S h and susceptible sheep S s become infected with the disease at rates β eh B/ðχ h + BÞ and β es B/ðχ s + BÞ, respectively, where B denotes the concentration of Echinococcus granulosus' eggs in the environment, β eh denotes the rate of ingestion of Echinococcus' egg from the environment by humans, β es denotes the rate of ingestion of Echinococcus' egg from the environment by sheep, and χ h and χ s are the half-saturation constants of parasite in the environment sufficient to infect humans and sheep, respectively. Susceptible dogs are infected by preying on the infected sheep. The disease transmission rate from sheep to dogs is denoted by β sd .
The model is developed based on the following assumptions.
(1) The total population of dogs, sheep, and humans is assumed constant, and at stable equilibrium,  (3) The infected human population could recover from the disease naturally, and their recovery rate is represented by α h , where as sheep and dogs cannot recover once they are infected. We assume that there is no Echinococcus-induced death. However, dogs, sheep, and humans die naturally at rates μ d , μ s , and μ h , respectively All variables and parameters are introduced in Tables 1  and 2, and the flowchart for the transmission dynamics of the disease is shown in Figure 1.
The flowchart in Figure 1 demonstrates the interactions between the three populations and the transition of individuals from one compartment to another. The solid arrows show progression of individuals from one compartment to another. The broken line from I s to the line between S d and E d tells us that dogs become infected when they feed on organs of an infected sheep. The broken lines from B to the line between E s and S s and from B to the line between E h and S h express the fact that sheep and humans are infected by accidentally ingesting an egg of Echinococcus granulosus.
Based on these assumptions and using the transitions among different classes of disease stages given in Figure 1, the transmission dynamics of cystic echinococcosis in the populations of dogs, humans, and sheep can be expressed by the following system of ordinary differential equations.

Well-Posedness of the Solutions
In this section, we test whether the model (7)- (17) is well posed epidemiologically and mathematically in a set or not.
Proof. The detailed proof of this theorem is presented in Appendix C.

Existence and Stability of Equilibria
The equilibrium points of the system (7)- (17) are found by equating each time derivative to zero. Thus, From equations (20) and (21), we, respectively, have From equations (18) and (19) and using (29), we have

Journal of Applied Mathematics
Substituting (30) in (18) yields Similarly, from equations (26)- (28), we obtain By equating (31) and (32), we then obtain a quadratic equation: Thus, we have two roots: All the remaining state variables S h , E h , I h , R h , S s , and E s obtained from equations (22)

The Basic Reproduction Number.
The basic reproduction number is one of the fundamental concepts in mathematical biology. It is a threshold parameter, intended to quantify the spread of disease by estimating the average number of secondary infections in a wholly susceptible population, giving an indication of the invasion strength of an epidemic. We determine the basic reproduction number using the nextgeneration matrix (NGM) approach [25]. See Appendix D for the detailed computation. Therefore, the basic reproduction number, denoted by R 0 , is given by We write it as corresponds to the average number of dogs in which one infectious sheep causes the disease over its expected infection period in a completely susceptible dog population, R s = γ s β es N * s /μ sγs corresponds to the number of sheep in which an infectious dog induces the disease over its expected infection period in a completely susceptible sheep population, and R e = δ/μ e χ s corresponds to the contribution of the environment to the sheep population as a result of one infectious dog subject during its infectious period.

Stability of the DFE
Theorem 4. If R 0 < 1, then the disease-free equilibrium X 0 is locally asymptotically stable.
Theorem 5. If R 0 < 1, then the disease-free equilibrium X 0 is globally asymptotically stable in D. If R 0 > 1, then the DFE is unstable, the system is persistent, and there is at least one endemic equilibrium in the interior of D.
Proof. To prove the global stability of the disease-free equilibrium X 0 , we use a matrix-theoretic method as explained in [26].
The disease compartments of the model (7)- (17) can be written as Journal of Applied Mathematics Since V −1 F is a reducible matrix [27], the condition of Theorem 2 in [26] fails. Instead, to establish the global stability of the DFE, we construct a Lyapunov function by using Theorem 1 of [26]. Let W T = ðw 1 , w 2 , w 3 , w 4 , w 5 , w 6 , w 7 Þ be the left eigenvector of V −1 F corresponding to the eigenvalue R 0 . Thus, As a result, we found that is a Lyapunov function for the model (7)- (17). Then, differentiating along the solutions of the system (7)- (17) gives LaSalle's invariance principle [28], the disease-free equilibrium X 0 is globally asymptotically stable if R 0 < 1. Furthermore, from (41), if R 0 > 1, then Q ′ > 0 in D when B > 0 and ðS d , S s Þ = ðN * d , N * s Þ. Thus, the disease-free equilibrium X 0 is unstable, and using Theorem 2 of [26], the system (7)-(17) is uniformly persistent and hence implies that there is at least one endemic equilibrium in the interior of D.

Existence and Stability of the Endemic Equilibrium (EE).
From the result in (34), an endemic equilibrium is obtained when B > 0. Thus, the endemic equilibrium point of the system in terms of the reproduction number R 0 is given by where Proof. The detailed proof of this theorem is presented in Appendix E.

Numerics: Elasticity Indices, Numerical
Simulations, and Control Strategies 5.1. Elasticity Indices. Our analysis in previous sections demonstrated that the quantity of R 0 plays a crucial role in determining the dynamic behavior of our model. In this section, we will identify which parameters have a high impact on the basic reproduction number R 0 using data from the literature and assumed (estimated) values as given in Table 3, and we use the result for intervention strategies.
To determine the parameter that contributes most to the disease transmission, we perform local sensitivity analysis by calculating the normalized sensitivity index (elasticity index) as introduced in [31,32]. The normalized forward sensitivity index (elasticity index) of a variable (R 0 ) with respect to a parameter p is the ratio of the relative change in the variable to the relative change in the parameter, given by Table 4 gives the elasticity indices of R 0 with respect to the key parameters of the model (7)-(17) at the baseline values indicated in Table 3 and arranged in descending order 7 Journal of Applied Mathematics of magnitudes. The sign of the elasticity index tells whether R 0 increases (positive sign) or decreases (negative sign) with the parameter, whereas the magnitude determines the relative importance of the parameter. From the magnitude of the elasticity index, we can notice that four parameters (β es , β sd , δ, and χ s ) have equal and the greatest influence for the transmission of the disease, followed by γ s and γ d .

Global Sensitivity Analysis.
From the local sensitivity analysis, we observed that it is impossible to differentiate explicitly the most influential parameter(s) of the model. In order to determine which parameter(s) among the six is (are) most influential in the dynamics of the disease, global sensitivity analysis is done. We employed the technique of Latin hypercube sampling (LHS) to test the sensitivity of the model to each input parameter, as described and implemented in [33], and partial rank correlation coefficients (PRCC) to assess the significance of each parameter with respect to each metric are used. Latin hypercube sampling is a stratified sampling technique that creates sets of parameters by sampling for each parameter according to a predefined probability distribution. To examine the dependence of R 0 on parameter variations, we determine the PRCC values by considering a range of parameters as given in Table 4, with sample size 1000. The result is depicted in Figure 2. The parameter with the PRCC value far away from zero indicates more strongly the parameter influence R 0 . The negative sign for PRCC indicates inverse proportionality.
From Figure 2, it is observed that the transmission rate of Echinococcus' eggs from the environment to sheep ðβ es Þ, the contamination rate of Echinococcus' eggs in the environment by infected dogs (δ), and the transmission rate from sheep to dogs (β sd ) are the most influential parameters among the six parameters in the disease dynamics.  Table 3, we obtain a reproductive number R 0 = 0:57 < 1. Figure 3 depicts the global stability of the disease-free equilibrium as proven in Theorem 5 with different initial conditions. We can notice that all disease com- h , E * s , and I * s converge asymptotically to zero, while the noninfected compartments S * d , S * h , and S * s converge to their respective total populations. Figure 4 shows the time evolution of human, sheep, and dog populations for the model (7)-(17) with parameter values given in Table 3 by increasing β es = 0:0001 to β es = 0:001. In this case, the reproductive number is R 0 = 1:8 > 1 and depicts the global stability of the endemic equilibrium as proven in Theorem 6. It can be noticed that all the compartments of the dog, human, and sheep populations converge asymptotically to their respective endemic equilibrium points irrespective of any initial conditions.

Control Strategies.
The global sensitivity analysis presented in Sections 5.1 and 5.2 demonstrates that cystic echinococcosis can be controlled by reducing the transmission rate of Echinococcus' eggs from the environment to sheep ðβ es Þ. In this section, we illustrate the effect of the transmission rate of Echinococcus' eggs from the environment to sheep ðβ es Þ.
The effect of the transmission rate of Echinococcus' eggs from the environment to sheep (β es ) using baseline parameter values in Table 3, and when β es varied from 0.01 to 0.0001, is displayed in Figure 5. As a result, the infectious sheep, dog, and human populations are, respectively, reduced from 123 to 0, 84 to 0, and 145 to 0, where the basic reproductive number is also reduced from R 0 = 5:72 to R 0 = 0:57. The result shows that increasing the transmission rate of Echinococcus' eggs from the environment to sheep ðβ es Þ will intensify the transmission of the disease. Thus, to control the disease transmission, we suggest that it is important to plan an intervention strategy to decrease the transmission rate of Echinococcus' eggs from the environment to sheep ðβ es Þ.

Conclusion
In this paper, we first derived a mathematical model of the predator-prey interaction of the sheep, dog, and human populations in order to ascertain the conditions for the coexistence of these populations in predator-prey relationship. We then determine the equilibrium points and study the stability of each of these points. We also formulated a compartmental model for transmission dynamics of cystic  Journal of Applied Mathematics echinococcosis. We calculated both the disease-free equilibrium (DFE) and the endemic equilibrium (EE) points of the model. To study the behavior of the disease, the basic reproduction number R 0 is derived. We proved that, when R 0 < 1, the DFE is locally asymptotically stable and globally asymptotically stable, which implies that the disease dies out, whereas when R 0 > 1, the EE is globally asymptotically stable, which implies that the disease persists in all the populations. To identify which parameter has a great impact on the disease transmission, we performed both local and global sensitivity analyses on the basic reproduction number, from which we have noticed that the most sensitive parameter is the transmission rate of Echinococcus' eggs from the environment to sheep (β es ). Furthermore, we have also observed that the variation of β es has a significant effect on the number of infectious sheep, dog, and human populations. Thus, the effective controlling strategy that controls the disease transmission is decreasing the transmission rate of Echinococcus' eggs from the environment to sheep (β es ). To this effect, we suggest that the transmission dynamics of the disease by considering different intervention strategies must be studied. Therefore, extension of our work by introducing intervention strategies is an interesting future study.
Using a standard theorem of the dynamical system [34], f ðXÞ is the Lipschitz continuous. Hence, there exists a unique solution of (1)-(3) for all times t ≥ 0.    Table 3, using different initial conditions, which gives R 0 = 0:57.    (1) is given by Since N s ð0Þ ≥ 0, we can notice that N s ðtÞ ≥ 0 for all t > 0. Equation (2) is a separable differential equation, whose integral solution is given by ðA:3Þ Thus, N d ð0Þ ≥ 0 implies that N d ðtÞ ≥ 0 for all t > 0. In a similar manner, an integral solution of (3) is given by Since N h ð0Þ ≥ 0, we obtain N h ðtÞ ≥ 0 for t > 0. Therefore, the solution ðN s ðtÞ, N d ðtÞ, N h ðtÞÞ of the model (1)-(3) is nonnegative for all t ≥ 0.
To prove the boundedness of the solution, let us assume that ðN s ðtÞ, N d ðtÞ, N h ðtÞÞ is the solution of the system (1)-(3) such that TðtÞ = N s ðtÞ + N d ðtÞ + N h ðtÞ. Then, ðA:5Þ For some ε > 0, we have Now choose ε such that 0 < ε ≤ μ d ; then, we have ðA:7Þ From the solution of the differential equation, we obtain  Table 3, with varying values of β es .

Journal of Applied Mathematics
Case 2. If ω − a > 0 and θ − c > 0, then we can observe that ðA:10Þ From the solution of the differential equation, we obtain ðA:11Þ In both cases, we have 0 < TðtÞ ≤ η/ε as t → ∞. This proves that all solutions of the system are uniformly bounded.

B.1. Local Stability of Equilibrium Points
Proof. Since we are working with a first-order nonlinear system of differential equations, we can analyze the stability of our model at its equilibrium points by linearizing the system using the Jacobian matrix. The Jacobian matrix for the sys- The eigenvalues of the Jacobian are determined from the characteristic equation: where a 1 = −TrðJÞ, a 2 = A 11 + A 22 + A 33 , and a 3 = −det ðJÞ, where A 11 , A 22 , and A 33 are the minors of the entries a 11 , a 22 , and a 33 for the Jacobian J = ða ij Þ 3 . The local stability conditions of the equilibrium points obtained in Section 2.1.1 are discussed as follows.
(1) The Jacobian matrix at the equilibrium point E 0 is The eigenvalues of JðE 0 Þ are λ 1 = r s − μ s , λ 2 = −μ d , and λ 3 = r h − μ h . Hence, E 0 is stable if r s < μ s and r h < μ h . Otherwise, it is unstable (2) The Jacobian matrix at the equilibrium point E 1 is The eigenvalues of JðE 1 Þ are λ 1 = μ s − r s , λ 2 = ωK s ð1 − ðμ s /r s ÞÞ − μ d , and θK s ð1 − ðμ s /r s ÞÞ + r h − μ h . Hence, E 1 is stable if r s > μ s , r h > μ h , and K s < r s μ d /ðr s − μ s Þ and unstable if otherwise (3) The Jacobian matrix at the equilibrium point E 2 is The eigenvalues of JðE 2 Þ are λ 1 = r s − μ s − cK h ð1 − ðμ h / r h ÞÞ, λ 2 = −μ d , and λ 3 = μ h − r h . Hence, E 2 is stable if μ h < r h and μ s > r s and unstable if otherwise (4) The Jacobian matrix at the equilibrium point E 3 is ðB:6Þ Hence, E 3 is unstable since the eigenvalue λ = ððθμ d Þ/ωÞ + r h > 0.

12
Journal of Applied Mathematics (5) The Jacobian matrix at the equilibrium point E 4 is ðB:9Þ ðB:10Þ In order to verify that VðN s , N d , N h Þ > 0 for all ðN s , N d , N h Þ ≠ E 5 , it suffices to show that VðE 5 Þ is a minimum (global minimum). In this case, we apply the second-order partial test for three variables. The Hessian matrix is shown as follows: ðB:11Þ 13 Journal of Applied Mathematics with D 1 = 1/ðN * s Þ 2 > 0, ðB:12Þ and D 3 = jHðE 5 Þj > 0. This indicates that VðE 5 Þ is the local minimum. Moreover, V is a convex function on a convex set. Therefore, VðE 5 Þ is minimum, and consequently, VðN s , Thus, V is a Lyapunov function. Therefore, the equilibrium point E 5 = ðN * s , N * d , N * h Þ is globally asymptotically stable.
C. Proof of Theorem 3

C.1. Existence and Uniqueness of Solutions
Proof. The model (7)- (17) with initial conditions can be expressed as  (7)- (17). Using a standard theorem of the dynamical system [34], f ðXÞ is the Lipschitz continuous. Hence, there exists a unique solution of (7)-(17) for all times t > 0. From equations (7), (11), and (15), we, respectively, obtain We prove that the remaining variables are positive by a method of contradiction. Suppose that the conclusion is not true. Then, there exists t 1 ∈ ½0, rÞ for some r > 0 such that ðC:3Þ If eðt 1 Þ = E d ðt 1 Þ, then from (8) and since S d ðtÞ > 0 for t > 0, we obtain dE d /dt > −γ d E d for all t ∈ ½0, t 1 Þ. It then follows that which leads to a contradiction. If eðt 1 Þ = I d ðt 1 Þ, then from (9), we have dI d /dt > −μ d I d for all t ∈ ½0, t 1 Þ. Thus, which leads to a contradiction.
From continuity of the functions (the state variables), any of the variables can never be negative. Therefore, the solution of (7)-(17) is positive for all t ≥ 0.
Secondly, we prove the boundedness of the solutions as follows.

D. Calculation of the Basic Reproduction Number
According to the concepts of the next-generation matrix and reproduction number presented in [25,36], we define The Jacobian matrix of the infection subsystem at X 0 can be decomposed as F − V, where F is a matrix of transmission rates given by  Obviously, the second and third terms of dV/dt are negative. To show the global stability of endemic equilibrium X E , it suffices to show that A is Volterra-Lyapunov stable in D/ fX E g. For this purpose, we show that matrixÃ is Volterra-Lyapunov stable and matrix U = g A −1 is Volterra-Lyapunov stable (or − g A −1 is diagonally stable) [37,38]. Clearly, the diagonal elements −ð1/β sd I s + μ d Þ < 0, −ð1/ γ d Þ < 0, and det ð− g D −1 Þ > 0. Thus, based on Lemma 2.4 of [37], − g D −1 is Volterra-Lyapunov stable. Therefore, D = −Ã is diagonally stable, and hence,Ã is Volterra-Lyapunov stable.

Condition 2. To show that g
A −1 is Volterra-Lyapunov stable, we consider Based on Lemma 2.4 of [37], we can also observe that − f E −1 is also Volterra-Lyapunov stable. Therefore, g A −1 is Volterra-Lyapunov stable. Based on Lemma 2.8 of [37], there exists a diagonal matrix W = diag fw 1 , w 2 , w 3 , w 4 g such that Wð−AÞ + ð−AÞ T W T > 0 or WA + A T W T < 0, which indicates that A is Volterra-Lyapunov stable.
Therefore, dV/dt < 0, and by LaSalle's invariance principle [28], X E is globally asymptotically stable in the interior of D and is unique provided that R 0 > 1.

Data Availability
The numerical data used in our research are obtained from published literature, which are cited therein. We also use reasonable estimate, for the data that are not available in literature.