Effect of a Vaccination against the Dengue Fever Epidemic in an Age Structure Population: From the Perspective of the Local and Global Stability Analysis

Abstract

The effect of vaccination on the dengue fever epidemic described by an age structured modified SIR (Susceptible-Infected-Retired) model is studied using standard stability analysis. The chimeric yellow fever dengue tetravalent dengue vaccine (CYD-TDV™) is a vaccine recently developed to control this epidemic in several Southeast Asian countries. The dengue vaccination program requires a total of three injections, 6 months apart at 0, 6, and 12 months. The ages of the recipients are nine years and above. In this paper, we analyze the mathematical dynamics SIR transmission model of the epidemic. The stability of the model is established using Routh–Hurwitz criteria to see if a Hopf Bifurcation occurs and see when the equilibrium states are local asymptotically stable or global asymptotically stable. We have determined the efficiency of CYD-TDV by simulating the optimal numerical solution for each age range for this model. The numerical results showed the optimal age for vaccination and significantly reduced the severity and severity of the disease.

1. Introduction

Dengue disease is an endemic disease in many tropical countries, especially in the urban areas in these countries. This disease is a mosquito-borne viral infection, caused by four serotypes of the dengue virus, namely DEN-1, DEN-2, DEN-3, and DEN-4 [1,2,3]. These dengue viruses are transmitted to human beings by the bite of an infected Aedes aegypti mosquito. When a parameter known as the basic reproduction number $R_0$ is equal or greater than one, it means that an outbreak has occurred in the population [4]. If the epidemic occurs in five continents at the same time, we have a pandemic.

General symptoms of dengue fever (DF) include acute high fevers, headache, pains around the eye socket and muscle pains. Dengue hemorrhagic fever (DHF) may have more severe symptoms which are nosebleeds, vomiting blood, and plasma leakage and in some cases could include hypotension, anuria, and shock, which is called dengue shock syndrome (DSS) [5,6]. Infection by one strain of the dengue virus provides immunity to further infection by that strain but not to the others. In fact, a sequential infection by a second strain leads to an allergic reaction to the strain and this may lead to DHF and then to DSS (most likely ending in death). To understand how any disease (or epidemic) develops, a mathematical model is developed to describe the dynamics of the epidemic. Having a mathematical model, we have a way to predict what would happen if we changed the parameters in the model. It should be remembered that the greatest contribution to the public health of humans was the realization by Ross [7] that malaria was a disease that kills over three million people a year, and occurs when a susceptible person is bitten by an infected (a different type from the one transmitting the dengue virus) mosquito and becomes infected by the parasite. When this newly infectious mosquito bites another healthy human, the bitten human in turn becomes an infectious human can infect a second mosquito and begin the cycle again. Under certain conditions, the basic reproduction number of the model will be equal or greater than one, upon which, the epidemic occurs.

To formulate the model, we note that the World Health Organization (WHO) reported that the lifetime of the Aedes mosquitoes is about 10 days, with a flying distance of about 100 m. The mosquitoes’ eggs hatch in 5–8 days [3,8]. They have also pointed out that an estimated 3.9 billion people in 128 counties are at risk of the disease. While DF does not lead to the death of the person infected with the disease or the likelihood of death is very low, at the same time, DSS does lead to the death of people It does make people sick who then cannot work. So, there is a great economic impact. Therefore, groups around the world have looked for a vaccine that can be used against this disease.

As we have pointed out, dengue fever has had a great economic impact on many countries in South East Asia (such as Thailand and the Philippines). In the late 80′s, Mahidol University (Bangkok) in Thailand began research on developing a tetravalent dengue vaccine. A tetravalent vaccine was sought because all four strains are present in Thailand. In reality, only three (DEN-1, DEN-2, and DEN-3) were prevalent (DEN-4 only arrived recently) [9]. In 2015, Sanofi Pasteur (a French Drug company) released a drug chimeric yellow fever dengue tetravalent dengue vaccine (CYD-TDV) as a dengue vaccine. It offered (i.e., it appeared to work in year 1 and 2) hope initially but problems arose in year 3 causing it now to be recommended only for people forty-five years or older. Chanprasopchai, Tang, and Pongsumpun [10] recently looked at the effects of vaccinations on the DF epidemic in a population consisting of a single population of humans divided into susceptible, infectious and recovered. There was no age structure in the human population. In that paper, the model predicted that without the vaccination, there would be a period of oscillatory behavior in the number of infected people, which would die out. With the vaccination, there would be two periods of oscillatory behaviors (See Figure 10a,b in ref. [10]). However, the manufacturer of the vaccine CYD-TDV emphasizes the importance of taking into account the ages of the people being vaccinated; recommending vaccination only for people between the ages 9–45 years.

Aguiar et al. [11] studied the impact of the CYD-TDV in endemic countries. The outcomes of two vaccination strategies were simulated. The two strategies were first for people in the ages group 9–45 years, who are seropositive for one of the strains and seronegative for the others, and second for people ages 9–45 years who were immune to an individual dengue strain. They divided the population into three groups, such as susceptible, infective, and recovered by considering four age ranges are less than 9 years with natural infections, between 9 to 45 years and natural infections, between 9 to 45 years and vaccine, and more than 45 to 65 years with natural infections. The models of Hadinegoro [12] and Coudeville [13] reported that the effectiveness of the vaccine can reduce hospitalization by 80.8%. The vaccine efficacies against dengue hemorrhagic fever were 92.9%. This vaccine will be effective in a population aged 9 years and older. The reported efficacy of CYD-TDV, which was tested in healthy children to receive three vaccinations, divided into months 0, 6, and 12 was in the third phase of the vaccine study. It was found that the appropriate age for injection was from the age of 9 years and older [14,15,16,17]. Chanprasochai et al. and Esteva and Vargas [10,18,19] studied the dynamical transmission using the Susceptible-Infected-Recovered (SIR) model, including vaccination, with standard dynamical systems analyses. Driessche and Watmough [20] have studied the reproductive number and its relation to the endemic equilibria of the disease transmission model. It was generally found that if $R_0<1$, the disease-free equilibrium point is locally stable, and the endemic equilibrium point is stable at $R_0>1$. In 2021, the authors [21] studied the optimal control strategies through the SEIRS model by considering only the pre-exposure infectious transmission status and the impact of post-reinfection vaccination strategies without regarding the serotype sequence in infection. Sanusi et al. [22] studied analyzing a SIRS model and also simulating the model to predict the number of dengue fever cases. The results of the SIRS model for dengue fever can predict the number of dengue cases in South Sulawesi. The authors [23] proposed a compartmental model of dengue disease that included the typical incidence relation between susceptible vectors and infected humans in order to examine the impact of model parameters on the basic reproduction number. The results demonstrate that the vector population’s bite rate has a greater impact on the basic reproduction number than the other parameters.

In the present model, we include the efficacies of the vaccination to the different strains of the dengue virus. The difference in efficacy takes into account that infection by one strain only gives partial immunity to the others. The presence of a weakened immunity allows for allergies to the other strains which might result in the more severe form of DF (DHF), as well as antibody dependent enhancement (ADE). According to data from the Ministry of Public Health of Thailand from 2003 to 2019, the number of DF patients in each year can be divided into age groups as shown in Figure 1. It is seen that the numbers of the DF patients are more pronounced in three age groups, namely, patients aged less than 10 years, aged between 10 and 44 years, and above 45. Figure 2 plots the number of infected cases each year from 2003–2019 for each of the three age groups [24]. The significant number of occurrences for the cluster of 10–44 years suggests that an age structural model would be appropriate in order to control the disease in Thailand [24]. Note, that the plotted data suggests that there is a rise in the DF cases in the last decade, which is attributable to the arrival of the DEN-4 strain in Thailand. The developed model, thus, takes into account the efficacy of the vaccination to the different strains, which have not been investigated before.

In this study, we consider the modified SIR mathematical model for dengue disease with a natural and chimeric yellow fever dengue tetravalent dengue vaccine (CYD-TDV) in humans. The total population was subdivided into three groups: patients aged less than 10 years, between 10–44 years of age, and greater than 44 years of age. Stability and bifurcation analyses are then given, followed by numerical simulations.

2. Methodology

2.1. Mathematical Model

In this section, a mathematical model governing the dynamics of the dengue disease is developed, taking into account the vector-host dynamics. The total human population is divided into three groups according to age. Each age group then incorporates a SIR structure, where each of the susceptible and the infected classes are sub-classified according to the efficacy of the applied vaccine, which is assumed for the purpose of our modeling to be constant. The state variables and parameters in Equations (1)–(7) are defined in

Table 1. The variable and parameters used in the system equation.

Variable and Parameters	References
$\overline{S_{H_{i_a}}}$	The number of the susceptible human populations
$\overline{I_{H_{i_a}}}$	The number of the infected human populations
$\overline{R_{H_{i_a}}}$	The number of the recovered human populations
$\overline{S_{H_{i_a}V}}$	The number of the susceptible human population with vaccines
$\overline{I_{H_{i_a}V}}$	The number of the infected human population with vaccines
$\overline{S_V}$	The number of the susceptible vector populations
$\overline{I_V}$	The number of the infected vector populations
$P_{i_a}$	$\mathrm{The}\mathrm{efficacy}\mathrm{of}\mathrm{vaccination}\mathrm{for}\mathrm{all}{i}_a=1,2,3$
${\gamma }_{i_a}$	$\mathrm{Recovery}\mathrm{rate}\mathrm{of}\mathrm{human}\mathrm{and}\mathrm{with}\mathrm{vaccine}\mathrm{populations}\mathrm{for}\mathrm{all}{i}_a=1,2,3$
${\alpha }_{i_a}$	$\mathrm{The}\mathrm{presence}\mathrm{of}\mathrm{the}\mathrm{efficacy}\mathrm{coefficient}\mathrm{for}\mathrm{all}{i}_a=1,2,3$
${\beta }_{H_{i_a}}$	$\mathrm{Transmission}\mathrm{rate}\mathrm{of}\mathrm{dengue}\mathrm{virus}\mathrm{from}\mathrm{vector}\mathrm{to}\mathrm{human}\mathrm{for}\mathrm{all}{i}_a=1,2,3$
${\beta }_V$	Transmission rate of dengue virus from human to vector
$A$	Constant recruitment rate of vector populations
$b_H$	Birth rate of the human populations
$d_H$	Natural mortality rate of the human populations
$d_V$	Natural mortality rate of the vector populations
$d_d$	Mortality rate from the infected of the human populations
$d_k$	Mortality rate from the infected of the vector populations
$N_H$	Total human populations
$N_V$	Total vector populations

.

The vector transmission admits a simple SI structure. Note, that the overall structure of the transmission model is now a modified SIR model. The schematics of the transmission are shown in Figure 3.

The transmission of the disease is shown in Figure 3 by the following system of differential equations which are the mathematical representations of the time rates of changes of the state variables (flow charts of the state variables):

$$\frac{d\overline{S_{H_{i_a}}}}{dt}=\left(1-P_{i_a}\right)b_H{N}_H-{\beta }_{H_{i_a}}\overline{S_{H_{i_a}}}\overline{I_V}-d_H\overline{S_{H_{i_a}}}$$

(1)

$$\frac{d\overline{I_{H_{i_a}}}}{dt}={\beta }_{H_{i_a}}\overline{S_{H_{i_a}}}\overline{I_V}-{\gamma }_{i_a}\overline{I_{H_{i_a}}}-\left(d_H+d_d\right)\overline{I_{H_{i_a}}}$$

(2)

$$\frac{d\overline{R_{H_{i_a}}}}{dt}={\gamma }_{i_a}\left(\overline{I_{H_{i_a}}}+\overline{I_{H_{i_a}V}}\right)-d_H\overline{R_{H_{i_a}}}$$

(3)

$$\frac{d\overline{S_{H_{i_a}V}}}{dt}=P_{i_a}{b}_H{N}_H-\left(1-{\alpha }_{i_a}\right){\beta }_{H_{i_a}}\overline{S_{H_{i_a}V}}\overline{I_V}-d_H\overline{S_{H_{i_a}V}}$$

(4)

$$\frac{d\overline{I_{H_{i_a}V}}}{dt}=\left(1-{\alpha }_{i_a}\right){\beta }_{H_{i_a}}\overline{S_{H_{i_a}V}}\overline{I_V}-{\gamma }_{i_a}\overline{I_{H_{i_a}V}}-\left(d_H+d_d\right)\overline{I_{H_{i_a}V}}$$

(5)

$$\frac{d\overline{S_V}}{dt}=A-{\beta }_V\overline{S_V}\left({\displaystyle \sum _{i_a=1}^3\overline{I_{H_{i_a}}}}+{\displaystyle \sum _{i_a=1}^3\overline{I_{H_{i_a}V}}}\right)-d_V\overline{S_V}$$

(6)

$$\frac{d\overline{I_V}}{dt}={\beta }_V\overline{S_V}\left({\displaystyle \sum _{i_a=1}^3\overline{I_{H_{i_a}}}}+{\displaystyle \sum _{i_a=1}^3\overline{I_{H_{i_a}V}}}\right)-\left(d_V+d_k\right)\overline{I_V}$$

(7)

For all $i_a=1,2,3$ where $i_a=1$: age between 0 to 9 years, $i_a=2$: age between 10 to 44 years, and $i_a=3$: age more than 45 years.

$${\displaystyle \sum _{i_a=1}^3\left(\overline{S_{H_{i_a}}}+\overline{I_{H_{i_a}}}+\overline{R_{H_{i_a}}}+\overline{S_{H_{i_a}V}}+\overline{I_{H_{i_a}V}}\right)}=N_H$$

(8)

$$\overline{S_V}+\overline{I_V}=N_V$$

(9)

Normalizing the equations by introducing the following normalized variables:

$$\begin{array}{l}\overline{S_{H_{i_a}}}=S_{H_{i_a}}{N}_H,\overline{I_{H_{i_a}}}=I_{H_{i_a}}{N}_H,\overline{R_{H_{i_a}}}=R_{H_{i_a}}{N}_H,\overline{S_{H_{i_a}V}}=S_{H_{i_a}V}{N}_H,\overline{I_{H_{i_a}V}}=I_{H_{i_a}V}{N}_H,\\ \overline{S_V}=S_V{N}_V,\overline{I_V}=I_V{N}_V\end{array}$$

(10)

Since the infection rate does not introduce exogenous deaths into our population, and with a minimal timespan, the population is thus assumed to be constant. In terms of the normalized variables, the assumption leads to:

$${\displaystyle \sum _{i_a=1}^3\left(S_{H_{i_a}}+I_{H_{i_a}}+R_{H_{i_a}}+S_{H_{i_a}V}+I_{H_{i_a}V}\right)}=1$$

(11)

$$S_V+I_V=1$$

(12)

The mathematical model of Equations (1)–(7) is now reduced to the following equations. In terms of the normalized variables, the original system of dynamical transmission becomes:

$$\frac{d{S}_{H_{i_a}}}{dt}=\left(1-P_{i_a}\right)b_H-{\beta }_{H_{i_a}}{S}_{H_{i_a}}{I}_V{N}_V-d_H{S}_{H_{i_a}}$$

(13)

$$\frac{d{I}_{H_{i_a}}}{dt}={\beta }_{H_{i_a}}{S}_{H_{i_a}}{I}_V{N}_V-{\gamma }_{i_a}{I}_{H_{i_a}}-\left(d_H+d_d\right)I_{H_{i_a}}$$

(14)

$$\frac{d{S}_{H_{i_a}V}}{dt}=P_{i_a}{b}_H-\left(1-{\alpha }_{i_a}\right){\beta }_{H_{i_a}}{S}_{H_{i_a}V}{I}_V{N}_V-d_H{S}_{H_{i_a}V}$$

(15)

$$\frac{d{I}_{H_{i_a}V}}{dt}=\left(1-{\alpha }_{i_a}\right){\beta }_{H_{i_a}}{S}_{H_{i_a}V}{I}_V{N}_V-{\gamma }_{i_a}{I}_{H_{i_a}V}-\left(d_H+d_d\right)I_{H_{i_a}V}$$

(16)

$$\frac{d{I}_V}{dt}={\beta }_V{S}_V{N}_H\left({\displaystyle \sum _{i_a=1}^3{I}_{H_{i_a}}}+{\displaystyle \sum _{i_a=1}^3{I}_{H_{i_a}V}}\right)-\left(d_V+d_k\right)I_V$$

(17)

For all $i_a=1,2,3$.

The Non-Negativeness of the Solutions

Proposition 1.

Let $\left(S_{H_{i_a}}(t),I_{H_{i_a}}(t),S_{H_{i_a}V}(t),I_{H_{i_a}V}(t),R_{H_{i_a}}(t),S_V(t),I_V(t)\right)$be the solutions of Equations (1)–(7) for all $i_a=1,2,3$. Denoting also invariant set $\varphi =\left\{\left(S_{H_{i_a}},I_{H_{i_a}},S_{H_{i_a}V},I_{H_{i_a}V},R_{H_{i_a}},S_V,I_V\right)\in R_{+}^7:N_T\le \frac{b_H}{d_H}+\frac{A}{d_V}\right\}$. Then the closed set is positive invariant.

Proof. 

We begin by setting:

$$N_T(t)=N_H(t)+N_V(t)={\displaystyle \sum _{i=1}^3\left(S_{H_{i_a}},I_{H_{i_a}},S_{H_{i_a}V},I_{H_{i_a}V},R_{H_{i_a}}\right)}+S_V+I_V,$$

and assume that $N_T=\frac{b_H}{d_H}+\frac{A}{d_V}$. The total populations’ $N_T$ is a non-negative definite on $R_{+}^7$. Differentiation yields:

$$\frac{d{N}_T}{dt}=\frac{d{N}_H}{dt}+\frac{d{N}_V}{dt}$$

$$=\frac{d}{dt}{\displaystyle \sum _{i_a=1}^3\left(S_{H_{i_a}}+I_{H_{i_a}}+S_{H_{i_a}V}+I_{H_{i_a}V}+R_{H_{i_a}}\right)}+\frac{d}{dt}\left(S_V+I_V\right)$$

$$\begin{array}{l}=\left(b_H-d_H{S}_{H_{i_a}}-d_H{I}_{H_{i_a}}-d_H{S}_{H_{i_a}V}-d_H{I}_{H_{i_a}V}-d_H{R}_{H_{i_a}}-d_d{I}_{H_{i_a}}-d_d{I}_{H_{i_a}V}\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\left(A-d_V{S}_V-d_V{I}_V-d_k{I}_V\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\left(b_H-d_H{N}_H-d_d\left(I_{H_{i_a}}+I_{H_{i_a}V}\right)\right)+\left(A-d_V{N}_V-d_k{I}_V\right)\end{array}$$

Integrating the above equation, it follows that $\frac{d{N}_T}{dt}=\frac{d{N}_H}{dt}+\frac{d{N}_V}{dt}\le 0$ on $0\le N_T(t)\le \left(N_H(0)e^{-d_{H}t}+\frac{b_H}{d_H}\left(1-e^{-d_{H}t}\right)\right)+\left(N_V(0)e^{-d_{V}t}+\frac{A}{d_V}\left(1-e^{-d_{V}t}\right)\right)$.

As $t\to \infty ,e^{-d_{H}t}\to 0,$ and $e^{-d_{V}t}\to 0$ we have $\underset{t\to \infty }{\lim }{N}_T(t)=\underset{t\to \infty }{\lim }{N}_H(t)+\underset{t\to \infty }{\lim }{N}_V(t)\le \frac{b_H}{d_H}+\frac{A}{d_V}$. $N_T(t)$ thus approaches $\frac{b_H{d}_V+A{d}_H}{d_H{d}_V}$. Since the region of all solutions of $\varphi$ is in $R_{+}^7$ [25,26]. □

2.2. Stability Analysis

2.2.1. Steady States and Basic Reproduction Number

The steady states are obtained by setting the right hand side of Equations (13)–(17) to zero. The resultant disease-free steady state $E_1^{\ast }$ is then given as:

$$E_1^{\ast }(t)=\left(E_{11},E_{12},E_{13},0,0,0,E_{21},E_{22},E_{23},0,0,0,0\right)$$

(18)

where

$$E_{1i_a}=\frac{b_H\left(1-P_{i_a}\right)}{d_H},\mathrm{and}{E}_{2i_a}=\frac{b_H{P}_{i_a}}{d_H}\mathrm{for}\mathrm{all}{i}_a=1,2,3$$

Likewise, the endemic steady state $E_2^{\ast }$ is given:

$$E_2^{\ast }(t)=\left(S_{H_1}^{\ast },S_{H_2}^{\ast },S_{H_3}^{\ast },I_{H_1}^{\ast },I_{H_2}^{\ast },I_{H_3}^{\ast },S_{H_{1}V}^{\ast },S_{H_{2}V}^{\ast },S_{H_{3}V}^{\ast },I_{H_{1}V}^{\ast },I_{H_{2}V}^{\ast },I_{H_{3}V}^{\ast },I_V^{\ast }\right)$$

(19)

where

$$S_{H_{i_a}}^{\ast }=\frac{{\epsilon }_5{\epsilon }_8{\epsilon }_{10}{\beta }_{H_{i_a}}+d_H{\epsilon }_{11}\left({\epsilon }_7+{\alpha }_{i_a}\left({\epsilon }_2{\epsilon }_4-{\epsilon }_7\left(P_{i_a}-2\right)\right)\right)+{\epsilon }_{12}}{2{\epsilon }_1\left(d_H+{\epsilon }_{11}\right)}$$

$$I_{H_{i_a}}^{\ast }=\frac{{\epsilon }_8\left(1+{\alpha }_{i_a}\left(1-2P_{i_a}\right)\right){\epsilon }_{11}{\beta }_{H_{i_a}}-d_H{\epsilon }_{11}\left({\alpha }_{i_a}{P}_{i_a}{\epsilon }_7+{\epsilon }_2{\epsilon }_4-{\epsilon }_7\right)-{\epsilon }_{12}}{2{\epsilon }_{10}{N}_H{N}_V{\alpha }_{i_a}\left(d_H+{\epsilon }_{11}\right){\epsilon }_4}$$

$$S_{H_{i_a}V}^{\ast }=\frac{d_H{\epsilon }_{11}\left({\epsilon }_7+{\alpha }_{i_a}\left({\epsilon }_7{P}_{i_a}+{\epsilon }_2{\epsilon }_4\right)\right)-{\epsilon }_5{\epsilon }_8{\epsilon }_{10}{\beta }_{H_{i_a}}-{\epsilon }_{12}}{2{\epsilon }_1\left(d_H-{\epsilon }_5{\epsilon }_{11}\right)}$$

$$I_{H_{i_a}V}^{\ast }=\frac{{\epsilon }_{12}-{\epsilon }_5{\epsilon }_8\left(2P_{i_a}{\alpha }_{i_a}-1\right){\epsilon }_{10}{\beta }_{H_{i_a}}-d_H{\epsilon }_{11}\left({\epsilon }_7+{\alpha }_{i_a}\left({\epsilon }_2{\epsilon }_4-{\epsilon }_7{P}_{i_a}\right)\right)}{2N_H{N}_V{\alpha }_{i_a}{\epsilon }_{10}\left(d_H-{\epsilon }_5{\epsilon }_{11}\right){\epsilon }_4}$$

$$I_V^{\ast }=\frac{{\epsilon }_5{\epsilon }_8{\epsilon }_{10}{\beta }_{H_{i_a}}+d_H{\epsilon }_{11}\left({\epsilon }_3{\epsilon }_7-d_d{\epsilon }_2{\epsilon }_6-d_H{\epsilon }_2{\epsilon }_6-{\gamma }_{i_a}{\epsilon }_2{\epsilon }_6\right)+{\epsilon }_{12}}{2{\epsilon }_5{\epsilon }_{10}^2\left({\epsilon }_7+{\epsilon }_2{\epsilon }_4\right)}$$

For all $i_a=1,2,3$ with:

$$\begin{array}{l}{\epsilon }_1=d_H{N}_H{N}_V{\alpha }_{i_a}{\beta }_{i_a}{\beta }_V,{\epsilon }_2=\left(d_V+d_k\right),{\epsilon }_3=\left(P_{i_a}{\alpha }_{i_a}-1\right),{\epsilon }_4=\left(d_d+d_H+{\gamma }_{i_a}\right),\\ {\epsilon }_5=\left({\alpha }_{i_a}-1\right),{\epsilon }_6=\left({\alpha }_{i_a}-2\right),{\epsilon }_7=b_H{N}_H{\beta }_V,{\epsilon }_8=b_H{N}_H{N}_V^2,{\epsilon }_9=b_H{N}_H{N}_V,\\ {\epsilon }_{10}={\beta }_V{\beta }_{i_a},{\epsilon }_{11}=N_V{\beta }_{i_a},{\epsilon }_{12}={\epsilon }_{11}^2\left(4d_H{\epsilon }_5\left({\epsilon }_7+{\epsilon }_2{\epsilon }_4\right)\left({\epsilon }_2{\epsilon }_9{\epsilon }_{10}+d_H{\epsilon }_2{\epsilon }_4\right)\right),\\ {\epsilon }_{13}={\epsilon }_{11}^2{\left(d_d{d}_H{\epsilon }_2{\epsilon }_6+d_H^2{\epsilon }_2{\epsilon }_6-{\epsilon }_5{\epsilon }_9{\epsilon }_{10}+d_H\left({\epsilon }_3{\epsilon }_7+{\epsilon }_2{\epsilon }_6{\gamma }_{i_a}\right)\right)}^2,{\epsilon }_{14}=\sqrt{{\epsilon }_{12}+{\epsilon }_{13}}.\end{array}$$

(20)

The basic reproductive number $(R_0)$ is computed using the next generation matrix method. To apply this method, we have to rewrite the normalized time rate of change equations,

$$\frac{d{I}_{H_{i_a}}}{dt}={\beta }_{H_{i_a}}{S}_{H_{i_a}}{I}_V{N}_V-{\gamma }_{i_a}{I}_{H_{i_a}}-\left(d_H+d_d\right)I_{H_{i_a}}$$

$$\frac{d{I}_{H_{i_a}V}}{dt}=\left(1-{\alpha }_{i_a}\right){\beta }_{H_{i_a}}{S}_{H_{i_a}V}{I}_V{N}_V-{\gamma }_{i_a}{I}_{H_{i_a}V}-\left(d_H+d_d\right)I_{H_{i_a}V}$$

(21)

$$\begin{gathered} \frac{d{I}_V}{dt}={\beta }_V{S}_V{N}_H\left({\displaystyle \sum _{i_a=1}^3{I}_{H_{i_a}}}+{\displaystyle \sum _{i_a=1}^3{I}_{H_{i_a}V}}\right)-\left(d_V+d_k\right)I_V \\ \mathrm{for}\mathrm{all}{i}_a=1,2,3. \end{gathered}$$

Writing Equation (21) as separate equations, we have

$$\frac{d{I}_{H_1}}{dt}={\beta }_{H_1}{S}_{H_1}{I}_V{N}_V-{\gamma }_1{I}_{H_1}-\left(d_H+d_d\right)I_{H_1}$$

(22)

$$\frac{d{I}_{H_2}}{dt}={\beta }_{H_2}{S}_{H_2}{I}_V{N}_V-{\gamma }_2{I}_{H_2}-\left(d_H+d_d\right)I_{H_2}$$

(23)

$$\frac{d{I}_{H_3}}{dt}={\beta }_{H_3}{S}_{H_3}{I}_V{N}_V-{\gamma }_3{I}_{H_3}-\left(d_H+d_d\right)I_{H_3}$$

(24)

$$\frac{d{I}_{H_{1}V}}{dt}=\left(1-{\alpha }_1\right){\beta }_{H_1}{S}_{H_{1}V}{I}_V{N}_V-{\gamma }_1{I}_{H_{1}V}-\left(d_H+d_d\right)I_{H_{1}V}$$

(25)

$$\frac{d{I}_{H_{2}V}}{dt}=\left(1-{\alpha }_2\right){\beta }_{H_2}{S}_{H_{2}V}{I}_V{N}_V-{\gamma }_2{I}_{H_{2}V}-\left(d_H+d_d\right)I_{H_{2}V}$$

(26)

$$\frac{d{I}_{H_{3}V}}{dt}=\left(1-{\alpha }_3\right){\beta }_{H_3}{S}_{H_{3}V}{I}_V{N}_V-{\gamma }_3{I}_{H_{3}V}-\left(d_H+d_d\right)I_{H_{3}V}$$

(27)

$$\frac{d{I}_V}{dt}={\beta }_V{S}_V{N}_H\left(I_{H_1}+I_{H_2}+I_{H_3}+I_{H_{1}V}+I_{H_{2}V}+I_{H_{3}V}\right)-\left(d_V+d_k\right)I_V$$

(28)

In the matrix form, they become $\det \left(F{V}^{-1}-\eta I\right)=0$ where the matrix $F$ for new infections and the matrix $V$ (or losses) are:

$$F=\left[\begin{array}{ccccccc}0& 0& 0& 0& 0& 0& F_{17}\\ 0& 0& 0& 0& 0& 0& F_{27}\\ 0& 0& 0& 0& 0& 0& F_{37}\\ 0& 0& 0& 0& 0& 0& F_{47}\\ 0& 0& 0& 0& 0& 0& F_{57}\\ 0& 0& 0& 0& 0& 0& F_{67}\\ F_{71}& F_{72}& F_{73}& F_{74}& F_{75}& F_{76}& 0\end{array}\right]$$

(29)

where

$$\begin{array}{l}{F}_{17}={\beta }_{H_1}{S}_{H_1}{N}_V,F_{27}={\beta }_{H_2}{S}_{H_2}{N}_V,F_{37}={\beta }_{H_3}{S}_{H_3}{N}_V,\\ F_{47}=\left(1-{\alpha }_1\right){\beta }_{H_1}{S}_{H_{1}V}{N}_V,F_{57}=\left(1-{\alpha }_2\right){\beta }_{H_2}{S}_{H_{2}V}{N}_V,F_{67}=\left(1-{\alpha }_3\right){\beta }_{H_3}{S}_{H_{3}V}{N}_V,\\ F_{71}=F_{72}=F_{73}=F_{74}=F_{75}=F_{76}={\beta }_V{S}_V{N}_H.\end{array}$$

$$V=\left[\begin{array}{ccccccc}{V}_{11}& 0& 0& 0& 0& 0& 0\\ 0& V_{22}& 0& 0& 0& 0& 0\\ 0& 0& V_{33}& 0& 0& 0& 0\\ 0& 0& 0& V_{44}& 0& 0& 0\\ 0& 0& 0& 0& V_{55}& 0& 0\\ 0& 0& 0& 0& 0& V_{66}& 0\\ 0& 0& 0& 0& 0& 0& V_{77}\end{array}\right]$$

(30)

where

$$\begin{array}{l}{V}_{11}=V_{44}={\gamma }_1+d_H+d_d,\\ V_{22}=V_{55}={\gamma }_2+d_H+d_d,\\ V_{33}=V_{66}={\gamma }_3+d_H+d_d,\\ V_{77}=d_V+d_k.\end{array}$$

Then

$$F{V}^{-1}=\left[\begin{array}{ccccccc}0& 0& 0& 0& 0& 0& G_{17}\\ 0& 0& 0& 0& 0& 0& G_{27}\\ 0& 0& 0& 0& 0& 0& G_{37}\\ 0& 0& 0& 0& 0& 0& G_{47}\\ 0& 0& 0& 0& 0& 0& G_{57}\\ 0& 0& 0& 0& 0& 0& G_{67}\\ G_{71}& G_{72}& G_{73}& G_{74}& G_{75}& G_{76}& 0\end{array}\right]$$

(31)

where

$$\begin{array}{l}{G}_{17}=\frac{{\beta }_{H_1}{S}_{H_1}{N}_V}{d_V+d_k},G_{27}=\frac{{\beta }_{H_2}{S}_{H_2}{N}_V}{d_V+d_k},G_{37}=\frac{{\beta }_{H_3}{S}_{H_3}{N}_V}{d_V+d_k},\\ G_{47}=\frac{\left(1-{\alpha }_1\right){\beta }_{H_1}{S}_{H_{1}V}{N}_V}{d_V+d_k},G_{57}=\frac{\left(1-{\alpha }_2\right){\beta }_{H_2}{S}_{H_{2}V}{N}_V}{d_V+d_k},G_{67}=\frac{\left(1-{\alpha }_3\right){\beta }_{H_3}{S}_{H_{3}V}{N}_V}{d_V+d_k},\\ G_{71}=G_{74}=\frac{{\beta }_V{S}_V{N}_H}{{\gamma }_1+d_H+d_d},G_{72}=G_{75}=\frac{{\beta }_V{S}_V{N}_H}{{\gamma }_2+d_H+d_d},G_{73}=G_{76}=\frac{{\beta }_V{S}_V{N}_H}{{\gamma }_3+d_H+d_d},.\end{array}$$

Solving the matrix equation $\det \left(F{V}^{-1}-\eta I\right)=0$, we obtain the eigenvalues:

$$\begin{array}{l}{\eta }_1={\eta }_2={\eta }_3={\eta }_4={\eta }_5=0,\\ {\eta }_6=-\frac{\sqrt{m_2{m}_3+m_4+m_5+m_6+m_7+m_8}}{m_1},\\ {\eta }_7=\frac{\sqrt{m_2{m}_3+m_4+m_5+m_6+m_7+m_8}}{m_1},\end{array}$$

where

$$\begin{array}{l}{d}_1=d_H,d_2=d_d,d_3=d_V,d_4=d_k,\\ s_1=S_{H_1},s_2=S_{H_2},s_3=S_{H_3},s_4=S_{H_{1}V},s_5=S_{H_{2}V},s_6=S_{H_{3}V},\\ k_1=\left(d_3+d_4\right),k_2=\left(d_1+d_2+{\gamma }_1\right),k_3=\left(d_1+d_2+{\gamma }_2\right),\\ k_4=\left(d_1+d_2+{\gamma }_3\right),k_5=s_4\left({\alpha }_1-1\right),k_6=s_5\left({\alpha }_2-1\right),k_7=s_6\left({\alpha }_3-1\right),\\ k_8=\left({\gamma }_1+{\gamma }_2\right),k_9=\left({\gamma }_1+{\gamma }_3\right),k_{10}=\left({\gamma }_2+{\gamma }_3\right),\\ m_1=\frac{1}{\sqrt{k_1{k}_2{k}_3{k}_4}},m_2=N_H{N}_V{S}_V{\beta }_V,\\ m_3=d_1^2\left(\left(s_1-k_5\right){\beta }_1+\left(s_2-k_6\right){\beta }_2+\left(s_3-k_7\right){\beta }_3\right),\\ m_4=d_2^2\left(\left(s_1-k_5\right){\beta }_1+\left(s_2-k_6\right){\beta }_2+\left(s_3-k_7\right){\beta }_3\right),\\ m_5=\left(s_3-k_7\right){\beta }_3{\gamma }_1{\gamma }_2,m_6=\left(\left(s_2-k_6\right){\beta }_2{\gamma }_1+\left(s_1-k_5\right){\beta }_1{\gamma }_2\right){\gamma }_3,\\ m_7=d_1\left(\begin{array}{l}2d_2\left(\left(s_1-k_5\right){\beta }_1+\left(s_2-k_6\right){\beta }_2+\left(s_3-k_7\right){\beta }_3\right)-\\ s_5{d}_2{\beta }_2{\gamma }_1+s_3{\beta }_3{\gamma }_1+s_6{\beta }_3{\gamma }_1-s_6{\alpha }_3{\beta }_3{\gamma }_1+s_3{\beta }_3{\gamma }_2+s_6{\beta }_3{\gamma }_2-\\ s_6{\alpha }_3{\beta }_3{\gamma }_2-k_6{\beta }_2{\gamma }_3+s_2{\beta }_2{k}_9+\left(s_1+s_4-s_4{\alpha }_1\right){\beta }_1{k}_{10}\end{array}\right),\\ m_8=d_2\left(\begin{array}{l}\left(s_3+s_6-s_6{\alpha }_3\right){\beta }_3{k}_8+s_2{\beta }_2{k}_9-\\ k_6{\beta }_2{k}_9+\left(s_1+s_4-s_4{\alpha }_1\right){\beta }_1{k}_{10}\end{array}\right).\end{array}$$

Hence, $R_0=\rho \left(F{V}^{-1}\right)$ where $\rho =\max \left\{{\eta }_1,{\eta }_2,{\eta }_3,\dots {\eta }_7\right\}.$ Therefore

$$R_0=\frac{\sqrt{m_2{m}_3+m_4+m_5+m_6+m_7+m_8}}{m_1}$$

(32)

In addition, the number $R_0$ is called the basic reproductive number [8].

2.2.2. Local and Global Stabilities

It is well known that the steady states are local asymptotically stable if all the eigenvalues have negative real parts. The eigenvalue is obtained by solving the characteristic equation:

$$\det \left(J_{E_K}-\lambda I\right)=0$$

(33)

where $J_{E_K^{\ast }}$ was the Jacobian matrix of steady state $E_K^{\ast }$ for $K=1,2$. $I$ is the identity matrix of dimension $5\times 5$ and $\lambda$ was the eigenvalue. The Jacobian matrix of the system in Equations (13)–(17) for $i_a=1,2,3$ are as follows:

$$J_{E_K}=\left[\begin{array}{ccccc}-J_{11}-d_H& 0& 0& 0& -J_{15}\\ J_{21}& J_{22}& 0& 0& J_{25}\\ 0& 0& -J_{33}-d_H& 0& -J_{35}\\ 0& 0& J_{43}& J_{44}& J_{45}\\ 0& J_{52}& 0& J_{54}& J_{55}\end{array}\right].$$

(34)

where

$$\begin{array}{l}{J}_{11}=J_{22}={\beta }_{H_{i_a}}{I}_V{N}_V,J_{15}=J_{25}={\beta }_{H_{i_a}}{S}_{H_{i_a}}{N}_V,J_{22}=J_{44}=-{\gamma }_{i_a}-d_H-d_d,\\ J_{35}=J_{45}=\left(1-{\alpha }_{i_a}\right){\beta }_{H_{i_a}}{S}_{H_{i_a}}{N}_V,J_{52}=J_{54}={\beta }_V\left(1-I_V\right)N_H,\\ J_{55}=-d_V-d_k,J_{33}=J_{43}=\left(1-{\alpha }_{i_a}\right){\beta }_{H_{i_a}}{I}_V{N}_V.\end{array}$$

Proposition 2.

The disease-free steady state $E_1^{\ast }$ is local asymptotically stable when $R_0<1$.

Proof. 

For disease-free steady state $E_1^{\ast }$ the characteristic equation is given by $\det (J_{H_1^{\ast }}-\lambda I)=0$, thus the eigenvalues are:

$$\begin{array}{l}{\lambda }_1={\lambda }_2=-d_H,\\ {\lambda }_3=-\left(d_H+d_d+{\gamma }_{i_a}\right),\\ {\lambda }_4={\lambda }_5=-\left(\frac{d_H\left({\epsilon }_2+{\epsilon }_4\right)+\sqrt{d_H\left({\left({\epsilon }_4-{\epsilon }_2\right)}^2-4{\epsilon }_3{\epsilon }_9{\epsilon }_{10}\right)}}{2d_H}\right).\end{array}$$

When ${\epsilon }_2,{\epsilon }_3,{\epsilon }_4,{\epsilon }_9,$ and ${\epsilon }_{10}$ are defined in Equation (20). Hence, $E_1^{\ast }$ is local asymptotically stable. □

Next, we consider the local property of the endemic steady state of Equations (13)–(17).

Proposition 3.

The endemic steady state $E_2^{\ast }$ is local asymptotically stable when $R_0>1$.

Proof. 

For endemic steady state $E_2^{\ast }$. The characteristic equation is defined by:

$$\left(-\lambda -d_H-d_d-{\gamma }_{i_a}\right)\left({\lambda }^4+a_1{\lambda }^3+a_2{\lambda }^2+a_3\lambda +a_4\right)=0$$

(35)

where

$$\begin{array}{l}{a}_1=4d_H+{\epsilon }_2+2{\epsilon }_{11}{I}_V^{\ast }-{\epsilon }_{11}{\alpha }_{i_a}{I}_V^{\ast }+{\gamma }_i\\ a_2=3d_H^2+d_d\left(2d_H+{\epsilon }_2-{\epsilon }_6{\epsilon }_{11}{I}_V^{\ast }\right)-{\epsilon }_{11}\left(d_k{\epsilon }_6{I}_V^{\ast }+d_V{\epsilon }_6{I}_V^{\ast }\right)+\\ {\epsilon }_{11}\left(N_H\left(I_V^{\ast }-1\right)\left(S_{H_{i_a}}^{\ast }-S_{H_{i_a}V}^{\ast }{\epsilon }_5\right){\beta }_{H_{i_a}}\right)+\\ {\epsilon }_{11}\left(\left({\epsilon }_2-{\epsilon }_6{\epsilon }_{11}{I}_V^{\ast }\right){\gamma }_{i_a}+d_H\left(3{\epsilon }_2-2{\epsilon }_6{\epsilon }_{11}{I}_V^{\ast }+2{\gamma }_{i_a}\right)\right)\\ a_3=d_H^3+d_d\left(d_H\left(d_H+2{\epsilon }_2\right)-\left(d_H+{\epsilon }_2\right){\epsilon }_6{\epsilon }_{11}{I}_V^{\ast }-{\epsilon }_{11}^2{\epsilon }_5{I}_V^{{\ast }^2}\right)+\\ d_H^2\left(3{\epsilon }_2-{\epsilon }_6{\epsilon }_{11}{I}_V^{\ast }+{\gamma }_{i_a}\right)+d_H\left(2{\epsilon }_2-{\epsilon }_6{\epsilon }_{11}{I}_V^{\ast }{\gamma }_{i_a}\right)\\ d_H\left({\epsilon }_{11}\left(-2{\epsilon }_2{\epsilon }_6{I}_V^{\ast }+2N_H\left(I_V^{\ast }-1\right)\left(S_{H_{i_a}}^{\ast }-S_{H_{i_a}V}^{\ast }{\epsilon }_5\right){\beta }_V-{\epsilon }_5{\epsilon }_{11}{I}_V^{{\ast }^2}\right)\right)+\\ {\epsilon }_{11}{I}_V^{\ast }\left(N_V{\epsilon }_2{\epsilon }_5{I}_V^{\ast }+N_H\left(I_V^{\ast }-1\right)\left(S_{H_{i_a}}^{\ast }-S_{H_{i_a}V}^{\ast }\right){\beta }_V\right){\beta }_{H_{i_a}}+\left({\epsilon }_2{\epsilon }_6+{\epsilon }_5{\epsilon }_{11}{I}_V^{\ast }\right){\gamma }_{i_a}\\ a_4=d_H^3{\epsilon }_2+d_H{\epsilon }_2\left(d_H+{\epsilon }_{11}{I}_V^{\ast }\right)\left(d_H-{\epsilon }_5{\epsilon }_{11}{I}_V^{\ast }\right)-{\epsilon }_2{\epsilon }_{11}^2{\epsilon }_5{I}_V^{{\ast }^2}{\gamma }_{i_a}+\\ d_H^2\left(-N_V\left({\epsilon }_2{\epsilon }_6{I}_V^{\ast }-N_H\left(I_V^{\ast }-1\right)\left(S_{H_{i_a}}^{\ast }-S_{H_{i_a}V}^{\ast }{\epsilon }_5\right){\beta }_V\right){\beta }_{H_{i_a}}+{\epsilon }_2{\gamma }_{i_a}\right)-\\ d_H{\epsilon }_{11}{I}_V^{\ast }\left(N_V{\epsilon }_5\left({\epsilon }_2{I}_V^{\ast }+N_H\left(I_V^{\ast }-1\right)\left(S_{H_{i_a}}^{\ast }-S_{H_{i_a}V}^{\ast }\right){\beta }_V\right){\beta }_{H_{i_a}}+{\epsilon }_2{\epsilon }_6{\gamma }_{i_a}\right)\end{array}$$

where $S_{{}_{H_{i_a}}}^{\ast },S_{{}_{H_{i_a}V}}^{\ast },I_{{}_V}^{\ast },{\epsilon }_2,{\epsilon }_5,{\epsilon }_6$, and ${\epsilon }_{11}$ are defined in Equations (19) and (20). □

The eigenvalue of the endemic steady state will have negative real parts when the coefficients of Equation (35) satisfy the Routh–Hurwitz criteria [10,18]:

$$\begin{array}{l}{a}_1>0,\\ a_3>0,\\ a_4>0,\\ a_1{a}_2{a}_3-a_3^2-a_1^2{a}_4>0.\end{array}$$

(36)

To ensure that the polynomial of Equation (35) is Hurwitz, a numerical computation of the criteria in Equation (36) is carried out. Results are shown in Figure 4, which clearly illustrates that the polynomial of Equation (35) is indeed Hurwitz, and thus the endemic steady state $E_2$ is locally asymptotically stable. □

The values of parameters used are given in Table 3. From the above figures, the Routh–Hurwitz criteria conditions are satisfied for $R_0>1$.

We now proceed to investigate the global stabilities of both fixed points. In this respect, two theorems are given which illustrate the global stabilities.

Theorem 1.

Let $E_1^{\ast }=\left(S_{H_{i_a}}^{\ast },I_{H_{i_a}}^{\ast },S_{H_{i_a}V}^{\ast },I_{H_{i_a}V}^{\ast },I_V^{\ast }\right)=\left(\frac{b_H\left(1-P_{i_a}\right)}{d_H},0,\frac{b_H{P}_{i_a}}{d_H},0,0\right)$. If$R_0<1$, then the disease-free$E_1^{\ast }$is globally asymptotically stable on$\Omega$for all$i_a=1,2,3$where

$$\Omega =\left\{\left(S_{H_{i_a}},I_{H_{i_a}},S_{H_{i_a}V},I_{H_{i_a}V},I_V\right)\in {ℝ}_{+}^5,N_H\le \frac{b_H}{d_H+d_d},N_V\le \frac{A}{d_V+d_k}\right\}$$

Proof. 

We consider a Lyapunov function

$$\begin{array}{l}\theta \left(S_{H_{i_a}}^{\ast },I_{H_{i_a}}^{\ast },S_{H_{i_a}V}^{\ast },I_{H_{i_a}}^{\ast },I_V^{\ast }\right)={\displaystyle \sum _{i_a=1}^3\left(S_{H_{i_a}}-S_{H_{i_a}}^{\ast }\ln S_{H_{i_a}}\right)}+{\displaystyle \sum _{i_a=1}^3{I}_{H_{i_a}}}+{\displaystyle \sum _{i_a=1}^3{R}_{H_{i_a}}}\\ +{\displaystyle \sum _{i_a=1}^3\left(S_{H_{i_a}V}-S_{H_{i_a}V}^{\ast }\ln S_{H_{i_a}V}\right)}+{\displaystyle \sum _{i_a=1}^3{I}_{H_{i_a}V}}+I_V\end{array}$$

(37)

Then we have:

$$\theta \left(S_{H_{i_a}}^{\ast },I_{H_{i_a}}^{\ast },S_{H_{i_a}V}^{\ast },I_{H_{i_a}}^{\ast },I_V^{\ast }\right)=\left(\frac{b_H\left(1-P_{i_a}\right)}{d_H},0,\frac{b_H{P}_{i_a}}{d_H},0,0\right)$$

(38)

And the derivative with respect to time yields:

$$\begin{array}{l}\frac{d\theta }{dt}={\displaystyle \sum _{i_a=1}^3\frac{d{S}_{H_{i_a}}}{dt}\left(1-\frac{S_{H_{i_a}}^{\ast }}{S_{H_{i_a}}}\right)}+{\displaystyle \sum _{i_a=1}^3\frac{d{I}_{H_{i_a}}}{dt}}+{\displaystyle \sum _{i_a=1}^3\frac{d{R}_{H_{i_a}}}{dt}}\\ +{\displaystyle \sum _{i_a=1}^3\frac{d{S}_{H_{i_a}V}}{dt}\left(1-\frac{S_{H_{i_a}V}^{\ast }}{S_{H_{i_a}V}}\right)}+{\displaystyle \sum _{i_a=1}^3\frac{d{I}_{H_{i_a}V}}{dt}}+\frac{d{I}_V}{dt}\end{array}$$

(39)

$$\begin{array}{l}\frac{d\theta }{dt}={\displaystyle \sum _{i_a=1}^3\left(\left(\left(1-P_{i_a}\right)b_H-{\beta }_{H_{i_a}}{S}_{H_{i_a}}{I}_V{N}_V-d_H{S}_{H_{i_a}}\right)\left(1-\frac{S_{H_{i_a}}^{\ast }}{S_{H_{i_a}}}\right)\right)}\\ +{\displaystyle \sum _{i_a=1}^3\left(\left({\beta }_{H_{i_a}}{S}_{H_{i_a}}{I}_V{N}_V-d_H{S}_{H_{i_a}}\right)-{\gamma }_{i_a}{I}_{H_{i_a}}-\left(d_H+d_d\right)I_{H_{i_a}}\right)}\\ +{\displaystyle \sum _{i_a=1}^3\left({\gamma }_{i_a}{I}_{H_{i_a}}+{\gamma }_{i_a}{I}_{H_{i_a}V}-d_H{R}_{H_{i_a}}\right)}\\ +{\displaystyle \sum _{i_a=1}^3\left(\left(P_{i_a}{b}_H-\left(1-{\alpha }_{i_a}\right){\beta }_{H_{i_a}}{S}_{H_{i_a}V}{I}_V{N}_V-d_H{S}_{H_{i_a}V}\right)\left(1-\frac{S_{H_{i_a}V}^{\ast }}{S_{H_{i_a}V}}\right)\right)}\\ +{\displaystyle \sum _{i_a=1}^3\left(\left(1-{\alpha }_{i_a}\right){\beta }_{H_{i_a}}{S}_{H_{i_a}V}{I}_V{N}_V-{\gamma }_{i_a}{I}_{H_{i_a}V}-\left(d_H+d_d\right)I_{H_{i_a}V}\right)}\\ +{\beta }_V{S}_V{N}_H\left({\displaystyle \sum _{i_a=1}^3{I}_{H_{i_a}}+}{\displaystyle \sum _{i_a=1}^3{I}_{H_{i_a}V}}\right)-\left(d_V+d_k\right)I_V\end{array}$$

Simplification yields:

$$\begin{array}{l}\frac{d\theta }{dt}={\displaystyle \sum _{i_a=1}^3\left(\left(1-P_{i_a}\right)b_H-d_H{S}_{H_{i_a}}-\left(1-P_{i_a}\right)b_H\frac{S_{H_{i_a}}^{\ast }}{S_{H_{i_a}}}+d_H{S}_{H_{i_a}}\frac{S_{H_{i_a}}^{\ast }}{S_{H_{i_a}}}\right)}\\ -{\displaystyle \sum _{i_a=1}^3\left(d_H+d_d\right)I_{H_{i_a}}}-{\displaystyle \sum _{i_a=1}^3{d}_H{R}_{H_{i_a}}}+{\displaystyle \sum _{i_a=1}^3\left(P_{i_a}{b}_H-d_H{S}_{H_{i_a}V}-P_{i_a}{b}_H\frac{S_{H_{i_a}V}^{\ast }}{S_{H_{i_a}V}}+d_H{S}_{H_{i_a}V}\frac{S_{H_{i_a}V}^{\ast }}{S_{H_{i_a}V}}\right)}\\ -{\displaystyle \sum _{i_a=1}^3\left(d_H+d_d\right)I_{H_{i_a}V}}-\left(d_V+d_k\right)I_V\end{array}$$

(40)

$$\begin{array}{l}\frac{d\theta }{dt}={\displaystyle \sum _{i_a=1}^3\left(\left(1-P_{i_a}\right)b_H\left(\left(1-\frac{S_{H_{i_a}}^{\ast }}{S_{H_{i_a}}}\right)+\left(1-\frac{S_{H_{i_a}}}{S_{H_{i_a}}^{\ast }}\right)\right)\right)}-{\displaystyle \sum _{i_a=1}^3\left(d_H+d_d\right)I_{H_{i_a}}}-{\displaystyle \sum _{i_a=1}^3{d}_H{R}_{H_{i_a}}}\\ +{\displaystyle \sum _{i_a=1}^3\left(P_{i_a}{b}_H\left(\left(1-\frac{S_{H_{i_a}V}^{\ast }}{S_{H_{i_a}V}}\right)+\left(1-\frac{S_{H_{i_a}V}}{S_{H_{i_a}V}^{\ast }}\right)\right)\right)}-{\displaystyle \sum _{i_a=1}^3\left(d_H+d_d\right)I_{H_{i_a}V}}-\left(d_V+d_k\right)I_V\\ \frac{d\theta }{dt}=-{\displaystyle \sum _{i_a=1}^3\left(1-P_{i_a}\right)b_H\frac{{\left(S_{H_{i_a}}^{\ast }-S_{H_{i_a}}\right)}^2}{S_{H_{i_a}}^{\ast }{S}_{H_{i_a}}}}-{\displaystyle \sum _{i_a=1}^3\left(d_H+d_d\right)I_{H_{i_a}}}-{\displaystyle \sum _{i_a=1}^3{d}_H{R}_{H_{i_a}}}\\ -{\displaystyle \sum _{i_a=1}^3{P}_{i_a}{b}_H\frac{{\left(S_{H_{i_a}V}^{\ast }-S_{H_{i_a}V}\right)}^2}{S_{H_{i_a}V}^{\ast }{S}_{H_{i_a}V}}}-{\displaystyle \sum _{i_a=1}^3\left(d_H+d_d\right)I_{H_{i_a}V}}-\left(d_V+d_k\right)I_V\end{array}$$

$$\frac{d\theta }{dt}=-\left(\begin{array}{l}{\displaystyle \sum _{i_a=1}^3\left(1-P_{i_a}\right)b_H\frac{{\left(S_{H_{i_a}}^{\ast }-S_{H_{i_a}}\right)}^2}{S_{H_{i_a}}^{\ast }{S}_{H_{i_a}}}}+{\displaystyle \sum _{i_a=1}^3\left(d_H+d_d\right)I_{H_{i_a}}}+{\displaystyle \sum _{i_a=1}^3{d}_H{R}_{H_{i_a}}}\\ +{\displaystyle \sum _{i_a=1}^3{P}_{i_a}{b}_H\frac{{\left(S_{H_{i_a}V}^{\ast }-S_{H_{i_a}V}\right)}^2}{S_{H_{i_a}V}^{\ast }{S}_{H_{i_a}V}}}+{\displaystyle \sum _{i_a=1}^3\left(d_H+d_d\right)I_{H_{i_a}V}}+\left(d_V+d_k\right)I_V\end{array}\right)$$

(41)

It is obvious that all the terms appearing in Equation (41) are always nonpositive. Now, using LaSalle’s extension to Lyapunov’s theorem [21,25,26,27,28,29], we have $\frac{d\theta }{dt}\le 0$ and so the function $\frac{d\theta }{dt}$ is to be negative definite. The limit set of each solution is contained in the largest invariant set for which $S_{H_{i_a}}=S_{H_{i_a}}^{\ast },S_{H_{i_a}V}=S_{H_{i_a}V}^{\ast }$ and $R_{H_{i_a}}=0$ which is the singleton $\left\{E_1{}^{\ast }\right\}$. Then LaSalle’s invariant principle implies that the disease-free steady state $E_1{}^{\ast }$ is globally asymptotically stable on $\Omega$. □

Next, we consider the global property of the endemic steady state of Equations (13)–(17).

Theorem 2.

Let $E_2^{\ast }=\left(S_{H_{i_a}}^{\ast },I_{H_{i_a}}^{\ast },S_{H_{i_a}V}^{\ast },I_{H_{i_a}V}^{\ast },I_V^{\ast }\right)$given by Equation (19) for all$i_a=1,2,3$. If$R_0>1$, then the endemic$E_2^{\ast }$is globally asymptotically stable on$\Omega$if and only if

$$\{\begin{array}{l}{d}_H={\beta }_V{S}_V{N}_H-d_d\\ d_V={\displaystyle \sum _{i_a=1}^3{\beta }_{H_{i_a}}{N}_V{S}_{H_{i_a}}^{\ast }}\\ d_k={\displaystyle \sum _{i_a=1}^3\left(1-{\alpha }_{i_a}\right){\beta }_{H_{i_a}}{N}_V{S}_{H_{i_a}V}^{\ast }}\end{array}$$

(42)

Proof. 

The Lyapunov function is in the form:

$$\begin{array}{l}\varphi \left(S_{H_{i_a}}^{\ast },I_{H_{i_a}}^{\ast },S_{H_{i_a}V}^{\ast },I_{H_{i_a}V}^{\ast },I_V^{\ast }\right)={\displaystyle \sum _{i_a=1}^3\left(S_{H_{i_a}}-S_{H_{i_a}}^{\ast }\ln S_{H_{i_a}}\right)}+{\displaystyle \sum _{i_a=1}^3{I}_{H_{i_a}}}\\ +{\displaystyle \sum _{i_a=1}^3\left(S_{H_{i_a}V}-S_{H_{i_a}V}^{\ast }\ln S_{H_{i_a}V}\right)}+{\displaystyle \sum _{i_a=1}^3{I}_{H_{i_a}V}}+I_V\end{array}$$

(43)

For all $i_a=1,2,3.$

Then the derivative with respect to time yields:

$$\frac{d\varphi }{dt}={\displaystyle \sum _{i_a=1}^3\frac{d{S}_{H_{i_a}}}{dt}\left(1-\frac{S_{H_{i_a}}^{\ast }}{S_{H_{i_a}}}\right)}+{\displaystyle \sum _{i_a=1}^3\frac{d{I}_{H_{i_a}}}{dt}}+{\displaystyle \sum _{i_a=1}^3\frac{d{S}_{H_{i_a}V}}{dt}\left(1-\frac{S_{H_{i_a}V}^{\ast }}{S_{H_{i_a}V}}\right)}+{\displaystyle \sum _{i_a=1}^3\frac{d{I}_{H_{i_a}V}}{dt}}+\frac{d{I}_V}{dt}$$

(44)

$$\begin{array}{l}\frac{d\varphi }{dt}={\displaystyle \sum _{i_a=1}^3\left(\begin{array}{l}\left(1-P_{i_a}\right)b_H-{\beta }_{H_{i_a}}{S}_{H_{i_a}}{I}_V{N}_V-d_H{S}_{H_{i_a}}-\\ \left(1-P_{i_a}\right)b_H\frac{S_{H_{i_a}}^{\ast }}{S_{H_{i_a}}}+{\beta }_{H_{i_a}}{S}_{H_{i_a}}{I}_V{N}_V\frac{S_{H_{i_a}}^{\ast }}{S_{H_{i_a}}}+d_H{S}_{H_{i_a}}\frac{S_{H_{i_a}}^{\ast }}{S_{H_{i_a}}}\end{array}\right)}\\ +{\displaystyle \sum _{i_a=1}^3\left(\left({\beta }_{H_{i_a}}{S}_{H_{i_a}}{I}_V{N}_V-d_H{S}_{H_{i_a}}\right)-{\gamma }_{i_a}{I}_{H_{i_a}}-\left(d_H+d_d\right)I_{H_{i_a}}\right)}\\ +{\displaystyle \sum _{i_a=1}^3\left(\begin{array}{l}{P}_{i_a}{b}_H-\left(1-{\alpha }_{i_a}\right){\beta }_{H_{i_a}}{S}_{H_{i_a}V}{I}_V{N}_V-d_H{S}_{H_{i_a}V}\\ -P_{i_a}{b}_H\frac{S_{H_{i_a}V}^{\ast }}{S_{H_{i_a}V}}+\left(1-{\alpha }_{i_a}\right){\beta }_{H_{i_a}}{S}_{H_{i_a}V}{I}_V{N}_V\frac{S_{H_{i_a}V}^{\ast }}{S_{H_{i_a}V}}+d_H{S}_{H_{i_a}V}\frac{S_{H_{i_a}V}^{\ast }}{S_{H_{i_a}V}}\end{array}\right)}\\ +{\displaystyle \sum _{i_a=1}^3\left(\left(1-{\alpha }_{i_a}\right){\beta }_{H_{i_a}}{S}_{H_{i_a}V}{I}_V{N}_V-{\gamma }_{i_a}{I}_{H_{i_a}V}-\left(d_H+d_d\right)I_{H_{i_a}V}\right)}\\ +{\beta }_V{S}_V{N}_H\left({\displaystyle \sum _{i_a=1}^3{I}_{H_i}+}{\displaystyle \sum _{i_a=1}^3{I}_{H_{i_a}V}}\right)-\left(d_V+d_k\right)I_V\end{array}$$

$$\begin{array}{l}\frac{d\varphi }{dt}={\displaystyle \sum _{i_a=1}^3\left(\begin{array}{l}\left(1-P_{i_a}\right)b_H-d_H{S}_{H_{i_a}}-\left(1-P_{i_a}\right)b_H\frac{S_{H_{i_a}}^{\ast }}{S_{H_{i_a}}}\\ +{\beta }_{H_{i_a}}{S}_{H_{i_a}}{I}_V{N}_V\frac{S_{H_{i_a}}^{\ast }}{S_{H_{i_a}}}+d_H{S}_{H_{i_a}}\frac{S_{H_{i_a}}^{\ast }}{S_{H_{i_a}}}\end{array}\right)}\\ +{\displaystyle \sum _{i_a=1}^3\left(\left(-d_H{S}_{H_{i_a}}\right)-{\gamma }_{i_a}{I}_{H_{i_a}}-\left(d_H+d_d\right)I_{H_{i_a}}\right)}\\ +{\displaystyle \sum _{i_a=1}^3\left(\begin{array}{l}{P}_{i_a}{b}_H-d_H{S}_{H_{i_a}V}-P_{i_a}{b}_H\frac{S_{H_{i_a}V}^{\ast }}{S_{H_{i_a}V}}\\ +\left(1-{\alpha }_{i_a}\right){\beta }_{H_{i_a}}{S}_{H_{i_a}V}{I}_V{N}_V\frac{S_{H_{i_a}V}^{\ast }}{S_{H_{i_a}V}}+d_H{S}_{H_{i_{a}V}}\frac{S_{H_{i_a}V}^{\ast }}{S_{H_{i_a}V}}\end{array}\right)}\\ +{\displaystyle \sum _{i_a=1}^3\left(-{\gamma }_{i_a}{I}_{H_{i_a}V}-\left(d_H+d_d\right)I_{H_{i_a}V}\right)}\\ +{\beta }_V{S}_V{N}_H\left({\displaystyle \sum _{i_a=1}^3{I}_{H_{i_a}}+}{\displaystyle \sum _{i_a=1}^3{I}_{H_{i_a}V}}\right)-\left(d_V+d_k\right)I_V\end{array}$$

(45)

Substituting the condition set out in Equation (42) into Equation (45), we obtain:

$$\begin{array}{l}\frac{d\varphi }{dt}={\displaystyle \sum _{i_a=1}^3\left(\left(1-P_{i_a}\right)b_H\left(\left(1-\frac{S_{H_{i_a}}^{\ast }}{S_{H_{i_a}}}\right)+\left(1-\frac{S_{H_{i_a}}}{S_{H_{i_a}}^{\ast }}\right)\right)\right)}-{\displaystyle \sum _{i_a=1}^3{\gamma }_{i_a}{I}_{H_{i_a}}}\\ +{\displaystyle \sum _{i_a=1}^3\left(P_{i_a}{b}_H\left(\left(1-\frac{S_{H_{i_a}V}^{\ast }}{S_{H_{i_a}V}}\right)+\left(1-\frac{S_{H_{i_a}V}}{S_{H_{i_a}V}^{\ast }}\right)\right)\right)}-{\displaystyle \sum _{i_a=1}^3{\gamma }_{i_a}{I}_{H_{i_a}V}}\end{array}$$

$$\frac{d\varphi }{dt}=-{\displaystyle \sum _{i_a=1}^3\left(1-P_{i_a}\right)b_H\frac{{\left(S_{H_{i_a}}-S_{H_{i_a}}^{\ast }\right)}^2}{S_{H_{i_a}}{S}_{H_{i_a}}^{\ast }}}-{\displaystyle \sum _{i_a=1}^3{\gamma }_{i_a}{I}_{H_{i_a}}}-{\displaystyle \sum _{i_a=1}^3{P}_{i_a}{b}_H\frac{{\left(S_{H_{i_a}V}-S_{H_{i_a}V}^{\ast }\right)}^2}{S_{H_{i_a}V}{S}_{H_{i_a}V}^{\ast }}}-{\displaystyle \sum _{i_a=1}^3{\gamma }_{i_a}{I}_{H_{i_a}V}}$$

$$\frac{d\varphi }{dt}=-\left({\displaystyle \sum _{i_a=1}^3\left(1-P_{i_a}\right)b_H\frac{{\left(S_{H_{i_a}}-S_{H_{i_a}}^{\ast }\right)}^2}{S_{H_{i_a}}{S}_{H_{i_a}}^{\ast }}}+{\displaystyle \sum _{i_a=1}^3{\gamma }_{i_a}{I}_{H_{i_a}}}+{\displaystyle \sum _{i_a=1}^3{P}_{i_a}{b}_H\frac{{\left(S_{H_{i_a}V}-S_{H_{i_a}V}^{\ast }\right)}^2}{S_{H_{i_a}V}{S}_{H_{i_a}V}^{\ast }}}+{\displaystyle \sum _{i_a=1}^3{\gamma }_{i_a}{I}_{H_{i_a}V}}\right)$$

(46)

Hence, the derivative $\frac{d\varphi }{dt}\le 0$ for all ${\left(S_{H_{i_a}},I_{H_{i_a}},S_{H_{i_a}V},I_{H_{i_a}V},I_V\right)}_{\forall i_a=1,2,3}\in \Omega$ with $\frac{d\varphi }{dt}=0$ if and only if $S_{H_{i_a}}=S_{H_{i_a}}^{\ast },S_{H_{i_a}V}=S_{H_{i_a}V}^{\ast },I_{H_{i_a}}=0,$ and $I_{H_{i_a}V}=0$. Therefore, by LaSalle’s extension to Lyapunov’s function theorem [21,25,26,27,28,29], the endemic steady state $E_2{}^{\ast }$ is globally asymptotically stable on $\Omega$. □

2.2.3. Bifurcation Analysis

The previous section investigates the local and global stabilities and established that the value of $R_0=1$ is indeed a bifurcation value since the disease-free equilibrium changes its local and global stability properties there. Now, to further investigate the exact kind of bifurcation the unity basic reproductive number contains, we will make use of the results below. Let us consider a general system of ODEs with a parameter $\xi$:

$$\stackrel{•}{{\displaystyle x}}=f(x,\xi );f:R^n\times R\to R^n,f\in C^2(R^n\times R)$$

(47)

Without loss of generally, we assume that $x=0$ is a steady state for Equation (47).

Theorem 3.

(See (Castillo-Chavez and Song [30])). Assume:

$\left(A1\right)$$A=D_{x}f(0,0)$is the linearization matrix of system (47) around the steady state$x=0$with$\xi$evaluated at 0. Zero is a simple eigenvalue of$A$and all other eigenvalues of$A$have negative real parts.

$\left(A2\right)$Matrix$A$has a right eigenvector$\omega$and a left eigenvector$\theta$corresponding to the zero eigenvalue. Let$f_k$denotes the$kth$component of$f$and,

$$a={\displaystyle \sum _{k,i,j}^n{\theta }_k{\omega }_i{\omega }_j\frac{{\partial }^2{f}_k}{\partial x_i\partial x_j}}(0,0),b={\displaystyle \sum _{k,i,j}^n{\theta }_k{\omega }_i\frac{{\partial }^2{f}_k}{\partial x_i\partial \xi }}(0,0).$$

(48)

Then the local dynamics of system Equation (47) around$x=0$are totally determined by$a$and$b$. The conditions and their meanings are:

$\left(i\right)$ $a>0,b>0$when$\xi <0,$with$\left|\xi \right|\ll 1,x=0$is locally asymptotically stable and there exists a positive unstable steady states; when$0<\left|\xi \right|\ll 1,x=0$is unstable and there exists a negative and locally asymptotically stable steady states;

$\left(ii\right)$ $a<0,b<0$when$\xi <0,$with$\left|\xi \right|\ll 1,x=0$is unstable; when$0<\left|\xi \right|\ll 1,x=0$is locally asymptotically stable and there exists a positive unstable steady states;

$\left(iii\right)$ $a>0,b<0$when$\xi <0,$with$\left|\xi \right|\ll 1,x=0$is unstable and there exists a locally asymptotically stable negative steady states; when$0<\left|\xi \right|\ll 1,x=0$is stable and a positive unstable steady states appears;

$\left(iv\right)$ $a<0,b>0$when$\xi$changes from negative to positive, $x=0$changes its stability from stable to unstable. Correspondingly, a negative unstable steady states becomes positive and locally asymptotically stable.

At$\xi =0$, a transcritical bifurcation takes place, that is, if$a<0$and$b>0$this bifurcation is forward, and if$a>0$and$b>0$the transcritical bifurcation is backward.

Proof. 

See Castillo-Chavez and Song [30]. □

Theorem 4.

At $R_0=1$direction of the bifurcation of system Equations (13)–(17) exhibits a backward bifurcation if$\mu >1$. If$\mu <1$, direction of the bifurcation of system Equations (13)–(17) exhibits a forward bifurcation when$R_0=1$.

Proof. 

We consider the disease-free steady states $E_1^{\ast }\left(t\right)=\left(\frac{b_H\left(1-P_{i_a}\right)}{d_H},0,\frac{b_H{P}_{i_a}}{d_H},0,0\right)$ for all $i_a=1,2,3$ and observe that the condition $R_0=1$ which gives ${\beta }_{H_{i_a}}={\beta }^{\ast }$ as

$${\beta }^{\ast }=-\left(\frac{d_H{d}_d{d}_k+d_H^2{d}_k+d_H{d}_d{d}_V+d_H^2{d}_V+d_H{d}_k{\gamma }_{i_a}-d_H{d}_V{\gamma }_{i_a}}{b_H{N}_H{N}_V(P_{i_a}{\alpha }_{i_a}-1){\beta }_V}\right)$$

Evaluating the Jacobian matrix at ${\beta }^{\ast }:$

$$J\left(E_1,{\beta }^{\ast }\right)=\left[\begin{array}{ccccc}-d_H& 0& 0& 0& -Q_1\\ 0& Q_4& 0& 0& Q_1\\ 0& 0& -d_H& 0& Q_2\\ 0& 0& 0& Q_4& Q_3\\ 0& {\beta }_V{N}_H& 0& {\beta }_V{N}_H& -d_V-d_k\end{array}\right]$$

(49)

where

$$Q_1=\frac{b_H{N}_V\left(1-P_{i_a}\right){\beta }^{\ast }}{d_H},Q_2=\frac{b_H{N}_V{P}_{i_a}\left({\alpha }_{i_a}-1\right){\beta }^{\ast }}{d_H},Q_3=\frac{b_H{N}_V{P}_{i_a}\left(1-{\alpha }_{i_a}\right){\beta }^{\ast }}{d_H},\mathrm{and}{Q}_4=-d_H-d_d-{\gamma }_{i_a}.$$

The eigenvalues of $J\left(E_1,{\beta }^{\ast }\right)$ are obtained again by solving for the zeros of the characteristic polynomial of $J\left(E_1,{\beta }^{\ast }\right)$, yielding: ${\lambda }_1={\lambda }_2=-d_H,{\lambda }_3=-d_H-d_d-{\gamma }_{i_a},{\lambda }_4=-d_H-d_d-{\gamma }_{i_a}-d_V-d_k,$ and ${\lambda }_5=0$.

Thus ${\lambda }_5=0$ is a simple zero eigenvalue of the matrix $J(E_1,{\beta }^{\ast })$. All other eigenvalues have negative real parts. Hence, when ${\beta }_{H_{i_a}}={\beta }^{\ast }$ (when $R_0=1$), the disease-free steady state $E_1$ is a non-hyperbolic steady state: the assumption $\left(A1\right)$ of Theorem 3 is satisfied. Representing the right eigenvector $\omega ={({\omega }_1,{\omega }_2,{\omega }_3,{\omega }_4,{\omega }_5)}^T$ with the associated eigenvalue ${\lambda }_5=0$. The following system of the equation then follows:

$$-d_H{\omega }_1-\frac{b_H{N}_V(1-P_{i_a}){\beta }^{\ast }}{d_H}{\omega }_5=0$$

(50)

$$-(d_H+d_d+{\gamma }_{i_a}){\omega }_2+\frac{b_H{N}_V(1-P_{i_a}){\beta }^{\ast }}{d_H}{\omega }_5=0$$

(51)

$$-d_H{\omega }_3+\frac{b_H{N}_V{P}_{i_a}({\alpha }_{i_a}-1){\beta }^{\ast }}{d_H}{\omega }_5=0$$

(52)

$$-(d_H+d_d+{\gamma }_{i_a}){\omega }_4+\frac{b_H{N}_V{P}_{i_a}(1-{\alpha }_{i_a}){\beta }^{\ast }}{d_H}{\omega }_5=0$$

(53)

$$N_H{\beta }_V{\omega }_2+N_H{\beta }_V{\omega }_4-(d_V+d_k){\omega }_5=0$$

(54)

Solving for the right eigenvector, we obtain:

$$\omega ={\left(-\frac{d_d+d_H+{\gamma }_{i_a}}{d_H},1,{\omega }_3,-\frac{P_{i_a}-P_{i_a}{\alpha }_{i_a}}{P_{i_a}-1},{\omega }_5\right)}^T$$

(55)

where

$${\omega }_3=-\frac{P_{i_a}({\alpha }_{i_a}-1)(d_d+d_H+{\gamma }_{i_a})}{d_H},\mathrm{and}{\omega }_5=\frac{N_H{\beta }_V(P_{i_a}{\alpha }_{i_a}-1)}{(d_k+d_V)(P_{i_a}-1)}.$$

Likewise, the system of equations for the left eigenvector $\theta =({\theta }_1,{\theta }_2,{\theta }_3,{\theta }_4,{\theta }_5)$ satisfying condition $\omega \cdot \theta =1$ is given by:

$$-d_H{\theta }_1=0$$

(56)

$$-(d_H+d_d+{\gamma }_{i_a}){\theta }_2+N_H{\beta }_V{\theta }_5=0$$

(57)

$$-d_H{\theta }_3=0$$

(58)

$$-(d_H+d_d+{\gamma }_{i_a}){\theta }_4+N_H{\beta }_V{\theta }_5=0$$

(59)

$$\left(\begin{array}{l}-\frac{b_H{N}_V(1-P_{i_a}){\beta }^{\ast }}{d_H}{\theta }_1+\frac{b_H{N}_V(1-P_{i_a}){\beta }^{\ast }}{d_H}{\theta }_2+\frac{b_H{N}_V{P}_{i_a}({\alpha }_{i_a}-1){\beta }^{\ast }}{d_H}{\theta }_3+\\ \frac{b_H{N}_V{P}_{i_a}(1-{\alpha }_{i_a}){\beta }^{\ast }}{d_H}{\theta }_4-(d_k-d_V){\theta }_5\end{array}\right)=0$$

(60)

The left eigenvector $\theta$ then turns out to be:

$$\theta =\left(0,{\theta }_2,0,{\theta }_4,{\theta }_5\right)$$

(61)

where

$${\theta }_2={\theta }_4=\frac{(d_k+d_V)(P_{i_a}-1)}{(P_{i_a}{\alpha }_{i_a}-1)(d_d+d_H+d_k+d_V+{\gamma }_{i_a})},$$

and

$${\theta }_5=\frac{(d_d+d_H+{\gamma }_{i_a})(d_k+d_V)(P_{i_a}-1)}{N_H{\beta }_V(P_{i_a}{\alpha }_{i_a}-1)(d_d+d_H+d_k+d_V+{\gamma }_{i_a})}.$$

Now we can thus compute the coefficient $a$ and $b$ given in Theorem 3, that is

$$a={\displaystyle \sum _{k,i,j=1}^5{\theta }_k{\omega }_i{\omega }_j\frac{{\partial }^2{f}_k}{\partial x_i\partial x_j}}(E_1,{\beta }^{\ast })$$

(62)

$$b={\displaystyle \sum _{k,i=1}^5{\theta }_k{\omega }_i\frac{{\partial }^2{f}_k}{\partial x_i\partial \xi }}(E_1,{\beta }^{\ast })$$

(63)

Taking into account the system Equations (13)–(17) and considering only the constituent terms in $a$ and $b$ for which the derivatives $\frac{{\partial }^2{f}_k}{\partial x_i\partial x_j}(E_1,{\beta }^{\ast })$ and $\frac{{\partial }^2{f}_k}{\partial x_i\partial \xi }(E_1,{\beta }^{\ast })$ are nonzero. Therefore, we have:

$$\begin{array}{l}a=2{\theta }_1{\omega }_1{\omega }_5\frac{{\partial }^2{f}_1}{\partial S_{H_{i_a}}\partial I_V}(E_1,{\beta }^{\ast })+2{\theta }_2{\omega }_1{\omega }_5\frac{{\partial }^2{f}_2}{\partial S_{H_{i_a}}\partial I_V}(E_1,{\beta }^{\ast })\\ +2{\theta }_3{\omega }_3{\omega }_5\frac{{\partial }^2{f}_3}{\partial S_{H_{i_a}V}\partial I_V}(E_1,{\beta }^{\ast })+2{\theta }_4{\omega }_3{\omega }_5\frac{{\partial }^2{f}_4}{\partial S_{H_{i_a}V}\partial I_V}(E_1,{\beta }^{\ast }),\end{array}$$

(64)

$$\begin{array}{l}b={\theta }_1{\omega }_1\frac{{\partial }^2{f}_1}{\partial S_{H_{i_a}}\partial {\beta }_{H_{i_a}}}(E_1,{\beta }^{\ast })+{\theta }_1{\omega }_5\frac{{\partial }^2{f}_1}{\partial I_V\partial {\beta }_{H_{i_a}}}(E_1,{\beta }^{\ast })\\ +{\theta }_2{\omega }_1\frac{{\partial }^2{f}_2}{\partial S_{H_{i_a}}\partial {\beta }_{H_{i_a}}}(E_1,{\beta }^{\ast })+{\theta }_2{\omega }_5\frac{{\partial }^2{f}_2}{\partial I_V\partial {\beta }_{H_{i_a}}}(E_1,{\beta }^{\ast })\\ +{\theta }_3{\omega }_3\frac{{\partial }^2{f}_3}{\partial S_{H_{i_a}V}\partial {\beta }_{H_{i_a}}}(E_1,{\beta }^{\ast })+{\theta }_3{\omega }_5\frac{{\partial }^2{f}_3}{\partial I_V\partial {\beta }_{H_{i_a}}}(E_1,{\beta }^{\ast })\\ +{\theta }_4{\omega }_3\frac{{\partial }^2{f}_4}{\partial S_{H_{i_a}V}\partial {\beta }_{H_{i_a}}}(E_1,{\beta }^{\ast })+{\theta }_4{\omega }_5\frac{{\partial }^2{f}_4}{\partial I_V\partial {\beta }_{H_{i_a}}}(E_1,{\beta }^{\ast }).\end{array}$$

(65)

Computing all the required derivatives of $f_1=S_{H_{i_a}},f_2=I_{H_{i_a}},f_3=S_{H_{i_a}V},f_4=I_{H_{i_a}V},$ and $f_5=I_V$ from Equations (13)–(17), and substituting all the values of ${\omega }_i$ from Equation (55) and ${\theta }_i$ of Equation (61) into system Equation (64) and Equation (65) we obtain,

$$a=\frac{2N_H{N}_V{\beta }_{H_{i_a}}{\beta }_V(d_d+d_H+{\gamma }_{i_a})({\alpha }_{i_a}-2)}{d_H(d_k+d_V+d_d+d_H+{\gamma }_{i_a})},$$

(66)

$$b=\frac{b_H{\beta }_V(1-P_{i_a}{\alpha }_{i_a})N_H{N}_V}{d_H(d_k+d_V+d_d+d_H+{\gamma }_{i_a})}.$$

(67)

Let us introduce:

$$\mu =\frac{{\alpha }_{i_a}}{2}.$$

(68)

It is clear from Equation (67) that the coefficient $b$ is positive because $(1-P_{i_a}{\alpha }_{i_a})$ cannot be greater than or equal to one, its constant rate. The coefficient $a$ is positive if and only if $\mu >1$ in this case. Hence, the direction of the bifurcation of system Equations (13)–(17) is backward at $R_0=1$. Note, that if $\mu <1$ then it is clear that $a<0$ and $b>0$. The direction of bifurcation is thus forward. □

2.2.4. Sensitivity Analysis

In this section, we carried out a sensitivity analysis [31,32] test of the effect of the parameters on the basic reproduction number R₀ to ascertain the impact of these parameters on dengue transmission.

Definition 1

[33]. The sensitivity index of $R_0$, which depends differentially on a parameter $\mu$, is defined as:

$${\psi }_{\mu }^{R_0}=\frac{\partial R_0}{\partial \mu }\times \frac{\mu }{R_0}$$

(69)

The basic reproduction number in order to control the spread of dengue disease transmission model is given by Equation (32). The sensitivity index of Equation (69) can now be used to evaluate the index of each parameter. The sensitivity indices of the basic reproduction number ($R_0$) with respect to the parameters are shown in

Table 2. Sensitivity indices.

Parameters	Sensitivity Indices
$N_H$	0.50000
$N_V$	0.50000
${\beta }_V$	0.50000
$d_H$	0.00010
$d_d$	0.18643
${\alpha }_1$	−0.01736
${\alpha }_2$	−0.03510
${\alpha }_3$	−0.00660
${\beta }_{H_1}$	0.17551
${\beta }_{H_2}$	0.11116
${\beta }_{H_3}$	0.21333
${\gamma }_1$	0.26206
${\gamma }_2$	0.34021
${\gamma }_3$	0.21120

. The parameters that have positive indices that are $N_H,N_V,d_H,d_d,{\beta }_{H_{i_a}},{\gamma }_{i_a}$, and ${\beta }_V$ for all $i_a=1,2,3$ have a positive effect on the basic reproduction number. It means the increase in the number of the infected human population ($I_{H_{i_a}}$) and infected human population with vaccination ($I_{H_{i_a}V}$). Furthermore, the parameters in which their sensitivity indices are negative is ${\alpha }_{i_a}$ for all $i_a=1,2,3$ has a negative effect to minimize the endemic of the disease. From Table 2, we give examples of parameters ${\beta }_{H_1},{\beta }_{H_2},{\beta }_{H_3}$ and ${\beta }_V$ effects on all classes as shown in Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9.

3. Numerical Results

This section presents the numerical results of the transmission of dengue disease in this study based on the SIR model with vaccination by considering the human population of all three age groups for all $i_a=1,2,3$. The numerical simulations assume the following values of the parameters given in

Table 3. Parameter and their values used.

Parameters	Value	References
${\beta }_V$	0.00015	assumed
$b_H$	1/(365 ∗ 70)	[10,18,24,31,34]
$d_H$	1/(365 ∗ 70)	[10,18,24,31,34]
$d_V$	1/12	[10,18,24,31,34]
$d_d$	1/14	assumed
$d_k$	1/7	assumed
$N_H$	100,000	[10,18,24,31,34]
$N_V$	100,000	[10,18,24,31,34]

. Note, that the transmission rates quoted per day, correspond to a life expectancy in Thailand of approximately 70 years. It is also assumed that the birth and death rates are equal and constant.

The parameter used for simulating the disease-free steady state are:

$$\begin{array}{l}{P}_1=0.2,P_2=0.6,P_3=0.2,{\gamma }_1=0.3,{\gamma }_2=0.5,{\gamma }_3=0.2,\\ {\alpha }_1=0.3,{\alpha }_2=0.4,{\alpha }_3=0.3,{\beta }_{H_1}=0.00030,{\beta }_{H_2}=0.00035,{\beta }_{H_3}=0.00025.\end{array}$$

These parameters give a basic reproduction number of $R_0=0.0014$. The vaccine efficacy parameters used for the endemic case are:

$$P_1=0.2,P_2=0.6,P_3=0.2,{\gamma }_1=0.3,{\gamma }_2=0.5,{\gamma }_3=0.2,{\alpha }_1=0.3,{\alpha }_2=0.4,{\alpha }_3=0.3$$

Note, that reports on the vaccine efficacies against the virus in the participant of the second age group has a pooled estimate of 65.5%, as compared to 44.6% for participants in the first age group [14,16,17,34,35]. This thus prompts the assumption that vaccination would work best for the second age group. Consequently, it is assumed that the vaccine effectiveness parameters $P_2$ for the second age group are 0.6, while the corresponding effectiveness parameters for the first and third groups are set to be significantly lowered at $P_1=P_3=0.2$. In addition, the following transmission rate vectors are used:

$${\beta }_{H_1}=\left[0.35,0.40,0.45,0.50,0.55\right],$$

(70)

$${\beta }_{H_2}=\left[0.60,0.65,0.70,0.75,0.80\right],$$

(71)

$${\beta }_{H_3}=\left[0.10,0.15,0.20,0.25,0.30\right],$$

(72)

$${\beta }_V=\left[0.5,0.6,0.7,0.8,0.9\right].$$

(73)

These vaccine efficacy parameters, along with the transmission rate vectors of Equations (70)–(73), yield $R_0>1$.

The dynamical model, after normalizing variables, has 13 compartments namely, $S_{H_1},S_{H_2},S_{H_3},I_{H_1},I_{H_2},I_{H_3},S_{H_{1}V},S_{H_{2}V},S_{H_{3}V},I_{H_{1}V},I_{H_{2}V},I_{H_{3}V},$ and $I_V$. For the disease-free case, Figure 5a shows the trajectories of $S_{H_1},S_{H_2}$, and $S_{H_3}$; Figure 5b shows the trajectories of $I_{H_1},I_{H_2},I_{H_3}$; Figure 5c shows the trajectories of $S_{H_{1}V},S_{H_{2}V},S_{H_{3}V}$; Figure 5d shows the trajectories of $I_{H_{1}V},I_{H_{2}V},I_{H_{3}V}$, and Figure 5e shows the trajectories of $I_V$. For the endemic state shown in Figure 6, Figure 7, Figure 8 and Figure 9, Figure 6 shows the trajectories of age group 1, namely, $S_{H_1},I_{H_1},S_{H_{1}V}$, and $I_{H_{1}V}$. Figure 7 shows the trajectories of age group 2: $S_{H_2},I_{H_2},S_{H_{2}V}$, and $I_{H_{2}V}$. Figure 8 shows the trajectories of age group 3: $S_{H_3},I_{H_3},S_{H_{3}V}$, and $I_{H_{3}V}$. Moreover, lastly, Figure 9 shows the trajectories of $I_V$.

The trajectories of the numerical solutions for disease-free state and endemic state of $S_{H_1},S_{H_2},S_{H_3},I_{H_1},I_{H_2},I_{H_3},S_{H_{1}V},S_{H_{2}V},S_{H_{3}V},I_{H_{1}V},I_{H_{2}V},I_{H_{3}V},$ and $I_V$ are shown in Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9. Note, that the time responses in some graphs appear to be congruent on a straight line for the first 5 days of the trajectory before splitting thereafter, notably for Figure 6a,c, Figure 7a,c and Figure 8a,c. However, upon zooming into the responses in question, it is found that such congruence is caused by the overlapping of multiple lines. The overall trend seen with the responses of Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9 is that of exponential decay to zero. The increments of the transmission rate of the dengue virus from vector to human ${\beta }_{H_1},{\beta }_{H_2},$ and ${\beta }_{H_3}$ serve to reduce the number of susceptible people, both without the vaccine and also with the vaccine; while increasing the infected classes (both with and without the vaccine), as the transmission rate is increased. However, the trajectories of Figure 7 are closer together than those of Figure 6 and Figure 8, thereby suggesting that the sensitivity of ${\beta }_{H_2}$ to changes in the responses is less than those of ${\beta }_{H_1}$ and ${\beta }_{H_3}$, whose trajectory exhibit greater spacing between one another as the transmission rate is changed.

These phase portraits all support the theoretical calculations that, for $R_0<1$, the disease-free equilibrium is stable, while for $R_0>1$, the endemic equilibrium is stable.

4. Discussion and Conclusions

In this paper, we analyzed the SIR model for human populations with vaccination and SI for vector populations. Specifically, the human population is divided into three age structured groups as follows: the first group is 0–9 years with $i_a=1$, the second is 10–44 years $i_a=2$ and the last group is order than 44 years $i_a=3$. The analysis is based on the Routh–Hurwitz criteria to establish the local asymptotically stability of the endemic steady state. Two steady states that we found are disease-free and endemic steady states. We obtained the basic reproductive number $R_0$ using the Next-generation matrix. If $R_0<1$, the disease-free steady state $E_1$ is local asymptotically stable. If $R_0>1$, $E_1$ is unstable and the endemic steady state $E_2$ is local asymptotically stable. The global is asymptotically stable under conditions in Theorem 1 and Theorem 2 using the Lyapunov function method. Conditions for backward and forward Hopf bifurcations are given in Theorem 3 and Theorem 4 when $R_0=1$.

Numerical simulations were also conducted for a set of different parameter values to obtain the different trajectories for different population groups in the model. These trajectories are shown in Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9. Note, that the differences in the time scales exist between the plots of the susceptible and infected individuals for the case of $R_0>1$ as shown in Figure 6, Figure 7 and Figure 8. These differences are due to the fact that, with $R_0>1$, an infected individual can infect on average a large number of people. As a result, the susceptible population takes a short time to decay to zero, while the infected population takes a lot longer to decay to zero. The infection for all the trajectories of Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9 is to depict the progression of the outbreak from its origin to its end without truncations in time. Note, that in the case of the endemic steady state transmission rate of dengue virus from vector to human, ${\beta }_{H_{i_a}}$ for all $i_a=1,2,3$ and the transmission rate of dengue virus from human to vector, ${\beta }_V$, affects the convergence into stability. If the transmission rate is large enough, the convergence will be faster. Similarly, smaller transmission rates slow the convergence into stability.

The dynamics concept of this model is appropriate for high-endemic areas in Thailand and should be vaccinated to prevent dengue severity; from the data in Figure 2, it can be seen that the percentage of dengue infection in the age group 9–44 years is more than all ages This is consistent with the World Health Organization’s recommendation for vaccination in this age group. Numerical solutions suggest that when the vaccine is highly effective in this age group, outbreaks can be controlled. This mathematical model may be an important guideline for Thailand’s Ministry of Public Health to apply to areas with severe outbreaks, such as the Northeast.

Though the presented model takes into account the age structure of the population, it does not consider the serostatus before vaccination because there is no testing of blood before injection [17,34,35]; nor does it consider the multiple coexisting serotypes of the dengue virus. Future work will attempt to consider such factors, as well as the environment and topography affecting the outbreak. Given the new limitations on whom the vaccination should be directed at, the fact that most of the people contracting the disease are below the age of twenty will mean that the dengue occurrences for this age group will be more difficult to manage. Further research could uncover ways of managing such constraints optimally.