INTRODUCTION
Zika virus (ZIKV) has been highly disruptive in Brazil in recent years, as widespread transmission has been reported in almost every state of the federation [Reference Brasil1, Reference Massad2]. The timing and origin of ZIKV introduction in Brazil, however, has been the subject of controversy.
First, strong evidence has been presented suggesting that ZIKV was imported into Brazil from French Polynesia (FP) [Reference Zanluca3]. The geographic history of ZIKV [Reference Enserink4] indicates that ZIKV resurged in Yap Island in 2007, a date that long preceded the 2015 Brazilian outbreak. Next, the virus spread to FP in 2013, which is the subject of this study. In 2014, ZIKV was introduced to New Caledonia, the Cook Islands and the Easter Islands, all sites with populations, numbers of cases and travel volumes to Brazil considered too small to explain ZIKV exportation. In a recent publication [Reference Lednicky5], Lednicky et al. reported the cases of three children infected with ZIKV in Haiti in 2014. The results of the phylogenetic analysis showed that these children were infected with a ZIKV strain that was similar to the strain identified in Brazil. In addition, the cases were identified before the introduction of ZIKV in Brazil, which suggests that Haiti may have served as an intermediary in the importation of ZIKV between FP and Brazil. However, the force of ZIKV infection in Haiti in 2014 was not sufficiently high to export the virus to Brazil, and it is likely that the virus was first introduced to Haiti and later introduced in Brazil, as both from cases were imported from FP without requiring any intermediate step. Although it is not possible to rule out other sources, our model suggests that FP was the most likely source of ZIKV importation into Brazil.
However, two possible time periods continue to be debated: (i) the 2014 World Cup soccer tournament (12 Jun–13 Jul) [Reference Musso, Nilles and Cao-Lormeau6] and (ii) the Va'a canoe race held in Rio de Janeiro (12–17 Aug 2014) [Reference Musso, Nilles and Cao-Lormeau6], time periods that were further supported by the phylogenetic study conducted by Zanluca et al. [Reference Zanluca3]. However, by August 2014, the epidemic in FP (population 287 679) was over, and although the phylogenetic data in the study conducted by [Reference Musso, Nilles and Cao-Lormeau6] may be correct with respect to origin, the proposed timeline is not compatible with the FP outbreak.
Phylogenetic, epidemiological and mobility data were used in another study to quantify ZIKV evolution and explore the introduction of the virus to the Americas [Reference Musso, Nilles and Cao-Lormeau6]. The results of that study suggest that the introduction of ZIKV to the Americas predated events (i) and (ii), and their molecular clock dates suggested that another event may also be possible: the 2013 Confederations Cup soccer tournament (15–30 Jun 2013). Notably, sports competitors from FP attended all three events. However, an argument against the Confederations Cup soccer tournament in June 2013 is that this event predated the ZIKV outbreak in FP. Based on their findings, the authors conclude that if the ZIKV epidemic in Brazil did indeed arise from a single point of introduction, then the virus must have been circulated in the country for at least 12 months prior to the first case being reported in May 2015.
Mass gathering events have been implicated in many outbreaks of infectious diseases, particularly those with a direct mode of transmission. Every year, millions of people travel around the world to attend a wide range of mass gatherings, including sport events and religious festivals. These events increase the risk of infectious diseases outbreaks, which may be facilitated by the large number of individuals staying in confined areas, sometimes for a short period of time [Reference Tabatabaei and Metanat7]. However, density-dependant transmission, characteristic of directly transmitted infections such as influenza or meningitis, is not a necessary condition for vector-borne infections. In fact, the introduction of ZIKV in Brazil may have occurred independent of the sporting events held in the country prior to the reporting of the first case in 2015, which the model combining the force of infection in likely sources and travel data seems to suggest may be the case.
The objective of this paper was to help resolve the controversy regarding the timing and origin of ZIKV introduction in Brazil. We were motivated to do so because Zika was a highly politically sensitive issue surrounding the 2016 Olympic Games. Furthermore, controversies regarding dengue risk during the FIFA World Cup and Olympic Games [Reference Massad8–Reference Massad, Coutinho and Wilder-Smith11] have demonstrated the importance of evaluating the timing and origin of infectious disease importation for diseases with the potential for international spread. To this end, we used previously published mathematical models to estimate the individual risk of ZIKV infection in the FP Islands, which, when combined with the volume of travellers on commercial flights from that region to Brazilian cities, allowed for the estimation of the number of ZIKV infections imported from FP to Brazil during the period between the end of 2013 and beginning of 2014.
METHODS
Estimating the force of ZIKV infection in the French Polynesian Islands
In the period between October 2013 to April 2014, a large-scale ZIKV outbreak occurred in the FP Islands, with 8598 reported cases [Reference Musso, Nilles and Cao-Lormeau6, Reference Faria12]; however, after considering the fact that 80% of ZIKV infections are subclinical or under-reported [13], it is likely that an estimated 43 000 cases (within a population of 280 000 people) occurred during this time period, representing an infection attack rate of 15%.
Figure 1 shows [Reference Massad8] the weekly number reported ZIKV cases in FP in 2013–2014 [Reference Musso, Cao-Lormeau and Gubler14].
Based on these data, we estimated the incidence of ZIKV infection by assuming an asymptomatic-to-symptomatic ratio of 5 : 1 [Reference Zanluca3], as described below (see [Reference Musso, Nilles and Cao-Lormeau6] and [Reference Amaku15–Reference Burattini17] for the mathematical details).
First, we fitted a continuous function to the number of actually reported infections, representing the time-dependant ZIKV incidence, as follows:
Equation (1) may be used to determine the number of new infections per time unit. In terms of the classical notation of vector-borne infections [ref], the time-dependant incidence of ZIKV is equal to the product of the force of infection, λ(t) multiplied by the number of susceptible humans, denoted S(t). As mentioned later in this paper, the force of infection in vector-borne infections is the product of the biting rate multiplied by the probability of transmission from infected mosquitoes to human hosts multiplied by the number of infected mosquitoes divided by the total number of humans. In equation (1), c 1 is a scale parameter that determines the maximum incidence; c 2 is the time at which the maximum incidence is reached, and c 3 represents the width of the time-dependent incidence function. Equation (1) is intended to reproduce a ‘Gaussian’ curve; therefore, c 1 is just a scale parameter, but c 2 represents the ‘mean’ (and mode or maximum)time, and c 3 represents the ‘variance’ of the time distribution of the cases. F(t) is an ad hoc function introduced to both improve model fit and to set the initial time of infection, c 5. This function has a ‘logistic’ form with the addition of the rect(t) function, as follows:
where c 4 indicates the rate at which the incidence increases; rect(t) is a rectangular function included to account for the slight increase in the risk at the beginning of the year, and c 6 is a parameter that represents the magnitude of this increase. The parameters c i , i = 4, …, 6, therefore, have no physical/biological meaning. However, by including this ad hoc rectangular function, the fit of the model increases by 30%. Without the rect(t) function, the values predicted using equation (1) decrease faster than the actual values. All parameters c i , i = 4, …, 6 were fitted to model (1) so that the force of infection, when applied to the dynamic model described below, reproduced the observed incidence of ZIKV in FP. This was accomplished by applying the prevalence of ZIKV (1) to simple models:
where S H and I H represent the susceptible and infected humans, respectively, and λ(t) is the force of infection (or incidence density rate), which, as mentioned above, represents the product of the mosquitoes biting rate, a; the probability of transmission from mosquitoes to humans, b; and the number of infected mosquitoes with respect to humans, I M/N H, and is normally denoted as λ = ab(I M/N H). However, we would like to note that neither of these parameters/variables have been estimated in this paper. Note that λ(t)S H is the ZIKV incidence. By numerically adjusting models (1) and (3) to the actual data, we determined the values of parameters c i , i = 4, …, 6 that generated λ(t)S(t), that is, the incidence data (reported cases); in other words, the fitted function IncidenceZIKV(t) (equation (1)) is simultaneously used with the set of equations (3) to obtain the incidence, λ(t)S(t). Note that recovery and mortality rates are not considered in model (3) because we are interested only in the notification of cases, irrespective of the eventual subsequent outcomes. This model served as a crude approximation of the risk of ZIKV infection, but the parameter of interest was the force of infection and its consequence in terms of notification rate. We recognise that recovery from infection would reduce the expected number of infected travellers, and hence, our risk estimate is likely an over-estimation of the true risk. However, in a previous publication, Burattini et al. [Reference Burattini17] proposed a complete model that considers recovery and lethality rates in the determination of risk of importation, but that model is much more complicated, and its application to the current set of data did not result in any significant difference in the calculated risk of importation. Thus, we decided to use the current simple model. The calculated incidence, λ(t)S(t), along with the total incidence (taking into account the 5 : 1 ratio mentioned above), are shown in Figure 2.
The small blip observed in January 2014 was the result of the rect(t) function, which was included to account for the slight increase in the risk observed at the beginning of the year, as the model does not consider the inevitable delays in notification. If these delays were considered in the model, they would cause the force of infection to be shifted in time.
Risk was expressed using the following equation:
(equation A7, see Appendix), which can also be interpreted as the probability that at least one individual randomly selected from a population is infected with ZIKV during time interval t 2 − t 1 in a particular region.
To determine the weekly number of travellers leaving FP airports during the period between December 2013 and February 2014, we used data from the IATA (International Air Transport Association), which includes passenger-level flight itineraries describing the full route of the traveller from point of origin to final destination.
Estimating the R 0 for Zika in Brazilian cities
The ability of ZIKV-infected travellers to catalyse autochthonous transmission in Brazilian cities would be dependant upon on the basic reproduction number, R 0 in the cities where the travellers disembarked. We calculated the value of R 0 for dengue between December 2013 and February 2014 for Rio de Janeiro, São Paulo, Fortaleza, Recife, Salvador, Goiania and Vitoria, which together were expected to have received 67% of the expected passengers estimated to be infected with ZIKV. We used the method described in [Reference Massad18] to calculate R 0 DENV. This method estimates R 0 by linearization of the differential equation for the proportion of infected individuals at the beginning of an epidemic. The equation for R 0, in terms of Λ, where Λ is the exponential growth rate at the beginning of an outbreak, is expressed as follows:
where μ M and μ H are the natural mortality rates of mosquitoes and humans, respectively, and γ H is the humans’ recovery rate from infection.
RESULTS
The weekly risk of ZIKV infection (i.e., the probability of identifying at least one infected individual in 1 week) in FP is shown in Figure 3.
Due to seasonal variations in mosquito density, the risk varied with time of the year. As Figure 3 shows, the risk was highest during the weeks with the highest number of reported ZIKV cases (again, as a consequence of the higher density of infected mosquitoes) and reached its maximum value (12%) at the end of December 2013. We would like to note that the potential lag between infection and reporting was neglected in these calculations.
After calculating the weekly risk of infection and given the number of travellers departing FP airports with final destinations in Brazil, it was possible to estimate the expected number of ZIKV-infected individuals arriving at Brazilian airports and the probability of ZIKV introduction.
Table 1 shows the number of airline travellers travelling from FP airports to Brazilian cities between October 2013 and March 2014.
* Source: International Air Transport Association (IATA).
Table 2 shows the individual risk of ZIKV infection in FP and the expected number of infections imported into Brazil by travellers arriving from the islands between December 2013 and February 2014. The number of infections was calculated by multiplying the individual risk of being infected (equation (A7) from the Appendix) by the number of air travellers arriving in Brazilian cities from FP.
* Numbers refers to people arriving at the end of the week.
The manner in which Table 2 has been constructed based on the calculated risk estimates and data on the volume of travel from FP to Brazil requires further clarification. Table 1 shows the weekly number of air passengers from October 2013 to March 2014. The FP ZIKV outbreak lasted from December 2013 until February 2014. A total number of 20 and 14 passengers travelled to São Paulo and Rio de Janeiro in the week ending on 2 December 2013, respectfully, for a total of 34 passengers, as shown in the third column of Table 2. By multiplying the 34 passengers by the individual probability of being infected (0·0979584), we obtain the expected number of travellers infected with ZIKV (3) and its corresponding confidence interval (CI) (2–4). During the next week, 11 individuals travelled from FP to Recife and five individuals travelled from FP to Porto Alegre, for a total of 16 travellers; this value was then multiplied by the risk (0·114501), and the expected number of travellers infected with ZIKV (2) and its corresponding CI (1–2) were calculated. All the remaining estimates provided in Table 2 were calculated using the same algorithm.
Based on these calculations, we estimated that between 11 and 21 ZIKV infections were imported to Brazil from FP via airline travellers between December 2013 and February 2014. While the ZIKV outbreak in FP started in November 2013 and lasted until March 2014, the estimated numbers of infected travellers before December 2013 and after February 2014 were zero, and hence, the time interval in Table 2 was restricted to the period between December 2013 and February 2014. This supports the hypothesis proposed by Faria et al. [Reference Faria12], which suggested that ZIKV may have been imported into Brazil from FP in 2013/2014.
As mentioned in the Methods section, whether ZIKV-infected travellers would trigger autochthonous transmission in Brazilian cities will depend on the basic reproduction number, R 0, in the cities where the travellers disembarked. We applied equation (4) and obtained the reproduction ratio R 0 ZIKV to R 0 DENV from [Reference Funk19] to estimate the likelihood of ZIKV introduction and establishment in those cities.
The result, assuming a (R 0 ZIKV/R 0 DENV) = 1·12 [13], is shown in Table 3.
DISCUSSION
The purpose of modelling the time and origin of ZIKV importation into Brazil is not to propose any causal model in the classical sense of the term but rather to provide more data-driven evidence [Reference Faria12] that indicates FP to be the most likely origin of ZIKV in Brazil. In addition, the results derived based on our model suggest that the exportation of ZIKV from FP to Brazil could have occurred between December 2013 and February 2014. In this sense, our mathematical model, together with published phylogenetic studies, suggests that the recent ZIKV outbreak in Brazil were most likely imported from FP, which was consistent with data from many other publications [Reference Musso, Nilles and Cao-Lormeau6, Reference Massad8–Reference Faria12]. With regard to the time window of introduction, it is highly unlikely that French Polynesian visitors to the FIFA World Cup in 2014, [Reference Zanluca3] or FP athletes [Reference Musso, Nilles and Cao-Lormeau6] competing in the canoe race in Rio de Janeiro in August 2014 could have introduced ZIKV, as the peak of the FP outbreak was already over at the time of those events. Based on the air passenger volume from FP to Brazil, the force of infection in the originating country (FP), and the higher vectorial capacity of Aedes mosquitoes in the summer months of Brazil, we concluded that the period between December 2013 and February 2014 was the most likely time interval during which ZIKV was introduced and established. Although ZIKV infections were reported in FP after March 2014, the estimated risk of infection was very low and resulted in the calculation of zero expected infected travellers. Hence, we restricted our calculations to the period between December 2013 and February 2013.
We recognise that other models [20–Reference Zinszer23] can be applied to test and eventually support the hypothesis that ZIKV was introduced in Brazil from FP. These models, however, are more complex and demand more information than the model we used in this paper. Another limitation of our model is that we did not consider other possible ZIKV sources. However, as mentioned before, the geographic history of ZIKV [Reference Enserink4] indicates that ZIKV resurged in Yap Island in 2007, a date that long preceded the 2015 Brazilian outbreak. Next, the virus spread to FP in 2013, and this outbreak was the subject of this study. In 2014, ZIKV was introduced to New Caledonia, the Cook Islands and the Easter Islands, all sites with populations, numbers of cases and travel volumes to Brazil considered too small to explain ZIKV exportation. Therefore, although it is not possible to rule out other sources, our model suggests that FP was the most likely source of the ZIKV outbreak to Brazil.
2During this time period, ZIKV could have been imported into any of the Brazilian cities mentioned above. Although the R 0 ZIKV value differed in value from city to city, it was greater than one in the seven cities analysed herein. Nevertheless, a recent study [Reference Zinszer23] suggested the presence of a southward pattern of introduction starting from the north-eastern coast, where the cities of Recife, Fortaleza and Salvador are located, and spreading towards the remainder of the country at an average speed of 42 km/day or 15 367 km/year.
If the value of R 0 ZIKV at the time when the infected travellers arrived was of the same order of magnitude as the estimates described in [Reference Fergusson24] for Brazilian cities (approximately 3·5 during the weeks considered), then the likelihood of ZIKV introduction and establishment in Brazil from FP would be even higher.
Finally, we believe it is important to emphasise that we considered only direct transportation from French Polynesia to Brazil. Hence, we cannot exclude the possibility that other indirect routes existed, and it is even possible that index cases were infected elsewhere or that ZIKV was initially transmitted from French Polynesia to another place and then from that place (and any number of additional intermediary locations) to Brazil. Moreover, our definition of ZIKV importation considered all individuals leaving FP and arriving in Brazil to be the sample; for example, an infected individual that departed from FP and travelled Australia and then to Chile prior finally arriving in Brazil was regarded as identical to an infected individual who departed from FP and flew directly to Brazil. Should this individual serve as the index case for an outbreak in any of these intermediary locations, then indirect importation would be important. The potential role of Haiti as an intermediary location that was commented on in the introduction section was not supported by the low force of ZIKV infection observed in that country in 2014. More detailed travel data should be necessary for further refinement of the model.
Therefore, although it was not possible to rule out other sources of ZIKV importation into Brazil because parallel or concurrent importation from other countries may have been possible, we must take into account the fact that surveillance is often variable and imperfect in many places, and only some data are available for analysis. Therefore, despite these limitations, our model suggests that FP was a source of ZIKV in Brazil, which is a possibility that is consistent with the available data.
APPENDIX
Estimating the probability of ZIKV infection among travellers from FP.
We assumed a closed population of size N and the absence of competitive risks, that is, recovery and mortality rates were neglected.
The probability of ZIKV infection can then be calculated using a two-state model without recovery, that is, an SI type model, in which S individuals are susceptible to ZIKV, and I individuals have been infected in the past and acquired lifelong immunity. The model is a finite continuous-time Markov chain (CTMC)Y = {Y(t):t ⩾ 0}, where Y(t) denotes the number of infectives at time t, and its infinitesimal transition probabilities are specified as follows:
where o(ϕ) is a mathematical notation that describes the limiting behaviour of a function when the argument ϕ demonstrates a tendency towards a particular value or infinity, and o(Δt) → 0 as Δt → 0. This means that a typical infective makes infectious contacts at the points of a Poisson process at a rate of λ(t) during an infectious period.
First, we calculated the probability that y individuals are in state I during the period of time between t and t + Δt, as follows (see reference [Reference Lopez25] for details):
In equation (A2), the first term refers to the probability that there were y individuals at time t in condition I (prevalence of notified cases) and that no susceptible individuals (x) acquired the infection in the period. The second term refers to the probability that there were (y − 1) individuals at time t in condition I and that one susceptible individual (x) acquired infection during this period.
Taking the limit of equation (A2) when Δt → 0, it is possible to obtain the Kolmogorov Forward Equation, as follows:
Note that in equation (A3), as mentioned before, we are assuming a closed population of size N = x + y.
The general equation for the Probability Generation Function (PGF), G(z,t), is expressed as follows:
For the specific model described in equation (A3), PGF was expressed in the following manner:
The intermediate steps between equations (A4) and (A5) are numerous, but the interested reader may obtain these details in [Reference Burattini17].
The average number of y individuals at time t can then be calculated by taking the first partial derivative of the PGF with respect to z at z = 1, as follows:
Hence, the average probability (πd(t)) of at least one ZIKV infection at time t may be estimated using the following equation:
The variance of the probability distribution (σ 2) for the number of infected individuals at time t may be calculated as follows:
which results in
To summarise, the risk of ZIKV acquisition at time t, or Risk(t) [Reference Massad, Rocklov and Wilder-Smith26], may be estimated using the following equation:
where the second term is the 95% CI (assuming a normally distributed error). In addition, we also considered the error propagation derived from the fitting procedures of the force of infection, which may be estimated as follows:
where, again, λ(t) is the force of infection or incidence density rate. Note that the risk expressed in equation (A10) indicates the risk of infection among travellers that remain in the ZIKV endemic area during the period between t 1 and t 2. For locals, t 2−t 1 is the time interval considered in the risk calculation (e.g., the week-by-week risk calculation).
ACKNOWLEDGEMENTS
This work has been funded by LIM01-HCFMUSP; CNPq; ZikaPLAN through the European Union's Horizon 2020 Research and Innovation Programme under grant agreement No. 734584; and Canadian Institutes of Health Research (CIHR), consortium Canada-Latin/America-Caribbean Zika Virus Program (Grant 372512).
DECLARATION OF INTEREST
None.