Mathematical modelling of echinococcosis in human, dogs and sheep with intervention

Abstract

In this study, a model for the spread of cyst echinococcosis with interventions is formulated. The disease-free and endemic equilibrium points of the model are calculated. The control reproduction number ${\mathcal{R}}_c$ for the model is derived, and the global dynamics are established by the values of ${\mathcal{R}}_c$. The disease-free equilibrium is globally asymptotically stable if and only if ${\mathcal{R}}_c<1$. For ${\mathcal{R}}_c>1$, using Volterra–Lyapunov stable matrices, it is proven that the endemic equilibrium is globally asymptotically stable. Sensitivity analysis to identify the most influential parameters in the dynamics of CE is carried out. To establish the long-term behaviour of the disease, numerical simulations are performed. The impact of control strategies is investigated. It is shown that, whenever vaccination of sheep is carried out solely or in combination with cleaning or disinfecting of the environment, cyst echinococcosis can be wiped out.

Keywords:

• Cystic echinococcosis
• equilibrium points
• global stability
• global sensitivity
• control strategies

1. Introduction

Cystic echinococcosis (CE) is caused by the tapeworm belonging to the family Taeniidae, called Echinococcus granulosus (for more detail, see [3,9,11,26]). Echinococcus granulosus organisms take part in both asexual reproduction and sexual reproduction. Asexual reproduction takes place via budding in the intermediate host, while sexual reproduction takes place by gamete fusion in the definitive host. Echinococcus granulosus is hermaphroditic, containing both male and female sex organs. The transmission dynamics of the disease depends on a number of factors. These include the parasite's biotic potential, stimulation of immunity in life cycle hosts, life expectancy and development time of the parasite. Social and ecological factors such as meat inspection practices, disposal of offal and casualty animals and populations of stray, feral or sylvatic hosts can all affect transmission of this parasite [26]. Consuming offal containing E. granulosus by definite hosts can lead to infection. The frequency of offal feedings and the prevalence of the parasites within the offal are factors that affect infection pressure within the definitive host. The immunity of both the definitive and intermediate host plays a large role in the transmission of the parasite, as well as the contact rate between the intermediate and the definitive host (such as herding dogs and pasture animals being kept in close proximity where dogs can contaminate grazing areas with fecal matter) [21]. The environment plays a powerful role in the transmission of infectious diseases. As that of many environmental disease, cystic echinococcosis is affected by various biological and environmental factors [12,23].

Infection with E. granulosus remains a major public health issue in several countries and regions, even in places where it was previously at low levels, as a result of a reduction of control programmes due to economic problems and lack of resources [13]. Although control programmes against human cystic (CE), caused by E. granulosus, have been established in some countries and effective control strategies are available, the parasite has still affected many countries of all continents. Thus human CE is persisting in many parts of the world with high incidences [9,13]. The human incidence can exceed 50 per 100,000 person-years in areas of endemicity, and prevalence rates as high as 5–10% can be found in some countries [14]. The incidence of human Hydatid disease in any country is closely related to the prevalence of the disease in domestic animals and is highest where there is a large dog population and high sheep production [7]. The average annual death rate from echinococcosis is 0.007 per 10,000 population, which is very low. The main causes of death are either complications of hepatic and pulmonary echinococcosis or echinococcosis of the heart. The complications of liver echinococcosis may develop due to the changes occurring not only in the parasitic cyst but also in the affected organ or in the patient's body [15].

According to WHO, cystic echinococcosis is a preventable disease [22]. To employ preventive measures, understanding of the transmission dynamics, both between dogs and sheep where the parasite maintains itself and from dogs to human is important. It is from this knowledge that effective control measures can be devised to reduce the prevalence of the parasite in animals and hence reduce the incidence of human disease. Understanding of the epidemiology of echinococcosis has been greatly improved, new diagnostic techniques for both humans and animals have been developed, new prevention strategies have emerged with the development of a vaccine against E. granulosus in intermediate hosts[9]. Since sheep has a substantial potential to transmission of parasite, vaccination of sheep with an E. granulosus recombinant antigen (EG95) offers encouraging prospects for prevention and control [5]. The vaccine is currently being produced commercially and is registered in China and Argentina. Trials in Argentina demonstrated the added value of vaccinating sheep, and in China the vaccine is being used extensively [1,7]. Currently there are no human vaccines against any form of echinococcosis [1].

Echinococcus eggs can be inactivated by disinfectants such as formalin, chlorine gas, certain freshly-prepared iodine solutions (but not most iodides) or lime can inhibit hatching of the embryo and reduce the number of viable eggs. Food safety precautions such as thorough washing of fruits and vegetables, combined with good hygiene, can reduce exposure to eggs on food. The hands should always be washed after handling pets or farming, gardening or preparing food, and before eating. Water from unsafe sources such as lakes should be boiled or filtered. Meat, particularly the intestinal tract of carnivores, should be thoroughly cooked before eating. PPE reduces the risk of infection when working with animal tissues or fecal samples and periodic surveillance of high-risk populations [10].

In [5], a mathematical model of the transmission dynamics of cyst echinococcosis without intervention was presented and analysed. From the result of sensitivity analysis, it was found that the transmission rate of E. granulosus eggs from the environment to sheep $\left({\text{β}}_{es}\right)$ is the most influential parameter that control the dynamics of disease. Thus the prevention of the disease depends on the interruption of the life cycle of E. granulosus. To break the parasite's life cycle and control the disease transmission, vaccination of sheep and/or cleaning of environment are practical alternatives.

In this paper, the model developed in [5] is extended by incorporating vaccination of sheep and disinfection or cleaning of the environment as control strategies. The objective of this paper is to determine the better control method from vaccination of sheep and disinfection or cleaning of environment.

The rest of this paper is organized as follows. In Section 2, mathematical model of cyst echinococcosis with intervention is presented. In Section 3, the positivity and boundedness of solutions are discussed. In Section 4, both the disease free and endemic equilibrium points are determined, local and global stability of these equilibrium points with the calculation of the control reproduction number are presented. Moreover, local and global sensitivity analyses, numerical simulations and effects of control strategies are presented in Section 5, and finally conclusion is drawn in Section 6.

2. Model formulation

We formulate a compartmental model to describe the transmission dynamics of the disease by considering the dog, sheep and human populations. The total populations of dog, sheep and human are assumed constant, and denoted by $N_d^{\ast }$, $N_s^{\ast }$ and $N_h^{\ast }$ respectively, so that the birth rate and death rate of each of the populations are equal. The dog population has three classes: the Susceptible ($S_d$), Exposed ($E_d$) and Infectious ($I_d$) classes. The human population has four classes: the Susceptible ($S_h$), Exposed ($E_h$), Infectious ($I_h$) and Removed ($R_h$) classes. The sheep population has four classes: the Susceptible ($S_s$), Exposed ($E_s$), Infectious ($I_s$) and Vaccinated ($V_s$) classes. Dog, sheep and human populations are recruited to susceptible class by birth at rates ${\mu }_d$, ${\mu }_s$ and ${\mu }_h$ respectively.

The transmission cycle of E. granulosus involves two hosts (dogs and sheep) and free living stages. The relative times spent in the different stages of the life cycle of the parasite and the hosts are of considerable importance to the analysis of model behaviour. The sexual reproduction of the parasites in the dog population is an important factor for the load of the parasite. The concentration of the parasite in the environment is increased by shedding from infected dog at a rate $\delta \left(P\right)$ and decreased by the natural death rate of E. granulosus eggs at rate ${\mu }_e$ and disinfection or cleaning of environment at rate μ, where $P\left(t\right)$ denote the mean number of parasite per dog host. Although modelling infectious diseases that are caused by parasites is best formulated using the intensity of infection which measures the number of parasites rather than the prevalence of infection [18], in this work, we considered the latter. Here, the incidence function represents sufficient number of parasites that can cause infection as reported in [30]. Thus the time evolution of the parasite egg's is represented by the following differential equation: ${\displaystyle \frac{\mathrm{d}B}{\mathrm{d}t}}=\delta \left(P\right)I_d\left(t\right)-{\mu }_{e}B\left(t\right),\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\delta \text{is constant.}$ 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 E. granulosus eggs from the environment. We introduced vaccination to the susceptible sheep population at a rate ν, so that the population of susceptible sheep is reduced through vaccination and moved to vaccinated $V_s$ class. We assume that vaccination is not lifelong. The sheep population may lose of vaccine-induced immunity and move back to susceptible class at a rate ρ. Rate of infection of susceptible sheep is $\frac{{\text{β}}_{es}B}{{\chi }_s+B}$, where ${\text{β}}_{es}$ denotes the rate of ingestion of Echinococcus egg from the environment by sheep and ${\chi }_s$ is the half-saturation constant of parasite in the environment sufficient to infect sheep. Rate of infection of susceptible humans is $\frac{{\text{β}}_{eh}B}{{\chi }_h+B}$, where ${\text{β}}_{eh}$ denotes the rate of ingestion of Echinococcus egg from the environment by human, and ${\chi }_h$ is the half-saturation constant of parasite in the environment sufficient to infect human. Susceptible dogs are infected by preying on the infected sheep. The disease transmission rate from sheep to dogs is denoted by ${\text{β}}_{sd}$. The rates at which exposed dog, sheep and human progress to infectious classes are denoted by ${\gamma }_d$, ${\gamma }_s$ and ${\gamma }_h$ respectively. Infected human population could recover from the disease naturally, at rate ${\alpha }_h$, whereas 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 ${\mu }_d$, ${\mu }_s$ and ${\mu }_h$ respectively. The density of E. granulosus eggs depends mainly on the number of infectious dogs. Its concentration in the environment is increased by shedding of a parasite from infected dog at a rate δ and decreased by the natural death rate of E. granulosus eggs at rate ${\mu }_e$ and disinfection or cleaning of environment at rate μ.

The general structure of the model is captured by the flow chart displayed in Figure 1.

Figure 1. The flow diagram for cyst echinococcosis transmission dynamics.

The transmission dynamics of the disease in the three populations can now be expressed by the following system of first-order differential equations: (1) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{S}_d}{\mathrm{d}t}}={\mu }_d{N}_d^{\ast }-{\text{β}}_{sd}{I}_s{S}_d-{\mu }_d{S}_d,\end{array}$(1) (2) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{E}_d}{\mathrm{d}t}}={\text{β}}_{sd}{I}_s{S}_d-({\mu }_d+{\gamma }_d)E_d,\end{array}$(2) (3) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{I}_d}{\mathrm{d}t}}={\gamma }_d{E}_d-{\mu }_d{I}_d,\end{array}$(3) (4) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}B}{\mathrm{d}t}}=\delta I_d-({\mu }_e+\mu )B,\end{array}$(4) (5) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{S}_h}{\mathrm{d}t}}={\mu }_h{N}_h^{\ast }-{\mu }_h{S}_h-{\displaystyle \frac{{\text{β}}_{eh}B}{{\chi }_h+B}}{S}_h,\end{array}$(5) (6) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{E}_h}{\mathrm{d}t}}={\displaystyle \frac{{\text{β}}_{eh}B}{{\chi }_h+B}}{S}_h-({\mu }_h+{\gamma }_h)E_h,\end{array}$(6) (7) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{I}_h}{\mathrm{d}t}}={\gamma }_h{E}_h-({\mu }_h+{\alpha }_h)I_h,\end{array}$(7) (8) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{R}_h}{\mathrm{d}t}}={\alpha }_h{I}_h-{\mu }_h{R}_h,\end{array}$(8) (9) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{S}_s}{\mathrm{d}t}}={\mu }_s{N}_s^{\ast }+\rho V_s-{\mu }_s{S}_s-{\displaystyle \frac{{\text{β}}_{es}B}{{\chi }_s+B}}{S}_s-\nu S_s\end{array}$(9) (10) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{V}_s}{\mathrm{d}t}}=\nu S_s-\rho V_s-{\mu }_{s}Vs,\end{array}$(10) (11) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{E}_s}{\mathrm{d}t}}={\displaystyle \frac{{\text{β}}_{es}B}{{\chi }_s+B}}{S}_s-({\mu }_s+{\gamma }_s)E_s,\end{array}$(11) (12) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{I}_s}{\mathrm{d}t}}={\gamma }_s{E}_s-{\mu }_s{I}_s,\end{array}$(12) with initial conditions $S_d\left(0\right)\ge 0$, $E_d\left(0\right)\ge 0$, $I_d\left(0\right)\ge 0$, $B\left(0\right)\ge 0$, $S_h\left(0\right)\ge 0$, $E_h\left(0\right)\ge 0$, $I_h\left(0\right)\ge 0$, $R_h\left(0\right)\ge 0$, $S_s\left(0\right)\ge 0$, $E_s\left(0\right)\ge 0$, $I_s\left(0\right)\ge 0$ and $V_s\left(0\right)\ge 0$.

For convenience, we make the following substitutions: $\widetilde{{\gamma }_d}={\mu }_d+{\gamma }_d$, $\widetilde{{\gamma }_h}={\mu }_h+{\gamma }_h$, $\widetilde{q_h}={\mu }_h+{\alpha }_h$ and $\widetilde{{\gamma }_s}={\mu }_s+{\gamma }_s$. Thus we can rewrite the system (1(1) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{S}_d}{\mathrm{d}t}}={\mu }_d{N}_d^{\ast }-{\text{β}}_{sd}{I}_s{S}_d-{\mu }_d{S}_d,\end{array}$(1) )–(12(12) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{I}_s}{\mathrm{d}t}}={\gamma }_s{E}_s-{\mu }_s{I}_s,\end{array}$(12) ) with initial conditions as (13) ${\displaystyle \frac{\mathrm{d}X}{\mathrm{d}t}}=g(t,X(t\left)\right),X\left(0\right)=X_0$(13) where $X\left(t\right)=(S_d,E_d,I_d,B,S_h,E_h,I_h,R_h,S_s,V_s,E_s,I_s{)}^T$, $g:{\mathbb{R}}_{+}\times {\mathbb{R}}_{+}^{12}\to {\mathbb{R}}_{+}^{12}$ defined by $g(t,X(t\left)\right)=\left(\begin{array}{c}{\mu }_{\mathrm{d}}{N}_d^{\ast }-{\text{β}}_{sd}{I}_s{S}_d-{\mu }_d{S}_d\\ {\text{β}}_{sd}{I}_s{S}_d-\widetilde{{\gamma }_d}{E}_d\\ {\gamma }_d{E}_d-{\mu }_d{I}_d\\ \delta I_d-({\mu }_e+\mu )B\\ {\mu }_h{N}_h^{\ast }-{\mu }_h{S}_h-{\displaystyle \frac{{\text{β}}_{eh}B}{{\chi }_h+B}}{S}_h\\ {\displaystyle \frac{{\text{β}}_{eh}B}{{\chi }_h+B}}{S}_h-\widetilde{{\gamma }_h}{E}_h\\ {\gamma }_h{E}_h-\widetilde{q_h}{I}_h\\ {\alpha }_h{I}_h-{\mu }_h{R}_h\\ {\mu }_s{N}_s^{\ast }+\rho V_s-{\mu }_s{S}_s-{\displaystyle \frac{{\text{β}}_{es}B}{{\chi }_s+B}}{S}_s-\nu S_s\\ \nu S_s-\rho V_s-{\mu }_{s}Vs\\ {\displaystyle \frac{{\text{β}}_{es}B}{{\chi }_s+B}}{S}_s-\widetilde{{\gamma }_s}{E}_s\\ {\gamma }_s{E}_s-{\mu }_s{I}_s\end{array}\right)$

3. Well-posedness of the solutions

Before we proceed with the mathematical analysis, we need to show that the models (1(1) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{S}_d}{\mathrm{d}t}}={\mu }_d{N}_d^{\ast }-{\text{β}}_{sd}{I}_s{S}_d-{\mu }_d{S}_d,\end{array}$(1) )–(12(12) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{I}_s}{\mathrm{d}t}}={\gamma }_s{E}_s-{\mu }_s{I}_s,\end{array}$(12) ) (alternatively model (13(13) ${\displaystyle \frac{\mathrm{d}X}{\mathrm{d}t}}=g(t,X(t\left)\right),X\left(0\right)=X_0$(13) )) is well-posed epidemiologically and mathematically in a feasible domain.

Theorem 3.1

The region $\mathbb{D}={\mathbb{D}}_1\times {\mathbb{D}}_2\times {\mathbb{D}}_3\times {\mathbb{D}}_4$, where ${\mathbb{D}}_1=\left\{\right(S_d,E_d,I_d)\subset {\mathbb{R}}_{+}^3:S_d+E_d+I_d\le N_d^{\ast }\}$, ${\mathbb{D}}_2=\left\{\right(S_h,E_h,I_h,R_h)\subset {\mathbb{R}}_{+}^4:S_h+E_h+I_h+R_h\le N_h^{\ast }\}$, ${\mathbb{D}}_3=\left\{\right(S_s,E_s,I_s)\subset {\mathbb{R}}_{+}^3:S_s+E_s+I_s+V_s\le N_s^{\ast }\}$ and ${\mathbb{D}}_4=\{B\subset {\mathbb{R}}_{+}:B\le \frac{\delta N_d^{\ast }}{{\mu }_e+\mu }\}$ is positively invariant for models (1(1) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{S}_d}{\mathrm{d}t}}={\mu }_d{N}_d^{\ast }-{\text{β}}_{sd}{I}_s{S}_d-{\mu }_d{S}_d,\end{array}$(1) )–(12(12) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{I}_s}{\mathrm{d}t}}={\gamma }_s{E}_s-{\mu }_s{I}_s,\end{array}$(12) ).

Proof.

Proof of this theorem is presented in Appendix 1

4. Existence and stability of equilibria

The equilibrium point(s) of the system (1(1) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{S}_d}{\mathrm{d}t}}={\mu }_d{N}_d^{\ast }-{\text{β}}_{sd}{I}_s{S}_d-{\mu }_d{S}_d,\end{array}$(1) )–(12(12) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{I}_s}{\mathrm{d}t}}={\gamma }_s{E}_s-{\mu }_s{I}_s,\end{array}$(12) ) are obtained by equating the right hand sides to zero as follows: (14a) $\begin{array}{rl}& {\mu }_{\mathrm{d}}{N}_d^{\ast }-{\text{β}}_{sd}{I}_s{S}_d-{\mu }_d{S}_d=0,\end{array}$(14a) (14b) $\begin{array}{rl}& {\text{β}}_{sd}{I}_s{S}_d-\widetilde{{\gamma }_d}{E}_d=0,\end{array}$(14b) (14c) $\begin{array}{rl}& {\gamma }_{\mathrm{d}}{E}_d-{\mu }_d{I}_d=0,\end{array}$(14c) (14d) $\begin{array}{rl}& \delta I_d-({\mu }_e+\mu )B=0,\end{array}$(14d) (14e) $\begin{array}{rl}& {\mu }_h{N}_h^{\ast }-{\mu }_h{S}_h-{\displaystyle \frac{{\text{β}}_{eh}B}{{\chi }_h+B}}{S}_h=0,\end{array}$(14e) (14f) $\begin{array}{rl}& {\displaystyle \frac{{\text{β}}_{eh}B}{{\chi }_h+B}}{S}_h-\widetilde{{\gamma }_h}{E}_h=0,\end{array}$(14f) (14g) $\begin{array}{rl}& {\gamma }_h{E}_h-\widetilde{q_h}{I}_h=0,\end{array}$(14g) (14h) $\begin{array}{rl}& {\alpha }_h{I}_h-{\mu }_h{R}_h=0,\end{array}$(14h) (14i) $\begin{array}{rl}& {\mu }_s{N}_s^{\ast }+\rho V_s-{\mu }_s{S}_s-{\displaystyle \frac{{\text{β}}_{es}B}{{\chi }_s+B}}{S}_s-\nu S_s=0,\end{array}$(14i) (14j) $\begin{array}{rl}& {\displaystyle \frac{{\text{β}}_{es}B}{{\chi }_s+B}}{S}_s-\widetilde{{\gamma }_s}{E}_s=0,\end{array}$(14j) (14k) $\begin{array}{rl}& {\gamma }_s{E}_s-{\mu }_s{I}_s=0,\end{array}$(14k) (14l) $\begin{array}{rl}& \nu S_s-\rho V_s-{\mu }_s{V}_s=0,\end{array}$(14l) From Equations (14c(14c) $\begin{array}{rl}& {\gamma }_{\mathrm{d}}{E}_d-{\mu }_d{I}_d=0,\end{array}$(14c) ) and (14d(14d) $\begin{array}{rl}& \delta I_d-({\mu }_e+\mu )B=0,\end{array}$(14d) ), we respectively have (15a) $I_d={\displaystyle \frac{({\mu }_e+\mu )B}{\delta }}\mathrm{a}\mathrm{n}\mathrm{d}{E}_d={\displaystyle \frac{{\mu }_d({\mu }_e+\mu )B}{\delta {\gamma }_d}}.$(15a) From Equations (14a(14a) $\begin{array}{rl}& {\mu }_{\mathrm{d}}{N}_d^{\ast }-{\text{β}}_{sd}{I}_s{S}_d-{\mu }_d{S}_d=0,\end{array}$(14a) ), (14b(14b) $\begin{array}{rl}& {\text{β}}_{sd}{I}_s{S}_d-\widetilde{{\gamma }_d}{E}_d=0,\end{array}$(14b) ) and using (15a(15a) $I_d={\displaystyle \frac{({\mu }_e+\mu )B}{\delta }}\mathrm{a}\mathrm{n}\mathrm{d}{E}_d={\displaystyle \frac{{\mu }_d({\mu }_e+\mu )B}{\delta {\gamma }_d}}.$(15a) ) we have (15b) $S_d=N_d^{\ast }-{\displaystyle \frac{({\mu }_e+\mu ))\widetilde{{\gamma }_d}B}{\delta {\gamma }_d}}.$(15b) Substituting (15b(15b) $S_d=N_d^{\ast }-{\displaystyle \frac{({\mu }_e+\mu ))\widetilde{{\gamma }_d}B}{\delta {\gamma }_d}}.$(15b) ) in (14a(14a) $\begin{array}{rl}& {\mu }_{\mathrm{d}}{N}_d^{\ast }-{\text{β}}_{sd}{I}_s{S}_d-{\mu }_d{S}_d=0,\end{array}$(14a) ) yields (15c) $I_s={\displaystyle \frac{1}{{\text{β}}_{sd}}}\left({\displaystyle \frac{({\mu }_e+\mu ){\mu }_d\widetilde{{\gamma }_d}B}{\delta {\gamma }_d{N}_d^{\ast }-\widetilde{{\gamma }_d}({\mu }_e+\mu )B}}\right).$(15c) Similarly, from Equations (14i(14i) $\begin{array}{rl}& {\mu }_s{N}_s^{\ast }+\rho V_s-{\mu }_s{S}_s-{\displaystyle \frac{{\text{β}}_{es}B}{{\chi }_s+B}}{S}_s-\nu S_s=0,\end{array}$(14i) )–(14l(14l) $\begin{array}{rl}& \nu S_s-\rho V_s-{\mu }_s{V}_s=0,\end{array}$(14l) ), we obtain (15d) $I_s={\displaystyle \frac{{\gamma }_s{\text{β}}_{es}(\rho +{\mu }_s)N_s^{\ast }B}{\widetilde{{\gamma }_s}\left[\right({\text{β}}_{es}(\rho +{\mu }_s)+{\mu }_s(\rho +{\mu }_s+\nu ))B+{\mu }_s{\chi }_s(\rho +{\mu }_s+\nu \left)\right]}}.$(15d) By Equating (15c(15c) $I_s={\displaystyle \frac{1}{{\text{β}}_{sd}}}\left({\displaystyle \frac{({\mu }_e+\mu ){\mu }_d\widetilde{{\gamma }_d}B}{\delta {\gamma }_d{N}_d^{\ast }-\widetilde{{\gamma }_d}({\mu }_e+\mu )B}}\right).$(15c) ) and (15d(15d) $I_s={\displaystyle \frac{{\gamma }_s{\text{β}}_{es}(\rho +{\mu }_s)N_s^{\ast }B}{\widetilde{{\gamma }_s}\left[\right({\text{β}}_{es}(\rho +{\mu }_s)+{\mu }_s(\rho +{\mu }_s+\nu ))B+{\mu }_s{\chi }_s(\rho +{\mu }_s+\nu \left)\right]}}.$(15d) ), we then obtain a quadratic equation (16) $d{B}^2-eB=0,$(16) where $d=\widetilde{{\gamma }_d}({\mu }_e+\mu )\left[{\mu }_d\widetilde{{\gamma }_s}\right({\text{β}}_{es}(\rho +{\mu }_s)+{\mu }_s(\rho +{\mu }_s+\nu ))+{\gamma }_s{\text{β}}_{es}{\text{β}}_{sd}{N}_s^{\ast }(\rho +{\mu }_s\left)\right]$ and $e=\delta {\gamma }_d{\gamma }_s(\rho +{\mu }_s){\text{β}}_{es}{\text{β}}_{sd}{N}_s^{\ast }{N}_d^{\ast }-({\mu }_e+\mu ){\mu }_s{\chi }_s{\mu }_d\widetilde{{\gamma }_d}\widetilde{{\gamma }_s}(\rho +{\mu }_s+\nu ).$ Thus we have two roots (17) $\begin{array}{rl}B& =0,\\ B& ={\displaystyle \frac{\delta {\gamma }_d{\gamma }_s(\rho +{\mu }_s){\text{β}}_{es}{\text{β}}_{sd}{N}_s^{\ast }{N}_d^{\ast }-{\mu }_s{\chi }_s{\mu }_d\widetilde{{\gamma }_d}\widetilde{{\gamma }_s}({\mu }_e+\mu )(\rho +{\mu }_s+\nu )}{\widetilde{{\gamma }_d}({\mu }_e+\mu )\left[{\mu }_d\widetilde{{\gamma }_s}\right({\text{β}}_{es}(\rho +{\mu }_s)+{\mu }_s(\rho +{\mu }_s+\nu ))+{\gamma }_s{\text{β}}_{es}{\text{β}}_{sd}(\rho +{\mu }_s\left)N_s^{\ast }\right]}}.\end{array}$(17) The results presented in Sections 4.1 and 4.4 follow from these two roots.

4.1. Disease-free equilibrium (DFE)

From algebraic computation when B = 0, the system (13(13) ${\displaystyle \frac{\mathrm{d}X}{\mathrm{d}t}}=g(t,X(t\left)\right),X\left(0\right)=X_0$(13) ) has the DFE given by $\begin{array}{rl}{X}^0& =(S_d^0,E_d^0,I_d^0,B^0,S_h^0,E_h^0,I_h^0,R_h^0,S_s^0,E_s^0,I_s^0,V_s^0)\\ & =(N_d^{\ast },0,0,0,N_h^{\ast },0,0,0,{\displaystyle \frac{(\rho +{\mu }_s)N_s^{\ast }}{\rho +{\mu }_s+\nu }},0,0,{\displaystyle \frac{\nu N_s^{\ast }}{\rho +{\mu }_s+\nu }}).\end{array}$

4.2. The control reproduction number

A key quantity in epidemiological models is the reproduction number. It is a useful threshold in the study of a disease for predicting outbreak and for evaluating the control strategies. Due to the presence of controls in the model (13(13) ${\displaystyle \frac{\mathrm{d}X}{\mathrm{d}t}}=g(t,X(t\left)\right),X\left(0\right)=X_0$(13) ), the term ‘the control reproduction number’ is used. The control reproduction number, denoted by ${\mathcal{R}}_c$ represents the average number of secondary infections caused by an infectious individual over the course of infectious period in a totally susceptible population under specified controls. We derive the control reproduction number ${\mathcal{R}}_c$ by using the Next Generation Matrix (NGM) approach [27] on the system (13(13) ${\displaystyle \frac{\mathrm{d}X}{\mathrm{d}t}}=g(t,X(t\left)\right),X\left(0\right)=X_0$(13) ). The detail computation is done in Appendix A. Thus the control reproduction number is given by ${\mathcal{R}}_c=\rho \left(F{V}^{-1}\right)=\sqrt{{\displaystyle \frac{\rho +{\mu }_s}{\rho +{\mu }_s+\nu }}{\displaystyle \frac{{\gamma }_s{\text{β}}_{es}{N}_s^{\ast }}{{\chi }_s({\mu }_e+\mu )}}{\displaystyle \frac{{\text{β}}_{sd}}{{\mu }_s\widetilde{{\gamma }_s}}}{\displaystyle \frac{\delta {\gamma }_d{N}_d^{\ast }}{{\mu }_d\widetilde{{\gamma }_d}}}}=\sqrt{{\displaystyle \frac{\rho +{\mu }_s}{\rho +{\mu }_s+\nu }}{\displaystyle \frac{{\mu }_e}{{\mu }_e+\mu }}}{\mathcal{R}}_0.$ where ${\mathcal{R}}_0=\sqrt{\frac{{\gamma }_s{\text{β}}_{es}{N}_s^{\ast }}{{\chi }_s{\mu }_e}{\displaystyle \frac{{\text{β}}_{sd}}{{\mu }_s\widetilde{{\gamma }_s}}}\frac{\delta {\gamma }_d{N}_d^{\ast }}{{\mu }_d\widetilde{{\gamma }_d}}}$ is the basic reproduction number as derived in [5], which represents the average number of secondary infections caused by an infectious individual over the course of infectious period in a totally susceptible population without vaccination and disinfection or cleaning of the environment. Here, we can notice that ${\mathcal{R}}_c<{\mathcal{R}}_0$.

4.3. Stability of the DFE

Theorem 4.1

If ${\mathcal{R}}_c\le 1$, then the disease free equilibrium $X^0$ is globally asymptotically stable in $\mathbb{D}$. If ${\mathcal{R}}_c>1$, then the DFE is unstable, the system is persistent and there is at least one equilibrium in the interior of $\mathbb{D}$.

Proof.

To prove the global stability of the disease free equilibrium $X^0$, we use a matrix-theoretic method as explained in [25].

The disease compartments of model (13(13) ${\displaystyle \frac{\mathrm{d}X}{\mathrm{d}t}}=g(t,X(t\left)\right),X\left(0\right)=X_0$(13) ) can be written as $x^{\prime }=(F-V)x-f(x,y)$ where $x=(E_d,I_d,B,E_h,I_h,E_s,I_s{)}^T$, $y=(S_d,S_h,R_h,S_s,V_s{)}^T$, F and V are matrices given in (A1a(A1a) $F=\left(\begin{array}{ccccccc}0& 0& 0& 0& 0& 0& {\text{β}}_{sd}{N}_d^{\ast }\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& {\displaystyle \frac{{\text{β}}_{eh}{N}_h^{\ast }}{{\chi }_h}}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& {\displaystyle \frac{{\text{β}}_{es}(\rho +{\mu }_s)N_s^{\ast }}{(\rho +{\mu }_s+\nu ){\chi }_s}}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\end{array}\right),$(A1a) ) and (A1b(A1b) $V=\left(\begin{array}{ccccccc}\widetilde{{\gamma }_d}& 0& 0& 0& 0& 0& 0\\ -{\gamma }_d& {\mu }_d& 0& 0& 0& 0& 0\\ 0& -\delta & {\mu }_e+\mu & 0& 0& 0& 0\\ 0& 0& 0& \widetilde{{\gamma }_h}& 0& 0& 0\\ 0& 0& 0& -{\gamma }_h& \widetilde{q_h}& 0& 0\\ 0& 0& 0& 0& 0& \widetilde{{\gamma }_s}& 0\\ 0& 0& 0& 0& 0& -{\gamma }_s& {\mu }_s\end{array}\right).$(A1b) ), and $f(x,y)=\left(\begin{array}{c}{\text{β}}_{sd}{I}_s(N_d^{\ast }-S_d)\\ 0\\ 0\\ {\displaystyle \frac{{\text{β}}_{eh}B}{{\chi }_h}}(N_h^{\ast }-{\displaystyle \frac{{\chi }_h{S}_h}{B+{\chi }_h}})\\ 0\\ {\displaystyle \frac{{\text{β}}_{es}B}{{\chi }_s({\chi }_s+B)}}(S_s^{0}B+{\chi }_s(S_s^0-S_s\left)\right)\\ 0\end{array}\right)\ge 0,$ since $S_d\le N_d^{\ast }$, $S_h\le N_h^{\ast }$, ${\chi }_h\le {\chi }_h+B$, and $S_s\le S_s^0$.

Matrices F and $V^{-1}$ are entry wise non negative with $\begin{array}{r}{V}^{-1}F=\left(\begin{array}{ccccccc}0& 0& 0& 0& 0& 0& {\displaystyle \frac{{\text{β}}_{sd}{N}_d^{\ast }}{\widetilde{{\gamma }_d}}}\\ 0& 0& 0& 0& 0& 0& {\displaystyle \frac{{\gamma }_d{\text{β}}_{sd}{N}_d^{\ast }}{{\mu }_d\widetilde{{\gamma }_d}}}\\ 0& 0& 0& 0& 0& 0& {\displaystyle \frac{\delta {\gamma }_d{\text{β}}_{sd}{N}_d^{\ast }}{({\mu }_e+\mu ){\mu }_d\widetilde{{\gamma }_d}}}\\ 0& 0& {\displaystyle \frac{{\text{β}}_{eh}{N}_h^{\ast }}{{\chi }_h\widetilde{{\gamma }_h}}}& 0& 0& 0& 0\\ 0& 0& {\displaystyle \frac{{\gamma }_h{\text{β}}_{eh}{N}_h^{\ast }}{{\chi }_h\widetilde{{\gamma }_h}\widetilde{q_h}}}& 0& 0& 0& 0\\ 0& 0& {\displaystyle \frac{{\text{β}}_{es}(\rho +{\mu }_s)N_s^{\ast }}{{\chi }_s\widetilde{{\gamma }_s}(\rho +{\mu }_s+\nu )}}& 0& 0& 0& 0\\ 0& 0& {\displaystyle \frac{{\gamma }_s{\text{β}}_{es}(\rho +{\mu }_s)N_s^{\ast }}{{\chi }_s\widetilde{{\gamma }_s}{\mu }_s(\rho +{\mu }_s+\nu )}}& 0& 0& 0& 0\end{array}\right)\end{array}$ Since $V^{-1}F$ is reducible matrix [4], the condition of Theorem 2.2 in [25] fails. Instead, to establish the global stability of the DFE, we construct a Lyapunov function by using Theorem 2.1 of [25].

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 ${\mathcal{R}}_c$. Thus $W^T{V}^{-1}F={\mathcal{R}}_{c}W.$ As a result, we found that $W=(0,0,1,0,0,0,\frac{\delta {\gamma }_d{\text{β}}_{sd}{N}_d^{\ast }}{{\mu }_d\widetilde{{\gamma }_d}({\mu }_e+\mu ){\mathcal{R}}_c})$ is the left eigenvector of $V^{-1}F$ corresponding to the eigenvalue ${\mathcal{R}}_c$. Thus by Theorem 2.1 of [25], $\begin{array}{rl}Q& =W^T{V}^{-1}x={\displaystyle \frac{\delta {\gamma }_d{E}_d}{{\mu }_d({\mu }_e+\mu )\widetilde{{\gamma }_d}}}+{\displaystyle \frac{\delta I_d}{{\mu }_d({\mu }_e+\mu )}}+{\displaystyle \frac{B}{({\mu }_e+\mu )}}\\ & +({\displaystyle \frac{{\gamma }_s}{{\mu }_s\widetilde{{\gamma }_s}}}{E}_s+{\displaystyle \frac{I_s}{{\mu }_s}}){\displaystyle \frac{\delta {\gamma }_d{\text{β}}_{sd}{N}_d^{\ast }}{{\mu }_d({\mu }_e+\mu )\widetilde{{\gamma }_d}{\mathcal{R}}_c}},\end{array}$ is a Lyapunov function for model (1(1) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{S}_d}{\mathrm{d}t}}={\mu }_d{N}_d^{\ast }-{\text{β}}_{sd}{I}_s{S}_d-{\mu }_d{S}_d,\end{array}$(1) )–(12(12) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{I}_s}{\mathrm{d}t}}={\gamma }_s{E}_s-{\mu }_s{I}_s,\end{array}$(12) ). Then differentiating along the solutions of the system (1(1) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{S}_d}{\mathrm{d}t}}={\mu }_d{N}_d^{\ast }-{\text{β}}_{sd}{I}_s{S}_d-{\mu }_d{S}_d,\end{array}$(1) )–(12(12) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{I}_s}{\mathrm{d}t}}={\gamma }_s{E}_s-{\mu }_s{I}_s,\end{array}$(12) ) gives (18) $\begin{array}{rl}{Q}^{\mathrm{\prime }}& =({\mathcal{R}}_c-1)({\displaystyle \frac{\delta {\gamma }_d{\text{β}}_{sd}{I}_s{N}_d^{\ast }}{{\mu }_d({\mu }_e+\mu )\widetilde{{\gamma }_d}{\mathcal{R}}_c}}+B)-\left({\displaystyle \frac{\delta {\gamma }_d{\text{β}}_{sd}{I}_s}{{\mu }_d({\mu }_e+\mu )\widetilde{{\gamma }_d}}}\right(N_d^{\ast }-S_d)\\ & +{\displaystyle \frac{{\gamma }_s{\gamma }_d{\text{β}}_{es}{\text{β}}_{sd}B{N}_d^{\ast }}{{\mu }_s({\mu }_e+\mu ){\mu }_d{\chi }_s\widetilde{{\gamma }_d}{\mathcal{R}}_c}}(S_s^0-{\displaystyle \frac{S_s{\chi }_s}{B+{\chi }_s}})).\end{array}$(18) Since $S_d\le N_d^{\ast }$, $S_s\le S_s^0$ and ${\chi }_s\le {\chi }_s+B$, $Q^{\mathrm{\prime }}<0$ if ${\mathcal{R}}_c\le 1$. If $Q^{\mathrm{\prime }}=0$, then B = 0 and $I_s=0$. Hence, the largest invariant set of the model where $Q^{\prime }=0$ in $int\left(\mathbb{D}\right)$ is the singleton $\left\{X^0\right\}$. Therefore, by Lasalle's invariance principle [16], the disease free equilibrium $X^0$ is globally asymptotically stable if ${\mathcal{R}}_c<1$. For ${\mathcal{R}}_c>1$, the first term in (18(18) $\begin{array}{rl}{Q}^{\mathrm{\prime }}& =({\mathcal{R}}_c-1)({\displaystyle \frac{\delta {\gamma }_d{\text{β}}_{sd}{I}_s{N}_d^{\ast }}{{\mu }_d({\mu }_e+\mu )\widetilde{{\gamma }_d}{\mathcal{R}}_c}}+B)-\left({\displaystyle \frac{\delta {\gamma }_d{\text{β}}_{sd}{I}_s}{{\mu }_d({\mu }_e+\mu )\widetilde{{\gamma }_d}}}\right(N_d^{\ast }-S_d)\\ & +{\displaystyle \frac{{\gamma }_s{\gamma }_d{\text{β}}_{es}{\text{β}}_{sd}B{N}_d^{\ast }}{{\mu }_s({\mu }_e+\mu ){\mu }_d{\chi }_s\widetilde{{\gamma }_d}{\mathcal{R}}_c}}(S_s^0-{\displaystyle \frac{S_s{\chi }_s}{B+{\chi }_s}})).\end{array}$(18) ) is positive in $\mathbb{D}$, and the remaining terms are equal to zero when B = 0 and $(S_d,S_s)=(S_d^0,S_s^0)$, where $S_d^0=N_d^{\ast }$ and $S_s^0=\frac{(\rho +{\mu }_s)N_s^{\ast }}{\rho +{\mu }_s+\nu }$. Consequently, $Q^{\mathrm{\prime }}>0$ in a neighbourhood of $X^0$. Thus the disease free equilibrium $X^0$ is unstable, and using Theorem 2.2 of [25], the system (13(13) ${\displaystyle \frac{\mathrm{d}X}{\mathrm{d}t}}=g(t,X(t\left)\right),X\left(0\right)=X_0$(13) ) is uniformly persistent and hence imply there is at least one endemic equilibrium in the interior of $\mathbb{D}$.

4.4. Existence and stability of the endemic equilibrium (EE)

From (17(17) $\begin{array}{rl}B& =0,\\ B& ={\displaystyle \frac{\delta {\gamma }_d{\gamma }_s(\rho +{\mu }_s){\text{β}}_{es}{\text{β}}_{sd}{N}_s^{\ast }{N}_d^{\ast }-{\mu }_s{\chi }_s{\mu }_d\widetilde{{\gamma }_d}\widetilde{{\gamma }_s}({\mu }_e+\mu )(\rho +{\mu }_s+\nu )}{\widetilde{{\gamma }_d}({\mu }_e+\mu )\left[{\mu }_d\widetilde{{\gamma }_s}\right({\text{β}}_{es}(\rho +{\mu }_s)+{\mu }_s(\rho +{\mu }_s+\nu ))+{\gamma }_s{\text{β}}_{es}{\text{β}}_{sd}(\rho +{\mu }_s\left)N_s^{\ast }\right]}}.\end{array}$(17) ), it follows that the model (13(13) ${\displaystyle \frac{\mathrm{d}X}{\mathrm{d}t}}=g(t,X(t\left)\right),X\left(0\right)=X_0$(13) ) admits an endemic equilibrium. Thus the endemic equilibrium point of the system in terms of the control reproduction number ${\mathcal{R}}_c$ is given by $X_E=(S_d^{\ast },E_d^{\ast },I_d^{\ast },B^{\ast },S_h^{\ast },E_h^{\ast },I_h^{\ast },R_h^{\ast },S_s^{\ast },E_s^{\ast },I_s^{\ast },V_s^{\ast }),$ where (19a) $\begin{array}{rl}{I}_d^{\ast }& ={\displaystyle \frac{({\mu }_e+\mu )B^{\ast }}{\delta }}\\ & ={\displaystyle \frac{{\chi }_s{\mu }_d{\mu }_s\widetilde{{\gamma }_s}({\mu }_e+\mu )(\rho +{\mu }_s+\nu )}{\delta {\mu }_d\widetilde{{\gamma }_s}\left[{\text{β}}_{es}\right(\rho +{\mu }_s)+{\mu }_s(\rho +{\mu }_s+\nu \left)\right]+\delta {\gamma }_s{\text{β}}_{sd}{\text{β}}_{es}(\rho +{\mu }_s)N_s^{\ast }}}({\mathcal{R}}_c^2-1)\end{array}$(19a) (19b) $\begin{array}{rl}& E_d^{\ast }={\displaystyle \frac{{\chi }_s{\mu }_d^2{\mu }_s\widetilde{{\gamma }_s}({\mu }_e+\mu )(\rho +{\mu }_s+\nu )}{\delta {\mu }_d{\gamma }_d\widetilde{{\gamma }_s}\left[{\text{β}}_{es}\right(\rho +{\mu }_s)+{\mu }_s(\rho +{\mu }_s+\nu \left)\right]+\delta {\gamma }_s{\gamma }_d{\text{β}}_{sd}{\text{β}}_{es}(\rho +{\mu }_s)N_s^{\ast }}}({\mathcal{R}}_c^2-1)\end{array}$(19b) (19c) $\begin{array}{rl}& S_d^{\ast }=N_d^{\ast }-{\displaystyle \frac{({\mu }_e+\mu )\widetilde{{\gamma }_d}B}{\delta {\gamma }_d}}\end{array}$(19c) (19d) $\begin{array}{rl}& S_h^{\ast }={\displaystyle \frac{{\mu }_h{N}_h^{\ast }({\chi }_h+B^{\ast })}{{\mu }_h({\chi }_h+B^{\ast })+{\text{β}}_{eh}{B}^{\ast }}}\end{array}$(19d) (19e) $\begin{array}{rl}& E_h^{\ast }={\displaystyle \frac{{\mu }_h{\text{β}}_{eh}{N}_h^{\ast }{B}^{\ast }}{\widetilde{{\gamma }_h}\left({\mu }_h\right({\chi }_h+B^{\ast })+{\text{β}}_{eh}{B}^{\ast })}}\end{array}$(19e) (19f) $\begin{array}{rl}& I_h^{\ast }={\displaystyle \frac{{\mu }_h{\gamma }_h{\text{β}}_{eh}{N}_h^{\ast }{B}^{\ast }}{\widetilde{q_h}\widetilde{{\gamma }_h}\left({\mu }_h\right({\chi }_h+B^{\ast })+{\text{β}}_{eh}{B}^{\ast })}}\end{array}$(19f) (19g) $\begin{array}{rl}& R_h^{\ast }={\displaystyle \frac{{\gamma }_h{\alpha }_h{\text{β}}_{eh}{N}_h^{\ast }{B}^{\ast }}{\widetilde{q_h}\widetilde{{\gamma }_h}\left({\mu }_h\right({\chi }_h+B^{\ast })+{\text{β}}_{eh}{B}^{\ast })}}\end{array}$(19g) (19h) $\begin{array}{rl}& S_s^{\ast }={\displaystyle \frac{{\mu }_s(\rho +{\mu }_s)({\chi }_s+B^{\ast })N_s^{\ast }}{\left({\text{β}}_{es}\right(\rho +{\mu }_s)+{\mu }_s(\rho +{\mu }_s+\nu \left)\right)B^{\ast }+{\mu }_s{\chi }_s(\rho +{\mu }_s+\nu )}}\end{array}$(19h) (19i) $\begin{array}{rl}& E_s^{\ast }={\displaystyle \frac{{\mu }_s{\text{β}}_{es}(\rho +{\mu }_s)N_s^{\ast }{B}^{\ast }}{\widetilde{{\gamma }_s}\left[\right({\text{β}}_{es}(\rho +{\mu }_s)+{\mu }_s(\rho +{\mu }_s+\nu ))B^{\ast }+{\mu }_s{\chi }_s(\rho +{\mu }_s+\nu \left)\right]}}\end{array}$(19i) (19j) $\begin{array}{rl}& I_s^{\ast }={\displaystyle \frac{{\gamma }_s{\text{β}}_{es}(\rho +{\mu }_s)N_s^{\ast }{B}^{\ast }}{\widetilde{{\gamma }_s}\left[\right({\text{β}}_{es}(\rho +{\mu }_s)+{\mu }_s(\rho +{\mu }_s+\nu ))B^{\ast }+{\mu }_s{\chi }_s(\rho +{\mu }_s+\nu \left)\right]}}\end{array}$(19j) (19k) $\begin{array}{rl}& V_s^{\ast }={\displaystyle \frac{\nu {\mu }_s({\chi }_s+B^{\ast })N_s^{\ast }}{\left({\text{β}}_{es}\right(\rho +{\mu }_s)+{\mu }_s(\rho +{\mu }_s+\nu \left)\right)B^{\ast }+{\mu }_s{\chi }_s(\rho +{\mu }_s+\nu )}}\end{array}$(19k) (19l) $\begin{array}{rl}& B^{\ast }={\displaystyle \frac{{\chi }_s{\mu }_d{\mu }_s\widetilde{{\gamma }_s}(\rho +{\mu }_s+\nu )}{{\mu }_d\widetilde{{\gamma }_s}\left[{\text{β}}_{es}\right(\rho +{\mu }_s)+{\mu }_s(\rho +{\mu }_s+\nu \left)\right]+{\gamma }_s{\text{β}}_{sd}{\text{β}}_{es}(\rho +{\mu }_s)N_s^{\ast }}}({\mathcal{R}}_c^2-1)\end{array}$(19l) For ${\mathcal{R}}_c=1$, it can be noted that the endemic equilibrium point reduces to disease free equilibrium point. Under the condition ${\mathcal{R}}_c<1$, the quadratic equation (17(17) $\begin{array}{rl}B& =0,\\ B& ={\displaystyle \frac{\delta {\gamma }_d{\gamma }_s(\rho +{\mu }_s){\text{β}}_{es}{\text{β}}_{sd}{N}_s^{\ast }{N}_d^{\ast }-{\mu }_s{\chi }_s{\mu }_d\widetilde{{\gamma }_d}\widetilde{{\gamma }_s}({\mu }_e+\mu )(\rho +{\mu }_s+\nu )}{\widetilde{{\gamma }_d}({\mu }_e+\mu )\left[{\mu }_d\widetilde{{\gamma }_s}\right({\text{β}}_{es}(\rho +{\mu }_s)+{\mu }_s(\rho +{\mu }_s+\nu ))+{\gamma }_s{\text{β}}_{es}{\text{β}}_{sd}(\rho +{\mu }_s\left)N_s^{\ast }\right]}}.\end{array}$(17) ) has no positive root. Hence, the model equation (13(13) ${\displaystyle \frac{\mathrm{d}X}{\mathrm{d}t}}=g(t,X(t\left)\right),X\left(0\right)=X_0$(13) ) has no positive endemic equilibrium whenever ${\mathcal{R}}_c<1$. This consequently indicates that the backward bifurcation phenomenon does not occur whenever ${\mathcal{R}}_c<1$. On the contrary, the disease will persist if ${\mathcal{R}}_c$ exceeds unity, where a stable endemic equilibrium exists. The phenomenon, where the disease-free equilibrium loses its stability and a stable endemic equilibrium appears as ${\mathcal{R}}_c$ increases through one, is known as forward bifurcation.

Theorem 4.2

If ${\mathcal{R}}_c>1$, then an endemic equilibrium $X_E=(S_d^{\ast },E_d^{\ast },I_d^{\ast },B^{\ast },S_h^{\ast },E_h^{\ast },I_h^{\ast },R_h^{\ast },S_s^{\ast },E_s^{\ast },I_s^{\ast },V_s^{\ast })$ defined by Equations (19a(19a) $\begin{array}{rl}{I}_d^{\ast }& ={\displaystyle \frac{({\mu }_e+\mu )B^{\ast }}{\delta }}\\ & ={\displaystyle \frac{{\chi }_s{\mu }_d{\mu }_s\widetilde{{\gamma }_s}({\mu }_e+\mu )(\rho +{\mu }_s+\nu )}{\delta {\mu }_d\widetilde{{\gamma }_s}\left[{\text{β}}_{es}\right(\rho +{\mu }_s)+{\mu }_s(\rho +{\mu }_s+\nu \left)\right]+\delta {\gamma }_s{\text{β}}_{sd}{\text{β}}_{es}(\rho +{\mu }_s)N_s^{\ast }}}({\mathcal{R}}_c^2-1)\end{array}$(19a) )–(19l(19l) $\begin{array}{rl}& B^{\ast }={\displaystyle \frac{{\chi }_s{\mu }_d{\mu }_s\widetilde{{\gamma }_s}(\rho +{\mu }_s+\nu )}{{\mu }_d\widetilde{{\gamma }_s}\left[{\text{β}}_{es}\right(\rho +{\mu }_s)+{\mu }_s(\rho +{\mu }_s+\nu \left)\right]+{\gamma }_s{\text{β}}_{sd}{\text{β}}_{es}(\rho +{\mu }_s)N_s^{\ast }}}({\mathcal{R}}_c^2-1)\end{array}$(19l) ) is globally asymptotically stable, and the bifurcation of endemic equilibrium point is forward when ${\mathcal{R}}_c>1$.

Proof.

Detailed proof of this theorem is presented in Appendix 3.

5. Numerics: elasticity indices, numerical simulations and control strategies

5.1. Elasticity indices

In this section, we carried out sensitivity analysis to determine the model robustness to parameter values. This is a tool to identify the most influential parameters in determining model dynamics. Sensitivity analysis is used to obtain the sensitivity index that is a measure of the relative change in a state variable when a parameter changes. We compute the sensitivity indices of ${\mathcal{R}}_c$ to the model parameters with the approach used by Chitnis et al. [6,20]. These indices could be computed numerically so as to figure out parameters that have high impact on control reproduction number ${\mathcal{R}}_c$, and the importance of each individual parameter in the disease transmission dynamics and prevalence. To perform the local sensitivity analysis, we use the normalized forward sensitivity index of a variable with respect to a parameter which is expressed as the ratio of the relative variation in the variable to the relative variation in the parameter. Thus the normalized forward sensitivity index (elasticity index) of a variable (${\mathcal{R}}_c$) 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 ${\mathrm{{\rm Y}}}_p^{{\mathcal{R}}_c}={\displaystyle \frac{\mathrm{\partial }{\mathcal{R}}_c}{\mathrm{\partial }p}}\times {\displaystyle \frac{p}{{\mathcal{R}}_c}}$ We use data from the literature, and due to the great variation of some parameter values from region to region, without deviation from the range of parameter values obtained from literature, assumed (estimated) values were used for sensitivity analysis, as given in Table 1. The parameter values were hypothetically chosen, because the intention for this work is not to validate the model results of the real situation obtained in a particular study, but for illustrative purposes only. With sensitivity analysis, we can get insight into the appropriate intervention strategies to prevent and control the spread of the disease.

Table 1. Parameters, baseline values and sources.

Parameters	Description	Baseline values	Sources
${\mu }_d$	Natural death rate=birth rate of dog	0.00088/day${}^{-1}$	[26]
${\text{β}}_{sd}$	The transmission rate from sheep to dogs	0.00001/day${}^{-1}$	Assumed
${\gamma }_d$	The rates at which exposed dog progress to infective classes	0.0001/day${}^{-1}$	Assumed
${\mu }_h$	Natural death rate=birth rate of human	0.00004/day${}^{-1}$	[29]
${\text{β}}_{eh}$	Transmission rate of Echinococcos eggs from the environment
	to human	0.0001/day${}^{-1}$	Assumed
${\gamma }_h$	The rates at which exposed human progress to infective classes	0.00019/day${}^{-1}$	[28]
${\alpha }_h$	Human recovery rate	0.0001/day${}^{-1}$	Assumed
${\mu }_s$	Natural death rate=birth rate of sheep	0.00114/day${}^{-1}$	[26]
${\text{β}}_{es}$	Transmission rate of Echinococcus eggs from the environment
	to sheep	0.001/day${}^{-1}$	Assumed
${\mu }_e$	Natural death rate of E. granulosus eggs	0.0039/day${}^{-1}$	[26]
${\gamma }_s$	The rates at which exposed sheep progress to infective classes	0.0001/day${}^{-1}$	[2]
δ	Eggs contamination rate of the environment by infected dogs	0.003 cells/ml.day	Assumed
${\chi }_h$	The concentration of E. granulosus eggs at which
	half of all contacts with human produce infection	40  cells/ml	Assumed
${\chi }_s$	The concentration of E. granulosus eggs at which
	half of all contacts with sheep produce infection	20 cells/ml	Assumed
ν	Rate of vaccination of sheep	0.005/day  [0,1]	Assumed
ρ	Rate of loss of vaccine-induced immunity of sheep	0.0001/day [0,1]	Assumed
μ	Rate of disinfection or cleaning of environment	0.001/day [0,1]	Assumed

Sensitivity indices of the control reproduction number ${\mathcal{R}}_c$ with respect to the model nine parameters are determined and presented in the table below. Ranges of these parameters are taken without deviation to the data obtained in the literature and assumed values. Natural birth rates and death rates of the populations are not considered for sensitivity analysis since these parameters are difficult to control.

Table 2 gives the elasticity indices of ${\mathcal{R}}_c$ with respect to key parameters of model (1(1) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{S}_d}{\mathrm{d}t}}={\mu }_d{N}_d^{\ast }-{\text{β}}_{sd}{I}_s{S}_d-{\mu }_d{S}_d,\end{array}$(1) )–(12(12) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{I}_s}{\mathrm{d}t}}={\gamma }_s{E}_s-{\mu }_s{I}_s,\end{array}$(12) ) at the baseline values indicated in Table 1 and arranged in descending order of magnitudes. The sign of the elasticity index tells whether ${\mathcal{R}}_c$ increases (positive sign) or decreases (negative sign) with the parameter, whereas the magnitude determines the relative importance of the parameter. From the magnitude of elasticity index, we can notice that four parameters $({\text{β}}_{es},{\text{β}}_{sd},\delta ,{\chi }_s)$ have equal and the greatest influence for the transmission of the disease, followed by ${\gamma }_s$, ${\gamma }_d$ and ν, and ρ has the least influence for the transmission of the disease.

Table 2. Elasticity indices of ${\mathcal{R}}_c$ relative to some model parameters.

Parameters	Formula$({\mathrm{{\rm Y}}}_p^{R_c}={\displaystyle \frac{\mathrm{\partial }{\mathcal{R}}_c}{\mathrm{\partial }p}}\times {\displaystyle \frac{p}{{\mathcal{R}}_c}})$	Baseline value	Range	Elasticity index
${\text{β}}_{es}$	${\displaystyle \frac{1}{2}}$	0.001	0.0001−0.002	0.5
${\text{β}}_{sd}$	${\displaystyle \frac{1}{2}}$	0.00001	0.000001−0.00002	0.5
δ	${\displaystyle \frac{1}{2}}$	0.003	0.0003−0.006	0.5
${\chi }_s$	$-{\displaystyle \frac{1}{2}}$	20	10−40	−0.5
${\gamma }_s$	${\displaystyle \frac{{\mu }_s}{2({\mu }_s+{\gamma }_s)}}$	0.0001	0.00001−0.0002	0.46
${\gamma }_d$	${\displaystyle \frac{{\mu }_d}{2({\mu }_d+{\gamma }_d)}}$	0.0001	0.00001−0.0002	0.45
ν	$-{\displaystyle \frac{\nu }{2(\rho +\nu +{\mu }_s)}}$	0.005	0.0005−0.01	−0.4
μ	$-{\displaystyle \frac{\mu }{2({\mu }_e+\mu )}}$	0.001	0.0001−0.002	−0.1
ρ	${\displaystyle \frac{\rho \nu }{2(\rho +{\mu }_s)(\rho +\nu +{\mu }_s)}}$	0.0001	0.00001−0.0002	0.032

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. To determine which parameter(s) among the nine 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 [19], and Partial Rank Correlation Coefficients (PRCCs) to assess the significance of each parameter with respect to each metric is 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 ${\mathcal{R}}_c$ on parameter variations, we determine the PRCC values by considering a range of parameters as given in Table 2, with sample size 1000. The result is depicted in Figure 2. The parameter with PRCC value far away from zero indicates the more strongly the parameter influence ${\mathcal{R}}_c$. The negative sign for PRCCs indicates inverse proportionality.

Figure 2. Global sensitivity analysis displaying the partial rank correlation coefficients (PRCC) of control reproduction number ${\mathcal{R}}_c$.

From Figure 2, it is observed that the transmission rate from sheep to dog (${\text{β}}_{sd}$) and Echinococcos eggs contamination rate of the environment by infected dogs (δ) are the most influential parameters among the eight parameters in the disease dynamics. On the other hand, rate of losing vaccine induced immunity of sheep are the least sensitive parameter for the dynamics of the disease.

5.3. Numerical simulations

In this section, we carry out numerical simulations for mathematical model of cyst echinococcosis in the populations of sheep, dogs and humans. We use the total populations $N_d^{\ast }=590$, $N_h^{\ast }=1014$, $N_s^{\ast }=1760$, and parameter values given in Table 1. This yields a control reproduction number ${\mathcal{R}}_c=0.72<1$. Using different initial conditions, the time evolution of human, sheep and dog populations for model (13(13) ${\displaystyle \frac{\mathrm{d}X}{\mathrm{d}t}}=g(t,X(t\left)\right),X\left(0\right)=X_0$(13) ) is displayed in Figure 3. We can notice that all disease compartments $E_d^{\ast }$, $I_d^{\ast }$, $B^{\ast }$, $E_h^{\ast }$, $I_h^{\ast }$, $R_h^{\ast }$, $E_s^{\ast }$ and $I_s^{\ast }$ converge asymptotically to zero, while the noninfected compartments $S_d^{\ast }$, $S_h^{\ast }$ and $S_s^{\ast }+V_s^{\ast }$ converge to their respective total populations. This asserts the global stability of the disease free equilibrium as proved in Theorem 4.1.

Figure 3. Time evolution of the dog, sheep and human populations with baseline parameter values as in Table 1, using different initial conditions which gives ${\mathcal{R}}_c=0.72$.

Figure 4 shows that the time evolution of human, sheep and dog populations for model (13(13) ${\displaystyle \frac{\mathrm{d}X}{\mathrm{d}t}}=g(t,X(t\left)\right),X\left(0\right)=X_0$(13) ) with parameter values given in Table 1 by increasing ${\text{β}}_{sd}=0.00001$ to ${\text{β}}_{sd}=0.0001$. In this case, the control reproductive number is ${\mathcal{R}}_c=2.28>1$, and depict the global stability of the endemic equilibrium as proved in Theorem 4.2. 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.

Figure 4. Time evolution of the dog, sheep and human populations with baseline parameter values as in Table 1, using different initial conditions, except for ${\text{β}}_{sd}=0.0001$ which gives ${\mathcal{R}}_c=2.28$, and with approximate equilibrium values $S_d^{\ast }=254$, $E_d^{\ast }=302$, $I_d^{\ast }=34$, $S_h^{\ast }=545$, $E_h^{\ast }=81$, $I_h^{\ast }=111$, $R_h^{\ast }=276$, $S_s^{\ast }=321$, $E_s^{\ast }=132$, $I_s^{\ast }=12$, $V_s^{\ast }=1295$.

In the case of endemicity, the prevalence rate of human is $\begin{array}{rl}& {\displaystyle \frac{\mathrm{N}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r}\mathrm{o}\mathrm{f}\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{e}\mathrm{a}\mathrm{s}\mathrm{e}\mathrm{d}\mathrm{u}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{s}\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{f}\mathrm{i}\mathrm{e}\mathrm{d}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{d}}{\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{g}\mathrm{e}\mathrm{p}\mathrm{o}\mathrm{p}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{z}\mathrm{e}\times \mathrm{d}\mathrm{u}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}\mathrm{f}\mathrm{f}\mathrm{o}\mathrm{l}\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{u}\mathrm{p}}}={\displaystyle \frac{111+81}{1014\times 1}}\times 100\mathrm{\%}=18.9\mathrm{\%}\end{array}$ The prevalence rate resulted in the numerical simulation is higher than the WHO published data, since the WHO report showed that the human prevalence rate is 5–10% (as indicated in [14]). This result has shown at least 8.9% discrepancy from WHO published data.

5.4. Effects of control strategies on ${\mathcal{R}}_c$

The numerical simulations are performed to illustrate the effect of vaccination of sheep and cleaning or disinfecting the environment in the dynamics of the disease transmission in the populations of sheep, dogs and humans while they are used alone or simultaneously.

The effect of vaccination of sheep using baseline parameter values in Table 1 except for $\mu =0$, and when ν is varied from $0.005/10$ to 0.005, is displayed in Figure 5. As a result the infectious sheep, dog and human populations are respectively reduced from 25 to 0, 13 to 0, and 80 to 0, where the control reproduction number is also reduced from ${\mathcal{R}}_c=1.53$ to ${\mathcal{R}}_c=0.8$. This result shows that increasing the rate of vaccination of sheep (ν) reduces the time evolution of infected human, sheep and dog populations. The effect of disinfection or cleaning the environment using baseline parameter values in Table 1 except for $\nu =0$, and when μ is varied from $0.001/{10}^2$ to 0.001, is displayed in Figure 6. As a result the infectious sheep, dog and human populations are respectively reduced from 37 to 32, 18 to 16 and 93 to 79, where the control reproduction number is also reduced from ${\mathcal{R}}_c=1.8$ to ${\mathcal{R}}_c=1.6$. This result shows that increasing the rate of disinfection or cleaning the environment (μ) has very less effect to eradicate the disease transmission in human, sheep and dog populations.

Figure 5. The numerical simulations displaying effects of vaccination of sheep only on cumulative number of infectious dog, human and sheep populations, using parameter values in Table 1, with varying values of ν ($\mu =0$).

Figure 6. The numerical simulations displaying effects of disinfection or cleaning the environment only on cumulative number of infectious dog, human and sheep populations, using parameter values in Table 1, with varying values of μ ($\nu =0$).

Numerical simulation is also performed to illustrate the effect of vaccination of sheep and disinfection or cleaning the environment when the two control strategies are administered simultaneously, as displayed in Figure 7. The number of infectious sheep, dog and human populations are respectively reduced from 26 to 0, 13 to 0 and 80 to 0, where the control reproduction number is also reduced from ${\mathcal{R}}_c=1.5$ to ${\mathcal{R}}_c=0.72$. One can observe that these combined effects allow to reduce the size of infected individuals. Thus increasing the vaccination rate of sheep alone or a simultaneous increase of the vaccination rate of sheep and rate of cleaning or disinfecting the environment is an effective control measure of cyst trichinosis.

Figure 7. The numerical simulations displaying effects of combining controlling strategies on cumulative number of infectious dog, human and sheep populations, using parameter values in Table 1.

Furthermore, we assess the impact of combined control strategies using contour plots of ${\mathcal{R}}_c$ as function of the control strategies, and with varying rate of transmission from sheep to dog (${\text{β}}_{sd}$), we estimate the least values of the two control parameters which will ensure the disease eradication in the populations.

Figure 8(a) shows contour curves of ${\mathcal{R}}_c$ as a function of ν and μ using baseline parameter values in Table 1. We can observe that, low control strategies are needed to ensure the eradication of the parasites, with range of ${\mathcal{R}}_c\in [0.076,0.96]$ and mean $0.518$. Contour plots of ${\mathcal{R}}_c$ as function of the control strategies, for rate of transmission from sheep to dog (${\text{β}}_{sd}=0.0001$) is displayed in Figure 8(b). The least values of ν and μ that will ensure parasites eradication are estimated to be 0.01 and 0.01 so that ${\mathcal{R}}_c=1$. In this case, the combined control strategies have effect in the disease transmission with range of ${\mathcal{R}}_c\in [0.24,3.03]$ and mean $1.635$. Contour plots of ${\mathcal{R}}_c$ as function of the control strategies, for rate of transmission from sheep to dog (${\text{β}}_{sd}=0.0004$) is displayed in Figure 8(c).The least values of ν and μ that will ensure parasites eradication are estimated to be 0.045 and 0.01 so that ${\mathcal{R}}_c=0.99$. In this case, the combined control strategies has effect in the disease transmission with range of ${\mathcal{R}}_c\in [0.48,6.06]$ and mean $3.27$. These results have shown that a remarkable increase in the control reproduction number is observed with an increase rate of transmission from sheep to dog (${\text{β}}_{sd}$). Hence, to ensure the eradication of parasites, we must introduce the least values of the controls that can bring the value ${\mathcal{R}}_c<1$ with respect to the rate of transmission from sheep to dog (${\text{β}}_{sd}$).

Figure 8. Contour curves of ${\mathcal{R}}_c$ as a function of ν and μ using different rate of transmission from sheep to dog (${\text{β}}_{sd}$) with (a) parameter values in Table 1, (b) ${\text{β}}_{sd}=0.0001$, (c) ${\text{β}}_{sd}=0.0004$.

6. Conclusion

In this paper, we proposed and analysed a deterministic model for transmission dynamics of cyst trichinosis that incorporates two control strategies namely vaccination of sheep and disinfection or cleaning the environment. The model has a Disease Free Equilibrium point (DEF) which is both locally and globally asymptotically stable whenever the control reproduction number ${\mathcal{R}}_c<1$. We also found the Endemic Equilibrium points (I) and proved that it is globally stable whenever the control reproduction number ${\mathcal{R}}_c>1$. Moreover, we have performed sensitivity analysis on the control reproduction number with the two control strategies, from which we have noted that the most sensitive parameters are the transmission rate from sheep to dog (${\text{β}}_{sd}$) and Echinococcus eggs contamination rate of the environment by infected dogs (δ). Numerical simulations of the model have shown that whenever the control strategies are carried out solely then vaccination of sheep is the better alternative to eradicate cyst echinococcosis, but when disinfection or cleaning the environment is carried out solely, its effect to eradicated the disease is very less. This indicates that to eradicate the disease from the three populations the model needs to incorporate other possible intervention strategies that can reduce ${\text{β}}_{sd}$. Our finding has shown that the two control strategies are not enough to control the disease. We suggest that more controls, which focus more on dog population, should be incorporated to eradicate the spread of cyst echinococcosis. In our future, we will extend the model by incorporating additional control(s) and investigate the effectiveness and cost effectiveness of the control measures.

Disclosure statement

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

Data availability

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

Appendices

Appendix 1. Proof of Theorem 3.1

Proof.

1. Existence and uniqueness of solutions:

   Since g has a continuous derivative on any compact subset of ${\mathbb{R}}^{12}$, $g(t,X(t\left)\right)$ is the Lipschitz continuous in some open ball containing $X_0$. Hence, by the Picard–Lindelof theorem, there exist a unique local solution of (13(13) ${\displaystyle \frac{\mathrm{d}X}{\mathrm{d}t}}=g(t,X(t\left)\right),X\left(0\right)=X_0$(13) ) for all times $t\ge 0$ [24].

2. Positivity and boundedness of solutions

   Suppose that $\left(S_d\right(t),E_d(t),I_d(t),B(t),S_h(t),E_h(t),I_h(t),R_h(t),S_s(t),E_s(t),I_s(t),V_s(t\left)\right)$ be the solution of the model (1(1) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{S}_d}{\mathrm{d}t}}={\mu }_d{N}_d^{\ast }-{\text{β}}_{sd}{I}_s{S}_d-{\mu }_d{S}_d,\end{array}$(1) )–(12(12) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{I}_s}{\mathrm{d}t}}={\gamma }_s{E}_s-{\mu }_s{I}_s,\end{array}$(12) ) defined for all $t\ge 0$.

   From Equations (1(1) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{S}_d}{\mathrm{d}t}}={\mu }_d{N}_d^{\ast }-{\text{β}}_{sd}{I}_s{S}_d-{\mu }_d{S}_d,\end{array}$(1) ) and (9(9) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{S}_s}{\mathrm{d}t}}={\mu }_s{N}_s^{\ast }+\rho V_s-{\mu }_s{S}_s-{\displaystyle \frac{{\text{β}}_{es}B}{{\chi }_s+B}}{S}_s-\nu S_s\end{array}$(9) ), we respectively obtain $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{S}_d}{\mathrm{d}t}}>-({\text{β}}_{sd}{I}_s+{\mu }_d)S_d\Rightarrow S_d\left(t\right)>S_d\left(0\right)e^{-{\int }_0^t({\text{β}}_{sd}{I}_s+{\mu }_d)\mathrm{d}\tau }>0,\\ & {\displaystyle \frac{\mathrm{d}{S}_h}{\mathrm{d}\tau }}>-({\mu }_h+{\displaystyle \frac{{\text{β}}_{eh}B}{{\chi }_h+B}})S_h\Rightarrow S_h\left(t\right)>S_h\left(0\right)e^{-{\int }_0^t({\mu }_h+{\displaystyle \frac{{\text{β}}_{eh}B}{{\chi }_h+B}})\mathrm{d}\tau }>0,\end{array}$ Thus the variables $S_d\left(t\right)$ and $S_h\left(t\right)$ are positive for all $t\ge 0$.

   We prove the remaining variables are positive by a method of contradiction. Suppose that the conclusion is not true. Then, there exists $t_1\in [0,r)$ for some r>0 such that $e\left(t_1\right)=min\left\{\right(E_d\left(t_1\right),I_d\left(t_1\right),B\left(t_1\right),E_h\left(t_1\right),I_h\left(t_1\right),R_h\left(t_1\right),S_s\left(t_1\right),E_s\left(t_1\right),I_s\left(t_1\right),V_s\left(t_1\right)\left)\right\}=0$ If $e\left(t_1\right)=E_d\left(t_1\right)$, then $I_s\left(t\right)\ge 0$ for $t\in [0,t_1]$. Thus from (2(2) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{E}_d}{\mathrm{d}t}}={\text{β}}_{sd}{I}_s{S}_d-({\mu }_d+{\gamma }_d)E_d,\end{array}$(2) ) we obtain $\frac{\mathrm{d}{E}_d}{\mathrm{d}t}>-\widetilde{{\gamma }_d}{E}_d$ for all $t\in [0,t_1)$. It then follows that $0=E_d\left(t_1\right)>E_d\left(0\right)e^{{\int }_0^{t_1}-\widetilde{{\gamma }_d}\mathrm{d}t}>0$ which leads to a contradiction.

   If $e\left(t_1\right)=I_d\left(t_1\right)$, then from (3(3) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{I}_d}{\mathrm{d}t}}={\gamma }_d{E}_d-{\mu }_d{I}_d,\end{array}$(3) ) we have $\frac{\mathrm{d}{I}_d}{\mathrm{d}t}>-{\mu }_d{I}_d$ for all $t\in [0,t_1)$. Thus $0=I_d\left(t_1\right)>I_d\left(0\right)e^{{\int }_0^{t_1}-{\mu }_d\mathrm{d}t}>0$ which leads to a contradiction.

   If $e\left(t_1\right)=B\left(t_1\right)$, then from (4(4) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}B}{\mathrm{d}t}}=\delta I_d-({\mu }_e+\mu )B,\end{array}$(4) ), we have $\frac{\mathrm{d}B}{\mathrm{d}t}>-({\mu }_e+\mu )B$ for all $t\in [0,t_1)$. Thus $0=B\left(t_1\right)>B\left(0\right)e^{{\int }_0^{t_1}-({\mu }_e+\mu )dt}>0$ which also leads to a contradiction.

   Similar contradictions can be obtained if $e\left(t_1\right)=E_h\left(t_1\right)$, $e\left(t_1\right)=I_h\left(t_1\right)$, $e\left(t_1\right)=R_h\left(t_1\right)$, $e\left(t_1\right)=S_s\left(t_1\right)$, $e\left(t_1\right)=E_s\left(t_1\right)$, $e\left(t_1\right)=I_s\left(t_1\right)$ or $e\left(t_1\right)=V_s\left(t_1\right)$.

   From continuity of the functions (the state variables), any of the variables can never be negative. Therefore, the solution of (1(1) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{S}_d}{\mathrm{d}t}}={\mu }_d{N}_d^{\ast }-{\text{β}}_{sd}{I}_s{S}_d-{\mu }_d{S}_d,\end{array}$(1) )–(12(12) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{I}_s}{\mathrm{d}t}}={\gamma }_s{E}_s-{\mu }_s{I}_s,\end{array}$(12) ) is positive for all $t\ge 0$.

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

   From Equations (1(1) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{S}_d}{\mathrm{d}t}}={\mu }_d{N}_d^{\ast }-{\text{β}}_{sd}{I}_s{S}_d-{\mu }_d{S}_d,\end{array}$(1) )–(3(3) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{I}_d}{\mathrm{d}t}}={\gamma }_d{E}_d-{\mu }_d{I}_d,\end{array}$(3) ), we have ${\displaystyle \frac{\mathrm{d}{N}_d}{\mathrm{d}t}}={\mu }_d{N}_d^{\ast }-{\mu }_d{N}_d,$ where $N_d=S_d+E_d+I_d$. After some algebraic manipulation, the solution of this differential equation results in $N_d\left(t\right)=N_d^{\ast }+\left(N_d\right(0)-N_d^{\ast })e^{-{\mu }_{d}t}$. This implies that $\underset{t\to \mathrm{\infty }}{lim sup}{N}_d\left(t\right)\le N_d^{\ast }$. Hence, $S_d\left(t\right)$, $E_d\left(t\right)$ and $I_d\left(t\right)$ are bounded.

   From Equations (5(5) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{S}_h}{\mathrm{d}t}}={\mu }_h{N}_h^{\ast }-{\mu }_h{S}_h-{\displaystyle \frac{{\text{β}}_{eh}B}{{\chi }_h+B}}{S}_h,\end{array}$(5) )–(8(8) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{R}_h}{\mathrm{d}t}}={\alpha }_h{I}_h-{\mu }_h{R}_h,\end{array}$(8) ), we have ${\displaystyle \frac{\mathrm{d}{N}_h}{\mathrm{d}t}}={\mu }_h{N}_h^{\ast }-{\mu }_h{N}_h,$ where $N_h=S_h+E_h+I_h+R_h$. Solving the differential equation results in $N_h\left(t\right)=N_h^{\ast }+\left(N_h\right(0)-N_h^{\ast })e^{-{\mu }_{h}t}$ Thus $\underset{t\to \mathrm{\infty }}{lim sup}{N}_h\left(t\right)\le N_h^{\ast }$. Hence, $S_h\left(t\right)$, $E_h\left(t\right)$ $R_h\left(t\right)$, and $I_h\left(t\right)$ are bounded.

   From Equations (9(9) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{S}_s}{\mathrm{d}t}}={\mu }_s{N}_s^{\ast }+\rho V_s-{\mu }_s{S}_s-{\displaystyle \frac{{\text{β}}_{es}B}{{\chi }_s+B}}{S}_s-\nu S_s\end{array}$(9) )–(12(12) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{I}_s}{\mathrm{d}t}}={\gamma }_s{E}_s-{\mu }_s{I}_s,\end{array}$(12) ), we have ${\displaystyle \frac{\mathrm{d}{N}_s}{\mathrm{d}t}}={\mu }_s{N}_s^{\ast }-{\mu }_s{N}_s,$ where $N_s=S_s+E_s+I_s+V_s$. In similar manner, $\underset{t\to \mathrm{\infty }}{lim sup}{N}_s\left(t\right)\le N_s^{\ast }$. Hence, $S_s\left(t\right)$, $E_s\left(t\right)$, $I_s\left(t\right)$ and $V_s\left(t\right)$ are bounded. Finally, we consider Equation (4(4) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}B}{\mathrm{d}t}}=\delta I_d-({\mu }_e+\mu )B,\end{array}$(4) ) of the model. Since $I_d\left(t\right)\le N_d\left(t\right)\le N_d^{\ast }$, we have $\begin{array}{rl}B\left(t\right)& \le e^{-{\int }_0^t({\mu }_e+\mu )\mathrm{d}\tau }\left({\int }_0^t\delta N_d^{\ast }{e}^{{\int }_0^{\tau }({\mu }_e+\mu )\mathrm{d}s}\mathrm{d}\tau \right)\\ & ={\displaystyle \frac{\delta N_d^{\ast }}{{\mu }_e+\mu }}-{\displaystyle \frac{\delta N_d^{\ast }}{{\mu }_e+\mu }}{e}^{-({\mu }_e+\mu )t}\\ & \le {\displaystyle \frac{\delta N_d^{\ast }}{{\mu }_e+\mu }}\end{array}$

Appendix 2. Calculation of the control reproduction number

According to the concepts of the next generation matrix and reproduction number presented in [8] and [27], we define $\mathcal{F}=\left(\begin{array}{c}{\text{β}}_{sd}{I}_s{S}_d\\ 0\\ 0\\ {\displaystyle \frac{{\text{β}}_{eh}B}{{\chi }_h+B}}{S}_h\\ 0\\ {\displaystyle \frac{{\text{β}}_{es}B}{{\chi }_s+B}}{S}_s\\ 0\end{array}\right)\mathrm{a}\mathrm{n}\mathrm{d}\mathcal{V}=\left(\begin{array}{c}\widetilde{{\gamma }_d}{E}_d\\ {\mu }_d{I}_d-{\gamma }_d{E}_d\\ ({\mu }_e+\mu )B-\delta I_d\\ \widetilde{{\gamma }_h}{E}_h\\ \widetilde{q_h}{I}_h-{\gamma }_h{E}_h\\ \widetilde{{\gamma }_s}{E}_s\\ {\mu }_s{I}_s-{\gamma }_s{E}_s\end{array}\right),$ 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 (A1a) $F=\left(\begin{array}{ccccccc}0& 0& 0& 0& 0& 0& {\text{β}}_{sd}{N}_d^{\ast }\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& {\displaystyle \frac{{\text{β}}_{eh}{N}_h^{\ast }}{{\chi }_h}}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& {\displaystyle \frac{{\text{β}}_{es}(\rho +{\mu }_s)N_s^{\ast }}{(\rho +{\mu }_s+\nu ){\chi }_s}}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\end{array}\right),$(A1a) and V is a matrix of transition rates given by (A1b) $V=\left(\begin{array}{ccccccc}\widetilde{{\gamma }_d}& 0& 0& 0& 0& 0& 0\\ -{\gamma }_d& {\mu }_d& 0& 0& 0& 0& 0\\ 0& -\delta & {\mu }_e+\mu & 0& 0& 0& 0\\ 0& 0& 0& \widetilde{{\gamma }_h}& 0& 0& 0\\ 0& 0& 0& -{\gamma }_h& \widetilde{q_h}& 0& 0\\ 0& 0& 0& 0& 0& \widetilde{{\gamma }_s}& 0\\ 0& 0& 0& 0& 0& -{\gamma }_s& {\mu }_s\end{array}\right).$(A1b) Thus $V^{-1}=\left(\begin{array}{ccccccc}{\displaystyle \frac{1}{\widetilde{{\gamma }_d}}}& 0& 0& 0& 0& 0& 0\\ {\displaystyle \frac{{\gamma }_d}{{\mu }_d\widetilde{{\gamma }_d}}}& {\displaystyle \frac{1}{{\mu }_d}}& 0& 0& 0& 0& 0\\ {\displaystyle \frac{\delta {\gamma }_d}{({\mu }_e+\mu ){\mu }_d\widetilde{{\gamma }_d}}}& {\displaystyle \frac{\delta }{({\mu }_e+\mu ){\mu }_d}}& {\displaystyle \frac{1}{{\mu }_e+\mu }}& 0& 0& 0& 0\\ 0& 0& 0& {\displaystyle \frac{1}{\widetilde{{\gamma }_h}}}& 0& 0& 0\\ 0& 0& 0& {\displaystyle \frac{{\gamma }_h}{\widetilde{{\gamma }_h}\widetilde{q_h}}}& {\displaystyle \frac{1}{\widetilde{q_h}}}& 0& 0\\ 0& 0& 0& 0& 0& {\displaystyle \frac{1}{\widetilde{{\gamma }_s}}}& 0\\ 0& 0& 0& 0& 0& {\displaystyle \frac{{\gamma }_s}{{\mu }_s\widetilde{{\gamma }_s}}}& {\displaystyle \frac{1}{{\mu }_s}}\end{array}\right),$ and the next generation matrix is $F{V}^{-1}=\left(\begin{array}{ccccccc}0& 0& 0& 0& 0& {\displaystyle \frac{{\gamma }_s{\text{β}}_{sd}{N}_d^{\ast }}{{\mu }_s\widetilde{{\gamma }_s}}}& {\displaystyle \frac{{\text{β}}_{sd}{N}_d^{\ast }}{{\mu }_s}}\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ {\displaystyle \frac{\delta r_d{\text{β}}_{eh}{N}_h^{\ast }}{{\chi }_h{\mu }_d\widetilde{{\gamma }_d}({\mu }_e+\mu )}}& {\displaystyle \frac{\delta {\text{β}}_{eh}{N}_h^{\ast }}{{\chi }_h{\mu }_d({\mu }_e+\mu )}}& {\displaystyle \frac{{\text{β}}_{eh}{N}_h^{\ast }}{{\chi }_h({\mu }_e+\mu )}}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ {\displaystyle \frac{\varphi \delta r_d{\text{β}}_{es}{N}_s^{\ast }}{{\chi }_s{\mu }_d\widetilde{{\gamma }_d}({\mu }_e+\mu )}}& {\displaystyle \frac{\varphi \delta {\text{β}}_{es}{N}_s^{\ast }}{{\chi }_s{\mu }_d({\mu }_e+\mu )}}& {\displaystyle \frac{\varphi {\text{β}}_{es}{N}_s^{\ast }}{{\chi }_s({\mu }_e+\mu )}}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\end{array}\right),$ where $\varphi =\frac{\rho +{\mu }_s}{\rho +{\mu }_s+\nu }$. Thus the control reproduction number is the spectral radius given by ${\mathcal{R}}_c=\rho \left(F{V}^{-1}\right)=\sqrt{{\displaystyle \frac{\rho +{\mu }_s}{\rho +{\mu }_s+\nu }}{\displaystyle \frac{{\gamma }_s{\text{β}}_{es}{N}_s^{\ast }}{{\chi }_s({\mu }_e+\mu )}}{\displaystyle \frac{{\text{β}}_{sd}}{{\mu }_s\widetilde{{\gamma }_s}}}{\displaystyle \frac{\delta {\gamma }_d{N}_d^{\ast }}{{\mu }_d\widetilde{{\gamma }_d}}}}$

Appendix 3. Proof of Theorem 4.2

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: $\begin{array}{rl}L& =w_1(S_d-S_d^{\ast }{)}^2+w_2(E_d-E_d^{\ast }{)}^2+w_3(I_d-I_d^{\ast }{)}^2+w_4(B-B^{\ast }{)}^2\\ & +{\displaystyle \frac{1}{2}}(S_h-S_h^{\ast }+E_h-E_h^{\ast }+I_h-I_h^{\ast }+R_h-R_h^{\ast }{)}^2\\ & +{\displaystyle \frac{1}{2}}(S_s-S_s^{\ast }+E_s-E_s^{\ast }+I_s-I_s^{\ast }+V_s-V_s^{\ast }{)}^2\end{array}$ Thus $\begin{array}{rl}{\displaystyle \frac{\mathrm{d}L}{\mathrm{d}t}}& =2w_1(S_d-S_d^{\ast })\left({\displaystyle \frac{\mathrm{d}{S}_d}{\mathrm{d}t}}\right)+2w_2(E_d-E_d^{\ast })\left({\displaystyle \frac{d{E}_d}{dt}}\right)\\ & +2w_3(I_d-I_d^{\ast })\left({\displaystyle \frac{\mathrm{d}{I}_d}{\mathrm{d}t}}\right)+2w_4(B-B^{\ast })\left({\displaystyle \frac{\mathrm{d}B}{\mathrm{d}t}}\right)\\ & +(S_s-S_s^{\ast }+E_s-E_s^{\ast }+I_s-I_s^{\ast }+V_s-V_s^{\ast })\left({\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\right(S_s+E_s+I_s+V_s\left)\right)\\ & +(S_h-S_h^{\ast }+E_h-E_h^{\ast }+I_h-I_h^{\ast }+R_h-R_h^{\ast })\left({\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\right(S_h+E_h+I_h+R_h\left)\right)\end{array}$ The time derivative of $L\left(t\right)$ along the solutions of model Equations (1(1) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{S}_d}{\mathrm{d}t}}={\mu }_d{N}_d^{\ast }-{\text{β}}_{sd}{I}_s{S}_d-{\mu }_d{S}_d,\end{array}$(1) )–(12(12) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{I}_s}{\mathrm{d}t}}={\gamma }_s{E}_s-{\mu }_s{I}_s,\end{array}$(12) ), we obtain $\begin{array}{rl}{\displaystyle \frac{dL}{dt}}=& Y(WA+A^T{W}^T)Y^T-{\mu }_s(S_s-S_s^{\ast }+E_s-E_s^{\ast }+I_s-I_s^{\ast }+V_s-V_s^{\ast }{)}^2\\ & -{\mu }_h(S_h-S_h^{\ast }+E_h-E_h^{\ast }+I_h-I_h^{\ast }+R_h-R_h^{\ast }{)}^2,\end{array}$ where $Y=(S_d-S_d^{\ast },E_d-E_d^{\ast },I_d-I_d^{\ast },B-B^{\ast })$, $W=diag(w_1,w_2,w_3,w_4)$ and $A=\left(\begin{array}{cccc}-{\text{β}}_{sd}{I}_s-{\mu }_d& 0& 0& 0\\ {\text{β}}_{sd}{I}_s& -\widetilde{{\gamma }_d}& 0& 0\\ 0& {\gamma }_d& -{\mu }_d& 0\\ 0& 0& \delta & -{\mu }_e-\mu \end{array}\right).$ Obviously, the second and third terms of ${\displaystyle \frac{dL}{dt}}$ are negative. To show the global stability of endemic equilibrium $X_E$, it suffice to show that A is Volterra–Lyapunov stable in $\mathbb{D}/\left\{X_E\right\}$. For this purpose, we show that matrix $\tilde{A}$ is Volterra–Lyapunov stable, and matrix $U=\widetilde{A^{-1}}$ is Volterra–Lyapunov stable (or $-\widetilde{A^{-1}}$ is diagonally stable) [17,31].

Condition 1:

To show $\tilde{A}$ is Volterra–Lyapunov stable, we consider $\begin{array}{r}D=-\tilde{A}=\left(\begin{array}{ccc}{\text{β}}_{sd}{I}_s+{\mu }_d& 0& 0\\ -{\text{β}}_{sd}{I}_s& \widetilde{{\gamma }_d}& 0\\ 0& -{\gamma }_d& {\mu }_d\end{array}\right).\end{array}$ The matrix $\tilde{D}$ is obtained by deleting its third row and third column, and results in $-\tilde{D}=\left(\begin{array}{cc}-{\text{β}}_{sd}{I}_s-{\mu }_d& 0\\ {\text{β}}_{sd}{I}_s& -\widetilde{{\gamma }_d}\end{array}\right).$ Clearly $a_{11}=-{\text{β}}_{sd}{I}_s-{\mu }_d<0$, $a_{22}=-\widetilde{{\gamma }_d}$ and $det(-\tilde{D})>0$. Based on Lemma 2.4 of [17], $-\tilde{D}$ is Volterra–Lyapunov stable. Moreover, we obtain $\begin{array}{r}-\widetilde{D^{-1}}=\left(\begin{array}{cc}-{\displaystyle \frac{1}{{\text{β}}_{sd}{I}_s+{\mu }_d}}& 0\\ -{\displaystyle \frac{{\text{β}}_{sd}{I}_s}{\widetilde{{\gamma }_d}({\mu }_d+{\text{β}}_{sd}{I}_s)}}& -{\displaystyle \frac{1}{\widetilde{{\gamma }_d}}}\end{array}\right).\end{array}$ Clearly, the diagonal elements $-\frac{1}{{\text{β}}_{sd}{I}_s+{\mu }_d}<0$, $-\frac{1}{\widetilde{{\gamma }_d}}<0$ and $det(-\widetilde{D^{-1}})>0$. Thus, based on Lemma 2.4 of [17], $-\widetilde{D^{-1}}$ is Volterra–Lyapunov stable. Therefore, $D=-\tilde{A}$ is diagonally stable, and hence $\tilde{A}$ is Volterra–Lyapunov stable.

Condition 2:

To show $\widetilde{A^{-1}}$ is Volterra–Lyapunov stable, we consider $E=-\widetilde{A^{-1}}=\left(\begin{array}{ccc}{\displaystyle \frac{1}{{\text{β}}_{sd}{I}_s+{\mu }_d}}& 0& 0\\ {\displaystyle \frac{{\text{β}}_{sd}{I}_s}{\widetilde{{\gamma }_d}({\mu }_d+{\text{β}}_{sd}{I}_s)}}& {\displaystyle \frac{1}{\widetilde{{\gamma }_d}}}& 0\\ {\displaystyle \frac{{\text{β}}_{sd}{I}_s{\gamma }_d}{\widetilde{{\gamma }_d}{\mu }_d({\mu }_d+{\text{β}}_{sd}{I}_s)}}& {\displaystyle \frac{{\gamma }_d}{\widetilde{{\gamma }_d}{\mu }_d}}& {\displaystyle \frac{1}{{\mu }_d}}\end{array}\right).$ This results in $-\tilde{E}=\left(\begin{array}{cc}-{\displaystyle \frac{1}{{\text{β}}_{sd}{I}_s+{\mu }_d}}& 0\\ -{\displaystyle \frac{{\text{β}}_{sd}{I}_s}{\widetilde{{\gamma }_d}({\mu }_d+{\text{β}}_{sd}{I}_s)}}& -{\displaystyle \frac{1}{\widetilde{{\gamma }_d}}}\end{array}\right).$ Using Lemma 2.4 of [17], it is easy to observe that $-\tilde{E}$ is Volterra–Lyapunov stable. Moreover, we have $-\widetilde{E^{-1}}=\left(\begin{array}{cc}-{\text{β}}_{sd}{I}_s-{\mu }_d& 0\\ {\text{β}}_{sd}{I}_s& -\widetilde{{\gamma }_d}\end{array}\right).$ Based on Lemma 2.4 of [17], we can also observe that $-\widetilde{E^{-1}}$ is also Volterra–Lyapunov stable. Therefore, $\widetilde{A^{-1}}$ is Volterra–Lyapunov stable.

Based on Lemma 2.8 of [17], there exist a diagonal matrix $W=diag\{w_1,w_2,w_3,w_4\}$ 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, $\frac{\mathrm{d}L}{\mathrm{d}t}<0$, and by LaSalle's invariance principle [16], $X_E$ is globally asymptotically stable in the interior of $\mathbb{D}$, and is unique provided that ${\mathcal{R}}_c>1$.