Mathematics of dengue transmission dynamics: Roles of vector vertical transmission and temperature fluctuations.

A new deterministic model is designed and used to gain insight into the effect of seasonal variations in temperature and vector vertical transmission on the transmission dynamics of dengue disease. The model, which incorporates (among many other features) the dynamics of the immature dengue-competent mosquitoes, vertical transmission in the vector population, density-dependent larval mortality and temperature effects, is rigorously analysed and simulated using data relevant to the disease dynamics in Chiang Mai province of Thailand. The non-trivial disease-free equilibrium of the model is shown to be globally-asymptotically stable when the associated basic reproduction number of the model is less than unity. Numerical simulations of the model, using data relevant to the disease dynamics in the Chiang Mai province of Thailand, show that vertical transmission in the vector population has only marginal impact on the disease dynamics, and that the effect of vertical transmission is temperature-dependent (in particular, the effect of vertical transmission on the disease dynamics increases for values of the mean monthly temperature in the range [ 16 - 28 ] ∘ C, and decreases with increasing mean monthly temperature thereafter). It is further shown that dengue burden (as measured in terms of disease incidence) is maximized when the mean monthly temperature is in the range [ 26 - 28 ] ∘ C (and dengue burden decreases for mean monthly temperature values above 28 ∘ C). Thus, this study suggests that anti-dengue control efforts should be intensified during the period when this temperature range is recorded in the Chiang Mai province (this occurs between June and August).

Introduction

Dengue, a viral disease caused by one of the four closely-related Flavivirus (DENV 1–4), is endemic in many tropical and subtropical regions of the world (with over 2.5 billion people at risk of acquiring dengue infection) (world health Organization, 2017). Annually, the disease account for approximately 50 million cases and fatalities (Ong, Sandar, Chen, & Sin, 2007; world health Organization, 2017). Dengue and dengue haemorrhagic fever are on the rise in the Americas (Gubler & Trent, 1994; Pan American Sanitary Bureau, 1995). In Latin America, about of the population (around 81 million people) live in urban areas, and the incidence of the rise has been on the increase in the past decade (Githeko, Lindsay, Confalonieri, & Patz, 2000; Pan American Sanitary Bureau, 1994). In Puerto Rico, for example, almost 10,000 dengue fever cases are reported annually (dengue outbreaks are recorded in almost all Caribbean countries and Mexico (Gubler & Trent, 1994; Pan American Sanitary Bureau, 1995)). Dengue has also been periodically endemic in Texas over the past 20 years (Gubler & Trent, 1994; Pan American Sanitary Bureau, 1995). The disease, which is transmitted to humans by female Aedes aegypti mosquito (following taking a blood meal, needed for eggs laying), threatens other non-endemic countries in Europe. For instance, the first local transmission of the disease in France and Croatia was recorded in 2010 (Novoselet al, 2015; Rogers & Hay, 2012) (outbreaks were also recorded in Madeira islands of Portugal in 2012, imported cases (mainly from Portugal) were also detected in three other European countries (Rogers & Hay, 2012; world health Organization, 2017)). Furthermore, in 2013, dengue outbreaks were recorded in Miami, USA and Yunnan province of China (Dick et al., 2012; Zhanget al., 2014). Dengue causes life-threatening complications (such as Dengue Haemorrhagic fever and Dengue Shock syndrome (Halide & Ridd, 2008)), often triggered by immune responses to secondary infections (Vaughn, 2000).

The incidence of dengue has significantly increased globally over the last few years (Dick et al., 2012; Rogers & Hay, 2012; world health Organization, 2017). This is due to a number of factors (Githeko et al., 2000) (notably the geographic expansion, enhanced transmission intensity in endemic areas, variability in local weather and habitat conditions). Furthermore, vertical (transovarial) transmission, which has been observed in dengue transmission dynamics (Cosner et al., 2009; Pacheco, Esteva, & Vargas, 2009), is believed to retain dengue viral disease in nature during inter-epidemic periods of dengue (Angel & Joshi, 2008).

Changes in local temperature is known to significantly affect the dynamics of vector-borne diseases, including dengue (Githeko et al., 2000; Watson, Zinyowera, & Moss, 1998). In particular, temperature variability affects the maturation, survival, biting rate and abundance of dengue-competent mosquitoes (Githeko et al., 2000). As the global temperature is increasing due to greenhouse-gas effects (daily average temperature in southern borders of USA have increased by C over the past 30 years (Karl & Plummer, 1995); and it is estimated that global temperature will rise by C over the next 100 years (Githeko et al., 2000; Watson et al., 1998)), it is imperative to carry out detailed modeling studies to analyse the potential impact of such increases on the dynamics of vector-borne diseases (it should be mentioned that health risks due to these climate changes differ between countries, depending on the level of infrastructure and economic development (Karl & Plummer, 1995)).

Vertical transmission (i.e., disease transmission from an infected mother to a child) is also another factor affecting the dynamics of many pathogens and diseases including dengue (Adams & Boots, 2010; Angel & Joshi, 2008; Coutinho, Burattini, Lopez, & Massad, 2006; Joshi et al., 2006; Kow, Koon, & Yin, 2001). For instance, evidence for vertical transmission has been established in the dynamics of vector-borne diseases such as, La Crosse virus (Miller, Defoliart, & Yuill, 1978), St Louis Encephalitis virus ((Nayar, Rosen, & Knight, 1986)), West Nile virus (Baqar, Hayes, Murphy, & Watts, 1993) and Yellow fever (Diallo, Thonnon, & Fontenille, 2000). Vertical transmission of dengue virus has been demonstrated in the lab in Aedes aegypti, Aedes albopictus and Aedes scutellaris mosquitoes (Freier & Rosen, 1987; Joshi, Mourya, & Sharma, 2002; Mitchell & Miller, 1990; Rosen, Shroyer, Tesh, Freier, & Lien, 1983; Shroyer, 1990) (including in the wild (Angel & Joshi, 2008; Joshi et al., 2006; Kow et al., 2001)). It should also be mentioned that local temperature affects vertical transmission in the vector population. In subtropical regions, for example, dengue disease shows a resurgent pattern with yearly epidemics (which starts typically in the months with heavy rains and heat, peaking some 3 or 4 months after the beginning of the rainy season) (Monath Heinzet al, 1996). In the dry months, the number of dengue cases typically drops essentially to zero (because the disease (vector) has virtually disappeared during this period) (Monath Heinzet al, 1996). Dengue continues to re-emerge for many years in some regions (Coutinho et al., 2006; Githeko et al., 2000). This (re-emergence) is attributed to many factors, such as the survival of long-lived infected adult female mosquitoes, infected eggs (laid by infected adult female mosquitoes, vertically) that remain infected during the dry season (and their hatching during the beginning of the raining season) (Coutinho et al., 2006; Githeko et al., 2000; Monath Heinzet al, 1996).

Although numerous climate variables (such as temperature, precipitation and humidity (Chenet al., 2012; Githeko et al., 2000; Karl & Plummer, 1995; Lafferty, 2009)) affect the transmission dynamics of vector-borne diseases, the current study focuses on the singular effects of temperature. Consequently, the purpose of the current study is to use mathematical modeling approaches to gain insight into the role of temperature variability, and vertical transmission in the vector population, on the transmission dynamics of dengue disease in the community. A new deterministic model, which includes the dynamics of both the immature (i.e., modeling the eggs-larvae-pupae lifecycle stages) and adult female Aedes aegypti mosquitoes, as well as the effect of vector vertical transmission and temperature variability, will be designed. The paper is organized as follows. The model is designed in Section 2. The case of the model where temperature is fixed (i.e., the autonomous equivalent of the non-autonomous model) is rigorously analysed (for the existence and asymptotic stability of some of its equilibria, as well as to characterize bifurcation types) in Section 3. Uncertainty and sensitivity analysis of the parameters of this version of the model, as well as numerical simulations of the effect of temperature variability, are also carried out. The (full) non-autonomous model, which accounts for daily fluctuations in local temperature, is rigorously analysed in Section 4.

Model formulation

This study is motivated by dengue transmission dynamics in Chiang Mai province of Thailand (City populationCHIANG MAI, 2017; Pongsiriet al., 2012; Thai Meteorological Department, 2017). Although there are four dengue serotypes (DENV 1–4), serological data from this province (for 2004 to 2010) shows that one of the four serotypes typically dominates the others each year (Pongsiriet al., 2012) (for example, in the year 2004, the percentage prevalence of DENV-1, DENV-2, DENV-3 and DENV-4 were , , and , respectively (Pongsiriet al., 2012)). Consequently, this study will consider only one dengue serotype in the community (this simplifying assumption allows for the tractability of the mathematical analysis to be carried out; a number of models for dengue transmission dynamics also use a single dengue serotype, such as those in (Adams & Boots, 2010; Coutinho et al., 2006; Esteva & Vargas, 2000; Garba, Gumel, & Bakar, 2008; Grunnill & Boots, 2016; Halide & Ridd, 2008)).

The model to be designed is based, first of all, on splitting the total immature mosquito population at time t (denoted by ) into compartments of susceptible eggs (), infected eggs (), susceptible larvae (), infected larvae (), susceptible pupae () and infected pupae (), so that

Furthermore, the total adult female mosquito population at time t (denoted by ) is split into the sub-populations of susceptible adult female mosquitoes () and infected adult female mosquitoes . Thus,

Finally, the total human population at time t (denoted by ) is sub-divided into susceptible (), exposed (), symptomatic () and recovered () humans. Hence,

The model to be designed incorporates the effect of variability in ambient (air) temperature (denoted by ) on the dynamics of the mosquito-borne disease in a community. It is assumed that is non-negative, continuous and bounded periodic functions of t (it is also assumed, for mathematical convenience, that air temperature and the temperature near the surface of the water are approximately the same; so that we can use to approximate the temperature near the surface of the water where the development process of the immature mosquitoes occurs (Angel & Joshi, 2008)).

The weather-driven model for the transmission dynamics of a vector-borne disease (dengue), with vertical transmission in the vector, is given by the following deterministic system of non-linear differential equations (the state variables and parameters of the model are described in Table 2, their ranges and values are given in Table 3; a flow diagram of the model is depicted in Fig. 2):where,

Table 2

Description of variables and parameters of the model (2.1).

Symbol	Description
Variables
S_E	Number of susceptible eggs
I_E	Number of infected eggs
S_L	Number of susceptible larvae
I_L	Number of infected larvae
S_P	Number of susceptible pupae
I_P	Number of infected pupae
S_M	Population of susceptible adult female mosquitoes
I_M	Population of infected adult female mosquitoes
S_H	Population of susceptible humans
E_H	Population of latently-exposed humans
I_H	Population of symptomatically-infected humans
R_H	Population of recovered humans
N_VI	Total population of immature mosquitoes
N_VA	Total population of adult female mosquitoes
N_H	Total population of humans
Parameters
σ_E(T)	Hatching rate of eggs into larvae
μ_E(T)	Natural mortality rate of eggs
σ_L(t)	Maturation rate of larvae to pupae
μ_L(T)	Natural mortality rate of larvae
δ_L(T)	Density-dependent mortality rate of larvae
σ_P(T)	Maturation rate of pupae to adult mosquito
μ_P(T)	Natural mortality rate of pupae
λ_HV(T,N_H,N_VA)	Transmission rate from infected humans to susceptible mosquitoes
λ_VH(T,N_H,N_VA)	Transmission rate from infected mosquitoes to susceptible humans
a_V(T)	Per capita contact (biting) rate of adult female mosquitoes on humans
ϕ_V(T)	Per capita egg oviposition rate
μ_V(T)	Natural mortality rate of adult mosquitoes
K_V(t)	Carrying capacity of breeding habitats for adult female mosquitoes to lay eggss
f_V	Proportion of new adult mosquitoes that are females
μ_H	Natural mortality rate of humans
σ_H	Rate of development of disease symptoms in humans
γ_H	Recovery rate for humans
Π_H	Recruitment rate (by birth and immigration) into the community
β_V	Probability of infection of a susceptible mosquito per bite on an infected human
β_H	Probability of infection of a susceptible human per bite by an infected mosquito
r	Proportion of infected eggs laid by infected adult female mosquitoes
	(due to vertical transmission)

Table 3

Values and ranges of the parameters of autonomous version of the model (2.1).

Parameter	Range	Baseline value	Reference
σ_E	(0.1,0.5)/day	0.4/day	(Clements, 1999; Feng & Jorge, 1997)
μ_E	(0.07,0.3)/day	0.2	(Clements, 1999; Feng & Jorge, 1997)
σ_L	(0.08,0.35)/day	0.14/day	(Clements, 1999; Feng & Jorge, 1997)
μ_L	(0.07,0.3)/day	0.18/day	(Clements, 1999; Feng & Jorge, 1997)
σ_P	(0.1,0.5)/day	0.3/day	(Clements, 1999; Feng & Jorge, 1997)
μ_P	(0.07,0.25)/day	0.17/day	(Clements, 1999; Feng & Jorge, 1997)
a_V	(0,1)/day	0.12/day	(Andraud, Hensand, Marais, & Beutels, 2012)
μ_V	(0.047,0.071)/day	0.05/day	(Agusto et al., 2015; Chitnis et al., 2006; Laperriere et al., 2011)
μ_H	(0.00003,0.000042)/day	0.00005/day	(United Nations and Department of Economic and Social Affairs and Population Division, 2011)
σ_H	(0,1)/day	0.15/day	(Garba et al., 2008)
γ_H	(0,1)/day	0.1428/day	(Garba et al., 2008)
Π_H	(60,300)/day	66/day	(City populationCHIANG MAI, 2017)
β_V	(0.3,0.75)	0.5	(Feng & Jorge, 1997; Garba et al., 2008)
β_H	(0.1,0.75)	0.4	(Feng & Jorge, 1997; Garba et al., 2008)
r	(0,0.3)	0.007	(Bosio, Thomas, Grimstad, & Rai, 1992; Freier & Rosen, 1987; Shroyer, 1990)
f_V	(0.4,0.6)	0.55	(Lounibos & Escher, 2008)
K_V	(10^4,10^6)	40,000	(Agusto et al., 2015; Laperriere et al., 2011)
ϕ_V	(1500)	1.84	(Agusto et al., 2015; Chitnis et al., 2006; Laperriere et al., 2011)

Fig. 2

Flow diagram of the model (2.1). Notation: Red arrow indicates infection route.

In the model (2.1), eggs are laid by adult female Aedes aegypti mosquitoes at a logistic growth rate , where is the temperature-dependent egg oviposition rate, is the carrying capacity of the breeding habitats for adult female mosquitoes to lay eggs. The notation , where with , is used to ensure that the term for all t. The parameter r (with ) represents the proportion of mosquito offsprings that are born infected (due to vertical transmission). Table 1 shows the average proportion of eggs laid that are infected, for various subtypes of dengue fever, based on laboratory experiments (Grunnill & Boots, 2016; Rosen et al., 1983).

Table 1

Proportion of infected eggs (r) for various dengue subtypes.

Species	Dengue subtype	Proportion of infected eggs (r) per1,000	Reference
Ae. aegypti	DENV-1-4	[0.61–15]	(Grunnill & Boots, 2016; Rosen et al., 1983)
Ae. albopictus	DENV-1	[1.7–14]	(Grunnill & Boots, 2016)
	DENV-2	[0.91–2.5]	(Grunnill & Boots, 2016)
	DENV-3	[0.23–0.78]	(Grunnill & Boots, 2016)
	DENV-4	[0.22–5.2]	(Grunnill & Boots, 2016)

Eggs hatch into larvae at a temperature-dependent rate , larvae mature into pupae at a temperature-dependent rate , and pupae become adult female mosquitoes at a temperature-dependent rate (where is the proportion of new adult mosquitoes that are females). Eggs, larvae and pupae suffer natural mortality loss at temperature-dependent rates and , respectively. Larvae suffer additional density-dependent mortality at a temperature-dependent rate (Lutambi, Penny, Smith, & Chitnis, 2013). A schematic diagram of the mosquito lifecycle is depicted in Fig. 1.

Fig. 1

Schematic description of the lifecycle of the Aedes aegypti mosquito (Entomology Department at Purdue University, 2017).

Susceptible adult female mosquitoes acquire dengue infection, following effective contact with an infected human (after taking a blood meal), at a rate (given in (2.2)), where is the temperature-dependent effective contact (biting) rate of adult female mosquitoes on humans (regardless of infectious status of the vector or the host), is the temperature-dependent probability that a bite from a susceptible adult female mosquito to an infected human results in an infection. Adult mosquitoes suffer natural death at a temperature-dependent rate (). The parameter represents the per capita recruitment rate of humans into the community (by birth or immigration). Susceptible humans acquire dengue infection, following effective contact with an infected adult female mosquito, at a rate (given in (2.2)), where is as defined above, and is the temperature-dependent probability of infection from an infected mosquito to a susceptible human per bite. Exposed humans develop clinical symptoms of the disease at a rate . Infectious humans recover at a rate (it is assumed that recovery induces permanent immunity against reinfection). Natural death occurs in all human compartments at a rate . Since dengue-induced mortality in humans is generally negligible (Agusto, Gumel, & Parham, 2015; Chitnis, Cushing, & Hyman, 2006; Garba et al., 2008; Laperriere, Brugger, & Rubel, 2011; Samui Times, 2017) (for instance, there were only 126 dengue-induced fatalities in Thailand in 2017 (Samui Times, 2017)), no human disease-induced mortality is assumed in the model. Furthermore, no disease-induced mortality is assumed in the adult mosquito population.

The model (2.1) is an extension of numerous dengue transmission models that include vertical transmission in the vector population (such as those in (Adams & Boots, 2010; Coutinho et al., 2006; Esteva & Vargas, 2000; Garba et al., 2008; Grunnill & Boots, 2016; Halide & Ridd, 2008)) by (inter alia):

adding the dynamics of immature mosquitoes (this was not included in (Esteva & Vargas, 2000; Garba et al., 2008; Halide & Ridd, 2008));

incorporating the effect of temperature variability (this was not considered in (Adams & Boots, 2010; Esteva & Vargas, 2000; Garba et al., 2008));

including the effects of temperature on vertical transmission in the vector population and on the dynamics of dengue disease (this was not considered in (Coutinho et al., 2006; Esteva & Vargas, 2000; Halide & Ridd, 2008)).

Functional forms of temperature-dependent parameters for adult mosquitoes

The functional forms of the temperature-dependent parameters of the model related to adult mosquitoes Aedes aegypti (namely, , , , and ) are defined as follows (for C C):

The biting rate of adult female mosquitoes on the human host () is given by (see (Scott, Amerasinghe, & Morrison, 2000), Fig. 5):

Fig. 5

Profile of the reproduction number () as a function of mean monthly temperature for Chiang Mai Thailand. Parameters values used as given by the baseline values in Table 3.

Table 5

Mean monthly temperature (in °C) for Chiang Mai province of Thailand for the period 2005–2016 (Thai Meteorological Department, 2017).

Month	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec
High	28.9	32.2	34.9	36.1	34.1	32.3	31.7	31.1	31.3	31.1	29.8	28.3
Mean	20.5	22.9	26.4	28.7	28.1	27.3	27.0	26.6	26.5	25.8	23.8	21.0
Low	13.7	14.9	18.2	21.8	23.4	23.7	23.6	23.4	23.0	21.8	19.0	15.0

Table 6

Average monthly dengue incidence data for Chiang Mai province of Thailand, for the period 2005–2016 (Centers for Disease Contorl and Prevention, 2016).

Month	Dengue cases per100,000
January	3.8
February	2.3
March	3.5
April	7.5
May	21.4
June	51.63
July	73.33
August	72.61
September	45.14
October	23.82
November	17.44
December	7.38

Table 7

Fitted values of monthly biting rates (obtained from fitting the model with data).

Month	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec
Biting rate (a_V)	0.05	0.01	0.01	0.09	0.15	0.12	0.15	0.1	0.05	0.03	0.02	0.02

The probability of infection from an infected mosquito to a susceptible human () per bite is given by (see (Lambrechts, Paaijmansb, & Fansiri, 2011), Supporting Information):

The probability of infection from an infected human to a susceptible mosquito per bite () is given by (see (Lambrechts et al., 2011), Supporting Information):

The oviposition rate () is given by ((Lambrechts et al., 2011)):

The mortality rate of the Aedes aegypti mosquito is given by (Yang, Macoris, & Galvani, 2009):

The aforementioned functional forms are depicted in Fig. 3. Furthermore, typically, a sinusoidal function of the following form is used to account for hourly fluctuations in local ambient temperature (Okuneye, Eikenberry, & Gumel, 2018):where is the mean daily air temperature, captures variation about the mean (i.e., is the diurnal temperature range), and denotes for time in hour for any given day.

Fig. 3

Profile of the functional forms of the temperature-dependent parameters of the model (2.1) related to the adult Aedes aegypti mosquito. (a) Probability of infection from an infected mosquito to a susceptible human per bite (). (b) The biting rate of adult female mosquitoes (). (c) Probability of infection from an infected human to a susceptible mosquito per bite (). (d) Egg oviposition rate (). (e) Mortality rate of adult female mosquito ().

If a formulation such as (2.8) is used for the temperature-dependent functional forms, then the parameters defined in Equations (2.5), (2.6), (2.7), (2.3), (2.4) are time-dependent. Hence, the model (2.1) is non-autonomous. However, if fixed temperature values are used (e.g., using the mean daily or mean monthly temperature), then each of the parameters defined in Equations (2.5), (2.6), (2.7), (2.3), (2.4) is constant. Hence, the model (2.1) is autonomous in this case. In the absence of good data to realistically derive the functional forms for the other temperature-dependent parameters related to the immature Aedes aegypti mosquitoes (i.e., , , , , and ), numerical simulations of the model (2.1) will be carried out using fixed (constant) values for these parameters (available in the literature).

Data fitting

The model (2.1) is, first of all, fitted using available epidemiological and weather data relevant (see Table 3, Table 5, Table 6) to dengue transmission dynamics in the Chiang Mai province of Thailand (Bureau of Epidemiology, 2015; Thai Meteorological Department, 2017) (using least square regression). In particular, both the temperature data (provided by Thai Meteorological Department (Thai Meteorological Department, 2017)) and incidence data (provided by Thailand Bureau of Epidemiology (Bureau of Epidemiology, 2015); see also Appendix A) are given for monthly periods between 2005 and 2016. The mean monthly temperature for Chiang Mai (given in Table 5) is plotted alongside the average monthly dengue incidence (given in Table 6) in Fig. 4 (a). This figure shows that peak dengue incidence is attained for temperatures between C and C, which are recorded in the Chiang Mai province during the period between June and August annually.

Fig. 4

Data fittings of the autonomous version of the model (2.1). (a) Plot of the average monthly temperature (in °C) for Chiang Mai (Table 5) superimposed with the monthly dengue incidence data (Table 6). (b) Data fitting of the model (2.1) using the average monthly incidence data in Table 6 (and the fitted monthly biting rates in Table 7). (c) Fitted biting rate () used to fit the model with the data. Plots are generated using the baseline parameter values in Table 3.

Furthermore, the model (2.1) is fitted using the aforementioned mean monthly temperature and incidence data (Table 5, Table 6) using the baseline parameter values given in Table 3, where the temperature-dependent biting rate () is chosen as a fitting parameter (and the remaining temperature-dependent parameters of the model are computed using their respective functional forms given in Section 2.1). The result obtained, depicted in Fig. 4 (b), shows a reasonably good fit. Plots of the fitted biting rates (using the data in Table 3) and the incidence data (given in Table 6) are depicted in Fig. 4 (c), from which it follows that the dengue incidence positively correlates with increasing biting rate.

Basic qualitative properties

Since the model (2.1) monitors the temporal dynamics of mosquitoes (immature and mature) and humans, all parameters of the model are assumed to be non-negative. Define the regionis positively-bounded and invariant for the model (2.1), where .

If the initial values of the system (2.1) lie in the region , then there exists a unique positive solution for (2.1), such that

Proof. The right-hand side of the expressions in the model (2.1) are continuous, with continuous partial derivatives in . Hence, the model has a unique solution that exists for all time t (Perko, 1991). It remains to be shown that the region is forward (positively)-invariant. Setting in the model (2.1) shows that for all . Similarly, it can be shown that , , , , , , , , , and for all (Theorem A4 of (Thieme, 2003)).

For the invariance of , it is convenient to let (for each initial data) . Adding the first eight equations of the model (2.1) gives:where, . Since , for all (and noting that is bounded for all (Abdelrazec & Gumel, 2016)), then:

Finally, adding the last four equations of the model (2.1) gives;so that,

Thus, and are positively-bounded for all . ∎

Analysis of autonomous version of the model

The model (2.1) will, first of all, be analysed for the case when fixed temperature values are used (i.e., the autonomous equivalent/version of the model (2.1) will be considered first).

Disease-free equilibria

It is convenient to define the quantity

The autonomous version of the model (2.1) has two disease-free equilibria, described below:where,with , and in . It follows from (3.2) that exists if and only if . The quantity represents the average number of new adult female mosquitoes generated by a single susceptible adult female mosquito that has successfully taken a blood meal. It is the product of the rate at which eggs are laid eggs by an adult female mosquito (), the proportion of eggs survived to become larvae (), proportion of larvae survived and matured into pupae (), proportion of pupae survived and become adult female mosquitoes (), and the average lifespan of a susceptible adult female mosquito ().

Trivial (mosquito-free) disease-free equilibrium (TDFE):

Non-trivial (mosquito-present) disease-free equilibrium (NDFE):

Local asymptotic stability of TDFE ()

The TDFE , which always exists, corresponds to the case without mosquitoes. It can be shown, by linearizing the autonomous version of the model (2.1) around , that the associated eigenvalues of the linearization have negative real part wheneverwhere, is given in (3.1). Since we assumed , it follows that in this case. The following result can be established using standard linearization of the autonomous version of the model (2.1) around the TDFE ().

The TDFE point () is locally-asymptotically stable (LAS) whenever , and unstable if .

It is worth noting that, since the quantity is the average number of new susceptible adult female mosquitoes generated by a single susceptible adult female mosquito that has successfully taken a blood meal, the quantity measures the average number of new infected adult female mosquitoes generated by a single infected adult female mosquito that has successfully taken a blood meal. It is worth noting, from the expression (3.3), the threshold quantity increases with increasing values of r. That is, as expected, increasing the vertical transmission rate () increase the disease burden (by increasing in ).

Local asymptotic stability of NDFE

The local asymptotic stability of can be established using the next generation operator method (Diekmann, Heesterbeek, & Metz, 1990; van den Driessche & Watmough, 2002). The non-negative matrix of new infection terms and the matrix of the transition terms associated with the autonomous case of the model (2.1) are, respectively, given by:where, , , , , , and . It follows, from (van den Driessche & Watmough, 2002), that the basic reproduction number () of autonomous case of the model (2.1) is given by (where ρ is the spectral radius):with, and . It is worth noting from (3.5) that, in the absence of vertical transmission, the reproduction threshold () reduces to

The result below follows from Theorem 2 of (van den Driessche & Watmough, 2002).

The NDFE () is LAS whenever , and unstable if .

The epidemiological implication of Theorem 3.2 is that a small influx of infected individuals or vectors into the population will not generate a large outbreak in the community if . Hence, the disease may be effectively-controlled if can be brought to (and maintained at) values less than unity. A plot of , as a function of fixed mean monthly temperature (for C) for the Chiang Mai province of Thailand (Thai Meteorological Department, 2017), is depicted in Fig. 5. This figure shows, that the profile of lies in the range . Furthermore, the values of increase with increasing temperature values in the range C (and decrease thereafter, for increasing temperatures above the peak temperature of C). Thus, disease burden increases with increasing temperature values in the range C, and decreases for increasing temperature values thereafter.

Effects of vertical transmission on the reproduction number

The effect of vertical transmission in the vector population (i.e., ) on the reproduction number () is assessed by using the approach in (Shuai, Heesterbeek, & den Driessche, 2013) as follows. Let (recall that the matrices and are given by Equation (3.4)) be the next-generation matrix of the autonomous version of the model (2.1) (Diekmann et al., 1990), given by:where, , , , , , , , , , . The notion of target reproduction number (as introduced by Shuai et al. (Diekmann et al., 1990; Shuai et al., 2013a)) will be used. Using the notation in (Shuai et al., 2013a), the entries () of the matrix K represent the average number of new cases of infections in humans (vectors) generated by an average infected (immature and mature) vector (humans). In particular, the entry is the effect of infected eggs () on the generation of new infected (exposed) humans (i.e., ). Furthermore, following (Shuai et al., 2013a), let be the set of target entries and and . Following (Shuai et al., 2013a), it is convenient to define the following projection matrices associated with the autonomous case of the model (2.1). Let I be the identity matrix, , and be matrices with entries if and otherwise, and and are projection matrices (e.g., if and otherwise) (Shuai et al., 2013a). It follows that the target reproduction number (denoted by ) with respect to the set S is given by (Shuai et al., 2013a):provided the spectral radius . Hence,provided . It is worth noting that, the expression in Equation (3.6) accounts for the average number of infected humans caused by an infected egg (after maturation to adulthood). Fig. 6 compares the target reproduction number () with the basic reproduction number () of the autonomous version of the model (2.1) for various values of r, from which it follows that is always less than for . Furthermore, for the range of the vertical transmission rate for Aedes aegypti mosquitoes given in Table 1 (where ), it follows from Fig. 6 that vertical transmission has very marginal population-level impact on the disease dynamics (vertical transmission becomes relevant far larger values of r, outside the aforementioned realistic range). This result is consistent with the numerical simulation results reported in (Adams & Boots, 2010).

Fig. 6

Plot of and as function of r. Parameters values used are as given by the baseline values in Table 3.

Fig. 7 shows that, in the absence of human-mosquito interaction (i.e., no mosquito bites on humans), vertical transmission (r) has very marginal effect on the abundance of infected adult female mosquitoes (Figure (7 a)). However, when the human-mosquito interaction is slightly increased (such as by setting the biting rate to ), the number of infected adult female mosquitoes significantly increases with increasing values of the proportion of new infected eggs (Figure (7 b)). Figure (7 b) clearly shows that vertical transmission in the vector population has little or no effect on the number of infected adult female mosquitoes (this result supports the finding in Fig. 6).

Fig. 7

Simulations of the model (2.1) showing the effect of the proportion of infected eggs (r) on the disease dynamics for (a) mosquito-human interaction set at , (b) mosquito-human interaction set at . Parameter values used are as given in Table 3.

Endemic equilibria

In this section, conditions for the existence of endemic equilibria will be derived. Let be any arbitrary endemic equilibrium of autonomous version of the model (2.1). Furthermore, let

It is convenient to definewhere, and are as defined in Section 3.1.2. It can be shown that () if and only if () (so that behaves like the target reproduction number discussed in (Shuai et al., 2013b)). It can be shown, by solving for the variables of the autonomous version of the model (2.1) at steady-state, that the solutions of the autonomous model satisfy the following linear equation in terms of :where, , . It follows from (3.9) that , (so that () if ()). The components of the positive equilibrium of autonomous version of the model (2.1) can then be obtained by solving for from (3.9), and substituting the result into the steady-state expressions for each of the state variables in (2.1). It follows that the autonomous case of the model (2.1) has a unique endemic equilibrium whenever . Thus, the autonomous version of the model (2.1) will not have an endemic equilibrium point if . Hence, the model will not undergo the usual phenomenon of backward bifurcation (Agusto et al., 2015; Chitnis et al., 2006; Garba et al., 2008; Laperriere et al., 2011). A global asymptotic stability result is given below for the NDFE (). It is convenient, first of all, to define the threshold quantity:where, and are as defined in Sections 3.1 and 3.1.2, respectively. We claim the following result.

The NDFE () of the autonomous version of the model (2.1) is GAS in whenever and .

The Proof of Theorem 3.3, based on using the approach in (Kamgang and Sallet, 2005, 2008), is given in Appendix B. The epidemiological significance of Theorem 3.3 is that, for the autonomous version of the model (2.1), reducing (and maintaining) the threshold quantity () to a value less than unity is necessary and sufficient for the effective control or elimination of the disease in the community.

Simulations: effect of temperature variability

The effect of temperature variability on the disease dynamics, as a function of vertical transmission in the vector (r), is monitored by simulating the model (2.1) with various fixed temperature values in the range C to C. In the absence of vertical transmission (i.e., ), Fig. 8 (a) shows that the disease burden, as measured in terms of the number of new infected humans, increases with increasing temperature values until C (where the maximum peak is attained). The disease burden then decreases for increasing temperatures thereafter. Similar pattern is observed when vertical transmission is assumed (Fig. 8 (b)), although an increase in the number of new cases in humans is observed as expected. Furthermore, the number of infected mosquitoes exhibit similar pattern (that is, the number of infected adult female mosquitoes increases with increasing temperature until C, and decreases for increasing temperatures thereafter), although oscillatory dynamics, due to the assumed logistic eggs oviposition rate (), was observed (Fig. 8 (c) and (d)).

Fig. 8

Simulations of the model (2.1), for the effect of temperature and vertical transmission on disease dynamics. (a) Total number of symptomatic humans () as a function of time for . (b) Total number of symptomatic humans () as a function of time for . (c) Total number of infected adult mosquitoes () as a function of time for . (d) Total number of infected adult mosquitoes () as a function of time for . Parameters values used are given in Table 3. The functional forms for the temperature-dependent parameters (, , , and ), given in Section 2, are used.

The overall effect of temperature on disease dynamics is assessed by simulating the model (2.1) using various values of temperature in the range C (for Chiang Mai province of Thailand (Thai Meteorological Department, 2017)). The results obtained (depicted in Fig. 9) show that the total number of new dengue cases (in both the human and mosquito populations) reaches a peak during the period June to August (which correspond to the temperature range of C– C).

Fig. 9

Simulation of the model 2.1 for the disease dynamics in Chiang Mai, Thailand. (a) Monthly total number of infected mosquitoes. (b) Monthly total number of infected humans. Parameters values used are given in Table 3. The functional forms for the temperature-dependent parameters (, , , and ), given in Section 2, are used, using mean monthly temperature (T) given in Table 5.

The potential impact of temperature variability on the burden of disease caused by vertical transmission (r) is monitored by simulating the non-autonomous model (2.1) using various value of r and temperature. Fig. 10 shows that the effect of r in generating new infected cases is more pronounced for temperature values in the range C (Fig. 10 (a) and (b)) and decreases thereafter (Fig. 10 (d) and (e)). Thus, this study shows that the ability of vertical transmission to cause significant increase in disease burden is temperature-dependent.

Fig. 10

Cumulative number of infected humans for various values of the proportion of infected eggs laid by an infected mosquito () and temperature (T): (a) , (b) , (c) , (d) , (e) . Color notation from blue () to gold () represent varying values of r, from 0 to 0.9, in steps of length 0.1.

Uncertainty and sensitivity analysis

The autonomous version of the model (2.1), with fixed temperature values, contains 19 parameters, and uncertainty in their estimates are expected to arise. The effect of such uncertainties is assessed using uncertainty and sensitivity analysis (Cariboni, Gatelli, Liska, & Saltelli, 2007). In particular, Latin Hypercube Sampling (LHS) and Partial Rank Correlation Coefficients (PRCC) is used for this model, as below. The purpose of sensitivity analysis is to assess the effects of parameters on the outcomes of the simulations of the model (Cariboni et al., 2007). A highly-sensitive parameter should be more carefully estimated, since a small change in that parameter can cause a large quantitative change in the result (Cariboni et al., 2007). On the other hand, a parameter that is not sensitive does not require as much attempt to estimate (because a small change in that parameter will not cause a large variation to the quantity of interest) (Blower & Dowlatabadi, 1994; Cariboni et al., 2007; Marino, Hogue, Ray, & Kirschner, 2008). The analyses will be carried out using data (temperature, demographic and epidemiological) relevant to dengue transmission dynamics in the Chiang Mai province of Thailand (City populationCHIANG MAI, 2017; Lambrechts et al., 2011; Scott et al., 2000; Yang et al., 2009). The total population of the Chiang Mai province is estimated to be 1.7 million, and the average lifespan is 65–72 years (City populationCHIANG MAI, 2017) (so that per day with a mean of ). Thus, million. Hence, per day. The analysis will be carried out using the baseline values and ranges tabulated in Table 3.

Using the population of infectious humans () as the response function, it is shown in Table 8 that the top PRCC-ranked parameters of the model (i.e., parameters with PRCC values greater or equal to 0.5) are the mosquito biting rate () and the transmission probability per contact for susceptible human (). Similarly, using the population of infectious mosquitoes () as the response function, the top PRCC-ranked parameters are the egg deposition rate (), maturation rate of pupae to adult mosquito ) and the environmental carrying capacity of immature mosquitoes ). Furthermore, using the population of infected pupae () as the response function, the top PRCC-ranked parameters of the model are the egg deposition rate (), the environmental carrying capacity () and the maturation rate of larvae to pupae (). Considering the population of infected larvae () as the response function, the top PRCC-ranked parameters of the model are egg deposition rate () and maturation rate of eggs to larvae (). Finally, using the population of infected eggs () as the response function, the top PRCC-ranked parameters of the model are egg deposition rate (), the environmental carrying capacity () and the natural death rate of adult female mosquitoes ().

Table 8

PRCC values for the parameters of autonomous case of the model (2.1) using the total number of infected eggs (), larvae (), pupae (), adult mosquitoes () and infectious humans () as response functions (with PRCC ). The top parameters that affect the model with respect to each of the six response functions are highlighted in bold font.

Parameters	I_E	I_L	I_P	I_M	I_H
σ_E	−0.21966	+0.8402	+0.5515	0.6804	−0.1074
μ_E	−0.0241	+0.1490	−0.1112	+0.0171	−0.0313
σ_L	−0.0334	−0.2195	+0.5821	+0.4081	+0.0458
μ_L	+0.2574	−0.3229	−0.0120	−0.0845	−0.0513
σ_P	+0.3572	+0.0673	−0.2042	+0.7500	−0.0245
μ_P	+0.2973	+1453	−0.3415	+0.0376	+0.0100
a_V	−0.1988	+0.3124	−0.1143	+0.281	+0.9463
μ_V	−0.4202	−0.0292	0.0854	+0.1701	−0.1462
μ_H	−0.3149	−0.3936	+0.2996	+0.2687	0.0856
σ_H	+0.0409	−0.1667	+0.0184	+0.2364	−0.0995
γ_H	+0.0475	+0.0176	−0.0262	+0.1769	−0.4349
Π_H	0.0506	+0.1099	+0.2162	−0.0028	−0.2039
β_V	−0.1364	−0.1257	+0.0843	+0.4012	+0.1832
β_H	−0.1176	−0.2754	−0.2405	−0.0476	+0.8009
r	−0.2010	+0.1766	−0.3528	+0.0543	−0.0186
f_V	+0.2057	+0.1517	+0.1094	+0.3173	−0.0903
K_V	+0.5157	+0.4100	+0.6001	+0.7446	+0.0310
ϕ_V	+0.8603	+0.8999	+0.7805	+0.8113	−0.0185

In summary, this study identifies four parameters that dominate the transmission dynamics of the autonomous version of the model (2.1), namely the environmental carrying capacity of immature mosquitoes (), biting rate (), probability of infection of a susceptible human () and the egg oviposition rate (). It is worth noting from Table 8 that the parameter related to vertical transmission in the vector (r) does not have significant PRCC values, suggesting that vertical transmission plays a marginal if at all role on the transmission dynamics of the disease.

Analysis of non-autonomous model

Consider, now, the full non-autonomous model (2.1), where the function (2.8), for the daily temperature fluctuations is used (so that the temperature-dependent parameters, given in Equations (2.5), (2.6), (2.7), (2.3), (2.4), are now functions of t). Although the concept of basic reproduction number has been extensively addressed for autonomous models for disease transmission over the decades, such a concept has only been recently extended to disease transmission models with periodic coefficients (see, for instance (Bacaër, 2007; Bacaër & Abdurahman, 2008; Bacaër & Ouifki, 2007; Wang & Zhao, 2008),). In this section, the methodology in (Wang & Zhao, 2008) will be used to compute the reproduction number associated with the non-autonomous model (2.1) with (2.8). Although the non-autonomous model (2.1) has two disease-free solutions, namely the trivial disease-free equilibrium and a non-trivial disease-free periodic solution, only the non-trivial disease-free periodic solution will be analysed (since the former, associated with the absence of mosquitoes in the population, is ecologically unrealistic). It is convenient to define the functional threshold quantity

The non-trivial disease-free solution (NDFS), obtained by setting in (2.1), has the form (recalling that ) with being the unique periodic solution (for for all ) satisfying:

The next generation matrix (of the new infection terms) and the M -Matrix (of the remaining transfer terms), associated with the non-autonomous model (2.1) with (2.8), are given, respectively, byand,where, , , , , , and , . Let be the vector field in the right-hand sides of Equation (4.1) and . Following (Wang & Zhao, 2008), let be the monodromy matrix of the linear ω-periodic systemwhere, . Further, let , be the evolution operator of the linear ω-periodic systemthat is, for each , the matrix satisfieswhere, is the identity matrix. Suppose that (ω - periodic in s) is the initial distribution of infectious individuals. Thus, is the rate at which new infections are produced by infected individuals who were introduced into the population at time s (Wang & Zhao, 2008). Since , it follows that is the distribution of those infected individuals who were newly-infected at time s, and remain infected at time t. Hence, the cumulative distribution of new infections at time t, produced by all infected individuals () introduced at a prior time , is given by (Wang & Zhao, 2008)

Let be the ordered Banach space of all ω-periodic functions from to , which is equipped with maximum norm and positive cone . Define a linear operator given by (Wang & Zhao, 2008)

The basic reproduction ratio for non-autonomous model (2.1) (denote by ) is then given by the spectral radius of the linear operator L, denoted by (Wang & Zhao, 2008). That is, . It can be verified that the assumptions in (Wang & Zhao, 2008) are valid for the model (2.1) with periodic parameters. Therefore, the result below is established from Theorem 2.2 in (Wang & Zhao, 2008).

Let for all . The NDFS (), of the non-autonomous model (2.1), is LAS if , and unstable if

We claim the following result.

Let . The NDFS () of the special case of the non-autonomous model (2.1) is GAS in if .

The Proof of Theorem 4.2 is given in Appendix C. Theorem 4.2 shows that the disease can be effectively-controlled or eliminated if the threshold quantity can be brought to and maintained at, a value less than unity. In other words, the prospects of such effective control in the Chiang Mai province is promising if the control measures implemented in the province can bring, and maintain, to a value less than unity.

Discussion and conclusions

This study is based on the design and analysis of a new deterministic model for assessing the impact of vertical transmission, in the vector population, and temperature variability on the transmission dynamics and control of dengue disease. Although dengue is primarily transmitted horizontally (via the vector-host-vector transmission cycle), vertical transmission has also been observed in the two main dengue-competent Aedes mosquitoes (namely Aedes albopictus and Aedes aegypti (Grunnill & Boots, 2016; Rosen et al., 1983)) and the human host population (Cosner et al., 2009). The consequence of the vertical transmission process is that infected vectors continue to emerge (during favorable temperature and habitat conditions) even when there are no infected hosts, since the infected eggs can survive the dry season and re-emerge as infected adult mosquitoes (Adams & Boots, 2010; Pacheco et al., 2009). Temperature affects the dynamics of both immature and adult stages of the dengue-competent mosquito lifecycle by generally affecting vector dynamics (e.g., survival, development, etc.) and mosquito-host interactions (e.g., biting) (Alto & Bettinardi, 2013).

The new model designed, which incorporates the dynamics of the aquatic stages of the mosquito (including logistic eggs oviposition and density-dependent larval mortality), vertical transmission effects in the vector, was rigorously analysed to gain insight into its dynamical features. It was further shown that, for small enough dengue mortality rate, the non-trivial disease-free equilibrium of autonomous version of the model is globally-asymptotically stable if the associated reproduction number of the model is less than unity. The epidemiological implication of this result is that the disease can be effectively-controlled if the control strategies implemented in the community can bring (and maintain) the reproduction number to a value less than unity. In other words, this result shows that bringing (and maintaining) the reproduction number to a value less than unity is necessary and sufficient for the effective control of the disease in the community.

The model is used to assess the population-level impact of vertical transmission. It should be recalled from Table 1 that proportion of dengue-competent vector born infected (r) is quite small (with ). Numerical simulations of the autonomous version of the model (Fig. 6(a)) show that vertical transmission has very marginal effect on the disease dynamics. However, when temperature effects are incorporated into the model, simulations of the resulting model show that the effect of vertical transmission is more pronounced for temperature values in the range C (Fig. 10). This effect decreases for temperature values greater than C.

The model designed in this study contains numerous parameters, and the effect of the associated uncertainties of the parameters on the numerical simulations of the model was assessed using Latin Hypercube Sampling (LHS) and Partial Rank Correlation Coefficients (PRCC) (Blower & Dowlatabadi, 1994; Cariboni et al., 2007; Marino et al., 2008), based on parameter values and ranges relevant to dengue transmission dynamics in the Chiang Mai province of Thailand (Thai Meteorological Department, 2017). These analyses reveal some of the parameters of the model that play a dominant role on the disease transmission dynamics, including the mosquito carrying capacity, biting rate and eggs oviposition rate. Hence, effective dengue control is dependent on the design of strategies that reduce the values of these parameters (e.g., using larviciding and adulticiding to reduce the egg oviposition rate and mosquito carrying capacity, using insecticide-treated bednets and insect repellents to minimize the biting rate). Furthermore, simulations of the model show that vertical transmission has very marginal (if at all) impact on the disease transmission dynamics in the community. Finally, it is shown that dengue-associated burden, as measured in terms of the total number of new dengue cases in humans, increases with increasing mean monthly temperature in the recorded range for Chiang Mai ( ). Further, such burden is maximized when the mean monthly temperature lie in the range . This range is recorded in Chiang Mai province during 3–4 months of the year (between June and August). Thus, this study suggests that anti-dengue control efforts should be intensified in the Chiang Mai province of Thailand during these months (this result supports the finding in (Abdelrazec, Bélair, Shan, & Zhu, 2008)). Furthermore, dengue-associated burden decreases with decreasing mean monthly temperature below C and above C and this result supports the finding in (City populationCHIANG MAI, 2017; Garba et al., 2008; Gumel, 2012; Lambrechts et al., 2011; Scott et al., 2000; Yang et al., 2009).

Table 9

Average monthly DENV incidence in Chiang Mai, Thailand, for the period of 2005–2016.

Month	2005	2006	2007	2008	2009	2010	2011	2012	2013	2014	2015	2016	Average (per 100,000)
January	4	8	1	21	87	44	29	19	138	12	17	80	3.8
February	12	4	3	21	44	30	7	23	90	9	13	28	2.3
March	9	8	3	27	34	44	11	5	175	2	7	42	3.5
April	21	18	12	74	54	45	12	20	573	5	25	61	7.5
May	164	95	32	227	151	158	87	65	1293	19	173	109	21.4
June	168	301	99	591	314	525	142	170	3120	89	400	277	51.63
July	148	217	160	987	322	1850	93	252	3146	150	400	1074	73.33
August	103	137	156	971	287	2304	96	335	1691	184	826	1624	72.61
September	79	47	98	561	177	1153	49	332	818	168	1067	868	45.14
October	56	35	48	383	144	283	22	373	240	83	889	303	23.82
November	35	16	43	270	133	69	40	215	106	40	911	215	17.44
December	9	9	10	128	44	41	13	129	42	25	363	73	7.38

Table 10

Full monthly DENV incidence in Chiang Mai, Thailand, for the period of 2005–2016 (Bureau of Epidemiology, 2015).

Month 2005–2016	Temp ( ∘C)	Rain (mm)	DENgue cases
2005
1	22.6	0	4
2	25.5	0	12
3	27.2	24.7	9
4	29.8	57.2	21
5	29.6	104.7	164
6	29	193.5	168
7	28.4	179.1	148
8	27	155.2	103
9	26.9	436.3	79
10	26.8	192	56
11	25.5	22.8	35
12	22.6	27.9	9
2006
1	22.4	0	8
2	25.1	0	4
3	28	18	8
4	29.2	206.7	18
5	26.4	219.5	95
6	28.6	180.4	301
7	26.3	269.3	217
8	26	341.4	137
9	26.6	194.8	47
10	25.7	69.9	35
11	23.9	0	16
12	21.8	0	9
2007
1	21	0	1
2	23.3	0	3
3	26.5	0	3
4	29.5	56	12
5	26.4	393.5	32
6	27.8	130.1	99
7	26.9	74.6	160
8	26.9	153.2	156
9	26.8	179.8	98
10	25.7	64.6	48
11	23	73.5	43
12	21.8	0	10
2008
1	22.2	16.6	21
2	24.5	13.8	21
3	27.7	9.4	27
4	29.8	57.2	74
5	27.3	158.7	227
6	27.9	147.1	591
7	27.7	101.6	987
8	27.2	170.9	971
9	26.9	236.4	561
10	26.6	188.1	383
11	24.2	34.1	270
12	21.5	7.1	128
2009
1	21.3	0	87
2	25.3	0	44
3	27	16.7	34
4	29.5	97.9	54
5	28.4	142	151
6	27.6	140.2	314
7	27.7	124	322
8	27.9	126.8	287
9	27.9	191.7	177
10	27.3	223.4	144
11	25	0	133
12	22.4	7.5	44
2010
1	24.1	21.7	44
2	24.4	0	30
3	26.9	0	44
4	31.5	3.9	45
5	31.3	46.4	158
6	29.6	122.7	525
7	28.5	114.5	1850
8	27	470.6	2304
9	27.4	196.2	1153
10	26.8	169.6	283
11	24.7	0	69
12	23.4	6.1	41
2011
1	22.6	2.6	29
2	24.4	0.8	7
3	25.3	60.4	11
4	27.2	92.6	12
5	27.2	292.7	87
6	27.7	216.8	142
7	27.6	191.2	93
8	26.8	260.5	96
9	27.1	254.9	49
10	26.4	69.7	22
11	24.9	6.7	40
12	23	0.6	13
2012
1	23	11	19
2	25.1	0	23
3	27.4	8.3	5
4	29.3	75.9	20
5	28.5	216.4	65
6	28.1	55.9	170
7	27.5	106	252
8	27.6	185.4	335
9	27.6	179.6	332
10	27.3	80.1	373
11	26.9	38.8	215
12	24.2	1	129
2013
1	23.2	25	138
2	26.9	31.6	90
3	27.6	17.1	175
4	31.2	1.2	573
5	29.8	89.9	1293
6	28.9	39.7	3120
7	27.9	272.9	3146
8	27.3	299.4	1691
9	27.4	275.6	818
10	26.2	123.4	240
11	26.4	85.4	106
12	21	26.8	42
2014
1	21.3	0	12
2	24.3	0	9
3	27.7	5.9	2
4	29.6	34.9	5
5	29.1	236.1	19
6	28.8	58.2	89
7	28	175.2	150
8	27.4	231.3	184
9	27.6	177.5	168
10	27.3	129.3	83
11	25.8	16	40
12	23.5	0	25
2015
1	22.3	78.9	17
2	24.3	0	13
3	27.9	27.5	7
4	29.7	53.8	25
5	30.4	76.5	173
6	29.9	15.2	400
7	28.2	120.2	613
8	28.2	143	826
9	28.2	139.4	1067
10	27.3	93.2	889
11	26.8	79.2	911
12	24.5	4.9	363
2016
1	21.6	34.2	80
2	24.2	45.3	28
3	29.4	0	42
4	32.4	17.7	61
5	31.1	85.7	109
6	28.1	236.1	277
7	27.6	162.1	1074
8	27.7	132.1	1624
9	27.6	213.1	868
10	27.6	141.7	303
11	26.3	105.3	215
12	24	6	73