# Structure in the variability of the basic reproductive number (R 0 ) for Zika epidemics in the Pacific islands

Clara Champagne, David Georges Salthouse, Richard Paul, Van-Mai

Cao-Lormeau, Benjamin Roche, Bernard Cazelles

► To cite this version:

Clara Champagne, David Georges Salthouse, Richard Paul, Van-Mai Cao-Lormeau, Benjamin Roche, et al.. Structure in the variability of the basic reproductive number (R 0 ) for Zika epidemics in the Pacific islands. eLife, 2016, 5, pp.e19874. 10.7554/eLife.19874.001 . hal-01487852

HAL Id: hal-01487852

https://hal.sorbonne-universite.fr/hal-01487852

Submitted on 13 Mar 2017

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Structure in the variability of the basic reproductive number ($R_0$) for Zika epidemics in the Pacific islands

Clara Champagne¹,²*, David Georges Salthouse¹, Richard Paul³,⁴, Van-Mai Cao-Lormeau⁵, Benjamin Roche⁶, Bernard Cazelles¹,⁶*

¹IBENS, UMR 8197 CNRS-ENS Ecole Normale Supérieure, Paris, France; ²CREST, ENSAE, Université Paris Saclay, France; ³Department of Genomes and Genetics, Institut Pasteur, Unité de Génétique Fonctionnelle des Maladies Infectieuses, Paris, France; ⁴Centre National de la Recherche Scientifique, URA 3012, Paris, France;

⁵Unit of Emerging Infectious Diseases, Institut Louis Malardé, Tahiti, France;

⁶International Center for Mathematical and Computational Modeling of Complex Systems (UMMISCO), UPMC/IRD, Bondy cedex, France

**Abstract** Before the outbreak that reached the Americas in 2015, Zika virus (ZIKV) circulated in Asia and the Pacific: these past epidemics can be highly informative on the key parameters driving virus transmission, such as the basic reproduction number ($R_0$). We compare two compartmental models with different mosquito representations, using surveillance and seroprevalence data for several ZIKV outbreaks in Pacific islands (Yap, Micronesia 2007, Tahiti and Moorea, French Polynesia 2013-2014, New Caledonia 2014). Models are estimated in a stochastic framework with recent Bayesian techniques. $R_0$ for the Pacific ZIKV epidemics is estimated between 1.5 and 4.1, the smallest islands displaying higher and more variable values. This relatively low range of $R_0$ suggests that intervention strategies developed for other flaviviruses should enable as, if not more effective control of ZIKV. Our study also highlights the importance of seroprevalence data for precise quantitative analysis of pathogen propagation, to design prevention and control strategies.

DOI: 10.7554/eLife.19874.001

## Introduction

In May 2015, the first local cases of Zika were recorded in Brazil and by December of the same year the number of cases had surpassed 1.5 million. On February 2016, the World Health Organization declared Zika as a public health emergency of international concern (Who, 2016) and in March 2016, local transmission of Zika was recognized in 34 countries. Previously the Zika virus had circulated in Africa and Asia but only sporadic human cases had been reported. In 2007 the outbreak on Yap (Micronesia) was the first Zika outbreak outside Africa and Asia (Duffy et al., 2009). Since, Zika outbreaks have been also reported in French Polynesia and in New Caledonia (Cao-Lormeau et al., 2014; Dupont-Rouzeyrol et al., 2015) between 2013 and 2014 and subsequently, there have been cases of Zika disease in the Cook Islands, the Solomon Islands, Samoa, Vanuatu, and Easter Island (Chile) (see Figure 1 in Petersen et al. [2016]).

Zika virus (ZIKV) is a flavivirus, mostly transmitted via the bites of infected Aedes mosquitoes, although non-vector-borne transmission has been documented (sexual and maternofoetal transmission, laboratory contamination, transmission through transfusion) (Musso and Gubler, 2016). The most common clinical manifestations include rash, fever, arthralgia, and conjunctivitis (Musso and Gubler, 2016) but a large proportion of infections are asymptomatic or trigger mild symptoms that can remain unnoticed. Nevertheless, the virus may be involved in many severe

*For correspondence: champagn@biologie.ens.fr (CC); cazelles@biologie.ens.fr (BC)

**Competing interests:** The authors declare that no competing interests exist.

**Funding:** See page 30

**Received:** 02 August 2016

**Accepted:** 22 November 2016

**Published:** 29 November 2016

**Reviewing editor:** Mark Jit, London School of Hygiene & Tropical Medicine, and Public Health England, United Kingdom

© Copyright Champagne et al. This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
**eLife digest** Zika virus is an infectious disease primarily transmitted between people by mosquitoes. While most people develop mild flu-like symptoms, infection during pregnancy can interfere with how the baby's head and brain develop. Until recently, the virus had only been seen sporadically in Africa and Asia, but since 2007, outbreaks have been recorded on several Pacific islands. In 2015, the Zika virus reached the Americas, and within six months over 1.5 million cases had been reported in Brazil alone.

There is an urgent need to understand how the Zika virus moves within a population in order to help policymakers, and public health professionals, plan treatment and control of outbreaks of the disease. Researchers often use predictive models to estimate how a disease will spread. A parameter commonly calculated by these models is the "basic reproductive number", or $R_0$, which represents the average number of additional cases of the disease caused by one infected individual.

Using models that incorporated data from Zika virus outbreaks that occurred on several Pacific islands, Champagne et al. have produced estimates of $R_0$ that range from 1.5-4.1.

The $R_0$ values are greater than one, indicating that infection will spread within a population, but in the same range as those obtained for dengue fever, another closely related mosquito-borne disease. This suggests that by taking appropriate measures, the spread of Zika and dengue can be controlled to similar extents.

A closer look at the relationship between the population size and the predicted $R_0$ value for each Pacific island revealed an unexpected inverse relationship: the smaller the population, the larger the value of $R_0$. Since other regional factors may also explain these large differences between settings, further work is needed to disentangle context-specific from disease-specific factors. In this respect, data about seroprevalence (the number of people whose blood shows evidence of a past infection) in different populations is crucial for precisely analyzing the spread of Zika virus.

DOI: 10.7554/eLife.19874.002

neurological complications, including Guillain-Barre syndrome (Cao-Lormeau et al., 2016) and microcephaly in newborns (Schuler-Faccini et al., 2015). These complications and the impressive speed of its geographically propagation make the Zika pandemic a public health threat (Who, 2016). This reinforces the urgent need to characterize the different facets of virus transmission and to evaluate its dispersal capacity. We address this here by estimating the key parameters of ZIKV transmission, including the basic reproduction number ($R_0$), based on previous epidemics in the Pacific islands.

Defined as the average number of secondary cases caused by one typical infected individual in an entirely susceptible population, the basic reproduction number ($R_0$) is a central parameter in epidemiology used to quantify the magnitude of ongoing outbreaks and it provides insight when designing control interventions (Diekmann et al., 2010). It is nevertheless complex to estimate (Diekmann et al., 2010; van den Driessche and Watmough, 2002), and therefore, care must be taken when extrapolating the results obtained in a specific setting, using a specific mathematical model. In the present study, we explore the variability of $R_0$ using two models in several settings that had Zika epidemics in different years and that vary in population size (Yap, Micronesia 2007, Tahiti and Moorea, French Polynesia 2013–2014, and New Caledonia 2014). These three countries were successively affected by the virus, resulting in the first significant human outbreaks and they differ in several ways, including population size and location specific features. Hence, the comparison of their parameter estimates can be highly informative on the intrinsic variability of $R_0$. For each setting, we compare two compartmental models using different parameters defining the mosquito populations. Both models are considered in a stochastic framework, a necessary layer of complexity given the small population size and recent Bayesian inference techniques (Andrieu et al., 2010) are used for parameter estimation.

## Results

We use mathematical transmission models and data from surveillance systems and seroprevalence surveys for several ZIKV outbreaks in Pacific islands (Yap, Micronesia 2007 (Duffy et al., 2009), Tahiti
and Moorea, French Polynesia 2013–2014 (CHSP, 2014; Mallet et al., 2015; Aubry et al., 2015b), New Caledonia 2014 [DASS, 2014]) to quantify the ZIKV transmission variability.

Two compartmental models with vector-borne transmission are compared (cf. Figure 1). Both models use a Susceptible-Exposed-Infected-Resistant (SEIR) framework to describe the virus transmission in the human population, but differ in their representation of the mosquito population. Figure 1a is a schematic representation derived from Pandey et al. (2013) and formulates explicitly the mosquito population, with a Susceptible-Exposed-Infected (SEI) dynamic to account for the extrinsic incubation period (time taken for viral dissemination within the mosquito).

By contrast, in the second model (Figure 1b) based on Laneri et al. (2010) the vector is modeled implicitly: the two compartments κ and λ do not represent the mosquito population but the force of infection for vector to human transmission. This force of infection passes through two successive stages in order to include the delay associated with the extrinsic incubation period: κ stands for this latent phase of the force of infection whereas λ corresponds directly to the rate at which susceptible humans become infected.

The basic reproduction number of the models ($R_0$) is calculated using the next Generation Matrix method (Diekmann et al., 2010):

$$R_0^{\text{Pandey}} = \sqrt{\frac{\beta_H \beta_V \tau}{\gamma \mu (\mu + \tau)}}$$

$$R_0^{\text{Laneri}} = \sqrt{\frac{\beta}{\gamma}}$$

In addition, we consider that only a fraction $\rho$ of the total population is involved in the epidemic, due to spatial heterogeneity, immuno-resistance, or cross-immunity. For both models we define $N = \rho \cdot H$ with H the total size of the population reported by census.

The dynamics of ZIKV transmission in these islands is highly influenced by several sources of uncertainties. In particular, the small population size (less than 7000 inhabitants in Yap) leads to high variability in transmission rates. Therefore all these models are simulated in a discrete stochastic framework (Poisson with stochastic rates [Bretó et al., 2009]), to take this phenomenon into account. Stochasticity requires specific inference techniques: thus estimations are performed with PMCMC algorithm (particle Markov Chain Monte Carlo [Andrieu et al., 2010]).

Using declared Zika cases from different settings, the two stochastic models (Figure 1) were fitted (Figures 2–3). These models allow us to describe the course of the observed number of cases and estimate the number of secondary cases generated, $R_0$. Our estimates of $R_0$ lie between 1.6 (1.5–1.7) and 3.2 (2.4–4.1) and vary notably with respect to settings and models (Figures 2–3 and Tables 1–2). Strikingly, Yap displays consistently higher values of $R_0$ in both models and in general, there is an inverse relationship between island size and both the value and variability of $R_0$. This phenomenon may be explained by the higher stochasticity and extinction probability associated with smaller populations and can also reflect the information contained in the available data. However, the two highly connected islands in French Polynesia, Tahiti and Moorea, display similar values despite their differing sizes.

Regarding model variability, $R_0$ estimates are always higher and coarser with the Pandey model than with the Laneri model (cf. Tables 1–2). The Pandey model has two additional estimated parameters (in particular, the mosquito lifespan), which can explain the higher variability of the output. It is worth noting that these parameters are very sensitive (see Materials and methods). The difference in $R_0$ may also be linked to the difference in the estimated initial number of infected individuals ($H_I(0)$), which is higher in the Laneri model than in the Pandey model. Because of the high proportion of asymptomatic cases (the ratio of asymptomatic:symptomatic is estimated to be 1:1.3, V.-M Cao-Lormeau personal communication), it is hard to determine which scenario is more realistic, the time between introduction of the disease into the island and the first reported symptomatic case being unknown in most settings.

For the durations of infectious and intrinsic incubation (in human) and extrinsic incubation (in mosquito) periods, the posterior density resembles the informative prior (cf. Figures 6–13), indicating the models' incapacity to identify properly these parameters without more informative data.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Graphical representation of compartmental models. 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. (a) Pandey model (Pandey et al., 2013). *H*<sub>S</sub> susceptible individuals; *H*<sub>E</sub> infected (not yet infectious) individuals; *H*<sub>I</sub> infectious individuals; *H*<sub>R</sub> recovered individuals; *σ* is the rate at which *H*<sub>E</sub>-individuals move to infectious class *H*<sub>I</sub>; infectious individuals (*H*<sub>I</sub>) then recover at rate *γ*; *V*<sub>S</sub> susceptible vectors; *V*<sub>E</sub> infected (not yet infectious) vectors; *V*<sub>I</sub> infectious vectors; *V* constant size of total mosquito population; *τ* is the rate at which *V*<sub>E</sub>-vectors move to infectious class *V*<sub>I</sub>; vectors die at rate *μ*. (b) Laneri model (Laneri et al., 2010). *H*<sub>S</sub> susceptible individuals; *H*<sub>E</sub> infected (not yet infectious) individuals; *H*<sub>I</sub> infectious individuals; *H*<sub>R</sub> recovered individuals; *σ* is the rate at which *H*<sub>E</sub>-individuals move to infectious class *H*<sub>I</sub>; infectious individuals (*H*<sub>I</sub>) then recover at rate *γ*; implicit vector-borne transmission is modeled with the compartments *κ* and *λ*; *λ* current force of infection; *κ* latent force of infection reflecting the exposed state for mosquitoes during the extrinsic incubation period; *τ* is the transition rate associated to the extrinsic incubation period.</figcaption></figure>

DOI: 10.7554/eLife.19874.003
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Results using the Pandey model. Posterior median number of observed Zika cases (solid line), 95% credible intervals (shaded blue area) and data points (black dots). First column: particle filter fit. Second column: Simulations from the posterior density. Third column: *R*₀ posterior distribution. (a) Yap. (b) Moorea. (c) Tahiti. (d) New Caledonia. The estimated seroprevalences at the end of the epidemic (with 95% credibility intervals) are: (a) 73% (CI95: 68–77, observed 73%); (b) 49% (CI95: 45–53, observed 49%); (c) 49% (CI95: 45–53, observed 49%); (d) 39% (CI95: 8–92). See **Figure 4**. DOI: 10.7554/eLife.19874.004</figcaption></figure>

Moreover, these parameters have a clear sensitivity (see Materials and methods) and precise field
measures are therefore crucial for reliable model predictions.

The fraction *ρ* of the population involved in the epidemic is well estimated when the seroprevalence is known (in Yap and French Polynesia). In these cases, the proportion of the population involved is slightly greater than the seroprevalence rate, indicating a very high infection rate among
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Results using the Laneri model. Posterior median number of observed Zika cases (solid line), 95% credible intervals (shaded blue area) and data points (black dots). First column: particle filter fit. Second column: Simulations from the posterior density. Third column: *R*₀ posterior distribution. (a) Yap. (b) Moorea. (c) Tahiti. (d) New Caledonia. The estimated seroprevalences at the end of the epidemic (with 95% credibility intervals) are: (a) 72% (CI95: 68–77, observed 73%); (b) 49% (CI95: 45–53, observed 49%); c) 49% (CI95: 45–53, observed 49%); d) 65% (CI95: 24–91). See **Figure 5**.
DOI: 10.7554/eLife.19874.005</figcaption></figure>

involved individuals. In New Caledonia, as no information on seroprevalence was available, the fraction of population involved displays very large confidence intervals (cf. Tables 1 and 2), indicating that the model did not manage to identify properly this parameter with the available data. In this case, this parameter is highly correlated to the observation rate *r* (cf Figures 17 and 21): *r* and *ρ* seem unidentifiable without precise information on seroprevalence.
<table><caption>**Table 1.** Parameter estimations for the Pandey model. Posterior median (95% credible intervals). All the posterior parameter distributions are presented in **Figures 6–9**.</caption><thead><tr><th>Pandey model</th><th></th><th>Yap</th><th>Moorea</th><th>Tahiti</th><th>New Caledonia</th></tr></thead><tbody><tr><td>Population size</td><td>H</td><td>6892</td><td>16,200</td><td>178,100</td><td>268,767</td></tr><tr><td>Basic reproduction number</td><td>R<sub>0</sub></td><td>3.2 (2.4–4.1)</td><td>2.6 (2.2–3.3)</td><td>2.4 (2.0–3.2)</td><td>2.0 (1.8–2.2)</td></tr><tr><td>Observation rate</td><td>r</td><td>0.024 (0.019–0.032)</td><td>0.058 (0.048–0.073)</td><td>0.060 (0.050–0.073)</td><td>0.024 (0.010–0.111)</td></tr><tr><td>Fraction of population involved</td><td>ρ</td><td>74% (69–81)</td><td>50% (48–54)</td><td>50% (46–54)</td><td>40% (9–96)</td></tr><tr><td>Initial number of infected individuals</td><td>H<sub>I</sub>(0)</td><td>2 (1–8)</td><td>5 (0–21)</td><td>329 (16–1047)</td><td>37 (1–386)</td></tr><tr><td>Infectious period in human (days)</td><td>γ<sup>-1</sup></td><td>5.2 (4.1–6.7)</td><td>5.2 (4.1–6.8)</td><td>5.2 (4.1–6.7)</td><td>5.5 (4.2–6.8)</td></tr><tr><td>Extrinsic incubation period in mosquito (days)</td><td>τ<sup>-1</sup></td><td>10.6 (8.7–12.5)</td><td>10.5 (8.6–12.4)</td><td>10.5 (8.6–12.6)</td><td>10.7 (8.9–12.5)</td></tr><tr><td>Mosquito lifespan (days)</td><td>μ<sup>-1</sup></td><td>15.6 (11.7–19.3)</td><td>15.3 (11.4–19.1)</td><td>15.1 (11.2–19.0)</td><td>15.4 (11.6–19.4)</td></tr></tbody></table>

DOI: 10.7554/eLife.19874.006

Discussion

The reproduction number $R_0$ is a key parameter in epidemiology that characterizes the epidemic dynamics and the initial spread of the pathogen at the start of an outbreak in a susceptible population. $R_0$ can be used to inform public health authorities on the level of risk posed by an infectious disease, vaccination strategy, and the potential effects of control interventions (Anderson and May, 1982). In the light of the potential public health crisis generated by the international propagation of ZIKV, characterization of the potential transmissibility of this pathogen is crucial for predicting epidemic size, rate of spread and efficacy of intervention.

Using data from both surveillance systems and seroprevalence surveys in four different geographical settings across the Pacific (Duffy et al., 2009; CHSP, 2014; Mallet et al., 2015; DASS, 2014; Aubry et al., 2015b), we have estimated the basic reproductive number $R_0$ (see Figures 2–3 and Tables 1–2). Our estimate of $R_0$ obtained by inference based on Particle MCMC (Andrieu et al., 2010) has values in the range 1.6 (1.5–1.7) – 3.2 (2.4–4.1). Our $R_0$ estimates vary notably across settings. Lower and finer $R_0$ values are found in larger islands. This phenomenon can at least in part be explained by large spatial heterogeneities and higher demographic stochasticity for islands with smaller populations, as well as the influence of stochasticity on biological and epidemiological processes linked to virus transmission. This phenomenon can also be specific to the selection of the studied islands or can reflect a highly clustered geographical pattern, the global incidence curve being the smoothed overview of a collection of more explosive small size outbreaks. However, it is notable that the two French Polynesian islands yield similar estimates of $R_0$ despite differing population sizes. Indeed, other important factors differ among French Polynesia, New Caledonia and Yap, such as the human genetic background and their immunological history linked to the circulation of others arboviruses. Moreover, whilst both New Caledonia and French Polynesia populations were infected by the same ZIKV lineage and transmitted by the same principle vector species, Aedes aegypti, the epidemic in Yap occurred

<table><caption>**Table 2.** Parameter estimations for the Laneri model. Posterior median (95% credible intervals). All the posterior parameter distributions are presented in **Figures 10–13**.</caption><thead><tr><th>Laneri model</th><th></th><th>Yap</th><th>Moorea</th><th>Tahiti</th><th>New Caledonia</th></tr></thead><tbody><tr><td>Population size</td><td>H</td><td>6892</td><td>16,200</td><td>178,100</td><td>268,767</td></tr><tr><td>Basic reproduction number</td><td>R<sub>0</sub></td><td>2.2 (1.9–2.6)</td><td>1.8 (1.6–2.0)</td><td>1.6 (1.5–1.7)</td><td>1.6 (1.5–1.7)</td></tr><tr><td>Observation rate</td><td>r</td><td>0.024 (0.019–0.033)</td><td>0.057 (0.047–0.07)</td><td>0.057 (0.049–0.069)</td><td>0.014 (0.010–0.037)</td></tr><tr><td>Fraction of population involved</td><td>ρ</td><td>73% (69–78)</td><td>51% (47–55)</td><td>54% (49–59)</td><td>71% (27–98)</td></tr><tr><td>Initial number of infected individuals</td><td>H<sub>I</sub>(0)</td><td>2 (1–10)</td><td>9 (1–28)</td><td>667 (22–1570)</td><td>82 (2–336)</td></tr><tr><td>Infectious period in human (days)</td><td>γ<sup>-1</sup></td><td>5.3 (4.1–6.6)</td><td>5.3 (4.1–6.7)</td><td>5.2 (4.1–6.7)</td><td>5.4 (4.1–6.8)</td></tr><tr><td>Extrinsic incubation period in mosquito (days)</td><td>τ<sup>-1</sup></td><td>10.7 (8.8–12.7)</td><td>10.6 (8.6–12.6)</td><td>10.5 (8.5–12.5)</td><td>10.8 (8.9–12.8)</td></tr></tbody></table>

DOI: 10.7554/eLife.19874.007
much earlier with a different ZIKV lineage (Wang et al., 2016) and vectored by a different mosquito species *Aedes hensilli* (Ledermann et al., 2014). In French Polynesia, the vector *Aedes polynesiensis* is also present and dominates in Moorea with higher densities than in Tahiti. Finally, different vector control measures have been conducted in the three countries.

To date, studies investigating Zika outbreaks in the Pacific have always estimated R₀ using a deterministic framework. Using a similar version of the Pandey model in French Polynesia, Kucharski et al. (Kucharski et al., 2016) estimated R₀ between 1.6 and 2.3 (after scaling to square root for comparison) for Tahiti and between 1.8 and 2.9 in Moorea. These estimates are slightly lower and less variable than ours. This difference can be explained firstly by the chosen priors on mosquito parameters and secondly because our model includes demographic stochasticity. Moreover, they predicted a seroprevalence rate at the end of the epidemic of 95-97%, far from the 49% measured. In Yap island, a study (Funk et al., 2016) used a very detailed deterministic mosquito model, and estimated an R₀ for Zika between 2.9 and 8. In this case, our lower and less variable estimates may come from the fact that our model is more parsimonious in the number of uncertain parameters, especially concerning the mosquito population. Finally, a third study (Nishiura et al., 2016a) relied on another method for R₀ calculation (based on the early exponential growth rate of the epidemic) in French Polynesia as a whole and in Yap. Again, the obtained parameters are lower than ours in French Polynesia and higher in Yap. The first estimations for South America using a similar methodology (Nishiura et al., 2016b; Towers et al., 2016; Gao et al., 2016) lead to similar R₀ values. In all these studies a deterministic framework is used excluding the possibility of accounting for the high variability of biological and epidemiological processes exacerbated by the small size of the population. In these three studies, like in ours, it is worth noting that little insight is obtained regarding mosquito parameters. Posterior distribution mimics the chosen prior (cf. Figures 6-13). Both the simulation of the epidemics and the estimated R₀ are highly sensitive to the choice of priors on mosquito parameters, for which precise field measures are rare.

In the absence of sufficient data, the modeling of mosquito-borne pathogen transmission is a difficult task due to non-linearity and non-stationarity of the involved processes (Cazelles and Hales, 2006). This work has then several limitations. First, our study is limited by the completeness and quality of the data, with regard to both incidence and seroprevalence, but, above all, by the scarcity of information available on mosquitoes. Incidence data is aggregated at the island scale and cannot disentangle the effects of geography and observation noise to explain bimodal curves observed in Yap and New Caledonia. Moreover, although all data came from national surveillance systems, we had very little information about the potential degree of under-reporting, especially due to the high proportion of mildly symptomatic cases, at a time when the dangerous complications associated with the virus were unknown. Moreover, some cases might have been misdiagnosed as other flaviviruses, due to similarity in clinical manifestations or cross-reactivity in assays. Seroprevalence data were gathered from small sample sizes and were also sensitive to cross reactivity in populations non naive to dengue. In addition, they were missing in New Caledonia, which leads to strong correlation between our estimation of the observation rate and the fraction of the population involved in the epidemic. Because of the high proportion of asymptomatic or mildly symptomatic cases, the magnitude of the outbreaks is difficult to evaluate without precise seroprevalence data (Metcalf et al., 2016) or detection of mild, asymptomatic or pre-symptomatic infections (Thompson et al., 2016). Considering vectors, no demographic data were available and this partly explains the large variability of our R₀ estimations. Secondly, incidence and seroprevalence data were difficult to reconcile; the use of incidence data led to higher infection rates than those observed in the seroprevalence data. This difficulty has been overcome by considering that only a fraction of the population (p) is involved in the epidemic and then our model manages to reproduce the observed seroprevalence rate. This exposed fraction could be the result of spatial heterogeneity and high clustering of cases and transmission, as observed for dengue. The small dispersal of the mosquito and the limited population inter-mingling likely leads to considerable spatial variation in the extent of exposure to the virus and pockets of refugia in Tahiti as elsewhere (Telle et al., 2016). For instance, previous dengue infection rates in French Polynesia display large spatial variations even within islands (Daudens et al., 2009). Finer scale incidence and seroprevalence data would be useful to explore this. Another explanation for higher predicted than observed infection rates could be due to interaction with other flaviviruses.
The Zika outbreak was concomitant with dengue outbreaks in French Polynesia (CHSP, 2014; Mallet et al., 2015) and New Caledonia (DASS, 2014). Examples of coinfection have been reported (Dupont-Rouzeyrol et al., 2015) but competition between these close pathogens may also have occurred. Finally, mathematical models with vectorial transmission may tend to estimate high attack rates, sometimes leading to a contradiction between observed incidence and observed seroprevalence. Assumptions on the proportionality between infected mosquitoes and the force of infection, as well as the density-dependence assumption in these models could be questioned. Indeed even if these assumptions are at the heart of the mathematical models of mosquito-borne pathogen transmission (Reiner et al., 2013; Smith et al., 2014), a recent review (Halstead, 2008) and recent experimental results (Bowman et al., 2014; Harrington et al., 2014) question these important points.

On a final note, the estimates of R₀ for ZIKV are similar to but generally on the lower side of estimates made for two other flaviviruses of medical importance, dengue and Yellow Fever viruses (Favier et al., 2006; Imai et al., 2015; Massad et al., 2003), even though caution is needed in the comparison of studies with differing models, methods and data sources. Interventions strategies developed for dengue should thus enable as, if not more effective control of ZIKV, with the caveat that ZIKV remains principally a mosquito-borne pathogen with little epidemiological significance of the sexual transmission route. Even though further work and data are needed to support this hypothesis (Brauer et al., 2016), two recent studies indicated that sexual transmission alone is not sufficient to sustain an epidemic (Gao et al., 2016; Towers et al., 2016).

In conclusion, using recent stochastic modeling methods, we have been able to determine estimates of R₀ for ZIKV with an unexpected relationship with population size. Further data from the current Zika epidemic in South America that is caused by the same lineage as French Polynesia lead to estimates in the same range of values (Nishiura et al., 2016b; Towers et al., 2016; Gao et al., 2016). Our study highlights the importance of gathering seroprevalence data, especially for a virus that often leads to an asymptomatic outcome and it would provide a key component for precise quantitative analysis of pathogen propagation to enable improved planning and implementation of prevention and control strategies.

# Materials and methods

## Data

During the 2007 outbreak that struck Yap, 108 suspected or confirmed Zika cases (16 per 1000 inhabitants) were reported by reviewing medical records and conducting prospective surveillance between April 1st and July 29th 2007 (Duffy et al., 2009). In French Polynesia, sentinel surveillance recorded more than 8700 suspected cases (32 per 1000 inhabitants) across the whole territory between October 2013 and April 2014 (CHSP, 2014; Mallet et al., 2015). In New Caledonia, the first Zika case was imported from French Polynesia on 2013 November 12th. Approximately 2500 cases (9 per 1000 inhabitants) were reported through surveillance between January (first autochthonous case) and August 2014 (DASS, 2014).

For Yap and French Polynesia, the post-epidemic seroprevalence was assessed. In Yap, a household survey was conducted after the epidemic, yielding an infection rate in the island of 73% (Duffy et al., 2009). In French Polynesia, three seroprevalence studies were conducted. The first one took place before the Zika outbreak, and concluded that most of the population was naive for Zika virus (Aubry et al., 2015a). The second seroprevalence survey was conducted between February and March 2014, at the end of the outbreak, and reported a seroprevalence rate around 49% (Aubry et al., 2015b). The third one concerned only schoolchildren in Tahiti and was therefore not included in the present study.

Demographic data on population size were based on censuses from Yap (Duffy et al., 2009), French Polynesia (Insee, 2012), and New Caledonia (Insee, 2014).

## Models and inference

### Model equations

Although the models are simulated in a stochastic framework, we present them with ordinary differential equations for clarity. The reactions involved in the stochastic models are the same as
those governed by the deterministic equations, but the simulation process differs through the use of discrete compartments. It is described in the next section.

The equations describing Pandey model are:

$$ 
\begin{aligned}
\frac{dH_S}{dt} &= -\beta_H v_I H_S \\
\frac{dH_E}{dt} &= \beta_H v_I H_S - \sigma H_E \\
\frac{dH_I}{dt} &= \sigma H_E - \gamma H_I \\
\frac{dH_R}{dt} &= \gamma H_I \\
\frac{dv_S}{dt} &= \mu - \frac{\beta_v H_I}{N} v_S - \mu v_S \\
\frac{dv_E}{dt} &= \frac{\beta_v H_I}{N} v_S - \tau v_E - \mu v_E \\
\frac{dv_I}{dt} &= \tau v_E - \mu v_I
\end{aligned}
$$

where $v_S = \frac{v_S}{V}$ is the proportion of susceptible mosquitoes, $v_E = \frac{v_E}{V}$ the proportion of exposed mosquitoes, and $v_I = \frac{v_I}{V}$ the proportion of infected mosquitoes. Since we are using a discrete model, we cannot use directly the proportions $v_S$, $v_E$ and $v_I$ whose values are smaller than one. Therefore, we rescale using $V = H$, which leads to $v'_S = v_S \cdot H$, $v'_E = v_E \cdot H$, and $v'_I = v_I \cdot H$. In this model, the force of infection for humans is $\lambda_H = \beta_H v_I$, and the force of infection for mosquitoes is $\lambda_V = \beta_V \frac{H_I}{N}$.

The equations describing Laneri model are:

$$ 
\begin{aligned}
\frac{dH_S}{dt} &= -\lambda H_S \\
\frac{dH_E}{dt} &= \lambda H_S - \sigma H_E \\
\frac{dH_I}{dt} &= \sigma H_E - \gamma H_I \\
\frac{dH_R}{dt} &= \gamma H_I \\
\frac{d\kappa}{dt} &= \frac{2\beta_H v_I \tau}{N} - 2\tau\kappa \\
\frac{d\lambda}{dt} &= 2\tau\kappa - 2\tau\lambda
\end{aligned}
$$

In this model, the role of mosquitoes in transmission is represented only through the delay they introduce during the extrinsic incubation period (EIP, incubation period in the mosquito). For modeling reasons, this delay is included by representing the force of infection from infected humans to susceptible humans with two compartments $\kappa$ and $\lambda$: in this formalism, the duration between the moment when an exposed individual becomes infectious and the moment when another susceptible individual acquires the infection has a gamma distribution of mean $\tau^{-1}$ (Laneri et al., 2010; Roy et al., 2013; Lloyd, 2001). Therefore, $\lambda$ represents the current force of infection for humans $\lambda_H = \lambda$. The compartment $\kappa$ represents the same force of infection but at a previous stage, reflecting the exposed phase for mosquitoes during the extrinsic incubation period. As an analogy to Pandey model, the force of infection for mosquitoes is $\lambda_V = \frac{\beta_H v_I \tau}{v_S N}$, and therefore, the parameter $\beta$ can be interpreted as the product of a transmission parameter $\beta'$ by the proportion of susceptible mosquitoes: $\beta = v_S \beta'$. The force of infection for mosquitoes is then similar to Pandey's: $\lambda_V = \beta' \tau \frac{H_I}{N}$.

Again, since we are using a discrete model, we cannot use directly the proportions $\lambda$ and $\kappa$ whose values are smaller than one. Therefore, we rescale up to a factor $N$, which leads to $L = \lambda N$ and $K = \kappa N$.

In both models, following the dominant paradigm that diseases transmitted by Aedes mosquitoes are highly clustered, we restricted the total population $H$ measured by census to a fraction $N = \rho \cdot H$, in which the parameter $\rho$ is estimated. This formulation makes the hypothesis that a fraction of the total population is not at risk from the epidemic, because of individual factors or because the individuals remain in areas where the virus is not present. Moreover as the vector's flying range is small, the clustering of ZIKV infection may be reinforced. This may result in escapees from infection within the population, even at a single island scale. The available data does not allow further exploration of
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Infected and recovered humans evolution during the outbreak with Pandey model. Simulations from the posterior density: posterior median (solid line), 95% and 50% credible intervals (shaded blue areas) and observed seroprevalence (black dots). First column: Infected humans ($H_i$). Second column: Recovered humans ($H_R$). (a) Yap. (b) Moorea. (c) Tahiti. (d) New Caledonia.</figcaption></figure>

DOI: 10.7554/eLife.19874.013

the mechanisms underlying these phenomena, which seem fundamental to understand ZIKV propa-
gation. At the very least, the restriction to a fraction $\rho$ enables the model to reproduce the observed
seroprevalence rates, and to provide coherent results with respect to both data sources (seropreva-
lence and surveillance data).
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Infected and recovered humans evolution during the outbreak with Laneri model. Simulations from the posterior density: posterior median (solid line), 95% and 50% credible intervals (shaded blue areas) and observed seroprevalence (black dots). First column: Infected humans ($H_i$). Second column: Recovered humans ($H_R$). (a) Yap. (b) Moorea. (c) Tahiti. (d) New Caledonia.</figcaption></figure>

DOI: 10.7554/eLife.19874.014

## Stochastic framework

Both models are simulated in a stochastic and discrete framework, the Poisson with stochastic rates formulation (Bretó et al., 2009), to include the uncertainties related to small population size. In this framework, the number of reactions occurring in a time interval $dt$ is approximated by a multinomial distribution. In a model with $m$ possible reactions and $c$ compartments, $z_i$ being the state of the
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Posterior distributions. Pandey model, Yap island.
DOI: 10.7554/eLife.19874.015</figcaption></figure>

system at time *t* and *θ* the model parameters, the probability that each reaction with rate *r*<sup>*k*</sup> occurs *n*<sub>*k*</sub> times in *dt* is calculated as follows (Dureau *et al.*, 2013):

$$p(n_1, \dots, n_m | z_t, \theta) = \prod_{i=1}^{c} \left\{ M_i \left(1 - \sum_{X(k)=i} p_k\right)^{n_i} \prod_{X(k)=i} (p_k)^{n_k} \right\} + o(dt)$$
with, $z_t^{(i)}$ being the number of individual in compartment i at time t,

• $p_k = (1 - \exp\{-\sum_{X(l)=i} r^l(z_t, \theta) z_t^{X(l)} dt\}) \frac{r^k(z_t, \theta)}{\sum_{X(l)=i} r^l(z_t, \theta)}$

• $\bar{n}_i = z_t^{(i)} - \sum_{X(k)=i} n_k$ (the number of individuals staying in compartment i in $dt$)

• $M_i = (z_t^{(i)} - \bar{n}_i)$ (multinomial coefficient)

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Posterior distributions. Pandey model, Moorea island.

DOI: 10.7554/eLife.19874.016</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Posterior distributions. Pandey model, Tahiti island.</figcaption></figure>

DOI: 10.7554/eLife.19874.017

## Observation models

The only observed compartments are the infected humans (incidence measured every week) and the recovered humans (seroprevalence at the end of the outbreak when data is available). In order to link the model to the data, two observation models, for both incidence and seroprevalence data, are needed.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Posterior distributions. Pandey model, New Caledonia.
DOI: 10.7554/eLife.19874.018</figcaption></figure>

## Observation model on incidence data

The observed weekly incidence is assumed to follow a negative binomial distribution (Bretó et al., 2009) whose mean equals the number of new cases predicted by the model times an estimated observation rate *r*.

The observation rate *r* accounts for non observed cases, due to non reporting from medical centers, mild symptoms unseen by health system, and asymptomatic infections. Without additional data, it is not possible to make a distinction between these three categories of cases. We also
implicitely make the assumption that these cases transmit the disease as much as reported symptomatic cases.

The observation model for incidence data is therefore :

$$Inc_{obs} = \text{NegBin}(\phi^{-1}, \frac{1}{1 + \phi r Inc})$$

$Inc_{obs}$ being the observed incidence, and $Inc$ the incidence predicted by the model. The dispersion parameter (Bretó et al., 2009) $\phi$ is fixed at 0.1.

## Observation model on seroprevalence data

Seroprevalence data is fitted for Tahiti, Moorea, and Yap settings. It is assumed that the observed seroprevalence at the end of the epidemic follows a normal distribution with fixed standard deviation, whose mean equals the number of individuals in the $H_R$ compartment predicted by the model.

The observation model for seroprevalence data is therefore :

$$H_R^{obs} = \text{Normal}(H_R, \Lambda)$$

at the last time step, with notations detailed for each model in Table 3.

## Prior distributions

Informative prior distributions are assumed for the mosquito lifespan, the duration of infectious period, and both intrinsic and extrinsic incubation periods. The initial numbers of infected mosquitoes and humans are estimated, and the initial number of exposed individuals is set to the initial number of infected to reduce parameter space. We assume that involved populations are naive to Zika virus prior to the epidemic and set the initial number of recovered humans to zero. The other priors and associated references are listed in Table 4.

The range for the prior on observation rate is reduced for Tahiti and New Caledonia, in order to reduce the parameter space and facilitate convergence. In both cases, we use the information provided with the data source. In French Polynesia, 8750 cases we reported, but according to local health authorities, more than 32,000 people would have attended health facilities for Zika (Mallet et al., 2015) (8750/32000 ≤ 0.3). In New Caledonia, approximately 2500 cases were reported but more than 11,000 cases were estimated by health authorities (DASS, 2014) (2500/ 11000 ≤ 0.23). In both cases, these extrapolations are lower bounds on the real number of cases (in particular, they do not estimate the number of asymptomatic infections), and therefore can be used as upper bounds on the observation rate.

## Estimations

### Inference with PMCMC

The complete model is represented using the state space framework, with two equation systems: the transition equations refer to the transmission models, and the measurement equations are given by the observation models.

In a deterministic framework, this model could be directly estimated using MCMC, with a Metropolis-Hastings algorithm targeting the posterior distribution of the parameters. This algorithm would require the calculation of the model likelihood at each iteration.

<table><caption>**Table 3.** Details of the observation models for seroprevalence</caption><thead><tr><th>Island</th><th>Date</th><th>Standard deviation</th><th>Observed seroprevalence</th></tr><tr><td></td><td></td><td>Λ</td><td>H<sub>R</sub><sup>obs</sup></td></tr></thead><tbody><tr><td>Yap</td><td>2007-07-29</td><td>150</td><td>5005 (Duffy et al., 2009)</td></tr><tr><td>Moorea</td><td>2014-03-28</td><td>325</td><td>0.49 × 16200 = 7938 (Aubry et al., 2015b)</td></tr><tr><td>Tahiti</td><td>2014-03-28</td><td>3562</td><td>0.49 × 178100 = 87269 (Aubry et al., 2015b)</td></tr></tbody></table>

DOI: 10.7554/eLife.19874.008
<table><caption>**Table 4.** Prior distributions of parameters. 'Uniform[0,20]' indicates a uniform distribution in the range [0,20]. 'Normal(5,1) in [4,7]' indicates a normal distribution with mean five and standard deviation 1, restricted to the range [4,7].</caption><thead><tr><th>Parameters</th><th>Pandey model</th><th>Laneri model</th><th>References</th></tr></thead><tbody><tr><td>R<sup>2</sup><sub>0</sub></td><td>Uniform[0, 20]</td><td>Uniform[0, 20]</td><td>assumed</td></tr><tr><td>β<sub>V</sub></td><td>Uniform[0,10]</td><td>.</td><td>assumed</td></tr><tr><td>γ<sup>-1</sup></td><td>Normal(5,1) in [4,7]</td><td>Normal(5,1) in [4,7]</td><td>(Mallet et al., 2015)</td></tr><tr><td>σ<sup>-1</sup></td><td>Normal(4,1) in [2,7]</td><td>Normal(4,1) in [2,7]</td><td>(Nishiura et al., 2016b; Bearcroft, 1956; Lessler et al., 2016)</td></tr><tr><td>τ<sup>-1</sup></td><td>Normal(10.5,1) in [4,20]</td><td>Normal(10.5,1) in [4,20]</td><td>(Hayes, 2009; Chouin-Carneiro et al., 2016)</td></tr><tr><td>μ<sup>-1</sup></td><td>Normal(15,2) in [4,30]</td><td>.</td><td>(Brady et al., 2013; Liu-Helmersson et al., 2014)</td></tr><tr><td>ρ</td><td>Uniform[0,1]</td><td>Uniform[0,1]</td><td></td></tr><tr><th>Initial conditions<br>(t=0)</th><th>Pandey model</th><th>Laneri model</th><td></td></tr><tr><td>H<sub>I</sub>(0)</td><td>infected humans</td><td>Uniform[10<sup>-6</sup>,1]N</td><td>Uniform[10<sup>-6</sup>,1]N</td></tr><tr><td>H<sub>E</sub>(0)</td><td>exposed humans</td><td>H<sub>I</sub>(0)</td><td>H<sub>I</sub>(0)</td></tr><tr><td>H<sub>R</sub>(0)</td><td>recovered humans</td><td>0</td><td>0</td></tr><tr><td></td><td>infected vectors</td><td>V<sub>I</sub>(0)=Uniform[10<sup>-6</sup>,1]N</td><td>L(0)=Uniform[10<sup>-6</sup>,1]N</td></tr><tr><td></td><td>exposed vectors</td><td>V<sub>E</sub>(0) = V<sub>I</sub>(0)</td><td>K(0)=L(0)</td></tr><tr><th>Local conditions</th><th>Yap</th><th>Moorea</th><th>Tahiti</th><th>New<br>Caledonia</th><th>References</th></tr><tr><td>r</td><td>observation rate</td><td>Uniform<br>[0,1]</td><td>Uniform<br>[0,1]</td><td>Uniform<br>[0,0.3]</td><td>Uniform<br>[0,0.23]</td><td>(Mallet et al., 2015; DASS, 2014)</td></tr><tr><td>H</td><td>population size</td><td>6,892</td><td>16,200</td><td>178,100</td><td>268,767</td><td>(Duffy et al., 2009; Insee, 2012, 2014)</td></tr></tbody></table>

DOI: 10.7554/eLife.19874.009

In our stochastic framework, the model output is given only through simulations and the likelihood is intractable. In consequence, estimations are performed with the PMCMC algorithm (particle Markov Chain Monte Carlo (Andrieu et al., 2010)), in the PMMH version (particle marginal Metropolis-Hastings). This algorithm uses the Metropolis-Hastings structure, but replaces the real likelihood by its estimation with Sequential Monte Carlo (SMC).

Algorithm 1 PMCMC (Andrieu et al., 2010) (PMMH version, as in SSM (Dureau et al., 2013))

In a model with *n* observations and *J* particles.

$q(\cdot|\theta^{(i)})$ is the transition kernel.

1: Initialize $\theta^{(0)}$.

2: Using SMC algorithm, compute $\hat{p}(y_{1:n}|\theta^{(0)})$ and sample $x_{0:n}^*$ from $\hat{p}(x_{0:n}|y_{1:n}, \theta^{(0)})$.

3: **for** *i* = 1...*N* **do**

4: &nbsp;&nbsp;&nbsp;&nbsp;Sample $\theta^*$ from $q(\cdot|\theta^{(i)})$

5: &nbsp;&nbsp;&nbsp;&nbsp;Using SMC algorithm, compute $L(\theta^*) = \hat{p}(y_{1:n}|\theta^*)$ and sample $x_{0:n}^*$ from $\hat{p}(x_{0:n}|y_{1:n}, \theta^*)$

6: &nbsp;&nbsp;&nbsp;&nbsp;Accept $\theta^*$ (et $x_{0:n}^*$) with probability $1 \wedge \frac{L(\theta^*)q(\theta^*)}{L(\theta^*)q(\theta^*)}$

7: &nbsp;&nbsp;&nbsp;&nbsp;If accepted, $\theta^{(i+1)} = \theta^*$ and $x_{0:n}^{(i+1)} = x_{0:n}^*$. Otherwise $\theta^{(i+1)} = \theta^{(i)}$ and $x_{0:n}^{(i+1)} = x_{0:n}^{(i)}$.

8: **end for**

SMC (Doucet et al., 2001) is a filtering method that enables to recover the latent variables and estimate the likelihood for a given set of parameters. The data is treated sequentially, by adding one more data point at each iteration. The initial distribution of the state variables is approximated by a sample particles, and from one iteration to the next, all the particles are projected according
to the dynamic given by the model. The particles receive a weight according to the quality of their prediction regarding the observations. Before the next iteration, all the particles are resampled using these weights, in order to eliminate low weight particles and concentrate the computational effort in high probability regions. Model likelihood is also computed sequentially at each iteration (Dureau et al., 2013; Doucet and Johansen, 2011).

<caption>Algorithm 2 SMC (Sequential Monte Carlo, as implemented in SSM [Dureau et al., 2013])</caption>

In a model with n observations and J particles.
L is the model likelihood p(y₁:T|θ). W<sub>k</sub><sup>(j)</sup> is the weight and x<sub>k</sub><sup>(j)</sup> is the state associated to particle j at iteration k.
1: Set L = 1, W<sub>0</sub><sup>(j)</sup> = 1/J.
2: Sample (x<sub>0</sub><sup>(j)</sup>)<sub>j=1:J</sub> from p(x|θ<sub>0</sub>).
3: **for** k = 0 : n − 1 **do**
4:     **for** j = 0 : jdo
5:         Sample (x<sub>k+1</sub><sup>(j)</sup>)<sub>j=1:J</sub> from p(x<sub>k+1</sub>|x<sub>k</sub>, θ)
6:         Set α<sup>(j)</sup> = p(y<sub>k+1</sub>|x<sub>k+1</sub><sup>(j)</sup>, θ)
7:     **end for**
8:     Set W<sub>k+1</sub><sup>(j)</sup> = <sup>g<sup>(j)</sup></sup>/<sub>∑<sub>i=1</sub> α<sup>(i)</sup></sub> and L = L<sub>1</sub> <sup>1</sup>/<sub>∑<sub>j</sub> α<sup>(j)</sup></sub>
9:     Resample (x<sub>0:k+1</sub><sup>(j)</sup>)<sub>j=1:J</sub> from W<sub>k+1</sub><sup>(j)</sup>
10: **end for**

A gaussian kernel $q(\cdot|\theta^{(i)})$ is used in the PMCMC algorithm, with mean $\theta^{(i)}$ and fixed variance $\Sigma^q$
(random walk Metropolis Hastings).

Initialization

PMCMC algorithm is very sensitive to initialization of both the parameter values $\theta^{(0)}$ and the covariance matrix $\Sigma^q$. Several steps of initialization are therefore used.

Firstly, parameter values are initialized by maximum likelihood through simplex algorithm on a
deterministic version of the model. We apply the simplex algorithm to a set of 1000 points sampled
in the prior distributions and we select the parameter set with the highest likelihood.

Secondly, in order to initialise the covariance matrix, an adaptative MCMC (Metropolis Hastings) framework is used (Roberts and Rosenthal, 2009; Dureau et al., 2013). It uses the empirical covariance of the chain $\Sigma^{(i)}$, and aims to calibrate the acceptance rate of the algorithm to an optimal value. The transition kernel is also mixed (with a probability $\alpha = 0.05$) with another gaussian using the identity matrix to improve mixing properties.

$$
q^A(\cdot|x^{(i)}) = \alpha N(x^{(i)}, \lambda \frac{2.38^2}{d} Id) + (1-\alpha) N(x^{(i)}, \lambda \frac{2.38^2}{d} \Sigma^{(i)})
$$

The parameter $\lambda$ is approximated by successive iterations using the empirical acceptance rate of
the chain.

$$
\lambda_{i+1} = \lambda_i a^i (AccRate_i - 0.234)
$$

The adaptative PMCMC algorithm itself may have poor mixing properties without initialization. A first estimation of the covariance matrix is computed using KMCMC algorithm (Dureau et al., 2013). In the KMCMC algorithm, the model is simulated with stochastic differential equations (intermediate between deterministic and Poisson with stochastic rates frameworks) and the SMC part of the adaptative PMCMC is replaced by the extended Kalman filter. When convergence is reached with KMCMC, then, adaptative PMCMC is used.

The PMCMC algorithm is finally applied on the output of the adaptative PMCMC, using 50,000 iterations and 10,000 particles. Calculations are performed with SSM software (Dureau et al., 2013) and R version 3.2.2.

$R_0$ Calculation

$R_0$ is calculated using the Next Generation Matrix approach (NGM) (19).
## R₀ Calculation in Pandey model

$$F = \begin{pmatrix} 0 & 0 & 0 & \beta_H \\ 0 & 0 & 0 & 0 \\ 0 & \beta_v & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix} \quad v = \begin{pmatrix} -\sigma & 0 & 0 & 0 \\ 0 & -\gamma & 0 & 0 \\ 0 & 0 & -(\mu + \tau) & 0 \\ 0 & 0 & 0 & -\mu \end{pmatrix}$$

Then we have,

$$V^{-1} = \begin{pmatrix} -1/\sigma & 0 & 0 & 0 \\ -1/\gamma & -1/\gamma & 0 & 0 \\ 0 & 0 & -1/(\mu + \tau) & 0 \\ 0 & 0 & -\tau/(\mu + \mu) & -1/\mu \end{pmatrix}$$

and

$$FV^{-1} = \begin{pmatrix} 0 & 0 & -\beta_H \tau / [\mu(\tau + \mu)] & -\beta_H / \mu \\ 0 & 0 & 0 & 0 \\ -\beta_v / \gamma & -\beta_v / \gamma & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}$$

We calculate the eigen values α of $-FV^{-1}$ :

$$\begin{vmatrix} -\alpha & 0 & \beta_H \tau / [\mu(\tau + \mu)] & \beta_H / \mu \\ 0 & -\alpha & 0 & 0 \\ \beta_v / \gamma & \beta_v / \gamma & -\alpha & 0 \\ 0 & 0 & 0 & -\alpha \end{vmatrix} = \alpha^2 \left( \alpha^2 - \frac{\beta_H \beta_v \tau}{\gamma \mu (\tau + \mu)} \right) = 0$$

Then $α = 0$ or $α = \pm \sqrt{\frac{\beta_H \beta_v \tau}{\gamma \mu (\tau + \mu)}}$ and the highest eigenvalue is $R_0 = \sqrt{\frac{\beta_H \beta_v \tau}{\gamma \mu (\tau + \mu)}}$.

This formula defines $R_0$ as "the number of secondary cases per generation" (Dietz, 1993): i.e $R_0$ can be written as the geometric mean $R_0 = \sqrt{R_0^v R_0^h}$, where $R_0^v$ is the number of infected mosquitoes after the introduction of one infected human in a naive population, and $R_0^h$ is the number of infected humans after the introduction of one infected mosquito in a naive population. With this definition, herd immunity is reached when $(1 - R_0^{-2})$ of the population is vaccinated (Dietz, 1993).

## R₀ Calculation in Laneri model

Following the analogy with Pandey model, we compute the spectral radius of the NGM for the Laneri model.

$$F = \begin{pmatrix} 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \\ 0 & \beta \tau & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix} \quad V = \begin{pmatrix} -\sigma & 0 & 0 & 0 \\ 0 & -\gamma & 0 & 0 \\ 0 & 0 & -\tau & 0 \\ 0 & 0 & \tau & -\tau \end{pmatrix}$$

Then we have,

$$V^{-1} = \begin{pmatrix} -1/\sigma & 0 & 0 & 0 \\ -1/\gamma & -1/\gamma & 0 & 0 \\ 0 & 0 & -1/\tau & 0 \\ 0 & 0 & -1/\tau & -1/\tau \end{pmatrix}$$

and

$$FV^{-1} = \begin{pmatrix} 0 & 0 & -1/\tau & -1/\tau \\ 0 & 0 & 0 & 0 \\ -\beta \tau / \gamma & -\beta \tau / \gamma & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}$$

We calculate the eigen values α of $-FV^{-1}$ :
$$ \left| \begin{array}{cccc} -\alpha & 0 & 1/\tau & 1/\tau \\ 0 & -\alpha & 0 & 0 \\ \beta\tau/\gamma & \beta\tau/\gamma & -\alpha & 0 \\ 0 & 0 & 0 & -\alpha \end{array} \right| = \alpha^2 \left( \alpha^2 - \frac{\beta\tau}{\gamma\tau} \right) = 0 $$

Then $\alpha = 0$ or $\alpha = \pm\sqrt{\frac{\beta}{\gamma}}$ and the highest eigenvalue is $\alpha_R = \sqrt{\frac{\beta}{\gamma}}$.

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Posterior distributions. Laneri model, Yap island.

DOI: 10.7554/eLife.19874.019</figcaption></figure>
<table><caption>**Table 5.** Square root of the number of secondary cases after the introduction of a single infected individual in a naive population. Median and 95% credible intervals of 1000 deterministic simulations using parameters from the posterior distribution.</caption><thead><tr><th></th><th>Pandey model</th><th>Laneri model</th></tr></thead><tbody><tr><td>Yap</td><td>3.1 (2.5–4.3)</td><td>2.2 (1.9–2.6)</td></tr><tr><td>Moorea</td><td>2.6 (2.2–3.3)</td><td>1.8 (1.6–2.0)</td></tr><tr><td>Tahiti</td><td>2.4 (2.0–3.2)</td><td>1.6 (1.5–1.7)</td></tr><tr><td>New Caledonia</td><td>2.0 (1.8–2.2)</td><td>1.6 (1.5–1.7)</td></tr></tbody></table>

DOI: 10.7554/eLife.19874.010

Because $\lambda$ and $\kappa$ can be seen as parameters rather than state variables, the interpretation of the spectral radius as the $R_0$ of the model is not straightforward. Therefore, we computed the $R_0$ of the model through simulations, by counting the number of secondary infections after the introduction of a single infected individual in a naive population. Since Laneri model is considered here as a vector model, the number of infected humans after the introduction of a single infected is considered as $R_0$. We simulated 1000 deterministic trajectories, using parameter values sampled in the posterior distributions for all parameters except initial conditions. With this method, the confidence intervals for number of infected humans ($R_0^2$) are similar to the ones of $\alpha_R^2$ estimated by the model. As a consequence, $R_0$ was approximated by the spectral radius of the NGM in our results with our stochastic framework (cf. Table 5).

As a robustness check, the same method was applied to Pandey model : the confidence intervals for the number of secondary cases in simulations are very similar to the ones of $R_0^2$ (cf. Table 5).

## Sensitivity analysis

In order to analyse the influence of parameter values on the model's outputs, a sensitivity analysis was performed, using LHS/PRCC technique (Blower and Dowlatabadi, 1994), on Tahiti example. Similar results were obtained for the other settings. Three criteria were retained as outputs for the analysis: the seroprevalence at the last time point, the intensity of the peak of the outbreak and the

<table><caption>**Table 6.** Sensitivity analysis in Pandey model. Tahiti island. 1000 parameter sets were sampled with latin hypercube sampling (LHS), using 'lhs' R package (Carnell, 2016). On each parameter set, the model was simulated deterministically in order to explore variability in parameters without interference with variations due to stochasticity. PRCC were computed using the 'sensitivity' R package (Pujol et al., 2016).</caption><thead><tr><th>Parameters</th><th>Distribution</th><th>Seroprevalence</th><th>Peak intensity</th><th>Peak date</th></tr></thead><tbody><tr><td colspan="5"><strong>Model parameters</strong></td></tr><tr><td>R<sup>2</sup><sub>0</sub></td><td>Uniform[0,20]</td><td>0.87</td><td>0.90</td><td>-0.55</td></tr><tr><td>β<sub>V</sub></td><td>Uniform[0,10]</td><td>-0.66</td><td>-0.73</td><td>0.35</td></tr><tr><td>γ<sup>-1</sup></td><td>Uniform[4,7]</td><td>-0.25</td><td>0.10</td><td>0.20</td></tr><tr><td>σ<sup>-1</sup></td><td>Uniform[2,7]</td><td>-0.03</td><td>-0.10</td><td>0.15</td></tr><tr><td>τ<sup>-1</sup></td><td>Uniform[4,20]</td><td>-0.05</td><td>-0.07</td><td>0.06</td></tr><tr><td>μ<sup>-1</sup></td><td>Uniform[4,30]</td><td>-0.56</td><td>-0.70</td><td>0.49</td></tr><tr><td colspan="5"><strong>Initial conditions</strong></td></tr><tr><td>H<sub>i</sub>(0)</td><td>Uniform[2.10<sup>-5</sup>,0.011]</td><td>0.05</td><td>-0.02</td><td>0.02</td></tr><tr><td>V<sub>i</sub>(0)</td><td>Uniform[10<sup>-4</sup>,0.034]</td><td>0.11</td><td>-0.00</td><td>-0.26</td></tr><tr><td colspan="5"><strong>Fraction involved and observation model</strong></td></tr><tr><td>ρ</td><td>Uniform[0.46,0.54]</td><td>0.47</td><td>0.15</td><td>-0.03</td></tr><tr><td>r</td><td>Uniform[0.048,0.072]</td><td>-0.04</td><td>0.03</td><td>0.05</td></tr></tbody></table>

DOI: 10.7554/eLife.19874.011
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Posterior distributions. Laneri model, Moorea island.
DOI: 10.7554/eLife.19874.020</figcaption></figure>

date of the peak. We used uniform distributions for all parameters, which are listed in Tables 6 and 7. For model parameters, we used the same range as for the prior distribution. For initial conditions, the observation rate *r* and the fraction involved in the epidemic *ρ*, we used the 95% confidence interval obtained by PMCMC, in order to avoid unrealistic scenarios.

For all criteria, the key parameters in both models are transmission parameters (*R*₀ and *β*_v). High values for *R*₀ are positively correlated with a large seroprevalence and a high and early peak. On the contrary, high values for the parameters introducing a delay in the model, the incubation periods in
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Posterior distributions. Laneri model, Tahiti island.
DOI: 10.7554/eLife.19874.021</figcaption></figure>

human ($\sigma^{-1}$) and in mosquito ($\tau^{-1}$), are associated with a lower and later peak, and have no significant effect on seroprevalence. Moreover, the simulations are clearly sensitive to the other model parameters, in particular the mosquito lifespan ($\mu^{-1}$) in Pandey model.

Concerning other parameters, the initial conditions have a noticeable effect on the date of the peak only. As expected, the fraction involved in the epidemic ($\rho$) influences the magnitude of the outbreak, by calibrating the proportion of people than can be infected, but it has no significant effect on the timing of the peak.
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Posterior distributions. Laneri model, New Caledonia.
DOI: 10.7554/eLife.19874.022</figcaption></figure>

## Complementary results

These complementary results include PMCMC results for both models in the four settings: the epidemic trajectories regarding the human compartments for infected and recovered individuals (Figures 4,5), the detailed posterior distributions for all parameters (Figures 6–13) and correlation plots for all models (Figures 14–21).
<table><caption>**Table 7.** Sensitivity analysis in Laneri model. Tahiti island. 1000 parameter sets were sampled with latin hypercube sampling (LHS), using 'lhs' R package (Carnell, 2016). On each parameter set, the model was simulated deterministically in order to explore variability in parameters without interference with variations due to stochasticity. PRCC were computed using the 'sensitivity' R package (Pujol et al., 2016).</caption><thead><tr><th>Parameters</th><th>Distribution</th><th>Seroprevalence</th><th>Peak intensity</th><th>Peak date</th></tr></thead><tbody><tr><td colspan="5"><strong>Model parameters</strong></td></tr><tr><td>R<sup>2</sup><sub>0</sub></td><td>Uniform[0,20]</td><td>0.62</td><td>0.93</td><td>-0.50</td></tr><tr><td>γ<sup>-1</sup></td><td>Uniform[4,7]</td><td>0.01</td><td>0.62</td><td>0.15</td></tr><tr><td>σ<sup>-1</sup></td><td>Uniform[2,7]</td><td>-0.03</td><td>-0.54</td><td>0.21</td></tr><tr><td>τ<sup>-1</sup></td><td>Uniform[4,20]</td><td>-0.03</td><td>-0.70</td><td>0.47</td></tr><tr><td colspan="5"><strong>Initial conditions</strong></td></tr><tr><td>H<sub>i</sub>(0)</td><td>Uniform[10<sup>-5</sup>,0.015]</td><td>0.05</td><td>0.02</td><td>-0.32</td></tr><tr><td>L(0)</td><td>Uniform[2.10<sup>-5</sup>,0.004]</td><td>0.05</td><td>0.00</td><td>-0.16</td></tr><tr><td colspan="5"><strong>Fraction involved and observation model</strong></td></tr><tr><td>ρ</td><td>Uniform[0.49,0.59]</td><td>0.80</td><td>0.34</td><td>0.02</td></tr><tr><td>r</td><td>Uniform[0.048,0.068]</td><td>-0.01</td><td>0.01</td><td>-0.02</td></tr></tbody></table>

DOI: 10.7554/eLife.19874.012

## Correlation between estimated parameters

The inference technique may fail to estimate some parameters due to identifiability issues. In particular, when two parameters are highly correlated to one another, the model manages to estimate the pair of parameters but not each one separately. The analysis of correlation between parameters' posterior distributions can reveal such cases. The following graphics display for each model the correlation coefficients between all pairs of parameters across the MCMC chain. For example, in models for New Caledonia, the observation rate and the fraction of the population involved in the epidemic

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Correlation plot of MCMC output. Pandey model, Yap island.
DOI: 10.7554/eLife.19874.023</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Correlation plot of MCMC output. Pandey model, Moorea island.
DOI: 10.7554/eLife.19874.024</figcaption></figure>

are strongly negatively correlated (Figures 17,21): the inference technique does not manage to estimate properly these two parameters, due to the lack of information on seroprevalence.

<figure><img src="image_18.png" /><figcaption>**Figure 16.** Correlation plot of MCMC output. Pandey model, Tahiti island.
DOI: 10.7554/eLife.19874.025</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 17.** Correlation plot of MCMC output. Pandey model, New Caledonia.
DOI: 10.7554/eLife.19874.026</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 18.** Correlation plot of MCMC output. Laneri model, Yap island.
DOI: 10.7554/eLife.19874.027</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure 19.** Correlation plot of MCMC output. Laneri model, Moorea island.
DOI: 10.7554/eLife.19874.028</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 20.** Correlation plot of MCMC output. Laneri model, Tahiti island.
DOI: 10.7554/eLife.19874.029</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Figure 21.** Correlation plot of MCMC output. Laneri model, New Caledonia.
DOI: 10.7554/eLife.19874.030</figcaption></figure>

## Code and source data files

The estimation tools are provided by the open source software SSM (Dureau *et al.*, 2013) (State Space Models, RRID:SCR_014647), available at https://github.com/JDureau/ssm and https://github.com/sballesteros/ssm-predict. The codes for the implementation of each model are provided as supplementary file.

## Acknowledgements

CC, DGS and BC are partially supported by the "Pépinière interdisciplinaire Eco-Evo-Devo" from the Centre National de la Recherche Scientifique (CNRS). The research leading to these results has also received funding from the European Commission Seventh Framework Program [FP7/2007–2013] for the DENFREE project under Grant Agreement 282 378. The funders played no role in the study design, data collection, analysis, or preparation of the manuscript.

## Additional information

### Funding

<table><thead><tr><td>Funder</td><td>Grant reference number</td><td>Author</td></tr></thead><tbody><tr><td>Centre National de la Recherche Scientifique</td><td>Pépinière interdisciplinaire Eco-Evo-Devo</td><td>Clara Champagne<br>David Georges Salthouse<br>Bernard Cazelles</td></tr><tr><td>European Commission</td><td>Seventh Framework Program [FP7 2007-2013] for the DENFREE project under Grant Agreement 282 348</td><td>Clara Champagne<br>David Georges Salthouse<br>Richard Paul<br>Van-Mai Cao-Lormeau<br>Bernard Cazelles</td></tr></tbody></table>

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
## Author contributions

CC, BC, Contributed to the study design, numerical part of the study, interpretation of the results and writing of the manuscript; DGS, Contributed to the numerical part of the study, interpretation of the results and writing of the manuscript; RP, V-MC-L, BR, Contributed to the interpretation of the results and writing of the manuscript

## Author ORCIDs

Clara Champagne, http://orcid.org/0000-0002-0369-6758

Bernard Cazelles, http://orcid.org/0000-0002-7972-361X

## Additional files

### Supplementary files

* Supplementary file 1. Codes for the implementation of each model.
  DOI: 10.7554/eLife.19874.031

* Supplementary file 2. Data file, Yap island. Number of cases per week and number of immune individuals at the end of the epidemic (Duffy et al., 2009).
  DOI: 10.7554/eLife.19874.032

* Supplementary file 3. Data file, Moorea island. Number of cases per week and number of immune individuals at the end of the epidemic (Mallet et al., 2015; Aubry et al., 2015b).
  DOI: 10.7554/eLife.19874.033

* Supplementary file 4. Data file, Tahiti island. Number of cases per week and number of immune individuals at the end of the epidemic (Mallet et al., 2015; Aubry et al., 2015b).
  DOI: 10.7554/eLife.19874.034

* Supplementary file 5. Data file, New Caledonia. Number of cases per week (DASS, 2014).
  DOI: 10.7554/eLife.19874.035

## References

* Anderson RM, May RM. 1982. Directly transmitted infections diseases: control by vaccination. *Science* **215**:1053-1060. doi: 10.1126/science.7063839, PMID: 7063839

* Andrieu C, Doucet A, Holenstein R. 2010. Particle Markov chain Monte Carlo methods. *Journal of the Royal Statistical Society: Series B* **72**:269-342. doi: 10.1111/j.1467-9868.2009.00736.x

* Aubry M, Finke J, Teissier A, Roche C, Broult J, Paulous S, Desprès P, Cao-Lormeau VM, Musso D. 2015a. Seroprevalence of arboviruses among blood donors in French Polynesia, 2011-2013. *International Journal of Infectious Diseases* **41**:11-12. doi: 10.1016/j.ijid.2015.10.005, PMID: 26482390

* Aubry M, Teissier A, Roche C, Teururai S, Desprès P, Mallet H, Merceron S, Huart M, Sicard S, Deparis X, Musso D, Cao-Lormeau V. 2015b. Serosurvey of Dengue, Zika and Other Mosquito-Borne Viruses in French Polynesia. *Poster 765*. Philadelphia, PA, USA: 64th Annual Meeting of the American Society of Tropical Medicine and Hygiene. Oct 25-29. http://www abstractsonline.com/Plan/ViewAbstract.aspx? sKey=2e5199c4-aceb-4d61-8769-342586917c5a& cKey=74ffe328-3a83-4f28-84c6-5656971501c8&mKey=%7bAB652FDF-0111-45C7-A5E5-0BA9D4AF5E12%7d

* Bearcroft WG. 1956. Zika virus infection experimentally induced in a human volunteer. *Transactions of the Royal Society of Tropical Medicine and Hygiene* **50**:442-448. doi: 10.1016/0035-9203(56)90090-6, PMID: 13380987

* Blower SM, Dowlatabadi H. 1994. Sensitivity and uncertainty analysis of complex models of disease transmission: An HIV Model, as an Example. *International Statistical Review / Revue Internationale de Statistique* **62**:229-243. doi: 10.2307/1403510

* Bowman LR, Runge-Ranzinger S, McCall PJ. 2014. Assessing the relationship between vector indices and dengue transmission: a systematic review of the evidence. *PLoS Neglected Tropical Diseases* **8**:e2848. doi: 10.1371/journal.pntd.0002848, PMID: 24810901

* Brady OJ, Johansson MA, Guerra CA, Bhatt S, Golding N, Pigott DM, Delatte H, Grech MG, Leishnam PT, Maciel-de-Freitas R, Styer LM, Smith DL, Scott TW, Gething PW, Hay SL. 2013. Modelling adult Aedes aegypti and Aedes albopictus survival at different temperatures in laboratory and field settings. *Parasites & Vectors* **6**: 351. doi: 10.1186/1756-3305-6-351, PMID: 24330720

* Brauer F, Castillo-Chavez C, Mubayi A, Towers S. 2016. Some models for epidemics of vector-transmitted diseases. *Infectious Disease Modelling*. doi: 10.1016/j.idm.2016.08.001

* Bretó C, He D, Ionides EL, King AA. 2009. Time series analysis via mechanistic models. *The Annals of Applied Statistics* **3**:319-348. doi: 10.1214/08-AOAS201
Cao-Lormeau V-M, Blake A, Mons S, Lastère S, Roche C, Vanhomwegen J, Dub T, Baudouin L, Teissier A, Larre P, Vial A-L, Decam C, Choumet V, Halstead SK, Willison HJ, Musset L, Manuguerra J-C, Despres P, Fournier E, Mallet H-P, et al. 2016. Guillain-Barré Syndrome outbreak associated with Zika Virus infection in French Polynesia: a case-control study. *The Lancet* **387**:1531–1539. doi: 10.1016/S0140-6736(16)00562-6

Cao-Lormeau V-M, Roche C, Teissier A, Robin E, Berry A-L, Mallet H-P, Sall AA, Musso D. 2014. Zika Virus, French Polynesia, South Pacific, 2013. *Emerging Infectious Diseases* **20**:1084–1086. doi: 10.3201/eid2006.140138

Carnell R. 2016. Package 'lhs'. 0.14. http://cran.stat.auckland.ac.nz/web/packages/lhs/lhs.pdf

Cazelles B, Hales S. 2006. Infectious diseases, climate influences, and nonstationarity. *PLoS Medicine* **3**:e328. doi: 10.1371/journal.pmed.0030328, PMID: 16903777

Chouin-Carneiro T, Vega-Rua A, Vazeille M, Yebakima A, Girod R, Goindin D, Dupont-Rouzeyrol M, Lourenço-de-Oliveira R, Failloux AB. 2016. Differential Susceptibilities of Aedes aegypti and Aedes albopictus from the Americas to Zika Virus. *PLoS Neglected Tropical Diseases* **10**:e0004543. doi: 10.1371/journal.pntd.0004543, PMID: 26938868

CHSP. 2014. Surveillance de la dengue et du zika en Polynésie française, données actualisées au 28 mars 2014. *Centre d'hygiène et de salubrité publique*. http://www.hygiene-publique.gov.pf/IMG/pdf/bulletin_dengue_28-03-14.pdf

DASS. 2014. Situation sanitaire en Nouvelle Calédonie. Les arboviroses : Dengue, Chikungunya, Zika. *Direction des Affaires Sanitaires et Sociales de Nouvelle Calédonie*. http://www.dass.gouv.nc/portal/page/portal/dass/librairie/fichiers/32142267.PDF

Daudens E, Lastere S, Hirschauer C, Cao-Lormeau V-M, Louette R, Roche C, Chee Ayee A, Goffard N, Vrousos E, Renou L, Wiegandt A. 2009. Épidémiologie de la dengue et stratégies de lutte en Polynésie française, 2006-2008. Numéro thématique. Polynésie française : une situation particulière. *Bulletin Épidémiologique Hebdomadaire*, N°48-49-50:499-503.

Diekmann O, Heesterbeek JAP, Roberts MG. 2010. The construction of next-generation matrices for compartmental epidemic models. *Journal of the Royal Society Interface* **7**:873–885. doi: 10.1098/rsif.2009.0386

Dietz K. 1993. The estimation of the basic reproduction number for infectious diseases. *Statistical Methods in Medical Research* **2**:23–41. doi: 10.1177/096228029300200103, PMID: 8261248

Doucet A, Freitas N, Gordon N. 2001. Sequential Monte Carlo Methods in Practice. New York: Springer. doi: 10.1007/978-1-4757-3437-9

Doucet A, Johansen AM. 2011. A tutorial on particle filtering and smoothing: fifteen years later.

Duffy MR, Chen TH, Hancock WT, Powers AM, Kool JL, Lanciotti RS, Pretrick M, Marfel M, Holzbauer S, Dubray C, Guillaumot L, Griggs A, Bel M, Lambert AJ, Laven J, Kosoy O, Panella A, Biggerstaff BJ, Fischer M, Hayes EB. 2009. Zika virus outbreak on yap island, Federated States of Micronesia. *New England Journal of Medicine* **360**:2536–2543. doi: 10.1056/NEJMoa0805715, PMID: 19516034

Dupont-Rouzeyrol M, O'Connor O, Calvez E, Daurès M, John M, Grangeon JP, Gourinat AC. 2015. Co-infection with Zika and dengue viruses in 2 patients, New Caledonia, 2014. *Emerging Infectious Diseases* **21**:381–382. doi: 10.3201/eid2102.141553, PMID: 25625687

Dureau J, Ballesteros S, Bogich T. 2013 SSM: Inference for time series analysis with State Space Models. https://github.com/JDureau/ssm/blob/master/doc/doc.pdf.

Favier C, Degallier N, Rosa-Freitas MG, Boulanger JP, Costa Lima JR, Luitgards-Moura JF, Menkès CE, Mondet B, Oliveira C, Weimann ET, Tsouris P. 2006. Early determination of the reproductive number for vector-borne diseases: the case of dengue in Brazil. *Tropical Medicine and International Health* **11**:332–340. doi: 10.1111/j.1365-3156.2006.01560.x, PMID: 16553913

Funk S, Kucharski AJ, Camacho A, Eggo RM, Yakob L, Edmunds WJ. 2016. Comparative analysis of dengue and Zika outbreaks reveals differences by setting and virus. *bioRxiv*. doi: 10.1101/043265

Gao D, Lou Y, He D, Porco TC, Kuang Y, Chowell G, Ruan S. 2016. Prevention and Control of Zika as a Mosquito-Borne and Sexually Transmitted Disease: A Mathematical Modeling Analysis. *Scientific Reports* **6**:28070. doi: 10.1038/srep28070, PMID: 27312324

Halstead SB. 2008. Dengue virus-mosquito interactions. *Annual Review of Entomology* **53**:273–291. doi: 10.1146/annurev.ento.53.103106.093326, PMID: 17803458

Harrington LC, Fleisher A, Ruiz-Moreno D, Vermeylen F, Wa CV, Poulson RL, Edman JD, Clark JM, Jones JW, Kitthawee S, Scott TW. 2014. Heterogeneous feeding patterns of the dengue vector, Aedes aegypti, on individual human hosts in rural Thailand. *PLoS Neglected Tropical Diseases* **8**:e3048. doi: 10.1371/journal.pntd.0003048, PMID: 25102306

Hayes EB. 2009. Zika Virus outside Africa. *Emerging Infectious Diseases* **15**:1347–1350. doi: 10.3201/eid1509.090442, PMID: 19788800

Imai N, Dorigatti I, Cauchemez S, Ferguson NM. 2015. Estimating dengue transmission intensity from sero-prevalence surveys in multiple countries. *PLoS Neglected Tropical Diseases* **9**:e0003719. doi: 10.1371/journal.pntd.0003719, PMID: 25881272

Insee. 2012. *Populations légales au recensement de la population 2012 De Polynésie française*: Institut National de la Statistique et des Etudes Economiques. https://www.insee.fr/fr/statistiques/2122859.

Insee. 2014. *Populations légales au recensement de la population 2014 de Nouvelle-Calédonie*: Institut National de la Statistique et des Etudes Economiques. https://www.insee.fr/fr/statistiques/2122859.

Kucharski AJ, Funk S, Eggo RM, Mallet HP, Edmunds WJ, Nilles EJ. 2016. Transmission Dynamics of Zika Virus in Island Populations: A Modelling Analysis of the 2013-14 French Polynesia Outbreak. *PLoS Neglected Tropical Diseases* **10**:e0004726. doi: 10.1371/journal.pntd.0004726, PMID: 27186984
Laneri K, Bhadra A, Ionides EL, Bouma M, Dhiman RC, Yadav RS, Pascual M. 2010. Forcing versus feedback: epidemic malaria and monsoon rains in northwest India. *PLoS Computational Biology* **6:e1000898. doi: 10.1371/journal.pcbi.1000898, PMID: 20824122**

Ledermann JP, Guillaumot L, Yug L, Saweyog SC, Tided M, Machieng P, Pretrick M, Marfel M, Griggs A, Bel M, Duffy MR, Hancock WT, Ho-Chen T, Powers AM. 2014. Aedes hensilli as a potential vector of Chikungunya and Zika Viruses. *PLoS Neglected Tropical Diseases* **8:e3188. doi: 10.1371/journal.pntd.0003188, PMID: 25299181**

Lessler J, Ott CT, Carcelen AC, Konikoff JM, Williamson J, Bi Q, Kucirka LM, Cummings DA, Reich NG, Chaisson LH. 2016. Times to key events in Zika Virus infection and implications for blood donation: a systematic review. *Bulletin of the World Health Organization* **94**:841–849. doi: 10.2471/BLT.16.174540, PMID: 27821887

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:e89783. doi: 10.1371/journal.pone.0089783, PMID: 24603439**

Lloyd AL. 2001. Destabilization of epidemic models with the inclusion of realistic distributions of infectious periods. *Proceedings of the Royal Society B: Biological Sciences* **268**:985–993. doi: 10.1098/rspb.2001.1599

Mallet H-P, Vial A, Musso D. 2015. Bilan de l'épidémie à virus Zika en Polynésie française. *Bulletin d'Informations Sanitaires Epidemiologiques Et Statistiques* **13** http://www.hygiene-publique.gov.pf/IMG/pdf/no13_-_mai_2015_-_zika.pdf.

Massad E, Burattini MN, Coutinho FA, Lopez LF. 2003. Dengue and the risk of urban yellow fever reintroduction in São Paulo State, Brazil. *Revista de Saúde Pública* **37**:477–484. doi: 10.1590/S0034-89102003000400013, PMID: 12937709

Metcalf CJE, Farrar J, Cutts FT, Basta NE, Graham AL, Lessler J, Ferguson NM, Burke DS, Grenfell BT. 2016. Use of serological surveys to generate key insights into the changing global landscape of infectious disease. *The Lancet* **388**:728–730. doi: 10.1016/S0140-6736(16)30164-7

Musso D, Gubler DJ. 2016. Zika virus. *Clinical Microbiology Reviews* **29**:487–524. doi: 10.1128/CMR.00072-15, PMID: 27029595

Nishiura H, Kinoshita R, Mizumoto K, Yasuda Y, Nah K. 2016a. Transmission potential of Zika Virus infection in the South Pacific. *International Journal of Infectious Diseases* **45**:95–97. doi: 10.1016/j.ijid.2016.02.017, PMID: 26923081

Nishiura H, Mizumoto K, Villamil-Gómez WE, Rodríguez-Morales AJ. 2016b. Preliminary estimation of the basic reproduction number of zika virus infection during Colombia epidemic, 2015-2016. *Travel Medicine and Infectious Disease* **14**:274–276. doi: 10.1016/j.tmaid.2016.03.016, PMID: 27060613

Pandey A, Mubayi A, Medlock J. 2013. Comparing vector-host and SIR models for dengue transmission. *Mathematical Biosciences* **246**:252–259. doi: 10.1016/j.mbs.2013.10.007, PMID: 24427785

Petersen LR, Jamieson DJ, Powers AM, Honein MA. 2016. Zika virus. *The New England Journal of Medicine* **374**: 1552–1563. doi: 10.1056/NEJMra1602113, PMID: 27028561

Pujol G, Iooss B, Janon A. 2016. Package 'sensitivity'. https://cran.r-project.org/web/packages/sensitivity/sensitivity.pdf

Reiner RC, Perkins TA, Barker CM, Niu T, Chaves LF, Ellis AM, George DB, Le Menach A, Pulliam JRC, Bisanzio D, Buckee C, Chiyaka C, Cummings DAT, Garcia AJ, Gatton ML, Gething PW, Hartley DM, Johnston G, Klein EY, Michael E, et al. 2013. A systematic review of mathematical models of mosquito-borne pathogen transmission: 1970-2010. *Journal of The Royal Society Interface* **10**:20120921 doi: 10.1098/rsif.2012.0921

Roberts GO, Rosenthal JS. 2009. Examples of adaptive MCMC. *Journal of Computational and Graphical Statistics* **18**:349–367. doi: 10.1198/jcgs.2009.06134

Roy M, Bouma MJ, Ionides EL, Dhiman RC, Pascual M. 2013. The potential elimination of Plasmodium vivax malaria by relapse treatment: insights from a transmission model and surveillance data from NW India. *PLoS Neglected Tropical Diseases* **7:e1979. doi: 10.1371/journal.pntd.0001979, PMID: 23326611**

Schuler-Faccini L, Ribeiro EM, Feitosa IML, Horovitz DDG, Cavalcanti DP, Pessoa A, Doriqui MJR, Neri JI, Neto Joao Monteiro de Pina, Wanderley HYC, Cernach M, El-Husny AS, Pone MVS, Serao CLC, Sanseverino MTV, Neto J, Neto JM, Brazilian Medical Genetics Society–Zika Embryopathy Task Force. 2015. Possible Association Between Zika Virus Infection and Microcephaly — Brazil, 2015. *MMWR. Morbidity and Mortality Weekly Report* **65**:59–62. doi: 10.15585/mmwr.mm6503e2

Smith DL, Perkins TA, Reiner RC, Barker CM, Niu T, Chaves LF, Ellis AM, George DB, Le Menach A, Pulliam JR, Bisanzio D, Buckee C, Chiyaka C, Cummings DA, Garcia AJ, Gatton ML, Gething PW, Hartley DM, Johnston G, Klein EY, et al. 2014. Recasting the theory of mosquito-borne pathogen transmission dynamics and control. *Transactions of the Royal Society of Tropical Medicine and Hygiene* **108**:185–197. doi: 10.1093/trstmh/tru026, PMID: 24591453

Telle O, Vaguet A, Yadav NK, Lefebvre B, Cebeillac A, Nagpal BN, Daudé E, Paul RE. 2016. The spread of dengue in an endemic urban milieu–The Case of Delhi, India. *PLoS One* **11:e0146539. doi: 10.1371/journal.pone.0146539, PMID: 26808518**

Thompson RN, Gilligan CA, Cunniffe NJ. 2016. Detecting presymptomatic infection is necessary to forecast major epidemics in the earliest stages of infectious disease outbreaks. *PLOS Computational Biology* **12:e1004836. doi: 10.1371/journal.pcbi.1004836, PMID: 27046030**

Towers S, Brauer F, Castillo-Chavez C, Falconar AK. I, Mubayi A, Romero-Vivas CME. 2016. Estimation of the reproduction number of the 2015 zika virus outbreak in Barranquilla, Colombia, and a first estimate of the relative role of sexual transmission. *arXiv*. https://arxiv.org/abs/1606.01422.
van den Driessche P, Watmough J. 2002. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. *Mathematical Biosciences* **180**:29-48. doi: 10.1016/S0025-5564 (02)00108-6, PMID: 12387915

Wang L, Valderramos SG, Wu A, Ouyang S, Li C, Brasil P, Bonaldo M, Coates T, Nielsen-Saines K, Jiang T, Aliyari R, Cheng G. 2016. From mosquitos to humans: Genetic Evolution of Zika Virus. *Cell Host & Microbe* **19**:561-565. doi: 10.1016/j.chom.2016.04.006, PMID: 27091703

**World Health Organisation**. 2016. WHO statement on the first meeting of the International Health Regulations (2005) (IHR 2005) Emergency Committee on Zika virus and observed increase in neurological disorders and neonatal malformations. http://www.who.int/mediacentre/news/statements/2016/1st-emergency-committee-zika/en/