Analysis of a stage-structured dengue model

In order to study the impact of control measures and limited resource on dengue transmission dynamics, we formulate a stage-structured dengue model. The basic investigation of the model, such as the existence of equilibria and their stability, have been proved. It is also shown that this model may undergo backward bifurcation, where the stable disease-free equilibrium co-exists with an endemic equilibrium. The backward bifurcation property can be removed by ignoring the disease-induced death in human population and the global stability of the unique endemic equilibrium has been proved. Sensitivity analysis with respect to \begin{document}$R_0$\end{document} has been carried out to explore the impact of model parameters. In addition, numerical analysis manifests that the more intensive control measures in targeting immature and adult mosquitoes are both effective in preventing dengue outbreaks. It is also shown that the earlier the control intervention begins, the less people would be infected and the earlier dengue would be eradicated. Even later epidemic prevention and control can also effectively reduce the severity of pandemic. Moreover, comprehensive control measures are more effective than a single measure.


1.
Introduction. Dengue, a mosquito-borne disease caused by any of four closelyrelated virus serotypes of the genus Flavivirus, is endemic in at least 100 countries in Africa, the Americas, the Eastern Mediterranean and subtropical regions of the world, inhibited by over 2.5 billion people( [6,29]). Dengue ranks second to malaria amongst deadly mosquito-borne diseases, each year claiming about 100 million infections and 20,000 deaths globally [29].
Dengue exhibits as many as four coexisting types of serotype (strains) in a region. Infection by any dengue virus strain produces long-lasting immunity but only temporary cross-immunity to other serotypes. A person infected by one of the four serotypes will lose immunity to the three other serotypes (heterologus immunity) in about 12 weeks and then becomes more susceptible to developing dengue haemorrhagic fever. Three of the vectors are Aedes aegypti Linnaeus, Aedes albopictus Skuse, and Aedes scutellaris Walk. Dengue infection follows the bite of a competent mosquito vector, principally Aedes aegypti (L.). Aedes aegypti mosquitoes acquire infection from infected individuals 6 to 18 h before onset of fever and then for the duration of the fever. A minimum extrinsic incubation period of 8 to 14 days is required after an infective blood meal before the mosquito becomes infectious [18].
The early models of dengue can be back to the models in [25]. Since then, more and more complicated models have been developed ( [3], [4], [7], [11], [12], [13], [14], [18], [19], [20], [32] etc). Feng et al.( [18]) studied a system that models the population dynamics of a SIR vector transmitted disease with two pathogen strains and argued that the existence of competitive exclusion in their system is a product of the interplay between the host superinfection process and frequency-dependent (vector to host) contact rates. Esteva et al.( [11]) proposed a one-serotype dengue model with a constant human population and variable vector population. With the assumption that the human population was supposed to grow exponentially, the same authors proposed another model considering only one serotype [12]. In [13], the impact of vertical transmission and interrupted feeding on the dynamics of the disease has been studied. Besides, a multi-serotype of dengue was also studied by them in [14]. They discussed conditions for the asymptotic stability of equilibria, supported by analytical and numerical methods and found that coexistence of both types of serotype was possible for a large range of parameters. Billings et al.([4]) formulated a multi-serotype disease model which did not include vector variables investigated the complex dynamics induced by Antibody-dependent Enhancement (ADE). Considering the effect of vaccination, a new multi-serotype disease model was studied in [3]. In addition, Shaw et al.( [28]) formulated and analyzed a compartmental model for multiple types of serotype exhibiting ADE. Using center manifold techniques, they showed how the dynamics rapidly collapsed to a lower dimensional system. In [10], authors studied a SEIR dengue by examining the role of temperature in driving vector dynamics.
Several experiments suggest that rate of mosquito development and the mosquito lifespan fluctuate with changes in temperature. Temperature also impacts the extrinsic incubation period of the virus within a mosquito host [5,34]. Moreover, the population dynamics of mosquitoes is closely related to the environment and resource in a region [2,17,22,33]. One can see that few models mentioned above include the immature mosquito stage which is more sensitive to climate change and resource. Recently, several mosquito-borne disease models incorporating stagestructured mosquitoes have been studied ( [1,22], etc.) Li ([22]) studied a malaria model with a aquatic-stage class of mosquitoes. A nonlinear maturation rate was assumed to model the intraspecific competition between larvae. In [1], Ai et al. formulated a mosquito-stage-structured malaria model which include four distinct metamorphic stages of mosquitoes. Nevertheless, the competition for reproduction resource (like blood meals and water reservoir) between adult female mosquitoes has not been considered.
Some control measures, like spraying insecticide and clearing standing water, would be implemented to control dengue [5]. While, mathematical models can provide useful strategic insights into control measures for dengue [23].
It is our aim to formulate a new dengue transmission model incorporating the immature mosquito stage explicitly and study the impact of the limited resource and some control measures on mosquito population and dengue transmission, as well as assess the effectiveness of different intervention strategies.
The rest of the paper is organized as follows. A dengue model incorporating some control measures is derived in Section 2. In Section 3, the basic investigation of the model, such as the existence and stability of equilibria, is performed in this part. Some numerical analysis and discussion are given in Section 4.
2. Model formulation. Mosquito life includes four distinct stages of development during a lifetime: egg, larva, pupa, and adult [2,22]. We group the three aquatic stages of mosquitoes into one class and divide the mosquito population into only two class, one of which consists of the first three stages, denoted by A(t), and the other one of which is the adult stage, denoted by M (t).
The total population of adult mosquitoes M (t) is subdivided into three classes: susceptible M S (t), exposed M E (t), infective M I (t). We ignore recovered mosquitoes in the model because of the short lifespan of mosquitoes. The total population of humans H(t) is subdivided into four classes: susceptible H S (t), exposed H E (t), infective H I (t), and recovered H R (t). All the state variables are listed in Table 1 Since the resource related to mosquito reproduction is limited, we assume that there is competition between female adult individuals for finding an appropriate water reservoir and blood meals to lay eggs. It's natural to adopt a saturated birth rate functionp A M/(1 + M ) to reflect the restriction of limited resource for reproduction of eggs, wherep A is the maximum value of the recruitment rate of viable mosquito eggs. The idea is inspired by [17,33]. γ a is the maturation rate of larvae.d a and κ are the density independent and dependent death rate, respectively, and the death rate of adults isd m .
It is worth emphasizing that it is different from the idea using in [22] in which the birth rate of mosquito eggs is not affected by limited resource and there is intraspecific competition between immature mosquitoes during the maturation process.
Let b be the number of bites on a human per mosquito per unit of time, β h is the probability of infection to a human host per bite and β m is the probability of infection per bite to a susceptible mosquito from an infective human host. Then, the infection rate of a human The major dengue control effort includes clearing standing water and spraying insecticide in targeting immature and adult mosquitoes, respectively. Clearing standing water decreases the probability that a female mosquito finds an appropriate water reservoir to lay eggs and kills some existing immature mosquitoes. While spraying insecticide can kill adult mosquitoes. We assume that the maximum of viable egg recruitment rate be (p A −α b ), the death rate of immature mosquitoes be (d a + α J ) and the death rate of adult mosquitoes be (d m + α v ), where α b , α J and α v are control intervention parameters.
We let Ψ h be the input flows of the susceptible humans including births, d h is the natural death rates of humans. δ h is disease-induced death rates for humans. 1/γ h and 1/γ m are the incubation period of humans and mosquitoes, respectively. η h is the recovery rate of humans. ψ h is the rate of loss of immunity for recovered humans. It is worth emphasizing that, unlike many of the published modeling studies on dengue transmission dynamics, the current study assumes that an infected human host can obtain temporary immunity, but he loses immunity some times later and then becomes susceptible again. (see, for instance, [8,9,32]). Fig. 1 shows a flow diagram to describe the transmission dynamics of dengue virus between humans and mosquitoes.
Hence, we can obtain the following system to describe the dynamic of the dengue transmission: Description of parameters used in the model is given in Table 2.
If we assume that (2) becomes the model studied in [16]. Let R m 0 = p A γa dm(γa+da) . According to the result of [16], we have Lemma 2.1. (Theorem in [16]) If R m 0 < 1, the trivial equilibrium (0,0) of the system (2) is a locally asymptotically stable, and there exists no positive equilibrium. If R m 0 > 1, the trivial equilibrium (0,0) of system (2) is unstable, and there exists a unique positive equilibrium (A 0 , M 0 ), which is globally asymptotically stable, where In the following, we assume that R m 0 > 1 and the positive equilibrium (A 0 , M 0 ) is globally asymptotically stable.
The continuity of the right side of the previous system and its derivatives implies that unique solutions exist in the maximal interval. Since solutions approach, enter or stay in Ω they are eventually bounded and hence exist for t ≥ 0. Therefore the model is mathematically and epidemiologically well posed.
3. Dynamical analysis. [31], the linear stability of E 0 can be established using the next generation operator method. The matrix F (the new infection terms) and V (the It follows then that the basic reproduction number, denoted by R 0 , is given by Therefore, using Theorem 2 of [31], we have established the following results:

ANALYSIS OF A STAGE-STRUCTURED DENGUE MODEL 4051
Theorem 3.1. For the model (1), the disease-free equilibrium E 0 is locally asymptotically stable if R 0 < 1, and unstable if R 0 > 1.
The epidemiological implication of Theorem 3.1 is that, generally speaking, when R 0 is less than unity, a small influx of infected mosquitoes into a completely susceptible community would not generate large outbreaks, and the disease dies out in time. However, we show in the next subsection that the disease may still persist even if R 0 < 1 owing to the backward bifurcation.
3.2. Existence of endemic equilibria and backward bifurcation. Using the method used in [22], we next explore the existence of endemic equilibria. Let the right-hand side of system (1) equal to zero and solve H S , H E , H I , H R in terms of λ h , then Therefore, where .It is easy to see that, at the equilibrium, Then, substituting (6) into the right-hand side of (1) yields Obviously, (7) has a unique positive root A 0 , if R m 0 > 1. Solving M S , M E , M I in terms of λ v , we have Substituting (10) into (9) yields where The solution λ h = 0 corresponds to the disease-free equilibrium E 0 . We are now trying to find a positive solution of (11). Let Then, (11) is equivalent to F (λ h ) = 0 when λ h > 0. Moreover,
Based on the above analysis, we have the following results. If ξ 2 +ξ 4 < ξ 3 , there exist a threshold value 0 < R c < 1, such that system (1) has no endemic equilibrium when R 0 < R c , system (1) has a unique endemic equilibrium when R 0 = R c and system (1) has two endemic equilibrium if R c < R 0 < 1, which implies that a backward bifurcation appears. When ξ 2 + ξ 4 ≥ ξ 3 , system (1) has a unique endemic equilibrium if and only if R 0 > 1 and backward bifurcation can not occur.
Epidemiologically, our results suggest that, if ξ 2 + ξ 4 < ξ 3 , even though the basic reproduction number is smaller than unity, there may be a stable endemic equilibrium due to backward bifurcation of the model and the basic reproductive number itself is not enough to describe whether dengue will prevail or not and a new threshold value R c is needed. If R c < R 0 < 1, we should pay more attention to the initial sizes of the involved populations (See Fig. 2). The associated bifurcation diagram is depicted in Fig. 3. 3.3. Stability of endemic equilibrium. If δ h ≤ d h , one can easily prove that Then, system (1) has no endemic equilibrium if R 0 < 1 and backward bifurcation can not occur. Therefore, disease-induced death rate is big enough is a necessary condition for backward bifurcation. The backward bifurcation property can be removed by ignoring the disease-induced death in the human population. It is difficult to determine the stability of the unique endemic equilibrium when R 0 > 1. In the following, we are going to prove that the unique endemic equilibrium where A * = A 0 , is globally asymptotically stable by ignoring the disease-induced death and temporal immunity in the human population. The method used here is inspired by [1]. The main results are given as follows.
. It follows from Lemma 2.1 that A(t) → A * as t → ∞. Then, the limited system of (1) is the same as the limited system of Model (2.1) in [1] when δ h = ψ h = 0. By the proof of Theorem 6.1 in [1], one can easily prove that E * is locally stable and φ 1 (t) → E 1 * as t → ∞, if φ(0) ∈ D\D 0 . Therefore, the unique endemic equilibrium E * is globally asymptotically stable. 4. Numerical analysis and discussion. In this paper, focusing on the impact of the limited resource and some control measures on the transmission dynamics of dengue, we formulated a dengue model with stage-structured mosquito population. The basic investigation of the model, such as the existence of equilibria and their stability, have been finished. It was also shown that the stage-structured dengue model may undergo backward bifurcation, where the stable disease-free equilibrium co-exists with an endemic equilibrium, which implies that even though the basic reproduction number is less than unity, there may be an endemic equilibrium due to the backward bifurcation of the model. The basic reproductive number itself is not enough to describe whether dengue will prevail or not and a new threshold value R c is needed. If R c < R 0 < 1, we should pay more attention to the initial sizes of the involved populations. Our results also suggested that the backward bifurcation property can be removed by ignoring the disease-induced death and the global stability of the unique endemic equilibrium has been proved.
In order to model the impact of limited resource, we assume the growth rate function of adult mosquitoes be a saturated function p A M 1+M . Fig. 4 (a) shows that, unlike the linear growth rate function or constant growth rate function, the growth rate increases and tends to a constant with the increasing of the number of adult mosquitoes M . Fig. 4 (b) depicts the curve of the growth rate function with respect to time t, which also tends to a constant. Therefore, the saturated growth rate function can depict the restriction of limited resource related to mosquito reproduction, like blood meal resource.
Temperature could cause the variation of the following six model parameters: the biting rate, b, the mosquito incubation time, 1/γ m , the duration of the whole cycle from egg laying to an adult mosquito eclosion, 1/γ a , the maximum value of the recruitment rate of viable mosquito eggs, p A , the death rate of immature mosquitoes, d a , and the death rate of adult mosquitoes, d m . Clear relationships between temperature and these parameters were noted with corresponding equations in [23] without control intervention. The temperature range is assumed to be [10. In the following, we are going to study the global and local sensitivity indexes of above six model parameters with respect to the reproductive number of infection, R 0 . Sensitivity analysis (SA) is defined as the study of how uncertainty in the quantity of interest(QoI) is attributed to different sources of uncertainties ( [26],  [27], [30], etc.). It helps to understand how the parameters of a model affect the quantity of interest(QoI), in our case the reproductive number of infection, R 0 . SA is often employed to rank the model parameters in order of their influence so that we can determine the most important model parameters. Fig. 5 demonstrates the global Sobol sensitivity analysis of these parameters. We observe that d m is the most global sensitive model parameter with respect to the basic reproductive number R 0 . b, γ a , p A and γ m are the second, third, forth and fifth most global sensitive model parameters, respectively. In addition, the death rate of immature mosquitoes, d a , is the least sensitive parameter. This global sensitivity study indicates that temperature variation brings the variation of R 0 and we should pay more attention to obtaining accurate data, especially those more sensitive parameters, like d m and b, if we want to have an accurate estimation of R 0 . In addition, local sensitivity analysis around the corresponding mean position is also performed to investigate the local sensitivity effect of above six model parameters on the basic reproductive number, R 0 . The local sensitivity analysis can provide the information on how a small change of a model parameter will cause the variation of R 0 at a given local setting of the six model parameters. Different from the global sensitivity analysis, the local sensitivity analysis does not depend on the range of the model parameters [27]. The magnitude of the local sensitivity indicates how much the variation of R 0 will be, given a small change of a model parameter. The sign of the local sensitivity indice indicates whether the variation increase of one model parameter will cause the variation increase or decrease of R 0 . As shown in Fig. 6, d m is most sensitive and the following are γ a , b and γ m , which is different from the global sensitivity study. Moreover, p A and d a have very small local sensitivity. Since the local sensitivity indexes of γ a , b, γ m , and p A are positive, the variation of the reproductive number of infection, R 0 , will increase as long as the variation of one of the four model parameters is increasing. R 0 will decrease as long as the variation of one of the two model parameters d a and d m are increasing according to the negative local sensitivity indexes. Therefore, temperature increase exacerbates dengue transmission because it brings the increase of the values of γ a , r, γ m , p A , as well as the decrease of the value of d a and d m .
The major dengue control effort includes clearing standing water and spraying insecticide in targeting immature and adult mosquitoes respectively. Clearing standing water decreases the probability that a female mosquito finds an appropriate water reservoir to lay eggs (decreases p A ) and kills some existing immature mosquitoes (increases d a ). While spraying insecticide can kill adult mosquitoes (increases d m ). We now focus on investigating the effectiveness of this two kinds of control measures and choose the intervention parameters α J , α b and α v , as critical parameters.  The contour of the basic reproduction number, R 0 as a function of α J , α b and α v are presented in Fig. 7 and Fig. 8. The parameter ranges have been extended for the sake of exploring the efficiency of control intervention and catching the trend of R 0 . The results manifest that more intensive control measures in targeting immature and adult mosquitoes are both effective in preventing dengue outbreaks.  Figure 11. The cumulative number of infected human hosts according to different control intervention time with single control intervention measure. Parameter values used are the same as that used in Fig. 9 except d m , p A and d a .
In order to study the impact of intervention time on dengue transmission, we carried out the numerical analysis on cumulative number of infected human hosts C(t) and dengue prevailing time according to three different control intervention strategies in Fig. 9 and Fig. 10. For the control intervention, we assume that α v = 0.04 and α J = 0.55 and α b = 4e + 10. Our results show that if the beginning time of intervention is at 30th day, the ultimate cumulative number of infected human hosts C(∞) will shrink from 407 to 206, which is a 49.4% decrease. And the prevailing time will be reduced from 460 days to 150 days. While, if the beginning time is at 10th day, C(∞) will shrink to 127, which is a 68.8% decrease. And the prevailing time will be further reduced to 136 days. Hence, the earlier the control intervention begins, the less people will be infected and the earlier dengue will be eradicated. Even later epidemic prevention and control can also effectively reduce the severity of pandemic. In addition, according to comparing Fig. 9 and Fig. 11, it can be asserted that comprehensive control measures are more effective than a single measure. Clearing standing water and spraying insecticide should be carried out simultaneously. The results could be helpful in guiding to use the most effective interventions and determine the intervention time in future dengue outbreaks with a limited budget.
Owing to the existence of backward bifurcation, even small fluctuation of parameters can bring a completely different result, especially as the basic reproduction number is near the subthreshold value R c . In order to control dengue transmission, people should concentrate on environmental and climate change, pay more attention to those more sensitive parameters, monitor the number of new cases, and take intervention measures as early as possible.
As an initial work, we are trying to model the impact of resource and some control measures on mosquito population and dengue transmission. Nevertheless, mosquito population may be affected by some other environmental and climate factors, like humidity and precipitation. Dengue transmission can also be affected by population mobility, vaccination, and other factors. We'll incorporate them into our model in the future work.