Modeling the Heterogeneity of Dengue Transmission in a City

Abstract

Dengue fever is one of the most important vector-borne diseases in the world, and modeling its transmission dynamics allows for determining the key influence factors and helps to perform interventions. The heterogeneity of mosquito bites of humans during the spread of dengue virus is an important factor that should be considered when modeling the dynamics. However, traditional models generally assumed homogeneous mixing between humans and vectors, which is inconsistent with reality. In this study, we proposed a compartmental model with negative binomial distribution transmission terms to model this heterogeneity at the population level. By including the aquatic stage of mosquitoes and incorporating the impacts of the environment and climate factors, an extended model was used to simulate the 2014 dengue outbreak in Guangzhou, China, and to simulate the spread of dengue in different scenarios. The results showed that a high level of heterogeneity can result in a small peak size in an outbreak. As the level of heterogeneity decreases, the transmission dynamics approximate the dynamics predicted by the corresponding homogeneous mixing model. The simulation results from different scenarios showed that performing interventions early and decreasing the carrying capacity for mosquitoes are necessary for preventing and controlling dengue epidemics. This study contributes to a better understanding of the impact of heterogeneity during the spread of dengue virus.

1. Introduction

Dengue fever (DF) is caused by four serotypes of dengue virus (DENV-1 to DENV-4), which can be transmitted by female mosquitoes. The major vector is Aedes aegypti, whereas Aedes albopictus acts as a secondary vector. Infection by one serotype of DENV induces lifelong immunity to the infecting serotype but only temporary partial immunity to the other three serotypes [1].

Dengue has become the most important arthropod-borne viral disease affecting humans [1], and its incidence has substantially increased worldwide in recent decades. It was estimated that there are 390 million dengue infections per year, of which 96 million manifest clinically (any level of disease severity) [2,3]. In another study, Brady et al. estimated that a total of 3.97 billion people in 128 countries are at risk of infection with DENV [3,4]. The drivers for the rapid dengue expansion include urbanization, globalization (such as travel and trade), lack of effective mosquito control, and climate change, among others [5,6].

Mathematical models have been used to model the transmission dynamics of dengue virus and to evaluate the effectiveness and/or cost effectiveness of interventions [7,8,9]. Most of these models were compartmental-type models, which generally assumed a homogeneous mixing for the mosquito bites of humans, i.e., each mosquito can bite any host with equal probability. However, this assumption was in conflict with real epidemics, which are highly heterogeneous with multiple levels: spatial heterogeneity and individual heterogeneity [10]. The spatial heterogeneity was largely governed by environmental conditions and ecological variations, causing poorly mixed mosquito-host encounters [9]. Specifically, the heterogeneous exposure of hosts to mosquito bites was due to factors such as proximity to the aquatic habitats of immature mosquitoes, the type of house [11], insecticide use [10], human movements [12,13], and so forth. The individual heterogeneity has been described in terms of heterogeneous biting, i.e., a given mosquito bites some hosts more than others [14], which was associated with human sex, age, size, pregnancy, defensive behavior, blood type, and variability in human odors, among others [10,11]. The impacts of heterogeneity on transmission dynamics have previously been addressed using meta-population models [15,16] and agent-based models [17,18,19]. However, this topic still deserves more attention [9].

The negative binomial distribution (NBD) has been widely used to represent count data, particularly for data exhibiting overdispersion. In the biological literature, NBD was used in different transmission dynamics models to explore the influence of transmission heterogeneity, including host–parasitoid models [20,21], insect–pathogen models [22], and so forth. Recently, the application of NBD was extended to epidemiological models. In a famous study, the authors integrated branching process theory and contact tracing data by assuming that the number of secondary cases resulting from each infection was described by an ’offspring distribution’ to analyze the influence of individual variations in infectiousness on disease emergence [23]. The results showed that the NBD ($R_0$, k), where the parameter k captures the skew in the transmission distribution, was favorable for a number of disease datasets. The values of k estimated from the data ranged from 0.01 to approximately 0.1, highlighting the large variance in individual infectiousness. Kucharski et al. used the same method as in [23] to examine the level of heterogeneity in the transmission of Middle East respiratory syndrome coronavirus in South Korea in 2015, and they found substantial potential for superspreading [24]. The NBD was also used to explore the spread of DF. Padmanabha et al. used a multinomial negative binomial model to describe the number of secondary infections to evaluate the overall effects of mosquito production and household human density on dengue transmission [25]. The simulation results showed that the intervention of super-productive containers was substantially more effective in reducing the risk of dengue at higher human densities. In addition, the NBD was also used in dynamics modeling for infectious diseases. Barlow used the NBD transmission function $kln(1+\frac{\beta I}{k})$ to model a possum-tuberculosis (TB) system, where k denotes the level of heterogeneity and $\beta$ is the transmission rate [26]. Kong et al. developed an SEIR model with an NBD transmission function, $kln(1+\frac{\beta I}{kN})S$, to model the heterogeneity of contact rate for direct infectious disease [27]. The influence of different transmission functions was studied by Hoch et al. [28].

In this paper, we modeled the heterogeneity of mosquito bites of humans using a compartmental model with NBD transmission functions. The impact of heterogeneity on the transmission dynamics of DF was explored through numerical simulations. Then, we extended the model by including the aquatic phase of mosquitoes and integrating the environment and climate impacts on them. The extended model was used to reproduce the 2014 DF outbreak in Guangzhou, China. Finally, we analyzed the impacts of different interventions on the transmission of DF through simulations in different scenarios.

2. Materials and Methods

2.1. Data

In 2014, an unanticipated severe DF outbreak occurred in Guangdong Province, with more than 45,000 cases of infection [29,30]. As the provincial capital of Guangdong Province, Guangzhou was the most affected city, with 37,420 reported local DF cases. In this paper, we chose Guangzhou as our study region, and its geographic location is shown in Figure 1. DF outbreaks in China were previously thought to be imported [29,31], although recent studies suggested that DF may be endemic to China [32,33]. Here, we only considered the local DF cases and considered them to be triggered by imported cases [31]. Figure 2 shows the daily new DF cases and cumulative cases. The data were collected based on the onset date in the case data obtained from the Chinese Center for Disease Control and Prevention (CDC, http://www.chinacdc.cn/). The first reported local DF case occurred on 11 June 2014, followed by sporadic cases until the middle of July; then, the number of cases gradually increased. A rapid increase in DF cases occurred in September, and the number of cases peaked in early October. After 11 October, the number of DF cases gradually decreased. The last case of DF was reported on 19 December 2014.

The main mosquito that transmitted dengue in Guangzhou, China, was Aedes albopictus [34,35], and the population of this species was influenced by climate. The temperature, precipitation and evaporation data from 2013 to 2014 were obtained from the China Meteorological Data Sharing Service System (http://data.cma.cn/). These data were monitored by meteorological stations located in Guangzhou and within 50 kilometers of the borders of Guangzhou (Figure 1). Then, the inverse distance weighted interpolation method [36] was employed to obtain the corresponding data for Guangzhou. The interpolation results (shown in Figure 3) were used as inputs to the following models.

2.2. The NBD SEIR-SEI Model

A number of mathematical models have been developed to understand the dynamics of dengue infection and to evaluate the effectiveness and/or cost effectiveness of control strategies [7,8]. Among these models, the SEIR(human)–SEI(vector) compartmental model has been widely used to model the transmission dynamics of a single serotype of DENV [7,8,37,38,39,40]. This type of model divides humans into susceptible, exposed, infectious and removed groups and divides mosquitoes into susceptible, exposed and infectious groups. DENV is passed to susceptible humans ($S_H$) through the bites of infectious female mosquitoes ($I_V$); the infected humans then enter the exposed compartment $E_H$. After an incubation period of 4 to 10 days [41], the patients enter compartment $I_H$ and can transmit the virus (for 4–5 days, maximum of 12 days [3]; in another study, 1–7 days, mean of 4.5 days [42]) via mosquito bites. Finally, the hosts become immune to that serotype and enter compartment $R_H$. Susceptible mosquitoes ($S_V$) are infected by biting infectious humans ($I_H$), and then they enter the exposed compartment $E_V$. After the extrinsic incubation period (EIP) of 8 to 12 days [41], infected mosquitoes are capable of transmitting the virus for the remainder of their life ($I_V$). Figure 4 presents the structure of this model, and

Table 1. Parameter notation, biological meaning, values and sources.

Parameter	Biological Meaning	Range	Value	Source
a	Average daily biting rate	0.3–1	1	[8]
$p_V$	Transmission probability from vector to human per bite	0.1–0.75	0.5	[8]
$p_H$	Transmission probability from human to vector per bite	0.5–1	0.75	[8]
$1/{\mu }_H$	Human life expectancy (years)	-	75	Assumed
${\mu }_V^{-1}$	Average lifespan of mosquitoes (days)	4–50	21	[8]
${\delta }_H^{-1}$	Intrinsic incubation period (IIP, days)	4–10	7	[41]
${\delta }_V^{-1}$	Extrinsic incubation period (EIP, days)	8–12	10	[8]
${\gamma }_H^{-1}$	Infectious period (days)	1–7	4.5	[42]

shows the biological meanings, ranges and values of the parameters used in the following simulation.

Most of these models assumed that the mosquito-host encounters were well mixed, i.e., any given host had an equal probability of being bitten by any given mosquito [7,8]. This assumption resulted in a well-mixed transmission function between infectious mosquitoes and susceptible human hosts, with the form:

$$\frac{a{p}_V{I}_V}{N_H}{S}_H$$

(1)

and that between infectious human hosts and susceptible mosquitoes:

$$\frac{a{p}_H{I}_H}{N_H}{S}_V.$$

(2)

However, this assumption is particularly unrealistic at large spatial scales (for example, a whole city): a given mosquito may only have the opportunity to bite a limited subset of hosts in the population, i.e., the numbers of bites of different hosts by any given mosquito differ. We assumed that the number of effective bites (bites causing the transmission of DENV) of the ith susceptible human from infectious mosquitoes, $X_i$, followed a Poisson distribution with parameter ${\theta }_i$. Unique values for all ${\theta }_i$ results in the homogeneous mixing model. Here, we assumed that the values of all ${\theta }_i$ were different from person to person and that it followed a Gamma distribution with shape parameter $k_1$ and a scale parameter. After some derivations [27], we obtained the NBD transmission function from infectious mosquitoes to susceptible humans:

$$k_{1}ln(1+\frac{a{p}_V{I}_V}{k_1{N}_H})S_H,$$

(3)

where $k_1$ characterized the level of heterogeneity of bites of susceptible humans by infectious mosquitoes.

Similarly, assuming that the number of effective bites of infectious hosts from the jth susceptible mosquito, $Y_j$, also followed a Poisson distribution with parameter ${\delta }_j$, which had a Gamma distribution with shape parameter $k_2$, we obtained the NBD transmission function from infectious humans to susceptible mosquitoes:

$$k_{2}ln(1+\frac{a{p}_H{I}_H}{k_2{N}_H})S_V,$$

(4)

where $k_2$ characterized the level of heterogeneity of bites of infectious humans by susceptible mosquitoes.

Because $k_1$ and $k_2$ both characterized the level of heterogeneity of mosquito bites of humans, we assumed that $k_1=k_2=k$ for simplicity. With the derived NBD transmission functions (Equations (3) and (4)), we constructed an NBD SEIR-SEI compartmental model to describe the transmission process of DENV. This model can be represented by a set of ordinary differential equations (Equation (5)), in which the human and vector populations were assumed to be constants:

$$\left\{\begin{array}{cc}\hfill \frac{\mathrm{d}{S}_H}{\mathrm{d}t}& ={\mu }_H{N}_H-kln(1+\frac{a{p}_V{I}_V}{k{N}_H})S_H-{\mu }_H{S}_H,\hfill \\ \hfill \frac{\mathrm{d}{E}_H}{\mathrm{d}t}& =kln(1+\frac{a{p}_V{I}_V}{k{N}_H})S_H-({\delta }_H+{\mu }_H)E_H,\hfill \\ \hfill \frac{\mathrm{d}{I}_H}{\mathrm{d}t}& ={\delta }_H{E}_H-({\gamma }_H+{\mu }_H)I_H,\hfill \\ \hfill \frac{\mathrm{d}{R}_H}{\mathrm{d}t}& ={\gamma }_H{I}_H-{\mu }_H{R}_H,\hfill \\ \hfill \frac{\mathrm{d}{S}_V}{\mathrm{d}t}& ={\mu }_V{N}_V-kln(1+\frac{a{p}_H{I}_H}{k{N}_H})S_V-{\mu }_V{S}_V,\hfill \\ \hfill \frac{\mathrm{d}{E}_V}{\mathrm{d}t}& =kln(1+\frac{a{p}_H{I}_H}{k{N}_H})S_V-({\delta }_V+{\mu }_V)E_V,\hfill \\ \hfill \frac{\mathrm{d}{I}_V}{\mathrm{d}t}& ={\delta }_V{E}_V-{\mu }_V{I}_V.\hfill \end{array}\right.$$

(5)

Because

$$\underset{k\to \infty }{lim}kln(1+\frac{a{p}_V{I}_V}{k{N}_H})S_H=\frac{a{p}_V{I}_V}{N_H}{S}_H,$$

$$\underset{k\to \infty }{lim}kln(1+\frac{a{p}_H{I}_H}{k{N}_H})S_V=\frac{a{p}_H{I}_H}{N_H}{S}_V,$$

the NBD transmission functions approximate that used in the corresponding homogeneous mixing models [37] (Equations (1) and (2)) when the parameter k approximates $+\infty$. Therefore, the homogeneous mixing models can be viewed as special cases of the NBD model (Equation (5)).

Based on the next-generation matrix operator approach [43], the basic reproductive number, $R_0$, was derived (see details in Appendix A):

$$R_0=\sqrt{\frac{N_V}{N_H}\frac{a^2{p}_H{p}_V{\delta }_H{\delta }_V}{{\mu }_V({\delta }_H+{\mu }_H)({\gamma }_H+{\mu }_H)({\delta }_V+{\mu }_V)}}.$$

It is identical to the $R_0$ in the corresponding homogeneous mixing models [38]. Note that it is independent of k, indicating that $R_0$ is not dependent on the level of heterogeneity. This result agrees with previous works [24,44], in which the authors inferred the maximum likelihood estimation of $R_0$ and the dispersion parameter k that characterized heterogeneity, and it was shown that the maximum likelihood estimation of $R_0$ was independent of k, although the associated confidence intervals do depend on it. The independence can be explained by $R_0$ being an average value fundamentally, neglecting heterogeneity. In this way, the disease outbreaks cannot be fully understood if individual variation in infectiousness is neglected [23], considering the transmission variance. Therefore, heterogeneity is indispensable for modeling transmission dynamics.

2.3. Extended Model

In reality, the number of mosquitoes is not constant but varies with the climate [45,46,47,48], environment [49], and human behavior [50]. Climate influences the transmission of DF by affecting mosquito development and mosquito/human interactions [45]. The entire lifecycle of the mosquito, including breeding, development and mortality, is influenced by temperature [45,46,47,48]. Temperature also influences virus replication and transmission. Studies have found that high temperature can increase the rate of viral replication within the vector and therefore shorten the EIP [41,51]. Precipitation also influences the mosquito population by providing habitats for their aquatic stages [45,47]. It has been shown that higher precipitation is associated with increased mosquito populations [47,52]. However, intense rainfall may wash out breeding sites and therefore have a negative effect on mosquito populations [45,52]. The mosquito population is also influenced by human behavior since man-made containers can provide habitats for the aquatic stages of mosquitoes [52,53,54,55].

When modeling dengue dynamics, it was necessary to consider the above factors. First, we extended the NBD SEIR-SEI model (Equation (5)) by explicitly modeling the aquatic immature life stages of mosquitoes. Four additional compartments, the eggs ($E_g$), larva ($L_a$), pupae ($P_u$) and the emerging females that are unable to reproduce ($A_e$), were added to the above model. Figure 5 presents the structure of the extended model, and the definition for each symbol in this figure is shown in

Table 2. Parameter notation, biological meaning and values used in the extended model.

Parameter	Biological Meaning	Values ${}^1$
$\sigma$	Diapause	1 in Mar. 15 to Oct. 25; 0 otherwise ([31])
$n_E$	Eggs per gonotrophic cycle (per female)	Temperature dependent
$f_{ag}$	1/duration for gonotrophic cycle (per day)	Temperature dependent
$f_E$	Egg development rate	Temperature and precipitation dependent
$f_L$	Larva development rate	Temperature and precipitation dependent
$f_P$	Development rate of pupae to emerging adults	Temperature dependent
${\mu }_E$	Egg mortality rate	0.05 ([47])
${\mu }_L$	Mortality for larva	Temperature and density dependent
${\mu }_P$	Mortality for pupa	Temperature dependent
${\mu }_A$	Mortality rate during adult emergence	0.1 ([47])
${\mu }_V$	Mortality rate of adult mosquitoes	Temperature dependent
$\xi$	Sex ratio of Aedes albopictus at emergence	0.5 ([61])
${\gamma }_A$	Development rate of emerging adults (day${}^{-1}$)	0.4 ([47])
$K_L$	Carrying capacity of mosquito larvae population	Precipitation and environment dependent
$K_P$	Carrying capacity of mosquito pupae population	Precipitation and environment dependent
${\pi }_{max}$	Maximum carrying capacity for immature mosquitoes	Environment and density dependent
$N_V$	Total population of adult female mosquitoes	$N_V=S_V+E_V+I_V$

${}^1$ The details for the climate- and/or environment-dependent functions are shown in Appendix B.

. Each life-stage was represented by an ordinary differential equation (the fifth to the eighth equation in Equation (6)). The adults that are able to mate and feed ($S_V$, $E_V$ and $I_V$) can oviposit and begin their life cycle [56]. Density-dependent mortality was assumed at the larval stage [47,57] because inter- and intra-specific competition (e.g., for food) at the larvae stage usually limits mosquito populations [58]. Because mosquito pupae do not feed and rely solely on energy stored from the larval stage [59], the proportion of emerging pupae that survives to emergence transits to stage ’emerging adults’. Pupa density-dependent success of adult emergence was assumed [47,57] as emergence success has been shown to be negatively correlated to pupa density [52,60]. Second, the constant parameters used in Equation (5), which were influenced by climate and the environment, were replaced by climate- and/or environment-dependent functions. These parameters included the oviposition rate, mortality and development rate of mosquitoes; the biting frequency; and EIP. The impacts of the environment and human activities on the mosquito population were characterized by the carrying capacity for immature mosquitoes (${\pi }_{max}$). It was related to the water level, which increased due to rainfall. However, if there was an extreme rainfall event and the water level was close to the maximum water level, then a spillover effect was triggered, in which a fraction of the immature mosquitoes were washed out from their breeding sites [19,31]. The form and coefficients of the climate- and/or environment-dependent functions were based on [19,31,47]. The detailed forms of these functions can be found in Appendix B. The extended model was represented by the set of differential equations shown in Equation (6):

$$\left\{\begin{array}{cc}\hfill \frac{\mathrm{d}{S}_H}{\mathrm{d}t}& ={\alpha }_H{N}_H-kln(1+\frac{b{p}_V{I}_V}{k{N}_H})S_H-{\mu }_H{S}_H,\hfill \\ \hfill \frac{\mathrm{d}{E}_H}{\mathrm{d}t}& =kln(1+\frac{a{p}_V{I}_V}{k{N}_H})S_H-({\delta }_H+{\mu }_H)E_H,\hfill \\ \hfill \frac{\mathrm{d}{I}_H}{\mathrm{d}t}& ={\delta }_H{E}_H-({\gamma }_H+{\mu }_H)I_H,\hfill \\ \hfill \frac{\mathrm{d}{R}_H}{\mathrm{d}t}& ={\gamma }_H{I}_H-{\mu }_H{R}_H,\hfill \\ \hfill \frac{\mathrm{d}{E}_g}{\mathrm{d}t}& =n_E{f}_{ag}{N}_V-\sigma f_E{E}_g-{\mu }_E{E}_g,\hfill \\ \hfill \frac{\mathrm{d}{L}_a}{\mathrm{d}t}& =\sigma f_E{E}_g-f_L{L}_a-{\mu }_L(1+\frac{L_a}{K_L})L_a,\hfill \\ \hfill \frac{\mathrm{d}{P}_u}{\mathrm{d}t}& =f_L{L}_a-f_P{P}_u-{\mu }_P{P}_u,\hfill \\ \hfill \frac{\mathrm{d}{A}_e}{\mathrm{d}t}& =\xi f_P{P}_{u}exp(-{\mu }_A(1+\frac{P_u}{K_P}))-{\gamma }_A{A}_e-{\mu }_V{A}_e,\hfill \\ \hfill \frac{\mathrm{d}{S}_V}{\mathrm{d}t}& ={\gamma }_A{A}_e-kln(1+\frac{b{p}_H{I}_H}{k{N}_H})S_V-{\mu }_V{S}_V,\hfill \\ \hfill \frac{\mathrm{d}{E}_V}{\mathrm{d}t}& =kln(1+\frac{a{p}_H{I}_H}{k{N}_H})S_V-({\delta }_V+{\mu }_V)E_V,\hfill \\ \hfill \frac{\mathrm{d}{I}_V}{\mathrm{d}t}& ={\delta }_V{E}_V-{\mu }_V{I}_V.\hfill \end{array}\right.$$

(6)

2.4. Model Calculation

To reproduce the 2014 DF outbreak in Guangzhou, we calibrated the extended model (Equation (6)) to the daily new DF cases data of that outbreak. A number of parameters were estimated using the GlobalSearch algorithm in the MATLAB Global Optimization Toolbox [62,63]. The normalized root mean square error (NRMSE) [64] was chosen as the goodness of fit between the model output and the referenced daily new DF cases, as well as the objective function of the model fitting procedure. NRMSE varied from $-\infty$ (poor fit) to 1 (perfect fit) [64].

The extended model developed here is a type of deterministic mathematical dynamics model. Therefore, the introduction of an infectious individual to the original disease-free population would result in an epidemic outbreak [65]. During the initial phase of the 2014 DF outbreak in Guangzhou, only sporadic cases occurred. Each of these cases had the chance to infect mosquitoes and then trigger the outbreak, while he or she also had the possibility of infecting no mosquitoes, thereby breaking the chain of transmission. To our best knowledge, no study showed which case triggered the overall outbreak. To better reproduce the overall outbreak, i.e., to better fit the daily new DF cases shown in Figure 2, we ignored this initial phase and only fit the data from 15 July 2014, to the end of the outbreak. We chose this period because from that day on, the daily new DF cases showed a clear increasing trend (Figure 2). To describe the initial trigger, we employed an unknown parameter, $t_0$, to denote the date on which a trigger individual was infected. This was realized by setting $I_H\left(t\right)=0$ when $t<t_0$ and $I_H\left(t_0\right)=1$.

The 2014 DF outbreak in Guangzhou peaked on October 1, with 1613 new cases. After that day, the daily new DF cases showed a notable decrease (Figure 2). This decrease in cases may have been caused by the interventions, including larval breeding eradication, killing of adult mosquitoes with pesticides, public health education and community involvement [66]. Here, we separated the entire transmission process into two phases. In the first phase, an exponential growth phase was observed, during which interventions played no or a negligible role. When the interventions began to take effect, the transmission process entered the second phase. The date separating the two phases was set to 23 September 2014, using the method presented in [67]. The mortality of the mosquitoes must be different in the two phases. However, it was difficult to estimate the mortality caused by interventions. Here, rather than estimating the mortality, we modeled the two phases with different heterogeneity levels, reflected by different parameters, $k^1$ and $k^2$ ($k^2<k^1$). The smaller $k^2$ reflected a higher degree of heterogeneity in the second phase. We can image that the mosquitoes eradicated by human interventions were still “alive”, but they lost the ability to bite humans. Therefore, fewer people had the chance of being bitten, and the level of heterogeneity increased. In the model, this was realized by setting $k=k^1$ before September 23, and after that day, setting $k=k^2$.

To solve Equation (6), initial values for all the compartments were needed. Since the mosquito population for only the first simulated year was affected by the initial value for eggs [31,47,68], we ran the model over the period 2013 to 2014, although our focus was to simulate the dengue outbreak in 2014. The year 2013 was used to achieve a mosquito population for 2014 that was not sensitive to the initial conditions. Here, we set $E_g\left(0\right)=1000$ and $L_a\left(0\right)=1$, and the initial values of all the other mosquito compartments were set to be zero. Vertical transmission was not considered here; thus, we ignored the DF cases in 2013 for simplicity. We set $S_H\left(0\right)=12.8389\times {10}^6$, which was the size of the permanent population at the end of the year 2012, and $E_H\left(0\right)=I_H\left(0\right)=R_H\left(0\right)=0$.

During the transmission process of DF, the mild or asymptomatic infections cannot be ignored [69]. It was estimated that the overall inapparent-to-symptomatic (I:S) ratio was 2.2:1 (95% CI: 1.1-4.2:1) for DF in Zhongshan City [69], which is adjacent to Guangzhou. For simplicity, we used the daily DF cases multiplied by 3.2 to fit the model output.

In addition to the three unknown parameters introduced above ($t_0$, $k^1$ and $k^2$), five more unknown parameters were estimated: (1) the transmission probability from vector to human per bite ($p_V$), (2) the transmission probability from human to vector per bite ($p_H$), (3) the reciprocal of the intrinsic incubation period (${\delta }_H$), (4) the reciprocal of the infectious period for humans (${\gamma }_H$), and (5) the carrying capacity for immature mosquitoes (${\pi }_{max}$). The ranges for parameters (1)–(4) are shown in Table 1. Eight parameters in total were estimated during the calibration process. Some other parameters that were used are shown in Appendix B.

3. Results

3.1. Transmission Dynamics with Different Heterogeneity Levels

Through numerical simulations, we explored the influence of the heterogeneity of mosquito bites of humans on the transmission dynamics. We simulated the NBD SEIR-SEI model (Equation (5)) with different values for k: $0.1$, $0.01$, $0.001$, and $0.0001$. As the value of k decreases, the heterogeneity level increases. For all simulations, we assumed an equal amount of humans and mosquitoes, i.e., $N_H=N_V$. We set the initial conditions as follows: $N_H={10}^6$, $I_H\left(0\right)=1$, $E_H\left(0\right)=R_H\left(0\right)=0$, $S_H\left(0\right)=N_H-E_H\left(0\right)-I_H\left(0\right)-R_H\left(0\right)$, $S_V\left(0\right)=N_V$, and $E_V\left(0\right)=I_V\left(0\right)=0$. The parameter values were set to the values in Table 1. The simulation results are shown in Figure 6. For comparison, the epidemic curve predicted by the corresponding homogeneous mixing model was also drawn.

As the level of heterogeneity increased (the value of k decreased), the gap between the infection curves from the NBD SEIR-SEI model and that from the homogeneous mixing model increased and the peak size decreased. We explained these results by imagining a scenario with a high level of heterogeneity, in which only a small proportion of people had the chance of being bitten by infectious mosquitoes and being infected while most of people avoided being bitten and infected. This resulted in a smaller transmission rate and a smaller peak size compared with the output from the homogeneous mixing model. In real DF outbreaks, human hosts and mosquitoes are heterogeneously mixed. A higher degree of heterogeneity resulted in a higher extinction probability of diseases [23].

3.2. Model Calibration Results

Figure 7 (A) shows the reported daily new DF cases and the fitted epidemic curve. The value of NRMSE, which determined the goodness of fit, was 0.7299. The optimized values of the unknown parameters are listed in

Table 3. Parameter notation, biological meaning, and optimal values.

Parameter	Biological Meaning	Optimized Value ± Std ${}^1$
$t_0$	The date on which an infectious human trigged the outbreak	17 June 2014 ± 7 days
$k^1$	Heterogeneity level in the 1st phase	$3.8263\times {10}^{-4}\pm 1.1316\times {10}^{-6}$
$k^2$	Heterogeneity level in the 2nd phase	$2.1379\times {10}^{-7}\pm 6.6802\times {10}^{-9}$
$p_V$	Transmission probability from vector to human per bite	0.414 $\pm 0.012$
$p_H$	Transmission probability from human to vector per bite	0.682 $\pm 0.019$
${\pi }_{max}$	Maximum carrying capacity of immature mosquitoes	$5.1354\times {10}^6$± 14665
${\delta }_H^{-1}$	Intrinsic incubation period (IIP, days)	7.853 ± 0.22
${\gamma }_H^{-1}$	Infectious period (days)	4.746 ± 0.135
	Estimated cumulative reported cases	37,733

${}^1$ std refers to the standard deviation of each parameter through 100 runs.

. Because of the intrinsic stochasticity of the GlobalSearch algorithm [62,63], slightly different results may be obtained when running the programs. With random values within the range of the parameters (Table 1 and assumed for the others), we ran the optimal program multiple times and obtained the standard derivation of each parameter.

Both $k^1$ and $k^2$ were small, indicating that there was high heterogeneity during the entire outbreak period. The large difference between $k^1$ and $k^2$ ($k^2<k^1$) indicated that the heterogeneity level in the second phase was considerably higher than that in the first stage. We believe that this result may be the combined effect of human interventions, climate conditions and enhancement of human protection.

3.3. Simulation of Dengue Transmission in Different Scenarios

We simulated the impact of interventions on the spread of DF in different scenarios. The epidemic curves were plotted by changing one parameter while keeping the others identical to the values fitted to the 2014 DF outbreak in Guangzhou.

First, we assumed that the carrying capacity for immature mosquitoes (${\pi }_{max}$) was increased or decreased by 10% and then simulated the DF spread in these two scenarios. The simulated epidemic curves are shown in Figure 7B. The simulation results showed that the outbreak would reach the peak at the same time but with different peak sizes: an increased carrying capacity resulted in a higher peak size, whereas a decreased carrying capacity resulted in a lower peak size. The total number of infections was also influenced, increased by 63.47% or decreased by 50.29%, corresponding to ${\pi }_{max}$ increasing or decreasing, respectively (shown in

Table 4. The total number of infections reported and simulated in assumed scenarios.

Scenario	Reported Total Infections ${}^1$	Change ${}^2$
Reported	37,420	-
Fitted	37,733	-
${\pi }_{max}$ decreased by 10%	18,758	Decreased by 50.29%
${\pi }_{max}$ increased by 10%	61,684	Increased by 63.47%
Intervention 10 days earlier	14,831	Decreased by 60.69%
Intervention 10 days delayed	84,746	Increased by 124.59%
Vaccinating 10%	21,232	Decreased by 43.73%
Vaccinating 20%	11,563	Decreased by 69.36%

${}^1$ The reported total infections were calculated by multiplying the total infections predicted by the models by a constant 3.2 according to the I:S ratio [69]; ${}^2$ The rates of change were calculated based on the fitted total infection.

).

Second, we assumed that the beginnings of the interventions were set earlier or delayed by 10 days. Here, we assumed that the much smaller $k^2$ compared with $k^1$ mainly resulted from human interventions. We estimated that the date that separated the two phases was 23 September for the 2014 DF outbreak in Guangzhou. To simulate the impact of the timing of interventions, we set the timing of the switch from $k^1$ to $k^2$ to be September 13 and October 3, respectively. The simulated epidemic curves are shown in Figure 7C. Both the peak time and the total number of infections were changed. The transmission was observed to peak early with a lower peak size when the timing of the switch from $k^1$ to $k^2$ was set earlier by 10 days. In contrast, if the timing was delayed by 10 days, the daily infected cases would continue to increase until it reached a much higher peak (Figure 7C). The total number of infections in the two scenarios decreased by 60.69% and increased by 124.59%, respectively (Table 4).

Third, we assessed the impact of dengue immunization programs if implemented in Guangzhou in the future. The dengue vaccine Dengvaxia (also referred to as CYD-TDV) has been licensed and registered in several countries [70]. Here, we assumed that a dengue vaccination program had been implemented in Guangzhou and that, among the population, proportions of 10% and 20% had been vaccinated and were immune to DENV. We simulated the transmission in these two scenarios by decreasing the initial value of the susceptible humans, $S_H\left(0\right)$, by 10% and 20%, respectively, and increased the initial value of the immune humans, $R_H\left(0\right)$, by the corresponding proportion. The simulated epidemic curves are shown in Figure 7D. The simulation results showed that, as the proportion of the effectively vaccinated population increased, the peak size became smaller and the total number of infections decreased (Table 4). In practice, it is necessary to consider the cost efficacy of the vaccination program.

4. Discussion

In this paper, a mathematical model with nonlinear NBD transmission functions was developed to analyze the impact of heterogeneity of mosquito bites of humans during the dengue transmission process. The simulation results suggested that a higher level of heterogeneity can result in a small peak size and a smaller number of total infections. The simulations in different scenarios suggested introducing interventions early and clearing the mosquito breeding environment to prevent and control the transmission of DF.

It has been emphasized in many previous studies that the heterogeneity should be taken into account to model the transmission of DF [9,71,72]. Studies have indicated that people, rather than mosquitoes, influence the spatial dynamics of dengue virus [71,73]. Favier et al. modeled the contact heterogeneity using a stochastic SEIR-SEI model and the probability for hosts visiting others in another house [71]. Cosner et al. used a discrete-time multi-patch model to model the effects of the mobility of people on DF transmission dynamics [15,16]. Agent-based, spatially explicit models were applied to understand the factors contributing to DF spread, to explore the impacts of mosquito movement and distribution on DF transmission, to find optimized dengue control strategies, and to estimate the effectiveness and cost effectiveness of future vaccination programs [11,18,19,74].

The NBD SEIR-SEI model developed in this study was novel in terms of the use of the NBD transmission functions, which characterized the heterogeneity of mosquito bites of humans by assuming that the encounters between humans and mosquitoes were dispersed. It is known that most female Ae. aegypti may spend their lifetime in or around the houses where they emerge as adults [73,75]. Thus, for a disease-free population, a person infected with DENV can only spread DENV to the mosquitoes near the position he or she had appeared. Similarly, only a small portion of humans had the chance of being infected by one specific infectious mosquito. Therefore, the rate of DF transmission was not as fast as that predicted by the homogeneous mixing models, nor was the transmission intensity. Whereas the homogeneous mixing model always predicted that almost all people were infected by the dengue virus after an epidemic [37], the NBD SEIR-SEI model can reproduce a real dengue epidemic well. The very small fitted $k^1$ and $k^2$ values indicated that it was indeed highly heterogeneous. Therefore, it was necessary to consider the influence of heterogeneity when modeling the dynamics of DF.

The extended model including the aquatic phase of immature mosquitoes and integrating the influence of the environment and climate factors on them was used to reproduce the 2014 DF outbreak in Guangzhou. The fitted results (Figure 7A) demonstrated the ability of the model to describe the real DF transmission processes. It was also used to simulate dengue transmission dynamics in different scenarios and to evaluate the effects of different interventions quantitatively. The simulation results showed that performing interventions as early as possible and decreasing the mosquito carrying capacity were necessary to prevent and control DF outbreaks. Therefore, two important issues should be addressed. The first is to develop an early warning system to predict the risk of a dengue epidemic. Once the risk is higher than a predefined ‘alert’ threshold, intervention strategies should be started and performed in a timely manner. The other issue is to decrease the carrying capacity of aquatic phase of mosquitoes, for example, by cleaning water containers and other egg-laying habitats.

As we know, the more sophisticated the model is, the more data that need to be collected and the higher the requirements on the parametrization. Compared to the meta-population and agent-based models [11,15,16,18,19,74], the NBD SEIR-SEI model was considerably more tractable in that it used only one parameter to characterize the heterogeneity of host-vector encounters. However, the modeling assumptions deserve further exploration. In future research, it would be interesting to integrate the NBD model with other methods for a greater understanding of the DF’s spatiotemporal diffusion.

5. Conclusions

DF in Guangdong Province in China has been a major public health concern. Understanding the spread of DF is key to risk assessment and performing preventive interventions. In this study, we proposed a dynamic mathematical model with NBD transmission functions to analyze the influence of the heterogeneity of mosquito bites of humans during the DF transmission process. The results showed a smaller peak size and a smaller number of total infections with a higher level of heterogeneity. These findings contribute to a better understanding of DF transmission dynamics. The presented model can also be modified to describe the transmission dynamics of other vector-borne diseases.