Dynamics of an epidemic model with imperfect vaccinations on complex networks

In recent times, we have been facing the threat of some infectious diseases such as influenza, HIV/AIDS, dengue, tuberculosis, Zika virus disease, Ebola virus disease, and the recent and rapidly emerging outbreak disease—coronavirus disease 2019 (COVID-19). These infectious diseases usually have a severe impact on human life and the economic system. Therefore, it is urgent to formulate a good epidemic prevention strategy. To prevent and control an infectious disease, we must understand the transmission dynamics of the disease, and try to predict the future evolution of the disease. Mathematical modeling is a powerful tool for analyzing the spread of infectious diseases and this approach has been extensively studied over a period of time [3, 24]. Among the various types of epidemic models, in this study, we are interested in the compartmental models. In a compartmental model, the population is assigned to different compartments according to epidemiological states (e.g., susceptible, infected, exposed, etc), and in a compartment, individuals have the same epidemiological state. Usually, we can formulate a system of differential equations to investigate how the compartments change over time.

Vaccination is recognized as an effective means of disease control and prevention. Therefore, vaccination should play an important role in the transmission dynamics of an infectious disease. There are several compartmental epidemic models that include such a factor to investigate the role of this control strategy [1, 2, 4, 5, 9, 10, 22, 27, 31, 32, 35, 37]. Although vaccination is a useful tool for reducing or eradicating the incidence of an infection, it is well known that not all vaccinations are completely effective. Therefore, it is important to study the effect of imperfect vaccinations on the transmission dynamics of infectious diseases. In general, there are three different types of imperfect vaccinations that are commonly considered in the literature. We use the terminology in [12]: (i) leaky vaccine is used to describe a vaccine that only reduces the probability of infection upon exposure; however, it does not eliminate it; (ii) all-or-nothing vaccine describes a vaccine that has no effect on some individuals but confers complete protection in others; (iii) waning vaccine is the one for which the protection conferred wanes over time. Different epidemic models including all three aforementioned imperfect vaccinations were investigated in [1, 2, 9, 22, 31]. A model that incorporates a nonlinear incidence rate and leaky vaccine was studied in [4]. The epidemic models including all-or-nothing and waning vaccines with a nonlinear incidence rate and a saturated treatment were considered in [5, 10] respectively. A model that incorporates the waning vaccine and general nonlinear incidence was analyzed in [37]. Epidemic models with imperfect vaccination under the impact of opinion dynamics were investigated in [29, 30]. A model that includes the element of vaccine failure and three types of infection awareness and immunization awareness was studied in [15].

Thus far, existing studies have mainly been based on the assumption that individuals in the population are homogeneously mixed. However, it is commonly accepted that contacts among individuals in a population are heterogeneous, and can have vital effects on the transmission of infectious diseases. To describe the effect of contact heterogeneity, the concept of complex networks is included in epidemic modeling. Because of the seminal study by Pastor-Satorras and Vespignani [28], several network-based epidemic models have been proposed and studied [6, 8, 14, 16, 18–20, 25, 33]. It should be noted that there have been relatively fewer studies on network epidemic models with imperfect vaccinations. In [7], two epidemic models incorporating variable population size, degree-related vaccination rate, and quarantine were introduced. Global stability and the conditions for the existence of multiple endemic equilibria were investigated. In [21], the authors studied an epidemic model with imperfect vaccination on scale-free networks. Their model included two degree-related transmission rates and a degree-related vaccination rate. Two threshold parameters were derived, and the existence conditions of multiple endemic equilibria were discussed. Furthermore, the effects of different immunization schemes and an optimal control strategy of vaccination were also presented. Generally, the results of network epidemic models with imperfect vaccinations are still rare. To gain more insight into the effects of imperfect vaccinations on the transmission of infectious diseases, in this study, we investigate a network-based epidemic model with three types of the aforementioned imperfect vaccinations. We first obtain a threshold value ${\mathcal{R}}_{v}$ that determines the stability of the disease-free equilibrium. Subsequently, we consider the existence of the endemic equilibrium by a bifurcation analysis. Furthermore, we study the persistence of the disease. Numerical simulations are conducted to validate our theoretical findings. The novelty and contributions of this work are as follows: (i) we not only performed the stability analysis of the disease-free equilibrium, but also studied the bifurcation behavior of the endemic equilibria; (ii) we obtained a necessary and sufficient condition that ensures the occurrence of a backward bifurcation and this condition is independent of network structure; (iii) from our analysis, we observed that the leaky vaccine plays the most important role in the dynamics of epidemic spreading. The leaky vaccine can not only cause the backward bifurcation, but also result in the largest region of backward bifurcation in comparison with the other two imperfect vaccines.

The remainder of this paper is organized as follows. In section 2, we describe the proposed model, and study the stability of the disease-free equilibrium. In section 3, we derive a condition that determines the direction of bifurcating the endemic equilibrium. In section 4, we study the persistence of the disease. In section 5, numerical experiments are detailed to illustrate the theoretical results. Finally, the conclusions are presented in section 6.

In this section, we present the network-based epidemic model with imperfect vaccinations inspired by the study by Magpantay et al [22]. We assume a population of individuals whose connections to one another form a network. Specifically, in this network, each individual is represented by a node, and the connections among the individuals are formed by edges. The degree of a node is the number of edges between this node and other ones. We classify all the nodes of this network into n groups, such that nodes in the same group have the same degree. This indicates that each node in the kth group has a degree k for k = 1, 2, ..., n, where n is the maximal degree of the network. Furthermore, in each group, we divide the nodes into three categories: susceptible, infected, and vaccinated. Let S_k (t), I_k (t), and V_k (t) represent the densities of susceptible, infected, and vaccinated nodes with node degree k at time t, respectively. The network-based epidemic model is governed by the following system of ordinary differential equations:

$$\begin{equation}\begin{cases}_{k}^{\prime }\left(t\right)=\left(1-\left(1-{\varepsilon }_{A}\right)p\right)\mu -\beta k{S}_{k}\left(t\right){\Theta}\left(t\right)+\gamma {I}_{k}\left(t\right)+\alpha {V}_{k}\left(t\right)-\mu {S}_{k}\left(t\right),\quad \hfill \\ {I}_{k}^{\prime }\left(t\right)=\beta k{S}_{k}\left(t\right){\Theta}\left(t\right)+{\varepsilon }_{L}\beta k{V}_{k}\left(t\right){\Theta}\left(t\right)-\left(\gamma +\mu \right){I}_{k}\left(t\right),\quad \hfill \\ {V}_{k}^{\prime }\left(t\right)=\left(1-{\varepsilon }_{A}\right)p\mu -{\varepsilon }_{L}\beta k{V}_{k}\left(t\right){\Theta}\left(t\right)-\left(\alpha +\mu \right){V}_{k}\left(t\right),\quad k=1,2,\dots ,n,\quad \hfill \end{cases}\end{equation} \tag{ 1 }$$

where the descriptions of all the parameters are given in table 1. In this model, we assume that the model population has equal birth and natural death rates of μ. The term (1 − (1 − ɛ_A )p)μ describes that a fraction p of the newborns are assumed to be vaccinated, and of this fraction, a proportion ɛ_A exhibit vaccine failure (i.e., all-or-nothing vaccine). The term βkS_k (t)Θ(t) represents the incidence, and the factor Θ(t) is the probability of a given edge to be connected to an infected node. Specifically, the probability that an edge is connected to a node with degree h is proportional to hP(h), where P(h) > 0 represents the probability of a node with degree h, and hence, ${\sum }_{h=1}^{n}P\left(h\right)=1$. Thus, the probability that a randomly chosen edge is connected to an infected node is given by

$$\begin{equation*}{\Theta}\left(t\right)=\frac{{\sum }_{h=1}^{n}hP\left(h\right){I}_{h}\left(t\right)}{{\sum }_{h=1}^{n}hP\left(h\right)}=\frac{1}{\langle k\rangle }\sum _{h=1}^{n}hP\left(h\right){I}_{h}\left(t\right),\end{equation*}$$

where $\langle k\rangle ={\sum }_{h=1}^{n}hP\left(h\right)$ denotes the average degree of the network. The parameter ɛ_L is introduced to model the phenomenon that vaccination may reduce the probability of infection upon exposure; however, it does not completely eliminate susceptibility to infection (i.e., leaky vaccine). In this study, we also assume the waning rate of vaccinated individuals by α = μɛ_W /(1 − ɛ_W ) (cf [22]) such that the parameter

$$\begin{equation*}{\varepsilon }_{W}=\frac{\alpha }{\alpha +\mu }=\frac{\mathrm{t}\mathrm{h}\mathrm{e}\enspace \mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}\enspace \mathrm{a}\mathrm{t}\enspace \mathrm{w}\mathrm{h}\mathrm{i}\mathrm{c}\mathrm{h}\enspace \mathrm{v}\mathrm{a}\mathrm{c}\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{e}\enspace \mathrm{p}\mathrm{r}\mathrm{o}\mathrm{t}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\enspace \mathrm{w}\mathrm{a}\mathrm{n}\mathrm{e}\mathrm{s}}{\mathrm{t}\mathrm{h}\mathrm{e}\enspace \mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}\enspace \mathrm{a}\mathrm{t}\enspace \mathrm{w}\mathrm{h}\mathrm{i}\mathrm{c}\mathrm{h}\enspace \mathrm{i}\mathrm{n}\mathrm{d}\mathrm{i}\mathrm{v}\mathrm{i}\mathrm{d}\mathrm{u}\mathrm{a}\mathrm{l}\mathrm{s}\enspace \mathrm{l}\mathrm{e}\mathrm{a}\mathrm{v}\mathrm{e}\enspace \mathrm{t}\mathrm{h}\mathrm{e}\enspace \mathrm{v}\mathrm{a}\mathrm{c}\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{d}\enspace \mathrm{c}\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{s}\enspace \mathrm{i}\mathrm{f}\enspace {\varepsilon }_{\text{L}}=0}.\end{equation*}$$

The parameter γ stands for the recovery rate of infected individuals. Let N_k (t) = S_k (t) + I_k (t) + V_k (t) for all k = 1, 2, ..., n. It follows from system (1) that N_k '(t) = μ − μN_k (t). Thus, we assume that the initial values satisfy the normalization condition, namely, N_k (0) = S_k (0) + I_k (0) + V_k (0) = 1 for all k, then, we have N_k (t) = 1 for all t ⩾ 0 and k = 1, 2, ..., n. Therefore, we consider the dynamics of the solutions of system (1) in the following biologically meaningful region:

$$\begin{equation*}{\Omega}{:=}\left\{\left({S}_{1},{I}_{1},{V}_{1},\dots ,{S}_{n},{I}_{n},{V}_{n}\right):0{\leqslant}{S}_{k},{I}_{k},{V}_{k}{\leqslant}1,{S}_{k}+{I}_{k}+{V}_{k}=1,\quad 1{\leqslant}k{\leqslant}n\right\}\end{equation*}$$

and the initial values are given in Ω. It is easy to verify that Ω is a positively invariant set for system (1), that is, every trajectory of (1) starting in Ω remains in Ω for all t ⩾ 0.

To find the equilibrium solutions for system (1), we can observe that there always exists a disease-free equilibrium ${\mathcal{E}}_{0}=\left({E}_{0},{E}_{0},\dots ,{E}_{0}\right)\in {\mathbb{R}}^{3n}$ with

$$\begin{equation}{E}_{0}=\left({S}_{0},{I}_{0},{V}_{0}\right){:=}\left(\frac{\alpha +\mu -\left(1-{\varepsilon }_{A}\right)p\mu }{\alpha +\mu },0,\frac{\left(1-{\varepsilon }_{A}\right)p\mu }{\alpha +\mu }\right).\end{equation} \tag{ 2 }$$

Next, we investigate the stability of the disease-free equilibrium ${\mathcal{E}}_{0}$. We first establish the local stability of ${\mathcal{E}}_{0}$, and subsequently, we consider the global stability of ${\mathcal{E}}_{0}$. The proofs of the stability results of ${\mathcal{E}}_{0}$ are given in appendix A.

Theorem 2.1. Define the threshold value

$$\begin{equation*}{\mathcal{R}}_{v}{:=}\frac{\alpha +\mu -\left(1-{\varepsilon }_{A}\right)p\mu +{\varepsilon }_{L}\left(1-{\varepsilon }_{A}\right)p\mu }{\alpha +\mu }\frac{\beta \langle {k}^{2}\rangle }{\left(\gamma +\mu \right)\langle k\rangle },\end{equation*}$$

where $\langle {k}^{2}\rangle ={\sum }_{h=1}^{n}{h}^{2}P\left(h\right)$. The disease-free equilibrium ${\mathcal{E}}_{0}$ is locally asymptotically stable if ${\mathcal{R}}_{v}{< }1$, and it is unstable if ${\mathcal{R}}_{v}{ >}1$.

Remark 1. Given that α = μɛ_W /(1 − ɛ_W ), one can derive that

$$\begin{align*}\hfill {\mathcal{R}}_{v}& =\left(1-\left(1-{\varepsilon }_{L}\right)\left(1-{\varepsilon }_{A}\right)\left(1-{\varepsilon }_{W}\right)p\right)\frac{\beta \langle {k}^{2}\rangle }{\left(\gamma +\mu \right)\langle k\rangle }\hfill \\ \hfill & =\left(1-\left(1-{\varepsilon }_{L}\right)\left(1-{\varepsilon }_{A}\right)\left(1-{\varepsilon }_{W}\right)p\right){\mathcal{R}}_{0},\hfill \end{align*}$$

where ${\mathcal{R}}_{0}{:=}\beta \langle {k}^{2}\rangle /\left(\left(\gamma +\mu \right)\langle k\rangle \right)$ is the threshold value in the absence of vaccination (p = 0). It can be seen that incorporating an effective vaccine can reduce the threshold value. Therefore, vaccination can prevent the outbreak of diseases.

From remark 1, we can observe that ${\mathcal{R}}_{v}{\leqslant}{\mathcal{R}}_{0}$. The following result shows that ${\mathcal{E}}_{0}$ is globally asymptotically stable when ${\mathcal{R}}_{0}{< }1$.

Theorem 2.2. If ${\mathcal{R}}_{0}{< }1$, then ${\mathcal{E}}_{0}$ is globally asymptotically stable.

From theorem 2.2, we only obtain the global stability of ${\mathcal{E}}_{0}$ when ${\mathcal{R}}_{0}{< }1$ instead of ${\mathcal{R}}_{v}{< }1$. This indicates that the condition ${\mathcal{R}}_{v}{< }1$ is not sufficient to eliminate the disease. Therefore, it is important to discuss whether the disease can persist even if ${\mathcal{R}}_{v}{< }1$. In this section, we characterize whether there exists a backward (or subcritical) bifurcation at ${\mathcal{R}}_{v}=1$. When such a bifurcation occurs, there exists a range of the threshold value ${\mathcal{R}}_{v}$ (i.e., ${\mathcal{R}}_{v}^{{\ast}}{< }{\mathcal{R}}_{v}{< }1$) that we refer to as the region of metastability. In this region, there are at least two endemic equilibria and typically, at least one of which is stable. This implies that ${\mathcal{E}}_{0}$ is not globally stable when ${\mathcal{R}}_{v}{< }1$, and the disease can persist even if ${\mathcal{R}}_{v}{< }1$. There are some mechanisms that have been identified to induce such a bifurcation. More precisely, epidemic models can exhibit backward bifurcation owing to mechanisms such as nonlinear treatment rates [10, 17, 19], vaccine failure [1, 2, 4, 9, 31], and host disease-induced mortality [11]. Besides, a model that considers a coupling of an epidemic spreading and an opinion dynamics has been shown to be able to undergo such a bifurcation [30].

Here, we study the bifurcation phenomenon at ${\mathcal{R}}_{v}=1$. More specifically, we will derive a condition that determines the direction of the bifurcation at ${\mathcal{R}}_{v}=1$. Therefore, let us consider the endemic equilibria of system (1). At an endemic equilibrium, we have S_k , I_k , V_k > 0 for all k = 1, 2, ..., n. Furthermore, it should satisfy the following equations:

$$\begin{equation}\begin{cases}\left(1-\left(1-{\varepsilon }_{A}\right)p\right)\mu -\beta k{S}_{k}{\Theta}+\gamma {I}_{k}+\alpha {V}_{k}-\mu {S}_{k}=0,\quad \hfill \\ \beta k{S}_{k}{\Theta}+{\varepsilon }_{L}\beta k{V}_{k}{\Theta}-\left(\gamma +\mu \right){I}_{k}=0,\quad \hfill \\ \left(1-{\varepsilon }_{A}\right)p\mu -{\varepsilon }_{L}\beta k{V}_{k}{\Theta}-\left(\alpha +\mu \right){V}_{k}=0.\quad \hfill \end{cases}\end{equation} \tag{ 3 }$$

It can be deduced from the third equation of system (3) that

$$\begin{equation}{V}_{k}=\frac{\left(1-{\varepsilon }_{A}\right)p\mu }{{\varepsilon }_{L}\beta k{\Theta}+\alpha +\mu }.\end{equation} \tag{ 4 }$$

Combining S_k + I_k + V_k = 1 with the first equation of system (3), we can derive

$$\begin{equation}{S}_{k}=\frac{\left(1-\left(1-{\varepsilon }_{A}\right)p\right)\mu +\alpha +\left(\gamma -\alpha \right){I}_{k}}{\beta k{\Theta}+\alpha +\mu }.\end{equation} \tag{ 5 }$$

Substituting the relations (4) and (5) into the second equation of system (3), we can obtain

$$\begin{equation}{I}_{k}=\frac{\beta k{\Theta}\left(\left(1-\left(1-{\varepsilon }_{A}\right)p\right)\mu +\alpha \right){{\Delta}}_{1}\left(k,{\Theta}\right)+{\varepsilon }_{L}\beta k{\Theta}\left(1-{\varepsilon }_{A}\right)p\mu {{\Delta}}_{2}\left(k,{\Theta}\right)}{\left(\beta k{\Theta}+\gamma +\mu \right)\left(\alpha +\mu \right){{\Delta}}_{1}\left(k,{\Theta}\right)},\end{equation} \tag{ 6 }$$

where Δ₁(k, Θ) = ɛ_L βkΘ + α + μ and Δ₂(k, Θ) = βkΘ + α + μ. Substituting relation (6) into the expression of Θ, the following equation can be derived:

$$\begin{equation}{\Theta}=\frac{1}{\langle k\rangle }\sum _{h=1}^{n}\frac{\beta {h}^{2}P\left(h\right){\Theta}\left(\left(1-\left(1-{\varepsilon }_{A}\right)p\right)\mu +\alpha \right){{\Delta}}_{1}\left(h,{\Theta}\right)+{\varepsilon }_{L}\beta {h}^{2}P\left(h\right){\Theta}\left(1-{\varepsilon }_{A}\right)p\mu {{\Delta}}_{2}\left(h,{\Theta}\right)}{\left(\beta h{\Theta}+\gamma +\mu \right)\left(\alpha +\mu \right){{\Delta}}_{1}\left(h,{\Theta}\right)}.\end{equation} \tag{ 7 }$$

Because we are looking for the endemic equilibria, we have I_k > 0 for all k. This implies that Θ > 0. Therefore, the endemic equilibria should obey the following equation of Θ from equation (7):

$$\begin{equation}\frac{1}{\langle k\rangle }\sum _{h=1}^{n}\frac{\beta {h}^{2}P\left(h\right)\left(\left(1-\left(1-{\varepsilon }_{A}\right)p\right)\mu +\alpha \right){{\Delta}}_{1}\left(h,{\Theta}\right)+{\varepsilon }_{L}\beta {h}^{2}P\left(h\right)\left(1-{\varepsilon }_{A}\right)p\mu {{\Delta}}_{2}\left(h,{\Theta}\right)}{\left(\beta h{\Theta}+\gamma +\mu \right)\left(\alpha +\mu \right){{\Delta}}_{1}\left(h,{\Theta}\right)}=1.\end{equation} \tag{ 8 }$$

Once we get Θ from equation (8), S_k , I_k , and V_k can be determined according to the relations (4)–(6). Consequently, the endemic equilibria can be found. We will study the relationship between ${\mathcal{R}}_{v}$ and Θ. Therefore, we try to express equation (8) in terms of ${\mathcal{R}}_{v}$ and Θ. From the definition of ${\mathcal{R}}_{v}$, we set

$$\begin{equation*}\beta =\frac{\alpha +\mu }{\alpha +\mu -\left(1-{\varepsilon }_{A}\right)p\mu +{\varepsilon }_{L}\left(1-{\varepsilon }_{A}\right)p\mu }\frac{\left(\gamma +\mu \right)\langle k\rangle }{\langle {k}^{2}\rangle }{\mathcal{R}}_{v}{:=}\tau {\mathcal{R}}_{v}.\end{equation*}$$

Therefore, equation (8) becomes

$$\begin{equation}\frac{1}{\langle k\rangle }\sum _{h=1}^{n}\frac{\tau {\mathcal{R}}_{v}{h}^{2}P\left(h\right)\left(\left(1-\left(1-{\varepsilon }_{A}\right)p\right)\mu +\alpha \right){{\Delta}}_{1}\left(h,{\Theta}\right)+{\varepsilon }_{L}\tau {\mathcal{R}}_{v}{h}^{2}P\left(h\right)\left(1-{\varepsilon }_{A}\right)p\mu {{\Delta}}_{2}\left(h,{\Theta}\right)}{\left(\tau {\mathcal{R}}_{v}h{\Theta}+\gamma +\mu \right)\left(\alpha +\mu \right){{\Delta}}_{1}\left(h,{\Theta}\right)}=1.\end{equation} \tag{ 9 }$$

Equation (9) can be used to define an endemic equilibrium curve in the first quadrant of the ${\mathcal{R}}_{v}{\Theta}$-plane. If we treat Θ as a function of ${\mathcal{R}}_{v}$, the direction of the bifurcation is determined by the slope of the curve at $\left({\mathcal{R}}_{v},{\Theta}\right)=\left(1,0\right)$. This indicates that if the derivative is negative at the critical value $\left({\mathcal{R}}_{v},{\Theta}\right)=\left(1,0\right)$, backward bifurcation occurs, and if the derivative at $\left({\mathcal{R}}_{v},{\Theta}\right)=\left(1,0\right)$ is positive, forward bifurcation takes place (figure 1). Therefore, we can obtain a necessary and sufficient criterion that decides the direction of the bifurcation by computing the derivative $\partial {\Theta}/\partial {\mathcal{R}}_{v}$ at $\left({\mathcal{R}}_{v},{\Theta}\right)=\left(1,0\right)$. To obtain this condition, we differentiate equation (9) implicitly with respect to R_v , and this leads to the following:

$$\begin{equation}\frac{1}{\langle k\rangle }\sum _{h=1}^{n}\frac{{\mathcal{I}}_{1}-{\mathcal{I}}_{2}}{{\left\{\left(\tau {\mathcal{R}}_{v}h{\Theta}+\gamma +\mu \right)\left(\alpha +\mu \right){{\Delta}}_{1}\left(h,{\Theta}\right)\right\}}^{2}}=0,\end{equation} \tag{ 10 }$$

where

$$\begin{align*}\hfill {\mathcal{I}}_{1}& =\left(\tau {\mathcal{R}}_{v}h{\Theta}+\gamma +\mu \right)\left(\alpha +\mu \right){{\Delta}}_{1}\left(h,{\Theta}\right){\times}\left\{\tau {h}^{2}P\left(h\right)\left(\left(1-\left(1-{\varepsilon }_{A}\right)p\right)\mu +\alpha \right)\right.\hfill \\ \hfill & \quad {\times}{{\Delta}}_{1}\left(h,{\Theta}\right)+\tau {\mathcal{R}}_{v}{h}^{2}P\left(h\right)\left(\left(1-\left(1-{\varepsilon }_{A}\right)p\right)\mu +\alpha \right){\varepsilon }_{L}\tau h\left({\Theta}+{\mathcal{R}}_{v}\frac{\partial {\Theta}}{\partial {\mathcal{R}}_{v}}\right)\hfill \\ \hfill & \quad +{\varepsilon }_{L}\tau {h}^{2}P\left(h\right)\left(1-{\varepsilon }_{A}\right)p\mu {{\Delta}}_{2}\left(h,{\Theta}\right)\hfill \\ \hfill & \left.\quad +{\varepsilon }_{L}\tau {\mathcal{R}}_{v}{h}^{2}P\left(h\right)\left(1-{\varepsilon }_{A}\right)p\mu \tau h\left({\Theta}+{\mathcal{R}}_{v}\frac{\partial {\Theta}}{\partial {\mathcal{R}}_{v}}\right)\right\},\hfill \end{align*}$$

and

$$\begin{align*}\hfill {\mathcal{I}}_{2}& =\tau {\mathcal{R}}_{v}{h}^{2}P\left(h\right)\left\{ \left(\left(1-\left(1-{\varepsilon }_{A}\right)p\right)\mu +\alpha \right){{\Delta}}_{1}\left(h,{\Theta}\right)+{\varepsilon }_{L}\left(1-{\varepsilon }_{A} \right)p\mu {{\Delta}}_{2}\left(h,{\Theta}\right) \right\}\hfill \\ \hfill & \quad {\times}\left\{\tau h\left({\Theta}+{\mathcal{R}}_{v}\frac{\partial {\Theta}}{\partial {\mathcal{R}}_{v}}\right)\left(\alpha +\mu \right){{\Delta}}_{1}\left(h,{\Theta}\right)\right.\hfill \\ \hfill & \quad +\left.\left(\tau {\mathcal{R}}_{v}h{\Theta}+\gamma +\mu \right)\left(\alpha +\mu \right){\varepsilon }_{L}\tau h\left({\Theta}+{\mathcal{R}}_{v}\frac{\partial {\Theta}}{\partial {\mathcal{R}}_{v}}\right)\right\}.\hfill \end{align*}$$

If ${\mathcal{R}}_{v}=1$ and Θ = 0, then

$$\begin{align*}\hfill {\mathcal{I}}_{1}& =\tau {h}^{2}P\left(h\right)\left(\gamma +\mu \right){\left(\alpha +\mu \right)}^{3}\left\{\left(1-\left(1-{\varepsilon }_{A}\right)p\right)\mu +\alpha +{\varepsilon }_{L}\left(1-{\varepsilon }_{A}\right)p\mu \right\}\hfill \\ \hfill & \quad +{\tau }^{2}{h}^{3}P\left(h\right)\left(\gamma +\mu \right){\left(\alpha +\mu \right)}^{2}\left\{\left(\left(1-\left(1-{\varepsilon }_{A}\right)p\right)\mu +\alpha \right){\varepsilon }_{L}+{\varepsilon }_{L}\left(1-{\varepsilon }_{A}\right)p\mu \right\}\frac{\partial {\Theta}}{\partial {\mathcal{R}}_{v}},\hfill \end{align*}$$

and

$$\begin{equation*}{\mathcal{I}}_{2}={\tau }^{2}{h}^{3}P\left(h\right){\left(\alpha +\mu \right)}^{2}\left\{\left(1-\left(1-{\varepsilon }_{A}\right)p\right)\mu +\alpha +{\varepsilon }_{L}\left(1-{\varepsilon }_{A}\right)p\mu \right\}\enspace \left\{\left(\alpha +\mu \right)+\left(\gamma +\mu \right){\varepsilon }_{L}\right\}\frac{\partial {\Theta}}{\partial {\mathcal{R}}_{v}}.\end{equation*}$$

Consequently, equation (10) becomes

$$\begin{equation*}\frac{1}{\langle k\rangle }\sum _{h=1}^{n}\frac{{\mathcal{I}}_{1}-{\mathcal{I}}_{2}}{{\left\{\left(\gamma +\mu \right){\left(\alpha +\mu \right)}^{2}\right\}}^{2}}=0.\end{equation*}$$

This implies that

$$\begin{align*}\hfill \sum _{h=1}^{n}\left({\mathcal{I}}_{1}-{\mathcal{I}}_{2}\right)& =\tau \langle {k}^{2}\rangle \left(\gamma +\mu \right){\left(\alpha +\mu \right)}^{3}\left\{\left(1-\left(1-{\varepsilon }_{A}\right)p\right)\mu +\alpha +{\varepsilon }_{L}\left(1-{\varepsilon }_{A}\right)p\mu \right\}\hfill \\ \hfill & \quad +{\tau }^{2}\langle {k}^{3}\rangle \left(\gamma +\mu \right){\left(\alpha +\mu \right)}^{2}\left\{\left(\left(1-\left(1-{\varepsilon }_{A}\right)p\right)\mu +\alpha \right){\varepsilon }_{L}\right.\hfill \\ \hfill & \quad +\left.{\varepsilon }_{L}\left(1-{\varepsilon }_{A}\right)p\mu \right\}\frac{\partial {\Theta}}{\partial {\mathcal{R}}_{v}}-{\tau }^{2}\langle {k}^{3}\rangle {\left(\alpha +\mu \right)}^{2}\left\{\left(1-\left(1-{\varepsilon }_{A}\right)p\right)\mu \right.\hfill \\ \hfill & \quad +\left.\alpha +{\varepsilon }_{L}\left(1-{\varepsilon }_{A}\right)p\mu \right\}\left\{\left(\alpha +\mu \right)+\left(\gamma +\mu \right){\varepsilon }_{L}\right\}\frac{\partial {\Theta}}{\partial {\mathcal{R}}_{v}}\hfill \\ \hfill & =\tau \langle {k}^{2}\rangle \left(\gamma +\mu \right){\left(\alpha +\mu \right)}^{3}\left\{\mu -\left(1-{\varepsilon }_{L}\right)\left(1-{\varepsilon }_{A}\right)p\mu +\alpha \right\}\hfill \\ \hfill & \quad +{\tau }^{2}\langle {k}^{3}\rangle \left(\gamma +\mu \right){\left(\alpha +\mu \right)}^{3}{\varepsilon }_{L}\frac{\partial {\Theta}}{\partial {\mathcal{R}}_{v}}-{\tau }^{2}\langle {k}^{3}\rangle {\left(\alpha +\mu \right)}^{3}\hfill \\ \hfill & \quad {\times}\left\{1-\left(1-{\varepsilon }_{L}\right)\left(1-{\varepsilon }_{A}\right)\left(1-{\varepsilon }_{W}\right)p\right\}\left\{\left(\alpha +\mu \right)+\left(\gamma +\mu \right){\varepsilon }_{L}\right\}\frac{\partial {\Theta}}{\partial {\mathcal{R}}_{v}}\hfill \\ \hfill & =0.\hfill \end{align*}$$

This leads to the following equality:

$$\begin{align}\hfill & \langle {k}^{2}\rangle \left(\gamma +\mu \right)\left\{\mu -\left(1-{\varepsilon }_{L}\right)\left(1-{\varepsilon }_{A}\right)p\mu +\alpha \right\}\hfill \\ \hfill & \qquad =\tau \langle {k}^{3}\rangle \left\{\left(1-\left(1-{\varepsilon }_{L}\right)\left(1-{\varepsilon }_{A}\right)\left(1-{\varepsilon }_{W}\right)p\right)\left(\left(\alpha +\mu \right)+\left(\gamma +\mu \right){\varepsilon }_{L}\right)\right.\hfill \\ \hfill & \qquad \quad -\left.\left(\gamma +\mu \right){\varepsilon }_{L}\right\}\frac{\partial {\Theta}}{\partial {\mathcal{R}}_{v}}.\hfill \end{align} \tag{ 11 }$$

It should be noted that the left-hand side of equation (11) is positive. Therefore, we can conclude that the sign of $\partial {\Theta}/\partial {\mathcal{R}}_{v}$ at $\left({\mathcal{R}}_{v},{\Theta}\right)=\left(1,0\right)$ is the same as the sign of

$$\begin{equation*}\left(1-\left(1-{\varepsilon }_{L}\right)\left(1-{\varepsilon }_{A}\right)\left(1-{\varepsilon }_{W}\right)p\right)\left(\left(\alpha +\mu \right)+\left(\gamma +\mu \right){\varepsilon }_{L}\right)-\left(\gamma +\mu \right){\varepsilon }_{L}.\end{equation*}$$

Consequently, the following necessary and sufficient condition for backward bifurcation can be made:

Zoom In Zoom Out Reset image size

Theorem 3.1. System (1) has backward bifurcation at ${\mathcal{R}}_{v}=1$ if and only if

$$\begin{equation}\left(1-\left(1-{\varepsilon }_{L}\right)\left(1-{\varepsilon }_{A}\right)\left(1-{\varepsilon }_{W}\right)p\right)\left(\left(\alpha +\mu \right)+\left(\gamma +\mu \right){\varepsilon }_{L}\right)-\left(\gamma +\mu \right){\varepsilon }_{L}{< }0.\end{equation} \tag{ 12 }$$

Remark 2. From theorem 3.1, we have some observations.

• (a)  
  If p = 0 or ɛ_L = 0, then condition (12) fails. Therefore, only forward bifurcation occurs at ${\mathcal{R}}_{v}=1$. This indicates that if there is no vaccine (i.e., p = 0), then there is only forward bifurcation in our model. However, when a vaccine is present and the protection of the vaccine is not 100% effective, we realize that leaky vaccine plays a vital role, leading to a backward bifurcation.
• (b)  
  The bifurcation condition (12) is independent of the network structure for our model.
• (c)  
  The bifurcation condition (12) is a necessary and sufficient condition while the conditions derived in [7, 21] are sufficient conditions. This is because we determine the slope of the endemic equilibrium curve at ${\mathcal{R}}_{v}=1$, and in [7, 21], the authors find the conditions to determine the number of endemic equilibria.

In this section, we study the persistence of the disease by using theorem 4.6 in [34]. It should be noted that the persistence of the disease has the same significance as the global stability of the endemic equilibrium in an epidemiological sense. We briefly sketch the idea of the proof. Define two subsets of Ω, namely,

$$\begin{equation*}{X}_{1}=\left\{\left({S}_{k}\left(t\right),{I}_{k}\left(t\right),{V}_{k}\left(t\right)\right)\in {\Omega}\vert I\left(t\right)=\sum _{h=1}^{n}P\left(h\right){I}_{h}\left(t\right){ >}0\right\}\end{equation*}$$

and

$$\begin{equation*}{X}_{2}={\Omega}{\backslash}{X}_{1}.\end{equation*}$$

The first step is to show that the set $\left\{{\mathcal{E}}_{0}\right\}$ is a covering of the union of $\left\{\omega \left({\mathcal{X}}_{0}\right)\vert {\mathcal{X}}_{0}\in {X}_{2},\enspace \phi \left(t,{\mathcal{X}}_{0}\right)\in {X}_{2},\quad \forall \enspace t{ >}0\right\}$, where $\phi \left(t,{\mathcal{X}}_{0}\right)$ and $\omega \left({\mathcal{X}}_{0}\right)$ are the solution and the omega limit set of the solution of system (1) with initial value ${\mathcal{X}}_{0}\in {\Omega}$, respectively. The second step is to claim that ${\mathcal{E}}_{0}$ is a weak repeller for X₁, i.e.,

$$\begin{equation*}\underset{t\to \infty }{\mathrm{lim sup}}\enspace d\left(\phi \left(t,{\mathcal{X}}_{0}\right),{\mathcal{E}}_{0}\right){ >}0\quad \mathrm{f}\mathrm{o}\mathrm{r}\enspace \mathrm{a}\mathrm{n}\mathrm{y}\enspace {\mathcal{X}}_{0}\in {X}_{1}.\end{equation*}$$

Here $d\left(\phi \left(t,{\mathcal{X}}_{0}\right),{\mathcal{E}}_{0}\right)$ denotes the distance between $\phi \left(t,{\mathcal{X}}_{0}\right)$ and ${\mathcal{E}}_{0}$. Owing to theorem 4.6 in [34], the conclusion of persistence can be made. The detailed proof is given in appendix B.

Theorem 4.1. If ${\mathcal{R}}_{v}{ >}1$, then there exists an η > 0, which is independent of the initial condition satisfying I(0) > 0, such that

$$\begin{equation*}\underset{t\to \infty }{\mathrm{lim inf}}\enspace I\left(t\right)=\underset{t\to \infty }{\mathrm{lim inf}}\sum _{h=1}^{n}P\left(h\right){I}_{h}\left(t\right){ >}\eta ,\end{equation*}$$

where $I\left(0\right)={\sum }_{h=1}^{n}P\left(h\right){I}_{h}\left(0\right)$.

In this section, we present some numerical simulations to illustrate our theoretical results. In the first example, we considered the stability of the disease-free equilibrium ${\mathcal{E}}_{0}$.

Example 5.1. In this example, we choose scale-free networks which satisfy a power-law degree distribution, namely, P(k) = η₁ k^−2.5 for k = 1, 2, ..., 100. The constant η₁ is chosen to make ${\sum }_{k=1}^{n}P\left(k\right)=1$. Therefore, we can get ⟨k⟩ = 1.7995 and ⟨k²⟩ = 13.8643. The parameters of system (1) are selected as β = 0.004, p = 0.85, μ = 0.015, γ = 0.05, ɛ_L = 0.2, and ɛ_A = ɛ_W = 0. Therefore, ${\mathcal{R}}_{v}=0.1517{< }1$ and ${\mathcal{R}}_{0}=0.4741{< }1$ can be established. It can be realized from theorem 2.2 that the disease-free equilibrium ${\mathcal{E}}_{0}$ is globally asymptotically stable. We use 30 different initial values to plot the trajectories of S_k (t), I_k (t), and V_k (t) for k = 30, 60, 90. As we can see in figure 2, all the trajectories approach to ${\mathcal{E}}_{0}$. This could support that the disease-free equilibrium ${\mathcal{E}}_{0}$ is globally stable.

Zoom In Zoom Out Reset image size

Example 5.2. In this example, we consider bifurcation. The following parameters are fixed: β = 0.001, μ = 0.015, and γ = 0.05. In figure 3, we plot phase diagrams with different values of ɛ_A and ɛ_W in the parameter domain (ɛ_L , p) to illustrate condition (12). In the red region, it means that condition (12) is satisfied, and hence, backward bifurcation occurs. On the contrary, forward bifurcation takes place when the parameters are in the blue region. Now, we focus on the left panel of figure 3. If we put p = 0.85 and ɛ_L = 0.2, then we can verify that

$$\begin{align*}\hfill & \left(1-\left(1-{\varepsilon }_{L}\right)\left(1-{\varepsilon }_{A}\right)\left(1-{\varepsilon }_{W}\right)p\right)\left(\left(\alpha +\mu \right)+\left(\gamma +\mu \right){\varepsilon }_{L}\right)-\left(\gamma +\mu \right){\varepsilon }_{L}\hfill \\ \hfill & \qquad =0.009-0.013{< }0.\hfill \end{align*}$$

This indicates that condition (12) is met. It can, therefore, be deduced from theorem 3.1 that backward bifurcation occurs at ${\mathcal{R}}_{v}=1$. However, if we change the value of ɛ_L from 0.2 to 0.05, then

$$\begin{align*}\hfill & \left(1-\left(1-{\varepsilon }_{L}\right)\left(1-{\varepsilon }_{A}\right)\left(1-{\varepsilon }_{W}\right)p\right)\left(\left(\alpha +\mu \right)+\left(\gamma +\mu \right){\varepsilon }_{L}\right)-\left(\gamma +\mu \right){\varepsilon }_{L}\hfill \\ \hfill & \qquad =0.0035-0.0033{ >}0.\hfill \end{align*}$$

Therefore, forward bifurcation occurs at ${\mathcal{R}}_{v}=1$. Because condition (12) is independent of the network structure, we will consider four different types of networks to demonstrate this result. The first one is the scale-free network that we have chosen in example 5.1; the second one is the network with an exponential degree of distribution, namely, P(k) = η₂ e^−k/5 for k = 1, 2, ..., 100; the third one is the Poisson network with degree distribution P(k) = η₃ e⁻¹⁵15^k /k! for k = 1, 2, ..., 100; the final one is a scientific collaboration network [26] whose degree of distribution is P(k) = η₄ k^−1.3 e^−k/10.7 for k = 1, 2, ..., 100. The constants η_i (i = 2, 3, 4) are selected to keep each degree of distribution satisfying ${\sum }_{k=1}^{\mathrm{100}}P\left(k\right)=1$. The bifurcation diagrams for these networks are given in figures 4–7. The numerical results confirm our theoretical analysis.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Example 5.3. In this example, let us focus on the scale-free networks with the degree of distribution described in example 5.1. The parameters are chosen as β = 0.0025, p = 0.85, μ = 0.015, γ = 0.05, ɛ_L = 0.2, and ɛ_A = ɛ_W = 0. Therefore, we have ${\mathcal{R}}_{v}=0.9482{< }1$ and ${\mathcal{R}}_{0}=2.9632{ >}1$. In addition, as we have verified in example 5.2, the condition (12) is satisfied. Therefore, backward bifurcation occurs at ${\mathcal{R}}_{v}=1$. Given that ${\mathcal{R}}_{v}=0.9482$, it can be observed in the right subfigure of figure 4 that there exist two positive values of Θ, and each of them determines an endemic equilibrium respectively. That is, under the aforementioned parameters, there are two endemic equilibria and one of them is stable. To confirm this observation, 30 different initial values are used to show the time evolutions of S_k (t), I_k (t), and V_k (t) for k = 30, 60, 90. A bistable phenomenon can be observed in figure 8.

Zoom In Zoom Out Reset image size

Example 5.4. This example is aimed at studying the effects of imperfect vaccines on the backward bifurcation behavior. We fix the parameters β = 0.01, μ = 0.015, γ = 0.2, p = 0.85, and ɛ_L = 0.2. We will consider the following four cases: (i) ɛ_A = ɛ_W = 0; (ii) ɛ_A = 0.2 and ɛ_W = 0; (iii) ɛ_W = 0.2 and ɛ_A = 0; (iv) ɛ_A = ɛ_W = 0.2. The bifurcation diagrams for the different types of network structures are depicted in figure 9, where the bifurcation curves corresponding to case (i)–(iv) are blue, red, black, and magenta respectively. We can observe that the region of metastability depends on the presence of the type of imperfect vaccine. More specifically, the region of backward bifurcation is the largest for case (i), and it is smallest for case (iv). Therefore, it seems that we need the strongest condition to eliminate the disease if there is only a leaky vaccine. However, in the presence of all the three imperfect vaccines, the weakest condition is needed.

Zoom In Zoom Out Reset image size

Example 5.5. In the final example, we will conduct simulations to illustrate the persistence of the result. We will still use the scale-free networks with the degree of distribution described in example 5.1. The model parameters are selected as β = 0.03, p = 0.85, μ = 0.015, γ = 0.05 with the 3-tuple (ɛ_L , ɛ_A , ɛ_W ) varying. We use 30 different initial values to depict the trajectories of I(t) with different settings of the 3-tuple (ɛ_L , ɛ_A , ɛ_W ) in figure 10. As be observed from each panel of figure 9, the threshold value is ${\mathcal{R}}_{v}{ >}1$. Therefore, the persistence of the disease can be assured from theorem 4.1. In fact, we can see that all the trajectories approach to a positive constant. This indicates that there exists a stable endemic equilibrium when ${\mathcal{R}}_{v}{ >}1$.

Zoom In Zoom Out Reset image size

In this study, we investigated the dynamics of a network-based epidemic model with imperfect vaccinations. We considered three types of imperfect vaccinations, namely, leaky, all-or-nothing, and waning vaccines. We proved that there exists a threshold value ${\mathcal{R}}_{v}$ that determines the stability of the disease-free equilibrium ${\mathcal{E}}_{0}$. More specifically, we proved that if ${\mathcal{R}}_{v}{< }1$, then ${\mathcal{E}}_{0}$ is locally asymptotically stable. Furthermore, because of the presence of imperfect vaccinations, we derived a necessary and sufficient condition that ensures the occurrence of a backward (or subcritical) bifurcation at ${\mathcal{R}}_{v}=1$. From our result, we found that a leaky vaccine plays a vital role in leading to a backward bifurcation at ${\mathcal{R}}_{v}=1$. When such a bifurcation occurs, ${\mathcal{R}}_{v}{< }1$ does not confirm the global stability of ${\mathcal{E}}_{0}$. This means that the condition ${\mathcal{R}}_{v}{< }1$ is insufficient to eradicate the disease. When ${\mathcal{R}}_{v}{ >}1$, we showed that the disease can persist. Several numerical experiments were performed to validate the theoretical results. We conclude that when a backward bifurcation occurs, we require a stronger condition to eliminate the diseases. In comparison with the three types of imperfect vaccinations, the leaky vaccine must be paid special attention as it explicitly affects the backward bifurcation phenomena, which are challenging from the perspective of disease control. Therefore, it is important to study how to control the spreading dynamics when such a bifurcation has been detected. Furthermore, in this work, we only considered a simple epidemic model on static networks. However, it is known that the real contact networks are time-varying. Hence, it is also necessary to study the spreading dynamics on real contact networks. Besides, we know that disease spread is a fundamentally spatiotemporal phenomenon. Thus, studying the effects of space-time heterogeneities on the spreading dynamics of infectious diseases is also an important and interesting issue. We defer these research topics for future investigation.

The authors would like to thank the anonymous referees for their valuable comments and suggestions that improved the presentation of this paper. This work was partially supported by the Ministry of Science and Technology of Taiwan under the Grant MOST 108-2115-M-017-002.

In this appendix, we give the detailed proofs of the stability results in section 2. The following lemma is used to facilitate the stability analysis.

Lemma A.1. ([36]) For a real n × n matrix A = (a_ij ) where a_ij = δ_ij σ_i + p_i q_j (p_i , q_j ⩾ 0, i, j = 1, 2, ..., n) and δ_ij denotes the Kronecker symbol. Then,

$$\begin{equation*}\mathrm{det}\left(A\right)={\sigma }_{1}{\sigma }_{2}\dots {\sigma }_{n}+{p}_{1}{q}_{1}{\sigma }_{2}{\sigma }_{3}\dots {\sigma }_{n}+{\sigma }_{1}{p}_{2}{q}_{2}{\sigma }_{3}\dots {\sigma }_{n}+\cdots +{\sigma }_{1}{\sigma }_{2}\dots {\sigma }_{n-1}{p}_{n}{q}_{n}.\end{equation*}$$

Specifically, if σ_k = σ for all k, then

$$\begin{equation}\mathrm{det}\left(A\right)={\sigma }^{n-1}\left(\sigma +\sum _{i=1}^{n}{p}_{i}{q}_{i}\right).\end{equation} \tag{ A.1 }$$

Proof of theorem 2.1. To determine the local stability of ${\mathcal{E}}_{0}$ by linearization, we rearrange the equations in system (1) as the following system:

$$\begin{equation}{\boldsymbol{y}}^{\prime }\left(t\right)=f\left(\boldsymbol{y}\left(t\right)\right),\end{equation} \tag{ A.2 }$$

where $\boldsymbol{y}\left(t\right)={\left({y}_{1}\left(t\right),{y}_{2}\left(t\right),\dots ,{y}_{3n}\left(t\right)\right)}^{\top }$ and $f\left(\boldsymbol{y}\left(t\right)\right)={\left({f}_{1}\left(\boldsymbol{y}\left(t\right)\right),{f}_{2}\left(\boldsymbol{y}\left(t\right)\right),\dots ,{f}_{3n}\left(\boldsymbol{y}\left(t\right)\right)\right)}^{\top }$ with y_k (t) = S_k (t), y_n+k (t) = I_k (t), y_2n+k (t) = V_k (t), and

$$\begin{align*}\hfill {f}_{k}\left(\boldsymbol{y}\left(t\right)\right)& =\left(1-\left(1-{\varepsilon }_{A}\right)p\right)\mu -\beta k{y}_{k}\left(t\right){\Theta}\left(t\right)+\gamma {y}_{n+k}\left(t\right)+\alpha {y}_{2n+k}\left(t\right)-\mu {y}_{k}\left(t\right),\hfill \\ \hfill {f}_{n+k}\left(\boldsymbol{y}\left(t\right)\right)& =\beta k{y}_{k}\left(t\right){\Theta}\left(t\right)+{\varepsilon }_{L}\beta k{y}_{2n+k}\left(t\right){\Theta}\left(t\right)-\left(\gamma +\mu \right){y}_{n+k}\left(t\right),\hfill \\ \hfill {f}_{2n+k}\left(\boldsymbol{y}\left(t\right)\right)& =\left(1-{\varepsilon }_{A}\right)p\mu -{\varepsilon }_{L}\beta k{y}_{2n+k}\left(t\right){\Theta}\left(t\right)-\left(\alpha +\mu \right){y}_{2n+k}\left(t\right),\quad k=1,2,\dots ,n.\hfill \end{align*}$$

Then, we can obtain the Jacobian matrix of system (A.2) at ${\mathcal{E}}_{0}$

$$\begin{equation*}J\left({\mathcal{E}}_{0}\right)=\left(\begin{matrix}\hfill {A}_{1}\hfill & \hfill \star \hfill & \hfill \star \hfill \\ \hfill O\hfill & \hfill {A}_{2}\hfill & \hfill O\hfill \\ \hfill O\hfill & \hfill \star \hfill & \hfill {A}_{3}\hfill \end{matrix}\right),\end{equation*}$$

where A₁ = diag{−μ, −μ, ..., −μ} and A₃ = diag{−(α + μ), −(α + μ), ..., −(α + μ)} are two n × n diagonal matrices; O is the n × n zero matrix; A₂ is an n × n matrix whose (i, j)-entry is given by

$$\begin{equation*}{\left({A}_{2}\right)}_{ij}=\begin{cases}\frac{1}{\langle k\rangle }\beta ijP\left(j\right)\left({S}_{0}+{\varepsilon }_{L}{V}_{0}\right)-\left(\gamma +\mu \right),\quad \hfill & i=j,\hfill \\ \frac{1}{\langle k\rangle }\beta ijP\left(j\right)\left({S}_{0}+{\varepsilon }_{L}{V}_{0}\right),\quad \hfill & i\ne j.\hfill \end{cases}\end{equation*}$$

It can be observed that matrix A₁ has a negative eigenvalue −μ with multiplicity n, and matrix A₃ has a negative eigenvalue −(α + μ) with multiplicity n. Therefore, the other n eigenvalues are determined by the matrix A₂. We define the (i, j)-entries of matrix A₂ − λI as

$$\begin{equation*}-\left(\gamma +\mu +\lambda \right){\delta }_{ij}+\frac{1}{\langle k\rangle }\beta ijP\left(j\right)\left({S}_{0}+{\varepsilon }_{L}{V}_{0}\right).\end{equation*}$$

Denote

$$\begin{equation*}{\sigma }_{i}=-\left(\gamma +\mu +\lambda \right),\quad {p}_{i}=\beta i\left({S}_{0}+{\varepsilon }_{L}{V}_{0}\right),\quad \mathrm{a}\mathrm{n}\mathrm{d}\quad {q}_{j}=\frac{jP\left(j\right)}{\langle k\rangle }.\end{equation*}$$

From equation (A.1), the following characteristic polynomial can be made:

$$\begin{equation*}\mathrm{det}\left({A}_{2}-\lambda I\right)={\left(-1\right)}^{n-1}{\left(\gamma +\mu +\lambda \right)}^{n-1}\left(-\left(\gamma +\mu +\lambda \right)+\frac{1}{\langle k\rangle }\left({S}_{0}+{\varepsilon }_{L}{V}_{0}\right)\sum _{i=1}^{n}\beta {i}^{2}P\left(i\right)\right).\end{equation*}$$

By solving the characteristic equation det(A₂ − λI) = 0, we can obtain a negative characteristic root λ = −(γ + d) that has multiplicity n − 1. Consequently, the stability of the disease-free equilibrium ${\mathcal{E}}_{0}$ completely depends on the sign of the root of

$$\begin{equation*}-\left(\gamma +\mu +\lambda \right)+\frac{1}{\langle k\rangle }\left({S}_{0}+{\varepsilon }_{L}{V}_{0}\right)\sum _{i=1}^{n}\beta {i}^{2}P\left(i\right)=0.\end{equation*}$$

This implies that

$$\begin{align*}\hfill \lambda & =-\left(\gamma +\mu \right)+\left({S}_{0}+{\varepsilon }_{L}{V}_{0}\right)\frac{\beta \langle {k}^{2}\rangle }{\langle k\rangle }\hfill \\ \hfill & =\left(\gamma +\mu \right)\left(\left({S}_{0}+{\varepsilon }_{L}{V}_{0}\right)\frac{\beta \langle {k}^{2}\rangle }{\left(\gamma +\mu \right)\langle k\rangle }-1\right)\hfill \\ \hfill & =\left(\gamma +\mu \right)\left({\mathcal{R}}_{v}-1\right).\hfill \end{align*}$$

Accordingly, if ${\mathcal{R}}_{v}{< }1$, then all the eigenvalues of $J\left({\mathcal{E}}_{0}\right)$ are negative. Therefore, ${\mathcal{E}}_{0}$ is asymptotically stable. Conversely, if ${\mathcal{R}}_{v}{ >}1$, then ${\mathcal{E}}_{0}$ is unstable because of the existence of a positive eigenvalue. This completes the proof. □

Proof of theorem 2.2. It suffices to show that ${\mathcal{E}}_{0}$ is globally attractive, and by combining the local stability result (theorem 2.1), the global stability of ${\mathcal{E}}_{0}$ can be achieved. To this end, first, we claim that I_k (t) → 0 when t → ∞ for all k. According to the second equation of system (1), we can obtain the derivative of Θ(t):

$$\begin{align*}\hfill {{\Theta}}^{\prime }\left(t\right)& =\frac{1}{\langle k\rangle }\sum _{h=1}^{n}hP\left(h\right){I}_{h}^{\prime }\left(t\right)\hfill \\ \hfill & =\frac{1}{\langle k\rangle }\sum _{h=1}^{n}hP\left(h\right)\left\{\beta h{S}_{h}\left(t\right){\Theta}\left(t\right)+{\varepsilon }_{L}\beta h{V}_{h}\left(t\right){\Theta}\left(t\right)-\left(\gamma +\mu \right){I}_{h}\left(t\right)\right\}.\hfill \end{align*}$$

Given that ɛ_L ∈ [0, 1) and Ω is a positive invariant set, we can derive that

$$\begin{align*}\hfill {{\Theta}}^{\prime }\left(t\right)& {\leqslant}\frac{1}{\langle k\rangle }\sum _{h=1}^{n}hP\left(h\right)\left\{\beta h\left(1-{I}_{h}\left(t\right)\right){\Theta}\left(t\right)-\left(\gamma +\mu \right){I}_{h}\left(t\right)\right\}\hfill \\ \hfill & {\leqslant}\frac{\beta \langle {k}^{2}\rangle }{\langle k\rangle }{\Theta}\left(t\right)-\left(\gamma +\mu \right){\Theta}\left(t\right)\hfill \\ \hfill & =\left(\gamma +\mu \right){\Theta}\left(t\right)\left({\mathcal{R}}_{0}-1\right).\hfill \end{align*}$$

Because ${\mathcal{R}}_{0}{< }1$, we can obtain Θ(t) → 0 as t → ∞. This implies that I_k (t) → 0 can be considered as t → ∞ for all k. Therefore, for any δ > 0, there exists T_δ > 0 such that 0 ⩽ I_k (t) ⩽ δ for t > T_δ . Next, we claim V_k (t) → V₀ as t → ∞ for all k. It follows from the third equation of system (1) that

$$\begin{equation*}{V}_{k}^{\prime }\left(t\right){\leqslant}\left(1-{\varepsilon }_{A}\right)p\mu -\left(\alpha +\mu \right){V}_{k}\left(t\right).\end{equation*}$$

This implies that

$$\begin{equation*}\underset{t\to \infty }{\mathrm{lim sup}}\enspace {V}_{k}\left(t\right){\leqslant}\frac{\left(1-{\varepsilon }_{A}\right)p\mu }{\alpha +\mu }={V}_{0}.\end{equation*}$$

In contrast, because 0 ⩽ I_k (t) ⩽ δ for t > T_δ , we have

$$\begin{align*}\hfill {V}_{k}^{\prime }\left(t\right)& {\geqslant}\left(1-{\varepsilon }_{A}\right)p\mu -\delta {\varepsilon }_{L}\beta k{V}_{k}\left(t\right)-\left(\alpha +\mu \right){V}_{k}\left(t\right)\hfill \\ \hfill & =\left(1-{\varepsilon }_{A}\right)p\mu -\left(\delta {\varepsilon }_{L}\beta k+\alpha +\mu \right){V}_{k}\left(t\right),\quad t{ >}{T}_{\delta }.\hfill \end{align*}$$

Considering δ → 0, we can obtain

$$\begin{equation*}\underset{t\to \infty }{\mathrm{lim inf}}\enspace {V}_{k}\left(t\right){\geqslant}\frac{\left(1-{\varepsilon }_{A}\right)p\mu }{\alpha +\mu }={V}_{0}.\end{equation*}$$

Therefore, V_k (t) → V₀ when t → ∞ for all k. In fact, S_k (t) = 1 − I_k (t) − V_k (t) leads to S_k (t) → S₀ when t → ∞ for all k. Consequently, ${\mathcal{E}}_{0}$ is globally attractive whenever ${\mathcal{R}}_{0}{< }1$, and the proof is completed. □

Proof. Define

$$\begin{equation*}{X}_{1}=\left\{\left({S}_{k}\left(t\right),{I}_{k}\left(t\right),{V}_{k}\left(t\right)\right)\in {\Omega}\vert I\left(t\right)=\sum _{h=1}^{n}P\left(h\right){I}_{h}\left(t\right){ >}0\right\}\end{equation*}$$

and

$$\begin{equation*}{X}_{2}={\Omega}{\backslash}{X}_{1}.\end{equation*}$$

First, we claim that X₁ is a positive invariant with respect to system (1). If (S_k (0), I_k (0), V_k (0)) ∈ X₁, then I(0) > 0 and

$$\begin{align*}\hfill {I}^{\prime }\left(t\right)& =\sum _{h=1}^{n}P\left(h\right){I}_{h}^{\prime }\left(t\right)\hfill \\ \hfill & =\sum _{h=1}^{n}P\left(h\right)\left(\beta h{S}_{h}\left(t\right){\Theta}\left(t\right)+{\varepsilon }_{L}\beta h{V}_{h}\left(t\right){\Theta}\left(t\right)-\left(\gamma +\mu \right){I}_{h}\left(t\right)\right)\hfill \\ \hfill & {\geqslant}-\left(\gamma +\mu \right)I\left(t\right).\hfill \end{align*}$$

Therefore, I(t) ⩾ I(0) exp(−(γ + μ)t) is considered for t > 0, making X₁ a positive invariant. Furthermore, there is a compact set B in which all the solutions of system (1) initiated in X₁ will enter and remain forever. The compactness condition (C_4.2) of theorem 4.6 in [34] is easily verified for set B.

Let $\phi \left(t,{\mathcal{X}}_{0}\right)$ be the solution of system (1) with the initial condition ${\mathcal{X}}_{0}\in {\Omega}$, and $\omega \left({\mathcal{X}}_{0}\right)$ is the ω-limit set of the solution of system (1) with an initial value of ${\mathcal{X}}_{0}\in {\Omega}$. Denote

$$\begin{equation*}Y=\cup \left\{\omega \left({\mathcal{X}}_{0}\right)\vert {\mathcal{X}}_{0}\in {X}_{2},\enspace \phi \left(t,{\mathcal{X}}_{0}\right)\in {X}_{2},\quad \forall \enspace t{ >}0\right\}.\end{equation*}$$

Restricting system (1) on ${M}_{\partial }=\left\{{\mathcal{X}}_{0}\vert \phi \left(t,{\mathcal{X}}_{0}\right)\in {X}_{2},\quad \forall \enspace t{\geqslant}0\right\}$ gives

$$\begin{equation}\begin{cases}_{k}^{\prime }\left(t\right)=\left(1-\left(1-{\varepsilon }_{A}\right)p\right)\mu +\alpha {V}_{k}\left(t\right)-\mu {S}_{k}\left(t\right),\quad \hfill \\ {I}_{k}^{\prime }\left(t\right)=-\left(\gamma +\mu \right){I}_{k}\left(t\right),\quad \hfill \\ {V}_{k}^{\prime }\left(t\right)=\left(1-{\varepsilon }_{A}\right)p\mu -\left(\alpha +\mu \right){V}_{k}\left(t\right),\quad k=1,2,\dots ,n.\quad \hfill \end{cases}\end{equation} \tag{ B.1 }$$

Clearly, the disease-free equilibrium ${\mathcal{E}}_{0}$ of system (1) is the only equilibrium of system (B.1) in M_∂. One can easily claim that ${\mathcal{E}}_{0}$ is globally asymptotically stable in M_∂ because (B.1) is a linear system. Therefore, $Y=\left\{{\mathcal{E}}_{0}\right\}$. Moreover, ${\mathcal{E}}_{0}$ is a covering of Y. This is isolated and acyclic. Finally, the proof will be completed if we show that ${\mathcal{E}}_{0}$ is a weak repeller for X₁, i.e.,

$$\begin{equation*}\underset{t\to \infty }{\mathrm{lim sup}}\enspace d\left(\phi \left(t,{\mathcal{X}}_{0}\right),{\mathcal{E}}_{0}\right){ >}0\quad \mathrm{f}\mathrm{o}\mathrm{r}\enspace \mathrm{a}\mathrm{n}\mathrm{y}\enspace {\mathcal{X}}_{0}\in {X}_{1}.\end{equation*}$$

Therefore, we need to verify that ${W}^{s}\left({\mathcal{E}}_{0}\right)\cap {X}_{1}=\varnothing$, where

$$\begin{equation*}{W}^{s}\left({\mathcal{E}}_{0}\right)=\left\{{\mathcal{X}}_{0}\in {X}_{1}\vert \enspace \underset{t\to \infty }{\mathrm{lim}}\enspace \phi \left(t,{\mathcal{X}}_{0}\right)={\mathcal{E}}_{0}\right\}.\end{equation*}$$

Supposing this is not true, then there exists a solution within the initial condition ${\mathcal{X}}_{0}\in {X}_{1}$ such that

$$\begin{equation*}{S}_{k}\left(t\right)\to {S}_{0},\enspace {I}_{k}\left(t\right)\to 0,\enspace {V}_{k}\left(t\right)\to {V}_{0}\quad \mathrm{a}\mathrm{s}\enspace t\to \infty ,\quad k=1,2,\dots ,n.\end{equation*}$$

Given that ${\mathcal{R}}_{v}{ >}1$, we can choose a small > 0 such that

$$\begin{equation*}\zeta {:=}\left({S}_{0}+{\varepsilon }_{L}{V}_{0}-\left(1+{\varepsilon }_{L}\right){\epsilon}\right)\frac{\beta \langle {k}^{2}\rangle }{\left(\gamma +\mu \right)\langle k\rangle }{ >}1.\end{equation*}$$

Considering > 0, there exists T > 0 such that for t > T, and we have

$$\begin{equation*}{S}_{0}-{\epsilon}{\leqslant}{S}_{k}\left(t\right){\leqslant}{S}_{0}+{\epsilon},\enspace 0{\leqslant}{I}_{k}\left(t\right){\leqslant}{\epsilon},{V}_{0}-{\epsilon}{\leqslant}{V}_{k}\left(t\right){\leqslant}{V}_{0}+{\epsilon},\quad k=1,2,\dots ,n.\end{equation*}$$

Consider the function $L\left(t\right)=\frac{1}{\left(\gamma +\mu \right)}{\Theta}\left(t\right)$. The derivative of L along solution (S_k (t), I_k (t), V_k (t)) for t > T is given by

$$\begin{align*}\hfill {L}^{\prime }\left(t\right)& =\frac{1}{\left(\gamma +\mu \right)\langle k\rangle }\sum _{h=1}^{n}hP\left(h\right){I}_{h}^{\prime }\left(t\right)\hfill \\ \hfill & =\frac{1}{\left(\gamma +\mu \right)\langle k\rangle }\sum _{h=1}^{n}hP\left(h\right)\left(\beta h{S}_{h}\left(t\right){\Theta}\left(t\right)+{\varepsilon }_{L}\beta h{V}_{h}\left(t\right){\Theta}\left(t\right)-\left(\gamma +\mu \right){I}_{h}\left(t\right)\right)\hfill \\ \hfill & {\geqslant}\frac{1}{\left(\gamma +\mu \right)\langle k\rangle }\sum _{h=1}^{n}hP\left(h\right)\left(\beta h\left({S}_{0}-{\epsilon}\right){\Theta}\left(t\right)+{\varepsilon }_{L}\beta h\left({V}_{0}-{\epsilon}\right){\Theta}\left(t\right)-\left(\gamma +\mu \right){I}_{h}\left(t\right)\right)\hfill \\ \hfill & =\frac{1}{\left(\gamma +\mu \right)\langle k\rangle }\sum _{h=1}^{n}\beta {h}^{2}P\left(h\right)\left({S}_{0}+{\varepsilon }_{L}{V}_{0}-\left(1+{\varepsilon }_{L}\right){\epsilon}\right){\Theta}\left(t\right)-\left(\gamma +\mu \right)L\left(t\right)\hfill \\ \hfill & =\left({S}_{0}+{\varepsilon }_{L}{V}_{0}-\left(1+{\varepsilon }_{L}\right){\epsilon}\right)\frac{\beta \langle {k}^{2}\rangle }{\left(\gamma +\mu \right)\langle k\rangle }{\Theta}\left(t\right)-\left(\gamma +\mu \right)L\left(t\right)\hfill \\ \hfill & =\left(\zeta -1\right){\Theta}\left(t\right)\hfill \\ \hfill & =\left(\zeta -1\right)\left(\gamma +\mu \right)L\left(t\right).\hfill \end{align*}$$

Because ζ > 1, we can conclude that L(t) → ∞ can be considered as t → ∞. This contradicts the boundedness of L(t). Therefore, ${W}^{s}\left({\mathcal{E}}_{0}\right)\cap {X}_{1}=\varnothing$. The proof is completed. □

In this appendix, we conduct parameters-related sensitivity and uncertainty analyses on ${\mathcal{R}}_{v}$ and on the bifurcation condition (12) since the dynamical behaviors of system (1) are based on these two quantities. We will use Latin hypercube sampling-based method to quantify uncertainty and sensitivity of ${\mathcal{R}}_{v}$ and the left-hand side of condition (12) as a function of 7 parameters (β, γ, μ, p, ɛ_L , ɛ_A , ɛ_W ) and a function of 6 parameters (γ, μ, p, ɛ_L , ɛ_A , ɛ_W ), respectively. The Latin hypercube sampling is a type of Monte Carlo sampling [13] as it densely stratifies the input parameters. The partial rank correlation coefficient (PRCC) is used as the sensitivity index to measure the impact of the parameters on the output variables, that is, ${\mathcal{R}}_{v}$ and the left-hand side of condition (12). More specifically, the PRCC measures the strength between the input parameters and outputs of the model (i.e., ${\mathcal{R}}_{v}$ and condition (12)) using correlation through the sampling done by the Latin hypercube sampling method [23]. Because we have no a priori data, a uniform distribution is assigned to each parameter [23]. After using the Latin hypercube sampling method, 2000 samples are taken and the histograms of the input parameter values are drawn in figure 11. Corresponding to the above mentioned distributions for input parameters, the resulting distributions for ${\mathcal{R}}_{v}$ and for the left-hand side of condition (12) are made in figure 12. The distribution of ${\mathcal{R}}_{v}$ is on the left panel of figure 12. In this figure, the probability of ${\mathcal{R}}_{v}{ >}1$ is about 0.328. The distribution for the left-hand side of condition (12) is on the right panel of figure 12. In this figure, the probability of negative values of condition (12) is about 0.045. The PRCC sensitivity results are carried out for ${\mathcal{R}}_{v}$ and for condition (12), respectively. The bar charts are given in figure 13. PRCC suggests that the parameters β and μ have the highest and the lowest influence on ${\mathcal{R}}_{v}$, respectively. The parameters ɛ_W and ɛ_L have the highest and the lowest influence on the bifurcation condition (12), respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size