Estimating risks of importation and local transmission of Zika virus infection

Background. An international spread of Zika virus (ZIKV) infection has attracted global attention. ZIKV is conveyed by a mosquito vector, Aedes species, which also acts as the vector species of dengue and chikungunya viruses. Methods. Arrival time of ZIKV importation (i.e., the time at which the first imported case was diagnosed) in each imported country was collected from publicly available data sources. Employing a survival analysis model in which the hazard is an inverse function of the effective distance as informed by the airline transportation network data, and using dengue and chikungunya virus transmission data, risks of importation and local transmission were estimated. Results. A total of 78 countries with imported case(s) have been identified, with the arrival time ranging from 1 to 44 weeks since the first ZIKV was identified in Brazil, 2015. Whereas the risk of importation was well explained by the airline transportation network data, the risk of local transmission appeared to be best captured by additionally accounting for the presence of dengue and chikungunya viruses. Discussion. The risk of importation may be high given continued global travel of mildly infected travelers but, considering that the public health concerns over ZIKV infection stems from microcephaly, it is more important to focus on the risk of local and widespread transmission that could involve pregnant women. The predicted risk of local transmission was frequently seen in tropical and subtropical countries with dengue or chikungunya epidemic experience.


INTRODUCTION
to have been caused by ZIKV infection among pregnant women (Staples et al., 2016). The causal relationship between ZIKV infection and microcephaly has yet to be fully established, but an increasing amount of scientific evidence strongly supports the epidemiological link of ZIKV with microcephaly (Mlakar et al., in press;Brasil et al., in press;Nishiura et al., 2016b). Considering the expected impact on neurological manifestations including microcephaly and Guillain-Barré Syndrome (GBS), the World Health Organization issued the Public Health Emergency of International Concern (PHEIC) on 1 February 2016.
ZIKV shares self-limiting clinical signs and symptoms with dengue virus (DENV) and chikungunya virus (CHIKV) (Ioos et al., 2014). Moreover, these viruses are also conveyed by the common mosquito vector, Aedes species. Technically, the infected cases could be observed in any part of the world in the presence of Aedes species, but the risk of transmission in Latin American and Caribbean countries has been anticipated to be high (Musso, Cao-Lormeau & Gubler, 2015;Bogoch et al., 2016). The transmission potential of ZIKV infection has been shown to be comparable to those of DENV and CHIKV (Nishiura et al., 2016a). Considering that there are no specific treatments for these arboviruses, it is of utmost importance to intervene transmission by preventing mosquito biting (e.g., the use of repellent) or controlling Aedes species, most notably by targeting highly invasive species including Aedes aegypti and Aedes albopictus. The control option of Aedes spp. includes source reduction (i.e., removal of water and covering water storage containers), the use of ovitrap, chemical control, bioinsecticides and biological control (e.g., mosquitofish) (Anonymous, 2016), while recent approaches include transgenic and symbiont-based approaches and the use of plant-borne molecules (Adelman & Tu, 2016;Benelli, 2016;Benelli & Mehlhorn, 2016).
To plan for forthcoming potential local transmission at each country level, it is fruitful to understand the actual risk of Zika virus epidemic in a quantitative manner. The present real-time study aimed to identify countries at high risk of importing ZIKV infection and also of local transmission as a part of risk assessment practice and for improved understanding of the predicted risk of the spread, with distinct differentiation of ZIKV importation alone from that of local transmission.

Epidemiological data
Arrival times of ZIKV importation, i.e., the time at which first ZIKV case was imported, in each imported country was collected from publicly available data sources (see Supplemental Information). The criterion to include a country as imported was the presence of a report based on confirmatory diagnosis made either by serological testing or detection of virus RNA. When ZIKV was isolated from mosquitoes or primates other than humans, it was considered that the country had already imported the virus, because the virus isolation implies that the local transmission had taken place. For each imported country, the arrival time was defined as the week in which the importation event had happened. When monthly or yearly data were only available or when a seroepidemiological study identified the local transmission, the median week was used as the week of importation. The latest time at which the importation week was systematically examined was 31 January 2016.
In addition to the importation data, we collected the list of countries with local transmission of ZIKV from the European Centre for Disease Prevention and Control (2016) and Gatherer & Kohl (2016). ECDC resource has specified countries or territories with reported confirmed autochthonous cases of ZIKV infection in the past 9 months. Since the European Centre for Disease Prevention and Control (2016) list misses countries that experienced local transmission in 2014 or earlier, we extracted the list of those earlier countries from Gatherer & Kohl (2016). Again, the virus isolation from mosquitoes or primates other than humans was regarded as a signature of local transmission.
As additional input of prediction models, an open source airline transportation network data were extracted. Using the Global Flights Network (2016) derived from the OpenFlights database as of 10 November 2014 (Contentshare, 2016), we obtained the total number of flight routes between each pair of countries (a total of 230 airports and 4,600 flight routes). Furthermore, the country-specific presence of Aedes species, especially focusing on Aedes albopictus, Aedes aegypti and Aedes africana, was explored (Kraemer et al., 2015;The Walter Reed Biosystematics Unit, USA, 2016). The list of countries with at least one of the abovementioned Aedes species was prepared, because the local human-vectorhuman transmission of ZIKV infection occurs only where the vector species are present. Similarly, we also retrieved the country-specific data regarding DENV and CHIKV epidemic (Centers for Disease Control and Prevention, 2016a;Centers for Disease Control and Prevention, 2016b), because they share major vector species with ZIKV and their epidemics indicate that the vector species have been actively involved in local transmission. Moreover, it has been explicitly discussed that international spread of Zika virus followed the geographic path of CHIKV (Musso, Cao-Lormeau & Gubler, 2015). The presence of Aedes species, DENV and CHIKV were all dealt with as dichotomous variable.

Prediction model
We estimated the risk of importing ZIKV and that of local transmission using a survival analysis model. Let T be a continuous random variable with probability density function of the time from importation in Brazil to importation in country j, f j (t ), and cumulative distribution function F j (t ) = Pr(T < t ). The function F j (t ) describes the probability that ZIKV has been already imported to a country j by time t . The time t = 0 corresponds to the time at which ZIKV infection was first recognized in Brazil and started to rapidly spread across countries (i.e., week 12 of 2015). We parameterize the hazard function λ(t ) using the effective distance, the metric derived from the airline transportation network (Contentshare, 2016). Let {n 1 ,n 2 ,...n l } be the sequence of transit countries at which a traveler starting from country n 1 (with a destination n l ) stops over. The length of path of that travel is l. The effective length of the path, dn 1 n l , is defined by where P ji denotes the conditional probability that an individual that left i moves to j.
Assuming that the number of passengers is identical among all international flights, the transition matrix is calculated as P ji = w ji k w ki , where w ki is the number of direct flights from country i to country k per unit time derived from airline transportation network data (Contentshare, 2016). Finally, the effective distance m j of a country j from the country of origin (i.e., Brazil) is calculated as the minimum of the all possible effective lengths of path that goes from the origin country to the country j. The effective distance has been known as an excellent predictor of the arrival time of SARS (severe acute respiratory syndrome) and influenza pH1N1 2009 (Brockmann & Helbing, 2013). Considering that the effective distance is a critical indicator of the risk of importation, we assume that the hazard function is an inverse of the effective distance, i.e., where k is a constant. Using the hazard function, the abovementioned density function is modeled as Modeling f j (t ) in this way, the mean arrival time of the ZIKV to a country j would be proportional to the effective distance from the origin of spread (i.e., Brazil in our case study).
Using the calculated risk of importation, we subsequently model the country specific risk of local transmission. Let p j be the conditional probability of ZIKV transmission given an importation event in a country j. Since our time scale of examining the risk of importation is much longer than the generation time of ZIKV infection, we ignore the time-lag that is required for observing the transmission cycle, and the risk of the local transmission of ZIKV in a country j is modeled as In order to estimate p j , the conditional probability of experiencing local transmission, we use datasets of the presence of Aedes species, CHIKV and DENV in each country. More specifically, we model p j by employing a logit model where a 0 is an intercept, x ji is the dichotomous variable that describes the presence of Aedes mosquitoes (for i = 1), CHIKV (for i = 2) and DENV (for i = 3) in country j, and β i represents the coefficient for corresponding variable i.

Parameter estimation and risk assessment
We estimate parameters k, a 0 and three β i using a maximum likelihood method. Since we start the clock from week 12, 2015 originated from Brazil (Campos, Bandeira & Sardi, 2015), countries that had imported ZIKV in advance (e.g., Africa and South Pacific) were removed before the implementation. In total, we predict the risks of importation and local transmission among a total of 189 countries. The model was fitted at week 46 since importation in Brazil (and the week 46 corresponds to the latest time at which our systematic survey of importation dates was completed in January 2016). Although five countries (Aruba, Trinidad and Tobago, Marshall Islands, Saint Vincent and the Grenadines, and American Samoa) were specified as countries with local transmission by ECDC data source (European Centre for Disease Prevention and Control , 2016), these countries have experienced importation and confirmation in February 2016, and were thus removed from the list of countries with importation as well as local transmission. The likelihood of importation risk adhered to survival model. That is, among countries that have already imported ZIKV by then, the arrival time t j was used for parameter estimation. Countries that have not imported ZIKV by week 46 were dealt with as the censored observation. In other words, the survival time for those countries is considered to be at least as long as the duration of our study. The likelihood function reads where A represents the set of countries which experienced local transmission of ZIKV by time of observation t m (31 January 2016), B represents countries which imported ZIKV before t m but without local transmission, and C is the set of countries which has not imported by t m . Once the parameter k is estimated, we are able to predict the probability of country j to experience local transmission of ZIKV by time t as We assessed the diagnostic performance of our models in predicting risks of importation and local transmission by employing the receiver operating characteristic (ROC) curve and measuring the area under the curve (AUC) (Greiner, Pfeiffer & Smith, 2000). The predictors of local transmission in Eq. (5) were varied by using 1-3 variables (by examining all possible combinations of Aedes, DENV and CHIKV), and we compared AUC to identify the best predictive model. To avoid serious collinearity among Aedes, DENV and CHIKV, Cohen's kappa, an agreement statistic, was computed to examine correlations between two dichotomous explanatory variables. For each model, the optimal cut-off value of estimated risk was calculated using Youden's index. While model fit was assessed in week 46, the prediction has been made in week 92 that corresponds to the end of 2016. week number of the report has been originally written and we used it as the week of importation. When the exact week was unavailable, the midpoint of the available time window was used as the week of importation. The first day of the year was set to be the first Sunday of January. One year was calculated to be exactly equal to 52 weeks in our analysis. Although three confirmed cases of ZIKV infection were reported in Chile on 2 February 2016, that data was ignored in the present study because the latest time was set at fourth week of 2016. We also did not count a case report from Easter Island, a Chilean island in the southeastern Pacific Ocean, because the island is geographically distant from other parts of Chile and has only one airport (World Health Organization, Western Pacific Region, 2016). In Australia, although one imported case is reported, there were no available sources of information to identify the date of confirmation. Instead, given that the corresponding article was received on 16 January 2013, we consider that the importation event took place in the previous year, 2012, and took the mid-point of the year as the week of importation (Kwong, Druce & Leder, 2013). In Zambia, a cross sectional study was conducted, but there was no available information to identify the date of importation or survey. Thus, the year of the acceptance of the article, 2015, was assumed as the year of importation (Babaniyi et al., 2015). (continued on next page)  Notes. a AUC, area under the curve. The confidence intervals were calculated using Mann-Whitney method (Gengsheng & Hotilovac, 2008). b CI, confidence interval.

Country
first identified ZIKV case in 1947, Uganda (Dick, 1952), but the country was not specified as the origin, because a long time has passed since the emergence and the path of global spread must not have followed a static network. Note that 39 countries which had already imported ZIKV prior to importation event in Brazil were excluded from the analysis.  Figure 2 shows the global distribution of the risk of importation. High risks of importation are identified in South America and the western part of European countries. Among the top 30 countries predicted at high risk of ZIKV importation in the end of 2016 (Fig. 3), 18 countries (60%) had already imported ZIKV before week 46. Table 2 shows the fitting results of various models to predict the risk of local transmission. Comparing predictive performance by AUC among models with all possible combinations of three explanatory variables, the model with CHIKV and DENV appeared to yield the greatest value of AUC (0.90 (95% CI [0.60-1.00])). Figures 1B and 1D show the distribution of the estimated risk and ROC curve, respectively. Using the best model, sensitivity was as high as 96.4% (95% CI [89.6-100.0]), while the specificity was estimated at 67.7% (95% CI [60.5-74.9]). The second best model was the one with CHIKV only.   Table 3 shows parameter estimates of the best fitted model. Compared with the absence of Chikungunya virus, the risk of local transmission in the presence of this virus was shown to be indicative of 22.9 times (95% CI [3.3-238.3]) higher. The presence of dengue was not statistically significant, but the adjusted odds ratio was estimated to be 7.7 (95% CI [1.0-73.6]). Table 4 shows correlations between two dichotomous variables as measured by Cohen's kappa. No particular correlation that could lead to multi-collinearity (e.g., kappa > 0.60) was identified. Figure 2B shows the country-specific global distribution of local transmission using the best model. The high risk of local transmission is seen mainly among countries in The importation risk of ZIKV by week 92 is colored by intensity (0-20%, 20-40%, 40-60%, 60-80%). The origin country, Brazil, is colored in grey. Other additional countries colored in grey were excluded, because they experienced importation of ZIKV infection prior to the event in Brazil (week 12, 2015). (B) The local transmission risk of ZIKV infection by week 92 accounting for dengue and chikungunya epidemic data. The local transmission risk of ZIKV is colored by intensity (0-15%, 15-30%, 30-45%, 45-60%). The origin country, Brazil, is colored in grey. Other additional countries colored in grey were excluded, because they experienced importation of ZIKV infection prior to the event in Brazil (week 12, 2015).  tropical and subtropical areas. Figure 3B show the top 30 countries with high risk of local transmission by week 92, using Brazil as the origin of spread. Among the total, 19 countries (63.3%) were predicted to have already allowed local transmission before week 46.

DISCUSSION
The present study estimated country-specific risk of importation and local transmission of ZIKV infection using a simple statistical model. As reported elsewhere, ZIKV infection was often internationally spread by mildly infected travellers (Tappe et al., 2014;Kutsuna et al., 2014;Korhonen et al., 2016). Potentially high importation risk in many temperate countries has motivated us to explore the risks of importation and local transmission fuelled by travellers. Our model is not as sophisticated as mapping precise risk of transmission using seasonal population dynamics of Aedes species and temperature/climatological data (Bogoch et al., 2016;Nah et al., in press), but the findings at country levels from our study are broadly consistent with what has been briefly described at finer scales (Bogoch et al., 2016). Without using finer scale spatial data (e.g., ecological data on Geographic Information Systems), our approach has crudely and clearly distinguished the risk of local transmission from importation risk at country levels using a more tractable approach. We have shown that the predicted risk of local transmission was frequently seen in tropical and subtropical countries with DENV or CHIKV epidemic experience, while the risk of importation was more scattered around the world. The diagnostic performance of risk model for local transmission was well supported by AUC value of 0.90. Our study contributes to the risk assessment practice at each country level using estimated risks of importation and local transmission expressed as probabilities. We have shown that the risk of importation may be high in several countries given continued global travel of mildly infected travellers. However, considering that the public health concerns over ZIKV infection stems from the presence of microcephaly, it is more important to focus on the risk of local and widespread transmission that could involve pregnant women. Compared with the risk of importation, the risk of local transmission was particularly highlighted in countries with DENV and CHIKV epidemic experience. Considering that the Olympic Game in Brazil 2016 will elevate the risk of ZIKV to spread over a wider spatial extent (Petersen et al., 2016), the distinction between importation and local transmission will be even more important than it has been recognized. Whereas our risk model of local transmission relied on the presence of transmission by DENV or CHIKV, the use of such a dataset should be deemed only as a proxy. The sensitivity of our best model was 96.4%, and thus there is limited case in which the absence of DENV or CHIKV lead to the local transmission of ZIKV. However, the specificity was limited to 67.7%, indicating that there would be a number of countries in which ZIKV transmission has yet to occur even in the presence of DENV and CHIKV. Our simple model must have missed important additional predictors of local transmission in this regard, and such predictors have to be sought in the future studies. Moreover, to further improve the prediction, ecological data of vector behaviour at finer spatial scale must be incorporated with validation using epidemic data.
A few limitations must be discussed. First, ascertainment bias cannot be ignored for ZIKV infection with substantial fraction of asymptomatic and mild infections. The empirical data in Table 1 must have missed several importation events. Presented risks of importation and local transmission were thus probably underestimated on a whole. Second, the network data we used were static, and the airline transportation data represented only the number of flight routes (i.e., not the number of passengers). Nevertheless, assuming that the geographic patterns of travel have been maintained and only the travel volume has changed over short period of time, our simplistic approach is still justifiable to predict the global spread from 2015. Third, country-specific modelling exercise suffers from heterogeneity within each single country. For instance, the United States, China and Australia have both high and low risk areas of transmission due to vast land with different climatological areas.
Finally, we estimated the importation risk of ZIKV in a country as a risk of importing ZIKV from Brazil to the country. Therefore, our model cannot capture the risk of importing ZIKV from other endemic area such as South Pacific.
Despite a clear need to improve predictions in the future, the present study successfully devised a simple global risk prediction of importation and local transmission. Countries with DENV and CHIKV epidemic experience are likely to be at particular high risk, and such countries should be prepared for vector control measures such as avoiding daytime biting and the use of mosquito repellent. To further improve model predictions, it is essential to have laboratory capacity built up in every single country at risk of local transmission.

CONCLUSIONS
Risks of importation and local transmission of Zika virus infection were estimated, analyzing epidemiological, ecological and mobility data. Whereas the risk of importation was well explained by the airline transportation network data, the risk of local transmission appeared to be best captured by additionally accounting for the presence of an epidemic of dengue and chikungunya viruses.