Ento-Epidemiological Framework

We devised an ordinary differential equation (ODE) model to capture the transmission dynamics of dengue between human and mosquito hosts. The human population is assumed to have constant size () and to be fully susceptible to the virus. Upon challenge with infectious mosquito bites (), individuals enter the incubation phase () with mean duration of days, later becoming infectious () for days and finally recovering () with life-long immunity. The dynamics of the human population are defined by the following set of ODEs:(1)(2)(3)(4)(5)

For the dynamics of the vector population we consider the model previously formulated by Yang and colleagues [16], in which individuals are divided into two pertinent life-stages: aquatic (eggs, larvae and pupae, ) and adult females (). We further extend the adult class by subdividing into the epidemiologically relevant stages for dengue transmission: susceptible (), incubating () for days and infectious (). For ease of reading, the temperature-dependent entomological factors are herein distinguished by a (dot) notation (further details in the following sections). The system of equations describing the vector population is:(6)(7)(8)(9)(10)

Here, the coefficients and are the fraction of eggs hatching to larvae and the fraction of female mosquitoes hatched from all eggs, respectively. For simplicity and lack of quantifications for the local mosquito population, we assume these to be 1 (see the original publication for a discussion [16]). Moreover, denotes the rate of transition from aquatic to adults, and are the mortality rates, and is the intrinsic oviposition rate. The logistic term can be understood as the physical/ecological available capacity to receive eggs, scaled by the carrying capacity term , used in the fitting approach to indirectly estimate the adult mosquito population size (see below). From the above system, the basic offspring number (), that is, the mean number of viable female offspring produced by one female adult during its entire time of survival (and in the absence of any density-dependent regulation), can be derived as:(11)

All parameters defining are temperature-dependent (see below). For a fixed temperature it is possible to derive expressions for the expected population sizes of each mosquito life-stage modelled. These are used to initialize the system, given the temperature present at the initial timepoint:(12)(13)

The vector-to-human () and human-to-vector () incidence rates are assumed to be density-dependent and frequency-dependent (respectively), in respect to the type of infected host being considered:(14)(15)

Here, is the biting rate and and are the vector-to-human and human-to-vector transmission probabilities per bite. This approach follows the recent framework from Althouse et al. which conforms to the expectations arising from the constant nature of the number of bites per mosquito [17]: conceptually, (i) an increase in the density of infectious vectors should directly raise the risk of infection to a single human; while (ii) an increase in the frequency of infected humans raises the risk of infection to a mosquito biting at a fixed rate.

With the two hosts, the expression for dengue's basic reproductive number is defined equally to previous modelling approaches [18], [19] but without human mortality:

Temperature-Dependent Parameters

In this section we summarize the methodologies used for each of the seven entomological parameters dependent on temperature (Table 1). Here, is temperature in Celsius, is temperature in Kelvin and is the universal gas constant in .

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Temperature-dependent parameters.

https://doi.org/10.1371/journal.pntd.0003083.t001

In the study by Yang et al., from where we base the developmental part of our vector dynamical system (see above), temperature-controlled experiments were performed on populations of Aedes aegypti to derive closed-form expressions (based on polynomials) for the model's rates (see Figures 2, 3, 4 and 5 of the original publication [16]). We integrate such solutions into our framework:(16)(17)(18)(19)

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. Model fitting to Madeira's dengue outbreak data.

(A,B) Reported cases (incidence and cumulative) per week (dotted, black) and example of model fitting (solid, purple). Coloured area (purple) is the standard deviation of all accepted steps in the MCMC chain. The dashed vertical line represents the date of the first reported clinical cases. The red dashed line represents the epidemic progression ignoring the first week in November, when a new surveillance method was introduced. (C) Stationary distributions of the estimated timepoint of first case for 30 independent MCMC runs with random initial conditions and 1 million steps.

https://doi.org/10.1371/journal.pntd.0003083.g003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 4. Model-derived epidemiological and entomological parameter estimates for 2012.

(A) Example of estimated values for 2012 (solid, red) together with the weekly minimum temperatures for 2012 (solid, blue) and long-term average of minimum temperatures (2001–2011, dashed green). The dashed red line marks the epidemic threshold . (B) Example of estimated number of mosquitoes per human (solid, black), incubation period (solid, cyan) and adult life-span (solid, orange) for 2012.

https://doi.org/10.1371/journal.pntd.0003083.g004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 5. Model-derived epidemic potential for the island of Madeira.

(A) Temperatures for the year of 2012 (red, solid line) and average temperatures for the past 10 years (2001–2011; blue, solid line). The points mark the mean outbreak size (number of cases) for 100 stochastic introductions at different timepoints using temperature data from 2012 (red) and the average over the past 10 years (blue). (B) Derived real-time (red, solid line) for 2012, with an annual mean of (dashed line). (C) Derived real-time (blue, solid line) for the past 10 year, with an annual mean of (dashed line). (B,C) Grey shaded areas are the frequency of simulations (in 100) achieving either more than 3 (light grey) or 1000 (dark grey) cases.

https://doi.org/10.1371/journal.pntd.0003083.g005

The relationship between the extrinsic incubation period and temperature has been formulated by Focks et al. [20] using an enzyme kinetics model previously proposed by other authors [21] and used in other dengue modelling approaches [22]. The model assumes that the rate of development is determined by a single rate-controlling enzyme. The expression used is:(20)

The probabilities of transmission per mosquito bite and are modelled as previously estimated by Lambrechts and colleagues [23]. The data used in their study was both sampled from several other studies and obtained from de novo experiments that measured the variations in proportion of infected and transmitting vectors according to changes in temperature. The analysis was done for a variety of arboviruses from the flavivirus family, including the Dengue virus, the West Nile virus, Murray Valley Encephalitis virus and St. Louis Encephalitis virus. The expressions used are:(21)(22)

Constant Parameters

The framework described above has only three fixed parameters that are neither temperature-dependent nor estimated in the MCMC approach. These can be found in Table 2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Constant parameters.

https://doi.org/10.1371/journal.pntd.0003083.t002

Data Series

The outbreak time series was compiled from the official weekly reports from the Direcção Geral de Saúde (Portuguese health ministry) [10] issued throughout 2012 and 2013 and the special report by the European Centre for Disease Prevention and Control (ECDC) [9]. Temperature data for the island of Madeira was assembled from Weather Underground, a Weather Channel's repository [24]. For this we chose a weather station located in the centre of Funchal, Madeira's capital city, where most cases took place. We resorted to the website of Aeroportos da Madeira (Madeira Airports) for the statistics on airline traffic [25]. Finally, the figures for yearly investment in tourism were obtained from the official local source, the Instituto de Desenvolvimento Regional (Institute for Regional Development) [26].

Markov Chain Monte Carlo Fitting

For the fitting process a Markov chain Monte Carlo approach [27] is used to find combinations of parameters that can describe qualitative properties of Madeira's outbreak. We define the jumping distribution as being symmetric (Gaussian), effectively defining a random walk Metropolis-Hastings algorithm:(23)(24)(25)(26)(27)

Here, the Markov chain state is generally denoted by M, the proposal of new parameters by Y and the ODE system (described above) output by O. In step 1, is the Markov chain state of parameter at step , the pre-defined variance for each jump of parameter and the resulting proposal for time . In step 2, is the probability of acceptance. For this, we calculate the least squares distance between the data series and the ODE output for both the proposal of parameters and the previously accepted parameters . The probability is assumed to decrease exponentially with increases in least squares distances to the data.

With this simple approach we explored all possible combinations of values from four open parameters (Table 3) that are able to closely describe the outbreak time series. Amongst these is the carrying capacity K, which we explore in order to indirectly estimate the number of adult mosquitoes per human, and , the timepoint of the first case. We also consider two linear coefficients, and , that scale the mortality rate and incubation period of adult mosquitoes - we argue that these entomological factors, as defined by Yang et al. in laboratory experiments [16], should be adjusted to possible biological/ecological local effects. For example, it has been previously demonstrated that mosquito and virus genotype can have an effect on both susceptibility and incubation [28], [29], while human and predator behaviour, as well as the local geospatial topology can affect adult mortality [30], [31]. By considering these linear effects, we do not change the relative effect of temperature variation on mortality and incubation per se, but rather allow the baselines to be different from the ones obtained from the laboratory, ideal conditions of Yang et al. study. For a discussion on how much field and laboratory entomological factors can differ, see the recent work by Brady and colleagues [32].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Estimated parameters.

https://doi.org/10.1371/journal.pntd.0003083.t003

We address MCMC convergence by visual inspection and also quantify it using , the Gelman-Rubin statistic, which compares the variance between and within M independent MCMC chains [33]. Consider that each chain has length N steps and that is the parameter value in chain Then, when defining B as the between-chain and W as the within-chain variances, can be obtained using:(28)(29)(30)(31)(32)(33)(34)

is expected to approximate 1 when the M chains have converged to the same stationary distribution. Values significantly larger than 1, for instance, indicate that the between-chain variance is greater that the within-chain variance, highlighting that the MCMC may need more time to converge or tuning of jump parameters is required [33]. In our approach we calculate and present for each estimated parameter (Table 3) using 30 independent chains started with random initial conditions. Our jump parameters are chosen to assure that all MCMC chains presented in this study have acceptance rates and chains are run for at least 1 million steps.

Stochastic Formulation

A stochastic version of the ento-epidemiological framework was developed by introducing demographic stochasticity in the transitions of the dynamical system. We used multinomial distributions to sample the effective number of individuals transitioning between classes per time step. Multinomial distributions are generalized binomials, here defined as , where equals the number of individuals in each class and equals the probability of the transition event (equal to the deterministic transition rate). This approach has been demonstrated elsewhere, see e.g. [34].

Observations and Model Fitting

The rainfall data displays three distinct peaks during this period, a small one coinciding with the start of the outbreak, one during its peak and another two months after the reported cases dropped close to zero (Figure 2). We first note that the actual amount of rainfall over the whole 6-months period under consideration was very small, and it is reasonable to question how much of an impact this could have made on the local mosquito population during the epidemic, especially given the year-round availability of breeding sites [5]. There are also contrasting observations in relationship to rainfall timing and case-response: the first peak was short and small but followed by a drastic increase in case numbers; in contrast, the second much longer and heavier rain episode coincided with the epidemic peak but was followed by a sudden decrease in case numbers; finally, the third peak in rainfall took place outside the time-range of interest. We therefore argued that rainfall was unlikely to have been a main player in the overall progression of this particular short-lived epidemic, although a small contribution cannot be ruled out.

Compared to rainfall, which mostly affects Aedes's habitat-quality and availability, temperature directly affects both the mosquito life-cycle and viral replication rates within the mosquito [16], [20], [23], [32], [38]. Accordingly, we note that the drop in minimum temperature towards the end of the year correlated with the decrease in case numbers (Figure 2), possibly delayed by the integrated length of the intrinsic cycles of human-vector-human transmission and aquatic-to-adult mosquito development.

We initially tested a variety of SIR-based model frameworks with constant parameters in time, but it became clear that these models were unable to fit the sudden decline in case numbers after the epidemic peak in November. That is, models that would match the steep exponential rise in incidence, for example, and which would therefore predict very high reproductive numbers (), would inevitably generate epidemics of significantly higher magnitude and longer duration than the 2012 outbreak. As human intervention can be ruled out [9], and given the apparent strong correlation between declining temperatures and fading case numbers, we instead focused on a temperature-driven human-mosquito transmission model (see Methods) to gain more insight into the dynamic progression of this outbreak.

In order to fit our ento-epidemiological model to the case data under temperature variations, or more specifically to derive particular parameter combinations that are able to reproduce the timing, shape and size of the epidemic, we employed a Markov chain Monte Carlo (MCMC) approach (see Methods). Figures 3A and B show the fitted model together with the weekly and cumulative incidence data, respectively. There is an overall close fit between model output and the data, although we notice two apparent deviations: one at the onset of the epidemic and one just after its peak. We argue that these are more likely artefacts of the data, however, rather than real discrepancies. That is, the sudden drop in incidence after the epidemic peak is likely due to the introduction of a new surveillance system in the first week of November [9], whereas the slight overestimation of cases during the initial phase of the outbreak could be due to the deterministic nature of our model and possible under-reporting at the onset of the epidemic. We also note that the new surveillance system did not change the clinical or laboratory-confirmation definitions per se (see annexe of [9] for case definition) but aimed at efficiently integrating data from all health care centres involved, including the private and public sectors. It is therefore not expected that this change in the system affected the sensitivity of case detection but instead increased the time and space resolution of the epidemic data from November onwards.

Using this method we estimated that the possible timepoint of introduction of dengue to the island occurred towards the end of August, with the first autochthonous human infection between the and of September, two to three weeks before the first reported clinical case. Convergence to this date (range) was confirmed by independent MCMC runs as demonstrated in Figure 3C. That is, given random initial conditions for the four free parameters (Table 3), the system robustly converged towards equally distributed parameter estimates (by design, this approach produces parameter distributions rather than point estimates); see Methods and Supplementary Figures S1 and S2 for the resulting distributions of other parameters and quantification of convergence.

At the time of the first local human infection, temperatures were already on a declining trajectory (Figure 2), which caused a significant reduction in the virus's reproductive potential due to a combination of shorter mosquito life-expectancy, smaller population size and an increase in the extrinsic incubation period. These effects are demonstrated in Figure 4A, which also illustrates the colder than usual winter and slightly warmer summer during 2012 (compared to the average temperatures over the past 10 years). Dengue's estimated reproductive potential, , is here given as a time-dependent quantity to highlight the temperature-driven dependencies of entomological factors. It is important to note that this differs from the often used time-varying effective reproduction number, , that takes into account the varying susceptibility levels in the population (see Methods for mathematical expressions).

Of particular interest is the increase in the length of the extrinsic incubation period beyond the mosquitoes' average life-expectancies (Figure 4B), substantially contributing to the sharp drop in from at its peak at the end of the summer, to during late autumn and winter. We thus believe that this temperature-driven phenomenon might explain the rapid decrease in dengue incidence and essentially the end of the outbreak, with the expected delay due to the total length of the transmission and developmental cycles.

Epidemic and Endemic Potential

Using these parameter insights we next investigated potential outcomes if the pathogen would have been introduced at different timepoints during 2012 and further considered introductory events during a ‘typical’ year using average temperatures of the past 10 years (2001–2011). In order to take into consideration the probabilistic nature of viral introduction and epidemic outcome, we expanded our framework to include demographic stochasticity and viral extinction (see Methods section).

As demonstrated in Figure 4B, there was a time window of several months during 2012 when adult mosquito counts where sufficiently high and, critically, the virus's temperature-dependent extrinsic incubation period was shorter than the mosquito's average life-expectancy, thus allowing for efficient vector-human transmission. When simulating introduction events, we found significant differences in the epidemic windows between 2012 and a typical year due to deviations in temperature trends throughout the studied periods. Notably, the winter in 2012 was unusually cold, which resulted in a shorter window during which outbreaks could be measured (Figure 5A). At the same time, slightly warmer temperatures during the summer months of 2012 increased the transmission potential, (Figures 5B and 5C), and resulted in bigger outbreaks when compared to a typical year (Figure 5A).

At first sight, differences of 2–4° Celsius in the summer months may seem insufficient to explain the differences in and consequently in outbreak sizes. However, according to experimental evidence, increases in temperature just above the critical point of 20° Celsius will strongly add to the overall vectorial capacity of Aedes mosquitoes [16], [20], [23], [32], [38]. This is a consequence of slight changes in the rates describing mortality, incubation and life-stage progression, which in concert have a cumulative effect and may be involved in positive feedback relationships. Hence, differences of a few degrees, especially when maintained over wide periods of time, can have significant and long lasting effects on the vector population and therefore on dengue's reproductive number.

We further investigated dengue's success of invasion into the island by quantifying the frequency of stochastic simulations that developed into outbreaks above certain sizes, differentiating between 2012 (Figure 5B) and ‘typical’ years (Figure 5C). Comparing the occurrence of any-size () or major-size () outbreaks, we found the risk for the latter to be strongly linked to introductions during the summer months. In fact, as demonstrated in Figures 5B and 5C, there is a substantial risk for major epidemic outbreaks for introductory events taking place weeks or even months before reaches its full potential. This is because any introduction during that period can enjoy the climate-driven ‘deterministic’ growth in until late summer. We can thus identify a key epidemic window, dictated by temperatures above Celsius, in which efforts to detect and control imported cases are crucial for public health planning in Madeira.

In agreement with the estimated differences in transmission potential between 2012 and average years, we also found the invasion success to be generally higher during 2012, which could potentially explain the success of the virus in that particular year. In fact, our results suggest that during a typical year a substantial proportion of introductions () are expected to go extinct before reaching epidemic potential, even during the peak in transmission potential (Figure 5C). Given the homogeneous assumptions of our modelling approach, we argue that these rates should be seen as ‘best’ case scenarios for the successful invasion of dengue in Madeira. In a more realistic scenario, in which heterogeneities in contacts and host and vector densities are present, we expect these rates to be potentially much lower, which could offer an explanation as to why dengue had failed to achieve sustained transmission on the island in the past.