Abstract Between October 2013 and April 2014, more than 30,000 cases of Zika virus (ZIKV) disease were estimated to have attended healthcare facilities in French Polynesia. ZIKV has also been reported in Africa and Asia, and in 2015 the virus spread to South America and the Caribbean. Infection with ZIKV has been associated with neurological complications including Guillain-Barré Syndrome (GBS) and microcephaly, which led the World Health Organization to declare a Public Health Emergency of International Concern in February 2015. To better understand the transmission dynamics of ZIKV, we used a mathematical model to examine the 2013–14 outbreak on the six major archipelagos of French Polynesia. Our median estimates for the basic reproduction number ranged from 2.6–4.8, with an estimated 11.5% (95% CI: 7.32–17.9%) of total infections reported. As a result, we estimated that 94% (95% CI: 91–97%) of the total population of the six archipelagos were infected during the outbreak. Based on the demography of French Polynesia, our results imply that if ZIKV infection provides complete protection against future infection, it would take 12–20 years before there are a sufficient number of susceptible individuals for ZIKV to re-emerge, which is on the same timescale as the circulation of dengue virus serotypes in the region. Our analysis suggests that ZIKV may exhibit similar dynamics to dengue virus in island populations, with transmission characterized by large, sporadic outbreaks with a high proportion of asymptomatic or unreported cases.

Author Summary Since the first reported major outbreak of Zika virus disease in Micronesia in 2007, the virus has caused outbreaks throughout the Pacific and South America. Transmitted by the Aedes genus of mosquitoes, the virus has been linked to possible neurological complications including Guillain-Barré Syndrome and microcephaly. To improve our understanding of the transmission dynamics of Zika virus in island populations, we analysed the 2013–14 outbreak on the six major archipelagos of French Polynesia. We found evidence that Zika virus infected the majority of the population, but only around 12% of total infections on the archipelagos were reported as cases. If infection with Zika virus generates lifelong immunity, we estimate that it would take at least 12–20 years before there are enough susceptible people for the virus to re-emerge. Our results suggest that Zika virus could exhibit similar dynamics to dengue virus in the Pacific, producing large but sporadic outbreaks in small island populations.

Citation: Kucharski AJ, Funk S, Eggo RM, Mallet H-P, Edmunds WJ, Nilles EJ (2016) Transmission Dynamics of Zika Virus in Island Populations: A Modelling Analysis of the 2013–14 French Polynesia Outbreak. PLoS Negl Trop Dis 10(5): e0004726. https://doi.org/10.1371/journal.pntd.0004726 Editor: Christopher M. Barker, University of California, Davis, UNITED STATES Received: February 5, 2016; Accepted: May 2, 2016; Published: May 17, 2016 Copyright: © 2016 Kucharski et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability: All relevant data are within the paper and its Supporting Information files. Funding: AJK and SF are supported by the Medical Research Council (fellowships MR/K021524/1 and MR/K021680/1). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist.

Introduction Originally identified in Africa [1], the first large reported outbreak of Zika virus (ZIKV) disease occurred in Yap, Micronesia during April–July 2007 [2], followed by an outbreak in French Polynesia between October 2013 and April 2014 [3], and cases in other Pacific countries [4, 5]. During 2015, local transmission was also reported in South American countries, including Brazil [6, 7] and Colombia [8]. Transmission of ZIKV is predominantly vector-borne, but can also occur via sexual contact and blood transfusions [9]. The virus is spread by the Aedes genus of mosquito [10], which is also the vector for dengue virus (DENV). ZIKV is therefore likely to be capable of sustained transmission in other tropical areas [11]. As well as causing symptoms such as fever and rash, ZIKV infection has also been linked to increased incidence of neurological sequelae, including Guillain-Barré Syndrome (GBS) [12, 13] and microcephaly in infants born to mothers who were infected with ZIKV during pregnancy [14]. On 1st February 2015, the World Health Organization declared a Public Health Emergency of International Concern in response to the clusters of microcephaly and other neurological disorders reported in Brazil, possibly linked to the recent rise in ZIKV incidence. The same phenomena were observed in French Polynesia, with 42 GBS cases reported during the outbreak [13, 15]. In addition to the GBS cluster, there were 18 fetal or newborn cases with unusual and severe neurological features reported between March 2014 and May 2015 in French Polynesia [16], including 8 cases with microcephaly [17]. Given the potential for ZIKV to spread globally, it is crucial to characterize the transmission dynamics of the infection. This includes estimates of key epidemiological parameters, such as the basic reproduction number, R 0 , defined as the average number of secondary cases generated by a typical infectious individual in a fully susceptible population, and how many individuals (including both symptomatic and asymptomatic) are typically infected during an outbreak. Such estimates could help assist with outbreak planning, assessment of potential countermeasures, and the design of studies to investigate putative associations between ZIKV infection and other conditions. Islands can be useful case studies for outbreak analysis. Small, centralized populations are less likely to sustain endemic transmission than a large, heterogeneous population [18], which means outbreaks are typically self-limiting after introduction from external sources [19]. Further, if individuals are immunologically naive to a particular pathogen, it is not necessary to consider the potential effect of pre-existing immunity on transmission dynamics [20]. Using a mathematical model of vector-borne infection, we examined the transmission dynamics of ZIKV on six archipelagos in French Polynesia during the 2013–14 outbreak. We inferred the basic reproduction number and the overall size of the outbreak, and hence how many individuals would still be susceptible to infection in coming years.

Methods Data We used weekly reported numbers of suspected ZIKV infections from the six main regions of French Polynesia between 11th October 2013 and 28th March 2014 (Table 1), as detailed in the Centre d’hygiène et de salubrité publique situation reports [21, 22]. Confirmed and suspected cases were reported from sentinel surveillance sites across the country; the number of such sentinel sites varied in number from 27–55 during the outbreak (raw data are provided in S1 Dataset). Clinical cases were defined as suspected cases if they presented to health practitioners with rash and/or mild fever and at least two of the following signs: conjunctivitis, arthralgia, or oedema. Suspected cases were defined as a confirmed case if they tested positive by RT-PCR on blood or saliva. In total, 8,744 suspected cases were reported from the sentinel sites. As there were 162 healthcare sites across all six regions, it has been estimated that around 30,000 suspected cases attended health facilities in total [21]. For each region, we calculated the proportion of total sites that acted as sentinels, to allow us to adjust for variation in reporting over time in the analysis. Population size data were taken from the 2012 French Polynesia Census [23]. In our analysis, the first week with at least one reported case was used as the first observation date. PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Table 1. Geographical breakdown of the 2013–14 French Polynesia ZIKV outbreak. https://doi.org/10.1371/journal.pntd.0004726.t001 Mathematical model We used a compartmental mathematical model to simulate vector-borne transmission [24, 25]. Both people and mosquitoes were modelled using a susceptible-exposed-infectious-removed (SEIR) framework. This model incorporated delays as a result of the intrinsic (human) and extrinsic (vector) incubation periods (Fig 1). Since there is evidence that asymptomatic DENV-infected individuals are capable of transmitting DENV to mosquitoes [26], we assumed the same for ZIKV: all people in the model transmitted the same, regardless of whether they displayed symptoms or were reported as cases. PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 1. Human-vector transmission model schematic. SH represents the number of susceptible people, EH the number of people incubating the virus, IH the number of infectious people, RH the number recovered people. Similarly, SV represents the proportion of mosquitoes currently susceptible, EV the proportion in their incubation period, and IV the proportion of mosquitoes infectious. Mosquitoes are assumed to remain infectious for life. β V is the transmission rate from humans to mosquitoes; β H is transmission from mosquitoes to humans; 1/α H and 1/α V are the mean latent periods for humans and mosquitoes respectively; 1/γ is the mean infectious period for humans; 1/δ is the mean lifespan of mosquitoes; and N is the human population size. https://doi.org/10.1371/journal.pntd.0004726.g001 The main vectors for ZIKV in French Polynesia are thought to be Ae. aegypti and Ae. polynesiensis [12]. In the southern islands, the extrinsic incubation period for Ae. polynesiensis is longer during the cooler period from May to September [27], which may act to reduce transmission. Moreover, temperature can also influence mosquito mortality, and hence vector infectious period [28]. However, climate data from French Polynesia [29] indicated that the ZIKV outbreaks on the six archipelagos ended before a decline in mean temperature or rainfall occurred (S1 Fig). Hence it is likely that transmission ceased as a result of depletion of susceptible humans rather than seasonal changes in vector transmission. Therefore we did not include seasonal effects in our analysis. In the model, SH represents the number of susceptible people, EH is the number of people currently in their incubation period, IH is the number of infectious people, RH is the number of people that have recovered, C denotes the cumulative number of people infected (used to fit the model), and N is the human population size. Similarly, SV represents the proportion of mosquitoes currently susceptible, EV the proportion in their incubation period, and IV the proportion of mosquitoes currently infectious. As the mean human lifespan is much longer than the outbreak duration, we omitted human births and deaths. The full model is as follows: (1) (2) (3) (4) (5) (6) (7) (8) Parameter definitions and values are given in Table 2. We used weakly informative prior distributions for the human latent period, 1/α H , infectious period, 1/γ, extrinsic latent period, 1/α v , and mosquito lifespan, 1/μ. For these prior distributions, we made the assumption that human latent period was equivalent to the intrinsic incubation period, i.e. that no transmission typically occurs before symptom onset. A systematic review of the incubation period for ZIKV in humans estimated a mean value of 5.9 days [30]; the infectious period, 1/γ, lasted for 4–7 days in clinical descriptions of 297 PCR-confirmed cases in French Polynesia [22]; the extrinsic latent period has been estimated at 10.5 days [1]; and mosquito lifespan in Tahiti was estimated at 7.8 days [31]. We therefore used these values for the respective means of 1/α H , 1/γ, 1/α v and 1/δ in our prior distributions. These parameters were estimated jointly across all six regions; as mentioned above, we assumed that the parameters remained fixed over time, as temperature and rainfall levels did not change substantially during the outbreak. The rest of the parameters were estimated for each region individually; we assumed uniform prior distributions for these. PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Table 2. Parameters used in the model. Prior distributions are given for all parameters, along with source if the prior incorporates a specific mean value. All rates are given in units of days−1. https://doi.org/10.1371/journal.pntd.0004726.t002 Serological analysis of samples from blood donors between July 2011 and October 2013 suggested that only 0.8% of the population of French Polynesia were seropositive to ZIKV [33]; we therefore assumed that the population was fully susceptible initially. We also assumed that the initial number of latent and infectious people were equal (i.e. ), and the same for mosquitoes ( ). The basic reproduction number was equal to the product of the average number of mosquitoes infected by the typical infectious human, and vice versa [24]: (9) Statistical inference We fitted the model using Markov chain Monte Carlo (MCMC), where incidence in week t, denoted c t , was the difference in the cumulative proportion of cases over the previous week i.e. c t = C(t) − C(t − 1). In the model, the total number of cases included asymptomatic and subclinical cases—which would not be detected at any site—as well as those that displayed symptoms. Hence there were two sources of potential underreporting: as a result of limited sentinel sites; and as a result of cases not seeking treatment. We adjusted for the first source of underreporting by defining κ t as the proportion of total sites that reported as sentinels in week t. We assumed that the population was uniformly distributed across the catchment areas of the healthcare sites. Under this assumption, the proportion of total sites that reported cases as sentinels in a particular week, κ t , was equivalent to the expected fraction of new cases that would be reported in that week if the reporting proportion, r, was equal to 1. The parameter r accounted for the second source of under-reporting, and represented the proportion of cases (both symptomatic and asymptomatic) that did not seek treatment. To calculate the likelihood of observing a particular number of cases in week t, y t , we assumed the number of confirmed and suspected cases in week t followed a negative binomial distribution with mean rκ t c t and dispersion parameter ϕ, to account for potential variability in reporting over time [34]. The dispersion parameter reflected variation in the overall proportion reported, as well as potential variation in size and catchment area of sentinel sites. Hence the log-likelihood for parameter set θ given data was L(θ|Y) = ∑ t log P(y t |c t ). As a sensitivity analysis (see Results), we also extended the model so the likelihood included the probability of observing 314/476 seropositive individuals in Tahiti after the outbreak, given that a proportion Z were infected in the model. Hence for Tahiti, L(θ|Y) = ∑ t log P(y t |c t ) + log P(X = 314), where X ∼ B(n = 476, p = Z). The joint posterior distribution of the parameters was obtained from eight replicates of 25,000 MCMC iterations, each with a burn-in period of 5,000 iterations (S2–S8 Figs). The model was implemented in R version 3.2.3 [35] using the deSolve package [36]. Demographic model We implemented a simple demographic model to examine the replacement of the number of susceptible individuals over time. In 2014, French Polynesia had a birth rate of b = 15.47 births/1,000 population, a death rate of d = 4.93 deaths/1,000 population, and net migration rate of m = −0.87 migrants/1,000 [37]. The number of susceptible individuals in year τ, S(τ), and total population size, N(τ), was therefore expressed as the following discrete process: (10) (11) We set S(2014) as the fraction of the population remaining in the S compartment at the end of the 2013–14 ZIKV outbreak, and propagated the model forward to estimate susceptibility in future years. The effective reproduction number, R eff (τ), in year τ was the product of the estimated basic reproduction number, and the proportion of the population susceptible: R eff (τ) = R 0 S(τ). We sampled 5,000 values from the estimated joint posterior distributions of S(2014) and R 0 to obtain the median and credible intervals.

Discussion Using a mathematical model of ZIKV transmission, we analysed the dynamics of infection during the 2013–14 outbreak in French Polynesia. In particular, we estimated key epidemiological parameters, such as the basic reproduction number, R 0 , and the proportion of infections that were reported. Across the six regions, our median estimates suggest that between 7–17% of infections were reported as suspected cases. This does not necessarily mean that the non-reported cases were asymptomatic; individuals may have had mild symptoms and hence did not enter the healthcare system. For example, although the attack rate for suspected ZIKV disease cases was 2.5% in the 2007 Yap ZIKV outbreak, a household study following the outbreak found that around 19% of individuals who were seropositive to ZIKV had experienced ZIKV disease-like symptoms during the outbreak period [2]. Our median estimates for R 0 ranged from 2.6–4.8 across the six main archipelagos of French Polynesia, and as a result the median estimates of the proportion of the populations that became infected in our model spanned 87–97%. This is more than the 66% (95% CI: 62–70%) of individuals who were found to be seropositive to ZIKV in a post-outbreak study in Tahiti [17]. When we constrained the model to reproduce this level of seroprevalence as well as the observed weekly reports, however, we obtained a much poorer fit to the case time series (S11 Fig). The discrepancy may be the result of population structure, which we did not include within each region; we used a homogeneous mixing model, in which all individuals had equal chance of contact. In reality, there will be spatial heterogeneity in transmission [41], potentially leading to a depletion of the susceptible human pool in some areas but not in others. Additionally, there is evidence that Ae. aegypti biting rate can vary between individual human hosts [42]. Whereas in the model everyone in regions with ZIKV infected mosquitoes had equal probability of infection, in reality there is likely to be individual-level heterogeneity in probability of infection, which could alter the proportion who seroconvert to ZIKV after the outbreak. As we used a deterministic model, differences in the estimate for the reporting dispersion parameter for different regions may to some extent reflect the limitations of the model in capturing observed transmission patterns, as well as true variability in reporting. The ZIKV outbreak in French Polynesia coincided with a significant increase in Guillain-Barré syndrome (GBS) incidence [13]. We found that although there was a raw incidence rate of 480 (95% CI: 346–648) GBS cases per 100,000 suspected ZIKV cases reported, the majority of the population was likely to have been infected during the outbreak, and therefore the rate per infected person was similar to the overall rate per capita. This could have implications for the design of epidemiological studies to examine the association between ZIKV infection and neurological complications in island populations. If infection with ZIKV confers lifelong immunity, we found it would take at least a decade before re-invasion were possible. In the Pacific island countries and territories, replacement of DENV serotypes occurs every 4–5 years [19, 40], and therefore each specific serotype re-emerges in a 12–15 year cycle. The similarity of this timescale to our results suggest that ZIKV may exhibit very similar dynamics to DENV in island populations, causing infrequent, explosive outbreaks with a high proportion of the population becoming infected. In September 2014, Chikungunya virus (CHIKV) caused a large outbreak in French Polynesia [43], and is another example of a self-limiting arbovirus epidemic in island populations [5]. However, it remains unclear whether ZIKV could become established as an endemic disease in larger populations, as DENV and CHIKV have. For immunising infectious diseases, there is typically a ‘critical community size’, below which random effects frequently lead to disease extinction, and endemic transmission cannot be sustained [18, 44]. Analysis of dengue fever outbreaks in Peru from 1994–2006 found that in populations of more than 500,000 people, dengue was reported in at least 70% of weekly records [41]. Large cities could have the potential to sustain other arboviruses too, and understanding which factors—from population to climate—influence whether ZIKV transmission can become endemic will be an important topic for future research. We did not consider seasonal variation in transmission as a result of climate factors in our analysis, because all six outbreaks ended before there was a substantial seasonal change in rainfall or temperature. Such changes could influence the extrinsic incubation period and mortality of mosquitoes, and hence disease transmission. If the outbreaks had ended as a result of seasonality, rather than depletion of susceptibles, it would reduce the estimated proportion of the population infected, and shorten the time interval before ZIKV would be expected to re-emerge. There are some additional limitations to our analysis. As we were only fitting to a single time series for each region, we also assumed prior distributions for the incubation and infectious periods in humans and mosquitoes. Sensitivity analysis on these prior distributions suggested it was not possible to fully identify these parameters from the available data. If seroprevalence data from each region were to become available in the future, it could provide an indication of how many people were infected, which may make it possible to constrain more of the model parameters, and evaluate the role of spatial heterogeneity discussed above. Such studies may require careful interpretation, though, because antibodies may cross-react between different flaviviruses [12]. Our results suggest that ZIKV transmission in island populations may follow similar patterns to DENV, generating large, sporadic outbreaks with a high degree of under-reporting. If a substantial proportion of such populations become infected during an outbreak, it may take several years for the infection to re-emerge in the same location. A high level of infection, combined with rarity of outbreaks, could also make it more challenging to investigate a potential causal link between infection and concurrent neurological complications.

Author Contributions Conceived and designed the experiments: AJK SF RME WJE. Performed the experiments: AJK. Analyzed the data: AJK SF RME HPM WJE EJN. Wrote the paper: AJK SF RME HPM WJE EJN.