# Dengue modeling in rural Cambodia: Statistical performance versus epidemiological relevance

Clara Champagne, Richard E. Paul, Sowath Ly, Veasna Duong, Rithea Leang, Bernard Cazelles

<figure><img src="image_1.png" /></figure>

## To cite this version:

Clara Champagne, Richard E. Paul, Sowath Ly, Veasna Duong, Rithea Leang, et al., Dengue modeling in rural Cambodia: Statistical performance versus epidemiological relevance. Epidemics, 2019, 26, pp.43-57. 10.1016/j.epidem.2018.08.004. hal-02747514

HAL Id: hal-02747514

https://hal.science/hal-02747514

Submitted on 19 Oct 2020

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dengue modeling in rural Cambodia: Statistical performance versus epidemiological relevance

<figure><img src="image_3.png" /></figure>

Check for
updates

Clara Champagne<sup>a,b,*</sup>, Richard Paul<sup>c,d</sup>, Sowath Ly<sup>e</sup>, Veasna Duong<sup>f</sup>, Rithea Leang<sup>g</sup>, Bernard Cazelles<sup>a,d,h,*</sup>

## ARTICLE INFO

### Keywords:

- Model selection

- Dengue

- Vector borne disease

- Mathematical model

- Cambodia

## ABSTRACT

Dengue dynamics are shaped by the complex interplay between several factors, including vector seasonality, interaction between four virus serotypes, and inapparent infections. However, paucity or quality of data do not allow for all of these to be taken into account in mathematical models. In order to explore separately the importance of these factors in models, we combined surveillance data with a local-scale cluster study in the rural province of Kampong Cham (Cambodia), in which serotypes and asymptomatic infections were documented. We formulate several mechanistic models, each one relying on a different set of hypotheses, such as explicit vector dynamics, transmission via asymptomatic infections and coexistence of several virus serotypes. Models are confronted with the observed time series using Bayesian inference, through Markov chain Monte Carlo. Model selection is then performed using statistical information criteria, and the coherence of epidemiological characteristics (reproduction numbers, incidence proportion, dynamics of the susceptible classes) is assessed in each model. Our analyses on transmission dynamics in a rural endemic setting highlight that two-strain models with interacting effects better reproduce the long term data, but they are difficult to parameterize when relying on incidence cases only. On the other hand, considering the available data, incorporating vector and asymptomatic components seems of limited added-value when seasonality and underreporting are already accounted for.

## 1. Introduction

Dengue is a vector-borne viral disease transmitted by *Aedes* spp. caused by any of four dengue virus (DENV) serotypes. Infection can result in a flu-like illness, and sometimes potentially lethal complications called Dengue Hemorrhagic Fever (DHF) and Dengue Shock Syndrome (DSS), although a significant proportion are subclinical or asymptomatic, causing insufficient discomfort for clinical presentation (Grange et al., 2014). Dengue is ubiquitous in the tropics and the subtropics, particularly in Southeast Asia, the Pacific and the Americas (Guzman et al., 2010). The World Health Organization (WHO) considers that dengue is a major public health issue worldwide, with four billion people in 128 countries exposed to the dengue virus (Messina

et al., 2013; WHO, 2017), an estimated 390 million infections every year and about 50–100 million symptomatic cases worldwide and a high disease burden (Bhatt et al., 2013; Shepard et al., 2016). Nowadays, there are more cases of dengue worldwide than any other arboviral disease (Halstead, 2007; Messina et al., 2015; Sharp et al., 2017).

The value of mathematical models and associated statistical tools for investigating public health policy questions has long been recognized and has provided insights into their transmission and control for more than one hundred years (Heesterbeek et al., 2015; Reiner et al., 2013). It is important, however, to adapt them as much as possible to a specific setting, in order to derive appropriate public health recommendations and accurately generate the key parameters using estimation tools, so that they can produce realistic conclusions, in

<footnote><sup>a</sup> Institut de biologie de l'Ecole normale supérieure (IBENS), Ecole normale supérieure, CNRS, INSERM, PSL Research University, 75005 Paris, France</footnote>

<footnote><sup>b</sup> CREST, ENSAE, Université Paris Saclay, 5 avenue Henry Le Chatelier, 91120 Palaiseau, France</footnote>

<footnote><sup>c</sup> Institut Pasteur, Unité de Génétique Fonctionnelle des Maladies Infectieuses, Department of Genomes and Genetics, F-75724 Paris cedex 15, France</footnote>

<footnote><sup>d</sup> Centre National de la Recherche Scientifique (CNRS), Génomique évolutive, modélisation et santé, UMR 2000, 75724 Paris cedex 15, France</footnote>

<footnote><sup>e</sup> Institut Pasteur in Cambodia, Epidemiology and Public Health Unit, Phnom Penh, Cambodia</footnote>

<footnote><sup>f</sup> Institut Pasteur in Cambodia, Virology Unit, Phnom Penh, Cambodia</footnote>

<footnote><sup>g</sup> National Dengue Control Program, Cambodia</footnote>

<footnote><sup>h</sup> International Center for Mathematical and Computational Modeling of Complex Systems (UMMISCO), UMI 209 UPMC/IRD, Bondy cedex, France</footnote>

<footnote>* Corresponding authors at: Institut de biologie de l'Ecole normale supérieure (IBENS), Ecole normale supérieure, CNRS, INSERM, PSL Research University, 75005 Paris, France.</footnote>

<footnote>E-mail addresses: champagn@biologie.ens.fr (C. Champagne), cazelles@biologie.ens.fr (B. Cazelles).</footnote>

<footnote>https://doi.org/10.1016/j.epidem.2018.08.004</footnote>

<footnote>Received 6 November 2017; Received in revised form 19 July 2018; Accepted 27 August 2018</footnote>

<footnote>Available online 29 August 2018</footnote>

<footnote>1755-4365/ © 2018 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license</footnote>

<footnote>(http://creativecommons.org/licenses/BY/4.0/).</footnote>
<table><caption>**Table 1**
Number of cases under 15 years old in the DENFREE study (index cases and community cases) and in the surveillance system NDSS, and associated theoretical quantities in models. One case in 2013 was coinfected with DENV-1 and DENV-2 and was not included in models with two strains (serotypes).</caption><thead><tr><th>Children under 15 years old</th><th>Total number of observed cases</th><th>Theoretical quantity</th><th>Observation rate</th><th>Model</th></tr></thead><tbody><tr><td>DENV3 data</td><td>year 2012</td><td>C<sub>I</sub></td><td>r<sub>D</sub></td><td>All</td></tr><tr><td>Confirmed symptomatic cases</td><td>236</td><td>C<sub>I1</sub></td><td>r<sub>D</sub></td><td>SEIR2, SEIR2psi</td></tr><tr><td>Denv-1</td><td>226</td><td>C<sub>I2</sub></td><td>r<sub>D</sub></td><td>SEIR2, SEIR2psi</td></tr><tr><td>Denv-2, Denv-3, Denv-4</td><td>5</td><td></td><td></td><td></td></tr><tr><td>Unknown serotype</td><td>5</td><td></td><td></td><td></td></tr><tr><td>Confirmed asymptomatic cases</td><td>5</td><td>C<sub>A</sub></td><td>r<sub>A</sub></td><td>SEIAR</td></tr><tr><td>Denv-1</td><td>5</td><td></td><td></td><td></td></tr><tr><td>Denv-2, Denv-3, Denv-4</td><td>0</td><td></td><td></td><td></td></tr><tr><td>NDSS data</td><td>2002–2015</td><td></td><td></td><td></td></tr><tr><td>Surveillance cases</td><td>10,780 (10,096 in 2002–2013)</td><td>C<sub>I</sub></td><td>r<sub>N</sub></td><td>All except SEIAR</td></tr><tr><td></td><td></td><td>C<sub>H</sub></td><td>r<sub>H</sub></td><td>SEIAR</td></tr></tbody></table>

accordance with the observed data.

Dengue dynamics are shaped by the complex interplay between many factors associated with the mosquito vector and human hosts and their interactions with the virus. Hitherto, the exploration of dengue dynamics has focused on the urban setting, where the incidence of dengue is highest (Clapham et al., 2015; Nisalak et al., 2016; Reich et al., 2013; Salje et al., 2012). Few studies have been carried out in rural settings (Aldstadt et al., 2012; Mammen et al., 2008; Strickman et al., 2000), despite growing evidence that rural dengue is an increasing problem. Guha-Sapir and Schimmer (2005) observed shifts in modal age, rural spread, and social determinants of dengue susceptibility, with major implications for health services. Muhammad Azami et al. (2011) observed similar dengue seroprevalence rates between urban and rural samples, showing that dengue is not confined to urban areas in Malaysia. Chareonsook et al. (1999) showed that DHF in Thailand, which was originally thought to be an urban disease, has spread to most areas of Thailand, and is now more common in rural than urban areas and studies suggest that rural dengue incidence can surpass urban and semi-urban communities within the same region (Reller et al., 2012; Vong et al., 2010). In addition, several studies have stressed that rural settings play an important role in the timing of dengue epidemics in Southeast Asia, with the seasonal dengue waves typically arriving later in major urban centers (Cazelles and Cazelles, 2014; Cuong et al., 2013; Teurlai et al., 2012).

In this study, we combine two datasets from rural Cambodia that provide information on different key factors. We contrast and compare several mechanistic models, incorporating differing levels of complexity with respect to vector dynamics, coexistence of several virus strains, and transmission via asymptomatic infections. Models are adapted to the observed time series using Bayesian inference, through Markov Chain Monte Carlo (MCMC) and compared in light of the data, using statistical indicators to identify the best model (Camacho et al., 2011; King et al., 2008; Pandey et al., 2013; Reich et al., 2013). In addition, we also analyze the epidemiological coherence of the estimated models in simulations. Critically, we do not merely focus on the observed infected individuals but also on other compartments, such as the susceptible class of individuals. By comparing these models, we try to find a realistic but parsimonious way of modeling dengue epidemics in rural Cambodia. The best model may then be used in the study of intervention scenarios or in comparative analyses with other settings. For instance, it could be readily expanded to understand the potential impact of different vaccination strategies in rural settings.

# 2. Methods

## 2.1. Data

### 2.1.1. Study area

Kampong Cham province is a densely populated rural province

120 km northeast from the capital Phnom Penh. Dengue is endemic and strongly seasonal, with outbreaks occurring every year from June to September, during the rainy season. The four virus serotypes co-circulate, even though one usually dominates the three others for about 3 to 5 years. We used two different datasets reporting dengue cases in the province: the results of a punctual study conducted in a 30 km radius around the city of Kampong Cham (DENFREE data), and the national surveillance data (NDSS data) in the four districts comprising the DENFREE study area (Kampong Cham, Kampong Siem, Prey Chhor and Tboung Khmum, with the administrative divisions of 2012–2013).

### 2.1.2. DENFREE data

The DENFREE study took place in the Kampong Cham region during the 2012 and 2013 outbreaks. Patients with acute dengue-like illness were enrolled in three hospitals in the Kampong Cham province. Positive DENV cases were considered as index cases, and an outbreak investigation was initiated in their neighbourhood, in order to detect additional symptomatic cases but also asymptomatic or mildly symptomatic cases. For both index and outbreak investigation cases, DENV infection was confirmed by qRT-PCR. The study protocol is extensively detailed in Duong et al. (2015).

We used the series of the total number of cases per week (index cases and outbreak investigation cases) between 18th June–5th November 2012 and 3rd June–30th September 2013 (Table 1). We also restricted the study to children under 15 years old for two major reasons: most of the reported dengue cases were in this age class (89.6% and 90.8% in 2012 and 2013 respectively, cf. Appendix B), and it allowed a comparison with other dengue reporting systems in Cambodia, which are mainly done at paediatric hospitals. Information on the serotype responsible for infection, and symptomatic/asymptomatic status of the patients were available (Table 2).

### 2.1.3. NDSS data

Because the DENFREE data covers only a relatively short period of time, surveillance data were added to improve the estimations. Surveillance of dengue is conducted at the national level in Cambodia, through the National Dengue Surveillance system (NDSS) (Huy et al., 2010; Teurlai et al., 2012), involving the paediatric departments of several hospitals throughout the country. Since surveillance is hospital based, mostly severe cases are observed. Diagnosis is done clinically and only a small fraction of the cases are confirmed serologically. Because of the co-circulation of other flaviviruses (Chikungunya, Japanese Encephalitis) and the relative non-specificity of symptoms, clinical misdiagnosis may be frequent. Since surveillance is carried out in paediatric departments, only cases among children under 16 years old are reported.

We selected all the cases under 15 years old in the four districts involved in the DENFREE study between January 2002 and December 2015 and aggregated them per week (cases under 15 years old
<table><caption>**Table 2**
Comparison of both datasets.</caption><thead><tr><th></th><th>DENFREE</th><th>NDSS</th></tr></thead><tbody><tr><td>Time window</td><td>Observations in 2012–2013</td><td>Observations in 2002–2015</td></tr><tr><td>Area</td><td>No observation during inter-epidemic period, dataset starts during the epidemic peak<br>30 km radius around Kampong Cham city<br>Clustered collecting process<br>An observation rate can be calculated</td><td>Observations all year round<br>4 districts comprising the DENFREE study<br>Stable reporting process over time<br>Unknown observation rate</td></tr><tr><td>Diagnosis</td><td>Laboratory confirmation of all cases (qRT-PCR)<br>Observation of symptomatic infections (index cases and outbreak investigation) and mildly symptomatic or asymptomatic infections (outbreak investigation)<br>Known serotype for both index and outbreak investigation cases<br>Primary and secondary infections are not distinguished</td><td>Clinical diagnosis<br>Symptoms ranging from dengue fever to DHF and DSS (mostly hospitalized)<br>Unknown serotype<br>Primary and secondary infections are not distinguished</td></tr></tbody></table>

represent 96.2% of reported cases over the period, cf. Appendix B). In this area, on average, 770 cases under 15 years old are reported per year (maximum 1985 cases in 2007, minimum 209 cases in 2014). We used data from 2002 to 2013 for estimations, and data for 2014 and 2015 as the test set.

### 2.1.4. Population

We take as the reference population (N = 161, 391) the number of children below 15 years old in four districts of the Kampong Cham province (Kampong Cham, Kampong Siem, Prey Chhor and Tboung Khmum, with the administrative divisions of 2012–2013) according to 2008 National Census (National Institute of Statistics, Ministry of Planning, 2009). Since the DENFREE study was conducted in a subpart of this area, we calculated the total population for the DENFREE study (n = 65, 208) as the sum of the population of children under 15 years old in all the villages investigated in either 2012 or 2013 (National Institute of Statistics, Ministry of Planning, 2009).

## 2.2. Models

All model parameters are defined in the figures captions and in Table 3.

### 2.2.1. One-strain models

We take a Susceptible-Exposed-Infected-Recovered (SEIR) model as the simplest model (cf. Fig. 1). $C_I$ represents the count of new cases and it is aggregated weekly to be compared with data on both NDSS and DENFREE symptomatic cases. In this model, the basic reproduction number, i.e. the number of secondary infections resulting from the introduction of a single infected in an entirely susceptible population, is
$$R_0^{SEIR}(t) = \frac{\beta(t)\sigma}{(\gamma + \mu_H)(\sigma + \mu_H)}$$

$$\begin{aligned} \frac{dH_S}{dt} &= \mu_H N - \beta(t) \frac{(H_I + i)H_S}{N} - \mu_H H_S \\ \frac{dH_E}{dt} &= \beta(t) \frac{(H_I + i)H_S}{N} - \sigma H_E - \mu_H H_E \\ \frac{dH_I}{dt} &= \sigma H_E - \gamma H_I - \mu_H H_I \\ \frac{dH_R}{dt} &= \gamma H_I - \mu_H H_R \\ \frac{dC_I}{dt} &= \sigma H_E \end{aligned} \quad (1)$$

This model is compared with two other models that include the mosquito vector transmission components. In the first one, which is a Ross-Mcdonald type model derived from Pandey et al. (2013), the vector is modelled explicitly with three compartments (Susceptible-Exposed-Infected) (cf. Fig. 2). In the second one, derived from Laneri et al. (2010), the vector is modelled implicitly as an external force of infection including two stages, latent (κ) and current (λ) (cf. Fig. 3). We derived $R_0$ for each model as $R_0^{Pandey}(t) = \frac{\beta_H \beta_V(t) \sigma}{(\gamma + \mu_H)(\sigma + \mu_H) \mu_V (\mu_V + \tau)}$ and the estimation $R_0^{Laneri}(t) = \frac{\beta(t) \sigma}{(\gamma + \mu_H)(\sigma + \mu_H)}$ (Champagne et al., 2016). In order to compare these models with the non-vector models, we considered the same definition (i.e. the number of secondary human infections

resulting from the introduction of a single infected human in a entirely susceptible population), and not the reproduction ratio per generation provided through the use of the next generation matrix.

The equations describing the Pandey model are:

$$\begin{aligned} \frac{dH_S}{dt} &= \mu_H N - \beta_H v_I H_S - \mu_H H_S \\ \frac{dH_E}{dt} &= \beta_H v_I H_S - \sigma H_E - \mu_H H_E \\ \frac{dH_I}{dt} &= \sigma H_E - \gamma H_I - \mu_H H_I \\ \frac{dH_R}{dt} &= \gamma H_I - \mu_H H_R \\ \frac{dv_S}{dt} &= \mu_V - \beta_V(t) \frac{(H_I + i)}{N} v_S - \mu_V v_S \\ \frac{dv_E}{dt} &= \beta_V(t) \frac{(H_I + i)}{N} v_S - \tau v_E - \mu_V v_E \\ \frac{dv_I}{dt} &= \tau v_E - \mu_V v_I \\ \frac{dC_I}{dt} &= \sigma H_E \end{aligned} \quad (2)$$

where $v_S$ is the proportion of susceptible mosquitoes, $v_E$ the proportion of exposed mosquitoes, and $v_I$ the proportion of infected mosquitoes.

The equations describing the Laneri model are:

$$\begin{aligned} \frac{dH_S}{dt} &= \mu_H N - \lambda H_S - \mu_H H_S \\ \frac{dH_E}{dt} &= \lambda H_S - \sigma H_E - \mu_H H_E \\ \frac{dH_I}{dt} &= \sigma H_E - \gamma H_I - \mu_H H_I \\ \frac{dH_R}{dt} &= \gamma H_I - \mu_H H_R \\ \frac{d\chi}{dt} &= \beta(t) \frac{2(H_I + i)\tau}{N} - 2\tau\kappa \\ \frac{d\lambda}{dt} &= 2\tau\kappa - 2\tau\lambda \\ \frac{dC_I}{dt} &= \sigma H_E \end{aligned} \quad (3)$$

### 2.2.2. Model with explicit asymptomatic individuals (SEIAR)

We also consider a model in which asymptomatic infections are explicitly taken into account in the transmission process (cf. Fig. 4). In this model, we assume that, after the incubation period, there are three possible manifestations of the disease: asymptomatic ($H_A$), mildly symptomatic not requiring hospitalization ($H_I$) and hospitalized cases ($H_H$). Asymptomatic cases are defined in the dengue study as asymptomatic or pauci-symptomatic (presence of other symptoms not being sufficient to classify as symptomatic). Hospital cases are defined as NDSS cases (reported by the surveillance system in hospitals). We assume that symptomatic DENFREE cases are either ($H_I$) or ($H_H$). We also assume that asymptomatic cases transmit the disease as much as symptomatic cases, as recently shown (Duong et al., 2015), and therefore, $R_0^{SEIAR}(t) = \frac{\beta(t)\sigma}{(\gamma + \mu_H)(\sigma + \mu_H)}$. $C_H$, $C_I$, $C_A$ represent respectively the count of new hospitalized, symptomatic and asymptomatic cases, and each is aggregated weekly to be compared respectively with NDSS data, DENFREE data on symptomatic cases, and DENFREE data on asymptomatic cases.
<caption>**Table 3**</caption>

**Prior distributions of parameters.** “Uniform[0,20)” indicates a uniform distribution in the range [0,20]. “Normal(0.44, 0.05) in [0.2, 1)” indicates a truncated normal distribution with mean 0.44 and standard deviation 0.05, restricted to the range [0.2,1].

<table><thead><tr><th>Parameter</th><th>Prior distribution</th><th>Reference</th><th>Models</th></tr></thead><tbody><tr><td colspan="4"><strong>Infectiousness, incubation and mortality rates</strong></td></tr><tr><td>γ<sup>-1</sup></td><td>4.5</td><td>WHO (2017)</td><td>All</td></tr><tr><td>σ<sup>-1</sup></td><td>5.9</td><td>Chan and Johansson (2012)</td><td>Pandey, Laneri</td></tr><tr><td>τ<sup>-1</sup></td><td>10</td><td>Chan and Johansson (2012)</td><td>Pandey, Laneri</td></tr><tr><td>σ<sup>-1</sup></td><td>15.9</td><td>Chan and Johansson (2012)</td><td>SEIR, SEIAR, SEIR2, SEIR2psi</td></tr><tr><td>μ<sub>H</sub><sup>-1</sup></td><td>15</td><td>Assumed</td><td>All</td></tr><tr><td>μ<sub>V</sub><sup>-1</sup></td><td>15</td><td>Liu-Helmersson et al. (2014)</td><td>Pandey</td></tr><tr><td>ρ<sub>A</sub></td><td>Uniform[0, 1]</td><td>Assumed</td><td>SEIAR</td></tr><tr><td>ρ<sub>H</sub></td><td>Uniform[0, 1]</td><td>Assumed</td><td>SEIAR</td></tr><tr><td colspan="4"><strong>Transmission parameters</strong></td></tr><tr><td>R<sub>0</sub></td><td>Uniform[0, 20]</td><td>Assumed</td><td>All</td></tr><tr><td>β<sub>V</sub></td><td>Uniform[0.1, 2]</td><td>Pandey et al. (2013)</td><td>Pandey</td></tr><tr><td>ψ</td><td>Uniform[0.5, 3]</td><td>Assumed</td><td>SEIR2psi</td></tr><tr><td colspan="4"><strong>Initial conditions</strong></td></tr><tr><td>H<sub>I</sub>(0)</td><td>Uniform[0, 100]</td><td>Assumed</td><td>SEIR, Pandey, Laneri</td></tr><tr><td>H<sub>E</sub>(0)</td><td>H<sub>I</sub>(0)</td><td>Assumed</td><td>SEIR, Pandey, Laneri, (3H<sub>I</sub>(0) in SEIAR)</td></tr><tr><td>H<sub>S</sub>(0)</td><td>N*Normal(0.44, 0.05) in [0.2, 1]</td><td>Thai et al. (2005)</td><td>All</td></tr><tr><td>v<sub>I</sub>(0) or λ<sub>I</sub>(0)</td><td>0</td><td>Assumed</td><td>Pandey or Laneri</td></tr><tr><td>v<sub>E</sub>(0) or κ<sub>I</sub>(0)</td><td>0</td><td>Assumed</td><td>Pandey or Laneri</td></tr><tr><td>H<sub>H</sub>(0), H<sub>A</sub>(0)</td><td>H<sub>I</sub>(0)</td><td>Assumed</td><td>SEIAR</td></tr><tr><td>H<sub>S1</sub>(0), H<sub>S2</sub>(0)</td><td>N*Uniform[0.01, 0.5]</td><td>Assumed</td><td>SEIR2, SEIR2psi</td></tr><tr><td>H<sub>I1</sub>(0), H<sub>I2</sub>(0)</td><td>Uniform[0, 100]</td><td>Assumed</td><td>SEIR2, SEIR2psi</td></tr><tr><td>H<sub>E1</sub>(0), H<sub>E12</sub>(0), H<sub>E21</sub>(0)</td><td>H<sub>I1</sub>(0)</td><td>Assumed</td><td>SEIR2, SEIR2psi</td></tr><tr><td>H<sub>E2</sub>(0), H<sub>I12</sub>(0), H<sub>E12</sub>(0)</td><td>H<sub>I2</sub>(0)</td><td>Assumed</td><td>SEIR2, SEIR2psi</td></tr><tr><td colspan="4"><strong>Observation process</strong></td></tr><tr><td>r<sub>N</sub></td><td>Uniform[0, 1]</td><td>Assumed</td><td>SEIR, Laneri, Pandey, SEIR2, SEIR2psi</td></tr><tr><td>r<sub>H</sub></td><td>1</td><td>Assumed</td><td>SEIAR</td></tr><tr><td>r<sub>D</sub></td><td>Fixed</td><td>cf. Table 4</td><td>All</td></tr><tr><td>r<sub>A</sub></td><td>Fixed</td><td>cf. Table 4</td><td>SEIAR</td></tr><tr><td>φ</td><td>Uniform[0, 1]</td><td>Assumed</td><td>All</td></tr><tr><td colspan="4"><strong>Seasonality parameters</strong></td></tr><tr><td>b</td><td>Uniform[0, 1]</td><td>Assumed</td><td>All</td></tr><tr><td>p</td><td>Uniform[-0.5, 0.5]</td><td>Assumed</td><td>All</td></tr><tr><td>i</td><td>Uniform[0, 10]</td><td>Assumed</td><td>All</td></tr><tr><td colspan="4"><strong>Total population</strong></td></tr><tr><td>N</td><td>161391</td><td>Assumed</td><td>All</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Graphical representation of SEIR model. $H_S$ susceptible individuals; $H_E$ infected (not yet infectious) individuals; $H_I$ infectious individuals; $H_R$ recovered individuals; $\beta(t)$ is the transmission parameter; $\sigma$ is the rate at which $H_E$-individuals move to the infectious class $H_I$; infectious individuals ($H_I$) then recover at rate $\gamma$; individuals leave the children population at rate $\mu_H$. $H_S + H_E + H_I + H_R = N$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Graphical representation of Pandey model (Pandey et al., 2013). Squared boxes and circles correspond respectively to human and vector compartments. Plain arrows represent transitions from one state to the next. Dashed arrows indicate interactions between humans and vectors. $H_S$ susceptible individuals; $H_E$ infected (not yet infectious) individuals; $H_I$ infectious individuals; $H_R$ recovered individuals; $\beta_H$ is the transmission parameter from vector to human; $\sigma$ is the rate at which $H_E$-individuals move to the infectious class $H_I$; infectious individuals ($H_I$) then recover at rate $\gamma$; individuals leave the children population at rate $\mu_H$; $H_S + H_E + H_I + H_R = N$; $v_S$ proportion of susceptible vectors; $v_E$ proportion of infected (not yet infectious) vectors; $v_I$ proportion of infectious vectors; $\beta_V(t)$ is the transmission parameter from human to vector; $\tau$ is the rate at which $v_E$-vectors move to the infectious class $v_I$; vectors die at rate $\mu_V$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 3. Graphical representation of Laneri model (Laneri et al., 2010). Squared boxes and circles correspond respectively to human and vector compartments. Plain arrows represent transitions from one state to the next. Dashed arrows indicate interactions between humans and vectors. Hs susceptible individuals; He infected (not yet infectious) individuals; Hj infectious individuals; HR recovered individuals; σ is the rate at which Hs-individuals move to the infectious class Hs; infectious individuals (Hj) then recover at rate γ; individuals leave the children population at rate μH; Hs + He + Hj + HR = N; implicit vector-borne transmission is modelled with the compartments κ and λ; λ current force of infection; κ latent force of infection reflecting the exposed state for mosquitoes during the extrinsic incubation period; β(t) is the transmission parameter; τ is the transition rate associated with the extrinsic incubation period.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. Graphical representation of SEIR model. Hs susceptible individuals; He infected (not yet infectious) individuals; HA asymptomatic infectious individuals; Hj mildly symptomatic infectious individuals; Hs hospitalized infectious individuals; HR recovered individuals; β(t) is the transmission parameter; σ is the rate at which He individuals move to the infectious classes Hs, HA and Hj; a proportion ρA of He individuals do not show symptoms during the infectious period; a proportion ρH of symptomatic individuals go to hospital; infectious individuals (Hs, HA, Hs) then recover at rate γ; individuals leave the children population at rate μH. Hs + He + HA + Hj + HR + Hs = N.</figcaption></figure>

$$
\begin{align}
\frac{dH_S}{dt} &= \mu_H N - \beta(t) \frac{(H_1 + H_A + H_2 + i)H_S}{N} - \mu_H H_S \\
\frac{dH_E}{dt} &= \beta(t) \frac{(H_1 + H_A + H_2 + i)H_S}{N} - \sigma H_E - \mu_H H_E \\
\frac{dH_A}{dt} &= \rho_A \sigma H_E - \gamma H_A - \mu_H H_A \\
\frac{dH_H}{dt} &= \rho_H (1 - \rho_A) \sigma H_E - \gamma H_H - \mu_H H_H \\
\frac{dH_I}{dt} &= (1 - \rho_H)(1 - \rho_A) \sigma H_E - \gamma H_I - \mu_H H_I \\
\frac{dH_R}{dt} &= \gamma (H_I + H_A + H_H) - \mu_H H_R \\
\frac{dC_I}{dt} &= (1 - \rho_A) \sigma H_E \\
\frac{dC_A}{dt} &= \rho_A \sigma H_E \\
\frac{dC_H}{dt} &= \rho_H (1 - \rho_A) \sigma H_E
\end{align}
\tag{4}
$$

2.2.3. Model with two virus serotypes

In the 2012 and 2013 epidemics, DENV-1 was highly dominant: the three other serotypes represented about 2% of the cases reported in the DENFREE study in 2012 and about 21% in 2013 (cf. Table 1). Therefore, a two-strain model is also studied, in which we separate DENV-1 cases from DENV-2, DENV-3 and DENV-4 combined (cf. Fig. 5). For simplicity and parsimony in the number of parameters, the two strains share the same parameter values. We first assume both strains to be independent (ψ = 1 in Eq. (5), called SEIR2 model). In this context, the

reproduction numbers for each strain are equal,
$R_0^{\text{SEIR21}}(t) = R_0^{\text{SEIR22}}(t) = \frac{\beta(t)\sigma}{(\gamma + \mu_H)(\sigma + \mu_H)}$.

$$
\begin{align}
\frac{dH_S}{dt} &= \mu_H N - \beta(t) \frac{(H_{11} + \psi H_{12} + i)H_S}{N} - \beta(t) \frac{(H_{12} + \psi H_{12} + i)H_S}{N} - \mu_H H_S \notag \\
\frac{dH_{E1}}{dt} &= \beta(t) \frac{(H_{11} + \psi H_{12} + i)H_S}{N} - \sigma H_{E1} - \mu_H H_{E1} \notag \\
\frac{dH_{I1}}{dt} &= \sigma H_{E1} - \gamma H_{I1} - \mu_H H_{I1} \notag \\
\frac{dH_{S1}}{dt} &= \gamma H_{I1} - \beta(t) \frac{(H_{12} + \psi H_{12} + i)H_{S1}}{N} - \mu_H H_{S1} \notag \\
\frac{dH_{E12}}{dt} &= \beta(t) \frac{(H_{12} + \psi H_{12} + i)H_{S1}}{N} - \sigma H_{E12} - \mu_H H_{E12} \notag \\
\frac{dH_{I12}}{dt} &= \sigma H_{E12} - \gamma H_{I12} - \mu_H H_{I12} \notag \\
\frac{dH_{E2}}{dt} &= \beta(t) \frac{(H_{12} + \psi H_{12} + i)H_S}{N} - \sigma H_{E2} - \mu_H H_{E2} \notag \\
\frac{dH_{I2}}{dt} &= \sigma H_{E2} - \gamma H_{I2} - \mu_H H_{I2} \notag \\
\frac{dH_{S2}}{dt} &= \gamma H_{I2} - \beta(t) \frac{(H_{11} + \psi H_{12} + i)H_{S2}}{N} - \mu_H H_{S2} \notag \\
\frac{dH_{E21}}{dt} &= \beta(t) \frac{(H_{11} + \psi H_{12} + i)H_{S2}}{N} - \sigma H_{E21} - \mu_H H_{E21} \notag \\
\frac{dH_{I21}}{dt} &= \sigma H_{E21} - \gamma H_{I21} - \mu_H H_{I21} \notag \\
\frac{dH_R}{dt} &= \gamma (H_{I12} + H_{I21}) - \mu_H H_R \notag \\
\frac{dC_I}{dt} &= \sigma (H_{E1} + H_{E2} + H_{E21} + H_{E12}) \notag \\
\frac{dC_{I1}}{dt} &= \sigma (H_{E1} + H_{E21}) \notag \\
\frac{dC_{I2}}{dt} &= \sigma (H_{E2} + H_{E12}) \tag{5}
\end{align}
$$

We also considered another version of the model including inter-
action between strains (ten Bosch et al., 2016), in order to reflect the
fact that secondary infection with a heterologous serotype leads more
often than primary infection to severe manifestations of the disease
(Halstead, 2007). In our model (called SEIR2psi model), primary and
secondary infections differ in infectiousness, through a parameter ψ
(Ferguson et al., 1999). This parameter is estimated between 0.5 and 3:
values superior to 1 correspond to transmission cross-enhancement
(because of higher virus titers during secondary infections (Ferguson
et al., 1999)) and values inferior to 1 suggest a lower infectivity for
secondary infected individuals (for example because they are hospita-
lized and less in contact with the population (Aguiar et al., 2011)). As in
Ferguson et al. (1999), we define R₀<sup>SEIR21</sup>(t) = R₀<sup>SEIR22</sup>(t) = <sup>β(t)σ</sup>&frasl;<sub>(γ + μ<sub>H</sub>)(σ + μ<sub>H</sub>)</sub>
the basic reproduction number for each strain.

C<sub>I</sub> represents the count of new cases for both serotypes and it is
aggregated weekly to be compared with NDSS data. C<sub>I1</sub> and C<sub>I2</sub> re-
present the count of new cases for strain 1 and strain 2, and are ag-
gregated weekly to be compared with symptomatic cases from the
DENFREE study for DENV-1 and DENV-2/ DENV-3/ DENV-4 respec-
tively.

2.2.4. Seasonality

All models include seasonality through the use of a time-varying transmission parameter $\beta(t) = \beta \left[ 1 + b \cdot \sin \left( 2\pi \left( \frac{t}{365} + p \right) \right) \right]$, according to a sinusoidal function whose phase $p$ and amplitude $b$ are estimated.

We also assume that a constant number of cases *i* are imported.

2.3. Prior distributions

The prior distributions are listed in Table 3.

Dirac priors based on the literature were used for the durations of infectiousness and incubation, as well as the mortality rates. In models without vectorial transmission, the incubation period is assumed to be the sum of the extrinsic (in mosquito) and intrinsic (in human) incubation periods, to reflect the generation time of the disease. For transmission parameters, we used wide weakly informative priors.

2.3.1. Initial conditions

The initial number of infected individuals is assumed to be equal to
<figure><img src="image_8.png" /></figure>

parameter; σ is the rate at which exposed individuals move to the infectious class; infectious individuals then recover at rate γ; ψ is the change in infectivity for secondary infected individuals in SEIR2psi model (in SEIR2 model, ψ = 1); individuals leave the children population at rate μH. Hs + HE1 + HE2 + H11 + H12 + Hs1 + HS2 + HE12 + HE21 + H112 + H121 + HR = N.

the number of exposed individuals and to be lower than 100, as the model starts in January, during the epidemic trough. Except for the initial proportion of susceptibles, all priors on initial conditions are uniform distributions.

The initial proportion of susceptibles is an influential parameter on the model outputs. It is highly correlated to the transmission parameter β (and therefore to the basic reproductive number), which makes it difficult to estimate them both. An informative gaussian prior was therefore used on Hs(0). To date, no large scale seroprevalence study is available for Cambodia, and we relied on a study conducted among schoolchildren in rural Vietnam (Thai et al., 2005), which we considered as the closest setting to be compared with Kampong Cham. We extrapolated their results on schoolchildren (7 to 14 years old) to a 1–15 years old population as follows, where Sλ is the proportion of susceptibles among 1–15 year-old children (using their estimation λ = 0.117):

$$S_{\lambda} = \frac{1}{15} \sum_{a=1}^{15} \exp(-\lambda a) = 0.44$$

Sλ is used as the mean of the gaussian prior, and the standard deviation is fixed at 0.05.

## 2.4. Estimation

### 2.4.1. Observation model

The propagation models are related to the observed data using a negative binomial observation model (Bretó et al., 2009). The observed number of cases during week k in each dataset, C(k)obs, is assumed to be drawn from a a negative binomial distribution with mean rC(k) and variance rC(k) + (rC(k))²φ, where C(k) is the total number of new cases simulated by the model during week k, r is the observation rate quantifying the amount of non reported cases and φ is an overdispersion parameter. The quantity C and observation rate r associated to each dataset are indicated in Table 1.

For each epidemiological model, estimations are performed simultaneously on NDSS and DENFREE data (DENFREE data being decomposed into two time series in SEIAR, SEIR2 and SEIR2psi models). The observations of all time series are assumed independent conditional on the underlying disease process, and the model likelihood is the product of the likelihoods calculated on each series.

The observation rate is estimated for NDSS data but fixed for DENFREE data, using informations from the sampling scheme (cf. Table 4), in order to account for the difference in coverage during outbreak investigation between 2012 and 2013. We assume that index cases are all reported and that the observation rate for community cases equals the ratio of people tested over the population of the area. We then extrapolate this observation rate to the total population of the four districts.

Different observation rates are used in the SEIAR model. As NDSS data are interpreted as hospitalized cases, the observation rate is fixed to 1 (which assumes that all hospitalized cases go through surveillance

**Fig. 5.** Graphical representation of SEIR2 models. Hs individuals susceptible to both strains; HE1 (resp. HE2) individuals infected (not yet infectious) to strain 1 (resp. strain 2); H11 (resp. H12) individuals infectious to strain 1 (resp. strain 2); Hs1 (resp. HS2) individuals immune to strain 1 only (resp. strain 2); HE12 (resp. HE21) individuals (not yet infectious) with a secondary infection to strain 2 (resp. strain 1); H112 (resp. H121) infectious individuals with a secondary infection to strain 2 (resp. strain 1); HR individuals immune to both strains; β(t) is the transmission

and neglects the presence of private hospitals or non reports from hospitals), but the proportion of hospitalized cases is estimated. For DENFREE cases, the observation rate for symptomatic and asymptomatic individuals is fixed as indicated in Table 4.

The overdispersion parameter φ is estimated along with the other parameters.

### 2.4.2. Estimation with Markov Chain Monte Carlo

Models are considered in deterministic framework and estimations are made using random walk Metropolis Hastings. SSM software (Dureau et al., 2013) is used for simulations and calculations. The covariance matrix of the proposal distribution was initialized using adaptive MCMC as in Dureau et al. (2013).

Due to unknown initial conditions, the posterior distribution is multimodal, especially in two strain models. This implies that the adaptive MCMC remains trapped in local maxima and global exploration of the parameter space was performed with latin hypercube sampling. We ran a simplex algorithm on 10,000 parameter sets sampled with latin hypercube sampling (with lhs R package (Carnell, n.d.)) and chose the one with the highest posterior value as initialization of the adaptive MCMC algorithm. The results associated to this best fitting equilibrium are displayed in the main text, but other examples of good fitting equilibria for two-strain models are indicated in Appendix E.

## 2.5. Model comparison

### 2.5.1. Statistical indicators

In order to identify the best model, the Deviance Information Criterion (DIC) (Spiegelhalter et al., 2002) is used. DIC is an indicator that combines a measure of model fit and a penalty on model complexity, commonly used with MCMC estimations. The best model is the one with the smallest DIC. As it does not enable comparison of models with differing number of observations, we also calculate the RMSE (root mean square error) between observations C(k)obs and simulations of the model with parameters sampled in the posterior distribution C(k,i), for weeks k = k0, ..., K and simulation i:

$$RMSE_i = \sqrt{\frac{1}{K - k_0} \sum_{k=k_0}^{K} (C^{(k,i)} - C_{\text{obs}}^{(k)})^2}$$

We then compute the mean and the 2.5% and 97.5% quantiles of this indicator among the simulated trajectories. We first calculated them on the data used for estimation, separating NDSS data and DENFREE data. For DENFREE data we considered the total number of symptomatic cases to enable comparison of all models: therefore, for this metric, asymptomatic DENFREE cases were not taken into account in SEIAR model, and cases of all serotypes were aggregated in two-strain models. Then we used these indicators to assess the predictive performance of the models, comparing projections of the model with NDSS observations for 2014 and 2015.
<table><caption>**Table 4**
Calculation of the observation rate for DENFREE data.</caption><thead><tr><th rowspan="2">Children &lt; 15 years old</th><th colspan="2">2012</th><th colspan="2">2013</th></tr><tr><th>Symptomatic</th><th>Asymptomatic</th><th>Symptomatic</th><th>Asymptomatic</th></tr></thead><tbody><tr><td>Population in investigated villages (n)</td><td>65,208</td><td></td><td>65,208</td><td></td></tr><tr><td>Population in the 4 districts (N)</td><td>161,391</td><td></td><td>161,391</td><td></td></tr><tr><td>Index cases (a)</td><td>151</td><td>0</td><td>376</td><td>0</td></tr><tr><td>Observation rate for index cases</td><td>1</td><td>.</td><td>1</td><td>.</td></tr><tr><td>Children tested in communities (b)</td><td>1722</td><td></td><td>4119</td><td></td></tr><tr><td>DENV positive in communities (c)</td><td>85</td><td>5</td><td>198</td><td>28</td></tr><tr><td>Observation rate for community cases (d = b/n)</td><td>0.0264</td><td></td><td>0.0632</td><td></td></tr><tr><td>Extrapolated number of community cases (e = c/d)</td><td>3219</td><td>189</td><td>3135</td><td>443</td></tr><tr><td>Extrapolated total number of cases (f = a + e)</td><td>3370</td><td>189</td><td>3511</td><td>443</td></tr><tr><td>Observed total number of cases (g = a + c)</td><td>236</td><td>5</td><td>574</td><td>28</td></tr><tr><td><strong>Observation rate (r<sub>D</sub> = (g/f) * (n/N) and r<sub>A</sub> = b/N)</strong></td><td><strong>0.0283</strong></td><td><strong>0.0107</strong></td><td><strong>0.0661</strong></td><td><strong>0.0255</strong></td></tr></tbody></table>

### 2.5.2. Epidemiological indicators

Models are also compared according to several indicators to describe their epidemiological behaviour. The basic reproduction number, the observation rate and the initial proportion of susceptibles are estimated using the MCMC chain. In the model with asymptomatic infections we report the estimated proportions of asymptomatic and hospitalized cases.

With parameters sampled in the posterior distribution, we can also re-simulate the model to study hidden states, such as the susceptible and infected classes. The effective reproduction number ($R_e$) is calculated as the seasonal basic reproduction number multiplied by the proportion of susceptibles at each time step, as indicated in Table 5. We then calculate the annual incidence proportion as the total number of infections over one year divided by the total population of susceptibles at the beginning of the year. In models with two strains, we separate the annual incidence of primary infections (as the total number of primary infections over one year divided by the total population of naïve individuals at the beginning of the year) and secondary infections (as the total number of secondary infections over one year divided by the total population of susceptibles to one strain only at the beginning of the year).

Calculations are made using R version 3.2.2 (Team, 2015), and graphics using ggplot2 (Wickham, 2009).

## 3. Results

The calculations are based on the MCMC chain associated with the highest posterior after LHS exploration. In multistrain models, the posterior distribution is multimodal, and examples of other local posterior modes are displayed in Appendix E.

### 3.1. Statistical comparison

For single strain models, the SEIR model proved the best with the DIC criterion (cf. Table 6). For two strain models, the SEIR2psi proved best. As regards simulation-based indicators on the 2002–2014 data, SEIR, Laneri, Pandey models have RMSE values in the same order of magnitude. Indeed, they produce a similar dynamic with respect to the 2002–2014 data (cf. Figs. 6–7), with a period of approximately six

years, and a large overestimation of the 2002 outbreak (the mean RMSE for 2002 is far higher (> 40) than the average for the other years (< 25)). Due to the small number of observed asymptomatic cases, SEIAR model also produces a similar dynamic and the SEIR2 model, despite a larger number of parameters, does not improve over single-strain models. The model with two interacting strains (SEIR2psi) outperforms the other models in terms of RMSE, but the difference is mainly explained by the first year of simulation (2002). However, the RMSE indicator does not penalize for the complexity of this more flexible model, and the differences are not very pronounced given the variability in RMSE values. When visualizing the simulations compared to the data (cf. Figs. 6–7), all models underestimate a large number of epidemic peaks and all models but the one with two interacting strains overpredict the first epidemic peak (as pointed out in Table 6). SEIR2psi is the model in which the large epidemic in 2007 is best reproduced and overall, SEIR2psi is the model that reproduces most accurately the observed data. Moreover, the estimated overdispersion is lower in this model (cf. Table 7), indicating that less observation noise is required to explain the data and that the model explains a larger part of the variability in the data.

Considering the predictive capacity of the models, all models overestimate the 2014 small epidemic, and to a lesser extent the one in 2015 (cf. Fig. 8). In the SEIR2psi model, this overestimation is less pronounced, and the associated RMSE is also smaller (cf. Table 6).

### 3.2. Epidemiological comparison

The average $R_0$ is estimated to be between 2 and 3 in most of the models and the maximum value between 3 and 4 (cf. Table 7), except with the Pandey model, in which it is higher (mean value above 3 and maximum value above 6). The estimated values are very close in the SEIR, Laneri and SEIAR models. These values relate to transmission among children only and we expect therefore the $R_0$ in the whole population to be higher, since children represent approximately one third of the total population, but more than 90% of the total cases observed. Estimates for $R_0$ in South East Asia in the literature range from 1 to 5 (Imai et al., 2015, 2016), with sometimes values above 10 (Johansson et al., 2011), displaying large variations in time and space, and in particular, the estimates for $R_0$ in Cambodia based on age-stratified

<table><caption>**Table 5**
Reproduction numbers calculation in each model. Calculations are based on the next generation matrix method (Diekmann et al., 2010).</caption><thead><tr><th></th><th>SEIR, Laneri, SEIAR</th><th>Pandey</th><th>SEIR2</th><th>SEIR2psi</th></tr></thead><tbody><tr><td>R<sub>0</sub>(t)</td><td><div style="text-align: center;"><span style="display: inline-block; border-bottom: 1px solid black; padding-bottom: 1px;">β(t)σ</span><span style="display: inline-block; padding-top: 1px;">(γ + μ<sub>H</sub>)(σ + μ<sub>H</sub>)</span></div></td><td><div style="text-align: center;"><span style="display: inline-block; border-bottom: 1px solid black; padding-bottom: 1px;">β<sub>H</sub>β̂<sub>V</sub>(t)σ</span><span style="display: inline-block; padding-top: 1px;">(γ + μ<sub>H</sub>)(σ + μ<sub>H</sub>)μ<sub>V</sub>(μ<sub>V</sub> + τ)</span></div></td><td><div style="text-align: center;"><span style="display: inline-block; border-bottom: 1px solid black; padding-bottom: 1px;">β(t)σ</span><span style="display: inline-block; padding-top: 1px;">(γ + μ<sub>H</sub>)(σ + μ<sub>H</sub>)</span></div></td><td><div style="text-align: center;"><span style="display: inline-block; border-bottom: 1px solid black; padding-bottom: 1px;">β(t)σ</span><span style="display: inline-block; padding-top: 1px;">(γ + μ<sub>H</sub>)(σ + μ<sub>H</sub>)</span></div></td></tr><tr><td>R<sub>e</sub>(t) (or R<sub>e</sub><sup>i</sup>, i = 1, 2)</td><td>R<sub>0</sub>(t) <div style="text-align: center;"><span style="display: inline-block; border-bottom: 1px solid black; padding-bottom: 1px;">H<sub>S</sub>(t)</span><span style="display: inline-block; padding-top: 1px;">N</span></div></td><td>R<sub>0</sub>(t) <div style="text-align: center;"><span style="display: inline-block; border-bottom: 1px solid black; padding-bottom: 1px;">H<sub>S</sub>(t)</span><span style="display: inline-block; padding-top: 1px;">N</span></div></td><td>R<sub>0</sub>(t) <div style="text-align: center;"><span style="display: inline-block; border-bottom: 1px solid black; padding-bottom: 1px;">H<sub>S</sub>(t) + H<sub>Sj</sub>(t)</span><span style="display: inline-block; padding-top: 1px;">N</span></div></td><td>R<sub>0</sub>(t) <div style="text-align: center;"><span style="display: inline-block; border-bottom: 1px solid black; padding-bottom: 1px;">H<sub>S</sub>(t) + ψH<sub>Sj</sub>(t)</span><span style="display: inline-block; padding-top: 1px;">N</span></div></td></tr></tbody></table>
<caption>**Table 6**</caption>

**Information criteria in deterministic models.** DIC is the Deviance Information Criterion (Spiegelhalter et al., 2002). RMSE is the root mean square error (mean and 2.5% and 97.5% quantiles) between simulations and observations: it is calculated separately on the datasets used for estimations (NDSS data for 2002–2013 and DENFREE data for 2012–2013) and on the test set (NDSS data for 2014–2015). It is also computed on separated years in order to highlight well or badly estimated years: for example, for each simulation $i$, $\frac{1}{12}(\text{RMSE}_i^{\text{NDSS 2002}})^2 + \frac{11}{12}(\text{RMSE}_i^{\text{NDSS 2003–2013}})^2 = (\text{RMSE}_i^{\text{NDSS}})^2$. Calculations are based on the MCMC chain associated with the highest posterior.

<table><thead><tr><th>Model</th><th>SEIR</th><th>Laneri</th><th>Pandey</th><th>SEIAR</th><th>SEIR2</th><th>SEIR2psi</th></tr></thead><tbody><tr><td>nb parameters</td><td>8</td><td>8</td><td>9</td><td>9</td><td>11</td><td>12</td></tr><tr><td>nb observations</td><td>665</td><td>665</td><td>665</td><td>704</td><td>704</td><td>704</td></tr><tr><td colspan="7"><strong>ESTIMATION SET</strong></td></tr><tr><td>DIC</td><td>4343</td><td>4360</td><td>4364</td><td>4472</td><td>4492</td><td>4446</td></tr><tr><td>RMSE NDSS</td><td>25 (22–28)</td><td>25 (22–29)</td><td>26 (23–29)</td><td>26 (23–29)</td><td>25 (22–29)</td><td>21 (19–24)</td></tr><tr><td>RMSE NDSS 2002</td><td>43 (28–59)</td><td>42 (27–58)</td><td>45 (29–63)</td><td>43 (27–61)</td><td>40 (24–61)</td><td>17 (12–25)</td></tr><tr><td>RMSE NDSS 2003–2013</td><td>23 (20–25)</td><td>23 (21–26)</td><td>23 (21–25)</td><td>23 (21–26)</td><td>24 (21–26)</td><td>22 (19–24)</td></tr><tr><td>RMSE DENFREE</td><td>24 (19–29)</td><td>24 (19–29)</td><td>24 (19–29)</td><td>24 (19–29)</td><td>23 (19–27)</td><td>25 (19–31)</td></tr><tr><td colspan="7"><strong>TEST SET</strong></td></tr><tr><td>RMSE 2014–2015</td><td>15 (11–19)</td><td>15 (11–19)</td><td>16 (12–20)</td><td>15 (11–20)</td><td>17 (13–22)</td><td>13 (10–16)</td></tr><tr><td>RMSE 2014</td><td>16 (11–23)</td><td>16 (11–23)</td><td>17 (11–24)</td><td>17 (11–24)</td><td>17 (11–24)</td><td>14 (10–19)</td></tr><tr><td>RMSE 2015</td><td>13 (9–17)</td><td>13 (9–17)</td><td>13 (10–18)</td><td>13 (10–17)</td><td>17 (12–23)</td><td>12 (9–16)</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Number of observed cases per week and NDSS data, 2002–2013. Simulations with negative binomial noise using parameters from the MCMC chain associated with the highest posterior, calculated using both NDSS and DENFREE datasets. Posterior median (solid line), 95% credible intervals (shaded blue area) and NDSS data points (black dots). A. SEIR model. B. Laneri model. C. Pandey model. D. SEIAR model. E. SEIR2 model. F. SEIR2psi model.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Number of observed cases per week and DENFREE data, 2012–2013. Simulations with negative binomial noise using parameters from the MCMC chain associated with the highest posterior, calculated using both NDSS and DENFREE datasets. Posterior median (solid line), 95% credible intervals (shaded blue area) and DENFREE data points (blue dots: symptomatic cases, green dots: asymptomatic cases, red dots: serotype specific symptomatic cases). A. SEIR model. B. Laneri model. C. Pandey model. D. SEIAR model, symptomatic cases. E. SEIAR model, asymptomatic cases. F. SEIR2 model, cases due to strain 1 and DENV1 observations. G. SEIR2 model, cases due to strain 2 and DENV2/DENV3/DENV4 observations. H. SEIR2 model with interactions, cases due to strain 1 and DENV1 observations. I. SEIR2 model with interactions, cases due to strain 2 and DENV2/DENV3/DENV4 observations. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

case-notification data are between 2 and 7 (Imai et al., 2016), so that our estimates are in the same order of magnitude. Nevertheless, our estimation of $R_0$ strongly depends on the estimation of the initial proportion of susceptibles ($H_S(0)$), which is unknown in the case of Cambodia and can bias the estimates. It is indeed very difficult to evaluate $R_0$ via incidence data (Imai et al., 2016), especially for endemic diseases without information on seroprevalence, and it may be more informative to study the effective reproductive number ($R_e$), which accounts for the proportion of immune individuals in the population. The estimated $R_e$ has a mean value around 1 and a maximum value around 1.7 in most models, but it is higher in the Pandey model (the time-varying behaviour of $R_e$ is displayed in Appendix D).

In SEIR2psi model, the parameter $\psi$ quantifying the interaction between strains is inferior to 1, suggesting a reduced infectivity of secondary infections on average, as in Aguiar et al. (2011), or in Coudeville and Garnett (2012). This suggests that, as far as

infectiousness is concerned, cross-protection among serotypes is more important than cross-enhancement to explain the results observed in the field.

It is estimated in the models that approximately half of the children are susceptible to the disease, which is close to the informative prior used, highlighting the difficulty to estimate this quantity due to identifiability issues. This proportion is lower in the Pandey model. In the models with two strains, the number of susceptibles differs between models, but as a whole, at least 50% of children are not totally immune to the disease and can still suffer infection. The proportion of children who are susceptible to one or both strains are however correlated in the MCMC chain (cf. Appendix), indicating that their relative proportions are not well identified. These values are in the range of the measures of seroprevalence in several Asian countries (L'Azou et al., 2016, 2018). As the measures reveal large differences between countries (L'Azou et al., 2016, 2018; Prayitno et al., 2017; Thai et al., 2005;
<table><caption>**Table 7**
Epidemiological criteria among children under 15 years old. Estimated parameters and indicators based on simulations over 2002–2015 (posterior median and 95% credible intervals). Calculations are based on the MCMC chain associated with the highest posterior.</caption><thead><tr><th>Model</th><th></th><th>SEIR</th><th>Laneri</th><th>Pandey</th><th>SEIAR</th><th>SEIR2</th><th>SEIR2 psi</th></tr></thead><tbody><tr><td>Mean R<sub>0</sub></td><td>Median (95%CI)</td><td>2.32 (2.25–2.41)</td><td>2.38 (2.3–2.46)</td><td>3.32 (3.16–3.47)</td><td>2.3 (2.23–2.37)</td><td>1.45 (1.42–1.48)</td><td>2.47 (2.35–2.57)</td></tr><tr><td>Max R<sub>0</sub></td><td>Median (95%CI)</td><td>3.64 (3.51–3.77)</td><td>3.66 (3.53–3.8)</td><td>5.87 (5.57–6.17)</td><td>3.6 (3.48–3.72)</td><td>2.29 (2.23–2.35)</td><td>3.78 (3.59–3.94)</td></tr><tr><td>ψ</td><td>Median (95%CI)</td><td></td><td></td><td></td><td></td><td></td><td>0.67 (0.62–0.73)</td></tr><tr><td>H<sub>5</sub>(0)/N (%)</td><td>Median (95%CI)</td><td>46 (45–47)</td><td>46 (45–47)</td><td>37 (36–38)</td><td>47 (46–47)</td><td>49 (45–56)</td><td>35 (30–38)</td></tr><tr><td>H<sub>51</sub>(0)/N (%)</td><td>Median (95%CI)</td><td></td><td></td><td></td><td></td><td>22 (14–25)</td><td>14 (10–21)</td></tr><tr><td>H<sub>52</sub>(0)/N (%)</td><td>Median (95%CI)</td><td></td><td></td><td></td><td></td><td>27 (19–30)</td><td>4 (1–11)</td></tr><tr><td>Observation rate (%)</td><td>Median (95%CI)</td><td>12 (11–13)</td><td>12 (11–13)</td><td>10 (9–11)</td><td></td><td>11 (10–12)</td><td>7 (6–7)</td></tr><tr><td>Hospitalized (%)</td><td>Median (95%CI)</td><td></td><td></td><td></td><td>14 (13–15)</td><td></td><td></td></tr><tr><td>Asymptomatic (%)</td><td>Median (95%CI)</td><td></td><td></td><td></td><td>15 (13–19)</td><td></td><td></td></tr><tr><td>Over-dispersion</td><td>Median (95%CI)</td><td>0.73 (0.59–0.91)</td><td>0.75 (0.61–0.93)</td><td>0.7 (0.57–0.87)</td><td>0.9 (0.74–0.99)</td><td>0.73 (0.59–0.91)</td><td>0.56 (0.45–0.76)</td></tr><tr><td>Median annual incidence proportion</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Primary infection (%)</td><td>Median 2002–2015 (min–max)</td><td>8 (5–16)</td><td>8 (4–16)</td><td>12 (7–25)</td><td>8 (4–16)</td><td>6 (3–10)</td><td>14 (7–27)</td></tr><tr><td>Secondary infection (%)</td><td>Median 2002–2015 (min–max)</td><td></td><td></td><td></td><td></td><td>3 (1–5)</td><td>7 (3–14)</td></tr><tr><td>Mean R<sub>e</sub></td><td>Median (95%CI)</td><td>1.03 (1.03–1.03)</td><td>1.06 (1.06–1.07)</td><td>1.14 (1.13–1.15)</td><td>1.03 (1.03–1.03)</td><td></td><td></td></tr><tr><td>Mean R<sub>e</sub> strain 1</td><td>Median (95%CI)</td><td></td><td></td><td></td><td></td><td>1.03 (1.02–1.03)</td><td>1.03 (1.02–1.03)</td></tr><tr><td>Mean R<sub>e</sub> strain 2</td><td>Median (95%CI)</td><td></td><td></td><td></td><td></td><td>1.04 (1.03–1.04)</td><td>1.02 (1.02–1.03)</td></tr><tr><td>Max R<sub>e</sub></td><td>Median (95%CI)</td><td>1.7 (1.66–1.74)</td><td>1.72 (1.69–1.77)</td><td>2.18 (2.12–2.27)</td><td>1.7 (1.67–1.74)</td><td></td><td></td></tr><tr><td>Max R<sub>e</sub> strain 1</td><td>Median (95%CI)</td><td></td><td></td><td></td><td></td><td>1.73 (1.69–1.77)</td><td>1.74 (1.71–1.77)</td></tr><tr><td>Max R<sub>e</sub> strain 2</td><td>Median (95%CI)</td><td></td><td></td><td></td><td></td><td>1.72 (1.69–1.75)</td><td>1.69 (1.65–1.73)</td></tr></tbody></table>

Vongpunsawad et al., 2017), a seroprevalence survey in Cambodia would be particularly useful to evaluate which scenario is more plausible.

The observation rate for NDSS data varies between models, from 6 to 13%. These values indicate that a large proportion of dengue cases are not reported in national surveillance, likely reflecting mild symptoms that do not require hospitalization or misdiagnosis or misreporting.

In the SEIR2psi model, we also plotted the current number of infected individuals for each strain (cf. Fig. 9). In our simulations, the first strain is responsible for large explosive outbreaks, whereas the second one has a more regular pattern over the years. Moreover, the two strains are asynchronous and each one dominates for two or three years. It is also qualitatively close to the dynamics observed in Thailand (Reich et al., 2013) or Singapore (Ong et al., 2017).

The proportion of susceptible individuals also displays asynchronous dynamics between strains, as they reflect the history of past epidemics. Despite the seasonality and the year-to-year variations, the total number of susceptibles remains high (cf. Fig. 9 for SEIR2psi and SEIR models), allowing the possibility for large outbreaks to occur in the future. However, multi-strain models are difficult to parameterize due to the uncertainty in the initial conditions, and models with differing trajectories and parameters may have close performances in terms of model fitting (cf. Appendix E).

We calculated the annual incidence proportion as the proportion of new infections over one year among the susceptibles at the beginning of that year (cf. Fig. 10). The values for primary infections are coherent with other studies in Vietnam or Indonesia who analyzed seroprevalence data or seroconversion data (Graham et al., 1999; Prayitno et al., 2017; Thai et al., 2005, 2007; Tien et al., 2010). The incidence proportion is highly variable from one year to the next, especially in models with two interacting strains.

## 4. Discussion

With two datasets reporting dengue cases in the Kampong Cham region in Cambodia, we compared several models to represent dengue transmission dynamics in a rural setting. In order to assess the quality of the models, we compared their statistical properties and analysed their

epidemiological features.

The best model describing the dengue trend over 14 years of data was the two strain model, with reduced infectivity for secondary infected individuals. Secondary infections being more prone to severe dengue (Halstead, 2007), these individuals may stay at home or at the hospital, and be less involved in the spread of the disease than the ones whose mild symptoms do not alter what Stoddard et al. (2013) call “house-to-house movement” (Yoon et al., 2012; Perkins et al., 2016). There is also some indication from studies on the infectiousness of infected people to mosquitoes, that symptomatic infections are less infectious than asymptomatic individuals (Duong et al., 2015). Likewise there was a small negative impact of IgG on transmission rates that would be indicative of secondary and potentially more symptomatic cases (Nguyen et al., 2013). This feature was previously analytically studied and coherent with dengue incidence time series in Thailand (Aguiar et al., 2011, 2014).

On the contrary, including vectorial transmission or a compartment for asymptomatic infections did not seem to improve the model fit despite the additional complexity. In our models, the role of mosquitoes in transmission seems well captured by the seasonal forcing implemented in the model. The limited added-value of including vectorial components has been observed by several authors previously, when studying one dengue season only (Pandey et al., 2013) or when a seasonal forcing was included in the model (Rocha et al., 2016). Mathematical analyses have suggested that because the time scale of the mosquito epidemiology is so fast compared to that in humans, it will be slaved by the slower human epidemiology (Rocha et al., 2013). Thus, for understanding human disease epidemiology, mainly the dynamics of the human time scale are essential and inclusion of mosquito dynamics results in an unnecessary increase in model complexity when vector data is not available (Dye and Williams, 1995; Rocha et al., 2013, 2016), and when potential seasonal variations are already accounted for.

Our models account for underreporting, and the lack of additional improvement when including explicitly the asymptomatic class is likely due to the very few asymptomatic infections observed, which may be due to a very strict definition of asymptomatic infections in the DENFREE protocol. Yoon et al. (2012) also observe many inapparent cases in their cohort but few strictly asymptomatic cases in their cluster
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Projections of the number of observed cases per week and NDSS data, 2014–2015. Simulations with negative binomial noise using parameters from the MCMC chain associated with the highest posterior, calculated using both NDSS and DENFREE datasets. Posterior median (solid line), 95% credible intervals (shaded blue area) and NDSS data points (black dots). A. SEIR model. B. Laneri model. C. Pandey model. D. SEIAR model. E. SEIR2 model. F. SEIR2psi model. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

study in Thailand. Therefore, in our model, including a compartment for asymptomatic individuals had only little influence on the overall transmission due to their small number.

However, even though the two strain model with interactions reproduces the observed data more adequately, it is difficult to parameterize. Differing model parameterizations may lead to differing dynamics that relate closely to the observed data, as shown in Appendix E, indicating identifiability issues (Miao et al., 2011). Identifiability is nonetheless a central problem in designing and estimating epidemiological models, that also affects models with vector-borne transmission (Kao and Eisenberg, 2018) and even simple models (Evans et al., 2002). In our case, this problem was related to two phenomena. On the one hand, some parameters values were correlated, and therefore difficult to distinguish from one another. On the other hand, in some models, the posterior distribution was multimodal, and therefore both hard to explore with many inference methods and hard to interpret given the presence of several possible trajectories. In order to limit these difficulties in all models, we fixed some parameters (in particular the durations of infection and incubation, or the observation rate in

DENFREE data), and relied on informative priors (in particular for the initial proportion of susceptibles). In addition, it is possible to rule out some trajectories that correctly fit the data but display epidemiological incoherences in the qualitative dynamics of their hidden states (such as the depletion of the susceptible population or the absence of serotype asynchrony), in the spirit of pattern matching methods (ten Bosch et al., 2016). We therefore studied the epidemiological characteristics of the selected models to ensure such coherence. Other trajectories may, nevertheless, remain equally plausible given the available information and require additional data to be further explored. These issues are even more important when evaluating public health measures, as models with similar fit may have different responses to control scenarios (Kao and Eisenberg, 2018; ten Bosch et al., 2016). Therefore, care is required when choosing the best model, and the uncertainty in the model structure must also be taken into account, especially when the observed data are scarce.

Over the six models that were explored, we also obtained some insight on the parameters describing transmission. The average annual $R_0$ in our paediatric population is estimated between 1 and 4 over all
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Number of infected and proportion of susceptible individuals. Simulations without negative binomial noise using parameters from the MCMC chain associated with the highest posterior, calculated using both NDSS and DENFREE datasets. Posterior median (solid line) and 95% credible intervals (shaded blue area). A. Individuals infected with strain 1 in SEIR2psi model ($H_{11} + H_{121}$). B. Individuals infected with strain 2 in SEIR2psi model ($H_{12} + H_{121}$). C. Proportion of susceptible individuals in SEIR model ($H_s/N$). D. Proportion of susceptible individuals with immunity to strain 1 in SEIR2psi model ($H_{s1}/N$). E. Proportion of susceptible individuals with immunity to strain 2 in SEIR2psi model ($H_{s2}/N$). F. Proportion of individuals susceptible to both strains in SEIR2psi model ($H_s/N$). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

models. These values are within the range observed in urban settings, suggesting that, despite very different population densities, the rural dynamics of dengue are not that dissimilar, or that dissimilarities are hidden by the variations between countries and populations and the uncertainties due to diverse estimation procedures, as this quantity is difficult to evaluate using incidence data. The median annual incidence at primary infection over the period is between 5 and 14, with large year-to-year variations. The estimated observation rate on surveillance data varies between models (9–13% in most models and 6–7% in the model with two interacting strains), indicating in both cases a high proportion of unreported infections, so that most of the transmission is due to unobserved infections. These values are in line with the large underrecognition highlighted in Wichmann et al. (2011), and also with other studies in South-East Asia (Graham et al., 1999; Thai et al., 2007;

Yoon et al., 2012). These results are also coherent with a recent study emphasizing that more than 80% of dengue transmission is due to inapparent infections (ten Bosch et al., 2018).

This work has however several limitations. Firstly, model comparison was not straightforward between one strain and two strain models, both statistically and epidemiologically. From the statistical point of view, the differing number of observations between models led us to use simulation based-indicators. From the epidemiological point of view, single strain models, two strain models and observations on four serotypes may be hard to compare because some indicators cover different interpretations. For example, in single strain models, there is no distinction in the susceptible class between individuals immune to one strain only and naive ones, and there is no strain specific reproduction number or incidence.

<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Annual incidence proportion of first dengue infection (%). Median and credible intervals per year, over 2002–2015, based on simulations without negative binomial noise using parameters from the MCMC chain associated with the highest posterior, calculated using both NDSS and DENFREE datasets.</figcaption></figure>
Secondly, the selected model formulations were restricted due to data availability. In particular, despite the endemicity of the four serotypes in rural Cambodia, we did not consider more than two dengue serotypes. This was done to limit model complexity, especially in the number of unknown initial conditions, but has also been previously shown to adequately describe dengue dynamics (Aguiar et al., 2013). When two serotypes were considered, we tested only interactions in terms of enhancement or restriction of infectiousness. We did not include models with (temporary) cross immunity, because of the too large increase in the number of parameters with respect to the data. We also did not include models with a finer spatial scale, even if small scale transmission plays a decisive role in dengue dynamics (Salje et al., 2012, 2017). On the one hand, NDSS data were only available at the district level, which was too large to follow transmission chains and too small for observing a sufficient number of cases. On the other hand, the clustered sampling protocol in the DENFREE study abnegated the interpretation of the spatial distribution of community cases. We also restricted the analysis to children under 15 years old and did not study the role of adults in transmission.

Thirdly, the projections are not completely able to describe the observed data, as most models overestimate the dengue epidemic in 2014. Nevertheless, 2014 was a particular year, with the lowest number of cases in the whole time series, potentially due to particular climatic conditions. In many countries in South-East Asia, except Malaysia, the reported incidence was lower than in 2013 (Cheng et al., 2017). Many provinces of Thailand also reported fewer cases than usual in 2014 (Reich et al., 2016). Our models are deterministic and do not take into account variations due to demographic stochasticity or environmental hazards such as climate.

Despite these limitations, combining two datasets permitted us to overcome some observation biases, such as the fact that surveillance data did not report serotype and DENFREE data did not reflect the long term dynamics. Nevertheless, some information is lacking in both datasets, in particular that on seroprevalence. Clearly the parameter estimations depend on the immunological status at the beginning of the simulations, especially in complex multi-strain models. As in our previous work (Champagne et al., 2016), our modeling study stresses the importance of seroprevalence data in order to more accurately estimate the initial conditions of our simulation and reduce identifiability problems. A seroprevalence survey in Cambodia would be of great value to evaluate the dengue burden, transmissibility potential and consider vaccination scenarios.

## 5. Conclusion

In conclusion, our analyses highlight that two-strain models with interacting effects better reproduce the long term dengue dynamics, but they are also difficult to parameterize relying on incidence data only. On the other hand, incorporating vector and asymptomatic components seems of limited added-value in this case when seasonality and underreporting are already accounted for. Although model complexity is framed by the scientific objectives, it must also be driven by the available data. The unavailability of mosquito data and the difficulty of observing asymptomatic infections questions their incorporation explicitly in the models. Another important aspect is related to the comparison of models considering the available data. In addition to model selection based on goodness of fit (Pandey et al., 2013), assessing the validity of model's outputs in terms of epidemiological features (such as reproduction numbers, annual incidence, dynamics of the susceptible classes) is an important step when data are scarce and identifiability issues are present.

## Acknowledgments

We thank the members of the Epidemiology and Public Health unit and the Virology unit in Institut Pasteur in Cambodia for their work

during the DENFREE project. All the authors have been supported by a grant from Agence Nationale de la Recherche for the PANIC project (Targeting PAthogen's NIChE: a new approach for infectious diseases control in low-income countries: ANR-14-CE02-0015-01). The research leading to these results has also received funding from the European Commission Seventh Framework Program FP7 for the DENFREE project under Grant Agreement 282 378. CC and BC are partially supported by the "Pépiniere interdisciplinaire Eco-Evo-Devo" from the Centre National de la Recherche Scientifique (CNRS). CC is supported by the PhD studentship of the Groupe des Ecoles Nationales d'Economie et Statistique (GENES). The funders played no role in the study design, data collection, analysis, or preparation of the manuscript.

## Appendix A. Supplementary data

Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.epidem.2018.08.004.

## References

Aguiar, M., Ballesteros, S., Kooi, B.W., Stollenwerk, N., 2011. The role of seasonality and import in a minimalistic multi-strain dengue model capturing differences between primary and secondary infections: complex dynamics and its implications for data analysis. J. Theor. Biol. 289, 181–196. https://doi.org/10.1016/j.jtbi.2011.08.043.

Aguiar, M., Kooi, B.W., Rocha, F., Ghaffari, P., Stollenwerk, N., 2013. How much complexity is needed to describe the fluctuations observed in dengue hemorrhagic fever incidence data? Ecol. Complex. 16, 31–40. https://doi.org/10.1016/j.ecocom.2012.09.001. http://www.sciencedirect.com/science/article/pii/S1476945X12000670.

Aguiar, M., Paul, R., Sakuntabhai, A., Stollenwerk, N., 2014. Are we modelling the correct dataset? Minimizing false predictions for dengue fever in Thailand. Epidemiol. Infect. 142 (11), 2447–2459. https://doi.org/10.1017/S09502688130003348. http://journals.cambridge.org/article/S09502688130003348.

Aldstadt, J., Yoon, I.-K., Tannitisupawong, D., Jarman, R.G., Thomas, S.J., Gibbons, R.V., Uppapong, A., Iamsirithaworn, S., Rothman, A.L., Scott, T.W., Endy, T., 2012. Space-time analysis of hospitalised dengue patients in rural Thailand reveals important temporal intervals in the pattern of dengue virus transmission. Trop. Med. Int. Health: TM & IH 17 (9), 1076–1085. https://doi.org/10.1111/j.1365-3156.2012.03040.x.

Bhatt, S., Gething, P.W., Brady, O.J., Messina, J.P., Farlow, A.W., Moyes, C.L., Drake, J.M., Brownstein, J.S., Hoen, A.G., Sankoh, O., Myers, M.F., George, D.B., Jaenisch, T., Wint, G.R.W., Simmons, C.P., Scott, T.W., Farrar, J.J., Hay, S.I., 2013. The global distribution and burden of dengue. Nature 496 (7446), 504–507. https://doi.org/10.1038/nature12060. https://www.nature.com/articles/nature12060.

Bretó, C., He, D., Ionides, E.L., King, A.A., 2009. Time series analysis via mechanistic models. Ann. Appl. Stat. 3 (1), 319–348. https://doi.org/10.1214/08-AOAS201. https://projecteuclid.org/euclid.aos/1239888373.

Camacho, A., Ballesteros, S., Graham, A.L., Carrat, F., Ratmann, O., Cazelles, B., 2011. Explaining rapid reinfections in multiple-wave influenza outbreaks: Tristan da Cunha 1971 epidemic as a case study. Proc. R. Soc. Lond. B: Biol. Sci. https://doi.org/10.1098/rspb.2011.0300.. rspb20110300, http://rspb.royalsocietypublishing.org/content/early/2011/04/22/rspb.2011.0300.

R. Carnell, lhS: Latin hypercube samples. https://CRAN.R-project.org/package-lhs.

Cazelles, B., Cazelles, K., 2014. Major urban centers have weak influence on the timing of dengue epidemics in Southeast Asia. Int. J. Infect. Dis. 21 (Suppl 1), 217. https://doi.org/10.1016/j.ijid.2014.03.873. http://www.sciencedirect.com/science/article/pii/S1201971214009321.

Champagne, C., Salthouse, D.G., Paul, R., Cao-Lormeau, V.-M., Roche, B., Cazelles, B., 2016. Structure in the variability of the basic reproductive number (R0) for Zika epidemics in the Pacific islands. eLife 5, e19874. https://doi.org/10.7554/eLife.19874. https://elifesciences.org/articles/19874.

Chan, M., Johansson, M.A., 2012. The incubation periods of Dengue viruses. PLOS ONE 7 (11), e50972. https://doi.org/10.1371/journal.pone.0050972.

Chareonsook, O., Foy, H.M., Teeraratkul, A., Silarug, N., 1999. Changing epidemiology of dengue hemorrhagic fever in Thailand. Epidemiol. Infect. 122 (1), 161–166.

Cheng, Q., Jing, Q., Spear, R.C., Marshall, J.M., Yang, Z., Gong, P., 2017. The interplay of climate, intervention and imported cases as determinants of the 2014 dengue outbreak in Guangzhou. PLOS Negl. Trop. Dis. 11 (6), e0005701. https://doi.org/10.1371/journal.pntd.0005701.

Clapham, H., Cummings, D.A.T., Nisalak, A., Kalayanarooj, S., Thaisomboonsuk, B., Klungthong, C., Fernandez, S., Srikiatkhachorn, A., Macareo, L.R., Lessler, J., Reiser, J., Yoon, I.-K., 2015. Epidemiology of infant dengue cases illuminates serotype-specificity in the interaction between immunity and disease, and changes in transmission dynamics. PLOS Negl. Trop. Dis. 9 (12), e0004262. https://doi.org/10.1371/journal.pntd.0004262.

Coudeville, L., Garnett, G.P., 2012. Transmission dynamics of the four dengue serotypes in Southern Vietnam and the potential impact of vaccination. PLOS ONE 7 (12), e51244. https://doi.org/10.1371/journal.pone.0051244.

Cuong, H.Q., Vu, N.T., Cazelles, B., Boni, M.F., Thai, K.T., Rabaa, M.A., Quang, L.C., Simmons, C.P., Huu, T.N., Anders, K.L., 2013. Spatiotemporal dynamics of dengue epidemics, Southern Vietnam. Emerg. Infect. Dis. 19 (6), 945–953. https://doi.org/
10.3201/eid1906.121232. http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3713821/.

Diekmann, O., Heesterbeek, J.A.P., Roberts, M.G., 2010. The construction of next-generation matrices for compartmental epidemic models. J. R. Soc. Interface 7 (47), 873–885. https://doi.org/10.1098/rsif.2009.0386. http://rsif.royalsocietypublishing.org/content/7/47/873.

Duong, V., Lambrechts, L., Paul, R.E., Ly, S., Lay, R.S., Long, K.C., Huy, R., Tarantola, A., Scott, T.W., Sakuntabhai, A., Buchy, P., 2015. Asymptomatic humans transmit dengue virus to mosquitoes. Proc. Natl. Acad. Sci. https://doi.org/10.1073/pnas.1508114112.

Dureau, J., Ballesteros, S., Bogich, T., 2013. ssm: Inference for State Space Models like playing with duplo blocks, original-date: 2014-08-29T13:15:10Z (Nov.). https://github.com/JDureau/ssm.

Dye, C., Williams, B., 1995. Nonlinearities in the dynamics of indirectly-transmitted infections (or, does having a vector make a difference? In: Grenfell, B.T., Grenfell, B.T., Dobson, A.P., Dobson, A.P. (Eds.), Ecology of Infectious Diseases in Natural Populations. Cambridge University Press Press, pp. 260–279 http://www.core/books/ecology-of-infectious-diseases-in-natural-populations/nonlinearities-in-the-dynamics-of-indirectly-transmitted-infections-or-does-having-a-vector-make-a-difference/1E571DC1187028C2262A1D6EFE5A73A4.

Evans, N., Chapman, M., Chappell, M., Godfrey, K., 2002. The structural identifiability of a general epidemic (SIR) model with seasonal forcing. IFAC Proc. Vol. 35 (1), 109–114. https://doi.org/10.3182/20020721-6-ES-1901.01327. http://linkinghub.elsevier.com/retrieve/pii/S1474667015397482.

Ferguson, N., Anderson, R., Gupta, S., 1999. The effect of antibody-dependent enhancement on the transmission dynamics and persistence of multiple-strain pathogens. Proc. Natl. Acad. Sci. USA 96 (2), 790–794. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC15215/.

Graham, R.R., Fuffrie, M., Tan, R., Hayes, C.G., Laksono, I., Ma'roef, C., Erlin, n., Sutaryo, n., Porter, K.R., Halstead, S.B., 1999. A prospective seroepidemiologic study on dengue in children four to nine years of age in Yogyakarta, Indonesia I. studies in 1995–1996. Am. J. Trop. Med. Hyg. 61 (3), 412–419.

Grange, L., Simon-Lorine, E., Sakuntabhai, A., Gresh, L., Paul, R., Harris, E., 2014. Epidemiological risk factors associated with high:global frequency of inapparent dengue virus: infections. Front. Immunol. 5(https://doi.org/10.3389/fimmu.2014.00280. http://www.ncbi.nlm.nih.gov/pmc/articles/PMC4052743/.

Guha-Sapir, D., Schimmer, B., 2005. Dengue fever: new paradigms for a changing epidemiology. Emerg. Themes Epidemiol. 2, 1. https://doi.org/10.1186/1742-7622-1-2.

Guzman, M.G., Halstead, S.B., Artsob, H., Buchy, P., Farrar, J., Gubler, D.J., Hunsperger, E., Kroeger, A., Margolis, H.S., Martínez, E., Nathan, M.B., Pelegroin, J.L., Simmons, C., Yoksan, S., Peeling, R.W., 2010. Dengue: a continuing global threat. Nat. Rev. Microbiol. 8, S7–S16. https://doi.org/10.1038/nrmicro2460. http://www.nature.com/nrmicro/journal/v8/n12/full/nrmicro2460.html.

Halstead, S.B., 2007. Dengue. Lancet (London, England) 370 (9599), 1644–1652. https://doi.org/10.1016/S0140-6736(07)61687-0.

Heesterbeek, H., Anderson, R.M., Andreasen, V., Bansal, S., De Angelis, D., Dye, C., Eames, K.T.D., Edmunds, W.J., Frost, S.D.W., Funk, S., Hollingsworth, T.D., House, T., Isham, V., Klepac, P., Lessler, J., Lloyd-Smith, J.O., Metcalf, C.J.E., Mollison, D., Pellis, L., Pulliam, J.R.C., Roberts, M.G., Viboud, C., Isaac Newton Institute IDD Collaboration, 2015. Modeling infectious disease dynamics in the complex landscape of global health. Science 347 (6227), aaa4339. https://doi.org/10.1126/science. aaa4339.

Huy, R., Buchy, P., Conan, A., Ngan, C., Ong, S., Ali, R., Duong, V., Yit, S., Ung, S., Te, V., Chroeung, N., Pheaktra, N.C., Uok, V., Vong, S., 2010. National dengue surveillance in Cambodia 1980–2008: epidemiological and virological trends and the impact of vector control. Bull. World Health Organ. 88 (9), 650–657. https://doi.org/10.2471/BLT.09.073908. http://www.who.int/bulletin/volumes/88/9/09-073908.pdf.

Imai, N., Dorigatti, I., Cauchemez, S., Ferguson, N.M., 2015. Estimating dengue transmission intensity from sero-prevalence surveys in multiple countries. PLoS Negl. Trop. Dis. 9 (4). https://doi.org/10.1371/journal.pntd.0003719. http://www.ncbi.nlm.nih.gov/pmc/articles/PMC4400108/.

Imai, N., Dorigatti, I., Cauchemez, S., Ferguson, N.M., 2016. Estimating dengue transmission intensity from case-notification data from multiple countries. PLOS Negl. Trop. Dis. 10 (7), e0004833. https://doi.org/10.1371/journal.pntd.0004833.

Johansson, M.A., Hombach, J., Cummings, D.A., 2011. Models of the impact of dengue vaccines: a review of current research and potential approaches. Vaccine 29 (35), 5860–5868. https://doi.org/10.1016/j.vaccine.2011.06.042. http://www.ncbi.nlm.nih.gov/pmc/articles/PMC4327892/.

Kao, Y.-H., Eisenberg, M.C., 2018. Practical unidentifiability of a simple vector-borne disease model: implications for parameter estimation and intervention assessment. Epidemics. https://doi.org/10.1016/j.epidem.2018.05.010. http://www.sciencedirect.com/science/article/pii/S1755436517301627.

King, A.A., Ionides, E.L., Pascual, M., Bouma, M.J., 2008. Inapparent infections and cholera dynamics. Nature 454 (7206), 877–880. https://doi.org/10.1038/nature07084. http://www.nature.com/doifinder/10.1038/nature07084.

Laneri, K., Bhadra, A., Ionides, E.L., Bouma, M., Dhiman, R.C., Yadav, R.S., Pascual, M., 2010. Forcing versus feedback: epidemic malaria and monsoon rains in Northwest India. PLoS Comput. Biol. 6 (9), e1000898. https://doi.org/10.1371/journal.pcbi.1000898.

L'Azou, M., Moureau, A., Sarti, E., Nealon, J., Zambrano, B., Wartel, T.A., Villar, L., Capeding, M.R., Ochiai, R.L., 2016. Symptomatic dengue in children in 10 Asian and Latin American Countries. New Engl. J. Med. 374 (12), 1155–1166. https://doi.org/10.1056/NEJMoa1503877.

L'Azou, M., Assoukpa, J., Fanouillere, K., Plennevaux, E., Bonaparte, M., Boukennooghe, A., Frago, C., Noriega, F., Zambrano, B., Ochiai, R.L., Guy, B., Jackson, N., 2018. Dengue seroprevalence: data from the clinical development of a tetravalent dengue

vaccine in 14 countries (2005–2014). Trans. R. Soc. Trop. Med. Hyg. https://doi.org/10.1093/trstmh/try037.

Liu-Helmersson, J., Stenlund, H., Wilder-Smith, A., Rocklöv, J., 2014. Vectorial Capacity of *Aedes aegypti*: effects of temperature and implications for global dengue epidemic potential. PLOS ONE 9 (3), e89783. https://doi.org/10.1371/journal.pone.0089783.

Mammen Jr., M.P., Pimgale, C., Koenraadt, C.J., Rothman, A.L., Aldstadt, J., Nisalak, A., Jarman, R.G., Jones, J.W., Srikiatkhachorn, A., Ypil-Butac, C.A., et al., 2008. Spatial and temporal clustering of dengue virus transmission in Thai villages. PLoS Med. 5 (11), e205. http://journals.plos.org/plosmedicine/article?id=10.1371/journal.pmed.0050205.

Messina, J.P., Brady, O.J., Scott, T.W., Zou, C., Pigott, D.M., Duda, K.A., Bhatt, S., Katzelnick, L., Howes, R.E., Battle, K.E., Simmons, C.P., Hay, S.I., 2014. Global spread of dengue virus types: mapping the 70 year history. Trends Microbiol. 22 (3), 138–146. https://doi.org/10.1016/j.tim.2013.12.011.

Messina, J.P., Brady, O.J., Pigott, D.M., Golding, N., Kraemer, M.U.G., Scott, T.W., Wint, G.R.W., Smith, D.L., Hay, S.I., 2015. The many predicted futures of dengue. Nat. Rev. Microbiol. 13 (4), 230–239. https://doi.org/10.1038/nrmicro3430.

Miao, H., Xia, X., Perelson, A.S., Wu, H., 2011. On identifiability of nonlinear ODE models and applications in viral dynamics. SIAM Rev. Soc. Ind. Appl. Math. 53 (1), 3–39. https://doi.org/10.1137/090757009. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3140286/.

Muhammad Azami, N.A., Salleh, S.A., Neoh, H.-m., Syed Zakaria, S.Z., Jamal, R., 2011. Dengue epidemic in Malaysia: not a predominantly urban disease anymore. BMC Res. Notes 4, 216. https://doi.org/10.1186/1756-0500-4-216.

Nguyen, N.M., Kien, D.T.H., Tuan, T.V., Quyen, N.T.H., Tran, C.N.B., Thi, L.V., Thi, D.L., Nguyen, H.L., Farrar, J.J., Holmes, E.C., Rabaa, M.A., Bryant, J.E., Nguyen, T.T., Nguyen, H.T.C., Nguyen, L.T.H., Pham, M.P., Nguyen, H.T., Luong, T.T.H., Wills, B., Nguyen, C.V.V., Wolbers, M., Simmons, C.P., 2013. Host and viral features of human dengue cases shape the population of infected and infectious *Aedes aegypti* mosquitoes. Proc. Natl. Acad. Sci. USA 110 (22), 9072–9077. https://doi.org/10.1073/pnas.1303395110. http://www.pnas.org/content/110/22/9072.

Nisalak, A., Clapham, H.E., Kalayanarooj, S., Klungthong, C., Thaisomboonsuk, T., Fernandez, S., Reiser, J., Srikiatkhachorn, A., Macroreo, L.R., Lessler, J.T., Cummings, D.A.T., Yoon, I.-K., 2016. Forty years of dengue surveillance at a Tertiary Pediatric Hospital in Bangkok, Thailand, 1973–2012. Am. J. Trop. Med. Hyg. 94 (6), 1342–1347. https://doi.org/10.4269/ajtmh.15-0337.

Ong, J., Yap, G., Hapuarachchi, C., Tang, C.S., Liang, S., Rajarethinam, J., Tan, S., Koo, C., Zeng, K., Ycasas, J., Tay, L., Chong, C.-S., Abdul Razak, M.A., Lee, C., Chiang, S., Lai, Y.-L., Ng, L.-C., 2017. Dengue review of 2016 and outlook for 2017. Singap. Epidemiol. News Bull 43 (2), 47–56.

National Institute of Statistics, Ministry of Planning, 2009. Phnom Penh, Cambodia: General Population Census of Cambodia 2008. National Report on Final Census Results.

Pandey, A., Mubayi, A., Medlock, J., 2013. Comparing vector-host and SIR models for dengue transmission. Math. Biosci. 246 (2), 252–259. https://doi.org/10.1016/j.mbs.2013.10.007. http://linkinghub.elsevier.com/retrieve/pii/S002556413002435.

Perkins, T.A., Paz-Soldan, V.A., Stoddard, S.T., Morrison, A.C., Forshey, B.M., Long, K.C., Halsey, E.S., Kochel, T.J., Elder, J.P., Kitron, U., Scott, T.W., Vazquez-Prokopec, G.M., 2016. Calling in sick: impacts of fever on intra-urban human mobility. Proc. R. Soc. B 283 (1834), 20160390. https://doi.org/10.1098/rspb.2016.0390. http://rspb.royalsocietypublishing.org/content/283/1834/20160390.

Prayitno, A., Taurel, A.-F., Nealon, J., Satari, H.I., Karyanti, M.R., Sekartini, R., Soedjatmiko, S., Gunardi, H., Medise, B.E., Sasmono, R.T., Simmerman, J.M., Boukennooghe, A., Hadinegoro, S.R., 2017. Dengue seroprevalence and force of primary infection in a representative population of urban dwelling Indonesian children. PLOS Negl. Trop. Dis. 11 (6), e0005621. https://doi.org/10.1371/journal.pntd.0005621.

Reich, N.G., Shrestha, S., King, A.A., Rohani, P., Lessler, J., Kalayanarooj, S., Yoon, I.-K., Gibbons, R.V., Burke, D.S., Cummings, D.A.T., 2013. Interactions between serotypes of dengue highlight epidemiological impact of cross-immunity. J. R. Soc. Interface 10 (86), 20130414. https://doi.org/10.1098/rsif.2013.0414.

Reich, N.G., Lauer, S.A., Sakrejda, K., Iamsirithaworn, S., Hinjoy, S., Suangtho, P., Suthachana, S., Clapham, H.E., Salje, H., Cummings, D.A.T., Lessler, J., 2016. Challenges in real-time prediction of infectious disease: a case study of dengue in Thailand. PLOS Negl. Trop. Dis. 10 (6), e0004761. https://doi.org/10.1371/journal.pntd.0004761.

Reiner, R.C., Perkins, T.A., Barker, C.M., Niu, T., Chaves, L.F., Ellis, A.M., George, D.B., Le Menach, A., Pulliam, J.R.C., Bisanzio, D., Buckee, C., Chiyaka, C., Cummings, D.A.T., Garcia, A.J., Gatton, M.L., Gething, P.W., Hartley, D.M., Johnston, G., Klein, E.Y., Michael, E., Lindsay, S.W., Lloyd, A.L., Pigott, D.M., Reisen, W.K., Ruktanonchai, N., Singh, B.K., Tatem, A.J., Kitron, U., Hay, S.I., Scott, T.W., Smith, D.L., 2013. A systematic review of mathematical models of mosquito-borne pathogen transmission: 1970–2010. J. R. Soc. Interface 10 (81), 20120921. https://doi.org/10.1098/rsif.2012.0921.

Reller, M.E., Bodinayake, C., Nagahawatte, A., Devasiri, V., Kodikara-Arachichi, W., Strouse, J.J., Broadwater, A., Östbye, T., de Silva, A., Woods, C.W., 2012.119062. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3310451/.

Rocha, F., Aguiar, M., Souza, M., Stollenwerk, N., 2013. Time-scale separation and centre manifold analysis describing vector-borne disease dynamics. Int. J. Comput. Math. 90 (10), 2105–2125. https://doi.org/10.1080/00207160.2013.783208.

Rocha, F., Mateus, L., Skwara, U., Aguiar, M., Stollenwerk, N., 2016. Understanding dengue fever dynamics: a study of seasonality in vector-borne disease models. Int. J. Comput. Math. 93 (8), 1405–1422. https://doi.org/10.1080/00207160.2015.
1050961.

Salje, H., Lessler, J., Endy, T.P., Curriero, F.C., Gibbons, R.V., Nisalak, A., Nimmanitya, S., Kalayanarooj, S., Jarman, R.G., Thomas, S.J., Burke, D.S., Cummings, D.A.T., 2012. Revealing the microscale spatial signature of dengue transmission and immunity in an urban population. Proc. Natl. Acad. Sci. USA 109 (24), 9535–9538. https://doi.org/10.1073/pnas.1120621109.

Salje, H., Lessler, J., Berry, I.M., Melendrez, M.C., Endy, T., Kalayanarooj, S., A-Nuegoonpipat, A., Chanama, S., Sangkijporn, S., Klungthong, C., Thaisomboonsuk, B., Nisalak, A., Gibbons, R.V., Iamsirithaworn, S., Macareo, L.R., Yoon, I.-K., Sangarsang, A., Jarman, R.G., Cummings, D.A.T., 2017. Dengue diversity across spatial and temporal scales: local structure and the effect of host population size. Science 355 (6331), 1302–1306. https://doi.org/10.1126/science.aaj9384. http://science.sciencemag.org/content/355/6331/1302.

Sharp, T.M., Tomashek, K.M., Read, J.S., Margolis, H.S., Waterman, S.H., 2017. A new look at an old disease: recent: insights into the global epidemiology of dengue. Curr. Epidemiol. Rep. 4 (1), 11–21. https://doi.org/10.1007/s40471-017-0095-y.

Shepard, D.S., Undurraga, E.A., Halasa, Y.A., Stanaway, J.D., 2016. The global economic burden of dengue: a systematic analysis. Lancet Infect. Dis. 16 (8), 935–941. https://doi.org/10.1016/S1473-3099(16)00146-8.

Spiegelhalter, D.J., Best, N.G., Carlin, B.P., Van Der Linde, A., 2002. Bayesian measures of model complexity and fit. J. R. Stat. Soc.: Ser. B (Stat. Methodol.) 64 (4), 583–639. https://doi.org/10.1111/1467-9868.00353.

Stoddard, S.T., Forshey, B.M., Morrison, A.C., Paz-Soldan, V.A., Vazquez-Prokopec, G.M., Astete, H., Reiner, R.C., Vilcarromero, S., Elder, J.P., Halsey, E.S., Kochel, T.J., Kitron, U., Scott, T.W., 2013. House-to-house human movement drives dengue virus transmission. Proc. Natl. Acad. Sci. USA 110 (3), 994–999. https://doi.org/10.1073/pnas.1213349110. http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3549073/.

Strickman, D., Sithiprasasna, R., Kittayapong, P., Innis, B.L., 2000. Distribution of dengue and Japanese encephalitis among children in rural and suburban Thai villages. Am. J. Trop. Med. Hyg. 63 (1–2), 27–35.

Team, R.D.C., 2015. R: A Language and Environment for Statistical Computing. http://www.R-project.org.

ten Bosch, Q.A., Singh, B.K., Hassan, M.R.A., Chadee, D.D., Michael, E., 2016. The role of serotype interactions and seasonality in dengue model selection and control: insights from a pattern matching approach. PLOS Negl. Trop. Dis. 10 (5), e0004680. https://doi.org/10.1371/journal.pntd.0004680.

ten Bosch, Q.A., Clapham, H.E., Lambrechts, L., Duong, V., Buchy, P., Althouse, B.M., Lloyd, A.L., Waller, L.A., Morrison, A.C., Kitron, U., Vazquez-Prokopec, G.M., Scott, T.W., Perkins, T.A., 2018. Contributions from the silent majority dominate dengue virus transmission. PLoS Pathogens 14 (5). https://doi.org/10.1371/journal.ppat.

1006965. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5933708/.

Teurlai, M., Huy, R., Cazelles, B., Duboz, R., Baehr, C., Vong, S., 2012. Can human movements explain heterogeneous propagation of dengue fever in Cambodia? PLOS Negl. Trop. Dis. 6 (12), e1957. https://doi.org/10.1371/journal.pntd.0001957.

Thai, K.T., Binh, T.Q., Giao, P.T., Phuong, H.L., Hung, L.Q., Nam, N.V., Nga, T.T., Groen, J., Nagelkerke, N., Vries, P.J., 2005. Seroprevalence of dengue antibodies, annual incidence and risk factors among children in southern Vietnam. Trop. Med. Int. Health 10 (4), 379–386. https://doi.org/10.1111/j.1365-3156.2005.01388.x/full.

Thai, K.T.D., Nga, T.T.T., Van Nam, N., Phuong, H.L., Giao, P.T., Hung, L.Q., Binh, T.Q., Van Doornum, G.J.J., De Vries, P.J., 2007. Incidence of primary dengue virus infections in Southern Vietnamese children and reactivity against other flaviviruses. Trop. Med. Int. Health 12 (12), 1553–1557. https://doi.org/10.1111/j.1365-3156.2007.01964.x.

Tien, N.T.K., Luxemburger, C., Toan, N.T., Pollissard-Gadroy, L., Huong, V.T.Q., Van Be, P., Rang, N.N., Wartel, T.-A., Lang, J., 2010. A prospective cohort study of dengue infection in schoolchildren in Long Xuyen, Viet Nam. Trans. R. Soc. Trop. Med. Hyg. 104 (9), 592–600. https://doi.org/10.1016/j.trstmh.2010.06.003.

Vong, S., Khieu, V., Glass, O., Ly, S., Duong, V., Huy, R., Ngan, C., Wichmann, O., Letson, G.W., Margolis, H.S., Buchy, P., 2010. Dengue incidence in urban and rural Cambodia: results from population-based active fever surveillance, 2006–2008. PLoS Negl. Trop. Dis. 4 (11), e903. https://doi.org/10.1371/journal.pntd.0000903.

Vongpun sawad, S., Intharasongkroh, D., Thongmee, T., Poovorawan, Y., 2017. Seroprevalence of antibodies to dengue and chikungunya viruses in Thailand. PLOS ONE 12 (6), e0180560. https://doi.org/10.1371/journal.pone.0180560.

WHO, 2017. Dengue – World Health Organisation. http://www.who.int/denguecontrol/en/.

Wichmann, O., Yoon, I.-K., Vong, S., Limkittikul, K., Gibbons, R.V., Mammen, M.P., Ly, S., Buchy, P., Sirivichayakul, C., Buathong, R., Huy, R., Letson, G.W., Sabchareon, A., 2011. Dengue in Thailand and Cambodia: an assessment of the degree of unrecognized disease burden based on reported cases. PLoS Negl. Trop. Dis. 5 (3), e996. https://doi.org/10.1371/journal.pntd.0000996.

Wickham, H., 2009. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. http://ggplot2.org.

Yoon, I.-K., Rothman, A.L., Tannitisupawong, D., Srikiatkhachorn, A., Jarman, R.G., Aldstadt, J., Nisalak, A., Mammen, M.P., Thammapalo, S., Green, S., Libraty, D.H., Gibbons, R.V., Getis, A., Endy, T., Jones, J.W., Koenraadt, C.J.M., Morrison, A.C., Fansiri, T., Pimgate, C., Scott, T.W., 2012. Underrecognized mildly symptomatic viremic dengue virus infections in rural Thai schools and villages. J. Infect. Dis. 206 (3), 389–398. https://doi.org/10.1093/infdis/jis357. http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3490697/.