1. Epidemiology and Global Health
  2. Microbiology and Infectious Disease
Download icon

Epidemiological and ecological determinants of Zika virus transmission in an urban setting

  1. José Lourenço Is a corresponding author
  2. Maricelia Maia de Lima
  3. Nuno Rodrigues Faria
  4. Andrew Walker
  5. Moritz UG Kraemer
  6. Christian Julian Villabona-Arenas
  7. Ben Lambert
  8. Erenilde Marques de Cerqueira
  9. Oliver G Pybus
  10. Luiz CJ Alcantara
  11. Mario Recker
  1. University of Oxford, United Kingdom
  2. FIOCRUZ, Brazil
  3. UMI 233, INSERM U1175 and Institut de Biologie Computationnelle, LIRMM, Université de Montpellier, France
  4. Universidade Estadual de Feira de Santana, Brazil
  5. University of Exeter, United Kingdom
Research Article
Cited
1
Views
1,563
Comments
0
Cite as: eLife 2017;6:e29820 doi: 10.7554/eLife.29820

Abstract

The Zika virus has emerged as a global public health concern. Its rapid geographic expansion is attributed to the success of Aedes mosquito vectors, but local epidemiological drivers are still poorly understood. Feira de Santana played a pivotal role in the Chikungunya epidemic in Brazil and was one of the first urban centres to report Zika infections. Using a climate-driven transmission model and notified Zika case data, we show that a low observation rate and high vectorial capacity translated into a significant attack rate during the 2015 outbreak, with a subsequent decline in 2016 and fade-out in 2017 due to herd-immunity. We find a potential Zika-related, low risk for microcephaly per pregnancy, but with significant public health impact given high attack rates. The balance between the loss of herd-immunity and viral re-importation will dictate future transmission potential of Zika in this urban setting.

https://doi.org/10.7554/eLife.29820.001

eLife digest

Mosquitoes can transmit viruses that cause Zika, dengue and several other tropical diseases that affect humans. Zika virus usually causes mild symptoms, but it is thought that infection during pregnancy can lead to brain abnormalities, including microcephaly, where babies are born with an abnormally small head. Recent studies have shed light on how the Zika virus spread from Africa to reach South America, the Caribbean and North America. However, much less is known about the ecological factors that contribute to the spread of the virus within towns, cities and other local areas.

In 2015, Brazil was struck by an outbreak of the Zika virus that led to an international public health emergency. Lourenço et al. used a mathematical model to explore the local conditions within Feira de Santana (a major urban center in Brazil) that contributed to the outbreak. The model took into account numerous factors, including temperature, humidity, rainfall and the mosquito life-cycle, which made it possible to reconstruct the history of the virus over the past three years and to make projections for the next decades.

It revealed that most of the infections occured during 2015, with approximately 65% of the population infected. The incidences of new infections declined in 2016, as increasing numbers of local people had already been exposed to the virus and became immune. Temperature and humidity appeared to have played a critical role in sustaining the mosquito population carrying the Zika virus.

Further analysis suggests that the risk of Zika virus causing microcephaly is very low – only 0.3–0.5% of the pregnant women in Feira de Santana who were infected with Zika gave birth to a baby with the condition. What therefore makes Zika a public health concern is the combination of a low risk with very high infection rates, which can affect a large number of pregnancies.

This study will help researchers and policy makers to predict how the Zika virus will behave in the coming years. It also highlights the limitations and successes of the current system of surveillance. Moreover, it will help to identify critical time periods in the year when mosquito control strategies should be implemented to limit the spread of this virus. In future, this could help shape new local strategies to control Zika virus, dengue and other diseases carried by mosquitoes.

https://doi.org/10.7554/eLife.29820.002

Introduction

The first cases of Zika virus (ZIKV) in Brazil were concurrently reported in March 2015 in Camaçari city in the state of Bahia (Campos et al., 2015) and in Natal, the state capital city of Rio Grande do Norte (Zanluca et al., 2015). During that year, the epidemic in Camaçari quickly spread to other municipalities of the Bahia state, including the capital city of Salvador, which together accounted for over 90% of all notified Zika cases in Brazil in 2015 (Faria et al., 2016a). During this period, many local Bahia health services were overwhelmed by an ongoing Chikungunya virus (CHIKV, East Central South African genotype) epidemic, that was first introduced in 2014 in the city of Feira de Santana (FSA) (Nunes et al., 2015; Faria et al., 2016b). The role of FSA in the establishment and subsequent spread of CHIKV highlights the importance of its socio-demographic and climatic setting, which may well be representative for the transmission dynamics of arboviral diseases in the context of many other urban centres in Brazil and around the world.

On the 1st February 2015 the first ZIKV cases were reported in FSA, followed by a large epidemic that continued into 2016. The rise in ZIKV incidence in FSA coincided temporally with an increase in cases of Guillain-Barré syndrome (GBS) and microcephaly (Faria et al., 2016a), with an unprecedented total of 21 confirmed cases of microcephaly in FSA between January 2015 and May 2017. There is wide statistical support for a causal link between ZIKV and severe manifestations such as microcephaly (Rubin et al., 2016; de Araújo et al., 2016; Soares de Araújo et al., 2016; Honein et al., 2017; Brasil et al., 2016; de Oliveira et al., 2017), and the proposed link in 2015 led to the declaration of the South American epidemic as an international public health emergency by the World Health Organization (WHO) in 2016; the response to which has been limited to vector control initiatives and advice to delay pregnancy in the affected countries (WHO, 2016b; WHO, 2016a). With few cohort studies published (Honein et al., 2017; Brasil et al., 2016) and the lack of an established experimental model for ZIKV infection (Aman and Kashanchi, 2016; Dowall et al., 2016), modelling efforts have taken a central role for advancing our understanding of the virus’s epidemiology (Chowell et al., 2016; Ferguson et al., 2016; Bogoch et al., 2016; Nishiura et al., 2016; Zhang et al., 2016; Perkins et al., 2016; Messina et al., 2016). In particular, our knowledege on parameters of public health importance, such as the basic reproduction number (R0), the duration of infection (Ferguson et al., 2016), attack and reporting rates (Kucharski et al., 2016), the risk of sexual transmission (Maxian et al., 2017; Gao et al., 2016; Moghadas et al., 2017) and birth-associated microcephaly (Bewick et al., 2016; Perkins et al., 2016) has advanced significantly from studies using transmission models. Climate variables are critical for the epidemiological dynamics of Zika and other arboviral diseases, such as dengue (Lourenço and Recker, 2014; Feldstein et al., 2015; Kraemer et al., 2015; van Panhuis et al., 2015) and chikungunya (Poletti et al., 2011; Mourya et al., 2004; Salje et al., 2016). Although these have also been previously addressed in mapping and/or modelling studies (e.g. (Bogoch et al., 2016; Zhang et al., 2016; Perkins et al., 2016; Messina et al., 2016)), their effects as ecological drivers for the emergence, transmission and endemic potential of the Zika virus, especially in the context of a well described outbreak, have not yet been addressed in detail.

In this study, focusing on an urban centre of Brazil (Feira de Santana), we explicitly model the mosquito-vector lifecycle under seasonal, weather-driven variations. Using notified case data of both the number of suspected Zika infections and confirmed microcephaly cases, we demonstrate how the combination of high suitability for viral transmission and a low detection rate resulted in an extremely high attack rate during the first epidemic wave in 2015. The rapid accumulation of herd-immunity significantly reduced the number of cases during the following year, when total ZIKV-associated disease was peaking at the level of the country. Projecting forward we find that the demographic loss of herd-immunity together with the frequency of reintroduction will dictate the risk of reemergence and endemic establishment of Zika in Feira de Santana. The conclusions of this study should be transferable to major urban centres of Brazil and elsewhere with similar climatic and demographic settings.

Methods summary

To model the transmission dynamics of ZIKV infections and estimate relevant epidemiological parameters, we fitted an ento-epidemiological, climate-driven transmission model to ZIKV incidence and climate data of FSA between 2015 and 2017 within a Bayesian framework, similar to our previous work on a dengue outbreak in the Island of Madeira (Lourenço and Recker, 2014).

The model is based on ordinary differential equations (ODE) describing the dynamics of viral infections within the human and mosquito populations (Equations 1-5 and 6-10, respectively). The human population is assumed to be fully susceptible before the introduction of ZIKV and is kept constant in size throughout the period of observation. After an infectious mosquito bite, individuals first enter an incubation phase, after which they become infectious to a mosquito for a limited period of time. Fully recovered individuals are assumed to retain life-long immunity. We assumed that sexual transmission did not significantly contribute to transmission dynamics and therefore ignored its effects (Yakob et al., 2016; Moghadas et al., 2017; Maxian et al., 2017).

For the dynamics of the vector populations we divided mosquitoes into two life-stages: aquatic and adult females. Adult mosquitoes were further divided into the epidemiologically relevant stages for arboviral transmission: susceptible, incubating and infectious. In contrast to human hosts, mosquitoes remain infectious for life. The ODE model comprised 8 climate-dependent entomological parameters (aquatic to adult transition rate, aquatic mortality rate, adult mortality rate, oviposition rate, incubation period, transmission probability to human, hatching success rate and biting rate), whose dependencies on temperature, rainfall and humidity were derived from other studies (see Table 1).

Table 1
Model climate-dependent parameters.
https://doi.org/10.7554/eLife.29820.003
NotationDescription
ϵAv(t)transition rate from aquatic to adult mosquito life-stages
μAv(t)mortality rate of aquatic mosquito life-stage
μVv(t)mortality rate of adult mosquito life-stage
θv(t)(human) intrinsic oviposition rate of adult mosquito life-stage
γv(t)(vector) extrinsic incubation period of adult mosquito life-stage
ϕvh(t)vector-to-human probability of transmission per infectious bite
cv(t)egg hatching success
av(t)adult vector biting rate

Four parameters (baseline mosquito biting rate, mosquito sex ratio, probability of transmission from human-to-vector and human lifespan) were fixed to their expected mean values, taken from the literature (see Table 2). To estimate the remaining parameters, alongside parameter distributions regarding the date of first infection, the human infectious and incubating periods, and the observation rate of notified ZIKV cases, we fitted the ODE model to weekly notified cases of ZIKV in FSA using a Bayesian Markov-chain Monte Carlo (MCMC) approach. The results are presented both in terms of mean dynamic behaviour of the ODE under the MCMC solutions and posterior distributions of key epidemiological parameters. A full description of the fitting approach and the estimated parameters can be found in the section Materials and methods.

Table 2
Model constant parameters.
https://doi.org/10.7554/eLife.29820.004
NotationValueDescriptionReferences
av0.25 per daymosquito biting rate[ 76, 88 ]
f0.5proportion of females (sex ratio)[ 52, 59 ]
ϕhv0.5human-to-vector probability of transmission per infectious bite
1/μh75 yearshuman mean lifespan[ 83 ]

Results

On the 1st February 2015 the first Zika virus (ZIKV) case was reported in Feira de Santana (FSA). Weekly cases remained very low for the following two months, adding up to just 10 notified cases by the end of March that year (Figure 1A). A rapid increase in the number of cases was observed in April, coinciding with Micareta, a local carnival-like festival that takes place across the urban centres of Bahia. The epidemic peaked in July 2015, which was followed by a sharp decline in notified cases over the next 1–2 months. This first epidemic wave was followed by a significantly smaller outbreak in 2016, peaking around March, and an even smaller outbreak in 2017 with no discernable epidemic peak.

Zika virus epidemics in Feira de Santana and Brazil (2015–2016).

(A) Comparison of weekly notified Zika cases (full red line) with monthly Microcephaly cases (blue bars) in Feira de Santana (FSA), overimposed with total Zika cases at the level of the country (BR, black dotted line). BR data for weeks 50–52 was missing. Green area highlights the time period for the Micareta festival and the dotted grey line the date of first notification. Incidence series is available as Dataset 3 and Microcephaly series as Dataset 4. (B) Age distribution and incidence rate ratio (IRR) for the 2015 (blue) and 2016 (green) FSA epidemics (data available as Dataset 2). The top panel shows the number of cases per age (full lines) and the proportion of total cases per age class (dashed lines), which peak at the age range 20–50. The bottom panel shows the age-stratified incidence risk ratio (IRR, plus 95% CI ), with the red dotted line indicating IRR=1. (C) Spatial distribution of cumulative notified cases in BR at the end of 2015 (left) and mid 2016 (right). Two largest urban centres in the Bahia state (Salvador, Feira de Santana) and at the country level (São Paulo, Rio de Janeiro) are highlighted.

https://doi.org/10.7554/eLife.29820.005

Confirmed (and monthly aggregated) microcephaly (MC) cases were absent by November 2015, after which a small epidemic was observed with peak counts in January 2016. We found a time lag of 5–6 months (20–24 weeks) between the first reported Zika epidemic wave and the MC peak in case counts. This coincides with previous observations suggesting a link between the development of neurological complications in newborns and ZIKV infection during the second trimester (Faria et al., 2016a). We note that our lag may be offset by around 1–4 weeks, however, since the date of MC cases in our dataset represents the date of diagnostic confirmation, which is usually done postpartum.

Overall, the epidemic behaviour in FSA was in sharp contrast with trends observed in notified cases across Brazil (BR) as a whole, for which the second epidemic in 2016 was approximately 6 times larger than the one in 2015 (Figure 1A), suggesting the Bahia state as a focus point in the emergence and initial spread of ZIKV in Brazil (Faria et al., 2016c; Faria et al., 2016a). Nonetheless, a clear temporal synchronization between country level and FSA case counts could be observed.

The age distribution of notified ZIKV cases in FSA suggested a higher proportion of cases between 20 and 50 years of age, but with no discernible differences between the two epidemics (Figure 1B, top panel). However, when corrected for the expected number of cases assuming an equal risk of infection per age class, we found the number of cases within this age group to be closer to most other groups (incidence rate ratio, IRR, close to 1, Figure 1B, bottom panel). The per capita case counts within the youngest age class (<1 years) appeared higher than expected, with an IRR significantly above 1 and also higher in 2016 (IRR = 4.4, 95% CI [2.8, 7.0]) than in 2015 (IRR = 1.95, 95% CI [1.5, 2.6]), possibly indicating biased reporting and/or health care seeking with increased awareness of the disease. There was also a consistent trend towards reduced IRR in the elderly (>65 years), although with significant uncertainties. Finally, a small increase in IRR could be detected in the 20–34 year olds, which could potentially be a signature of sexual transmission in this age group (Gao et al., 2016; Carlson et al., 2016; Foy et al., 2011; Turmel et al., 2016; Yakob et al., 2016; Maxian et al., 2017; Moghadas et al., 2017). At this stage and without more detailed data it was not possible to ascertain whether these findings indicated age-related risk of disease, age-dependent exposure risk or simply notification biases in particular age groups, however.

The spatial distribution of total notified cases for BR highlighted the expected clustering of ZIKV cases within the Bahia state by the end of 2015 as well as the wider geographical range by July 2016 (Figure 1C). We speculate that the difference in geographical range could explain the higher number of cases observed during the 2016 epidemic at the country level. This, on the other hand, did not explain why the second epidemic in FSA was nearly 7 times smaller than the first and with only sporadic cases in 2017. To answer this question and to obtain robust parameter estimates of ZIKV epidemiological relevance we utilised a dynamic transmission model, which we fitted to notified case data and local climate variables of FSA within a Bayesian framework (see Materials and methods).

Climate-driven vectorial capacity

The reliance on Aedes mosquitoes for transmission implies that the transmission potential of ZIKV is crucially dependent on temporal trends in the local climate. We therefore investigated daily rainfall, humidity and mean temperature data in FSA between 2013 and May 2017 (Figure 2A). The data showed erratic fluctuations in rainfall with sporadic episodes of intense rain but without a clear seasonal trend. Temperature, on the other hand, presented a much clearer seasonal signature with fixed amplitudes between 22 and 27 degree Celsius, peaking between December and May. Humidity showed an intermediate scenario and appeared correlated with periods of intense rainfall but negatively correlated with temperature.

Figure 2 with 3 supplements see all
Eco-epidemiological factors and model fit to notified cases.

(A) Zika case data (black) and daily climatic series for rainfall (gold), humidity (blue) and mean temperature (green) for Feira de Santana (FSA). Climate data available as Dataset 1. (B) Resulting Bayesian MCMC fit to weekly (black line: data, purple line: model fit) and cumulative incidence (black line: data, grey line: model fit). (A,B) The grey areas highlight the period before the Zika outbreak, the white areas highlight the period for which notified case data was available, and the yellow shaded areas highlight the period for which mean climatic data was used (see Materials and Methods). (C) Climatic series as in A and estimated R0 for the period of the outbreak (2015–2017) (R0 absolute values in Figure 2—figure supplement 3). (D) Correlations between the estimated R0 and climatic variables (intercepts: 0.839 for humidity, 0.067 for rainfall and 0.658 for temperature). (E) Correlations between the case counts and climatic variables (intercepts: 0.487 for humidity, 0.024 for rainfall and 0.862 for temperature). (D,E) Points presented are from timepoints (weeks) for which incidence was notified. (A–E) Y-axis normalised between 0 and 1 for visualisation purposes.

https://doi.org/10.7554/eLife.29820.006

By fitting our climate-driven transmission model to the local climate and ZIKV case data (see Material and Methods and Figure 2B) we obtained parameter estimates for the mosquito lifespan as well as the viral extrinsic incubation period (EIP) for the same period. Mosquito lifespan and EIP are main drivers of vectorial capacity and both showed seasonal oscillations with median values of around 9 and 5 days, respectively (Figure 2—figure supplement 3), which are in line with ranges found in the literature (Trpis et al., 1995; Trpis and Hausermann, 1986; Andraud et al., 2012; Hugo et al., 2014 and Table 3). Importantly, there was a strong negative temporal correlation between these two variables, with periods of longer EIP coinciding with shorter lifespans and vice-versa. This negative relationship resulted in large temporal variations in vectorial capacity and thus seasonal oscillations in the daily reproductive numbers, R0, with a median value of 2.7 in the period 2015–2017 (range 1.0–4.3, Figure 2—figure supplement 3), and 2.2 before 2015, peaking in the local summer months between December and April (Figure 2C). Importantly, R0 remained above 1 for the entire period, indicating a high suitability for ZIKV in FSA. It should be noted that R0 in this context is a time-dependent variable, i.e. R0(t), but out of convenience we simply refer to it as R0.

Table 3
Literature-based reports on key ZIKV epidemiological and entomological parameters.
https://doi.org/10.7554/eLife.29820.010
Parameter/FunctionValues and ranges reportedReferences
Intrinsic incubation period6.5, 5.9 days[ 34, 50 ]
Human infectious period4.7, 9.9 days[ 34, 50 ]
Extrinsic incubation period8.2, <10, <7 days[ 34, 51, 84 ]
Attack rates74, 50, 73, 94, 52%[ 17, 26, 47, 28 ]
R03.2, 2.5, 4.8, 2.05, 2.6–4.8, 4.3–5.8, 1.8–2.0[ 17, 37, 38, 47, 64 ]
Observation rate0.024, 0.06, 0.03, 0.11[ 17, 37, 47 ]

We also looked at the relationship between each climatic variable and R0 and case counts (Figure 2D and E, respectively). The transmission potential was strongly and positively correlated with temperature (r2=0.728) and negatively with humidity (r2=0.26). As expected, from the highly random patterns in the climate series, there was no correlation between R0 and rainfall (r2=0.008). In contrast, there was an opposite trend in the relationship between the climatic variables and case counts, with a positive correlation with humidity (r2=0.28) and a negative correlation with temperature (r2=0.23). As with R0 there was only a weak observable trend in the relationship between rainfall and the number of Zika cases. It should be understood that this macroscopic analysis does not take into account the expected temporal lags due to mosquito development, incubation periods etc., so the purpose here was simply to identify a general qualitative relationship between climate, vectorial capacity and disease incidence.

Model fit and parameter estimates

Four parameters of public health importance were estimated by our MCMC framework: the date of introduction, the human infectious period, the human (intrinsic) incubation period, and the case observation rate (Table 4). The posterior for the introduction date showed a strong support for an introduction into FSA in early-mid December 2014 (estimated median: 10th of December), i.e. around 7–8 weeks before the first notified case (Figure 3A). The estimated human infectious period was 6 days (Figure 3C, median = 5.9, 95% CI [5.47–6.14]), which was very similar to the estimated incubation period (Figure 3D, median = 5.8, 95% CI [5.6–6.15]) and in line with previously estimated ranges for ZIKV (Table 3). In this context it is important to note that informative priors had been used for these 2 parameters (Figure 2—figure supplement 2.), and the posterior for the incubation period presented an adjustment of -0.5 days relative to the proposed distribution from the literature.

Table 4
Model estimated parameters.
https://doi.org/10.7554/eLife.29820.011
NotationDescriptionRanges / priors
t0time point of first case (in a human)(, )
Kaquatic carrying capacity(0, )
ηlinear factor for mosquito mortality(0, )
αlinear factor for extrinsic incubation period(0, )
ρnon-linear factor for effects of humidity and rainfall(0, )
σhhuman infectious period(0, 15)
γhhuman (intrinsic) incubation period(0, 15)
ζobservation rate(0, 1)
Figure 3 with 2 supplements see all
Estimated epidemiological and ecological parameters.

MCMC posterior distributions, based on model fitting to notified case data between 2015–2017 and obtained from sampling 1 million MCMC steps after burn-in. (A) Posterior of the introduction date with median 10th December 2014 (95% CI [01–16 Dec]). (B) Posterior of the observation rate with median 0.0039 (95% CI [0.0038–0.0041])). (C) Posterior of the human infectious period with median 5.9 days (95% CI [5.47–6.14]). (D) Posterior of the human (intrinsic) incubation period with median 5.8 days (95% CI [5.6–6.15]). Representative samples of 500 MCMC chain states are available in Supplementary files 16. See Figure 3—figure supplement 1 for sample chain behaviour.

https://doi.org/10.7554/eLife.29820.012

Of particular interest here was the very low observation rate (Figure 3B), with a median of just under 0.004 (median = 0.0039, 95% CI [0.0038–0.0041]), which equates to less than 4 in 1000 infections having been notified during the epidemic in FSA. Although lower than other previously reported estimates, this would explain the relatively long period of low viral circulation before the epidemic took off in April 2015. That is, based on our estimates, there were around 2,700 Zika infections during the first 2 months, of which only 10 were notified. More importantly, when applying this rate to the total number of cases we found that by the end of the first epidemic wave around 65% (95% CI [57.0–72.9]) of the population in FSA had been infected by the virus. This high attack rate is not unusual for Zika, however, and is in general agreement with observations elsewhere (Table 3).

Future transmission potential for Zika virus

As illustrated by the cumulative attack rate in Figure 4A, and similar to estimates from other regions in the world (Table 3), nearly 65% of the population got infected by ZIKV by the end of 2015, which rose to over 75% (95% CI [76.9–84.3]) by the end of 2016. During the first wave most cases occurred off-season, here defined by our estimated daily reproductive number (R0), while the second wave appeared much more synchronized with the period of high transmission potential. Notably, this temporal phenomenon has also been observed for the chikungunya virus (CHKV) when it was first introduced into FSA in 2014 (Faria et al., 2016b).

Projected Zika virus dynamics and transmission potential.

(A) Fitted and projected epidemic attack rate (% population infected, or herd-immunity, green), basic reproduction number (R0, red) and effective reproduction number (Re, blue).(B) Colourmap showing the projected total number of annual cases depending on rate of external introduction of infectious individuals.The black arrow in the color scale marks the total number of real cases necessary for 1 notified case to be reported in FSA. (C) Projected incidence dynamics when considering less than 1 (green), 2 (blue) and 12 (red) external introductions per year. Grey and white shaded areas delineate different years. The Y-axes are normalised to 1 in each subplot for visualisation purposes. In (B, C) results are based on 1000 stochastic simulations with parameters sampled from the posterior distributions (Figure 3). Representative model solutions for incidence, R0 and Re from 500 MCMC chain samples are available in Supplementary files 16 (both deterministic and stochastic).

https://doi.org/10.7554/eLife.29820.015

The amassed accumulation of herd-immunity during the first wave resulted in a marked difference between the estimated basic reproductive number, R0, and the effective reproductive ratio (Re) by the end of 2015 (Figure 4A). This in turn might explain the marked reduction in Zika cases in FSA in 2016, at a time when the virus was infecting large numbers of individuals elsewhere in the country (Figure 1A,C). At the start of 2016, Re was estimated to be more than 3 times smaller than R0, which increased to 5 by the beginning of 2017. Projecting into the future using average climate data for this region showed that the mean effective reproductive number is expected to remain low and close to 1 for the next few years, suggesting a very weak potential for ZIKV endemicity in the near future. In fact, the sporadic nature of Zika cases in 2017 strongly suggest that herd immunity in this region is at a sufficiently high level to prevent sustained transmission. Furthermore, during 2017, Re was on average less than 1 (mean: 0.62, range: 0.25–1.06), and we would therefore argue that the small number of cases (1.4% of 2015–2017) were mostly a result of small transmission chains, either from resonant transmission from the previous year, or from introduction events from nearby locations. Crucially, this would also explain why our ODE model matched both the dynamics and the sizes of the first two epidemic waves in FSA between 2015 and 2016 but failed to capture the small number of cases during 2017 (Figure 2B).

Without external introductions of infectious individuals (human or vector) our results predicted an epidemic fade-out by 2017, in accordance with the lack of notified cases after March 2017 (Figure 4A). We therefore projected ZIKV’s epidemic potential over the next two decades (until 2040) using stochastic simulations (see Material and Methods) while assuming different rates of viral introduction (Figure 4B,C). Our results showed that the potential for ZIKV to cause another outbreak or to establish itself endemically in FSA is strongly dependent on the frequency of re-introductions, whereby higher rates of external introductions might in fact help to sustain high levels of herd immunity, whereas infrequent introductions are more likely to result in notable outbreaks. That is, semi-endemic behaviour was only observed in simulations with low introduction rates (Figure 4B–C), as these scenarios strike a fine balance between a low number of new cases affecting herd-immunity levels and population turnover. In contrast, high introduction rates quickly exhaust the remaining susceptible pool, resulting in very long periods without epidemic behaviours.

Sensitivity to reporting and microcephaly risk

In effect, our estimated observation rate entails the proportion of real infections that would have been notified if symptomatic and correctly diagnosed as Zika. Based on the previously reported Yap Island epidemic of 2007 (Duffy et al., 2009), the percentage of symptomatic infections can be assumed to be close to 18%. Unfortunately, measures of the proportion of individuals seeking medical attention and being correctly diagnosed do not exist for FSA, although it is well known that correct diagnosis for DENV is imperfect in Brazil (Silva et al., 2016). We therefore performed a sensitivity analysis by varying both the proportions of infected symptomatic individuals seeking medical attention and the proportion of those being correctly diagnosed for Zika. Figure 5A shows that if any of these proportions is less than 10%, or both between 15–20%, our observation rate of 3.9 per 1000 infections can easily be explained.

Sensitivity to reporting and microcephaly risk in Feira de Santana (FSA).

(A) The observation rate (OR) can be expressed as the product of the proportion of cases that are symptomatic (0.18 [Duffy et al., 2009]), with the proportion of symptomatic that seek medical attention, and the proportion of symptomatic that upon medical attention get correctly diagnosed with Zika. In the white area the expected number of notified cases is the range obtained from fitting FSA case data (OR = 0.0039, 95% CI [0.0038–0.0041], Figure 3). (B) Expected number of cases of microcephaly (MC) for theoretical ranges of birth rate (per 1000 females) and risk of MC assuming 65% exposure of all pregnancies as estimated by our model for 2015 in FSA. (C) Expected number of MC per 100,000 individuals under the same conditions as in B. The symbols in B and C represent the total confirmed MC cases (21, red diamond), and the 21 MC plus 3 fetal deaths with confirmed Zika infection (24, white circle); the dashed horizontal line marks the number of births for FSA in 2015, and the vertical lines are the estimated risks per pregnancy.

https://doi.org/10.7554/eLife.29820.016

Finally we investigated the sensitivity of our results with regards to the expected number of newborns presenting microcephaly (MC). Following the observation that virtually all reported MC cases were issued before the summer of 2016 and with a lag of 5–6 months (Figure 1A), we assumed that the vast majority of Zika-associated MC cases would have been a consequence of the first epidemic wave in 2015. We used the estimated attack rate of approximately 65% from 2015 (Figure 4A) and varied the local birth rate and the theoretical risk of MC to obtain an expected number of cases. In agreement with other reports (de Araújo et al., 2016; Cauchemez et al., 2016; Jaenisch et al., 2016; Johansson et al., 2016), our model predicted a relatively low risk for MC given ZIKV infection during pregnancy (Figure 5B,C). In particular, using a conservative total of 21 confirmed MC cases in FSA, i.e. rejecting suspected or other complications, we estimate an average risk of approximately 0.35% of pregnancies experiencing ZIKV infection. Including the 3 foetal deaths where ZIKV infections were confirmed during pregnancy, i.e. using a total of 24 cases, only increased the risk to 0.39%. More generally, based on the results from our fitting approach and using the average birth rates of FSA as guideline, we estimate that on average 3–4 MC cases are expected per 100 k individuals at 65% exposure to the virus.

Discussion

Using an ento-epidemiological transmission model, driven by temporal climate data and fitted to notified case data, we analysed the 2015–2017 Zika outbreak in the city of Feira de Santana (FSA), in the Bahia state of Brazil and determined the conditions that led to the rapid spread of the virus as well as its future endemic and epidemic potential in this region. Given FSA’s high suitability for ZIKV mosquito-vectors and its particular geographical setting as a state commerce and transport hub, our results should have major implications for other urban centres in Brazil and elsewhere.

The pattern of reported ZIKV infections in FSA was characterized by a large epidemic in 2015, in clear contrast to total reports at the country-level, peaking during 2016. Most notably for FSA was the epidemic decay in 2016 and fadeout in 2017. In order to resolve whether this was due to a lower transmission potential of ZIKV in 2016/2017 in FSA, we calculated the daily reproductive number (R0) between 2013 and 2017 but found no notable decrease in 2016. Interestingly, the maximum R0 in that period was observed in the season 2015/2016, coinciding with El Niño (Golden Gate Weather Services, 2017) and thus in line with the hypothesis that this phenomenon may temporary boost arboviral potential (Caminade et al., 2017; van Panhuis et al., 2015). By fitting our model to weekly case data we also estimated the observation rate, i.e. the fraction of cases that were notified as Zika out of the estimated total number of infections. It has previously been reported that the vast majority of Zika infections go unnoticed (Table 3), which is in agreement with our estimates of an observation rate below 1%. Based on this, around 65% of the local population were predicted to have been infected by ZIKV during the first wave in 2015, which is in the same range as the reported Zika outbreaks in French Polynesia (66%) (Cauchemez et al., 2016) and Yap Island (73%) (Duffy et al., 2009). The accumulation of herd-immunity caused a substantial drop in the virus’s effective reproductive number (Re) and hence a significantly lower number of cases during the second wave in 2016 and subsequent demise in 2017. In the context of FSA, it is possible that the high similarity of case definition to DENV, the concurrent CHIKV epidemic, and the low awareness of ZIKV at that time could have resulted in a significant number of ZIKV infections being classified as either dengue or chikungunya. Furthermore, based on our analysis, we would argue that the percentage of correctly diagnosed ZIKV infections and infected individuals seeking medical attention must have been exceptionally low (both lower than 20%).

The age structure of notified cases showed a higher than expected incidence risk ratio (IRR) for individuals under the age of 4 years and a lower than expected risk for individuals aged + 50 years. This contrasts the observation during the Zika outbreak on Yap Island in 2007, where all age classes, except the elderly, presented similar attack rates (Duffy et al., 2009). We note here, however, that the Yap Island analysis was based on both a retrospective analysis of historical hospital records and prospective surveillance (serology, surveys). It is therefore possible that the signatures amongst the youngest and oldest individuals in FSA may reflect deficiencies and/or biases in local notified data. Such signatures could emerge by both a rush of parents seeking medical services driven by a hyped media coverage or prioritization of child-care due to the emergence of microcephaly during the Zika epidemic and a small proportion of the elderly seeking or having access to medical attention. In fact, the increased risk in young children in 2016 may have been a result of increased awareness as well as the interventions by the WHO in the second year. We also found a small increase in IRR in the 20–34 years age group, particularly during 2016, which could be indicative of the small contribution of sexual transmission (Moghadas et al., 2017; Maxian et al., 2017). Most of these observations are speculative, however, and more detailed data will be required to fully understand these age-related risk patterns. For instance, initiatives such as the ZiBRA Project (ZIBRA, 2016; Faria et al., 2016c; Faria et al., 2017), which perform mobile and real-time sampling with portable genome sequencing, could prove to be essential for a retrospective and future analysis of the ZIKV epidemic in Brazil, especially in areas where high levels of herd-immunity will prevent large-scale circulation in the coming years (Ferguson et al., 2016).

The implicit consideration of climate variables as drivers of vector biology allowed us to ascertain the relative roles of temperature, humidity and rainfall for Zika’s basic and effective reproductive potentials (R0 and Re, respectively). Similar to other studies in temperate and tropical settings, we found that temperature, with its direct influence on mosquito lifespan, aquatic development and extrinsic incubation period, was the key driver of seasonal oscillations in the transmission potential (Lourenço and Recker, 2014; Mordecai et al., 2016; Mourya et al., 2004; Feldstein et al., 2015). Rainfall, on the other hand, only seemed to play a marginal role and we argue that it may be a relevant player for arboviral transmission mainly in tropical regions subject to intense rain seasons, such as areas in South East Asia (Cuong et al., 2011; Hii et al., 2012; Xuan et al., 2014). We also noted that the correlations between climatic variables and case counts were inverted when addressed against the transmission potential. For instance, while temperature was positively correlated with R0 it was negatively correlated with Zika cases. This implies that the transmission potential is readily responsive to climatic variation but that the Zika epidemics in FSA showed a slight but expected delay in relation to the peak in transmission potential, with case numbers generally increasing after a stable period of maximum R0, followed by epidemic peaks that tended to coincide with declining R0. An interesting observation is that the 2015 epidemic peaked approximately 3 months after the estimated peak in the virus’s transmission potential, whereas there was much higher synchrony during the second wave in 2016. The same behaviour has been described for the CHIKV outbreak in FSA in 2014–2015 and which has been linked to highly discordant spatial distributions between the first two epidemics (Faria et al., 2016b). It is likely that similar spatial effects (Kraemer et al., 2017) were present in FSA’s ZIKV outbreaks. Unfortunately we did not have access to sufficiently detailed spatial data to explore this hypothesis further.

A phylogenetic analysis has proposed that the introduction of ZIKV into Bahia took place between March and September 2014, although without direct evidence for its circulation in FSA at that time (Naccache et al., 2016). Our estimated date of introduction showed support for a date in early-mid December 2014, a few months after the proposed introduction into Bahia and just over 7 weeks before the first case of Zika was notified in FSA. Similar periods between the first notification and estimated introduction often represent the time taken to complete one or more full transmission cycles (human-mosquito-human) before a cluster of cases is generated of sufficient size for detection by passive surveillance systems (Lourenço and Recker, 2014). The case data also shows a 2 months period after the first notification during which weekly case numbers remained extremely low. This long period was unexpected as persistent circulation of ZIKV could hardly be justified by the observed total of only 10 cases. Given our estimated observation rate, however, the number of ZIKV infections during this time could have amounted to over 2700 actual cases. In April, the number of cases increased rapidly, coinciding with the Micareta festival, which we argue may have played a role in igniting the exponential phase of the epidemic by facilitating human-vector mixing as well as a more rapid geographical expansion.

After calibrating our model to the 2015–2017 epidemic, we projected the transmission of ZIKV beyond 2017 using stochastic simulations and average climatic variables for this region. Without the possibility of externally acquired infections, local extinction was very likely by 2017 due to the high levels of herd-immunity. According to our study, Zika’s reproductive potential (Re) reached its lowest point in 2017, and it is expected to remain low for the next couple of years, given the slow replenishment of susceptibles in the population through births. When explicitly modelling the importation of infectious cases our projections for the next two decades corroborated the conclusions of previous modelling studies that suggest a weak endemic potential for ZIKV after the initial exhaustion of the susceptible pool (Ferguson et al., 2016; Kucharski et al., 2016). However, our simulations also showed that the future epidemic behaviour is strongly dependent on the frequency of re-introductions, where sporadic and unpredictable epidemics could still be in the order of hundreds of cases. Furthermore, given our estimated observation rate for the 2015–2017 epidemic, passive surveillance systems are unlikely to fully detect the scale and occurrence of such small epidemics, missing their actual public health impact, and as such efforts should thus be placed to improve ZIKV detection and diagnosis in order to optimize the local reporting rates and potential for control.

Human sexual and vertical transmission of ZIKV is an important public health concern, especially within the context of potential Zika-associated microcephaly (MC) and other neurological complications in pre- and neonatals. With a total of over 10,000 live births in 2015 in FSA, our crude estimate for the risk of Zika-associated MC per pregnancy was below 4 cases per 100,000 individuals in a generalized population under an attack rate of 65%. As discussed elsewhere (Cauchemez et al., 2016), this risk is extremely low when compared to other known viral-associated complications, such as those caused by infections by cytomegalovirus (CMV) and the rubella virus (RV) (De Santis et al., 2006; Naing et al., 2016). It is therefore crucial to reiterate that what makes the ZIKV a public health concern is not necessarily the per pregnancy risk of neurological complications, but rather the combination of low risk with very high attack rates. Other studies have reported that the risk for complications during the 1st trimester of gestation is higher than the one estimated here. For example, in the French Polynesia (FP) outbreak (Cauchemez et al., 2016), the risk associated with ZIKV infection during the 1st trimester was 1%, while the overall, full pregnancy risk was 0.42%, similar to our FSA estimates. For the Yap Island epidemic, no microcephaly cases have been reported. With an estimated 24 births per 1000 females (census 2000 as in (Duffy et al., 2009)) and using an overall risk of approximately 0.4% per pregnancy, only 0–3 cases per 100,000 individuals would have been expected. However, the island’s small population size (7391 individuals (Duffy et al., 2009)) together with a general baseline of 0–2 microcephaly cases per 100,000 in many areas of the world (Johansson et al., 2016; Butler, 2016; EUROCAT, 2003) would explain the absence of reported cases. It is also important to consider that a variety of birth defects have been found to be statistically associated with Zika virus infection during pregnancy, of which MC is one possible outcome. While the risk for birth defects per pregnancy is consistently reported to be high, estimations for the risk of MC vary considerably. For example, recent clinical trials (Honein et al., 2017; Brasil et al., 2016) suggested that the risk of Zika-associated MC could be an order of magnitude higher than the estimate reported in this or other previous studies (Cauchemez et al., 2016; Duffy et al., 2009). At this stage it is not possible to explain these differences, but it is tempting to speculate that other factors must influence either the actual or estimated risk. For example, there could be diagnostic biases or differences between epidemiological and clinical studies. Alternatively, viral or host genetic background, as well as the pre-exposition to other arboviruses may influence the absolute risk experienced by local populations or cohorts.

Official notification of Zika infections in Brazil started on the 1st of January 2016, although cases were reported in many other regions in Brazil during 2015. It is therefore plausible that the observation rate changed upon official guidelines and that the capacity to accurately diagnose and report Zika infections could have been lower in 2015 compared subsequent years. To explore this, we reran our fitting approach allowing for a possible change in the observation rate for 2016 and onwards (Figure 3—figure supplement 2) and found a similar observation rate for 2015 (0.0039 versus 0.0034) as well as a similar attack rate between the two model variants. However, the estimated observation rate for 2016 and beyond was 4 times larger than for 2015, implying a positive change due to changes in the surveillance system. Nevertheless, only about 13–14 out of 1000 Zika cases were reported after the 1st of January 2016. It is hard to discern where the positive changes took place, but we suggest the revised diagnosis guidelines may have increased the proportion correctly diagnosed while the proportion of symptomatic individuals visiting medical facilities did not change. It is also tempting to speculate that the 2015/2016 imbalance in reporting may have been a general phenomenon across Brazil. As described elsewhere, it is thus possible that FSA is a good example of states and urban centres that may have witnessed larger epidemics than reported in 2015 (de Oliveira et al., 2017). This, together with our conclusion that low MC risk with very high attack rates makes ZIKV a public health concern, could explain why most MC reports at the level of the country were in 2015 (de Oliveira et al., 2017), although for many regions the total reported number of ZIKV cases may have been surprisingly small that year.

There are certain limitations to our approach, many of which could be revisited when more detailed data becomes available. For example, we assumed homogeneous mixing between human and mosquito hosts but it is possible that spatio-temporal heterogeneities may have played a role in FSA. Furthermore, we have curated and integrated functional responses of key entomological parameters to temperature, rainfall and humidity variation, which were originally reported for dengue viruses. Our fitting approach is also dependent on notified case data and it is possible that the reported cases are not representative of the initial expansion of the virus, which may have thwarted the obtained posterior of the introduction date. Finally, our future projections for the endemic and epidemic potential of ZIKV are based on average climatic trends of past years and do not capture the occurrence of natural variation between years, in particular for years affected by major Southern American climate events, such as the El Niño (Caminade et al., 2017).

In this study we have addressed the local determinants of ZIKV epidemiology in the context of a major urban centre of Brazil. Our results imply that control and surveillance of ZIKV should be boosted and focused in periods of high temperature and during major social events. These factors could identify windows of opportunity for local interventions to mitigate ZIKV introduction and transmission and should be transferable to other areas for which both temperature data and community event schedules are available. We further confirm that the high transmission potential of ZIKV in urban centres can lead to the exhaustion of the local susceptible pool, which will in turn dictate the long-term epidemic and endemic behaviour of the virus. Depending on the rate of re-introduction, sporadic outbreaks are to be expected, although these will be unlikely to result in a notable increase in the number of microcephaly cases due to their limited sizes and low risk per pregnancy. Nonetheless, these local sporadic occurrences could still have important public health consequences, and we argue that much better diagnostics and reporting rates are required for local authorities to detect and respond to such events in the near future. Our integrated mathematical framework is capable of deriving key insights into the past and future determinants of ZIKV epidemiology and its findings should be applicable to other major urban centres of Brazil and elsewhere.

Materials and methods

Demographic and socio-economic setting

Feira de Santana (FSA) is a major urban centre of Bahia, located within the state’s largest traffic junction, serving as way points to the South, the Southeast and central regions of the country. The city has a population of approximately 620.000 individuals (2015) and serves a greater geographical setting composed of 80 municipalities (municipios) summing up to a population of 2.5 million. Although major improvements in water supply have been accomplished in recent decades, with about 90% of the population having direct access to piped water, supply is unstable and is common practice to resort to household storage. Together with an ideal (tropical) local climate, these are favourable breeding conditions for species of the Aedes genus of mosquitoes, which are the main transmission vectors of ZIKV, CHIKV and the dengue virus (DENV) that are all co-circulating in the region (Kraemer et al., 2015; Carlson et al., 2016). FSA’s population is generally young, with approximately 30% of individuals under the age of 20% and 60% under the age of 34. In the year of 2015, the female:male sex ratio in FSA was 0.53 and the number of registered births was 10352, leading to a birth rate standard measure of 31 new-borns per 1000 females in the population.

Climate data

Local climatic data (rainfall, humidity, temperature) for the period between January 2013 and May 2017 was collected from the Brazilian open repository for education and research (BDMEP, Banco de Dados Meteorológicos para Ensino e Pesquisa) (Brazil BDMEP, 1961). The climate in FSA is defined as semi-arid (warm but dry), with sporadic periods of rain concetrated within the months of April and July. Between 2013 and 2015, mean yearly temperature was 24.6 celsius (range 22.5–26.6), total precipitation was 856 mm (range 571–1141), and mean humidity levels 79.5% (range 70.1–88.9%). Temperature, humidity and precipitation per day is available as Dataset 1.

Zika virus notified case data

ZIKV surveillance in Brazil is conducted through the national notifiable diseases information system (Sistema de Informação de Agravos de Notificação, SINAN), which relies on passive case detection. Suspected cases are notified given the presence of pruritic maculopapular rash (flat, red area on the skin that is covered with small bumps) together with two or more symptoms among: low fever, or polyarthralgia (joint pain), or periarticular edema (joint swelling), or conjunctival hyperemia (eye blood vessel dilation) without secretion and pruritus (itching) (Brazil SINAN, 2016; Brazil, 2016). The main differences to case definition of DENV and CHIKV are the particular type of pruritic maculopapular rash and low fever (as applied during the Yap Island ZIKV epidemic (Duffy et al., 2009)). The data presented in Figure 1 for both Brazil and FSA represents notified suspected cases and is available as Dataset 3 (please refer to the Acknowledgement section for sources). Here, we use the terms epidemic wave and outbreak interchangeably (but see (Perkins et al., 2016)).

Microcephaly and severe neurological complications case data

A total of 53 suspected cases with microcephaly (MC) or other neurological complications were reported in FSA between January 2015 and February 2017. Using guidelines for microcephaly diagnosis provided in March 2016 by the WHO (as in (Faria et al., 2016c)), a total of 21 cases were confirmed after birth and follow-up. A total of 3 fetal deaths were reported for mothers with confirmed ZIKV infection during gestation but for which no microcephaly assessment was available. The first confirmed microcephaly case was reported on the 24th of November 2015 and virtually all subsequent cases were notified before August 2016 (with the exception of 2). The microcephaly case series can be found in Dataset 4.

Ento-epidemiological dynamic model

The ordinary differential equations (ODE) model and the Markov-chain Monte Carlo (MCMC) fitting approach herein used are based on the framework previously proposed to study the introduction of dengue into the Island of Madeira in 2012 (Lourenço and Recker, 2014). We have changed this framework to relax major modelling assumptions on the mosquito sex ratio and success of egg hatching, have included humidity and rainfall as critical climate variables, and have also transformed the original least squares based MCMC into a Bayesian MCMC. The resulting framework is described in the following sections, in which extra figures are added for completeness.

The dynamics of infection within the human population are defined in Equations 1-5. In summary, the human population is assumed to have constant size (N) with mean life-expectancy of μh years, and to be fully susceptible before introduction of the virus. Upon challenge with infectious mosquito bites (λvh), individuals enter the incubation phase (Eh) with mean duration of 1/γh days, later becoming infectious (Ih) for 1/σh days and finally recovering (Rh) with life-long immunity.

(1)dShdt=μhNλvhμhSh(2)dEhdt=λvhγhEhμhEh(3)dIhdt=γhEhσhIhμhIh(4)dRhdt=σhIhμhRh(5)N=Sh+Eh+Ih+Rh

For the dynamics of the mosquito population (Equations 6-10), individuals are divided into two pertinent life-stages: aquatic (eggs, larvae and pupae, A) and adult females (V) as in (Yang et al., 2009). The adults are further divided into the epidemiologically relevant stages for arboviral transmission: susceptible (Sv), incubating (Ev) for 1/γ˙v days and infectious (Iv) for life. The ˙ (dot) notation is here adopted to distinguish climate-dependent entomological factors (further details in the following sections).

(6)dAdt=c˙vfθ˙v(1AK(R+1))V(ϵ˙Av+μ˙Av)A(7)dSvdt=ϵ˙AvAλhvμ˙VvSv(8)dEvdt=λhvγ˙vEvμ˙VvEv(9)dIvdt=γ˙vEvμ˙VvEv(10)V=Sv+Ev+Iv

Here, the coefficient c˙v is the fraction of eggs hatching to larvae and f the resulting female proportion. For simplicity and lack of quantifications for local mosquito populations, it is assumed that the sex ratio remains at 1:1 (i.e. f=0.5). Moreover, ϵ˙Av denotes the rate of transition from aquatic to adult stages, μ˙Av the aquatic mortality, μ˙Vv the adult mortality, and θ˙v is the success rate of oviposition. The logistic term (1-AK(R+1)) can be understood as the ecological capacity to receive aquatic individuals (Tran et al., 2013), scaled by a carrying capacity term K(R+1) in which K determines the maximum capacity and R is the local rainfall contribution (further details on following sections).

From Equations 6-10, the mean number of viable female offspring produced by one female adult during its life-time, i.e. the basic offspring number Q, was derived (Equation 11). Most parameters defining Q are climate-dependent, and for fixed mean values of the climate variables (ex. mean rainfall R¯), expressions were derived for the expected population sizes of each mosquito life-stage modelled (A0,V0) which are used to initialize the vector population (Equations 12-13).

(11)Q=ϵ˙Avϵ˙Av+μ˙Avc˙fθ˙vμ˙Vv(12)A0=K(R¯+1)(11Q)(13)V0=K(R¯+1)(11Q)ϵAv˙μVv˙

Viral transmission

In respect to the infected host-type being considered, the vector-to-human (λvh) and human-to-vector (λhv) incidence rates are assumed to be, respectively, density-dependent and frequency-dependent (Equations 14-15). Here, av˙ is the biting rate and ϕ˙vh and ϕhv are the vector-to-human and human-to-vector transmission probabilities per bite. Conceptually, this implies that (i) an increase in the density of infectious vectors should directly raise the risk of infection to a single human, while (ii) an increase in the frequency of infected humans raises the risk of infection to a mosquito biting at a fixed rate. The basic reproductive number (R0) is defined similarly to previous modelling approaches (Equation 16) (Wearing and Rohani, 2006; Lourenço and Recker, 2013). We further derived an expression for the effective reproductive ratio (Re, Equation 17), taking into account the susceptible proportion of the population in real-time.

(14)λvh=(av˙ϕ˙vhIvSh/N)Iv(15)λhv=(av˙ϕ˙hvIhSv/N)Ih/N(16)R0=(V/N) av˙ av˙ ϕ˙vh ϕhv γ˙v γhμ˙Vv(σh+μh)(γh+μh)(γ˙v+μ˙Vv)(17)Re=(Sh/N)×(Sv/N)×R0/(V/N)

Markov chain monte carlo fitting approach

For the fitting process, the MCMC algorithm by Lourenco et al. is here altered to a Bayesian approach by formalising a likelihood and parameter priors (Lourenço and Recker, 2014). For this, the proposal distributions (q) of each parameter were kept as Gaussian (symmetric), effectively retaining a random walk Metropolis kernel. We define our acceptance probability α of a parameter set Θ, given model ODE output y as:

(18)α=min{1,π(y|Θ)p(Θ)q(Θo|Θ)π(y|Θo)p(Θo)q(Θ|Θo)}

where Θ and Θo are the proposed and current (accepted) parameter sets (respectively); π(y|Θ) and π(y|Θo) are the likelihoods of the ODE output representing the epidemic data given each parameter set; p(Θo) and p(Θ) are the prior-related probabilities given each parameter set. We fit the Zika virus cumulative case counts per week, for which no age-related or geographical data is taken into consideration.

For computational reasons and based on a previous approach (Dorigatti et al., 2013), the likelihoods π were calculated as the product of the conditional Poisson probabilities of each epidemic data (di) and ODE (yi) data point:

(19)π(y|Θ)=i=1N[Pr{yi=di}]

Note, in this case where we have low cases numbers in a large population, the Poisson likelihood represents a reasonable approximation to the Binomial process, which is expected to underlie the observed data.

Fitted parameters

With the MCMC approach described above, all combinations of the open parameters in the ODE system that most likely represent the outbreak are explored (Table 4). In summary, the MCMC estimates the distributions for: (1) the carrying capacity K, an indirect estimate of the number of adult mosquitoes per human; (2) time point of the first case t0, assumed to be in a human; (3) a linear coefficient η that scales the effect of temperature on aquatic and adult mortality rates; (4) a linear coefficient α that scales the effect of temperature on the extrinsic incubation period; (5) a non-linear coefficient ρ that scales the effects of humidity and rainfall on entoi

mological parameters; (6) the human infectious period 1/σh; and (7) the human incubation period 1/γh.

By introducing the linear coefficients η and α, the relative effect of temperature variation on mortality and incubation is not changed per se, but instead the baselines are allowed to be different from the laboratory conditions used by Yang et al. (Yang et al., 2009). For solutions in which η,α1, the laboratory-based relationships are kept. For a discussion on possible biological factors that may justify η and α please refer to the original description of the method in (Lourenço and Recker, 2014) and (Brady et al., 2013). Finally, the introduction of ρ allows the MCMC to vary the strength by which entomological parameters react to deviations from local humidity and rainfall means. In practice, the effect of rainfall and humidity can be switched off when ρ0 and made stronger when ρ+ (details below).

Initial analysis of the MCMC output raised an identifiability issue between the human infectious period (1/σh) and the linear coefficient (η) that scales the effect of temperature on vector mortality (η scales the baseline mortality without changes to the response of mortality to temperature). Hence, changes in both η and 1/σh result in similar scaling effects on the transmission potential R0 (Equation 16) and thus unstable MCMC chains for η and 1/σh, with the resulting posteriors appearing to be bimodal (for which there was no biological support). We addressed this issue by using informative priors for four parameters for which biological support exists in the literature: η, 1/σh, 1/γh, and α. Gaussian priors were used with means and standard deviations taken from the literature (see Figure 2—figure supplement 2).

Constant parameters

The framework described above has only 4 fixed parameters that are neither climate-dependent nor estimated in the MCMC approach (Table 2). Amongst these, ϕhv is the per bite probability of transmission from human-to-mosquito, which we assume to be 0.5 (Lounibos and Escher, 2008; Mohammed and Chadee, 2011); the sex ratio of the adult mosquito population f is assumed to be 1:1 (Lounibos and Escher, 2008; Mohammed and Chadee, 2011); the life-expectancy of the human population is assumed to be an average of 75 years (WHO, 2016c); and the biting rate is taken to be on average 0.25 although with the potential to vary dependent on humidity levels (details below) (Trpis and Hausermann, 1986; Yasuno and Tonn, 1970).

Climate-Dependent parameters

For each of the temperature-dependent entomological parameters, polynomial expressions are found de novo or taken from previous studies fitting laboratory entomological data with temperature (T) values used in Celsius. For rainfall (R) and humidity (U), positive or negative relationships to entomological parameters are introduced using simple expressions, with values used after normalization to [0,1]. We assume that some parameters are affected by a combination of temperature with either rainfal or humidity, but take their effects to be independent. A list of climate-dependent parameters and references is found in Table 1.

Polynomials of 4th degree for the mortality (μAv,μVv) and success ovipositon (θv) rates are taken from the study by Yang and colleagues under temperature-controlled experiments on populations of Aedes aegypti (Equations 19-21) (Yang et al., 2009). For aquatic to adult (ϵAv) rate we use the 7th degree polynomial of the same study (Equation 20). For the relationship between the extrinsic incubation period (1/γv) and temperature we apply the formulation by Focks et al. which assumes that replication is determined by a single rate-controlling enzyme (Focks et al., 1995; Schoolfield et al., 1981; Otero et al., 2006) (Equation 24). The probability of transmission per mosquito bite (ϕvh) is here modelled (Equation 25) as estimated by Lambrechts and colleagues (Lambrechts et al., 2011). Finally, the relationship between temperature and the fraction of eggs that successfully hatch (cv) is estimated de novo (Equation 26) by fitting a 3rd degree polynomial to Aedes aegypti and albopictus empirical data described by Dickerson et al. (see Figure 2—figure supplement 1) (Dickerson, 2007; Mohammed and Chadee, 2011).

(20)ϵAv(T)=0.1310.05723T+0.01164T20.001341T3+0.00008723T40.000003017T5+5.153×108T63.42×1010T7(21)μAv(T)=2.130.3797T+0.02457T20.0006778T3+0.000006794T4(22)μVv(T)=0.86920.1599T+0.01116T20.0003408T3+0.000003809T4(23)θv(T)=5.4+1.8T0.2124T2+0.01015T30.0001515T4(24)γv(T)=0.003359Tk298×exp(15000R(12981Tk))1+exp(6.203×1021R(12.176×10301Tk))(25)ϕvh(T)=0.001044T×(T12.286)×(32.461T)1/2(26)cv(T)=(184.8+27.94T0.9254T2+0.009226T3)/100.0

We normalise the time series of rainfall (R) and humidity (U), further using the mean normalised values (R¯,U¯) as reference for extreme deviations from the expected local tendencies (Bicout and Sabatier, 2004; Tran et al., 2013). Rainfall is assumed to affect positively the fraction of eggs that successfully hatch (cv) (Alto and Juliano, 2001; Rossi et al., 2015; Tran et al., 2013; Madeira et al., 2002). A similar positive relationship is taken for the vector biting rate (av) and humidity levels (Yasuno and Tonn, 1970), in contrast to a negative effect on the adult mosquito mortality rate (μVv) (Alto and Juliano, 2001).

(27)cv(R)=(RR¯)/1+(RR¯)2(28)av(U)=(UU¯)/1+(UU¯)2(29)μVv(U)=U¯(UU¯)/1+(UU¯)2

Below is the complete formulation for each entomological parameter in time (t), depending on the climatic variables for which relationships are assumed to exist, including the MCMC fitted linear (α,η) and non-linear (ρ) factors described above.

(30)ϵAv(t)=ϵAv(T)(31)μAv(t)=ημAv(T)(32)μVv(t)=ημVv(T)[1+μVv(U)]ρ(33)θv(t)=θv(T)(34)γv(t)=αγv(T)(35)ϕvh(t)=ϕvh(T)(36)cv(t)=cv(T)[1+cv(R)]ρ(37)av(t)=av[1+av(U)]ρ

Stochastic formulation of the ento-epidemiological model

A stochastic version of the ento-epidemiological framework was developed by introducing demographic stochasticity in the transitions of the dynamic system. This followed the original strategy described in (Lourenço and Recker, 2014), in which multinomial distributions are used to sample the effective number of individuals transitioning between classes per time step. Multinomial distributions are generalized binomials - Binomial(n,p) - where n equals the number of individuals in each class and p the probability of the transition event (equal to the deterministic transition rate). This approach has also been demonstrated elsewhere (Lampoudi et al., 2009).

Source code

The approach used in this study uses code in C/C++, bash and R scripts and is available at https://github.com/lourencoj/ArboWeD2/tree/ArboWeD2_V1. A copy is archived at https://github.com/elifesciences-publications/ArboWeD2.

References

  1. 1
  2. 2
  3. 3
  4. 4
    Zika Virus: Endemic Versus Epidemic Dynamics and Implications for Disease Spread in the Americas
    1. S Bewick
    2. WF Fagan
    3. JM Calabrese
    4. F Agusto
    (2016)
    bioRxiv 041897.
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
    Using phenomenological models to characterize transmissibility and forecast patterns and final burden of zika epidemics
    1. G Chowell
    2. D Hincapie-Palacio
    3. J Ospina
    4. B Pell
    5. A Tariq
    6. S Dahal
    7. S Moghadas
    8. A Smirnova
    9. L Simonsen
    10. C Viboud
    (2016)
    PLoS Currents, 8, 10.1371/currents.outbreaks.f14b2217c902f453d9320a43a35b9583, 27366586.
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
    Epidemiology of Chikungunya Virus in Bahia, Brazil, 2014-2015
    1. NR Faria
    2. J Lourenço
    3. E Marques de Cerqueira
    4. M Maia de Lima
    5. L Carlos Junior Alcantara
    (2016)
    PLoS Currents, 8, 10.1371/currents.outbreaks.c97507e3e48efb946401755d468c28b2.
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
    Zika virus dynamics: When does sexual transmission matter?
    1. O Maxian
    2. A Neufeld
    3. EJ Talis
    4. LM Childs
    5. JC Blackwood
    (2017)
    Epidemics, 10.1016/j.epidem.2017.06.003, 28688996.
  56. 56
  57. 57
  58. 58
  59. 59
  60. 60
    Effect of temperature stress on immature stages and susceptibility of Aedes aegypti mosquitoes to chikungunya virus
    1. DT Mourya
    2. P Yadav
    3. AC Mishra
    (2004)
    The American journal of tropical medicine and hygiene 70:346–350.
  61. 61
  62. 62
  63. 63
  64. 64
  65. 65
  66. 66
  67. 67
  68. 68
    The dynamics of temperature- and rainfall-dependent dengue transmission in tropical regions
    1. MM Rossi
    2. LF Lopez
    3. E Massad
    (2015)
    Annals of Biometrics & Biostatistics 2:1–6.
  69. 69
  70. 70
  71. 71
  72. 72
  73. 73
  74. 74
    Doctor of Philosophy
    1. CZ Dickerson
    (2007)
    The effects of temperature and humidy on the eggs of Aedes aegypti and Aedes albopictus in Texas, Doctor of Philosophy, College Station, Texas, Texas A and M University.
  75. 75
  76. 76
  77. 77
  78. 78
  79. 79
  80. 80
  81. 81
    Mosquito (vector) control emergency response and preparedness for Zika virus
    1. WHO
    (2016a)
    World Health Organization. Accessed June 2017.
  82. 82
  83. 83
    World Health Organization - Brazil - Statistics
    1. WHO
    (2016c)
    World Health Organization. Accessed June 2017.
  84. 84
  85. 85
  86. 86
  87. 87
  88. 88
    A study of biting habits of Aedes aegypti in Bangkok, Thailand
    1. M Yasuno
    2. RJ Tonn
    (1970)
    Bulletin of the World Health Organization 43:319–325.
  89. 89
  90. 90
  91. 91

Decision letter

  1. Mark Jit
    Reviewing Editor; London School of Hygiene & Tropical Medicine, and Public Health England, United Kingdom

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

[Editors’ note: a previous version of this study was rejected after peer review, but the authors submitted for reconsideration. The first decision letter after peer review is shown below.]

Thank you for submitting your work entitled "Epidemiological and ecological determinants of Zika virus transmission in an urban setting" for consideration by eLife. Your article has been favorably evaluated by a Senior Editor and three reviewers, one of whom, Mark Jit, is a member of our Board of Reviewing Editors.

Our decision has been reached after consultation between the reviewers. Based on these discussions and the individual reviews below, we regret to inform you that your work in its current form will not be considered further for publication in eLife.

The three reviewers agreed that this is an interesting paper that has the potential to set a standard for ecological-epidemiological analysis of Zika outbreaks in these settings. However, we had serious misgivings about the model fitting process, interpretation of input data and lack of detail about some outputs. Partly as a result of this, we were not persuaded by the findings either.

If these methodological issues could be addressed, and the results can be used to substantiate a much more compelling and convincing story about Zika in Brazil then we would be willing to consider a revised manuscript as a new submission, which may be sent to the same reviewers. However, at this stage we cannot guarantee that we will review a revised manuscript.

If the technical issues can be addressed but the main hypotheses cannot be sustained, then we still think the manuscript has merit and would encourage submitting it to a more specialised journal.

The specific technical issues that would need to be addressed are the following:

a) Issues around the model fit

- The bimodal posterior distribution for infectious and incubation period in Figures 3C-D suggest there could be an issue with the model fitting. Two possibilities come to mind: either a lack of parameter identifiability in the model itself, or poor mixing of the MCMC chains. We would suggest testing for each of these. What do the pairwise correlations between the posteriors look like? What is the effective sample size of the MCMC outputs for each parameter? According to Table 4, there are 8 parameters estimated – what do the posteriors look like for the other 4?

- In Figure 4, the paper states that a stochastic model was used, but this wasn't mentioned in the Materials and methods. How was stochasticity incorporated? In addition, what time of year were the infections introduced in Figure 3C, and how was this value chosen? It seems to me timing would have a big effect on the number of cases, depending on whether introduction co-coincided with a high value of Re.

- In Figure 5A, we did not understand why the entire region of 0-8 cases was shaded blue, rather than just a line representing 4/1000 infections (or perhaps a boundary region to represent the posterior distribution of the estimate).

- Figure 5B is hard to interpret. The colour gradient seems to have been selected so that the line appears to go through the central microcephaly data point of 27, which makes it difficult to identify which regions produce high and low case numbers. It is also not clear what Figure 5C adds, other than normalising the results by the population size – in which case, should the numbers not be 6.2 times smaller (as the population is 620,000)?

- Why was a Poisson likelihood used for the observation process (equation 19), rather than, say, a binomial distribution?

- Why was the polynomial simplified to a 3rd degree one in equation 20? What effect did this have on the model in practice? Similarly, why was a 3rd degree polynomial used to fit the data in Figure 2—figure supplement 1? What impact could this assumption have had on model results?

- Should the observation rate be time-independent? One would expect that surveillance (and health care seeking) would improve as awareness about Zika increased.

b) Issues around interpretation of data and results

- We would have liked to see more discussion around estimates for the α and rho parameters, which control the extent to which environmental factors influence entomological dynamics. What contribution did humidity and temperature have? What are the implications for analysis in other settings, e.g. with stronger or weaker seasonal effects?

- In the fifth paragraph of the Discussion, it seems a stretch to suggest that the model estimates could be consistent with an autumn 2014 introduction. The lower 95% credible interval in Figure 3 is given as 2nd Jan 2015. What proportion of the posterior density falls within the range of dates implied by phylogenetic data?

- The context of the results is seriously misleading about the epidemiology of ZIKV in the rest of Brazil, specifically comparing 2015 to 2016. The state of Bahia is unique in that substantial surveillance was done in 2015. Zika did not become nationally reportable until January 2016, when reported case numbers increased substantially. From both microcephaly reports and anecdotal information, however, it is clear that massive outbreaks occurred in other states in 2015 (e.g. Pernambuco). This is not a problem with the model per se, but it is a big limitation to the conclusion that FSA was different from the rest of Brazil, (e.g. Figure 1A; Results, second paragraph and fourth paragraphs; Discussion, second paragraph). These sections should be rewritten to specifically address the uncertainty in national reporting (almost complete), remove national versus FSA comparison of the epidemiology, and highlight how this model and other models could tell us something about what likely happened in other places. We see that as one of the strengths of this approach and it is not explored at all.

- The finding about infants is notable, but also susceptible to a clear bias. It should be clearly noted that there may be an increased probability of reporting infants for whom care is often prioritized, both by families who would seek it and institutions who would provide it. This likely increased even more after recognition of the association with microcephaly. Furthermore, this dataset may offer a unique opportunity to assess changes in reporting as that became clear. This would be true for both infants and women of reproductive age.

- The information on microcephaly and GBS was insufficient. In the Introduction, it is stated that both were coincident with ZIKV incidence, but that seems unlikely given that both tend to lag behind incidence. These curves should by shown and discussed more specifically as this is a key component to understanding generalizability and the reliability of the data being used.

- There should be a bit more context of other work on ZIKV and climate; it's not accurate to say "the effects of local climate variables, such as temperature and rainfall, have not yet been explored in relation to Zika transmission." That's true in some ways, but not that generalizable. A number of studies are already cited, e.g. Bogoch et al., 2016; Zhang et al., 2016; Perkins et al., 2016; Messina et al., 2016. This manuscript should point out what is unique here.

- The observation rate estimate is very low. Lower than both Yap and French Polynesia. Is there other evidence that supports such a low rate? Limited surveillance? Is it alternatively possible that the epidemic was spatial heterogeneity actually resulted in a smaller epidemic that had a higher reporting rate, lower attack rate, yet nonetheless produced herd immunity effects?

- There should be more discussion about the risk of microcephaly. The comparisons to French Polynesia and Yap are great but there is a lot of other work that has been done, especially clinical studies: https://www.ncbi.nlm.nih.gov/pubmed/26943629 and https://www.ncbi.nlm.nih.gov/pubmed/27960197. It is especially important to understand why the estimates in this manuscript may be on the very low end of what is being reported elsewhere in studies specifically aimed at measuring that risk.

- The analysis suggests that most susceptibles become infected and then immune soon after the first wave of the epidemic in 2015. The second wave in 2016 has a much lower attack rate with a higher proportion of infants. However, there is potential for a new outbreak some years in the future (the exact time is difficult to determine because the x-axes in Figure 4 are incorrectly labelled I think). It would be useful to show the age distribution and predicted microcephaly incidence related to the later outbreaks. If these occur mainly in young children born after 2015 then the public health relevance may be minimal. This has wider implications – does it imply that the long-term public health impact of Zika is minimal once the virus has been established as an endemic childhood infection? These are obviously very large claims that are probably unsustainable from the model in its current state, but without further clarity about results they are obvious extrapolations that readers may make.

c) Issues around reproducibility

- As is normally the rule with eLife modelling papers, the model code, input data and results (including MCMC samples from the converged joint posterior distribution) needed to reproduce the figure should be included as supplementary data files. Public data from cited online sources may be moved, edited or removed in future, so it is important to include everything required to reproduce the descriptive and modelling analysis with the paper itself.

Reviewer #1:

This manuscript fits a model with entomological, epidemiological and climactic variables to data on Zika cases during the 2015/16 outbreaks in one city in Brazil. The model suggests that most susceptibles were infected during the 2015 wave which led to lower incidence in the 2016 wave, and would prevent further epidemics till some years in the future.

Some questions:

1) The analysis suggests that most susceptibles become infected and then immune soon after the first wave of the epidemic in 2015. The second wave in 2016 has a much lower attack rate with a higher proportion of infants. However, there is potential for a new outbreak some years in the future (the exact time is difficult to determine because the x-axes in Figure 4 are incorrectly labelled I think). It would be useful to show the age distribution and predicted microcephaly incidence related to the later outbreaks. If these occur mainly in young children born after 2015 then the public health relevance may be minimal. This has wider implications – does it imply that the long-term public health impact of Zika is minimal once the virus has been established as an endemic childhood infection?

2) However, it is not clear to me exactly how the model fit works, e.g. is the age dependent notification data even used or just the aggregated counts? It would be useful to give the actual likelihood function being used as equation (Cuong et al., 2011) in the appendix is too general (e.g. we aren't told exactly what yi or di are).

3) The relationship between transmission and climactic variables is established via a set of mechanistic equations linking variables governing vector life cycle with climate. While this is sophisticated, it would be useful to see a more conventional multi-variable regression approach, just to ensure that some obvious relationship has not been lost in the detail.

4) Should the observation rate be time-independent? One would expect that surveillance (and health care seeking) would improve as awareness about Zika increased.

5) In Figure 4, it is not clear whether the x-axis in panels A and C are in days or years.

Reviewer #2:

The authors present a transmission modelling analysis of Zika in Feira de Santana, Brazil. I think their broad approach is an important one – combining environmental data in a mechanistic model has the potential to reveal some valuable insights into Zika epidemiology. However, I had some concerns about the robustness of the model estimates, and interpretation of the results.

I have the following comments:

- The bimodal posterior distribution for infectious and incubation period in Figures 3C-D suggest there could be an issue with the model fitting. Two possibilities come to mind: either a lack of parameter identifiability in the model itself, or poor mixing of the MCMC chains. I would suggest testing for each of these. What do the pairwise correlations between the posteriors look like? What is the effective sample size of the MCMC outputs for each parameter? According to Table 4, there are 8 parameters estimated – what do the posteriors look like for the other 4?

- In Figure 4, the authors state they use a stochastic model, but this wasn't mentioned in the Materials and methods. How was stochasticity incorporated? In addition, what time of year were the infections introduced in Figure 3C, and how was this value chosen? It seems to me timing would have a big effect on the number of cases, depending on whether introduction co-coincided with a high value of Re.

- In Figure 5A, I did not understand why the entire region of 0-8 cases was shaded blue, rather than just a line representing 4/1000 infections (or perhaps a boundary region to represent the posterior distribution of the estimate).

- I found Figure 5B hard to interpret. It seems the authors have selected a colour gradient so the line appears to go through the central microcephaly data point of 27, which makes it difficult to identify which regions produce high and low case numbers. It is also not clear to me what Figure 5C adds, other than normalising the results by the population size – in which case, should the numbers not be 6.2 times smaller (as the population is 620,000)?

- I would have liked to see more discussed of estimates for the α and rho parameters, which control the extent to which environmental factors influence entomological dynamics. What contribution did humidity and temperature have? What are the implications for analysis in other settings, e.g. with stronger or weaker seasonal effects?

- In the fifth paragraph of the Discussion, it seems a stretch to suggest that the model estimates could be consistent with an autumn 2014 introduction. The lower 95% credible interval in Figure 3 is given as 2nd Jan 2015. What proportion of the posterior density falls within the range of dates implied by phylogenetic data?

- In the fifth paragraph of the Discussion, the authors suggest they do not have access to spatial data, but Figure 1 indicates they do, at least at some level of resolution. Could they clarify why this is not suitable for exploring heterogeneities to support their discussion point?

- In the subsection “Viral Transmission”, what was the motivation for have density and frequency dependent transmission for vector-human and H-V transmission?

- Why was a Poisson likelihood used for the observation process (equation 19), rather than, say, a binomial distribution?

- Why was the polynomial simplified to a 3rd degree one in equation 20? What effect did this have on the model in practice? Similarly, why was a 3rd degree polynomial used to fit the data in Figure 2—figure supplement 1? What impact could this assumption have had on model results?

Reviewer #3:

The manuscript describes a detailed transmission model of the ZIKV epidemic in the city of Feira de Santa, Brazil. The work is well done and generally interesting, but there were several important shortcomings.

First, the context of the results is seriously misleading about the epidemiology of ZIKV in the rest of Brazil, specifically comparing 2015 to 2016. The state of Bahia is unique in that substantial surveillance was done in 2015. Zika did not become nationally reportable until January 2016, when reported case numbers increased substantially. From both microcephaly reports and anecdotal information, however, it is clear that massive outbreaks occurred in other states in 2015 (e.g. Pernambuco). This is not a problem with the model per se, but it is a big limitation to the conclusion that FSA was different from the rest of Brazil, (e.g. Figure 1A; Results, second and fourth paragraphs; Discussion, second paragraph). These sections should be rewritten to specifically address the uncertainty in national reporting (almost complete), remove national versus FSA comparison of the epidemiology, and highlight how this model and other models could tell us something about what likely happened in other places. I see that as one of the strengths of this approach and it is not explored at all.

The finding about infants is notable, but also susceptible to a clear bias. It should be clearly noted that there may be an increased probability of reporting infants for whom care is often prioritized, both by families who would seek it and institutions who would provide it. This likely increased even more after recognition of the association with microcephaly. Furthermore, this dataset may offer a unique opportunity to assess changes in reporting as that became clear. This would be true for both infants and women of reproductive age.

I also felt the information on microcephaly and GBS was insufficient. In the Introduction, it is stated that both were coincident with ZIKV incidence, but that seems unlikely given that both tend to lag behind incidence. These curves should by shown and discussed more specifically as this is a key component to understanding generalizability and the reliability of the data being used.

There should be a bit more context of other work on ZIKV and climate; it's not accurate to say "the effects of local climate variables, such as temperature and rainfall, have not yet been explored in relation to Zika transmission." That's true in some ways, but not that generalizable. A number of studies are already cited, e.g. Bogoch et al., 2016; Zhang et al., 2016; Perkins et al., 2016; Messina et al., 2016. This manuscript should point out what is unique here.

The observation rate estimate is very low. Lower than both Yap and French Polynesia. Is there other evidence that supports such a low rate? Limited surveillance? Is it alternatively possible that the epidemic was spatial heterogeneity actually resulted in a smaller epidemic that had a higher reporting rate, lower attack rate, yet nonetheless produced herd immunity effects?

Lastly, there should be more discussion about the risk of microcephaly. The comparisons to French Polynesia and Yap are great but there is a lot of other work that has been done, especially clinical studies: https://www.ncbi.nlm.nih.gov/pubmed/26943629 and https://www.ncbi.nlm.nih.gov/pubmed/27960197. It is especially important to understand why the estimates in this manuscript may be on the very low end of what is being reported elsewhere in studies specifically aimed at measuring that risk.

[Editors’ note: what now follows is the decision letter after the authors submitted for further consideration.]

Thank you for resubmitting your work entitled "Epidemiological and ecological determinants of Zika virus transmission in an urban setting" for further consideration at eLife. Your revised article has been favorably evaluated by Prabhat Jha (Senior Editor) and three reviewers, one of whom is a member of our Board of Reviewing Editors.

The work and manuscript have greatly improved, and the reviewers are satisfied with the approach of the model. However, there are some remaining issues about the way aspects of the model and the results are described that need to be addressed before acceptance, as outlined below:

1) The relationship between the Zika epidemic in FSA vs. the rest of Brazil is still not very clear, although it is improved from before. For example, in the last paragraph of the Introduction it states that Zika peaked elsewhere in the country in 2016. That may be true for some locations, but data suggest otherwise for many of the most affected states. As the authors now note, Zika case surveillance changed substantially in 2016. Clearly higher case numbers are associated with this, but microcephaly numbers were much higher for many states as a result of the 2015 epidemics, so in fact many of the NE states in particular likely had much bigger epidemics in 2015, and FSA may be a very good representation of what happened in those states, not an anomaly as implied in the manuscript. The authors even found evidence of a 4-fold increase in reported FSA between 2015 and 2016. In my view, this is a very important finding and emphasizes how little we know about 2015, even in area with relatively strong surveillance.

While we agree that it is not very helpful to speculate about exactly when the Zika epidemic peaked in different parts of Brazil, we think the overall message is still misleading – it currently appears to suggest that the epidemic in FSA was earlier than the rest of the country when the quality of surveillance in 2015 is probably not good enough to make such conclusions. A more appropriate message could be that detailed analysis of FSA indicated a large epidemic in 2015, which could possibly have occurred in other places although this was not picked up by the limited surveillance at that time.

2) The addition of informative priors is sensible. It would also be helpful to have the posterior outputs so the main figures and analysis can be reproduced. This will enable groups working in other parts of Brazil (or indeed other countries) to make use of the analyses to compare to their findings. References should be given for the priors shown in Figure 2—figure supplement 2 – at the moment only the human incubation and infectious period are given priors.

It would also be helpful to have a brief descriptive file that states what data are used in which figure (if not the actual plotting code itself) – e.g. it is not clear where the age distribution in Figure 1B came from.

3) With regards to the actual code, we take the comment that this would run to many thousands of lines. We actually have facilities to host this on eLife, but we would be happy as an alternative if this was archived e.g. in a suitable GitHub repository.

4) The causal link between ZIKV infection and at least some congenital outcomes (e.g. microcephaly) is now clear from multiple studies in multiple locations. While much is left to be learned, we recommend being more direct in the second paragraph of the Introduction.

5) The term 'herd-immunity' seems a little imprecise as immunity has unlikely reached a level that is truly protective against invasion. Perhaps 'population-level immunity' would be more appropriate?

https://doi.org/10.7554/eLife.29820.028

Author response

[Editors’ note: the author responses to the first round of peer review follow.]

[…] The specific technical issues that would need to be addressed are the following:

a) Issues around the model fit

- The bimodal posterior distribution for infectious and incubation period in Figures 3C-D suggest there could be an issue with the model fitting. Two possibilities come to mind: either a lack of parameter identifiability in the model itself, or poor mixing of the MCMC chains. We would suggest testing for each of these. What do the pairwise correlations between the posteriors look like? What is the effective sample size of the MCMC outputs for each parameter? According to Table 4, there are 8 parameters estimated – what do the posteriors look like for the other 4?

It is true that our ODE model is very complex, such that the fitting procedure is likely to suffer from a parameter identifiability issue, which may have led to the seemingly bi-modal posterior distribution. It is also true that there is a certain degree of covariance between the parameters, which again can cause non-normal posteriors. In the original manuscript we used uninformative priors for all of our parameters and relied on empirical observations as a validation of our dynamic model and fitting output – e.g. mean rates of latency and recovery, and R0 values close to reported ranges in the literature. In order to address the possible identifiability issue and shape of the posteriors we now used informed (Gaussian) priors for 4 out of 8 estimated parameters. This was based on relevant biological / epidemiological observations, including human infectious and latency periods and the eta and α factors that determine the vector life-span and incubation periods according to climatic variables.

The inclusion of informative priors resulted in better-behaved posteriors and removed the semi-bimodal shapes from the posteriors for the human infectious and incubation periods (Figure 3). It is important to note, however, that the updated outputs resulting from the addition of the informative priors, as well as the addition of a further 9 months of notified data, are consistent with the ones presented in the original manuscript (see other responses below for comparison of new results under the revised changes to the methodology as proposed by the Editors and Reviewers).

We have updated our Materials and methods section accordingly to accommodate for the inclusion of informative priors. Figure 2—figure supplement 2 and Figure 3—figure supplement 1 now illustrate the posteriors distributions and output from our MCMC chains, demonstrating their stability both in terms of sampling and the levels of correlation between each pair of parameters.

- In Figure 4, the paper states that a stochastic model was used, but this wasn't mentioned in the Materials and methods. How was stochasticity incorporated? In addition, what time of year were the infections introduced in Figure 3C, and how was this value chosen? It seems to me timing would have a big effect on the number of cases, depending on whether introduction co-coincided with a high value of Re.

We apologize for not providing any description about the stochastic implementation, which we have now added to our revised manuscript. In summary, we used multinomial distributions to sample the effective number of individuals transitioning between classes per time step, as implemented in our previous work and elsewhere. A detailed description of this approach can now be found in the Materials and methods section.

With regards to the introduction of infected hosts, this is also modelled stochastically with a fixed migration rate (of infected individuals) per year, which is varied between the sensitivity scenarios presented in Figure 4. In effect, different simulations with the same migration rate will experience introduction events at different time points. The mean and variation presented in the stochastic model output (Figure 4) represent the outcome of thousands of simulations and therefore include all possible scenarios in which introductions may have occurred, covering periods of both high and low Re/ R0.

- In Figure 5A, we did not understand why the entire region of 0-8 cases was shaded blue, rather than just a line representing 4/1000 infections (or perhaps a boundary region to represent the posterior distribution of the estimate).

We agree that the blue area in Figure 5A was not the clearest way to present our results. As suggested by the reviewer, we have now revised the figure to illustrate the exact area that represents Feira de Santana in the background of the sensitivity output (using the 95% CI of the estimated observation rate, taken from the posterior).

- Figure 5B is hard to interpret. The colour gradient seems to have been selected so that the line appears to go through the central microcephaly data point of 27, which makes it difficult to identify which regions produce high and low case numbers. It is also not clear what Figure 5C adds, other than normalising the results by the population size – in which case, should the numbers not be 6.2 times smaller (as the population is 620,000)?

The original Figure 5B was intended to illustrate the sensitivity between the total of MC cases reported and the model estimations, in which the colours were chosen to have the mean at the total of 27 cases (orange line) and variation around this value fading away. The plot would show that under a very strict MC case definition (21, red diamond) or including all suspected counts (36, blue square) the risk of MC per pregnancy is similarly low. This implies that the public health impact of Zika is not necessarily due to high risk for complications per se but rather the potentially explosive nature of Zika epidemics than can affect a significant number of pregnancies. The idea behind Figure 5C was just to complement Figure 5B, showing normalised (per 100,000) MC case numbers. In this context, we confirm that the expected number was indeed 6.2 times smaller.

Following the comments from the reviewers, we have revised Figure 5. We now include a uniform colour scale in Figure 5B to help “identify which regions produce high and low case numbers”, and with the updated MC data set, which has been further curated to include solely confirmed and rejected MC reports, we have also updated the estimated risk per pregnancy based on confirmed cases only. With regards to Figure 5C, we would still argue that it complements Figure 5B, but we would be happy to leave it out of the manuscript if the reviewers / Editors feel it is superfluous.

- Why was a Poisson likelihood used for the observation process (equation 19), rather than, say, a binomial distribution?

The decision of using a Poisson likelihood was based on previous experience, computational cost and the fact that we fit a small number of cases per time step in a large population size, in which case the Poisson distribution is a reasonable approximation for the Binomial. Poisson-based likelihoods have also been used successfully in other modelling approaches, and we have now referenced and justified our decision accordingly in the main text.

- Why was the polynomial simplified to a 3rd degree one in equation 20? What effect did this have on the model in practice? Similarly, why was a 3rd degree polynomial used to fit the data in Figure 2—figure supplement 1? What impact could this assumption have had on model results?

The simplifications in equation 20 and 24 had been done purely for computational reasons but we appreciate that this could raise doubts about the sensitivity of our results on these approximations. To avoid such problems we have revised the equations, now using their original formulations, and have updated the Materials and methods section accordingly. We would like to note, however, that this did not have an effect on our results.

- Should the observation rate be time-independent? One would expect that surveillance (and health care seeking) would improve as awareness about Zika increased.

We agree that the observation rate is likely to have changed in time between 2015 and 2017. In fact, it has now become clear that the Brazilian authorities changed their surveillance system on the 1st January 2016. We therefore explored the effect of using a semi time-dependent observation rate: one for 2015 (zeta) and a different one (zeta’) from 1st January 2016 onwards.

We found that the inclusion of a second observation rate did not change the model fit and results for 2015 and only slightly for 2016. In particular, the posterior of the observation rate in 2015 had virtually the same distribution (mean zeta = 0.0039, mean zeta’ = 0.0034). This resulted in very similar attack rates for 2015 between the two model variants. Consequently, our general conclusions for the epidemic behaviour in 2016-2017 were the same, whereby high levels of herd-immunity generated by the epidemic outbreak in 2015 significantly reduced Zika’s effective reproductive number and hence incidence. In contrast, the posterior mean of the second observation rate (zeta’) showed a 4-fold increase, suggesting that the implementation of the new surveillance system helped to detect / correctly diagnose a higher number of Zika cases.

We have decided to keep the simpler model with one observation rate in the main results of the manuscript. We took this decision on the basis that the observation and attack rates for the year 2015 did not change and because these are the main drivers for the projected behaviour of Zika in the near future. We fully recognize the public health importance of increasing the observation rates due to changes in the surveillance system, but we believe that it should be mainly explored in the Discussion section and supported by Figure 3—figure supplement 2. We are naturally open for guidance and comments on this decision by the Editors and reviewers, however.

b) Issues around interpretation of data and results

- We would have liked to see more discussion around estimates for the α and rho parameters, which control the extent to which environmental factors influence entomological dynamics. What contribution did humidity and temperature have? What are the implications for analysis in other settings, e.g. with stronger or weaker seasonal effects?

We agree that these topics are of great interest and currently lacking in the literature. We have significantly changed Figure 2 and now present two subplots on the relationship between the climatic variables and the observed case counts and transmission potential. We have further added a paragraph discussing these findings in our revised manuscript. However, we believe that a more in-depth exploration would require well-curated arboviral data (as the one presented here) from more than one location, potentially both urban and rural, and that this topic goes beyond the scope of our manuscript. (It is indeed something that we are planning to do as a separate research study in the future).

- In the fifth paragraph of the Discussion, it seems a stretch to suggest that the model estimates could be consistent with an autumn 2014 introduction. The lower 95% credible interval in Figure 3 is given as 2nd Jan 2015. What proportion of the posterior density falls within the range of dates implied by phylogenetic data?

We fully agree that the discussion about the introduction of Zika was highly speculative. In the context of the revised manuscript with changes in the updated dataset and addition of informative priors it is clear that an introduction dating back to early autumn 2014 is not supported. We have therefore removed this point from our revised manuscript.

- The context of the results is seriously misleading about the epidemiology of ZIKV in the rest of Brazil, specifically comparing 2015 to 2016. The state of Bahia is unique in that substantial surveillance was done in 2015. Zika did not become nationally reportable until January 2016, when reported case numbers increased substantially. From both microcephaly reports and anecdotal information, however, it is clear that massive outbreaks occurred in other states in 2015 (e.g. Pernambuco). This is not a problem with the model per se, but it is a big limitation to the conclusion that FSA was different from the rest of Brazil, (e.g. Figure 1A; Results, second paragraph and fourth paragraphs; Discussion, second paragraph). These sections should be rewritten to specifically address the uncertainty in national reporting (almost complete), remove national versus FSA comparison of the epidemiology, and highlight how this model and other models could tell us something about what likely happened in other places. We see that as one of the strengths of this approach and it is not explored at all.

We had no intention to mark Feira de Santana (FSA) as a special case of the Brazilian Zika epidemic. The purpose of the graph in Figure 1 was simply to compare the epidemiological timeseries and to show that the epidemic in FSA peaked in 2015, in contrast to Brazil as a whole, where Zika peaked the following year; this does obviously not exclude other places where the epidemic could also have peaked in 2015. Furthermore, we believe that the behaviour described in this study is potentially ubiquitous to many major urban centres in Brazil and elsewhere. Hence, we do not agree that this is in any way misleading. In our revised version we have tried to clarify that this is not a comparison between FSA and other places in Brazil, which we agree would be very interesting to investigate in more detail, but which would require a different approach and more detailed data from different locations across the country.

- The finding about infants is notable, but also susceptible to a clear bias. It should be clearly noted that there may be an increased probability of reporting infants for whom care is often prioritized, both by families who would seek it and institutions who would provide it. This likely increased even more after recognition of the association with microcephaly. Furthermore, this dataset may offer a unique opportunity to assess changes in reporting as that became clear. This would be true for both infants and women of reproductive age.

We fully agree that the apparent increases or decreases in risks in some of the age groups can potentially be due to reporting bias. We tried to extract as much information out of the data as possible but unfortunately it was too limited to go into further details about the apparent risk differences. In our manuscript we explicitly state that “With the lack of more detailed data it was not possible to ascertain whether these findings indicated age-related risks of disease, age-dependent exposure risks or simply notification biases”. Furthermore, and akin to the reviewer’s comment, in the original Discussion we mentioned that “… such signatures could emerge by both a rush of parents seeking medical services driven by a hyped media coverage during the ZIKV epidemic, and a very small proportion of the elderly seeking or having access to medical attention”. The wording of this paragraph has now been revised to reflect better the reviewer’s comment.

- The information on microcephaly and GBS was insufficient. In the Introduction, it is stated that both were coincident with ZIKV incidence, but that seems unlikely given that both tend to lag behind incidence. These curves should by shown and discussed more specifically as this is a key component to understanding generalizability and the reliability of the data being used.

We agree that the text relating to microcephaly (MC) and Zika infections may have been misleading and potentially incomplete, due to not having had access to the time scales in the MC dataset. To clarify, we did not intend to state that MC and Zika infections should coincide in time, and in fact our work from March 2016 was one of the first to report a lag of several months between the two (Faria et al. (2016) Science 352:aaf5036)). We have now obtained MC counts in time up to May 2017 and have added this data to Figure 1A (unfortunately, as before, data on GBS was not available). By comparing the Zika and MC time series for Feira de Santana we can now detect a lag of about 5.5-6 months, which is slightly longer than our previous report. This, however, can be explained by the fact that the dates for MC cases are the actual dates of confirmation, which is usually done postpartum (30-60 days after delivery). Crucially, as shown in the revised Figure 1A, we are now able to validate our original conclusion that most MC cases reported were a consequence of the high attack rate in 2015 and consequently the high number of pregnancies at risk during that year. We have changed the manuscript to include results and discussion on the new time scale of the MC epidemic in Feira de Santana.

- There should be a bit more context of other work on ZIKV and climate; it's not accurate to say "the effects of local climate variables, such as temperature and rainfall, have not yet been explored in relation to Zika transmission." That's true in some ways, but not that generalizable. A number of studies are already cited, e.g. Bogoch et al., 2016; Zhang et al., 2016; Perkins et al., 2016; Messina et al., 2016. This manuscript should point out what is unique here.

It is true that the studies we reference include some weather proxies or climate variables. However, in our manuscript we present an integrated, dynamic framework that explicitly models the life-cycle of the vector under climatic drivers, alongside viral transmission between mosquitoes and humans. Studies Bogoch et al., 2016, Perkins et al., 2016 and Messina et al., 2016 on the other hand, do not explore a dynamic transmission model and use instead climatic variables for geo-mapping and suitability indexation. The study Zhang et al., 2016 does include a dynamic model and climate-dependent parameters, and although the proposed framework offers a great opportunity to study the impact of climate on Zika transmission, this has not been explored. In order to not overstate the focus of our study we have changed the statement to: “Climate variables are critical for the epidemiological dynamics of Zika and other arboviral diseases, such as dengue [Gao et al., 2016; Bewick et al., 2016; Lourenço and Recker, 2014; Feldstein et al., 2015] and chikungunya [Kraemer et al., 2015; van Punhuis et al., 2011; Poletti et al., 2011]. Although these have also been previously addressed in mapping and / or modelling studies (e.g. [Bogoch et al., 2016; Zhang et al., 2016; Perkins et al., 2016; Messina et al., 2016]), their effects as ecological drivers for the emergence, transmission and endemic potential of the Zika virus, especially in the context of a well described outbreak, have not yet been addressed in detail.”

- The observation rate estimate is very low. Lower than both Yap and French Polynesia. Is there other evidence that supports such a low rate? Limited surveillance? Is it alternatively possible that the epidemic was spatial heterogeneity actually resulted in a smaller epidemic that had a higher reporting rate, lower attack rate, yet nonetheless produced herd immunity effects?

At this point we can only speculate but would argue that it is likely a combination of various factors. The first one is the already mentioned similarity in symptoms to other arboviral diseases, which, given the chikungunya outbreak that took place in 2014-2015, one season before the introduction of the Zika virus, could have led to a high number of misdiagnosed cases. Limited surveillance is certainly another big factor, and as shown in this study could have contributed to the large number of missed cases in 2015-2016. With regards to spatial effects, or rather the lack thereof in our ODE model approach, this is an interesting aspect and we agree that this could have led to an overall overestimation of the attack rate and hence underestimation of the observation rate (as is usually the case with ODE models). Unfortunately we do not have access to higher spatially resolved case data to explore this in more detail but have added this as a Discussion point in our revised version. Because of these uncertainties we present sensitivity analyses over the observation rate (Figure 5A) as well as the number of microcephaly (MC) cases and the associated risk. In this context, we would like to note that model validation comes partially from estimated parameters. That is, although we leave certain parameters free in the model, these are obtained by the fitting procedure and crucially match expectations from the literature (e.g. human infectious period, mosquito life-span, etc.). Furthermore, we obtain a risk of MC per pregnancy that matches previous reports. We therefore argue that although it is possible that several factors are missing from our modelling approach, such small validations are reassuring in the sense that they match expectations and dictate that model calibration by the MCMC approach, based on reported data, fits biological and epidemiological expectations.

- There should be more discussion about the risk of microcephaly. The comparisons to French Polynesia and Yap are great but there is a lot of other work that has been done, especially clinical studies: https://www.ncbi.nlm.nih.gov/pubmed/26943629 and https://www.ncbi.nlm.nih.gov/pubmed/27960197. It is especially important to understand why the estimates in this manuscript may be on the very low end of what is being reported elsewhere in studies specifically aimed at measuring that risk.

We would like to note that our submission preceded the publication of these clinical trials, and thank the reviewer for their reference. In our study we estimate the risk of MC per pregnancy and find it to be consistent with observations in Yap Island and the French Polynesia. The clinical studies mentioned above do not focus on MC alone, but report instead a wide variety of birth defects and associated risks. The general impression from clinical studies is that the risk for birth defects upon Zika infection is indeed high. For instance, Honein et al. report that “Birth defects were reported in 9 of 85 (11%; 95% CI, 6%-19%) completed pregnancies with maternal symptoms or exposure”, while Brasil et al. found that “Among 117 live infants born to 116 ZIKV-positive women, 42% were found to have grossly abnormal clinical or brain imaging findings or both”. Notably, the particular risk of MC in both studies is much lower than the risk for general birth defects (4/117 for Brasil et al., and 18/442 for Honein et al.) but is still ~10 times larger than the risk reported in our manuscript and for the French Polynesia.

At this stage it is difficult, if not impossible, to resolve these differences, but it is tempting to speculate that other (known) factors must influence either the actual or the estimated risk. For example, there may be significant differences in the estimation of Zika attack rates in pregnant women versus the general population due to cross-reactivity of serological assays. There may also be diagnostic biases or differences between epidemiological and clinical studies. Viral or host genetic background, as well as the pre-exposition to other arboviruses, may also influence the absolute risk experienced by local populations, although it is hard to see how this could explain a ten-fold difference in risk. In the absence of a definite explanation we had to resolve to simply discuss this issue in the Discussion of our manuscript, but we would be interested to hear from the reviewer(s) if they had another, plausible explanation.

- The analysis suggests that most susceptibles become infected and then immune soon after the first wave of the epidemic in 2015. The second wave in 2016 has a much lower attack rate with a higher proportion of infants. However, there is potential for a new outbreak some years in the future (the exact time is difficult to determine because the x-axes in Figure 4 are incorrectly labelled I think). It would be useful to show the age distribution and predicted microcephaly incidence related to the later outbreaks. If these occur mainly in young children born after 2015 then the public health relevance may be minimal. This has wider implications – does it imply that the long-term public health impact of Zika is minimal once the virus has been established as an endemic childhood infection? These are obviously very large claims that are probably unsustainable from the model in its current state, but without further clarity about results they are obvious extrapolations that readers may make.

Our estimated cumulative attack rate of over 65% by the end of 2016 would suggest that major outbreaks are unlikely in the near future (the x-axes in Figure 4 have been updated). As our updated dataset shows, there have been a number of Zika cases in 2017, possibly through external introductions or through low-level background transmission events. The nature of these mini-outbreaks will have a big influence on herd-immunity levels in the population, which in turn will dictate, to a certain degree, when a new epidemic would become likely (as shown in Figures 4B, C). However, the potential timing of a new outbreak is too speculative and depends on too many factors for any model to make accurate predictions.

If Zika was to become an endemic disease then the average age of infection would mostly be determined by its reproductive number, as with most other endemic diseases. One could even speculate that the number of MC cases would be relatively low, if the average age of infection was significantly below child bearing age, but that would depend on a much better knowledge about host immunity to Zika as well as the risk of MC in pre-exposed mothers (see discussion above). Given our current knowledge we would feel uncomfortable to make a strong statement about the future occurrence of MC. For these reasons we decided not to enter further into this discussion in our manuscript apart from mentioning that what makes the ZIKV a public health concern is not the per pregnancy risk of neurological complications, but rather the combination of low risk with very high attack rates. High numbers of MC cases should therefore mainly be a phenomenon of major epidemics observed in fully susceptible populations but unlikely in the context of endemic circulation.

c) Issues around reproducibility

- As is normally the rule with eLife modelling papers, the model code, input data and results (including MCMC samples from the converged joint posterior distribution) needed to reproduce the figure should be included as supplementary data files. Public data from cited online sources may be moved, edited or removed in future, so it is important to include everything required to reproduce the descriptive and modelling analysis with the paper itself.

We are happy to provide the data and MCMC samples we analysed in this work as supplementary data files (done in the revised version), and we hope that our new Materials and methods section will be sufficient to allow the interested reader to reproduce the figures under the assumption they have the required computing skills. However, we would question the necessity, or even usefulness, to upload our entire computer code plus script files that contain thousands of lines of code and make use of various programming languages, and therefore like to ask the editor for clarification / guidance.

Reviewer #1:

[…] Some questions:

1) The analysis suggests that most susceptibles become infected and then immune soon after the first wave of the epidemic in 2015. The second wave in 2016 has a much lower attack rate with a higher proportion of infants. However, there is potential for a new outbreak some years in the future (the exact time is difficult to determine because the x-axes in Figure 4 are incorrectly labelled I think). It would be useful to show the age distribution and predicted microcephaly incidence related to the later outbreaks. If these occur mainly in young children born after 2015 then the public health relevance may be minimal. This has wider implications – does it imply that the long-term public health impact of Zika is minimal once the virus has been established as an endemic childhood infection?

Please see our response above.

2) However, it is not clear to me exactly how the model fit works, e.g. is the age dependent notification data even used or just the aggregated counts? It would be useful to give the actual likelihood function being used as equation (Cuong et al., 2011) in the appendix is too general (e.g. we aren't told exactly what yi or di are).

We apologize for not being clear on which data was used. The MCMC fits the ODE model output to the aggregated case counts and no age-related or spatial information is used. In fact, age-related information is not available in time and only as aggregated counts over an entire year (Figure 1B). We have now made this clear in the Materials and methods section and explicitly mention that we fit case counts only, with no consideration of age or spatial information. We have further clarified that yi relates to an ODE data point and di to an observed data point.

3) The relationship between transmission and climactic variables is established via a set of mechanistic equations linking variables governing vector life cycle with climate. While this is sophisticated, it would be useful to see a more conventional multi-variable regression approach, just to ensure that some obvious relationship has not been lost in the detail.

The relationship between climate variables and vector life cycle is taken from the published literature and based on experimental data, and we are unsure how a multi-variable regression approach, especially given the strong non-linear and non-monotonic relationships between many of the variables, could add to our methods / robustness of our results.

4) Should the observation rate be time-independent? One would expect that surveillance (and health care seeking) would improve as awareness about Zika increased.

Please see our response above.

5) In Figure 4, it is not clear whether the x-axis in panels A and C are in days or years.

We apologise for the poor labelling, this has been updated in a revised figure.

Reviewer #2:

[…] I have the following comments:

- The bimodal posterior distribution for infectious and incubation period in Figures 3C-D suggest there could be an issue with the model fitting. Two possibilities come to mind: either a lack of parameter identifiability in the model itself, or poor mixing of the MCMC chains. I would suggest testing for each of these. What do the pairwise correlations between the posteriors look like? What is the effective sample size of the MCMC outputs for each parameter? According to Table 4, there are 8 parameters estimated – what do the posteriors look like for the other 4?

Please see our response above.

- In Figure 4, the authors state they use a stochastic model, but this wasn't mentioned in the Materials and methods. How was stochasticity incorporated? In addition, what time of year were the infections introduced in Figure 3C, and how was this value chosen? It seems to me timing would have a big effect on the number of cases, depending on whether introduction co-coincided with a high value of Re.

Please see our response above.

- In Figure 5A, I did not understand why the entire region of 0-8 cases was shaded blue, rather than just a line representing 4/1000 infections (or perhaps a boundary region to represent the posterior distribution of the estimate).

Please see our response above.

- I found Figure 5B hard to interpret. It seems the authors have selected a colour gradient so the line appears to go through the central microcephaly data point of 27, which makes it difficult to identify which regions produce high and low case numbers. It is also not clear to me what Figure 5C adds, other than normalising the results by the population size – in which case, should the numbers not be 6.2 times smaller (as the population is 620,000)?

Please see our response above.

- I would have liked to see more discussed of estimates for the α and rho parameters, which control the extent to which environmental factors influence entomological dynamics. What contribution did humidity and temperature have? What are the implications for analysis in other settings, e.g. with stronger or weaker seasonal effects?

Please see our response above.

- In the fifth paragraph of the Discussion, it seems a stretch to suggest that the model estimates could be consistent with an autumn 2014 introduction. The lower 95% credible interval in Figure 3 is given as 2nd Jan 2015. What proportion of the posterior density falls within the range of dates implied by phylogenetic data?

Please see our response above.

- In the fifth paragraph of the Discussion, the authors suggest they do not have access to spatial data, but Figure 1 indicates they do, at least at some level of resolution. Could they clarify why this is not suitable for exploring heterogeneities to support their discussion point?

The spatial data presented in Figure 1 is at the level of the municipio, equivalent to councils in the United Kingdom. The areas covered by municipios are large and extend far beyond the area of a city such Feira de Santana. Hence, while we agree that the data would be suitable for a macro-spatial analysis, it cannot be used to explore spatial heterogeneities at the level of an urban centre. Furthermore, the geographic data is aggregated by month, whereas in this study we are exploring the ecological and epidemiological determinant of Zika transmission at a much finer temporal scale.

- In the subsection “Viral Transmission”, what was the motivation for have density and frequency dependent transmission for vector-human and H-V transmission?

This is based on the following observations: for a vector-transmitted disease, an increase in the density of infectious vectors directly raises the risk of infection for a single human host. On the other hand, an increase in the frequency of infected humans raises the risk of infection to a single mosquito host, assuming fixed biting rates.

- Why was a Poisson likelihood used for the observation process (equation 19), rather than, say, a binomial distribution?

Please see our response above.

- Why was the polynomial simplified to a 3rd degree one in equation 20? What effect did this have on the model in practice? Similarly, why was a 3rd degree polynomial used to fit the data in Figure 2—figure supplement 1? What impact could this assumption have had on model results?

Please see our response above.

Reviewer #3:

[…] Lastly, there should be more discussion about the risk of microcephaly. The comparisons to French Polynesia and Yap are great but there is a lot of other work that has been done, especially clinical studies: https://www.ncbi.nlm.nih.gov/pubmed/26943629 and https://www.ncbi.nlm.nih.gov/pubmed/27960197. It is especially important to understand why the estimates in this manuscript may be on the very low end of what is being reported elsewhere in studies specifically aimed at measuring that risk.

Please see our response above to the issues raised by this reviewer.

[Editors' note: the author responses to the re-review follow.]

The work and manuscript have greatly improved, and the reviewers are satisfied with the approach of the model. However, there are some remaining issues about the way aspects of the model and the results are described that need to be addressed before acceptance, as outlined below:

1) The relationship between the Zika epidemic in FSA vs. the rest of Brazil is still not very clear, although it is improved from before. For example, in the last paragraph of the Introduction it states that Zika peaked elsewhere in the country in 2016. That may be true for some locations, but data suggest otherwise for many of the most affected states. As the authors now note, Zika case surveillance changed substantially in 2016. Clearly higher case numbers are associated with this, but microcephaly numbers were much higher for many states as a result of the 2015 epidemics, so in fact many of the NE states in particular likely had much bigger epidemics in 2015, and FSA may be a very good representation of what happened in those states, not an anomaly as implied in the manuscript. The authors even found evidence of a 4-fold increase in reported FSA between 2015 and 2016. In my view, this is a very important finding and emphasizes how little we know about 2015, even in area with relatively strong surveillance.

While we agree that it is not very helpful to speculate about exactly when the Zika epidemic peaked in different parts of Brazil, we think the overall message is still misleading – it currently appears to suggest that the epidemic in FSA was earlier than the rest of the country when the quality of surveillance in 2015 is probably not good enough to make such conclusions. A more appropriate message could be that detailed analysis of FSA indicated a large epidemic in 2015, which could possibly have occurred in other places although this was not picked up by the limited surveillance at that time.

It was indeed our intention to speculate that FSA had a particular behaviour of peaking in 2015, when overall case counts in Brazil peaked in 2016, but we did not intend to focus on FSA being an anomaly – we understand, however, how our statements could be misleading. We also agree that relating our results of low observation rates and potential (estimated) changes of the surveillance system between 2015 and 2016 should be better contextualized for observations across the country. We have therefore followed the Editor’s advice and have now changed our statements in the main text, as highlighted below.

Changes in the Introduction section:

“The rapic accumulation of herd-immunity significantly reduced the number of cases during the following year, when total ZIKV-associated disease was peaking at the level of the country.”

Changes in the Discussion section:

The pattern of reported ZIKV infections in FSA was characterized by a large epidemic in 2015, in clear contrast to total reports at the country-lveel, peaking during 2016. Most notably for FSA was the epidemic decay in 2016 and fadeout in 2017.”

Changes in the Discussion section:

“It is also tempting to speculate that the 2015/2016 imbalance in reporting may have been a general phenomenon across Brazil. […] This, together with our conclusion that low MC risk with very high attack rates makes ZIKV a public health concern, could explain why most MC reports at the level of the country were in 2015 [de Oliveira et al., 2017], although for many regions the total reported number of ZIKV cases may have been surprisingly small that year.”

2) The addition of informative priors is sensible. It would also be helpful to have the posterior outputs so the main figures and analysis can be reproduced. This will enable groups working in other parts of Brazil (or indeed other countries) to make use of the analyses to compare to their findings. References should be given for the priors shown in Figure 2—figure supplement 2 – at the moment only the human incubation and infectious period are given priors.

It would also be helpful to have a brief descriptive file that states what data are used in which figure (if not the actual plotting code itself) – e.g. it is not clear where the age distribution in Figure 1B came from.

A) We would like to note that while we completely agree that making the posteriors available is essential, we argue that the millions of data points used in the figures would make the supplementary files unnecessarily large.

As part of the previous revision, we had submitted supplementary tables (csv files) with 500 samples of the posteriors and corresponding model solutions for R0 and Re (both deterministic and stochastic in a total of 4 files). These samples are representative of the posteriors presented in both Figures 3 and 7. We apologize for not making it clear in the main text that this information had been made available, and have now changed the figure legends for Figures 3 and Figure 2—figure supplement 1. to include reference to the supplementary files.

B)References for eta and α priors were indeed missing – we apologize for this. We have now added these references to the legend of Figure 2—figure supplement 2..

C) We had not included the age-related data in the supplementary files submitted in the last revision. We apologize for this. A new supplementary file has now been added (Dataset 2), which has been listed in the Supplementary files section and mentioned in the legend of Figure 1.

3) With regards to the actual code, we take the comment that this would run to many thousands of lines. We actually have facilities to host this on eLife, but we would be happy as an alternative if this was archived e.g. in a suitable GitHub repository.

We are totally committed for source code to be made public for transparency and research purposes (something the main author has already some background, in terms of publishing R-packages). Indeed, since submission of the last revision we have planned and offered a student project to develop an R-package with the goal of simulating the model described in this manuscript. We believe such project is a more useful way of making the model and methods described in this manuscript useful for the community and into the public domain. We would like to propose that a statement be added to the main text, on the grounds that the source code can be made available by the authors upon request.

“The approach used in this study uses code in C/C++, bash and R scripts. The code can be made available by the authors upon request.”

We note that once the R-package is developed and made public, this current manuscript submitted to eLife will serve as the main case study and first demonstration of its potential.

We are hopeful that the Editors and reviewers understand our position and that this does not affect the final decision from eLife to publish our manuscript.

4) The causal link between ZIKV infection and at least some congenital outcomes (e.g. microcephaly) is now clear from multiple studies in multiple locations. While much is left to be learned, we recommend being more direct in the second paragraph of the Introduction.

We have changed our statement to be more direct about the general consensus of a causal link:

“There is wide statistical support for a causal link between ZIKV and severe manifestations such as microcephaly [Rubin, Greene and Baden, 2016; de Araújo et al., 2016a; de Araújo et al., 2016b; Honein et al., 2017; Brasil et al., 2016; de Oliveira et al., 2017], and the proposed link in 2015 led to the declaration of the South American epidemic…”

5) The term 'herd-immunity' seems a little imprecise as immunity has unlikely reached a level that is truly protective against invasion. Perhaps 'population-level immunity' would be more appropriate?

We are unsure if this comment refers to a particular sentence or the general use of the term ‘herd-immunity’. We note that ‘herd-immunity’ is synonymous to 'population-level immunity' (and not a particular level of population-level immunity that offers protection). We have therefore applied no changes in regards to the term herd-immunity, but are happy to address this further if the editors/reviewers feel strongly about it.

https://doi.org/10.7554/eLife.29820.029

Article and author information

Author details

  1. José Lourenço

    Department of Zoology, University of Oxford, Oxford, United Kingdom
    Contribution
    Conceptualization, Data curation, Formal analysis, Methodology
    For correspondence
    jose.lourenco@zoo.ox.ac.uk
    Competing interests
    No competing interests declared
    ORCID icon 0000-0002-9318-2581
  2. Maricelia Maia de Lima

    Laboratory of Haematology, Genetics and Computational Biology, FIOCRUZ, SalvadorBahia, Brazil
    Contribution
    Conceptualization, Data curation, Investigation, Writing—original draft, Writing—review and editing
    Competing interests
    No competing interests declared
  3. Nuno Rodrigues Faria

    Department of Zoology, University of Oxford, Oxford, United Kingdom
    Contribution
    Conceptualization, Resources, Supervision, Investigation, Writing—original draft, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon 0000-0001-8839-2798
  4. Andrew Walker

    Department of Zoology, University of Oxford, Oxford, United Kingdom
    Contribution
    Data curation, Funding acquisition, Project administration
    Competing interests
    No competing interests declared
  5. Moritz UG Kraemer

    Department of Zoology, University of Oxford, Oxford, United Kingdom
    Contribution
    Conceptualization, Data curation
    Competing interests
    No competing interests declared
    ORCID icon 0000-0001-8838-7147
  6. Christian Julian Villabona-Arenas

    Institut de Recherche pour le Développement, UMI 233, INSERM U1175 and Institut de Biologie Computationnelle, LIRMM, Université de Montpellier, Montpellier, France
    Contribution
    Conceptualization, Investigation, Methodology
    Competing interests
    No competing interests declared
    ORCID icon 0000-0001-9928-3968
  7. Ben Lambert

    Department of Zoology, University of Oxford, Oxford, United Kingdom
    Contribution
    Conceptualization, Investigation, Methodology
    Competing interests
    No competing interests declared
  8. Erenilde Marques de Cerqueira

    Centre of PostGraduation in Collective Health, Department of Health, Universidade Estadual de Feira de Santana, Feira de SantanaBahia, Brazil
    Contribution
    Conceptualization, Data curation
    Competing interests
    No competing interests declared
  9. Oliver G Pybus

    Department of Zoology, University of Oxford, Oxford, United Kingdom
    Contribution
    Resources, Data curation, Supervision
    Competing interests
    No competing interests declared
  10. Luiz CJ Alcantara

    Laboratory of Haematology, Genetics and Computational Biology, FIOCRUZ, SalvadorBahia, Brazil
    Contribution
    Conceptualization, Data curation, Funding acquisition, Investigation, Writing—original draft, Project administration, Writing—review and editing
    Competing interests
    No competing interests declared
  11. Mario Recker

    Centre for Mathematics and the Environment, University of Exeter, Penryn, United Kingdom
    Contribution
    Conceptualization, Data curation, Investigation, Writing—original draft, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon 0000-0001-9489-1315

Funding

European Research Council (268904 - DIVERSITY)

  • José Lourenço
  • Andrew Walker

Wellcome Trust and Royal Society (204311/Z/16/Z)

  • Nuno Rodrigues Faria

International Development Emerging Pandemic Threats Program-2 (AID-OAA-A-14-00102)

  • Moritz UG Kraemer

Labex EpiGenMed (ANR-10-LABX-12-01)

  • Christian Julian Villabona-Arenas

Engineering and Physical Sciences Research Council

  • Ben Lambert

European Research Council (614725-PATHPHYLODYN)

  • Oliver G Pybus

Royal Society

  • Mario Recker

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

Acknowledgements

We are most grateful to Wanderson Klebeler de Oliveira, Livia Carla Vinhal, Mariana Pastorello Verotti, Giovanini Evelim Coelho and Claudio Maierovitch Pessanha Henrique from the Brazilian Ministry of Health for providing epidemiological data regarding Zika virus notified cases in Brazil. MML and EMC curated the Zika virus notified cases in Feira de Santana. JL and ASW received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement no. 268904 – DIVERSITY. MR was supported by a Royal Society University Research Fellowship. The European Research Council under the European Union’s Seventh Framework Program (FP7/2007-2013)/ERC grant agreement no. 614725-PATHPHYLODYN funded OP. MUGK’s contribution was made possible by the generous support of the American people through the United States Agency for International Development Emerging Pandemic Threats Program-2 PREDICT-2 (Cooperative Agreement No. AID-OAA-A-14–00102). CJVA was supported by a fellowship from the Labex EpiGenMed, via the National Research Agency, Program for Future Investment and University of Montpellier [ANR-10-LA-12–01]. BL received funding from the Engineering and Physical Sciences Research Council (EPSRC) in the UK. NRF was supported by a Sir Henry Dale Fellowship jointly funded by the Wellcome Trust and the Royal Society (grant number 204311/Z/16/Z). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Reviewing Editor

  1. Mark Jit, Reviewing Editor, London School of Hygiene & Tropical Medicine, and Public Health England, United Kingdom

Publication history

  1. Received: June 21, 2017
  2. Accepted: September 4, 2017
  3. Accepted Manuscript published: September 9, 2017 (version 1)
  4. Version of Record published: October 12, 2017 (version 2)

Copyright

© 2017, Lourenço 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.

Metrics

  • 1,563
    Page views
  • 403
    Downloads
  • 1
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Comments

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)