Modeling the effect of temperature on dengue virus transmission with periodic delay differential equations

1.   Introduction

Dengue fever is a rapidly spreading mosquito-borne disease all over the world, and dengue virus (DENV) is transmitted to human by Aedes (stegomyia) mosquito, a genus of mosquitoes. Aedes aegypti and Aedes albopictus are the vectors for transmitting dengue virus ^[1], and dengue fever can be triggered by dengue virus serotypes: DENV-1, DENV-2, DENV-3, DENV-4 ^[2,3]. When human infects the DENV, the symptoms of illness can range from mild forms such as a sudden onset of fever, headache, myalgia and a skin erythema ^[4], to severe forms including dengue hemorrhagic fever (DHF) and dengue shock syndrome (DSS) ^[5]. Since severe outbreaks of dengue fever in the 1950s in Philippines and Tahiland, it has spread throughout the Africa, Southeast Asia, Western Pacific Region, South America and Eastern Mediterranean Region, which leads to roughly 12,500 to 25,000 deaths each year ^[6,7].

According to the World Health Organization (WHO), about 50–100 million dengue cases are reported every year over 100 countries ^[7,8], and more than 2.5 billion people meet the risk of dengue virus transmission. There are very serious infection of dengue fever in Guangdong province, China ^[9]. It has been reported that there were 45230 dengue cases in Guangdong Province in 2014, which is a big outbreak of dengue fever in Guangdong. In recent decades, dengue fever has leaded a serious health and social problem in the world, and it has turned into a worldwide public health problem.

Mathematical models are used to study the transmission dynamics of dengue fever. It is the first time that Esteva et al. ^[10] proposed a dengue model and established the global stability of the endemic equilibrium. Later on, they analyzed the effect of mechanical and vertical transmission on the transmission dynamics of dengue virus in ^[11]. Chowell et al. ^[12] evaluated the basic reproduction number of dengue with data, and Khan et al. ^[13] studied a single-strain dengue model and the basic reproduction number was analyzed. Tewa et al. ^[14] developed the single-strain dengue model and proved the global stability of the equilibria. Yang and Ferreira ^[15] assessed the effect of vector control on dengue virus transmission, and obtained that the dengue fever can be eradicated if mosquito population dies out. Garba et al. ^[16] studied the backward bifurcation in dengue virus transmission model, and revealed that the implement of an imperfect vaccine could result in the effective control of dengue in a community provided that the vaccine efficacy and vaccine coverage are high enough. Abdelrazec et al. ^[17] investigated the effect of available resources of health system on the spread and control of dengue fever. Although a few studies have explored the dynamics of dengue virus transmission, there are few work addressing on the effect of climatic factors on dengue virus transmission.

Transmission of dengue virus is greatly influenced by temperature and precipitation ^[18,19]. Rueda et al. ^[20] showed that circadian rhythm of mosquito population is hugely affected by temperature and rainfall. Tun-Lin et al. ^[21] studied the impact of temperature on maturation period of mosquitoes, and they showed that the duration of development from egg to adult is closely related to temperature through the experiment methods, and found that higher temperature can shorten the duration of virus replication and increase mosquito reproduction and contacts with human ^[22]. Fan et al. ^[23] analyzed a delay differential equation model, and explored the effect of temperature on the mosquito abundance. Li et al. ^[24] explored the driving force of the dengue outbreak in Guangdong province in 2014 through a seasonally-driven model. Bai et al. ^[25] studied the effect of of climate change on the transmission of mosquito-borne diseases in China, and found that global warming extends the land area suitable for vector mosquitoes, and expands distribution of dengue. Therefore, it is necessary and important to study the effect of precipitation and temperature on the transmission of dengue virus.

In this paper, we propose a dengue model with maturation delay for mosquito production, where the birth rate and biting rate of mosquitoes depend on the average temperature. The basic reproduction number is computed by the next generation matrix. We prove the global stability of the disease-free steady state and uniform persistence of the system which are determined entirely by the basic reproduction number. Furthermore, we investigate the outbreak of dengue in Guangdong Province in 2014 through our model and data.

This paper is organized as follows. In section 2, we formulate a dengue virus transmission model, and compute the basic reproduction number $\mathbb{R}_0$. In section 3, the dynamics of the dengue model are completely determined by the basic reproduction number $\mathbb{R}_0$. If $\mathbb{R}_0 < 1$, the disease-free steady state is globally asymptotically stable and dengue fever dies out; if $\mathbb{R}_0 > 1$, there is at least one positive periodic solution and dengue fever persists. Section 4 provides a case study on dengue fever in Guangdong to illustrate our results which can be applied to other mosquito-borne diseases. Finally, a brief discussion completes the paper.

2.   Model formulation

The mosquito population $N_m(t)$ is divided into two compartments: The susceptible female adult mosquitoes $M_s(t)$ and the infectious female adult mosquitoes $M_i(t)$. Let $N_m(t) = M_s(t)+M_i(t)$ represent the total number of female adult mosquitoes. Here we ignore the recovered class of mosquitoes, because a mosquito never recovers from the infection in the short life span ^[26]. Let $N_h(t)$ denote the number of the total human population at time $t$, which is composed of the number of susceptible individuals $H_s(t)$, the number of infectious individuals $H_i(t)$, and the number of recovered individuals $H_r(t)$. Then $N_h(t) = H_s(t)+H_i(t)+H_r(t)$. The flow chart (Figure 1) depicts the transmission cycle of dengue virus between mosquitoes and humans.

Based on the above assumptions, and extending the ideas in ^{[23,27,28,29,30]}, our dengue virus transmission model can be written as

where all parameters are positive, and $r_m(t)$ and $k(t)$ are positive, continuous $\omega$-periodic functions with $\omega = 12$ months. Due to the density dependence of mosquito development from egg to adult, a Ricker type function is chosen as the birth rate of mosquitoes ^[29,31,32]. Let $\tau$ be the maturation period of mosquitoes, and $r_m$ is the intrinsic oviposition rate, and $d_j$ is the death rate of juvenile mosquitoes. We assume that $d_m$ is the natural death rate of adult female mosquitoes, and $d_h$ is the natural death rate of humans. Then $r_m(t)N_m(t-\tau)$ denotes the number of mosquito's eggs laid $\tau$ units ago, and $e^{-d_j\tau}$ denotes the survival probability during immature stages. Let $k(t)$ be the per capita biting rate of mosquitoes on humans, and $\beta_m(\text{resp. }\beta_h)$ represents the transmission probability from an infectious human to a susceptible mosquito (resp. from an infectious mosquito to a susceptible human). Since the birth rate $r_m(t)$ depending on the temperature affects the demographic of vector mosquitoes, and biting rate $k(t)$ depending on the temperature influences the virus transmission between humans and mosquitoes, then we assume that $r_m(t)$ and $k(t)$ are $\omega$-periodic functions. The parameter $r_h$ is the recruitment rate of human and $\mu_h$ denotes the recovery rate of the infectious human.

It follows from model (2.1) that the mosquito population $N_m(t) = M_s(t)+M_i(t)$ contents the following equation:

From ^[33,34,35], we know that there exists a globally asymptotically stable positive $\omega$-periodic solution $N^*_m(t)$ in model (2.1). In this paper, we consider the dengue virus transmission model with periodic environment. Then we assume

$($H$)$ A function $r_m(t)$ is periodic in $C^1(R_+)$ with period $\omega$, and then we can choose a $l_0 > 0$ which gives $r_m(t)e^{-d_j\tau}e^{-\alpha L}-d_mL < 0, \forall L > l_0$. Besides, there is a globally asymptotically stable positive $\omega$-periodic solution $N^*_m(t)$ in $C([-\tau, 0], R_{+})\backslash \{0\}$ for model (2.1).

Thus, the following results are established.

Theorem 2.1. For $\forall \phi \in C([-\tau, 0], R_{+}^{5})$, there is a unique nonnegative solution of model (2.1) through $\phi$, and all solutions of model (2.1) are uniformly bounded and ultimately bounded.

Proof. For $\forall \phi \in C([-\tau, 0], R_{+}^{5})$, we denote $f(t, \phi): =$

For $\forall \phi \in C([-\tau, 0], R_{+}^{5})$, it is clear that $f(t, \phi)$ is continuous and Lipschitzian with respect to $\phi$. Then model (2.1) has a unique solution through $(0, \phi)$. Since $f_i(t, \varphi)\geq 0$ for $\varphi \geq 0$ and $\varphi_i(0) = 0$, then $C([-\tau, 0], R_{+}^{5})$ is positively invariant set (see Remark 5.2.1 in ^[36]).

For human population $N_h(t)$,

It follows from the comparison principle that the solution of model (2.1) always exists for any time $t\geq 0$. In addition,

and

ensure that all solutions are ultimately bounded. Besides, when $N_h(t)\geq \max\{l_0, r_h/d_h\}$ and $N_v(t)\geq \max\{l_0, r_h/d_h\}$, the inequalities

hold, which means that all solutions of model (2.1) are uniformly bounded.

For system (2.1), there always exist the trivial equilibrium $(0, 0, N^*_h, 0, 0)$ and disease-free state $(N^*_m(t), 0, N^*_h, 0, 0)$, where $N^*_h = r_h/d_h$ and $N^*_m(t)$ is the positive periodic solution of model (2.2). For the trivial equilibrium, the disease is not transmitted when there is no mosquito, and then trivial equilibrium $(0, 0, N^*_h, 0, 0)$ will not be studied in this paper. Using the theory in ^[34], we compute the basic reproduction number of model (2.1). Linearizing model (2.1) at the disease-free periodic state $(N^*_m(t), 0, N^*_h, 0, 0)$, we have

Denote

and

Then model (2.3) becomes

with $x(t) = (M_i(t), H_i(t), H_r(t))^T$.

For periodic system

$Y(t, s)(t\geq s)$ represents the evolution operator. Then we have

For all $\omega-$periodic functions from $R$ to $R^5$, we define the Banach space $C_{\omega}$ and initial distribution of infectious mosquitoes and humans. Then the rate of new infectious mosquitoes and humans is $F(s)\phi(s)$, where $s$ is the time for infected mosquitoes and humans introduced. The distribution of infected individuals is $Y(t, s)F(s)\phi(s)$, where those infected mosquitoes and humans were infected at time $s$ and stay the infected classes at time $t\geq s$. Therefore, the distribution of accumulative new infectious mosquitoes and humans at time $t$ produced by infected mosquitoes and humans $\phi(s)$ brought at the previous time is

Linear operator $G: G_{\omega}\rightarrow G_{\omega}$ satisfies

Denote $\rho(G)$ as the spectral radius of $G$. According to ^[34], we obtain the next generation operator $G$ and ${\mathbb{R}_0} = \rho(G)$ is the basic reproduction number.

For system

we know that $W(t, \lambda)$ is the monodromy matrix and parameter $\lambda > 0$. Obviously, $-V(t)$ is cooperative and $F(t)$ is nonnegative, and then for $\lambda > 0$, $\rho(W(\omega, \lambda))$ is continuous and nonincreasing and $\lim_{\lambda\rightarrow \infty}\rho(W(\omega, \lambda)) < 1$. It follows from Theorem 2.1 in ^[34] that we have the following results.

Lemma 2.2. For the basic reproduction number, we have the following results:

$\bf(I)$ If there is $\lambda_0 > 0$ for $\rho(W(\omega, \lambda)) = 1$, then $\lambda_0$ is an eigenvalue of $G$ and $\mathbb{R}_0 > 0$.

$\bf(II)$ If $\mathbb{R}_0 > 0$, then $\rho(W(\omega, \lambda)) = 1$ has a unique solution $\lambda = \mathbb{R}_0$

$\bf(III)$ $\mathbb{R}_0 = 0$ if and only if $\rho(W(\omega, \lambda)) < 1,$ $\forall \lambda > 0$.

The Lemma 2.2 can be applied to compute the basic reproduction number $\mathbb{R}_0$ by numerical simulations ensures that $\mathbb{R}_0$ is in scale with $k(t)$.

3.   Threshold dynamics

In this section, we prove the global asymptotically stability of the disease-free periodic state when $\mathbb{R}_0 < 1$ and the uniform persistence of system (2.1) when $\mathbb{R}_0 > 1$.

Define

and

Obviously, $X_0$ is an open set. Suppose that (2.1) has an unique solution $u(t, \phi)$ with $u_0(\phi) = \phi$ and $\Phi(t)\varphi = u_t(\varphi)$. In addition, we assume that $P: X\rightarrow X$ is the Poincaré map corresponding to the model (2.1).

3.1.   Global stability of the disease-free periodic state

It follows from Theorem 2.2 in ^[34] that the disease-free periodic state is locally stable if $\mathbb{R}_0 < 1$ and unstable if $\mathbb{R}_0 > 1$. Let $X_K = C([-\tau, 0], [0, K]^5)$ for some $K\geq \max\{l_0, r_h/{d_h}\}$. Then we will show that the disease dies out when $\mathbb{R}_0 < 1$.

Theorem 3.1. If $\mathbb{R}_0 < 1$ and $(H)$ holds, then the disease-free periodic state of model (2.1) is globally asymptotically stable in $C([-\tau, 0], R_{+}^{5})\backslash \{(0, 0, N^*(h), 0, 0\}$.

Proof. It follows from Theorem 2.1 that $X_K$ is positively invariant. Then

Denote

From Theorem 2.2 in ^[34], $\mathbb{R}_0 < 1$ if and only if $\rho(\Phi_{F-V}(\omega)) < 1$. Since that the solutions are continuous in $\sigma$, then $\lim_{\sigma\rightarrow 0^+}\Phi_{M_{\sigma}}(\omega) = \Phi_{F-V}(\omega)$. Besides, we know that the spectrum for matrices is continuous ^[37], which implies that $\lim_{\sigma\rightarrow 0^+}\rho(\Phi_{M_{\sigma}}(\omega)) = \rho(\Phi_{F-V}(\omega))$. Therefore, we have $\rho(\Phi_{M_{\sigma}}(\omega)) < 1$ for some $\sigma$. From mosquito and human population, there is a $T_1 = T(\sigma)$ such that

Subsequently, we obtain

Using the Lemma 2.1 in ^[38], we obtain a positive $\omega$-periodic function $\hbar(t)$ satisfying that $h'(t) = M_{\sigma}(t)h(t)$ has a solution $h(t) = e^{\theta t}\hbar(t)$ and $\theta = \frac{1}{\omega}\ln\rho(\Phi_{M_{\sigma}}(\omega)) < 0$. The boundedness of the positive $\omega$-periodic function $\hbar(t)$ implies that $\lim_{t\rightarrow\infty}h(t)\rightarrow0$. Applying the comparison principle in ^[36], we obtain that $M_i(t)\rightarrow0, H_i(t)\rightarrow0$ and $H_r(t)\rightarrow0$ as $t\rightarrow\infty$. The theory of asymptotically periodic semiflow ^[39] gives that $M_s(t)\rightarrow N_m^{*}(t)$ and $H_s(t)\rightarrow N^*_h$ as $t\rightarrow\infty$. Therefore, the disease-free periodic state of model (2.1) is globally asymptotically stable in $C([-\tau, 0], R_{+}^{5})\backslash \{(0, 0, N^*_h, 0, 0\}$ when $\mathbb{R}_0 < 1$.

3.2.   Uniform persistence when $\mathbb{R}_0 > 1$

In this section, we show the uniform persistence of the disease when $\mathbb{R}_0 > 1$.

Theorem 3.2. When the assumption $(H)$ holds and $\mathbb{R}_0 > 1$, there exists an $\varepsilon > 0$ such that any solution of the system (2.1) in $C([-\tau, 0], R_{+}^{5})$ with $M_i(0) > 0, H_i(0) > 0$ and $H_r(0) > 0$ satisfies

In addition, system (2.1) has at least one positive periodic solution.

Proof. Since $P: X\rightarrow X$ is the Poincaré map corresponding to system (2.1), then $P(\phi) = u_{\omega}(\phi), \forall \phi \in X$, and then $\Phi(t)(X_0)\subset X_0, \forall t\geq 0.$ From Theorem 2.1, we obtain that $P^{n_0}$ is compact for all $n_0\omega > \tau$ and $P$ is point dissipative. According to the Theorem 2.9 in ^[40], $P$ generates a global attractor which is denoted as $A$.

Denote $E_1 = \{(0, 0, N^*_h, 0, 0)\}$ and $E_2 = \{(N^*_{m0}, 0, N^*_h, 0, 0)\}$, where $N^*_{m0}(\theta) = N^*_{m}(\theta), \forall \theta\in [-\tau, 0]$. There exists an $\delta_0 > 0$ such that $\inf_{t\geq 0}N^*_{m}(t) > 2\delta_0$. $\lim_{\phi\rightarrow E_1}(\Phi(t)\phi-E_1) = 0$ uniformly for $t\in [0, \omega]$ gives that we can choose an $\delta_1$ such that

We claim that

Suppose not, then there are some $\phi \in X_0$ such that $\limsup_{n\rightarrow\infty}\|\Phi(n\omega)\varphi-E_1\| < \delta_1$. Thus we choose an $T_2\geq 1$ such that $\|\Phi(t)\varphi-E_1\| < \delta_1, \forall n\geq T_2$. For $t-\tau\geq T_2\omega$, $t = n\omega+\bar{t}$, $n\geq T_2, \bar{t}\in [0, \omega]$ and $\|\Phi(t)\varphi-E_1\| = \|\Phi(\bar{t})\Phi(n\omega)\varphi-E_1\|\leq \delta_0$ hold. Then $M_s(t)\leq \delta_0, M_i(t)\leq \delta_0$ and $N_m(t)\leq \delta_0$ for $t-\tau\geq T_2\omega$. $N_m(0) = \varphi_1(0)+\varphi_2(0)$ implies that $\lim_{t\rightarrow\infty}(N_m(t)-N^*_m(t)) = 0$, which is a contradiction with $\inf_{t\geq 0}N^*_{m}(t) > 2\delta_0$.

Denote

Using the same discussion in Theorem 3.1, there is an $\sigma_1 > 0$ such that $\rho(\Phi_{\bar{M}_{\sigma}}(\omega)) > 1, \forall \sigma\in [0, \sigma_1]$. $\lim_{\phi\rightarrow E_2}(\Phi(t)\phi-\Phi(t)E_2) = 0$ uniformly for $t\in [0, \omega]$ implies that we can choose an $\varepsilon_1$ such that

The we claim that

Suppose not, then there are some $\varphi \in X_0$ satisfying $\limsup_{n\rightarrow\infty}\|\Phi(n\omega)\varphi-E_2\| < \varepsilon_1$. We choose an $T_3\geq 1$ such that $\|\Phi(n\omega)\varphi-E_2\| < \varepsilon_1, \forall n\geq T_3$. For $t-\tau\geq T_3\omega, t = n\omega+\bar{t}, n\geq T_3, \bar{t}\in [0, \omega]$, then $\|\Phi(t)\varphi-\Phi(t)E_2\| = \|\Phi(\bar{t})\Phi(n\omega)\varphi-\Phi(\bar{t})\Phi(n\omega)E_2\| = \|\Phi(\bar{t})\Phi(n\omega)\varphi-\Phi(\bar{t})E_2\|$, and

Lemma 2.1 in ^[38] implies that $\bar{h}(t) = e^{\theta t}\hbar(t)$ is a solution of equation $\bar{h}'(t) = M_{\sigma}(t)\bar{h}(t)$ for $\hbar(t)$ with period $\omega$, where $\theta = \frac{1}{\omega}\ln\rho(\Phi_{M_{\sigma}}(\omega)) > 0$. The boundedness of the positive $\omega-$periodic function $\hbar(t)$ gives that $\lim_{t\rightarrow\infty}\bar{h}(t)\rightarrow+\infty$. From $\Phi(t)\varphi\in X_0, \forall t\geq 0$, we can choose an $T_0\geq T_3$ and $\bar{\varepsilon}$ such that

The comparison principle ^[36] ensures that

Then $\lim_{t\rightarrow\infty}|(M_i(t), H_i(t), H_r(t))| = \infty$, which is a contradiction with Theorem 2.1.

Denote

Next we prove that $M_{\partial} = M_1\cup M_2$. $u_i(t, \varphi) = 0$, $\forall \varphi \in M_2$, $i = 1, 2$. Then $M_2\subset M_{\partial}$. For $\forall \varphi \in M_1$, $Z(t)\in (R_+, R^5_+)$ ensures that $Z_i(t) = 0, \forall t\geq 0$ for $i = 2, 4, 5$. Then

where $Z_1(0) = \varphi_1(0)$. Assume that $Z_3(t)$ satisfies

Then system (2.1) through $\varphi$ has a solution $Z(t)$. The uniqueness of solution implies that $u(t, \varphi) = Z_1(t), \forall t\geq 0$ and $M_1\subset M_{\partial}$. Therefore, $M_1\cup M_2 \subset M_{\partial}$. Next, we will prove $M_{\partial}\subset M_1\cup M_2$. Since $\sum_{j = 1}^2 \varphi_j(0) > 0$ for $\forall \varphi\in \partial X_0\backslash (M_1\cup M_2)$, then $\lim_{t\rightarrow\infty}|\sum_{j = 1}^2u_j(t, \varphi)-N^*_m(t)| = 0$. We can choose an $t_0 > 0$ such that $u_1(t, \varphi) > 0, \forall t > t_0$. Clearly, $u_3(t, \varphi) > 0, \forall t > 0$. If $\varphi_2(0) > 0$, then $u_j(t, \varphi) > 0$ for $i\in \{2, 4, 5\}, \forall t > 0$. Furthermore, if $\varphi_4(0) > 0$, then $u_2(t, \varphi) > 0, \forall t > t_0$. Therefore, $u(t, \varphi)\in X_0, \forall t > t_0$, which means that $\forall \varphi \in \partial X_0\backslash (M_1\cup M_2)$. There are some $n\omega > t_0$ such that $P^n(\varphi) \not\in \partial X_0$, and then $M_{\partial}\subset M_1\cup M_2$. Based on the above discussions, we obtain that for $P\in M_{\partial}$, $E_1$ and $E_2$ are disjoint, compact and isolated invariant sets, and then $\bar{A}_{M_{\partial}} = \bigcup_{\phi\in M_{\partial}}\omega(\phi) = \{E_1, E_2\}$. Besides, there is no cycle in $M_{\partial}$ for subset of $\{E_1, E_2\}$. Therefore, for $P\in X$, $E_1$ and $E_2$ are isolated invariant sets and $W^s(E_i)\cup X_0 = \emptyset (i = 1, 2)$, where $W^s(E_i)$ is the stable set of $E_i$. Note that acyclicity theorem on uniform persistence ^[39] gives the uniform persistence of $P$ and periodic semiflow $\Phi(t)$. According to Theorem 4.5 in ^[40], we obtain that there is an $\omega-$periodic solution $\Phi(t)\phi^*$ with $\phi^*\in X_0$ in model (2.1).

Using Theorem 4.5 in ^[40] with $\rho(x) = d(x, \partial X_0)$, we have that there is a compact global attractor $A_0$ for $P: X_0\rightarrow X_0$, and $A_0 = P(A_0) = \Phi(\omega)A_0$ ensures that $\varphi_i(0) > 0$ for $i = 2, 4, 5$, which means that $\phi_i(0) > 0$ for $i = 1, 3$ due to the invariance of $A_0$. Denote $D_0 = \bigcup_{t\in[0, \omega]}\Phi(t)A_0$, then $\varphi_i(0) > 0, \forall \varphi\in D_0, i = 1, 2$. Also, $D_0 \subset X_0$ and $\lim_{t\rightarrow \infty}d(\Phi(t)\phi, D_0) = 0, \forall \phi\in X_0$. We define $q: X\rightarrow R_+$ with $q(\phi) = \min_{1\leq i \leq 5}\{\phi_i(0)\}, \forall \phi \in X$ which is continuous. The compactness of $D_0 \in X_0$ implies that $\inf_{\phi\in D_0} q(\phi) = \min_{\phi\in D_0}q(\phi) > 0$. Then we can choose an $\varepsilon > 0$ such that

In addition, $\liminf_{t\rightarrow\infty}\min(\Phi(t)\phi^*)\geq \varepsilon$, and then $u_i(t, \phi^*) > 0, 1\leq i \leq 5, \forall t \geq 0$. Consequently, system (2.1) has at least one positive periodic solution $u(t, \phi^*)$.

When the basic reproduction number is greater than one, the disease persists in both mosquito and human population, and there has always infected mosquitoes and humans if they exist at the initial moment. In addition, periodic oscillations appear in both mosquito and human population.

4.   A case study

In this section, we study the outbreak of dengue fever in Guangdong Province, China in 2014. We verify our theretical results and explore the effect of maturation delay on the dengue virus transmission. The simulations are implemented by Matlab software.

China is hit by the outbreak of dengue fever since 1978, and then it spreads over Southern regions. There are serious infections of dengue in Guangdong after 1997 ^[9]. Cumulative number of reported cases of Guangdong province accounts for $94.99\%$ of the total reported cases in China, and all death cases are in Guangdong Province. Especially in 2014, Guangdong Province suffered from the worst dengue outbreak. This can be reflected by reported dengue fever cases in Guangzhou city. To carry out numerical simulations on the seasonal transmission pattern in Guangdong Province, parameter values will be estimated in model (2.1).

4.1.   Parameter estimates

From the National Bureau of Statistics of China in 2010, the number of people in Guangdong Province is 104,410,000. Here we assume that all people in Guangdong Province are susceptible. The average life expectancy of Chinese people was 74.83 years in 2010, and then the natural death rate of human is $d_h = \frac{1}{74.83\times12} = 0.001114$ $Month^{-1}$. Hence, the recruitment rate $r_h$ of human for Guangdong province can be calculated by $r_h = d_h\times104410000 = 116313$. The details of parameters are shown in Table 1.

Using the monthly mean temperature (National Centers for Environmental Information) observed in six monitoring stations in Guangdong Province from January 2003 to December 2012, we can calculate the average monthly temperature for Guangdong Province which is shown in Table 2. By the relationship between temperature, recruitment rate and biting pattern, we can estimate the parameters $r_m(t)$ and $k(t)$.

Based on the experiments in ^[44] on Aedes aegypti mosquitoes over the range of $10.54\; ^\circ$C, we choose the expressions of the intrinsic oviposition rate

where $C$ is the temperature in Celsius. In order to make sure the positivity of the intrinsic oviposition rate, we choose $C\geq 12\; ^\circ$C. The intrinsic oviposition rate is zero when $C < 12\; ^\circ$C. The intrinsic oviposition rate in Guangdong Province can be fitted by

The best fitted curve for the duration of the mosquito gonotrophic cycle is given by

Therefore, the temperature-dependent contact rate per unit time in Guangdong Province can be written as

We fit the biting rate in Guangdong province, and then we have

4.2.   Model validation

In this subsection, we use MATLAB software to simulate the dengue virus transmission in Guangdong Province based on the data from the Health Department of Guangdong Province, and obtain the monthly numbers of newly reported dengue cases from June 2015 to May 2017. The initial values are chosen as $M_s(0) = 241280, M_i(0) = 300, H_s(0) = 10000, H_i(0) = 15000$ and $H_r(0) = 98139900$. We take June 1 as the starting point. The simulations from June 2015 to May 2017 are presented in Figure 2 which shows the comparison between the monthly reported dengue cases and the model predicting cases for Guangdong Province. Our simulations are in accordance with the dengue cases reported in Guangdong Province.

4.3.   Long term behaviors

In this subsection, we show the long term behaviors of dengue virus transmission in Guangdong Province. Setting $M_s(\theta) = 24120, M_i(\theta) = 1800, H_s(\theta) = 3000, H_i(\theta) = 150, H_r(\theta) = 30$ for all $\theta\in[-\tau, 0]$. Firstly, we choose $\beta_m = 0.75, \beta_h = 0.75, \alpha = 0.0001, d_h = 0.001114, r_h = 116313, d_m = 0.07, d_j = 0.46, \tau = 0.8224, \mu = 0.1428$. Using Matlab software and Lemma 2.2, we get the basic reproduction number $\mathbb{R}_0 = 0.9981 < 1$. In this case, it follows from Figures 3 and 4 that both infectious human and infectious mosquito approach to zero eventually, which means that the dengue fever dies out. This demonstrates the theoretical results in Theorem 3.1.

Similarly, we take $d_m = 0.06, d_j = 0.4$, $\tau = 0.3947$, and all other parameters are the same as these in Figures 3 and 4. We can get $\mathbb{R}_0 = 1.1296 > 1$. In this case, it follows from Figures 5 and 6 that a positive periodic solution is observed and the dengue fever persists uniformly. The periodic fluctuations of the infectious humans (resp. mosquitoes) compartments are shown in Figure 5 (resp. Figure 6). Numerical simulations demonstrate that our simulations are consistent with the Theorem 3.2.

5.   Discussions and conclusions

Climatic factors contribute to the dengue virus transmission. Scientific evidences have shown that higher temperature will increase the proportion of infected mosquitoes, which is dangerous for human ^[45]. Hale et al. ^[18,46] found that climate change could increase the risk of dengue virus transmission, and about half of the global population will suffer from the risk of dengue infection in 2085 if the climate gets worse. Thus it is reasonable to study the impact of climatic factors such as temperature on dengue virus transmission.

Mathematical models can offer us an important method to explore the transmission and risk of dengue virus and help us to control the disease in time-varying environments. In this paper, a dengue model with maturation delay and seasonality is formulated to study the impact of temperature on the transmission of dengue virus, where the birth rate and biting rate of mosquitoes depend on the average temperature. This model gives a baseline against that climate change can be estimated in the long term. Applying the methods in ^[34], we have figured out the basic reproduction number $\mathbb{R}_0$. From the theoretical point of view, the disease-free periodic state $E_0$ is globally asymptotically stable when $\mathbb{R}_0 < 1$, which means that the infection can be cleared from the population if there is a small invasion. System admits at least one positive periodic solution when $\mathbb{R}_0 > 1$, that is, the dengue fever will persist and show seasonal fluctuations (Figures 5 and 6). For our model, we chose seasonable parameter values for modeling dengue fever cases in Guangdong province in China. The basic reproduction number of the model is computed numerically. For the period of June 2016 to May 2017, since $\mathbb{R}_0 = 0.9981 < 1$, then the disease will be cleared out eventually in Guangdong province (Figures 3 and 4), which is consistent with the data reported from June 2016 to May 2017.

If we only consider the delay effect and assume that mosquito growth rate $r_m(t)$ and biting rate $k(t)$ are time-independent, then the total mosquito population satisfies

For scalar Eq (5.1), it has been proven by Fan et al. ^[23] that $N_m(t) = 0$ is globally asymptotically stable if $r_me^{-d_j\tau}\leq d_m$, and $N_m(t) = \frac{1}{\alpha}\ln\left(\frac{r_m}{d_m e^{d_j\tau}}\right)$ is globally asymptotically stable if $r_me^{-d_j\tau} > d_m$. Therefore, the total mosquito population will not oscillate for any $\tau > 0$. For the full system (2.1), Eq (5.1) and the last four equations of system (2.1) are decoupled. Notice that the standard incidence rate and linear death rate are adopted in system (2.1). It has been proven by Shan et al. ^[47] that such incidence interaction and linear death will not generate oscillations. Hence, the delay alone will not induce any periodicity in our system.

In the present work, we only consider maturation delay as a constant, it should be more reasonable to consider maturation delay $\tau$ as a time-dependent function ^[48,49]. In fact, there are four development stages during a mosquito's lifetime which are egg, larva, pupa and adult. The dynamics of the first three aquatic stages have great impact on the spread of mosquito-borne diseases, and it will be more interesting to develop a structured dengue model with seasonality. We can also consider some effective measures to prevent and control the spread of dengue. These problems are left as our future work.

Our results will help us to understand dengue virus transmission between human and mosquitoes, and also be applied to other mosquito-borne disease such as West Nile virus, Zika and Chikungunya.

Acknowledgments

The authors thank the handling editor and anonymous referee for their very valuable suggestions, which greatly improve the manuscript. The research of H. Song was supported by the National Natural Science Foundation of China (11601291), and Program for the Outstanding Innovative Teams (OIT) of Higher Learning Institutions of Shanxi, and Shanxi Scholarship Council of China and Scientific and Technological Innovation Programs (STIP) of Higher Education Institutions in Shanxi. The research of C. Shan was supported by the Simons Foundation-Mathematics and Physical Sciences 523360.

Conflict of interest

All authors declare no conflicts of interest in this paper.