Abstract

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 and prove that the disease-free equilibrium is globally asymptotically stable, and hence, the disease dies out if . We further show that the endemic equilibrium is globally asymptotically stable, and hence, the disease persists if . 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 (). Thus, the effective controlling strategies are associated with this parameter.

1. 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 [24].

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 transmitted 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 [57].

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 [911].

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 [1517]. 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, 1823]. 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.

2. Mathematical Modeling

2.1. 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 () would exponentially decrease and becomes extinct. However, is increased by a rate proportional to the number of encounters between the sheep population () and dog population (). The growth rate of the dog population due to consumption of sheep is , 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 () is increased by a rate proportional to the number of encounters between the sheep population () and human population (). The growth rate of the human population due to consumption of sheep is , where denotes the conversion efficiency of consumed sheep into the human reproduction rate. Without sheep as the source of food, the human population grows exponentially at a per-capita growth rate but eventually increases up to the carrying capacity of the environment . 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 and , respectively. However, in the absence of dog and human populations, the sheep population increases logistically. It initially increases exponentially at a per-capita growth rate , and it is self-limited by the carrying capacity of the environment . Furthermore, the sheep, dog, and human populations die naturally at rates represented by , , and , 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 , , and , where , , , , , , , , and .

Theorem 1. The region , where and , is positively invariant (mathematically well-posed) for the model (1)–(3).

Proof. The proof of this theorem is presented in Appendix A.

2.1.1. Equilibrium Points

Our next result concerns the existence of equilibrium points of the system (1)–(3) that are biologically feasible and determines the conditions for the existence of the equilibrium of the system in the feasible region. At equilibrium, equations of the model (1)–(3) are

Equilibrium points are as follows: (1)The trivial fixed point , which represents the extinction of the three populations(2), which represent the coexistence of the sheep(3), which represent the coexistence of humans(4) and , 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 are presented in the following theorem.

Theorem 2. (1) is locally asymptotically stable provided that , , and , where(2)The equilibrium point is globally asymptotically stable in the

Proof. The proof of this theorem is presented in Appendix B.

2.2. 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 (), exposed (), infectious (), and removed () classes. The dog population has three classes: the susceptible (), exposed (), and infectious () classes. The sheep population has also three classes: the susceptible (), exposed (), and infectious () 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 and susceptible sheep become infected with the disease at rates and , respectively, where denotes the concentration of Echinococcus granulosus’ eggs in the environment, denotes the rate of ingestion of Echinococcus’ egg from the environment by humans, denotes the rate of ingestion of Echinococcus’ egg from the environment by sheep, and and 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 .

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, , so that the birth rate and death rate of each of the populations are equal(2)Dog, sheep, and human populations are recruited to the susceptible class by birth at rates , , and , respectively(3)The infected human population could recover from the disease naturally, and their recovery rate is represented by , 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 , , and , 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 to the line between and tells us that dogs become infected when they feed on organs of an infected sheep. The broken lines from to the line between and and from to the line between and 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. where , , , and .

3. 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.

Theorem 3. The region , where , , , and , is positively invariant for the model (7)–(17).

Proof. The detailed proof of this theorem is presented in Appendix C.

4. 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

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 , , , , , and obtained from equations (22)–(27) are expressed in terms of , where . Thus, we have two cases, and the results are presented in Sections 4.1 and 4.4.

4.1. Disease-Free Equilibrium (DFE)

From algebraic computation when , the system (7)–(17) has the DFE .

4.2. 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 next-generation matrix (NGM) approach [25]. See Appendix D for the detailed computation. Therefore, the basic reproduction number, denoted by , is given by

We write it as , where 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, 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 corresponds to the contribution of the environment to the sheep population as a result of one infectious dog subject during its infectious period.

4.3. Stability of the DFE

Theorem 4. If , then the disease-free equilibrium is locally asymptotically stable.

Proof. See Theorem 2 in [25].

Theorem 5. If , then the disease-free equilibrium is globally asymptotically stable in . If , then the DFE is unstable, the system is persistent, and there is at least one endemic equilibrium in the interior of .

Proof. To prove the global stability of the disease-free equilibrium , we use a matrix-theoretic method as explained in [26].
The disease compartments of the model (7)–(17) can be written as where , , and are matrices given in (D.2) and (D.3), and since , , , , and .
Matrices and are entry-wise nonnegative with Since 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 be the left eigenvector of corresponding to the eigenvalue . Thus, As a result, we found that is the left eigenvector of corresponding to the eigenvalue . Thus, by Theorem 1 of [26], is a Lyapunov function for the model (7)–(17). Then, differentiating along the solutions of the system (7)–(17) gives Since , , and , if . Furthermore, implies that and . Thus, using the results obtained in Section 4, we have , , , , , , , , and . Hence, the largest invariant set of the model where in is the singleton . Therefore, by LaSalle’s invariance principle [28], the disease-free equilibrium is globally asymptotically stable if . Furthermore, from (41), if , then in when and . Thus, the disease-free equilibrium 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 .

4.4. Existence and Stability of the Endemic Equilibrium (EE)

From the result in (34), an endemic equilibrium is obtained when . Thus, the endemic equilibrium point of the system in terms of the reproduction number is given by where

Theorem 6. If , then an endemic equilibrium defined by equations (43)–(53) is globally asymptotically stable.

Proof. The detailed proof of this theorem is presented in Appendix E.

5. Numerics: Elasticity Indices, Numerical Simulations, and Control Strategies

5.1. Elasticity Indices

Our analysis in previous sections demonstrated that the quantity of 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 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 () with respect to a parameter 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 with respect to the key parameters of the model (7)–(17) at the baseline values indicated in Table 3 and arranged in descending order of magnitudes. The sign of the elasticity index tells whether 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 (, and ) have equal and the greatest influence for the transmission of the disease, followed by and .

5.2. 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 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 . 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 , the contamination rate of Echinococcus’ eggs in the environment by infected dogs (), and the transmission rate from sheep to dogs () are the most influential parameters among the six parameters in the disease dynamics.

5.3. Numerical Simulations

In this section, numerical simulations to gain insight into some of quantitative features of the model (7)–(17) are presented.

For the model (1)–(3), we estimate the parameter values: , , , , , , , , , and . These parameter values result in an asymptotically stable equilibrium point . Thus, using the total populations , , and , an initial concentration of Echinococcus’ eggs , and the baseline parameter values given in Table 3, we obtain a reproductive number . 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 compartments , , , , , , , and converge asymptotically to zero, while the noninfected compartments , , and 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 to . In this case, the reproductive number is 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.

5.4. 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 . In this section, we illustrate the effect of the transmission rate of Echinococcus’ eggs from the environment to sheep .

The effect of the transmission rate of Echinococcus’ eggs from the environment to sheep () using baseline parameter values in Table 3, and when 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 to . The result shows that increasing the transmission rate of Echinococcus’ eggs from the environment to sheep 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 .

6. 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 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 is derived. We proved that, when , the DFE is locally asymptotically stable and globally asymptotically stable, which implies that the disease dies out, whereas when , 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 (). Furthermore, we have also observed that the variation of 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 (). 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.

Appendix

A. Proof of Theorem 1

Proof. Here, we show the existence and uniqueness of solutions, positivity, and boundedness of the solution of the model (1)–(3), since the state variables , , and represent population densities. The positivity of these variables is to conform with the reality for biological populations, while boundedness of the solution is to conform with the natural restriction to growth due to limited resources.

A.1. Existence and Uniqueness of Solutions

The system (1)–(3) with initial conditions can be expressed as where is a vector in and is the vector field in such that , , and . Using a standard theorem of the dynamical system [34], is the Lipschitz continuous. Hence, there exists a unique solution of (1)–(3) for all times .

A.2. Positivity and Boundedness of Solutions

Equation (1) is a Bernoulli type. Thus, the integral solution of equation (1) is given by

Since , we can notice that for all .

Equation (2) is a separable differential equation, whose integral solution is given by

Thus, implies that for all . In a similar manner, an integral solution of (3) is given by

Since , we obtain for . Therefore, the solution of the model (1)–(3) is nonnegative for all .

To prove the boundedness of the solution, let us assume that is the solution of the system (1)–(3) such that . Then,

For some , we have

Now choose such that ; then, we have

Case 1. If and , then From the solution of the differential equation, we obtain

Case 2. If and , then we can observe that Thus, , where .
From the solution of the differential equation, we obtain In both cases, we have as . This proves that all solutions of the system are uniformly bounded.

B. Proof of Theorem 2

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 system (1)–(3) is The eigenvalues of the Jacobian are determined from the characteristic equation: where , , and , where , , and are the minors of the entries , , and for the Jacobian .
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 is The eigenvalues of are , , and . Hence, is stable if and . Otherwise, it is unstable (2)The Jacobian matrix at the equilibrium point is The eigenvalues of are , , and . Hence, is stable if , , and and unstable if otherwise (3)The Jacobian matrix at the equilibrium point is The eigenvalues of are , , and . Hence, is stable if and and unstable if otherwise (4)The Jacobian matrix at the equilibrium point is Hence, is unstable since the eigenvalue . (5)The Jacobian matrix at the equilibrium point is where , , and . By the necessary and sufficient conditions of the Routh-Hurwitz criteria, is stable provided that , , and (6)The Jacobian matrix at the equilibrium point is By the necessary and sufficient conditions of the Routh-Hurwitz criteria, is asymptotically stable provided that , , and , where

B.2. Global Stability of

The method of Lyapunov functions enables the analysis to be extended beyond only a small region near the equilibrium point (global analysis). To study the global stability of the coexistence equilibrium point, we construct a Lyapunov function for a predator-prey model (1)–(3) (Theorem 4.2 of [35]). For an equilibrium point , let us define a function: (i)

The time derivative of is

Thus, .

In order to verify that for all , it suffices to show that 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: with , and . This indicates that is the local minimum. Moreover, is a convex function on a convex set. Therefore, is minimum, and consequently, for all . (ii)From the result in (i), we can see that whenever (iii)Clearly if

Thus, is a Lyapunov function. Therefore, the equilibrium point 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 where is a vector in and is the vector field in such that are right sides of the model (7)–(17). Using a standard theorem of the dynamical system [34], is the Lipschitz continuous. Hence, there exists a unique solution of (7)–(17) for all times .

C.2. Positivity and Boundedness of Solutions

Since the model deals with the human, sheep, and dog populations, we need to show that the state variables remain positive for all times. Suppose that is the solution of the model (7)–(17) defined for all .

From equations (7), (11), and (15), we, respectively, obtain

Thus, the variables , , and are positive for all .

We prove that the remaining variables are positive by a method of contradiction. Suppose that the conclusion is not true. Then, there exists for some such that

If , then from (8) and since for , we obtain for all . It then follows that which leads to a contradiction.

If , then from (9), we have for all . Thus, which leads to a contradiction.

If , then from (10), we have for all . Thus, which also leads to a contradiction.

Similar contradictions can be obtained if , , , , or .

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 .

Secondly, we prove the boundedness of the solutions as follows.

From equations (7)–(9), we have where . After some algebraic manipulation, the solution of this differential equation results in . This implies that . Hence, , , and are bounded.

From equations (11)–(14), we have where . Solving the differential equation results in . Thus, . Hence, , , , and are bounded.

From equations (15)–(17), we have where . In a similar manner, . Hence, , , and are bounded. Finally, from equation (10) of the model and since , we have

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 can be decomposed as , where is a matrix of transmission rates given by and is a matrix of transition rates given by

Thus, and the next-generation matrix is

The basic reproduction number is the spectral radius. Thus,

E. Proof of Theorem 6

Proof. To prove the global asymptotic stability of the endemic equilibria, we use the method of Lyapunov functions combined with the theory of Volterra-Lyapunov stable matrices. To do this, we define a Lyapunov function: Thus, The time derivative of along the solutions of model equations (7)–(17), we obtain where , , and Obviously, the second and third terms of are negative. To show the global stability of endemic equilibrium , it suffices to show that is Volterra-Lyapunov stable in . For this purpose, we show that matrix is Volterra-Lyapunov stable and matrix is Volterra-Lyapunov stable (or is diagonally stable) [37, 38].

Condition 1. To show that is Volterra-Lyapunov stable, we consider The matrix is obtained by deleting its row and column and results in Clearly, , , and . Based on Lemma 2.4 of [37], is Volterra-Lyapunov stable.
Moreover, we obtain Clearly, the diagonal elements , , and . Thus, based on Lemma 2.4 of [37], is Volterra-Lyapunov stable.
Therefore, is diagonally stable, and hence, is Volterra-Lyapunov stable.

Condition 2. To show that is Volterra-Lyapunov stable, we consider This results in Using Lemma 2.4 of [37], it is easy to observe that is Volterra-Lyapunov stable.
Moreover, we have Based on Lemma 2.4 of [37], we can also observe that is also Volterra-Lyapunov stable.
Therefore, is Volterra-Lyapunov stable.
Based on Lemma 2.8 of [37], there exists a diagonal matrix such that or , which indicates that is Volterra-Lyapunov stable.
Therefore, , and by LaSalle’s invariance principle [28], is globally asymptotically stable in the interior of and is unique provided that .

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.

Conflicts of Interest

The authors declare that there is no conflict of interests regarding the publication of this paper.