A mathematical model to study the 2014–2015 large-scale dengue epidemics in Kaohsiung and Tainan cities in Taiwan, China

1.   Introduction

Dengue virus (DENV), a flavivirus primarily transmitted to human via a bites of an infected mosquitoes of the genus Aedes (i.e., Aedes aegypti or Aedes albopictus) ^{[1,2,3,4,5,6,7,8,9,10]}. This species is capable of transmitting other mosquito-borne diseases such as Chikungunya, Zika, Yellow fever ^{[3,11,12,13,14]}. The DENV is endemic in at least 100 countries in the tropical and subtropical regions of the world, with nearly 4 billion people at risk ^[7,15]. It is second to malaria amongst the deadly mosquito-borne diseases with fast-spreading ability ^[7], with an estimate of 390 million new infections of which about 96 million are clinical, up to 1 million are severe cases (i.e., dengue hemorrhagic fever/dengue shock syndrome also known as DHF/DSS), and more than 20 thousand deaths every year globally ^[6,7,16]. The DHF/DSS have a high mortality risk with death occurring within 24hrs after the onset of shock ^[17,18,19]).

The DENV are of four serotypes (i.e., DENV-1, DENV-2, DENV-3 and DENV-4) with low cross-immunity among each other which can result in secondary infection after an infection with one serotype has happened ^[7]. An individual infected with one serotype can be infected with different serotypes about a half year later ^[3], yet there is no proof of reinfection with the same serotype. DEN infection of all serotypes can cause three distinct syndromes: classic dengue fever, DHF and DSS ^[3].

The epidemiological cycle occurs as follows. Female Aedes aegypti mosquitoes are infected by the blood meal from biting an infectious human during the viremic stage of the infection. After a period of time called the extrinsic incubation period (EIP), the mosquitoes become infectious ^[20]. The range of the extrinsic incubation period lies between 4 to 10 days ^[20]. The infected mosquito does not have immunity and can transmit the virus throughout its lifespan ^[3,7,20]. An infectious mosquito can transmit the virus to a susceptible individual via bites ^[21], after a period of time called the intrinsic incubation period (IIP), which varies from 4 to 7 days ^[20], the virus can evolve from classical dengue fever to the secondary cases of dengue infection which may lead to DHF or DSS ^[17]. The manifestation of the symptoms varies between 3 to 7 days ^[3], which corresponds roughly to the infectious period. Thereafter, the individual develops the life-long immunity against the same serotype of the virus ^[3].

The prevention against DENV is limited to controlling mechanisms applied on the mosquito, because the vaccine is not yet available (even though there is a vaccine called dengvaxia, which is still under clinical trial, released by Sanofi Pasteur in 2015 which was approved by morethan 10 countries ^[7,22]. The vaccine efficacy varies by the serotypes of the DENV (i.e., 54.7 % for DENV-1, 43.0 % for DENV-2, 71.6 % for DENV-3, and 76.9 % for DENV-4 ^[7,23]). To describe the overall transmission dynamics of the DENV, we developed a mathematical model that takes into account the female adult mosquito coupled with the human population to analyses and explain the transmission dynamics of the DENV epidemics in Taiwan in 2014–2015.

A number of mathematical modeling studies have been carried out to gain insight into the DENV transmission dynamics in human and/or mosquito populations ^{[3,4,5,6,11,24,25]}. Garba et al. ^[6] developed a deterministic model for the DENV transmission within human and mosquito populations, and extended the model to incorporate an imperfect vaccine against the virus, and find that both models exhibited backward bifurcation phenomenon, and the courses have been explained in detail by Gumel ^[26]. Yang and Ferreiri ^[3] proposed a model to describe the transmission dynamics of the DENV in human and mosquito populations, they incorporated aquatic stages (i.e., eggs, larvae, pupae) in the mosquito population, and calculated the basic reproduction number to be less than one (which indicate that the disease can dies out in time). Yang et al. ^[24] designed a simple model in order to assess the effect of temperature on the population of the mosquito by using the temperature controlled experiments. Yang et al. ^[25] extended Yang et al. ^[24] by including human compartments and updated most of the parameters to assess the effect of temperature on the risk of the DENV outbreaks. Yang ^[5] developed and analyzed a model encompassing four quiescence stages of the mosquito, and used the model to assess the influence of the quiescent eggs on the transmission dynamics of the mosquito.

In the current study, we focus on the design and analysis of a model, which extends some of the aforementioned studies, for assessing the transmission dynamics of the DENV in human and mosquito populations. Our model includes both the human and mosquito populations, and incorporates asymptomatic ($A_h$) and severe cases (those with DHF/DSS) of the DENV infection ($I_{h_2}$) in human population. We performed model fitting by using the plug-and-play inference framework for the human population, and our results contribute to better understanding of the transmission dynamics of the DENV epidemics in Taiwan in 2014-15 and provide useful guidelines for the design of control strategies in the future.

The paper is organized as follows. The model is formulated in section 2 and analyzed in section 3. Sensitivity analysis results of the full model is also given in section 4. Finally numerical analysis of the human only sub-model is performed with a plug-and-play inference framework in section 5.

2.   Model formulation

2.1.   Dengue epidemic model

The total human population at time $t$, denoted by $N_{h}(t)$, is divided into six compartments of susceptible, $S_h(t)$, exposed (infected but not infectious individuals that can progress to $I_{h1}$ class), $E_h(t)$, asymptomatically infected, $A_h(t)$ (i.e., expose individuals that can progress and becomes infectious without showing any clinical symptoms, which are lumped together for computational conveniences), symptomatically infected, $I_{h_1}(t)$, symptomatically infected with severe DENV (DHF/DSS), $I_{h_2}(t)$, and recovered, $R_h(t)$, humans. Hence,

The total mosquito population at time $t$, denoted by $N_m(t)$, is split into susceptible $S_m(t)$, exposed (infected but not infectious mosquitoes) $E_m(t)$ and symptomatic (infectious mosquitoes) $I_m(t)$ classes. Hence,

The model for the DENV transmission in the human and mosquito populations is given by the following deterministic ordinary differential equations (ODE) systems (2.1). The flow diagram of the model (2.1) is depicted in Figure 1. The state variables and parameters of the model are described in Tables 1 and 2, respectively.

Here, the infection rate for humans ($\lambda_{h}$) and vectors ($\lambda_{m}$) are given by

In model (2.1), $\Pi_h$ is the recruitment rate for human by birth, $\lambda_h$ is the infection rate of susceptible human from infectious mosquitoes, where $a$ is the mosquito biting rate, $b$ is the maximum number of bites a human can receive per unit time, $\beta_h$ is the transmission probability from infected mosquitoes to the susceptible humans while $\beta_m$ is transmission probability from an infectious humans to the susceptible mosquitoes per bite, and $\mu_h$ is the natural death rate of human. The term $\theta$ is the symptomatic ratio, and $(1-\theta)$ remains the fraction of asymptomatic ratio. The parameter $v_h$ measures the rate at which an individuals in the $E_h$ class develop a primary clinical symptoms of the virus, while $\eta$ is a fraction of humans who will not get a severe DENV, $\delta_h$ accounts for the disease-induced death rate of humans from the $I_{h_2}$ class, and $\gamma_i(i = 0, 1, 2)$ measures the rate of recovery of individual from $A_h$, $I_{h_1}$, and $I_{h_2}$, respectively.

The susceptible adult mosquitoes acquire the DENV at a rate $\lambda_m$, this population is decreased by natural death at a rate $\mu_m$. Finally, $\sigma_m$ is a progression rate of the DENV from the $E_m(t)$ class to the $I_m(t)$. It is also worth stating that the model (2.1) accounts for the conservation law of mosquito bites. Thus, the infection rates, $\lambda_h$ and $\lambda_m$, are normalized by the total host population, that is $N_h(t)$ ^[27].

2.2.   Basic properties

The basic properties of the model (2.1) will now be explored. Consider the following equations for the rate of change of the total human $N_h^{\prime}(t)$ and mosquito $N_m^{\prime}(t)$ populations

and

here, the prime, $^\prime$, represent a differentiation with respect to time, $t$.

Furthermore, consider the region,

It can be shown by solving for $N_h$ and $N_m$ in Eqs (2.3)-(2.4) so that all solutions of the system starting in the region $\Omega$ will remain in $\Omega$ for all time $t$ with $t \geqslant 0$. Thus, the region $\Omega$ is positively-invariant, and it is sufficient to consider solutions restricted in $\Omega$. In this region, the usual existence, uniqueness and continuation results hold for the system (2.1) ^[28,29].

3.   Mathematical analysis

3.1.   Disease-free equilibrium and its stability

The disease-free equilibrium (DFE) of the model (2.1) obtained at steady state is given by

Using the next generation matrix method ^[30], we obtained the associated reproduction number of the model (2.1), denoted by ${\mathcal{R}}_{0} = \rho(FV^{-1})$. The term $\rho$ represents the function to find the spectral radius of the next generation matrix, ${\bf{G}} = FV^{-1}$. The matrices $F$ (for the new infection terms) and $V$ (for the remaining transition terms), associated with the model (2.1), are given by

where $n_1 = v_h+\mu_h$, $n_2 = \gamma_h+\mu_h$, $n_3 = \gamma_1+\mu_h$, $n_4 = \delta_h+\gamma_2+\mu_h$, and $n_5 = \sigma_m+\mu_m$. Therefore, the basic reproduction number, ${\mathcal{R}}_{0}$, is given by

where

and

${\mathcal{R}}_{0}$ interpretation: The threshold quantity ${\mathcal{R}}_{0}$ is ecologically and epidemiologically interpreted as follows.

ⅰ. $\frac{ab\beta_h}{N_h}$ accounts for the number of new infected human hosts caused by an infected mosquito over its expected infectious period,

ⅱ. $\frac{v_h}{n_1}$ is the probability that an exposed human survives the exposed stage and move to the infectious stage, $I_{h_1}$,

ⅲ. $\frac{1}{n_3}$ is the average duration in the infectious stage,

ⅳ. $\frac{(1-\eta)\gamma_1}{n_4}$ is the probability that an infectious human survives the infectious stage (including severe stage of the DENV) and move to the recovered class,

ⅴ. $\frac{(1-\theta)}{n_2}$ is the probability that an infected human move to the asymptomatic class,

ⅵ. $ab\beta_m\frac{\sigma_m}{n_5}\frac{S_m^{0}}{\mu_m}$ accounts for the possible number of new infected mosquito caused by an infected human hosts over its exposed infectious period.

The theorem 3.1 below follows the Theorem 2 of Ref ^[30].

Theorem 3.1. The DFE, $E^0$, of the model (2.1), is locally-asymptotically stable (LAS) in $\Omega$ if ${\mathcal{R}}_{0} < 1$, but unstable if ${\mathcal{R}}_{0} > 1$.

3.2.   Existence of endemic equilibrium and backward bifurcation

3.2.1.   Endemic equilibrium and backward bifurcation

The endemic equilibrium (EE) for the system (2.1), which is

in terms of the forces of infection, $\lambda_h^{*}$ and $\lambda_m^{*}$, is given by

Substituting the above equilibrium points (equation (3.6)) into equation (2.2), we have

where

Therefore, one can obtain the positive EE of the model (2.1) by simplifying the quadratic equation (3.7), and substituting only the positive values of the $\lambda_h$ into the EE points. The occurrence of the quadratic equation (3.7) shows that the backward bifurcation (BB) for the model (2.1) exists.

Hence, the following Theorem 3.2 is established.

Theorem 3.2. The model (2.1) has

$i)$ a EE if $C < 0\Leftrightarrow {\mathcal{R}}_{0} > 1$;

$ii)$ a unique EE if $B < 0$ and $C = 0$ or $B^2-4AB = 0$;

$iii)$ two EEs if $B < 0$, $C > 0$ and $B^2-4AC > 0$; or

$iv)$ no EE otherwise.

The existence of the phenomenon of the BB in the model (2.1) highlights a co-existence between a stable disease free-equilibrium (DFE) and a stable endemic equilibrium (EE) even if the basic reproduction number (${\mathcal{R}}_{0}$) is less than unity, this makes the disease control more difficult. Furthermore, the disease eradication would no longer depends on the basic concept of ${\mathcal{R}}_{0}$, i.e., when ${\mathcal{R}}_{0} < 1$, the disease dies out in time; and while ${\mathcal{R}}_{0} > 1$, the disease persists in the community. Although the BB has been firstly shown to exist in the DENV from the work of Garba et al. ^[6] and its courses has been explained extensively by Gumel ^[26]. In this study, our analysis suggests that the BB can be removed if the ${\mathcal{R}}_{0}$ can be equal to one (i.e., ${\mathcal{R}}_{0} = 1$), so that equation (3.7) will be $A\lambda_h+B = 0$ which is linear equation.

3.2.2.   Global stability analysis of the endemic equilibrium

The following Theorem 3.3 is claimed.

Theorem 3.3. The endemic equilibrium (EE), $E^{*}$, is globally-asymptotically stable (GAS) in $\Omega$ when ${\mathcal{R}}_{0} > 1$ provided that

The proof of the above Theorem 3.3 is given in Appendix 7, and the Figure 2 is used to demonstrate the analytical result obtained.

4.   Sensitivity analysis of the model (2.1)

Following the previous studies ^{[13,14,34,35]}, we adopted the partial ranked correlation coefficient (PRCC) for sensitivity analysis. The PRCC of the ${\mathcal{R}}_{0}$ and IAR of the model (2.1) is estimated. The sensitivity analysis results in Figure 3 suggest the top three ranked parameters (i.e., $a$, $b$ and $\Pi_m$) are prioritized in controlling the DENV epidemics in Taiwan.

5.   Numerical simulations: a special case

5.1.   Plug-and-play inference framework

The analyses in this section will be carried out for the special scenario of the model (2.1) in the absence of mosquito compartments (i.e., a human only sub-model). We assumed that there is no birth and death rate ($\Pi_h = \mu_h = 0$), and no disease-induced mortality rate ($\delta_h = 0$) in a short period. We also assumed that the effective contact rate of the simplified version of the model to be $\beta = \beta(t) = ab \cdot m(t)$, where $m(t)$ is the time-dependent mosquito to human ratio. This setting allows us to fit the "SIR" based model to a vector-borne disease under the "absence of vector" model structure. The term of the product of $(a b m)$ is equivalent to the standard Ross-Macdonald malaria model ^[36].

The equivalent (simplified) model without mosquito compartment together with the above assumptions is given by

where the basic reproduction ratio of the simplified model (5.1) is given by ${\mathcal{R}}_{0}(t) = \frac{\beta\theta}{\gamma_1}$.

The weekly number of cases of the $i$-th week is

where the term $\rho$ denotes a constant reporting rate of the DENV cases. The $Z_i$ denotes the theoretical weekly DENV cases yield from the model (5.1).

We model the observed (reported) DENV cases, $C_i$ for the $i$-th week during the study period, as a partially observed Markov process (POMP, also know as hidden Markov model, HMM). Instead of implementing Poisson-distributed priors ^[34], all $C_i$s are assumed to follow over-dispersed Poisson distributions according to the theoretical modelling (Eq (5.2)) outputs, $Z_i$s. Since the rate of Poisson distribution is a Gamma random variable, the observed weekly number of people who are confirmed DENV infections ($C_i$ of the $i$-th week) is a random sample from a Negative-binomial (NB) distribution. Therefore,

where $\psi$ is the over-dispersion parameter for NB distribution under estimation. In this model, we set $L_i (\cdot)$ be the likelihood function for the week, which is the "probability" of $C_i$, given the real cases from simulations $Z_i$ under the NB distribution ^[37,38]. In this section, the model simulations are conducted by using the R (version 3.4.1) package "POMP" ^[39].

The overall log-likelihood, $l$, for the whole-time series is

where $\Theta$ is the parameter vector under estimation. The term $T$ denotes the total number of weeks during the study period. We apply the iterated filtering algorithm with the plug-and-play likelihood-based inference framework to estimate the maximum likelihood estimates (MLE) of $\Theta$, (see ^{[37,40,41,42,43,44,45,46]}. We use the fixed-time-step Euler-multinomial algorithm to simulate the ODE system (5.1) ^{[14,34,37,38,40]}, We compared different models using the small-sample-size corrected Akaike's Information Criterion (AICc) ^[44] as a measurement of the trade-off between model complexity and the goodness-of-fit. The AICc is given by

where $N$ is the number of data points and $k$ is the number of free parameters.

5.2.   Fitting results

The (simplified) epidemic model (5.1) is fitted to the dengue cases time series in Kaohsiung and Tainan cities in 2014–2016. We proposed the time-dependent basic reproduction number (i.e., ${\mathcal{R}}_{0}(t)$) driven by variable mosquito abundance and dengue transmission rate (e.g., mosquito biting rate, etc). Here, we compared two numerical approaches to fit the epidemic data and reconstruct the ${\mathcal{R}}_{0}(t)$ series, namely,

● the cubic spline reconstruction approach, and

● the square wave reconstruction approach.

The cubic spline approach is defined from the previous studies ^[14,47]. We used $p_q$ to denote the number of nodes in the cubic spline of ${\mathcal{R}}_{0}(t)$ to be reconstructed. We tested different possible $p_q$'s with the aim to find the value of $p_q$ that leads to the lowest AICc. We further restrict that Kaohsiung and Tainan shared exactly the same cubic spline function of ${\mathcal{R}}_{0}(t)$. Hence, the two cities were fitted simultaneously with the same set of parameter estimates. For the cubic spline approach, the fitting is conducted with completely the same time series of ${\mathcal{R}}_{0}(t)$ (i.e., same timing, scale, and all settings are the same) for both cities. And we find the goodness-of-fit (in term of the MLL) can still be improved for both cities. Even though the spline can yield a relatively flexible (in term of the degree of freedom of the parameters to be estimated) ${\mathcal{R}}_{0}(t)$ reconstruction outcome, the fittings are still not satisfied. Thus, we speculate a non-synchronized ${\mathcal{R}}_{0}(t)$ would be more reasonable. The fitting results of the ${\mathcal{R}}_{0}(t)$ reconstructed by the cubic spline function are shown in Figure 4, and the parameter estimates are summarized in Table 5. The square wave is defined as a baseline transmission rate plus a constant additional value for each year. The additional values are different for both Kaohsiung and Tainan in the same year, and also vary in different years. To explore the possible non-synchronization of ${\mathcal{R}}_{0}(t)$'s in Kaohsiung and Tainan cities, we further allow the changing time of the square wave to have different estimates in the two cities. Thus, the reconstructed ${\mathcal{R}}_{0}(t)$s in Kaohsiung and Tainan are expected to have different shape which indicates that the ${\mathcal{R}}_{0}(t)$'s of the two cities are not perfectly synchronized. In Tainan, we fix the time of the occurrence of the ${\mathcal{R}}_{0}(t)$ high values (i.e., with the additional value) to be May 1 and October 31 in both 2014 and 2015. In Kaohsiung, our estimates are that the high values of ${\mathcal{R}}_{0}(t)$ occurred between March 16 and October 31 in 2014, and between May 1 and November 25 in 2015. These time shifts are crucial to achieve good fitting outcomes in Kaohsiung city. Thus, this finding indicates that the DENV outbreaks in Kaohsiung turned on early in 2014, and turned off late in 2015. For the square-wave approach, the fitting is conducted with similar series of ${\mathcal{R}}_{0}(t)$'s for the two cities. The timing (of change) and the scale of ${\mathcal{R}}_{0}(t)$'s were slightly different. Although this approach is less flexible as the spline approach (i.e., with 4 degrees of freedom decreased), it leads to better fitting results with MLL improved by 10 units. And together, the AIC improves (i.e., decreases) by 24 units. This indicates that the non-synchronized square-wave approach is a strong and significant improvement from the synchronized spline approach. The fitting results of ${\mathcal{R}}_{0}(t)$ reconstructed by the square wave function are shown in Figure 5, and the parameter estimates are summarized in Table 6. Therefore, for the spline versus square-wave approach, the non-synchronized square-wave approach has better goodness-of-fit (i.e., higher MLL) with less flexibility (i.e., less degree of freedom). In other words, the non-synchronized square-wave approach has better fitting performance (in term of the AIC). Hence, we conclude that non-synchronized ${\mathcal{R}}_{0}(t)$'s were likely to occur in the two cities.

Our (simplified) DENV epidemic model well-explained the temporal patterns of the DENV outbreaks in the two cities in Taiwan under biologically reasonable conditions (i.e., set of model parameters). We discovered that the epidemics in the two cites can be reconstructed by similar transmission series (i.e., ${\mathcal{R}}_{0}(t)$), and also estimated with the same model structure. The square wave approach not only reduced the number of free parameters by 4, but also increased the log likelihood by 10 units. Hence, its AICc decreases roughly by 24 units from the cubic spline reconstruction, which is a significant improvement.

Although we find that the estimated transmission rates exhibited similar seasonal patterns in Kaohsiung and Tainan, they are neither completely synchronized between the two cities despite their proximity, nor periodically in-phase perfectly. If they were completely synchronized or perfectly periodic, the cubic spline approach would identify such a pattern reflected by an equivalent level of goodness-of-fit as in the square wave approach. We conclude that the time lags between the seasonal waves of the transmission rates (${\mathcal{R}}_{0}(t)$) in the two cities are likely to occurred. Other than the pathogenic factors that can causes the phase difference in the transmissibility ^[48], these time lags could be due to the heterogeneity in the local conditions in Kaohsiung and Tainan cities.

6.   Conclusions

We formulated and analyzed a mathematical model to study the transmission dynamics of the dengue outbreaks in Kaohsiung and Tainan cities in Taiwan in 2014–2015. We showed that the model (2.1) has two equilibria, of which the disease-free equilibrium (DFE) is locally asymptotically stable (LAS) whenever the basic reproduction number (${\mathcal{R}}_{0} < 1$), and unstable if the ${\mathcal{R}}_{0} > 1$. And the endemic equilibrium (EE) is globally asymptotically stable (GAS) in the region of attraction, $\Omega$, whenever the ${\mathcal{R}}_{0} > 1$ (numerical examples, see Figure 2, is used to demonstrate the analytical results obtained). Our analysis shows the existence of backward bifurcation (BB) phenomenon of the model (2.1), a situation where a stable disease free-equilibrium co-exists with a stable endemic equilibrium even when the basic reproduction number (${\mathcal{R}}_{0}$) is less than one, which makes the disease control even more difficult and is no longer dependent on the basic concept of the ${\mathcal{R}}_{0}$ (i.e., ${\mathcal{R}}_{0} < 1$ the disease dies out in time, while ${\mathcal{R}}_{0} > 1$ the disease persist to the community). In this paper, the analysis suggest that the backward bifurcation can be removed when the basic reproduction number is equal to one (i.e., ${\mathcal{R}}_{0} = 1$).

We employed the plug-and-play statistical inference framework on the special case of the model (2.1) in the absence of mosquito. We fitted the simplified model (5.1) to the weekly number of cases in the Kaohsiung and Tainan Cities in Taiwan (obtained from the National Infectious Disease Statistics System ^[49], see Figures. 4 and 5). Our (simplified) dengue epidemic model is able to explained the temporal patterns of the dengue outbreaks in the two cities in Taiwan under biologically reasonable conditions (i.e., set of model parameters). We obtained that the epidemics in the two cites can be reconstructed by the same transmission series (i.e., ${\mathcal{R}}_{0}(t)$) and also estimated with the same model structure. This findings suggests the similarity of the dengue outbreak in a two different places during the same period of time, which could be due to the similar demographic and/or meteorological conditions. The proposed two reconstruction approaches presents (roughly) equivalent goodness-of-fit in terms of the values of the log-likelihood. However, the AICc improves (decreases) roughly by 24 units from the cubic spline to the square wave reconstruction. We also found that despite the proximity in Kaohsiung and Tainan, the estimated transmission rates were neither completely synchronized, nor periodically in-phase perfectly in the two cities. The time lags between the seasonal waves of the transmission rates in the two cities were likely to occur. These time lags could be due to the differences in their local conditions. Furthermore, we show that the non-synchronized square-wave approach has better goodness-of-fit (i.e., higher MLL) with less flexibility (i.e., less degree of freedom). In other words, the non-synchronized square-wave approach has a better fitting performance (in term of the AIC) than the cubic spline approach. Hence, we concluded that non-synchronized ${\mathcal{R}}_{0}(t)$s were likely to occur in the two cities.

The sensitivity analysis results shows the top three ranked parameters (i.e., $a$, $b$ and $\Pi_m$), that is, the mosquito biting rate, the maximum number of bites a human can receive per unit time, and the mosquito recruitment rate are the key parameters to be prioritized for controlling the Taiwan dengue outbreaks. This further suggests that proper sanitation of mosquito breeding sites and avoiding the mosquito bites are the key control measures to future dengue outbreaks in Taiwan.

List of Abbreviations

We listed the abbreviations used in Table 4.

Acknowledgments

Zhen Jin is partly supported by Shanxi Scientific and Technology Innovation Team (201705D15111172 & 201805D131012-1). We thank reviewers for helpful comments.

Conflict of interest

The authors declare that they have no conflict of interest in this paper.

Appendices

A1.   The proof of theorem 3.3

Proof. Define a Lyapunov function as follows

where $J_1 = 2$, $J_2 = \frac{1}{1-\theta}$, $J_3 = \frac{1}{(\theta)}$, $J_4 = \frac{\lambda_{h}^* S_h^*}{v_h E_{h}^{*}}$, $J_5 = \frac{\lambda_{h}^* S_h^*}{(1-\eta)\gamma_1 I_{h1}^*}$, $K_1 = K_2 = 1$, and $K_3 = \frac{\lambda_{m}^*S_m^*}{\sigma_m E_m^*}$.

Thus, the Lyapunov derivative computed along solutions of the system (2.1) is given by

Hence,

Following ^[50,51], the function $v(x) = 1-x+\ln x$, if $x > 0$, it leads to $v(x) \leqslant 0$. And if $x = 1$, we have $v(x) = 0$. Thus, $x-1 \geqslant \ln(x)$ for any $x > 0$. Then, we also have that

Similarly, we have

and

Furthermore, we have that

and also

as well as

Hence, the equations (A1-1)-(A1-9) together with condition (3.8) ensure that $\frac{dV}{dt}\leqslant 0$. Furthermore, the strict inequality $\frac{dV}{dt} = 0$ holds only for $S_{h} = S_{h}^{*}$, $E_h = E_{h}^{*}$, $A_h = A_{h}^{*}$, $I_{h_1} = I_{h_1}^{*}$, $I_{h_2} = I_{h_2}^{*}$, $R_h = R_{h}^{*}$, $S_m = S_{m}^{*}$ $E_m = E_{m}^{*}$ and $I_m = I_{m}^{*}$. Thus, the endemic equilibrium state $E^{*}$ is the only positive invariant set to the system (2.1) contained entirely in $\biggl\{(S_h, E_h, A_h, I_{h_1}, I_{h_2}, R_h, S_m, E_m, I_m) \in \Omega : S_h = S_{h}^{*}, E_h = E_{h}^{*}, A_h = A_{h}^{*}, I_{h_1} = I_{h_1}^{*}, I_{h_2} = I_{h_2}^{*}, R_h = R_{h}^{*}, S_m = S_{m}^{*}, E_m = E_{m}^{*}, I_m = I_{m}^{*} \bigamma\}$. Therefore, it follows from the LaSalle's invariance principle ^[52] that every solutions to the equations in (A1-2) with initial conditions in $\Omega$ converge to stable endemic equilibrium point, $E^{*}$, as $t\rightarrow\infty$. Hence, the positive endemic equilibrium is globally asymptotically stable.

A2.   Summary table of ition / Descrnitial state values and parameters estimation results of the model (5.1)