Mathematical evaluation of the role of cross immunity and nonlinear incidence rate on the transmission dynamics of two dengue serotypes

Abstract

Dengue fever is a common disease which can cause shock, internal bleeding, and death in patients if a second infection is involved. In this paper, a multi-serotype dengue model with nonlinear incidence rate is formulated to study the transmission of two dengue serotypes. The dynamical behaviors of the proposed model depend on the threshold value $R_{{0}}^{{n}}$ known as the reproductive number which depends on the associated reproductive numbers with serotype-1 and serotype-2. The value of $R_{{0}}^{{n}}$ is used to reflect whether the disease dies out or becomes endemic. It is found that the proposed model has a globally stable disease-free equilibrium if $R_{{0}}^{{n}}\leq 1$, which indicates that if public health measures that make (and keep) the threshold to a value less than unity are carried out, the strategy in disease control is effective in the sense that the number of infected human and mosquito populations in the community will be brought to zero irrespective of the initial sizes of sub-populations. When $R_{{0}}^{{n}}>1$, the endemic equilibria called the co-existence primary and secondary infection equilibria are locally asymptotically stable. The effects of cross immunity and nonlinear incidence rate are explored using data from Thailand to determine the effective strategy in controlling and preventing dengue transmission and reinfection.

1 Introduction

Dengue is a viral mosquito-borne infection which in recent years has become a major international public health concern, a leading cause of illness and death in the tropics and subtropics with more than 50 million dengue fever cases per year [1, 2]. There are four distinct serotypes of the dengue virus (DEN 1, DEN 2, DEN 3, and DEN 4) that coexist in many endemic areas [3, 4], although recently a potential fifth strain has been discovered [5]. Infection with one serotype affords life-long immunity to that serotype [6], but only cross immunity or none at all to the other three serotypes [6, 7]. Patients reinfected with a different serotype, called secondary infection, face an increased risk of developing DHF and DSS [8]. Recently, cross protection, or cross immunity between serotypes, has been conjectured to play a role in the dynamics of dengue in the sense that while a primary dengue infection with a particular serotype may confer life-long immunity to that serotype, it may also confer temporary cross immunity to the other serotypes. It means that the cross immunity may be complete (individuals cannot contract a secondary infection during the cross immune period) or partial (cross immune individuals have a reduced but not zero probability to contract the disease) [9, 10]. Cross immunity may be the result of an immunological response to the disease. It acts to reduce the susceptibility to a secondary infection lowering the effective probability for reinfection to happen [11]. In general, the length of the cross immunity period may vary depending on the disease. In case of dengue fever, the cross immunity may last from two to nine months [12] after the antibodies have dropped to sufficiently low levels that allow infection with other serotypes and subsequent ADE. The study in Bianco et al. [13] showed that the main risk factor associated with DHF/DSS in secondary infection is the presence of pre-existing dengue antibodies at sub-neutralizing level after the cross immunity period. This increasing risk has been attributed to the effect of cross immunity [14]. These indicate that the cross immunity, thus, plays a crucial role in the co-circulation of serotypes and the pathogen diversity [15, 16].

Today, dengue is regarded as the most prevalent and rapidly spreading mosquito-borne viral disease of human beings. This leads to an upsurge in research on dengue virology, pathogenesis, and immunology and in development of antivirals and vaccines. Especially, mathematical models are one useful tool to investigate the cause of epidemic and to suggest the best way to control and prevent dengue [9, 16–27]. Kooi et al. [22] studied an asymmetric two-strain dengue model for predicting characteristic dynamic behavior and chaos occurring for smaller parameter regimes. Zheng et al. [24] formulated a two-strain dengue model and suggested that the measures of enhancing awareness of the infected and susceptible human self-protection should be taken, and the mosquito control measure is necessary in order to prevent the transmission of dengue virus from mosquitoes to humans. Souto-Maior [25] studied the transmission of dengue virus using multiple-serotype models and showed that the proposed model can sustain oscillations in the absence of any of the latter effects or stochasticity. This model was applied to simulate both epidemiological and viral evolution for epidemiological surveillance. These research works constructed the epidemiological model and employed bilinear and standard incidence rates (given by βSI and $\frac{\beta \mathrm{SI}}{N}$, respectively), which may be a good approximation if the number of available partners is large enough and everybody could not make more contacts than is practically feasible [28, 29]. Woodall et al. [16] presented a new modeling framework based on SIR model with enhancement to study cross enhancement between dengue serotypes which may be influencing the epidemic oscillations. González Morales et al. [23] studied the asymptotic and dynamical behavior of a two-strain dengue model under the application of a vaccine and indicated that vaccination and cross immunity period are seen to decrease the frequency and magnitude of outbreaks but in a differentiated manner with specific effects depending upon the interaction vaccine and the strain type. Meanwhile, Anggriani et al.[26] studied the effect of reinfection with the same serotype on dengue transmission dynamics by developing a multi-strain dengue mathematical model, which suggested that reinfection with the same serotype may be one of the underlying factors causing an increase in the number of secondary infections. Further, the dengue model with standard incidence rate proposed by Janreung et al. [27] was applied in order to investigate only the effect of vaccine parameters in controlling reinfection dengue by vaccination. In addition, empirical data in Thailand 2018 [30] showed that secondary infection leads to severe disease which is a major cause of death from dengue. Therefore, the reinfection between two serotypes is a major cause of yearly outbreaks of dengue [8, 30]. However, a number of authors have pointed out that a nonlinear incidence rate that includes the behavioral change and crowding effect of the infected individuals and also prevents unboundedness of the contact rate may be more realistic than others during the disease transmission process [9, 17, 31–41]. In addition, studies have found that epidemic models with nonlinear incidences have more complicated dynamics than those with bilinear or standard incidence. For example, Capasso and Serio [31] suggested a saturated nonlinear incidence $\frac{\beta I S}{1+\alpha I}$ to better model the cholera epidemic spread in Bari in 1973. This is important because the number of effective contacts between infective and susceptible individuals may saturate at high infective level due to overcrowding of infective individuals or due to protective measures enclosed by susceptible. Since then Cunningham [32] introduced an incidence rate $k(\mathrm{SI})^{p}\ (p > 1)$ to point out that there may exist periodic solution in a deterministic model for measles. Li et al. [33] studied the influence of nonlinear incidence rates $\frac{\beta I^{l} S}{1+\alpha I^{h}}$ upon the behavior of SIRS epidemiological models. Roop-O et al. [38] studied the influence of nonlinear incidence rates $\frac{\beta I S}{1+\alpha I}$ upon the backward bifurcation for a malaria model with temporary immunity. The study results suggested the different control strategies of infectious diseases. For dengue fever, infection by one serotype confers life-long immunity to that serotype and a period of temporary cross immunity to other serotypes. Infection with other kinds of dengue viruses in a person who already had a dengue fever in the past is called secondary dengue and has been found to be more severe disease than the primary infection due to cross immunity between serotypes [7, 9, 17]. In this paper, therefore, a mathematical model of dengue fever epidemiology with two serotypes is constructed by incorporating cross immunity and nonlinear incidence rate to explore the risk of secondary infection.

The paper is organized as follows: In Sect. 2, a mathematical model of two serotype dengue fever is constructed based on the biological aspects of dengue fever epidemiology. In Sect. 3, the formulated model is analyzed for the existence and stability of its equilibria. The reproductive number is derived using the next generation. By constructing a Lyapunov function, the disease-free equilibrium of the formulated model is globally asymptotically stable if the reproductive number is less than or equal to unity. The existence of endemic equilibrium is determined and its local stability is proved by applying the center manifold theory [42]. In Sect. 4, numerical simulations are illustrated to determine the appropriate parameter values used for the formulated model and to determine the parameter which is sensitive to disease prevalence. The effect of temporary cross immunity and nonlinear incidence rate for disease prevention and secondary infection control is investigated. Finally, discussion and conclusion are presented in Sect. 5.

2 Multi-serotype dengue model

In this section, the dengue model is formulated based on the situation of dengue transmission as follows. Dengue fever is caused by any one of four types of dengue viruses spread by mosquitoes. When a mosquito bites a person infected with a dengue virus, the virus enters the mosquito. When the infected mosquito then bites another person, the virus enters that person’s bloodstream. Thus, the study populations are human population and mosquito population. Although, there are four types of dengue virus, infection with any one of the viruses provides lifetime immunity to future infections from the same virus but it will have cross protection (cross immunity) against other types for a short period which leads to the possibility of secondary infection. Moreover, secondary infection will increase the risk of developing dengue hemorrhagic fever. Thus, temporary cross immunity can lead to subsequent infection with a different serotype. To study the influential factors to secondary dengue infection, the infection is divided into two parts: primary dengue infection and secondary dengue infection, respectively. The total human population is divided into two groups: human populations who are infected with dengue with serotype-i for the first time (the primary infection) and populations who are infected with dengue with serotype-j for the second time. In the primary infection, human populations consist of humans susceptible to both serotypes ($S_{{_{H}}}$), i.e., susceptible to dengue serotype-1 and dengue serotype-2; exposed humans with serotype-i ($E_{{H{{i}}}}$), infected humans with serotype-i ($I_{{H{{i}}}}$), and recovered humans with serotype-i ($R_{{H{{i}}}}$) who recover from the first infection with serotype-i, that is, if individuals are infected with the same dengue virus serotype they become immune to future infections. However, if individuals are infected subsequently with a different serotype, immunity wanes over time, which increases the risk of developing dengue hemorrhagic fever [43]. Thus, these sub-populations become susceptible to infection with a different serotype. For this reasons, the sub-populations in the second infection consist of recovered humans with serotype-i ($R_{{H{{i}}}}$) and then becoming susceptible to other serotypes (serotype-j), exposed with serotype-j when the first infection was caused by serotype-i ($E_{{H{{ij}}}}$), infected with serotype-j when the first infection was caused by serotype-i ($I_{{H{{ij}}}}$); and recovered humans from the secondary infection ($R_{{H{{22}}}}$) who then have life-long immunity to both serotypes, respectively. Thus, the total population is given by

$$\begin{aligned} N_{{H}}(t)={}&S_{{H}}(t)+\sum _{i=1}^{2} \bigl(E_{{Hi}}(t)+I_{{Hi}}(t)+R_{{Hi}}(t) \bigr) \\ &{}+\sum_{i=1}^{2}\sum _{\substack{j=1 \\ i\neq j}}^{2} \bigl(E_{{Hij}}(t)+I_{{Hij}}(t) \bigr)+R_{{H22}}(t). \end{aligned}$$

(2.1)

For mosquito population, there is a pool of mosquitoes that are exposed and yet incapable of passing on dengue because virus incubation is 4–10 days and after this period an infected mosquito is capable of transmitting the virus for the rest of its life. The total mosquito population at time t, therefore, denoted by $N_{{V}}$, is split into five sub-populations: susceptible mosquito $(S_{{V}}(t))$, exposed mosquito with serotype-i$(E_{{Vi}})$, and infected mosquito with serotype-i$(I_{{Vi}}(t))$, so that

$$\begin{aligned} N_{{V}}(t)=S_{{V}}(t)+\sum _{i=1}^{2} \bigl(E_{{Vi}}(t)+I_{{Vi}}(t) \bigr). \end{aligned}$$

(2.2)

It is noted that the population of recovered mosquitoes is not consider, because once a mosquito is infected with one serotype it never recovers from the infection and it cannot be reinfected with a different serotype since the life cycle of Aedes aegypti is short (approximated one and a half to three weeks) [19]. The description of variables in (2.1) and (2.2) are summarized in Table 1.

Table 1 Description of variables in model (2.4)–(2.5)

2.1 Nonlinear incidence rate

The incidence rate in an epidemiological model is the rate at which susceptible become infectious, which plays a very important role in analyzing the transmission of diseases. Several researchers thus used the standard incidence for studying dengue transmission, (see for instance [7, 9, 17, 20] and the references therein). However, dengue is a viral infection caused by four types of viruses (DENV-1, DENV-2, DENV-3, DENV-4), infection with one serotype will provide life-long protection to future infections from the same virus but only short-term cross immunity to the other types, leading to the possibility of secondary infections by other serotypes and developing into severe dengue. In addition, the dengue viruses are transmitted through the bite of infected female mosquitoes that feed both indoors and outdoors during the daytime (from dawn to dusk). These mosquitoes thrive in areas with standing water, including puddles, water tanks, containers, and old tires. Lack of reliable sanitation and regular garbage collection are also the causes of contribution to the spread of the mosquitoes. These are due to influence on nonlinearities in dengue transmission. For these reasons, the standard incidence rate should be modified into a nonlinear incidence rate order to be more realistic (see [31–33, 38, 40, 41]). Moreover, Bartley et al. [44] suggested that the strongest influences on the pattern of dengue infection and its seasonal variation were duration of infectiousness of the host, vector mortality, and biting rate. Let $\varphi _{{i}}$ and $\varphi _{{Vi}}$ represent the nonlinear incidence rate with serotype-i from an infected mosquito to a susceptible human and the nonlinear incidence rate with serotype-i from an exposed human and an infected human to a susceptible mosquito, respectively. Motivated by Capasso and Serio [31], the nonlinear incidence rates for humans and mosquitoes are defined by

$$\begin{aligned} \varphi _{{i}}= \frac{b\beta _{{i}}I_{{Vi}}}{1+\alpha _{{V}}I_{{Vi}}} \quad\text{and}\quad \varphi _{{Vi}}= \frac{b\beta _{{Vi}} (\eta _{{Hi}} (E_{{Hi}}+E_{{Hji}} )+I_{{Hi}}+I_{{Hji}} )}{1+\alpha _{{H}}(I_{{Hi}}+I_{{Hji}})} \end{aligned}$$

(2.3)

for $i,j=1,2, i\neq j$. The modification parameter $\eta _{{Hi}}$ accounts for the relative infectiousness of exposed humans with serotype-i in relation to infected humans with serotype-i, and $\alpha _{{V}}$, $\alpha _{{H}}$ denote the parameters measuring the inhibitory effect for human and mosquito, respectively. Notice that when $\alpha _{{V}}=0$ and $\alpha _{{H}}=0$, the nonlinear incidence rate (2.3) becomes the bilinear incidence rate.

2.2 Cross immunity

Cross immunity introduces significant challenges to researchers looking to create an accurate and identifiable epidemiological model of the disease. Temporary cross immunity can give rise to complex temporal dynamics in disease incidence, creating oscillating time series from which it is difficult to elicit information on or draw conclusions about parameters that govern the underlying epidemiological processes. It is known that a person currently infected with dengue or recovered from one serotype may experience a double infection or reinfection by a different serotype, but at a rate reduced by cross immunity. The phylogenetic data further suggest that there is a clade-specific immunological reaction between two serotypes. This cross immunity is encapsulated by the parameter λ, which modifies the probability of secondary infection. Values of λ between 0 and 1 represent cross protection with $\lambda =0$ giving complete protection against secondary infection and $\lambda =1$ no protection [45].

2.3 Multi-serotype dengue model with a nonlinear incidence rate

A flow diagram of the transmission of multi-serotype dengue in the primary and secondary infections is shown in Fig. 1. The model of dengue with two serotypes and nonlinear incidence rates, from Fig. 1, is described by the following system of differential equations: for $i,j = 1,2$ and $i\neq j$.

Figure 1

Schematic diagram of multi-serotype dengue model

Human population:

$$\begin{aligned} \left. \textstyle\begin{array}{l} {\frac{dS_{{H}}}{dt}} = {\varPi _{{H}}-\sum_{i=1}^{2} \varphi _{i} S_{{H}}-\mu _{{H}}S_{{H}},} \\ {\frac{dE_{{Hi}}}{dt}} = {\varphi _{i} S_{{H}}-K_{{i}}E_{{Hi}},} \\ {\frac{dI_{{Hi}}}{dt}} = {\sigma _{{i}}E_{{Hi}}-K_{{i+2}}I_{{Hi}},} \\ {\frac{dR_{{Hi}}}{dt}} = {\gamma _{{i}}I_{{Hi}}-\lambda \varphi _{j} R_{{Hi}} -\mu _{{H}}R_{{Hi}},} \\ {\frac{dE_{{Hij}}}{dt}} = { \lambda \varphi _{j} R_{{Hi}}-K_{{j}}E_{{Hij}},} \\ {\frac{dI_{{Hij}}}{dt}} = {\sigma _{{j}}E_{{Hij}}-K_{{j+2}}I_{{Hij}},} \\ {\frac{dR_{{H22}}}{dt}} = {\sum_{i=1}^{2}\sum_{\substack{j=1 \\ i\neq j}}^{2}\gamma _{{i}}I_{{Hji}}-\mu _{{H}}R_{{H22}},} \end{array}\displaystyle \right\} \end{aligned}$$

(2.4)

Mosquito population:

$$\begin{aligned} \left. \textstyle\begin{array}{l} {\frac{dS_{{V}}}{dt}} = \varPi _{{V}}- {\sum_{i=1}^{2} \varphi _{{Vi}} S_{{V}} -\mu _{{V}}S_{{V}}}, \\ {\frac{dE_{{Vi}}}{dt}} = {\sum_{i=1}^{2} \varphi _{{Vi}} S_{{V}}-K_{{i+4}}E_{{Vi}},} \\ {\frac{dI_{{Vi}}}{dt}} = {\sigma _{{Vi}}E_{{Vi}}-\mu _{{V}}I_{{Vi}}}, \end{array}\displaystyle \right\} \end{aligned}$$

(2.5)

where $K_{{i}} = \sigma _{{i}}+\mu _{{H}}, {}K_{{i+2}} = \gamma _{{i}}+ \delta _{{i}}+\mu _{{H}}, {}K_{{i+4}} = \sigma {_{{Vi}}}+\mu _{{V}}$. System (2.4)–(2.5) is called a multi-serotype dengue model. The dynamics of human populations in (2.4) are described as follows. The susceptible human population is generated via recruitment of humans (by birth or immigration) into the community at a constant rate $\varPi _{{H}}$. This population is decreased due to the primary infection with serotype-i, which can be acquired via effective contact with an infectious mosquito with serotype-i at a rate $\varphi _{{i}}$ and the natural death at a rate $\mu _{{H}}$. It is also assumed that all classes of human population die at the same natural death rate $\mu _{{H}}$. In the primary infection by serotype-i, exposed humans ($E_{{Hi}}$) develop clinical symptoms of dengue and move to the infectious class at a rate $\sigma _{{i}}$. Infected humans ($I_{{Hi}}$) recover and move into the recovered class, which has life-long immunity to serotype-i, but partial cross immunity to serotype-j, at a rate $\gamma _{{i}}$ and suffer disease-induced death by infection with serotype-i at a rate $\delta _{{i}}$. In the secondary infection by serotype-j, recovered humans with serotype-i ($R_{{Hi}}$) in the primary infection are susceptible to dengue with the second serotype-j. Due to the cross immunity, this population can be infected with serotype-j and becomes exposed human population at the reduced rate $\lambda \varphi _{{j}}$, $0<\lambda <1$. Exposed humans ($E_{{Hij}}$) develop clinical symptoms of dengue by the secondary infection with serotype-j and move to the infectious class at a rate $\sigma _{{j}}$. Infected humans ($I_{{Hij}}$) recover and move into the recovered class ($R_{{H22}}$) with life-long immunity to both serotypes at the rate $\gamma _{{j}}$. This population suffers the disease-induced death by infection with serotype-j at a rate $\delta _{{j}}$.

For the mosquito population in (2.5), the susceptible mosquito population is generated by birth at a constant rate $\varPi _{{V}}$. Susceptible mosquitoes are infected by both serotypes by biting the exposed and infected humans with serotype-i at the rate $\varphi _{{Vi}}$. After virus incubation for $4-10$ days, the exposed mosquito with serotype-i develops symptoms of disease and becomes the infected mosquito with serotype-i at the rate $\sigma _{{Vi}}$. It is assumed that all mosquito populations die due to their finite life span (natural death) at the same rate $\mu _{{V}}$. It is noted that the index i refers to the virus serotype (DEN-1 or DEN-2 or DEN-3 or DEN-4), and the new index j refers to different virus serotypes. The state variables and parameters of the multi-serotype dengue model are described in Tables 1 and 2, respectively.

Table 2 Description of parameters and values used for the multi-serotype dengue model (2.4)–(2.5)

Since multi-serotype model (2.4)–(2.5) monitors human and vector populations, all the associated parameters and state variables are nonnegative. The region of biological interest, then, is given in the following theorem.

Theorem 2.1

The closed set

$$\begin{aligned} \varOmega = \varOmega _{{1}}\cup \varOmega _{{2}}\subset \Re ^{12}_{+} \times \Re ^{5}_{+}, \end{aligned}$$

where

$$\begin{aligned} &\varOmega _{{1}} = \biggl\{ (S_{{H}},E_{{H1}},I_{{H1}},R_{{H1}},E_{{H12}},I_{{H12}},E_{{H2}},I_{{H2}},R_{{H2}},E_{{H21}},I_{{H21}}, R_{{H22}} )\in \Re ^{12}_{+}: \\ &\phantom{\varOmega _{{1}} =} N_{{H}}\leq \frac{\varPi _{{H}}}{\mu _{{H}}} \biggr\} , \\ &\varOmega _{{2}}= \biggl\{ (S_{{V}},E_{{V1}},I_{{V1}},E_{{V2}},I_{{V2}} )\in \Re ^{5}_{+}:N_{{V}}=\frac{\varPi _{{V}}}{\mu _{{V}}} \biggr\} , \end{aligned}$$

is positively invariant and attracting.

Proof

Let $\mathcal{P_{{H}}} = (S_{{H}},E_{{H1}},I_{{H1}},R_{{H1}},E_{{H12}},I_{{H12}},E_{{H2}},I_{{H2}},R_{{H2}},E_{{H21}},I_{{H21}}, R_{{H22}} )\in \Re ^{12}_{+}$ be a solution of human system (2.4) with nonnegative initial conditions. Adding all equations in (2.4) yields

$$\begin{aligned} \frac{dN_{{H}}}{dt}=\varPi _{{H}}-\mu _{{H}}N_{{H}}- (\delta _{{1}}I_{{H1}}+ \delta _{{1}}I_{{H21}}+ \delta _{{2}}I_{{H2}}+\delta _{{2}}I_{{H12}} ). \end{aligned}$$

(2.6)

It follows from (2.6) that

$$\begin{aligned} \varPi _{{H}}- (\mu _{{H}}+\delta _{{1}}+\delta _{{2}} )N_{{H}} \leq \frac{dN_{{H}}}{dt} \leq \varPi _{{H}}-\mu _{{H}}N_{{H}}. \end{aligned}$$

(2.7)

Using Gronwall’s inequality and solutions of two linear first order equations in (2.7), we get

$$\begin{aligned} \frac{\varPi _{{H}}}{\mu _{{H}}+\delta _{{1}}+\delta _{{2}}}+N_{{H}}(0)e^{- (\mu _{{H}}+\delta _{{1}}+\delta _{{2}} )t}\leq N_{{H}}(t) \leq N_{{H}}(0)e^{-\mu _{{H}}t}+\frac{\varPi _{{H}}}{\mu _{{H}}}, \end{aligned}$$

(2.8)

where $N_{{H}}(0)$ represents the initial value of a total human population (2.1) at time $t=0$. Thus, as $t\rightarrow \infty$,

$$\begin{aligned} \frac{\varPi _{{H}}}{\mu _{{H}}+\delta _{{1}}+\delta _{{2}}} \leq N_{{H}}(t) \leq \frac{\varPi _{{H}}}{\mu _{{H}}}. \end{aligned}$$

Therefore, all feasible solutions of the human population (2.4) enter the region $\varOmega _{{1}}$:

$$\begin{aligned} \varOmega _{{1}} ={}& \biggl\{ (S_{{H}},E_{{H1}},I_{{H1}},R_{{H1}},E_{{H12}},I_{{H12}},E_{{H2}},I_{{H2}},R_{{H2}},E_{{H21}},I_{{H21}}, R_{{H22}} )\in \Re ^{12}_{+}: \\ &{}\frac{\varPi _{{H}}}{\mu _{{H}}+\delta _{{1}}+\delta _{{2}}}+ \varepsilon \leq N_{{H}}\leq \frac{\varPi _{{H}}}{\mu _{{H}}}+ \varepsilon , \forall \varepsilon >0 \biggr\} . \end{aligned}$$

Further, let $\mathcal{P_{{V}}} = (S_{{V}},E_{{V1}},I_{{V1}},E_{{V2}},I_{{V2}} )\in \Re ^{5}_{+}$ be a solution of the mosquito system (2.5) with nonnegative initial conditions. Adding all equations in (2.5) yields

$$\begin{aligned} \frac{dN_{{V}}}{dt} = \varPi _{{V}}-\mu _{{V}}N_{{V}}. \end{aligned}$$

(2.9)

Obviously, the total mosquito population $N_{{V}}(t)$ approaches $\varPi _{{V}}/\mu _{{V}}$ as $t\rightarrow \infty$. Then, all feasible solutions of the mosquito population (2.2) enter the region $\varOmega _{{2}}$:

$$\begin{aligned} \varOmega _{{2}}= \biggl\{ (S_{{V}},E_{{V1}},I_{{V1}},E_{{V2}},I_{{V2}} )\in \Re ^{5}_{+}:N_{{V}}=\frac{\varPi _{{V}}}{\mu _{{V}}} \biggr\} . \end{aligned}$$

Taking $N=\operatorname{min} \{N_{{H}},N_{{V}} \}$, it follows that all possible solutions of system (2.4)–(2.5) will enter the region $\varOmega =\varOmega _{{1}}\cup \varOmega _{{2}}$. □

Hence, the region of biological interest Ω is positively invariant under the flow induced by system (2.4)–(2.5). Furthermore, the existence and continuation results of system (2.4)–(2.5) hold in Ω. Thus, model (2.4)–(2.5) is mathematically and epidemiology well posed. It is, therefore, sufficient to consider the dynamics of the flow generated by system (2.4)–(2.5) in Ω.

3 Analysis of multi-serotype dengue model

3.1 Stability of disease-free equilibrium

Multi-serotype dengue model (2.4)–(2.5) has a disease-free equilibrium (DFE) denoted by $P_{{n}}^{{o}}$ which is given by

$$\begin{aligned} P_{{n}}^{{o}} & = \bigl(S_{{H}}^{{o}},E_{{H1}}^{{o}},I_{{H1}}^{{o}},R_{{H1}}^{{o}},E_{{H12}}^{{o}}, I_{{H12}}^{{o}},E_{{H2}}^{{o}},I_{{H2}}^{{o}},R_{{H2}}^{{o}},E_{{H21}}^{{o}},I_{{H21}}^{{o}},R_{{H22}}^{{o}}, S_{{V}}^{{o}},E_{{V1}}^{{o}},I_{{V1}}^{{o}},E_{{V2}}^{{o}},I_{{V2}}^{{o}} \bigr) \\ & = \biggl(\frac{\varPi _{{H}}}{\mu _{{H}}},0,0,0,0,0,0,0,0,0,0,0, \frac{\varPi _{{V}}}{\mu _{{V}}},0,0,0,0 \biggr). \end{aligned}$$

(3.1)

According to the next generation method and the notations proposed by van den Driessche and Watmough [46], the $12\times 12$ matrices F (containing the new infection terms) and V (containing transfer terms) are given by

$$\begin{aligned} F=\left[ \textstyle\begin{array}{@{}c@{\quad }c@{\quad}c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{}} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & B_{1} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & B_{2} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \eta _{{H1}}A_{{1}} & A_{{1}} & 0 & 0 & 0 & 0 & \eta _{H1}A_{1} & A_{{1}} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \eta _{{H2}}A_{{2}} & A_{{2}} & \eta _{{H2}}A_{{2}} & A_{{2}} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{array}\displaystyle \right] \end{aligned}$$

and

$$\begin{aligned} V=\left[ \textstyle\begin{array}{@{}c@{\quad }c@{\quad}c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{}} K_{{1}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ -\sigma _{{1}} & K_{{3}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & K_{{2}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & -\sigma _{{2}} & K_{{4}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & K_{{2}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & -\sigma _{{2}} & K_{{4}} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & K_{{1}} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & -\sigma _{{1}} & K_{{3}} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & K_{{5}} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -\sigma _{{V1}} & \mu _{{V}} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & K_{{6}} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -\sigma _{{V2}} & \mu _{{V}} \end{array}\displaystyle \right], \end{aligned}$$

respectively, where $A_{{i}}=\frac{\varPi _{{V}}b\beta _{{Vi}}}{\mu _{{V}}}$ and $B_{{i}}=\frac{\varPi _{{H}}b\beta _{{i}}}{\mu _{{H}}}$ for $i=1,2$. The matrix $G:=FV^{{-1}}=[g_{{ij}}]_{{12 \times 12}}$ is a sparse matrix with nonzero elements as specified below:

$$\begin{aligned} \left .\textstyle\begin{array}{l} { g_{{1,9}} = \frac{\varPi _{{H}}b\beta _{{1}}\sigma _{{V1}}}{\mu _{{H}}\mu _{{V}}K_{{5}}}},\qquad {g_{{1,10}} = \frac{\varPi _{{H}}b\beta _{{1}}}{\mu _{{H}}\mu _{{V}}}}, \qquad {g_{{5,11}} = \frac{\varPi _{{H}}b\beta _{{2}}\sigma _{{V2}}}{\mu _{{H}}\mu _{{V}}K_{{6}}}}, \qquad {g_{{5,12}} = \frac{\varPi _{{H}}b\beta _{{2}}}{\mu _{{H}}\mu _{{V}}}} \\ { g_{{9,1}} = g_{{9,7}} = \frac{\varPi _{{V}}b\beta _{{V1}}(\eta _{{H1}}K_{{3}}+\sigma _{{1}})}{\mu _{{V}}K_{{1}}K_{{3}}}}, \qquad {g_{{9,2}} = g_{{9,4}} = \frac{\varPi _{{V}}b\beta _{{V1}}}{\mu _{{V}}K_{{3}}}} \\ {g_{{11,3}} = g_{{11,5}} = \frac{\varPi _{{V}}b\beta _{{V2}}(\eta _{{H2}}K_{{4}}+\sigma _{{2}})}{\mu _{{V}}K_{{2}}K_{{4}}}}, \qquad {g_{{11,4}} = g_{{11,6}} = \frac{\varPi _{{V}}b\beta _{{V2}}}{\mu _{{V}}K_{{4}}}}. \end{array}\displaystyle \right \} \end{aligned}$$

(3.2)

The entries of G in (3.2) are interpreted as the number of secondary infections produced by infected mosquitoes and humans during the course of their infection. The entries $g_{{9,1}}, g_{{9,7}}, g_{{11,3}}$, and $g_{{11,5}}$ are defined as the expected number of infectious humans with serotype-i generated by the infectious period of mosquitoes with serotype-i. The entries $g_{{1,9}}$ and $g_{{5,11}}$ are defined as the expected number of infectious mosquitoes with serotype-i generated by incubation period and infectious period, respectively. It is noted that infected humans produce infected mosquitoes and vise versa. Thus, the entries $g_{{9,2}}, g_{{9,4}}, g_{{11,4}}$, and $g_{{11,6}}$ display that humans cannot infect humans. The entries $g_{{1,10}}$ and $g_{{5,12}}$ show that mosquitoes cannot infect mosquitoes. The other entries are zero, which means no secondary infection in mosquitoes. Let $R_{{0}}^{{n}}=\rho (G)$, where ρ is the spectral radius (dominant eigenvalue in magnitude) of G. It is easy to show that

$$\begin{aligned} R_{{0}}^{{n}}= \operatorname{max} \{ R_{{n1}}, R_{{n2}} \}, \end{aligned}$$

(3.3)

where

$$\begin{aligned} &R_{{ni}} =\sqrt{ {\frac{\varPi _{{H}}\varPi _{{V}}b^{2}\beta _{{i}}\beta _{{Vi}}\sigma _{{Vi}} (\eta _{{Hi}}K_{{i+2}}+\sigma _{{i}})}{\mu _{{H}}\mu _{{V}}^{2}K_{{i}}K_{{i+2}}K_{{i+4}}}}}\quad \text{for } i=1,2. \end{aligned}$$

(3.4)

According to Theorem 2 of [46], the local stability of disease-free equilibrium (DFE) $P^{{o}}_{{n}}$ is established in the following lemma.

Lemma 3.1

The DFE $P^{{o}}_{{n}}$, given in (3.1), of system (2.4)–(2.5) is locally asymptotically stable (LAS) inΩwhenever$R_{{0}}^{{n}}<1$and unstable if$R_{{0}}^{{n}}>1$.

This lemma verifies that $R_{{0}}^{{n}}$ is a threshold value called the basic reproduction number of multi-serotype dengue model and $R_{{ni}}$ is the basic reproduction number associated with infection serotype-i for $i=1,2$. In epidemiology, the reproductive number $R_{{0}}^{{n}}$ represents the number of cases one case generates in a completely susceptible population [47]. It is found that $R_{{0}}^{{n}}$, see (3.3), depends on the values of $R_{{n1}}$ and $R_{{n2}}$, see (3.4). Further, it is found that if $R_{{0}}^{{n}}<1$ whenever $R_{{n1}} <1$ and $R_{{n2}}<1$, a small influx of infected cases into the community would not generate large outbreaks in both primary and secondary infections, and the disease dies out in time because the DFE is LAS as guaranteed by Lemma 3.1. However, to ensure the effective disease control of both primary and secondary dengue infections, it is imperative to show that the DFE is globally asymptotically stable (GAS) in order to make sure that the number of infected cases in the community at steady state are independent of the initial sizes of the sub-populations of a multi-serotype dengue model. The following theorem is established and its proof is given in the Appendix.

Theorem 3.1

The DFE, $P^{{o}}_{{n}}$, of multi-serotype dengue model (2.4)–(2.5) is globally asymptotically stable (GAS) inΩwhenever$R_{{0}}^{{n}}\leq 1$.

Lemma 3.1 and Theorem 3.1 verify that the asymptotic behavior of multi-serotype dengue model (2.4)–(2.5) is determined by the basic reproduction number $R_{{0}}^{{n}}$. As evident in Fig. 2(a)–(m), all profile solutions of multi-serotype dengue model (2.4)–(2.5) always converge to $P^{{o}}_{{n}}$ for all initial sizes of sub-populations used whenever $R_{{0}}^{{n}} \leq 1$ in the sense that the numbers of infectious human and mosquito populations decline to zero, see Fig. 2(c),(d),(e),(f), and (j)–(m). These study results suggest that when $R_{{0}}^{{n}} \leq 1$ the disease will be eradicated from the population irrespective of the initial sizes of sub-populations, which should be the great public health interest.

Figure 2

Time series plots of the multi-serotype dengue model (2.4)–(2.5) with different initial sizes of the sub-populations. The parameter values are given in Table 2 at which $(R_{n1}=0.9747, R_{n2}=0.9851)$ and $R_{{0}}^{{n}}=\max \{R_{n1}, R_{n2}\} = R_{n2} < 1$. (a)–(h) Showing the dynamics of human sub-populations. (i)–(m) Showing the dynamics of mosquito sub-populations

3.2 Existence of endemic equilibrium

Let

$$\begin{aligned} P^{{**}}_{{n}}= {}&\bigl(S_{{H}}^{{**}},E_{{H1}}^{{**}},I_{{H1}}^{{**}},R_{{H1}}^{{**}},E_{{H12}}^{{**}}, I_{{H12}}^{{**}},E_{{H2}}^{{**}},I_{{H2}}^{{**}}, R_{{H2}}^{{**}},E_{{H21}}^{{**}},I_{{H21}}^{{**}}, \\ &R_{{H22}}^{{**}}, S_{{V}}^{{**}},E_{{V1}}^{{**}},I_{{V1}}^{{**}},E_{{V2}}^{{**}},I_{{V2}}^{{**}} \bigr) \end{aligned}$$

(3.5)

be an arbitrary endemic equilibrium. Solving the endemic (positive) equilibrium of multi-serotype dengue model (2.4)–(2.5) in terms of state variables at steady state is laborious owing to its high dimensionality. For simplicity, we start from defining the force of infection with serotype-i of human and the force of infection with serotype-i of mosquito at steady state as follows: for $i,j=1,2$, $(i \neq j)$,

$$\begin{aligned} \varphi _{{i}}^{{**}}= \frac{b\beta _{{i}}I_{{Vi}}^{{**}}}{1+\alpha _{{Vi}}I_{{Vi}}^{{**}}}\quad \text{and}\quad \varphi _{{Vi}}^{{**}}= \frac{b\beta _{{Vi}}(\eta _{{Hi}}(E_{{Hi}}^{{**}}+E_{{Hij}}^{{**}})+I_{{Hi}}^{{**}}+I_{{Hji}}^{{**}})}{1+\alpha _{{Hi}}(I_{{Hi}}^{{**}}+I_{{Hji}}^{{**}})}, \end{aligned}$$

(3.6)

respectively, where ^∗∗ represents the component of the positive equilibrium at steady state. Then all state variables of system (2.4)–(2.5) at steady state can be expressed in terms of $\varphi _{{i}}^{{**}}$ and $\varphi _{{Vi}}^{{**}}$ as given in the following: for $i,j =1,2$$(i \neq j)$,

$$\begin{aligned} \left . \textstyle\begin{array}{l} {S_{{H}}^{{**}}=\frac{\varPi _{{H}}}{G_{\mu \varphi }^{**}}},\qquad E_{{Hi}}^{{**}}= {\frac{\varPi _{{H}}\varphi _{{i}}^{{**}}}{K_{{i}}G_{\mu \varphi }^{**}}},\qquad {I_{{Hi}}^{{**}}=\frac{\varPi _{{H}}\sigma _{{i}}\varphi _{{i}}^{{**}}}{K_{{i}}K_{{i+2}}G_{\mu \varphi }^{**}}},\\ R_{{Hi}}^{{**}}= {\frac{\varPi _{{H}}\sigma _{{i}}\gamma _{{i}}\varphi _{{i}}^{{**}}}{K_{{i}}K_{{i+2}}G_{\mu \varphi \lambda }^{**} G_{\mu \varphi }^{**}}}, \\ {E_{{Hij}}^{{**}}=\frac{\varPi _{{H}}\sigma _{{i}}\gamma _{{i}}\lambda \varphi _{{i}}^{{**}}\varphi _{{j}}^{{**}}}{K_{{i}}K_{{j}}K_{{i+2}} G_{\mu \varphi \lambda }^{**} G_{\mu \varphi }^{**}},\qquad} {I_{{Hij}}^{{**}}=\frac{\varPi _{{H}}\sigma _{{i}}\sigma _{{j}}\gamma _{{i}}\lambda \varphi _{{i}}^{{**}}\varphi _{{j}}^{{**}}}{K_{{i}}K_{{j}}K_{{i+2}}K_{{j+2}} G_{\mu \varphi \lambda }^{**} G_{\mu \varphi }^{**}},} \\ {R_{{H22}}^{{**}}=\sum_{i=1}^{2}\sum_{\substack{j=1 \\ i\neq j}}^{2}\frac{\gamma _{{i}}I_{{Hji}}}{\mu _{{H}}}},\qquad {S_{{V}}^{{**}}=\frac{\varPi _{{V}}}{G_{\mu \varphi _{{V}} }^{**}},\qquad E_{{Vi}}^{{*}}=\frac{\varPi _{{V}}\varphi _{{Vi}}^{{**}}}{K_{{i+4}}G_{\mu \varphi _{{V}} }^{**}}},\\ {I_{{Vi}}^{{**}}=\frac{\varPi _{{V}}\sigma _{{Vi}}\varphi _{{Vi}}^{{**}}}{\mu _{{V}}K_{{i+4}} G_{\mu \varphi _{{V}} }^{**}}}, \end{array}\displaystyle \right \} \end{aligned}$$

(3.7)

where $G_{\mu \varphi \lambda }^{**} = \lambda \varphi _{{j}}^{{**}}+\mu _{{H}}$, $G_{\mu \varphi }^{**} = \varphi _{{i}}^{**}+\varphi _{{j}}^{**}+ \mu _{{H}}$ and $G_{\mu \varphi _{{V}} }^{**} = \varphi _{{Vi}}^{{**}}+\varphi _{{Vj}}^{{**}}+ \mu _{{V}}$.

Substituting expressions (3.7) into (3.6) and after simplifying, we get the following system:

$$\begin{aligned} & \left .\textstyle\begin{array}{l} \varphi _{{1}}^{**} = {\frac{\varPi _{{V}}b\beta _{{1}}\varphi _{{V1}}^{**}}{\mu _{{V}}K_{{5}}(\varphi _{{V1}}^{**}+\varphi _{{V2}}^{**}+\mu _{{V}})+\alpha _{{V1}}\varPi _{{V}}\sigma _{{V1}}\varphi _{{V1}}^{**}}}, \\ \varphi _{{2}}^{**} = {\frac{\varPi _{{V}}b\beta _{{2}}\varphi _{{V2}}^{**}}{\mu _{{V}}K_{{6}}(\varphi _{{V1}}^{**}+\varphi _{{V2}}^{**}+\mu _{{V}})+\alpha _{{V2}}\varPi _{{V}}\sigma _{{V2}}\varphi _{{V2}}^{**}}}, \\ \varphi _{{V1}}^{{**}} = {\frac{\varPi _{{H}}b\beta _{{V1}}(\eta _{{H1}}K_{{2}}+\sigma _{{1}}) (K_{{2}}K_{{4}} G_{\lambda _{1}}^{**}+\lambda _{{1}}\gamma _{{2}}\sigma _{{2}}\varphi _{{2}}^{**}])\varphi _{{1}}^{**}}{K_{{1}}K_{{2}}K_{{3}}K_{{4}}G_{\lambda _{1}}^{**}G^{**} +\varPi _{{H}}\alpha _{{H1}}\sigma _{{1}}(K_{{2}}K_{{4}}G_{\lambda }^{**}+\lambda _{{1}}\gamma _{{2}}\sigma _{{2}}\varphi _{{2}}^{**})\varphi _{{1}}^{**}}}, \\ \varphi _{{V2}}^{{**}} = {\frac{\varPi _{{H}}b\beta _{{V2}}(\eta _{{H2}}K_{{3}}+\sigma _{{2}}) (K_{{1}}K_{{3}}G_{\lambda _{2}}^{**}+\lambda \gamma _{{1}}\sigma _{{1}}\varphi _{{1}}^{**})\varphi _{{2}}^{**}}{K_{{1}}K_{{2}}K_{{3}}K_{{4}}G_{\lambda _{2}}^{**}G^{**} +\varPi _{{H}}\alpha _{{H2}}\sigma _{{2}}(K_{{1}}K_{{3}}G_{\lambda _{2}}^{**}+\lambda _{{2}}\gamma _{{1}}\sigma _{{1}}\varphi _{{1}}^{**})\varphi _{{2}}^{**}}}, \end{array}\displaystyle \right \} \end{aligned}$$

(3.8)

where $G^{**} = \varphi _{{1}}^{**}+\varphi _{{2}}^{**}+\mu _{{H}}$, $G_{\lambda _{1}}^{**} =\lambda \varphi _{{1}}^{**}+\mu _{{H}}$, and $G_{\lambda _{2}}^{**} =\lambda \varphi _{{2}}^{**}+\mu _{{H}}$.

In the case when $\varphi _{{1}}^{**}= \varphi _{{2}}^{**} = 0$, it follows from (3.8) that $\varphi _{{V1}}^{{**}}=\varphi _{{V2}}^{{**}}=0$. Substituting these values into expressions (3.7), we get the disease-free equilibrium.

In the case when $\varphi _{{i}}^{**},\varphi _{{j}}^{**} > 0$ for $i,j=1,2$ ($i \neq j$), system (3.8) is rearranged to the system

$$\begin{aligned} &\biggl( \mu _{{V}} K_{{i+4}} + \varPi _{{V}} \sigma _{{Vi}} \biggl( \alpha _{{Vi}}-\frac{b\beta _{{i}}}{\varphi _{{i}}^{**}} \biggr) \biggr) \varphi _{{Vi}}^{**} +\mu _{{V}}K_{{i+4}}\bigl( \varphi _{{Vj}}^{**}+ \mu _{{V}} \bigr) =0, \end{aligned}$$

(3.9)

$$\begin{aligned} &R_{{ni}}^{{2}} = R_{{nj}}^{{2}} \biggl[ \frac{ (1+{\frac{\lambda \varphi _{{i}}^{**}}{\lambda \varphi _{{j}}^{**}+\mu _{{H}}}} M_{{i}}R_{{ni}}^{{2}} ) (1-{\frac{\alpha _{{Vj}}\varphi _{{j}}^{**}}{b\beta _{{j}}}} ) P_{{i}}}{ (1+{\frac{\lambda \varphi _{{j}}^{**}}{\lambda \varphi _{{i}}^{**}+\mu _{{H}}}}M_{{j}}R_{{nj}}^{{2}} ) (1-\frac{\alpha _{{Vi}}}{b\beta _{{i}}}\varphi _{{i}}^{**} ) P_{{j}}} \biggr], \end{aligned}$$

(3.10)

where

$$\begin{aligned} \begin{aligned}& P_{{i}} = \varphi _{{i}}^{**}+\varphi _{{j}}^{**}+\mu _{{H}} + R_{{ni}}^{{2}}M_{{i+2}} \varphi _{{i}}^{**} \biggl( 1+ \frac{\lambda \varphi _{{j}}^{**}M_{{j}}R_{{nj}}^{{2}}}{\lambda \varphi _{{i}}^{**}+\mu _{{H}}} \biggr), \\ &M_{i}= {\frac{\mu _{{H}}\mu _{{V}}^{{2}}\sigma _{{i}}\gamma _{{i}}K_{{i+4}}}{\varPi _{{H}} \varPi _{{V}}b^{{2}}\beta _{{i}}\beta _{{Vi}}\sigma _{{Vi}}(\eta _{{Hi}}K_{{i+2}}+\sigma _{{i}})}},\qquad M_{{i+2}}=\frac{\varPi _{{H}}\alpha _{{Hi}}}{\gamma _{{i}}}M_{{i}}. \end{aligned} \end{aligned}$$

(3.11)

Clearly, system (3.9)–(3.10) is consistent if and only if $R_{{ni}} /R_{{nj}} > 1$, which implies that system (3.8) is consistent whenever $R_{{ni}}\geq R_{{nj}} > 1$. Therefore, the endemic equilibrium $P^{{**}}_{{n}}$ given in (3.5) of system (2.4)–(2.5) is obtained by solving the fixed point problem (3.8) for positive $\varphi _{i}^{**}$ and $\varphi _{{Vi}}^{{**}}$, $i=1,2$, and then substituting the results obtained into (3.7). Further, the existence of endemic equilibrium of system (2.4)–(2.5) in the primary infection and secondary infection is investigated based on the cross immunity when $0 \leq \lambda <1$ as follows.

Case i: $\lambda =0$

In this case, system (3.9)–(3.10) gives the endemic equilibrium which has only the primary infection whenever $R_{{0}}^{{n}}= \operatorname{max} \{ R_{{n1}}, R_{{n2}} \} > 1$, that is, $R_{{n1}} = R_{{n2}}>1$. This equilibrium is called the co-existence primary infection equilibrium and is given below.

Theorem 3.2

If$\lambda =0$, the multi-serotype dengue model has the co-existence primary infection equilibrium, denoted by$P_{{H}}^{{*}}$:

$$\begin{aligned} P_{{H}}^{{*}}= \bigl(S_{{H}}^{{**}},E_{{H1}}^{{**}},I_{{H1}}^{{**}},R_{{H1}}^{{**}},0, 0,E_{{H2}}^{{**}},I_{{H2}}^{{**}}, R_{{H2}}^{{**}},0,0,0, S_{{V}}^{{**}},E_{{V1}}^{{**}},I_{{V1}}^{{**}},E_{{V2}}^{{**}},I_{{V2}}^{{**}} \bigr), \end{aligned}$$

whenever$R_{{n1}} = R_{{n2}}>1$.

This theorem indicates that when $\lambda =0$, the infection with one serotype produces complete immunity against the other serotypes in the sense that there is no reinfection dengue.

Case ii: $0<\lambda <1$

In the case when $R_{{0}}^{{n}}= \operatorname{max} \{ R_{{n1}}, R_{{n2}} \} >1$, that is, $R_{{ni}} \geq R_{{nj}}>1$, system (3.9)–(3.10) has the endemic equilibrium which is called the co-existence secondary infection equilibrium in the sense that the second infection occurs following the primary infection after some time interval. Thus, the following result is claimed.

Theorem 3.3

If$0<\lambda <1$, the multi-serotype dengue model has the co-existence secondary infection equilibrium, denoted by$P_{{n}}^{{**}}$, whenever$R_{{ni}} \geq R_{{nj}}>1$, $i \neq j$.

This theorem indicates that the infection with one serotype produces partial immunity against another serotypes when $0<\lambda <1$.

3.3 Local stabilities of co-existence equilibria

The local stabilities of co-existence equilibria of the multi-serotype dengue model (see Theorem 3.2 or Theorem 3.3) are analyzed based on the use of the center manifold theory [42]. To simplify, the state variables are changed first of all. Let $S_{{H}}=x_{{1}}, E_{{H1}}=x_{{2}}, I_{{H1}}=x_{{3}}, R_{{H1}}=x_{{4}},E_{{H12}}=x_{{5}}, I_{{H12}}=x_{{6}}, E_{{H2}}=x_{{7}}, I_{{H2}}=x_{{8}}, R_{{H2}}=x_{{9}}, E_{{H21}}=x_{{10}}, I_{{H21}}=x_{{11}}, R_{{H22}}=x_{{12}},S_{{V}}=x_{{13}}, E_{{V1}}=x_{{14}}, I_{{V1}}=x_{{15}}, E_{{V2}}=x_{{16}}$, and $I_{{V2}}=x_{{17}}$, so that $N_{{H}}=\sum_{n=1}^{12}{x_{{n}}}$ and $N_{{V}}=\sum_{n=1}^{5}{x_{{n+12}}}$, respectively. Further, by using the vector notations $x= (x_{{1}},x_{{2}},x_{{3}},\ldots,x_{{17}} )^{T}$ and $F(x)= (f_{{1}}(x),f_{{2}}(x),\ldots,f_{{17}}(x) )^{T}$, multi-serotype dengue model (2.4)–(2.5) is written in the form

$$\begin{aligned} {\frac{dx}{dt}}=F(x), \end{aligned}$$

(3.12)

where

$$\begin{aligned} \left .\textstyle\begin{array}{l} {f_{{1}}}= {\varPi _{{H}}-(\varphi _{{1}}+\varphi _{{2}})x_{{1}}-\mu _{{H}}x_{{1}}}, \qquad {f_{{2}}}= {\varphi _{{1}}x_{{1}}-(\sigma _{{1}}+\mu _{{H}})x_{{2}}}, \\ {f_{{3}}}= {\sigma _{{1}}x_{{2}}-(\gamma _{{1}}+\delta _{{1}}+\mu _{{H}})x_{{3}}}, \qquad {f_{{4}}}= {\gamma _{{1}}x_{{3}}-(\lambda \varphi _{{2}}+\mu _{{H}})x_{{4}}}, \\ {f_{{5}}}= {\lambda \varphi _{{2}}x_{{4}}-(\sigma _{{2}}+\mu _{{H}})x_{{5}}} \qquad {f_{{6}}}= {\sigma _{{2}}x_{{5}}-(\gamma _{{2}}+\delta _{{2}}+\mu _{{H}})x_{{6}}}, \\ {f_{{7}}}= {\varphi _{{2}}x_{{1}}-(\sigma _{{2}}+\mu _{{H}})x_{{7}}}, \qquad {f_{{8}}}= {\sigma _{{2}}x_{{7}}-(\gamma _{{2}}+\delta _{{2}}+\mu _{{H}})x_{{8}}}, \\ {f_{{9}}}= {\gamma _{{2}}x_{{8}}-(\lambda \varphi _{{1}}+\mu _{{H}})x_{{9}}}, \qquad {f_{{10}}}= {\lambda \varphi _{{1}}x_{{9}}-(\sigma _{{1}}+\mu _{{H}})x_{{10}}}, \\ {f_{{11}}}= {\sigma _{{1}}x_{{10}}-(\gamma _{{1}}+\delta _{{1}}+\mu _{{H}})x_{{11}}}, \qquad {f_{{12}}}= {\gamma _{{1}}x_{{11}}+\gamma _{{2}}x_{{6}}-\mu _{{H}}x_{{12}}}, \\ {f_{{13}}}= {\varPi _{{V}}-(\varphi _{{V1}}+\varphi _{{V2}})x_{{13}}-\mu _{{V}}x_{{13}}}, \qquad {f_{{14}}}= {\varphi _{{V1}}x_{{13}}-(\sigma _{{V1}}+\mu _{{V}})x_{{14}}}, \\ {f_{{15}}}= {\sigma _{{V1}}x_{{14}}-\mu _{{V}}x_{{15}}}, \qquad {f_{{16}}}= {\varphi _{{V2}}x_{{13}}-(\sigma _{{V2}}+\mu _{{V}})x_{{16}}}, \\ {f_{{17}}}= {\sigma _{{V2}}x_{{16}}-\mu _{{V}}x_{{17}}}, \qquad \varphi _{{1}}= {\frac{b\beta _{{1}}x_{{15}}}{1+\alpha _{{V1}}x_{{15}}}}, \varphi _{{2}}= {\frac{b\beta _{{2}}x_{{17}}}{1+\alpha _{{V2}}x_{{17}}}}, \\ \varphi _{{V1}}= {\frac{b\beta _{{V1}} (\eta _{{H1}}(x_{{2}}+x_{{10}})+x_{{3}}+x_{{11}} )}{1+\alpha _{{H1}}(x_{{3}}+x_{{11}})}}, \qquad\varphi _{{V2}}= {\frac{b\beta _{{V2}} (\eta _{{H2}}(x_{{5}}+x_{{7}})+x_{{6}}+x_{{8}} )}{1+\alpha _{{H2}}(x_{{6}}+x_{{8}})}}. \end{array}\displaystyle \right \} \end{aligned}$$

(3.13)

Since system (3.12) is identical to system (2.4)–(2.5), the disease-free equilibrium $P^{{o}}_{{n}}$ and the reproduction number $R_{{0}}^{{n}}$ of system (3.12) are given in (3.1) and (3.3), respectively. Consider $R_{{0}}^{{n}}= \operatorname{max}\{R_{{n1}}, R_{{n2}}\} =1$, then

$$\begin{aligned} \beta _{{2}}=\beta _{{2}}^{{**}}= \frac{\mu _{{H}}\mu _{{V}}^{{2}}K_{{2}}K_{{4}}K_{{6}}}{\varPi _{{H}}\varPi _{{V}}b^{{2}}\beta _{{V2}}\sigma _{{V2}}(\eta _{{H2}}K_{{4}}+\sigma _{{2}})} \end{aligned}$$

is chosen to be a bifurcation parameter with

$$\begin{aligned} \beta _{{1}}= \frac{\beta ^{{**}}_{{2}}\beta _{{V2}}\sigma _{{V2}}(\eta _{{H2}}K_{{4}}+\sigma _{{2}})K_{{1}}K_{{3}}K_{{5}}}{\beta _{{V1}}\sigma _{{V1}}(\eta _{{H1}}K_{{3}}+\sigma _{{1}})K_{{2}}K_{{4}}K_{{6}}}. \end{aligned}$$

The Jacobian of system (3.12) evaluated at $P^{{o}}_{{n}}$ with $\beta _{{2}}=\beta _{{2}}^{{**}}$ is given by

$$\begin{aligned} J\bigl(\beta _{{2}}^{{**}}\bigr) = \begin{bmatrix} J_{{1}} & 0_{{6\times 5}} & J_{{2}} &\\ 0_{{5\times 6}} & J_{{3}} & J_{{4}} &\\ J_{{5}} & J_{{6}} & J_{{7}} \end{bmatrix}, \end{aligned}$$

where $0_{{m \times n}}$ is a zero matrix of m rows by n column and the matrices $J_{{i}}, i=1,\ldots,7$, are given by

$$\begin{aligned} & J_{{1}} = \begin{bmatrix} -\mu _{{H}} & 0 & 0 & 0 & 0 & 0 \\ 0 & -K_{{1}} & 0 & 0 & 0 & 0 \\ 0 & \sigma _{{1}} & -K_{{3}} & 0 & 0 & 0 \\ 0 & 0 & \gamma _{{1}} & -\mu _{{H}} & 0 & 0 \\ 0 & 0 & 0 & 0 & -K_{{2}} & 0 \\ 0 & 0 & 0 & 0 & \sigma _{{2}} & -K_{{4}} \end{bmatrix},\\ &J_{{2}} = \begin{bmatrix} 0 & 0 & 0 & -C\beta _{{2}}^{{**}} & 0 & -C\beta _{{2}}^{{**}} \\ 0 & 0 & 0 & C\beta _{{2}}^{{**}} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix}, \\ & J_{{3}} = \begin{bmatrix} -K_{{2}} & 0 & 0 & 0 & 0 \\ \sigma _{{2}} & -K_{{4}} & 0 & 0 & 0 \\ 0 & \gamma _{{2}} & -\mu _{{H}} & 0 & 0 \\ 0 & 0 & 0 & -K_{{1}} & 0 \\ 0 & 0 & 0 & \sigma _{{1}} & -K_{{3}} \end{bmatrix},\qquad J_{{4}} = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & C\beta _{{2}}^{{**}} \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix}, \\ & J_{{5}} = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & \gamma _{{2}} \\ 0 & -\eta _{{H1}}A_{{1}} & -A_{{1}} & 0 & -\eta _{{H2}}A_{{2}} & -A_{{2}} \\ 0 & \eta _{{H1}}A_{{1}} & A_{{1}} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \eta _{{H2}}A_{{2}} & A_{{2}} \\ 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix},\\ &J_{{6}} = \begin{bmatrix} 0 & 0 & 0 & 0 & \gamma _{{1}} \\ -\eta _{{H2}}A_{{2}} & -A_{{2}} & 0 & -\eta _{{H1}}A_{{1}} & -A_{{1}} \\ 0 & 0 & 0 & \eta _{{H1}}A_{{1}}& A_{{1}} \\ 0 & 0 & 0 & 0 & 0 \\ \eta _{{H2}}A_{{2}} & A_{{2}} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{bmatrix},\\& J_{{7}} = \begin{bmatrix} -\mu _{{H}} & 0 & 0 & 0 & 0 & 0 \\ 0 & -\mu _{{V}} & 0 & 0 & 0 & 0 \\ 0 & 0 & -K_{{5}} & 0 & 0 & 0 \\ 0 & 0 & \sigma _{{V1}} & -\mu _{{V}} & 0 & 0 \\ 0 & 0 & 0 & 0 & -K_{{6}} & 0 \\ 0 & 0 & 0 & 0 & \sigma _{{V2}} & -\mu _{{V}} \end{bmatrix}, \end{aligned}$$

with $A_{{i}}={\frac{\varPi _{{V}}b\beta _{{Vi}}}{\mu _{{V}}}}$, $i=1,2$, and $C= {\frac{\varPi _{{H}}b\beta _{{V2}}\sigma _{{V2}}(\eta _{{H2}}K_{{4}}+\sigma _{{2}})K_{{1}}K_{{3}}K_{{5}}}{\mu _{{H}}\beta _{{V1}}\sigma _{{V1}}(\eta _{{H1}}K_{{3}}+\sigma 1)K_{{2}}K_{{4}}K_{{6}}}}$, respectively. It can be shown that $J(\beta _{{2}}^{{**}})$ has at least one eigenvalue with zero real part. Hence, the center manifold theory [42] can be used to analyze the dynamics of (3.12) with $\beta _{{2}}=\beta _{{2}}^{{**}}$.

Next, the right eigenvector of $J(\beta _{{2}}^{{**}})$ associated with the zero eigenvalue at $\beta _{{2}}=\beta _{{2}}^{{**}}$ is given by $w = [w_{{1}},w_{{2}},w_{{3}},\ldots,w_{{17}}]^{{T}}$, where

$$\begin{aligned} \left .\textstyle\begin{array}{l} { w_{{1}}= -\frac{K_{{1}}w_{{2}}+K_{{2}}w_{{7}}}{\mu _{{H}}},\qquad w_{{3}} = \frac{\sigma _{{1}}w_{{2}}}{K_{{3}}},\qquad w_{{4}} = \frac{\gamma _{{1}}\sigma _{{1}}w_{{2}}}{\mu _{{H}}K_{{3}}},\qquad w_{{2}} = w_{{2}}>0}, \\ w_{{5}} = w_{{6}}= 0,\qquad w_{{8}}= \frac{\sigma _{{2}} w_{{7}}}{K_{{4}}}, \qquad w_{{9}}=\frac{\gamma _{{2}} \sigma _{{2}} w_{{7}}}{\mu _{{H}}K_{{4}}},\\ w_{{10}} = w_{{11}}=w_{{12}} = 0, \qquad w_{{7}} = w_{{7}} > 0, \\ {w_{{13}} = -\frac{b\varPi _{{V}}\beta _{{V1}}K_{{4}} (\eta _{{H1}}K_{{3}}+\sigma _{{1}} )}{\mu _{{V}}^{{2}}K_{{3}}K_{{4}}}w_{{2}} -\frac{b\varPi _{{V}}\beta _{{V2}}K_{{3}}(\eta _{{H2}}K_{{4}}+\sigma _{{2}})}{\mu _{{V}}^{{2}}K_{{3}}K_{{4}}}w_{{7}}}, \\ { w_{{14}} = \frac{b\varPi _{{V}}\beta _{{V1}}(\eta _{{H1}}K_{{3}}+\sigma _{{1}})}{\mu _{{V}} K_{{3}}K_{{5}}}w_{{2}},\qquad {} w_{{15}} = \frac{b\varPi _{{V}}\beta _{{V1}}\sigma _{{V1}}(\eta _{{H1}}K_{{3}}+\sigma _{{1}})}{\mu _{{V}}^{{2}}K_{{3}}K_{{5}}}w_{{2}}}, \\ { w_{{16}} = \frac{b\varPi _{{V}}\beta _{{V2}}(\eta _{{H2}}K_{{4}}+\sigma _{{2}})}{\mu _{{V}}K_{{4}}K_{{6}}}w_{{7}}, \qquad w_{{17}} = \frac{b\varPi _{{V}}\beta _{{V2}}\sigma _{{V2}}(\eta _{{H2}}K_{{4}}+\sigma _{{2}})}{\mu _{{V}}^{{2}}K_{{4}}K_{{6}}}w_{{7}}}. \end{array}\displaystyle \right \} \end{aligned}$$

(3.14)

Similarly, the components of the left eigenvector of $J(\beta _{{2}}^{{**}})$ associated with the zero eigenvalues at $\beta _{{2}}=\beta _{{2}}^{{**}}$, denoted by $v=[v_{{1}},v_{{2}},v_{{3}},\ldots,v_{{17}}]$, are given by

$$\begin{aligned} \left .\textstyle\begin{array}{l} v_{{1}}=v_{{4}}=v_{{9}}= v_{{12}} = v_{{13}} = 0,\qquad v_{{14}} = v_{{14}}>0, \qquad v_{{16}} = v_{{16}}>0, \\ {v_{{2}}=\frac{b\varPi _{{V}}\beta _{{V1}}(\eta _{{H1}}K_{{3}}+\sigma _{{1}})}{\mu _{{V}}K_{{1}}K_{{3}}}v_{{14}}, \qquad v_{{3}}=\frac{b\varPi _{{V}}\beta _{{V1}}}{\mu _{{V}}K_{{3}}}v_{{14}},\qquad v_{{5}} = \frac{b\varPi _{{V}}\beta _{{V2}}(\eta _{{H2}}K_{{4}}+\sigma _{{2}})}{\mu _{{V}}K_{{2}}K_{{4}}}v_{{16}}}, \\ {v_{{6}} = \frac{b\varPi _{{V}}\beta _{{V2}}}{\mu _{{V}}K_{{4}}}v_{{16}},\qquad v_{{7}} = \frac{b\varPi _{{V}}\beta _{{V2}}(\eta _{{H2}}K_{{4}}+\sigma _{{2}})}{\mu _{{V}}K_{{2}}K_{{4}}}v_{{16}},\qquad v_{{8}} = \frac{b\varPi _{{V}}\beta _{{V2}}}{\mu _{{V}}K_{{4}}}v_{{16}}}, \\ v_{{10}}=\frac{b\varPi _{{V}}\beta _{{V1}}(\eta _{{H1}}K_{{3}}+\sigma _{{1}})}{\mu _{{V}}K_{{1}}K_{{3}}}v_{{14}},\qquad v_{{11}} = \frac{b\varPi _{{V}}\beta _{{V1}}}{\mu _{{V}}K_{{3}}}v_{{14}},\\ v_{{15}} =\frac{K_{{5}}}{\sigma _{{V1}}}v_{{14}} \quad \text{and}\quad v_{{17}} = \frac{K_{{6}}}{\sigma _{{V2}}}, \end{array}\displaystyle \right \} \end{aligned}$$

(3.15)

respectively.

Computing the associated nonzero partial derivatives of $F(x)$ at $P^{{o}}_{{n}}$ and using the expression in (3.14)–(3.15), the coefficients a and b defined in Theorem 4.1 [48] are given by

$$\begin{aligned} a ={}& \sum_{{k,n,m=1}}^{{17}}v_{{k}}w_{{n}}w_{{m}} \frac{\partial ^{{2}}f_{{k}}}{\partial x_{{n}}\partial x_{{m}}}(0,0) \\ ={}& {-} \frac{2b\varPi _{{V}}\beta _{{V1}}K_{{7}}w_{{2}}w_{{14}}}{\mu _{{V}}\varPi _{{H}}K_{{3}}K_{{4}}} \bigl( K_{{1}}K_{{4}} w_{{2}}+ ( K_{{2}}K_{{4}}-\gamma _{{2}} \sigma _{{2}} \lambda ) w_{{7}} \bigr) \\ &{}- \frac{2b\varPi _{{V}}\beta _{{V2}}K_{{8}}w_{{7}}w_{{16}}}{\mu _{{V}}\varPi _{{H}}K_{{3}}K_{{4}}} \bigl( K_{{1}}K_{{4}} w_{{2}}+ ( K_{{2}}K_{{4}}-\gamma _{{2}} \sigma _{{2}} \lambda ) w_{{7}} \bigr) \\ &{}- \frac{2b^{{2}}\varPi _{{V}}^{{2}}}{\mu _{{V}}^{{3}}K_{{3}}^{{2}}K_{{4}}^{{2}}K_{{5}}K_{{6}}} \bigl( \beta _{{V1}}^{{2}} \sigma _{{V1}} \alpha _{{V1}} K_{{7}}^{{2}} K_{{4}}^{{2}}K_{{6}} w_{{2}}^{{2}}v_{{14}} +\beta _{{V2}}^{{2}} \sigma _{{V2}} \alpha _{{V2}} K_{{3}}^{{2}} K_{{5}} K_{{8}}^{{2}} w_{{7}}^{{2}}v_{{16}} \bigr) \\ &{}- \frac{2b^{{2}}\varPi _{{V}} (\beta _{{V1}}K_{{4}}K_{{7}}w_{{2}}+\beta _{{V2}}K_{{3}}K_{{8}}w_{{7}} )}{\mu _{{V}}^{{2}}K_{{3}}^{{2}}K_{{4}}^{{2}}} ( \beta _{{V1}} K_{{4}} K_{{7}} w_{{2}}v_{{14}} +\beta _{{V2}} K_{{3}} K_{{8}} w_{{7}}v_{{16}} ) \end{aligned}$$

(3.16)

and

$$\begin{aligned} b &= \sum_{{k,n=1}}^{{17}}v_{{k}}w_{{n}} \frac{\partial ^{2}f_{{k}}}{\partial x_{{n}}\partial \beta _{{2}}^{{**}}}(0,0) \\ & = \frac{\varPi _{{H}}\varPi _{{V}}^{{2}}b^{{3}}\beta _{{V2}}\sigma _{{V2}}K_{{8}} ( \beta _{{V1}}K_{{4}}K_{{7}}w_{{2}}v_{{14}} + \beta _{{V2}}K_{{3}}K_{{8}}w_{{7}}v_{{16}} )}{\mu _{{V}}^{{3}}\mu _{{H}}K_{{2}}K_{{3}}K_{{4}}^{{2}}K_{{6}}}, \end{aligned}$$

(3.17)

where $K_{{7}}=\eta _{{H1}}K_{{3}}+\sigma _{{1}}$ and $K_{{8}}=\eta _{{H2}}K_{{4}}+\sigma _{{2}}$.

Clearly, the coefficient b is always positive and the coefficient a is positive whenever $0 \leq \lambda <1$. Similar results are obtained when the bifurcation is chosen to be $\beta _{{1}}=\beta _{{1}}^{{**}}$. According to Theorems 3.2–3.3 and Theorem 4.1 of [48], the following lemmas are established.

Lemma 3.2

If$\lambda =0$, the co-existence primary infection equilibrium$P_{{H}}^{{*}}$of the multi-serotype dengue model exists and is locally asymptotically stable whenever$R_{{0}}^{{n}} > 1$ (that is, $R_{{ni}} = R_{{nj}}>1$ for $i,j=1,2$ and $i \neq j$) and is close to 1.

Lemma 3.3

If$0 < \lambda <1$, the co-existence secondary infection equilibrium$P^{{**}}_{{n}}$of multi-serotype dengue model exists and is locally asymptotically stable whenever$R_{{0}}^{{n}} > 1$ (that is, $R_{{ni}} \geq R_{{nj}}>1$ for $i,j=1,2$ and $i \neq j$) and is close to 1.

4 Numerical simulations

4.1 Estimation of model parameters

The Ministry of Public Health Thailand reported the dengue outbreak during 2014–2018 and indicated that dengue is epidemic every year, especially the number of dengue cases increased in the year 2018 [30]. Thus, the multi-serotype dengue model is applied for predicting the number of infectious populations and computing the cumulative number of dengue cases by solving the following differential equation [49]:

$$\begin{aligned} \frac{dC}{at}=\kappa _{{1}}(I_{{H1}}+I_{{H21}})+ \kappa _{{2}}(I_{{H2}}+I_{{H12}}), \end{aligned}$$

(4.1)

where C denotes the cumulative number of dengue cases, $\kappa _{{1}}$ and $\kappa _{{2}}$ are the rates of progression from infective to diagnosed. The parameters used in simulation are $\lambda = 0.0111$, $\alpha _{{H}} = 8.648 \times 10^{-5}$, $\alpha _{{V}} =1.3404 \times 10^{-6}$, $\kappa _{{1}} =0.0361$, $\kappa _{{2}} =0.0283$, and the other parameter values are given in Table 2. The values of λ, $\alpha _{{H}}$, and $\alpha _{{V}}$ [30] correspond to three months of cross immunity, the prevalence rates per one million people, and mosquitoes’ density, respectively. The other model parameters chosen are reliable data corresponding to the dengue transmission in Thailand and those in the literature.

The cumulative numbers of dengue cases produced by the multi-serotype dengue model are compared with the real data in the year 2018 [30], see Fig. 3. The accuracy of the approximated data obtained are tested by using the coefficient of determination denoted by $R^{2}$ [50]. It is found that $R^{{2}} = 0.9997$ verifies that the approximated data are closer to the real data, see Fig. 3. This study indicates that the values of model parameter in Table 2 are appropriated values for the multi-serotype dengue model. These values will be used to investigate the sensitivity of model parameters that impact the transmission of dengue and then to analyze the effects of cross immunity and nonlinear incidence rate on the occurrence of reinfection dengue.

Figure 3

Comparison of the cumulative number of infected number produced by the multi-serotype dengue model and the cumulative number of reported dengue cases in Thailand, 2018

4.2 Sensitivity analysis

The parameters that have high impact on the transmission and should be targeted by interference strategies can be discovered by the sensitivity analysis [51]. When a parameter changes, the relative change in a variable can be measured by the normalized forward sensitivity index of a variable with respect to a parameter. The normalized forward sensitivity index of the reproduction number of a multi-dengue model $R_{{0}}^{{n}}$ with respect to the parameter p is given by

$$\begin{aligned} \varUpsilon ^{{R_{{0}}^{{n}}}}_{{p}} = \frac{\partial R_{{0}}^{{n}}}{\partial p} \times \frac{ p }{ R_{{0}}^{{n}}}. \end{aligned}$$

(4.2)

It is found, from predicted cumulative number of dengue cases, that $R_{{0}}^{{n}} =R_{{2}}^{{n}}>1$. Then the sensitivity indices of $R_{{2}}^{{n}}$ to the twelve different parameters at the baseline parameter values (see Table 2) are shown in Fig. 4. The study results show that the most sensitive parameter for mosquito is the natural death rate $\mu _{{V}}$ followed by the mosquito biting rate b, the probability of disease transmission from infectious humans to susceptible mosquitoes $\beta _{{V1}}$, and the recruitment rate by mosquito birth $\varPi _{{V}}$, respectively. Meanwhile, the most sensitive parameter for human population is the mosquito biting rate b followed by the natural death rate $\mu _{{H}}$, the recruitment rate $\varPi _{{H}}$, the transmission probability rate from infected mosquito to susceptible human $\beta _{{1}}$, the human recovery rate $\gamma _{{1}}$, and the modification parameters $\eta _{{H1}}$, respectively. It is also found that the sensitivity indices of $R_{{0}}^{{n}}$ with respect to $\beta _{{1}}$, $\beta _{{V1}}$, $\varPi _{{H}}$, and $\varPi _{{V}}$ do not depend on any parameter values, which implies that increasing these model parameters will lead to increasing $R_{{0}}^{{n}} =R_{{1}}^{{n}}$ and vice versa. Similarly, the same results are found when $R_{{0}}^{{n}} =R_{{1}}^{{n}}>1$ by changing index 2 to index 1. These study results indicate that the model parameters $\beta _{{V2}}$, b, $\beta _{{2}}$, and $\eta _{{H2}}$ affect the transmission of dengue. On the other hand, the model parameters $\mu _{{V}}$ and $\gamma _{{1}}$ impact its prevention.

Figure 4

Sensitivity indexes of $R_{{0}}^{{n}}=R_{{n2}}^{{n}}$ to the twelve different parameters for dengue serotype-2: $\mu _{{V}}$, b, $\mu _{{H}}$, $\varPi _{{V}}$, $\varPi _{{H}}$, $\beta _{{2}}$, $\beta _{{V2}}$, $\gamma _{{2}}$, $\sigma _{{V2}}$, $\eta _{{H2}}$, $\sigma _{{2}}$, $\delta _{{2}}$ at the baseline parameter values given in Table 2

4.3 Effect of cross immunity and nonlinear incidence rates on secondary infections

It is found that, from (3.3)–(3.4), the reproductive number of the multi-serotype dengue model $R_{{0}}^{{n}}$ does not depend on the cross immunity and nonlinear incidence rates. However, the co-existence equilibria of the multi-serotype dengue model exist whenever the ratio of $R_{{ni}}^{{2}}$ and $R_{{nj}}^{{2}}$ is greater than unity (see system (3.9)–(3.10)) and is guaranteed by Lemmas 3.2–3.3, respectively. To investigate the effect of cross immunity and nonlinear incidence rates on the secondary infections, the relation of $R_{{ni}}^{{2}}$ and $R_{{nj}}^{{2}}$ given in (3.10) is compared with the following functions:

$$\begin{aligned} R_{{ni}}^{{2}} > f_{i}\bigl(R_{{nj}}^{{2}} \bigr) \quad\text{if } R_{{ni}} \geq R_{{nj}} >1, \end{aligned}$$

(4.3)

where

$$\begin{aligned} f_{i}\bigl(R_{{nj}}^{{2}}\bigr):= { \frac{R_{{nj}}^{{2}}}{1+\lambda (1-\alpha _{{V}})(1-\frac{\alpha _{{H}}\sigma _{{i}}}{K_{{i}}K_{{i+2}}})M_{{j}}[R_{{nj}}^{{2}}-1]}}, \end{aligned}$$

(4.4)

and $M_{{j}}= {\frac{\varPi _{{H}}\varPi _{{V}}b^{{2}}\beta _{{j}}\beta _{{Vj}}\sigma _{{Vj}}(\eta _{{Hj}}K_{{j+2}}+\sigma _{{j}})}{\mu _{{H}}\mu _{{V}}^{{2}}\sigma _{{j}}\gamma _{{j}}K_{{j+4}}}}$ for $i,j=1,2$ and $i\neq j$, respectively.

Let $\lambda = \lambda _{{i}}^{*}$, $\alpha _{{H}}=\alpha ^{*}_{{Hi}}$, and $\alpha _{{V}}=\alpha ^{*}_{{Vi}}$ be the bifurcation parameters of $f_{i}$, $i=1,2$. The first and second derivatives of $f_{i}(R_{{nj}}^{{2}})$ with respect to $R_{{nj}}^{{2}}$ are given by

$$\begin{aligned} f_{i}'\bigl(R_{{nj}}^{{2}} \bigr) = { \frac{1-A}{ (1+ A(R_{{nj}}^{{2}}-1) )^{2}}} \quad\text{and}\quad f_{i}'' \bigl(R_{{nj}}^{{2}}\bigr) = { \frac{-2A(1-A)}{ (1+A(R_{{nj}}^{{2}}-1) )^{3}}}, \end{aligned}$$

(4.5)

respectively, where $A_{i}= \lambda (1-\alpha _{{V}}) (1-{ \frac{\alpha _{{H}}\sigma _{{i}}}{K_{{i}}K_{{i+2}}}} )M_{{j}}$ for $j=1,2$ and $i\neq j$. Further, solving $f'_{{i}}(R_{{nj}}^{{2}})=0$ for $\lambda ^{*}_{{i}}$, $\alpha ^{*}_{{Hi}}$, and $\alpha ^{*}_{{Vi}}$, we get

$$\begin{aligned} &\lambda ^{*}_{{i}} = \frac{1}{M_{j}(1-\alpha _{{V}}) (1- {\frac{\alpha _{{H}}\sigma _{{i}}}{K_{{i}}K_{{i+2}}}} )},\qquad \alpha _{{Vi}}^{*} =1- \frac{1}{\lambda M_{{j}} (1- {\frac{\alpha _{{H}}\sigma _{{i}}}{K_{{i}}K_{{i}+2}}} )}, \end{aligned}$$

(4.6)

$$\begin{aligned} & \alpha _{{Hi}}^{*} =\frac{K_{{i}}K_{{i+2}}}{\sigma _{{i}}} \biggl( 1- \frac{1}{\lambda M_{{j}}(1-\alpha _{{V}})} \biggr). \end{aligned}$$

(4.7)

According to (4.4)–(4.7), the following lemma is established.

Lemma 4.1

For all$R_{{ni}}^{{2}}>R_{{nj}}^{{2}}>1$, $i,j=1,2$and$i\neq j$.

• (i) The function$f_{i}(R_{{nj}}^{{2}})$is increasing and concave down if one of the following conditions holds:

  • (a) $\lambda < \lambda _{{i}}^{*}$;

  • (b) $\alpha _{{H}} > \alpha _{{Hi}}^{*}$;

  • (c) $\alpha _{{V}} > \alpha _{{Vi}}^{*}$.

• (ii) The function$f_{i}(R_{{nj}}^{{2}})$is decreasing and concave up if one of the following conditions holds:

  • (a) $\lambda > \lambda _{{i}}^{*}$;

  • (b) $\alpha _{{H}} < \alpha _{{Hi}}^{*}$;

  • (c) $\alpha _{{V}} < \alpha _{{Vi}}^{*}$.

• (iii) The functions$f_{{1}}(R_{{n2}}^{{2}})=f_{{2}}(R_{{n1}}^{{2}})=1$if one of the following conditions holds:

  • (a) $\lambda = \lambda _{{i}}^{*}$;

  • (b) $\alpha _{{H}} = \alpha _{{Hi}}^{*}$;

  • (c) $\alpha _{{V}} = \alpha _{{Vi}}^{*}$.

• (iv) The functions$f_{{1}}(R_{{n2}}^{{2}})=f_{{2}}(R_{{n1}}^{{2}})>1$if one of the following conditions holds:

  • (a) $\lambda = 0$;

  • (b) $\alpha _{{H}} = { \frac{K_{{i}}K_{{i+2}}}{\sigma _{{i}}}}$;

  • (c) $\alpha _{{V}} = 1$.

4.4 Theoretical results

The results in Lemma 4.1 are applied to monitor the effect of cross immunity (λ) and the inhibitory effect for humans ($\alpha _{{H}}$) and mosquitoes ($\alpha _{{V}}$) on the secondary infection as follows. First, the conditions in Lemma 4.1 are tested to determine the area of secondary infection by plotting the functions $f_{{1}}(R_{{n2}}^{{2}})$ versus $R_{{n2}}^{{2}}$ and $f_{{2}}(R_{{n1}}^{{2}})$ versus $R_{{n1}}^{{2}}$. The model parameter values/ranges used for simulation are given in Table 2. The results obtained are shown in Fig. 5. In this figure, the secondary infection is the area bounded by the profiles of $f_{{1}}(R_{{n2}}^{{2}})$ and $f_{{2}}(R_{{n1}}^{{2}})$. The primary infection with dengue serotype-i represents the stability region of co-existence primary infection equilibrium $P_{{H}}^{{*}}$ (in line with Theorem 3.2), and the secondary infection represents the stability region of co-existence secondary infection equilibrium $P_{{n}}^{{*}}$ (in line with Theorem 3.3) whenever $R_{{ni}}^{{2}} \geq R_{{nj}}^{{2}}>1$ for $i,j=1,2$ and $i\neq j$. It is found that the second infection occurs whenever conditions (i)–(iii) in Lemma 4.1 hold, see Fig. 5(b)–(d). On the other hand, the second infection is reduced if the condition (i) holds, see Fig. 5(b), and then it is eliminated if condition (iv) holds, see Fig. 5(a). These study results reveal that the cross immunity and nonlinear incidence rates for human and mosquito are important factors in controlling and preventing the secondary infection.

Figure 5

Showing the region of primary and secondary infections which are separated by the functions $f_{{1}}$ and $f_{{2}}$ given in (4.4) under the condition in Lemma 4.1

4.5 Effect of cross immunity

The effect of cross immunity is investigated with various λ. The values of λ are chosen to be 0.0714, 0.0110, 0.0037, and 0.0014 which correspond to the period of cross immunity as two weeks [52], three months, nine months, and two years [53, 54], respectively. With these parameters and the parameter uses given in Table 2, the secondary infection, which is the area bounded by the profiles of $f_{{1}}(R_{{n2}}^{{2}})$ and $f_{{2}}(R_{{n1}}^{{2}})$, is displayed in Fig. 6. It is found that when $\lambda =0.0714$, the area of secondary infection is widest and would reduce when λ increases to 0.0110, 0.0037, and 0.0014, respectively. This study interprets that the cross immunity affects the decrease in the occurrence of dengue reinfection after getting the primary infection.

Figure 6

Showing the region of secondary infection covered by the functions $f_{{1}}$ and $f_{{2}}$ given in (4.4) with $\lambda =0.0714$ (brown solid), 0.0110 (blue dashed), 0.0037 (red dashdot), and 0.0014 (green solid)

The number of infected populations at steady state is predicted by simulating multi-serotype dengue model (2.4)–(2.5) using the values of model parameters given in Table 2 and various λ. With the parameter values used in Table 2, the reproductive numbers associated with infection dengue serotype-1 and dengue serotype-2 are given by $R_{{n1}}=3.8987$ and $R_{{2}}=3.9402$ so that the reproductive number of multi-serotype dengue model (2.4)–(2.5) is $R_{{0}}^{{n}}= \operatorname{max} \{ R_{{n1}}, R_{{n2}} \} =R_{{2}}=3.9402$. The results obtained are tabulated in Table 3. It is found from Table 3 that when $R_{{n2}} > R_{{n1}} > 1$, increasing the period of cross immunity would decrease the number of infected humans in primary infection, which results in decreasing the number of mosquitoes carrying dengue serotype-1 and serotype-2. However, increasing the period of cross immunity would decrease the number of infected humans in secondary infection, which exchanges the stabilities of two equilibria ($P^{{**}}_{{n}}$ and $P_{{H}}^{{*}}$). This, of course, is in line with Lemmas 3.2–3.3. Further, Figs. 7–8 display the number of infected humans and mosquitoes produced by multi-serotype dengue model (2.4)–(2.5) using the values of model parameters given in Table 2 and various λ ($\lambda = 0.0111$ and $\lambda = 0.0014$). The value of λ corresponds to the period of cross immunity as three months and two years, respectively. The initial conditions in simulation vary, which corresponds to the dynamics of dengue transmission. The numbers of infected humans and mosquitoes for ten years are shown in Figs. 7–8. It is found that if the period of cross immunity is short (approximated three months), the number of infected humans in secondary infection increases and fluctuates, see Fig. 7(c)–(d). On the other hand, increasing the period of cross immunity (approximated two years) would decrease the number of infected humans in secondary infection but it would not die out as shown in Fig. 8(c)–(d). This is a strong evidence that the immunological interactions between serotypes are of central importance in understanding epidemiological dynamics and anticipating the impact of strategy in controlling and preventing dengue.

Figure 7

Time series plots of the number of infected humans and mosquitoes predicted by the multi-serotype dengue model (2.4)–(2.5) with different initial sizes of the sub-populations. The values of model parameters used are given in Table 2 and $\lambda = 0.0111$ which corresponds to cross immunity lasting as long as three months after primary experience infection

Figure 8

Time series plots of the number of infected humans and mosquitoes predicted by the multi-serotype dengue model (2.4)–(2.5) with different initial sizes of the sub-populations. The values of model parameters used are given in Table 2 and $\lambda = 1.3699\times 10^{-3}$ which corresponds to two years of cross immunity after primary experience infection

Table 3 Effect of the period of cross immunity ($1/\lambda$) on the number of infected populations: $I_{{H1}}$, $I_{{H2}}$, $I_{{H12}}$, $I_{{H21}}$, $I_{{V1}}$, and $I_{{V2}}$ at steady state using the multi-serotype dengue model (2.4)–(2.5)

4.6 Effect of nonlinear incidence rates

It is evident from Table 3 that only the number of infected humans in secondary infection die out as the period of cross immunity is increased, while the infected humans in primary infection and the infected mosquitoes still be in the community. This may cause the dengue outbreak again. For these reasons, in the case when the period of cross immunity is three months, the effect of nonlinear incidence rates for human and mosquito is investigated with various the inhibitory effect of human $\alpha _{{H}}$ and the inhibitory effect of mosquito $\alpha _{{V}}$.

When the inhibitory effect of human $\alpha _{{H}}$ varies as $\alpha _{{H}}= 8.6480\times 10^{-5}$, 0.05, 0.125, respectively. With these values and the other parameter values given in Table 2, the area of secondary infection is shown in Fig. 9(a). It is found that the area of secondary infection decreases as $\alpha _{{H}}$ increases. A similar result is found in the case when the inhibitory effect of mosquito $\alpha _{{V}}$ is chosen to be $\alpha {_{V}} =1.3405\times 10^{-6}$, 0.125, 0.85, respectively, see Fig. 9(b). By comparing the region of secondary infection, it is observed that slight increase in $\alpha _{{H}}$ would reduce the region of secondary infection. It is also found that in order to reduce the area of secondary infection while increasing $\alpha _{{H}}$, the value of $\alpha _{{V}}$ must be close to unity. These results indicate that the inhibitory effect of human is sensitive to decreasing the region of secondary infection more than the inhibitory effect of mosquito.

Figure 9

Showing the region of secondary infection covered by the functions $f_{{1}}$ and $f_{{2}}$ given in (4.4) with various $\alpha _{{H}}$ and $\alpha _{{V}}$. (a) $\alpha _{{V}} =1.3405\times 10^{-6}$ and $\alpha _{{H}}$ varies. (b) $\alpha _{{H}} = 8.6480\times 10^{-5}$ and $\alpha _{{V}}$ varies

Further simulations are carried out using various values of $\alpha _{{H}}$ and $\alpha _{{V}}$ and are tabulated in Tables 4–6. It is evident from Table 4 that, as the effectiveness of $\alpha _{{H}}$ increases, the number of infected humans in primary infection decreases, which causes the numbers of infected humans in secondary infection and infected mosquitoes to greatly reduce and finally die out. Meanwhile, increasing the efficiency level of $\alpha _{{V}}$, the number of infected humans in primary infection decreases, which causes the number of infected humans in secondary infection to greatly reduce and finally die out. The results also show that although the effectiveness of $\alpha _{{V}}$ would decrease the number of infected mosquitoes, it would not eliminate the number of infected mosquitoes. This study suggests that the strategy of vector control only may induce the risk of reinfection dengue in these populations.

Table 4 Effect of $\alpha _{{H}}$ on the number of infected populations: $I_{{H1}}$, $I_{{H2}}$, $I_{{H12}}$, $I_{{H21}}$, $I_{{V1}}$, and $I_{{V2}}$ at steady state using the multi-serotype dengue model (2.4)–(2.5)

The study results also show that at 5% effectiveness of inhibitory effects (that is, $\alpha _{{H}}=\alpha _{{V}}=0.05$), the steady states $P_{{H}}^{{*}}$ and $P_{{n}}^{{**}}$ exchange their stability, see Tables 4–5. This new result verifies that even nonlinear incidence rates $\alpha _{{H}},\alpha _{{V}}$ are not parameters in the associated reproductive number with serotype-i, but they are important factors that influence exchange of the stability region of co-existence primary infection equilibrium $P_{{H}}^{{*}}$ and the stability region of co-existence secondary infection equilibrium $P_{{n}}^{{**}}$, which leads to eliminating the infected populations and decreasing the risk of reinfection dengue. Thus, the combination of inhibitory effect of human and mosquito is explored by varying $\alpha _{{H}}$ and $\alpha _{{V}}$. The obtained result is shown in Table 6. The results in Table 6 verify that if the efficiency level of $\alpha _{{H}}$ is 85%, the efficiency level of $\alpha _{{V}}$ must be 95% in order to eliminate the infected human in secondary infection and the infected mosquito. It is also shown that when $\alpha _{{H}}$ is close to unity and $\alpha _{{V}}$ slightly increases, that is, $\alpha _{{H}}=0.9$ and $\alpha _{{V}}=0.3$, these shall reduce the number of infected humans in primary infection, which leads to eliminating the infected humans in secondary infection and the infected mosquitoes. Therefore, these study results suggest that the best way of dengue prevention and control should be the effective vector control measures and personal protection.

Table 5 Effect of $\alpha _{{V}}$ on the number of infected populations: $I_{{H1}}$, $I_{{H2}}$, $I_{{H12}}$, $I_{{H21}}$, $I_{{V1}}$, and $I_{{V2}}$ at steady state using the multi-serotype dengue model (2.4)–(2.5)

Table 6 Effect of $\alpha _{{H}}$ and $\alpha _{{V}}$ on the number of infected populations: $I_{{H1}}$, $I_{{H2}}$, $I_{{H12}}$, $I_{{H21}}$, $I_{{V1}}$, and $I_{{V2}}$ at steady state using the multi-serotype dengue model (2.4)–(2.5)

5 Discussion and conclusion

In this paper, we formulate the multi-serotype dengue model with nonlinear incidence rate (2.4)–(2.5) for studying the influence of cross immunity and nonlinear incidence rate on the secondary infection. The main results from this study are as follows. The reproduction number of the multi-serotype dengue model derived by using the next generation method is given by $R_{{0}}^{{n}} = \operatorname{max}\{R_{{n1}}, R_{{n2}}\}$. This number depends on the reproduction number associated with infection serotype-i ($R_{{ni}}$, $i=1,2$) in the sense that the transmission of dengue serotype-1 or dengue serotype-2 in primary infection has strong influence on the possibility of secondary infection in some period when the cross immunity is weak. Mathematical analysis has shown that the multi-serotype dengue model has a globally stable disease-free equilibrium if $R_{{0}}^{{n}}\leq 1$, as guaranteed by Theorem 3.2. This result indicates that if public health measures that make (and keep) the threshold to a value less than unity are carried out, the number of infected human and mosquito populations in the community will be brought to zero irrespective of the initial sizes of sub-populations, see Fig. 2. The local stability of co-existence infection equilibria of the multi-serotype dengue model is analyzed by applying the center manifold theory. When $R_{{0}}^{{n}} > 1$, the endemic equilibria of the multi-serotype dengue model called the coexistence primary and secondary equilibria exist and are locally asymptotically stable as guaranteed by Theorems 3.2–3.3 and Lemmas 3.2–3.3, respectively.

The multi-serotype dengue model is applied to predict the number of infected humans and compare it with the real data reported by the Ministry of Public Health Thailand [30]. This study gives the appropriated values of model parameters of the formulated model. Sensitivity analysis shows that the model parameters $\beta _{{V2}}$, b, $\beta _{{2}}$, and $\eta _{{H2}}$ influence the transmission of dengue. On the other hand, the model parameters $\mu _{{V}}$ and $\gamma _{{1}}$ impact the prevention. In biological meaning, the sensitivity of parameter $\mu _{V}$ indicates that dengue prevention and control depends on effective vector control measures. Meanwhile, the sensitivities of parameters $\beta _{{V2}}$, b, $\beta _{{2}}$, and $\eta _{{H2}}$ indicate that reducing the number of contacts between humans and mosquitoes, through a reduction in either or both, the frequency of mosquito blood meals, and the number of bites that a human will tolerate would have the largest effect on disease transmission. The sensitivity of parameter $\eta _{{H}}$ also shows that the infected asymptomatic humans are the main carriers and multipliers of the virus, serving as a source of the virus for uninfected mosquitoes. In addition, the sensitivity of recovery rate in both primary and secondary infection humans indicates that early detection and access to proper medical care lowers fatality rates due to no specific treatment for dengue/severe dengue [55].

The appropriate values of model parameters relevant to dengue transmission dynamics in Thailand are used to investigate the effect of cross immunity and nonlinear incidence rate on the transmission of dengue and reinfection. The study results reveal that the cross immunity affects reinfection of dengue transmission. These study results suggest that increasing immunity to people who have previously had a dengue infection is the best way to prevent dengue, which is in line with the developing dengue vaccine [56], including the manufacturer recommendation that the vaccine only be used in people who have previously had a dengue infection, as outcomes may be worsened in those who have not been previously infected [57].

Further, the results are shown that the nonlinear incidence rates $\alpha _{{H}}$ and $\alpha _{{V}}$ affect exchange of the stabilities of the coexistence primary and secondary equilibria. The study results verify that increasing the effective level of $\alpha _{{H}}$ would decrease the numbers of the infected humans in secondary infection and infected mosquitoes and they would finally die out. On the other hand, increasing the effective level of $\alpha _{{V}}$ leads to elimination of the infected humans in secondary infection but it does not eliminate the number of infected mosquitoes, which may induce the risk of reinfection dengue in these populations. In addition, increasing both the effective levels of $\alpha _{{H}}$ and $\alpha _{{V}}$ shall eliminate both the infected humans in secondary infection and infected mosquitoes. Therefore, these study results suggest that the best way of dengue prevention and control should be the effective vector control measures and personal protection, especially boosting immunity should be provided to people who have previously had a dengue infection.

In summary, the multi-serotype dengue model considered here is refined to incorporate cross immunity between serotypes and nonlinear incidence rate. These findings will impact strategies for designing dengue vaccine studies and future multi-strain modeling efforts in order to understand the evolutionary pressures in multi-strain disease systems which will be the future work.

Appendix: Proof of Theorem 3.1

Proof

Consider the Lyapunov function for model (2.4)–(2.5)

$$\begin{aligned} \mathcal{F} ={}& h_{{1}}(E_{{H1}}+E_{{H21}})+h_{{2}}(I_{{H1}}+I_{{H21}})+h_{{3}}(E_{{H2}}+E_{{H12}}) +h_{{4}}(I_{{H2}}+I_{{H12}}) \\ &{}+h_{{5}}E_{{V1}}+h_{{6}}I_{{V1}}+h_{{7}}E_{{V2}}+h_{{8}}I_{{V2}}, \end{aligned}$$

where

$$\begin{aligned} &h_{{1}}=\varPi _{{V}}b\beta _{{V1}}(\eta _{{H1}}K_{{3}}+\sigma _{{1}}) \sigma _{{V1}} \gamma _{{1}},\qquad {} h_{{2}}=\varPi _{{V}}b\beta _{{V1}} \sigma _{{V1}}K_{{1}}\gamma _{{1}},\\ & h_{{3}}=\varPi _{{V}}b\beta _{{V2}}( \eta _{{H2}}K_{{4}}+\sigma _{{2}})\sigma _{{V2}} \gamma _{{2}}, \\ &h_{{4}}=\varPi _{{V}}b\beta _{{V2}}\sigma _{{V2}}K_{{2}}\gamma _{{2}}, \qquad{} h_{{5}}=\mu _{{V}}K_{{1}}K_{{3}}\sigma _{{V1}}\gamma _{1}R_{{n1}}^{{2}},\qquad {} h_{{6}}=\mu _{{V}}K_{{1}}K_{{3}}K_{{5}}\gamma _{1}R_{{n1}}^{{2}}, \\ &h_{{7}}=\mu _{{V}}K_{{2}}K_{{4}}\sigma _{{V2}}\gamma _{2}R_{{n2}}^{{2}},\qquad {} h_{{8}}=\mu _{{V}}K_{{2}}K_{{4}}K_{{6}} \gamma _{2}R_{{n2}}^{{2}}. \end{aligned}$$

The derivative of $\mathcal{F}$ along the solution of system (2.4)–(2.5) is given by

$$\begin{aligned} \frac{d \mathcal{F}(t)}{dt} ={}& \varPi _{{V}}b\beta _{{V1}}( \eta _{{H}1}K_{{3}}+ \sigma _{{1}})\sigma _{{V1}}\gamma _{{1}} \biggl[ \frac{b\beta _{{1}}I_{{V1}}}{1+\alpha _{{V1}}I_{{V1}}}S_{{H}}-K_{{1}}E_{{H1}} \biggr] \\ &{}+ \varPi _{{V}}b\beta _{{V1}}\sigma _{{V1}}K_{{1}} \gamma _{{1}} [\sigma _{{1}}E_{{H1}}-K_{{3}}I_{{H1}} ] \\ &{}+\varPi _{{V}}b\beta _{{V2}}(\eta _{{H}2}K_{{4}}+ \sigma _{{2}}) \sigma _{{V2}}\gamma _{{2}} \biggl[ \frac{\lambda _{{2}}b\beta _{{2}}I_{{V2}}}{1+\alpha _{{V2}}I_{{V2}}}R_{{H1}}-K_{{2}}E_{{H12}} \biggr] \\ &{}+ \varPi _{{V}}b\beta _{{V2}}\sigma _{{V2}}K_{{2}} \gamma _{{2}} [\sigma _{{2}}E_{{H12}}-K_{{4}}I_{{H12}} ] \\ &{}+ \varPi _{{V}}b\beta _{{V2}}(\eta _{{H}2}K_{{4}}+ \sigma _{{2}}) \sigma _{{V2}}\gamma _{{2}} \biggl[ \frac{b\beta _{{2}}I_{{V2}}}{1+\alpha _{{V2}}I_{{V2}}}S_{{H}}-K_{{2}}E_{{H2}} \biggr] \\ &{}+\varPi _{{V}}b\beta _{{V2}}\sigma _{{V2}}K_{{2}} \gamma _{{2}} [\sigma _{{2}}E_{{H2}}-K_{{4}}I_{{H2}} ] \\ &{}+ \varPi _{{V}}b\beta _{{V1}}(\eta _{{H}1}K_{{3}}+ \sigma _{{1}}) \sigma _{{V1}}\gamma _{{1}} \biggl[ \frac{\lambda _{{1}}b\beta _{{1}}I_{{V1}}}{1+\alpha _{{V1}}I_{{V1}}}R_{{H2}}-K_{{1}}E_{{H21}} \biggr] \\ &{}+\varPi _{{V}}b\beta _{{V1}}\sigma _{{V1}}K_{{1}} \gamma _{{1}} [\sigma _{{1}}E_{{H21}}-K_{{3}}I_{{H21}} ] \\ &{}+\mu _{{V}}K_{{1}}K_{{3}}\sigma _{{V1}} \gamma _{1}R_{{1}}^{{n}} \biggl[ \frac{b\beta _{{V1}} (\eta _{{H1}} (E_{{H1}}+E_{{H21}} )+I_{{H1}}+I_{{H21}} )}{1+\alpha _{{H1}}(I_{{H1}}+I_{{H21}})}S_{{V}}-K_{{5}}E_{{V1}} \biggr] \\ &{}+\mu _{{V}}K_{{1}}K_{{3}}K_{{5}}\gamma _{1}R_{{1}}^{{n}} [\sigma _{{V1}}E_{{V1}}- \mu _{{V}}I_{{V1}} ] \\ &{}+\mu _{{V}}K_{{2}}K_{{4}}\sigma _{{V2}} \gamma _{2},R_{{2}}^{{n}} \biggl[ \frac{b\beta _{{V2}} (\eta _{{H2}} (E_{{H2}}+E_{{H12}} )+I_{{H2}}+I_{{H12}} )}{1+\alpha _{{H2}}(I_{{H2}}+I_{{H12}})}S_{{V}}-K_{{6}}E_{{V2}} \biggr] \\ &{}+\mu _{{V}}K_{{2}}K_{{4}}K_{{6}}\gamma _{2}R_{{2}}^{{n}} [\sigma _{{V2}}E_{{V2}}- \mu _{{V}}I_{{V2}} ] \\ \leq{} & \varPi _{{V}}b\beta _{{V1}}\sigma _{{V1}}\eta _{{H1}}K_{{1}}K_{{3}} \gamma _{{1}} \bigl[R_{{n1}}^{{2}}-1 \bigr]E_{{H1}} + \varPi _{{V}}b \beta _{{V1}}\sigma _{{V1}}K_{{1}}K_{{3}} \gamma _{{1}} \bigl[R_{{n1}}^{{2}}-1 \bigr]I_{{H1}} \\ &{}+ \varPi _{{V}}b\beta _{{V2}}\sigma _{{V2}}\eta _{{H2}}K_{{2}}K_{{4}} \gamma _{{2}} \bigl[R_{{n2}}^{{2}}-1 \bigr]E_{{H12}} +\varPi _{{V}}b \beta _{{V2}}\sigma _{{V2}}K_{{2}}K_{{4}} \gamma _{{2}} \bigl[R_{{n2}}^{{2}}-1 \bigr]I_{{H12}} \\ &{}+ \varPi _{{V}}b\beta _{{V2}}\sigma _{{V2}}\eta _{{H2}}K_{{2}}K_{{4}} \gamma _{{2}} \bigl[R_{{n2}}^{{2}}-1 \bigr]E_{{H2}} +\varPi _{{V}}b \beta _{{V2}}\sigma _{{V2}}K_{{2}}K_{{4}} \gamma _{{2}} \bigl[R_{{n2}}^{{2}}-1 \bigr]I_{{H2}} \\ &{}+ \varPi _{{V}}b\beta _{{V1}}\sigma _{{V1}}\eta _{{H1}}K_{{1}}K_{{3}} \gamma _{{1}} \bigl[R_{{n1}}^{{2}}-1 \bigr]E_{{H21}} +\varPi _{{V}}b \beta _{{V1}}\sigma _{{V1}}K_{{1}}K_{{3}} \gamma _{{1}} \bigl[R_{{n1}}^{{2}}-1 \bigr]I_{{H21}} \\ &{}+ \mu _{{V}}^{{2}}K_{{1}}K_{{3}}K_{{5}} \gamma _{{1}}R_{{n1}}^{{2}} \bigl[R_{{n1}}^{{2}}-1 \bigr]I_{{V1}} + \mu _{{V}}^{{2}}K_{{2}}K_{{4}}K_{{6}} \gamma _{{2}}R_{{n2}}^{{2}} \bigl[R_{{n2}}^{{2}}-1 \bigr]I_{{V2}} \\ = {}& \varPi _{{V}}b\beta _{{V1}}\sigma _{{V1}}K_{{1}}K_{{3}} \gamma _{{1}} \bigl[\eta _{{H1}}(E_{{H1}}+E_{{H21}})+I_{{H1}}+I_{{H21}} \bigr] \bigl[R_{{n1}}^{{2}}-1 \bigr] \\ &{}+\varPi _{{V}}b\beta _{{V2}}\sigma _{{V2}}K_{{2}}K_{{4}} \gamma _{{2}} \bigl[\eta _{{H2}}(E_{{H2}}+E_{{H12}})+I_{{H2}}+I_{{H12}} \bigr] \bigl[R_{{n2}}^{{2}}-1 \bigr] \\ &{}+ \mu _{{V}}^{{2}}K_{{1}}K_{{3}}K_{{5}} \gamma _{{1}}R_{{n1}}^{{2}} \bigl[R_{{n1}}^{{2}}-1 \bigr]I_{{V1}} + \mu _{{V}}^{{2}}K_{{2}}K_{{4}}K_{{6}} \gamma _{{2}}R_{{n2}}^{{2}} \bigl[R_{{n2}}^{{2}}-1 \bigr]I_{{V2}}. \end{aligned}$$

It is obvious that ${\frac{d \mathcal{F}(t)}{dt} < 0}$ whenever $R_{{0}}^{{n}}= \operatorname{max} \{ R_{{n1}}, R_{{n2}} \} <1$ and ${\frac{d \mathcal{F}(t)}{dt}=0}$ if and only if $E_{{H1}}=I_{{H1}}=E_{{H12}}=I_{{H12}}=E_{{H2}}=I_{{H2}}=E_{{H21}}=I_{{H21}}=I_{{V1}}=I_{{V2}}=0$. This indicates that the maximum invariant set in $\{ (S_{{H}},E_{{H1}},I_{{H1}},R_{{H1}},E_{{H12}},I_{{H12}},E_{{H2}},I_{{H2}}, R_{{H2}},E_{{H21}},I_{{H21}}, R_{{H22}},S_{{V}}, E_{{V1}},I_{{V1}},E_{{V2}},I_{{V2}} ) \in \varOmega: {\frac{d \mathcal{F}(t)}{dt}=0} \}$ is the singleton ${P_{{n}}^{{o}}}$. Hence, $\mathcal{F}$ is a Lyapunov function in Ω. Thus, it follows from LaSalle’s invariance principle [58] that

$$\begin{aligned} \bigl(E_{{H1}}(t),I_{{H1}}(t),E_{{H12}}(t),I_{{H12}}(t), E_{{H2}}(t),I_{{H2}}(t),E_{{H21}}(t),I_{{H21}}(t) \bigr) \rightarrow (0,0,0,0,0,0,0,0) \end{aligned}$$

as $t\rightarrow \infty$.

Substituting $E_{{H1}}=I_{{H1}}=E_{{H12}}=I_{{H12}}=E_{{H2}}=I_{{H2}}=E_{{H21}}=I_{{H21}}=I_{{V1}}=I_{{V2}}=0$ into the first and the thirteenth of model (2.4)–(2.5) gives $S_{{H}}(t) \rightarrow S_{{H}}^{{*}}$ and $S_{{V}}(t) \rightarrow S_{{V}}^{{*}}$ as $t\rightarrow \infty$. Therefore, every solution of multi-serotype dengue model (2.4)–(2.5) approaches the DFE $(P_{{n}}^{{o}})$, as $t\rightarrow \infty$ for $R_{{0}}^{{n}}= \operatorname{max} \{ R_{{n1}}, R_{{n2}} \} <1$, so that $P_{{n}}^{{o}}$ is GAS in Ω if $R_{{0}}^{{n}}= \operatorname{max} \{ R_{{n1}}, R_{{n2}} \} <1$. □