Traveling waves in a free boundary problem for the spread of ecosystem engineers

1.   Introduction

Zika virus (ZIKV) is a member of the family Flaviviridae, and is the viral cause of Zika fever, which, when symptomatic, involves only mild symptoms including fever, red eyes, joint pain, headache, and a maculopapular rash; Zika resembles a mild form of dengue fever ^[1]. The virus is named for the Zika Forest, in Uganda, where the virus was first isolated in 1947 from the serum of a rhesus monkey ^[2]. Aedes species, especially Aedes aegypti and Aedes albopictus, serve as vectors for ZIKV, allowing transmission between Aedes and humans. The transmission is autochthonous in many South American countries ^[3,4,5]. Historically, human ZIKV infection is less common than dengue fever, but international spread accelerated the global epidemic, and epidemics were observed in 2007 in Micronesia, in 201314 in French Polynesia, and, since 2015, in Brazil and across the world ^[6,7,8,9]. Following the detection of a large number of cases in Brazil in late 2015, an excessive number of microcephaly cases was reported approximately 30 weeks later in 2016, in northeastern Brazil, which attracted global attention to this virus ^[10,11,12]. Although ZIKV infections tend to be mild, infected pregnant women are at risk of delivering babies with microcephaly or other congenital disorders. Additionally, the infection can result in Guillain-Barre syndrome in adults ^[13,14].

The 2015 global epidemic was initially concentrated in northeastern Brazil; however, the infection rapidly spread throughout Latin America and the Caribbean and subsequently to many countries around the world ^[15,16,17]. More than 45 American countries reported ZIKV importation ^[18,19], and, of these, all countries except Canada and the US reported ZIKV importation after the Brazil epidemic. While many European countries have reported ZIKV importation from Brazil since the first imported case was reported in 2015, only a few African countries reported importation after the 2015 Brazil epidemic ^[20]. There are two possible explanations for the timing and the presence of imported cases in a given country, including why some countries did not experience importation. First, due to a countrys effective distance from Brazil or any other epidemic locations, importation may not actually occur. Second, due to limited laboratory capacity, imported cases could have gone undetected, resulting in reporting delays and underreporting. The status and quality of ZIKV diagnostic capacity are known to be highly variable, even among European countries ^[21], to say nothing of the limitations healthcare workers navigate in low-income countries.

Previous mathematical models were developed to describe vector-borne disease transmission dynamics ^[22,23], and, after the 2015 Brazil epidemic, mathematical models were applied to ZIKV to describe both local ^{[24,25,26,27]} and international transmission ^[28,29,30]. Both deterministic and stochastic models for ZIKV were formulated to fit data from Colombia, El Salvador, and Suriname in ^[5]. Massad et.al ^[31] identified high-risk countries in Europe by using mathematical models to estimate the importation risk from Brazil according to important factors, such as travel volume. In fact, there is an associated question with respect to the predictability of importation time for ZIKV: real-time forecasting for emerging infectious disease epidemics can be achieved by employing meta-population models ^[32,33], and, as a possible short cut, the so-called "effective distance, " which can be computed from airline network data, has been proposed ^[34,35]. Use of effective distance indicates that the relative arrival time for imported cases is determined by the path length and adjacency matrix from the origin of the epidemic to the destination. However, in the case of ZIKV, effective distance alone is not necessarily an adequate predictor. Figure 1A shows arrival time as a function of the effective distance from Brazil. While a positive correlation is certainly observed, the resulting importation times were highly variable. Such variation was not identified during severe acute respiratory syndrome (SARS) epidemic and the pandemic of H1N1 influenza in 2009 ^[34]. Figure 1B classifies countries into three different groups by gross domestic product (GDP) per capita. It seems that the variation in arrival time and deviation from the linear predictor in Figure 1 were partly explainable by GDP or other variables that influence the underlying diagnostic capacity in a given country.

To explore the mechanisms of ZIKV importation, we argued that it is crucial to account for possible inter-country variations in testing, diagnostics, and reporting ZIKV. In this study, we statistically estimated the actual arrival time of ZIKV importation in each country using airline network data, while adjusting for reporting delays. Describing the data generating process of imported cases, we mathematically explored mechanisms behind observed variations in ZIKV arrival time.

2.   Materials and method

2.1.   Epidemiological data

The time from Brazil to "reported" importation of ZIKV is defined as the time interval from the first notification of ZIKV infection in Brazil to the time at which each country reported the first ZIKV-infected case. The observed (reported) arrival time was collected from publicly available data sources for a total of 219 countries. Following a recent similar study that was conducted elsewhere ^[36], we updated our arrival time data. As of September 25, 2018,110 countries reported ZIKV importation (see below for non-imported countries). Of these, we excluded 40 countries that experienced importation before the 2015 Brazil epidemic. In addition to arrival time data, airline transportation data were obtained from the OpenFlights database ^[37], which uses the Global Flights Network (2016) that includes 230 airports and 4600 flight routes. We selected datasets that can be considered to reflect the diagnostic and reporting capacity of affected countries ^[38,39]. Specifically, we collected datasets of GDP per capita and several types of health expenditure metrics, including government expenditure on health, per capita government expenditure on health, and private expenditure on health. Moreover, we categorized countries into geographic regions using two different classification systems: first, by World Health Organization (WHO) region, and, second, by continent. Other variables that we examined included religion (Christian, Muslim, and others) and language (English, Spanish, and others). See Otsuki and Nishiura ^[40] for data collection methods for those explanatory variables.

2.2.   Importation risk

To determine actual ZIKV arrival times, we constructed a mathematical model that convolutes two functions (Figure 2). For the first function, we modeled the time from Brazil to actual importation, in which we exploited the relationship between arrival time and effective distance. We derived the effective distance, which was created by Brockmann and Helbing ^[34], from the abovementioned flight network, and we calculated it as the minimum of the summation of all possible path lengths and the logarithm of the product of transition probabilities from the origin country (i.e., Brazil) to each destination country $i$. A key feature of effective distance is that, as long as we handle emerging infectious diseases, it has shown a very strong correlation with arrival time, demonstrating itself to be an excellent predictor of arrival time. The validity of the predictive performance of effective distance has been described elsewhere ^[34], and has been used in studies exploring data from the 2003 SARS epidemic and the influenza A(H1N1-2009) pandemic, and it has been widely used to analyze a variety of global infectious disease pandemics ^[36,40,41]. Let $m_i$ be the effective distance from Brazil to country $i = 1, 2\cdots, 170$. Due to the linear relationship between arrival time and effective distance, the hazard function of importation for country $i$ is modeled as an inverse of the effective distance:

where $k$ is a constant parameter (to be estimated below). The country-specific hazard function yields the probability density function for ZIKV importation from Brazil to a country $i$ at time $t$, written as

where $t_0$ is the time at which the global epidemic started at the origin, i.e., Brazil. In this study, we set $t_0$ as 48 weeks prior to the date that Brazil reported their first ZIKV importation, referring to Faria et al.s molecular clock analysis ^[18].

2.3.   Time from importation to reporting

For the second function, we modeled the time from importation to reporting, accounting for reporting delays that differ by country (Figure 2). The country-specific mean reporting delay $\mu_{i}$ was modeled with a linear predictor using a combination of two variables: (1) financial capacity for diagnosis and reporting ($x_{1}$), which was either modeled by GDP per capita or health expenditure (per capita government expenditure on health), and (2) region $x_{2}$ according to the WHO classification. For $x_{1}$, we divided 170 countries into 3 groups i.e., $G_1 = \left\lbrace A, B, and \ C \right\rbrace$ based on tertiles of high, intermediate, and low, respectively. Three-group discretization was adopted because we did not observe GDP per capita as having precise predictor performance when handled as a continuous variable. Additionally, we classified countries into three groups based on WHO region for the variable $x_2$: region of the Americas (AMR), African region (AFR), and countries which belong to neither AMR nor AFR ($Others$), i.e., $G_2 = \left\lbrace AMR, AFR, and \ Others \right\rbrace$. Again, we employed discrete grouping to approximately capture the observed patterns, including the tendency that countries in the AMR reported ZIKV early, while those in the AFR reported relatively late. The linear predictor of $\mu_{i}$ for country $i$ is

where $\beta$s are coefficients to be estimated, and $x_{y, G_y(j)}$ represents the dummy variable $x_y$ that belongs to group $j$, defined as

That is, for $x_{1, G_1(j)}$, $j = 1$ represents high income group, $j = 2$ represents intermediate, and $j = 3$ low income. For $x_{2, G_2(j)}$, $j = 1$ represents countries that belong to AMR, $j = 2$ indicates countries that belong to AFR, and $j = 3$ represents other countries. The country grouping $j$ is determined according to country $i$. Let $J_y (i)$ be the group of variables $x_y$ for country $i$. Then, $j$ corresponds to $J_y (i)$. It shoud be noted that, while AMR is a mix of North America and other countries, the USA and Canada imported ZIKV case prior to Brazil and were excluded from our analysis. We calculated the correlations between $x_{1, G_1}$ and $x_{2, G_2}$ using the Kendall rank correlation coefficient. If significant correlations $p < 0.05$ were identified, we attempted to address the interaction using the interaction term from the linear predictor. However, if the interaction was not adjusted by the interaction-term alone, we removed the candidate model with two correlating variables prior to model comparisons.

For a country $i$, groups of two variables can be described as $g_i = \lbrace G_1(J_1(i)), G_2(J_2(i)) \rbrace$. We assumed that the reporting delay followed a gamma distribution. Supposing that country i belongs to group ${g_i}$, the reporting delay was assumed to be the result of an independent random sampling from the probability density function $g(t; \mu_{i}, \sigma_{i, g_i}^{2})$ with a country-specific mean $\mu_{i}$ and variance $\sigma_{i, g_i}^{2}$. To obtain group specific variance $\sigma_{i, g_i}^{2}$, we first estimated the coefficient of variation in two different wayseither using a group-specific standard deviation $\sigma_{g_i}$ (i.e., $cv_{g_i} = \frac{\sigma_{g_i}}{\mu_{g_i}}$) or a constant standard deviation over all the groups $\sigma$ (i.e., $cv_{g_i}^{'} = \frac{\sigma}{\mu_{g_i}}$). Consequently, group-specific variance for gamma distribution obeyed the formula $\sigma_{i, {g_i}}^{2} = cv_{g_i}\mu_i$ or $\sigma_{i, g_i}^{2} = cv_{g_i}^{'}\mu_i$.

2.4.   Likelihood functions

A total of 219 countries were divided into three different groups:

(ⅰ) 56 countries that reported ZIKV importation after the 2015 Brazil epidemic; the country $i$ in group $g_i$ is defined by $i \in C^{RP}_{g_i}$ where RP stands for "reported."

(ⅱ) 89 countries that reported no ZIKV importations by the end of 2017; the country $i$ in group $g_i$ is defined by $i \in C^{NRP}_{g_i}$ where NRP stands for "never reported."

(ⅲ) 74 excluded countries, including 40 countries that reported ZIKV importation before the 2015 Brazil epidemic and 34 countries for which there are no available data for either GDP per capita or health expenditure.

For the countries in group (ⅰ), the probability density function of the reporting time $h_{i} (t, t_0;k, \mu_i, \sigma_{i, g_i}^{2})$ obeys the formula:

For simplicity, let $h_{i}(t)$ be the simpler notation of $h_i(t, t_0;k, \mu_i, \sigma_{i, g_i}^{2})$. The likelihood for observing $t_i$ in this group is

where $t_i$ is the time of reported ZIKV importation for each country $i$. For the countries in group (ⅱ), the probability density function of the reporting time $w_{i} (t, t_0;k, \mu_i, \sigma_{i, g_i}^{2})$ is calculated as

As was done for group (ⅱ), let $w_i (t)$ be the simpler notation of $w_{i}(t, t_0;k, \mu_i, \sigma_{i, g_i}^{2})$. The likelihood function for observing censoring for this group is

where $t_c$ is the substituted time of reported ZIKV importation into country $i$, which we assumed was long enough to represent no reported cases. Finally the full likelihood function is

Profile-likelihood based confidence intervals were computed to obtain the 95% confidence interval (CI). We computed the second order Akaike information criterion (AICc) for model comparison.

2.5.   Ethical considerations

Herein, we only analyzed publicly available data. As such, the datasets were de-identified and fully anonymized in advance, and the analysis of publicly available data without identifying information does not require ethical approval.

2.6.   Data sharing policy

This study fully relies on published data, essential components of which we made available as supplementary material.

3.   Results

Models with variable combinations of explanatory parameters for describing the time from importation to reporting were optimized using the same datasets, yielding a total of 56 imported and 89 non-imported countries during the study period. We compared 10 models and selected the best according to which generated the lowest AICc value. Among the financial variables that influenced diagnostic and reporting capacities, GDP per capita and health expenditure were correlated, indicating that they could not jointly explain the reporting delay. Accordingly, when either GDP per capita or health expenditure was included, the other variable was excluded from the final model. Table 1 compares the goodness-of-fit of those 10 different models. Of these, the model with GDP and WHO region yielded the lowest AICc value of 681.7 (Table 1). We compared the mean reporting delay according to GDP per capita tertiles (Figure 3) using the best model (i.e. Model 1 in Table 1). We found that countries with the highest GDPs yielded the shortest reporting delays. It can be interpreted that countries with relatively high GDPs were more likely to report ZIKV infection earlier than countries with lower GDPs. The median reporting delays for first, second, and third GDP tertiles were 12 weeks (95% CI: 1222), 30 weeks (95% CI: 2030), and 35 weeks (95% CI: 2435), respectively. On average, approximately 4 months after the countries in group A reported a first case, countries in group B reported ZIKV importation, and countries in group C subsequently reported ZIKV importation approximately 1 month later. Figure 4 shows the linear relationship between the effective distance and the time from actual importation in Brazil to actual arrival in each country, according to GDP tertile grouping. As the effective distance became longer, the time required for importation increased. Moreover, comparing the time since actual arrival by GDP groups, the countries in the highest GDP tertile (A) experienced a longer time to actual arrival, even up to 75 weeks, since initial importation in Brazil. The majority of countries in the lowest tertile group experienced importation within 25 weeks of actual importation in Brazil.

4.   Discussion and conclusion

We statistically estimated ZIKV arrival time around the world after it first appeared in Brazil in 2015. Taking a similar approach to our previous study that we conducted in real time ^[36], we modeled inter-country infection spread using effective distance, but, as Figure 1 indicates, there was substantial variation among observed importation times, and we determined that effective distance alone was not an adequate predictor of ZIKV arrival time. To better decipher the importation mechanisms, we also accounted for reporting delay, which we regressed by plausible indicators of country-specific laboratory capacity and reporting, including GDP per capita, health expenditure per capita, and geographic regions. We found that high GDP is a good predictor of short reporting delays. Additionally, reporting delay was dependent on WHO geographic region. ZIKV infection is generally mild and, without substantial laboratory capacity, cases can be underestimated. Herein, we successfully highlighted this feature and identified several variables as important to the data generating process of time from Brazil to reported importation around the world.

There are two major findings from our regression models of reporting delay, which yielded smaller AICs than models using a constant to explain delay. First, the reporting delays were shorter in countries with higher GDP per capita. We chose GDP per capita to reflect the capacity of laboratory testing and surveillance. For similar reasons, we also evaluated health expenditure per capita. Because GDP and health expenditure were correlated, the weaker predictorhealth expenditurewas not included in the final model. Because of this finding, the estimated actual arrival time from the actual start of the Brazil epidemic was longest in countries in the highest tertile of GDP per capita. To the best of our knowledge, this study is the first to identify and decipher country-dependent mechanisms behind reporting delays.

Second, we found that South American countries had shorter reporting delays than other countries. Owing to their geographic proximity to Brazil and an elevated awareness of the virus among South American countries, this is an intuitive finding. We found longer reporting delays among the African countries. This is also in line with our expectation that the laboratory capacity in the African region may be lower relative to other countries. Moreover, owing to the mild nature of ZIKV infection, which is historically endemic in African countries ^[13], it is probable that there is a high rate of underreporting in this region.

As shown in Figure 4, our study does not fully explain the global variation in arrival time. Rather, owing to the discrete grouping of GDP per capita and other variables, the estimated variation in actual arrival time was escalated compared with the observed (reported) arrival time. As a response to country-specific variations, herein, we showed that laboratory capacity for testing and surveillance, as well as elevated illness awareness and geographic proximity to the origin, are likely to drive diagnosis and reporting rates. Rather than reducing the arrival-time variance, our contribution with this study was to identify the effect of reporting delays on country-specific observed arrival times of the ZIKV epidemic. A critical lesson learned from modeling the 201516 ZIKV pandemic is that the use of effective distance alone is not sufficiently precise to capture the observed patterns of country-specific arrival times; rather, we need to additionally account for case ascertainment, laboratory capacity, and virus awareness to address possible variations in reporting delays and laboratory coverage.

Several technical limitations of our study should be noted. First, even though our linear regression model explored reporting delay by country, ZIKV epidemics include asymptomatic infections; thus, it is possible that several importation events were missed. Second, even though country-specific variation was regressed, our approach adopted discrete classification of GDP per capita into three groups, but this approach could not fully explain variations in observed arrival times. Additional mechanisms need to be identified via country-specific analyses. Third, we used static network data to capture human mobility patterns, but it is possible that global spread in 2016 was enhanced due to the Olympic Games, which were held in Rio de Janeiro that year ^[42]. Lastly, we estimated the ZIKV importation risk in each country as a risk of importation from Brazil. Importations from other areas that experienced a contemporaneous ZIKV epidemic (e.g., South Pacific) were ignored.

Despite these limitations, we argue that, with this study, we have successfully elucidated key importation mechanisms, including diffusion by human migration and reporting delay. The time required for ZIKV to arrive and be reported in a country was determined not only by global airline travel patterns, but also by other factors, most notably heterogeneous laboratory and reporting capacities.

Acknowledgments

HN received funding support from the Japan Agency for Medical Research and Development (JP18fk0108050); Japan Society for the Promotion of Science KAKENHI (Grant Numbers 16KT0130, 17H04701, 17H05808 and 18H04895); Health and Labour Sciences Research Grant (H28-AIDS-General-001); the Inamori Foundation; the Telecommunication Advancement Foundation; and the Japan Science and Technology Agency (JST) CREST program (JPMJCR1413). HL received financial support from the Japan Society for the Promotion of Science, Program for Advancing Strategic International Networks to Accelerate the Circulation of Talented Researchers; and the Japan Society for the Promotion of Science KAKENHI (Grant number 18H06385). The funders had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. We also thank Ashley Hazel, PhD, from Edanz Group (www.edanzediting.com/ac) for editing a draft of this manuscript.

Conflict of interest

The authors declare that they have no conflict of interest.