key: cord-1013248-cvuvont8
authors: Wang, Lin; Wu, Joseph T.
title: Characterizing the dynamics underlying global spread of epidemics
date: 2018-01-15
journal: Nat Commun
DOI: 10.1038/s41467-017-02344-z
sha: fed39ba33ece570bf3f7e8879cf448bf93cd069e
doc_id: 1013248
cord_uid: cvuvont8

Over the past few decades, global metapopulation epidemic simulations built with worldwide air-transportation data have been the main tool for studying how epidemics spread from the origin to other parts of the world (e.g., for pandemic influenza, SARS, and Ebola). However, it remains unclear how disease epidemiology and the air-transportation network structure determine epidemic arrivals for different populations around the globe. Here, we fill this knowledge gap by developing and validating an analytical framework that requires only basic analytics from stochastic processes. We apply this framework retrospectively to the 2009 influenza pandemic and 2014 Ebola epidemic to show that key epidemic parameters could be robustly estimated in real-time from public data on local and global spread at very low computational cost. Our framework not only elucidates the dynamics underlying global spread of epidemics but also advances our capability in nowcasting and forecasting epidemics.

S ince the 1980s, metapopulation epidemic models built with worldwide air-transportation network (WAN) data have been the main tool for studying global spread of epidemics, such as pandemic influenza [1] [2] [3] [4] , SARS 5, 6 , MERS-CoV 7 , Ebola 8 , and Zika 9, 10 . The complexity of these models has substantially grown over the past few decades, advancing from 55 populations in the Rvachev-Longini model in 1985 1 to more than 3500 populations in the state-of-the-art simulator GLEAM powered by supercomputer 11, 12 . Despite the long history and widespread use of these models [13] [14] [15] [16] [17] , most studies on global spread of epidemics have relied on computationally intensive simulations that provide limited epidemiologic insights, whereas an analytical understanding of the underlying epidemic dynamics has only been partially elucidated in recent years [18] [19] [20] . Here, we build on these recent advancements and develop a novel framework for analytically characterizing how epidemic arrivals for different populations around the world depend on the epidemiologic parameters and structure of the WAN. We first validate this framework using global epidemic simulations. We then illustrate its potential to enhance our ability to nowcast and forecast epidemics by applying it retrospectively to the 2009 influenza A/H1N1 pandemic and the 2014 West African Ebola epidemic in Liberia.

Major assumptions in the framework. Throughout this paper, we consider only global spread of epidemics with relatively fast timescales in which epidemics in each population peak within 300 days after establishment (e.g., pandemic influenza, MERS, Ebola) such that changes in demographics (e.g., births, aging) is negligible. In metapopulation epidemic models, populations (e.g., cities) around the world are connected through the travel of individuals via the WAN (see WAN metapopulation epidemic model in Methods for details). We designate population i as the epidemic origin which is seeded with s i infections at time 0. For any given population j, we denote its population size by N j and initial epidemic growth rate by λ j . If populations j and k are directly connected in the WAN, the mobility rate from population j to k is defined as w jk = F jk /N j , where F jk is the direct airtraffic (passengers per day) and w jk ranges mostly between 10 -6 and 10 -3 per day in the WAN ( Supplementary Fig. 1 ). Let T n ij be the time at which population j receives its nth imported infection. The epidemic arrival time (EAT) for population j is defined as T 1 ij . Our framework is built upon the following assumption 19, 21 .

Assumption 1: Suppose populations j and k are directly connected in the WAN and only population j is infected. Exportation of infections from population j to k is a nonhomogeneous Poisson process (NPP) 22 with intensity function w jk I j (t) where I j (t) is the disease prevalence (number of infectives) in population j at time t (see Details on assumption 1 in Methods for details).

Supplementary Figure 2 shows that assumption 1 is very accurate for a wide range of plausible epidemic scenarios. We will show that the dynamics of global spread is largely analytically tractable because the following assumption is also accurate across these same scenarios.

Assumption 2: After the epidemic has established in a given population j, the first few exportations occur while disease prevalence is still growing exponentially, i.e., I j t ð Þ ¼ s j exp λ j t À Á . To this end, we progressively build up our framework by characterizing the probability distribution of EATs for all populations in three metapopulation networks with increasingly complex structure: (i) The two-population network which has the simplest metapopulation structure; (ii) the shortest-path-tree of the WAN (WAN-SPT hereafter) which is the dominant subnetwork driving global spread of epidemics as described by the seminal study by Brockmann and Helbing 20 ; and (iii) the WAN.

The two-population network. In the two-population network, the origin population i is only connected to population j. Under assumption 2, the probability density function (pdf) of T n ij can be expressed in closed-form:

where α ij = w ij s i , which we term adjusted mobility rate. Figure 1 shows that if n is smaller than 10, Eq. 1 is accurate across a wide range of realistic scenarios (e.g., the percent error in expected EAT is uniformly below 2%), which correspond to epidemics ranging from pandemic influenza (with doubling time around 4 −5 days) to Ebola (with doubling time longer than 20 days). This result leads to the following corollaries for the WAN-SPT and WAN analysis: (i) Exportation of the first n infections is essentially an NPP with intensity function α ij exp λ i t ð Þ; and (ii) the expected time of the nth exportation is given by

which corresponds to the (n−1)th exportation for an epidemic that starts at time T 1 ij with seed size s i exp λ i T 1 ij . These analytics can be used to formulate closed-form likelihood functions for inferring parameters from disease surveillance and global spread data (see Methods).

For the WAN-SPT and WAN analysis, we use the worldwide passenger booking data from the Official Airline Guide (OAG) and the Gridded Population of the World data set (Version 4) from the NASA Socioeconomic Data and Applications Center at Columbia University to build a stochastic metapopulation global epidemic simulator with 2309 populations and 54,106 connections (see The global epidemic simulator in Methods). This simulator is similar to GLEAM but without the effect of local commuting which has negligible impact on global spread 23 . Brockmann and Helbing 20 suggested that global spread of epidemics is primarily driven by the WAN-SPT subnetwork in which each population is connected to the epidemic origin via only one path. We will show that for each population k in the WAN-SPT, the time at which the nth importation occurs, namely T n ik , can be well characterized by f n (t|λ,α) (Eq. 1), where λ and α are specifically parameterized to account for the hub-effect and continuous seeding (explained in the next two sections and Fig. 2a, b) . This provides a profound insight: the epidemic arrival process for each population k in the WAN-SPT can be approximated as an NPP with intensity function in the form of α exp λt ð Þ. Hub-effect: Hubs are populations that have direct connections to many populations in the WAN, e.g., Hong Kong, Beijing, New York. If the epidemic origin is a hub, the growth of local disease prevalence can be substantially reduced if a significant proportion of infections travel outward as the epidemic unfolds 18 . Let D i,c be the set of populations that are c degrees of separation from the epidemic origin in the WAN-SPT 24 . From the perspective of the importation process for a given population j ∈ D i,1 (i.e., directly connected to the epidemic origin), the prevalence in population i grows exponentially at rate λ ij = λ i − ∑ k ≠ j w ik (see Fig. 2a and The WAN-SPT analysis in Methods for details). Figure 2d and Supplementary Fig. 3 show that with this hub-effect adjustment, f n (t|λ ij ,α ij ) accurately characterizes the probability distribution of T n ij for all populations in D i,1 (e.g., the percent error in expected EAT is uniformly below 4%). Continuous seeding: Unlike the epidemic origin which has a single seeding event at time 0, all the other populations in the WAN-SPT are continuously seeded by infections coming from their upstream populations. Suppose population k ∈ D i,2 is connected to the epidemic origin via population j along the path ψ: i → j → k. After the epidemic has arrived at population j at time T 1 ij , population i continues to export infections to population j before the epidemic arrives at population k at time T 1 ik (illustrated in Fig. 2b ). Under assumption 2, each imported infection in population j (arriving at times T 1 ij , T 2 ij , …) spawns an infection tree that grows exponentially at the hub-adjusted rate λ jk . Therefore, the prevalence in population j, namely I j (t), is simply the sum of the prevalence for all these infection trees. As such, assumption 1 warrants that the exportation of infections from population j to k is an NPP with intensity function w jk I j (t), which is itself a stochastic process because of its dependence on the random variables T 1 ij , T 2 ij , … (see The WAN-SPT analysis in Methods). We conjecture that this highly complex stochastic process can be greatly simplified with little loss of accuracy by assuming that conditional on T 1 ij (the EAT for population j),

2). In other words, the major source of stochasticity in I j (t) comes from T 1 ij , which is characterized by f 1 (t|λ ij ,α ij ) (Eq. 1). Figure 2e and Supplementary  Fig. 4 show that our conjecture is valid. The resulting approximate pdf of T n ik is accurate for D i,2 populations for all realistic epidemic scenarios. Furthermore, this pdf can in turn be well approximated with f 1 (t|λ ψ ,α ψ ), where λ ψ and α ψ are obtained by minimizing the relative entropy 25 

2.5% 5% 0% 5% 10% Fig. 1 Validating the framework in the two-population model. a-c The analytical (red dashed lines) and simulated (gray lines) pdf of T 1 ij , T 5 ij , and T 10 ij for an exemplary influenza pandemic, where the mean generation time T g is 3.5 days and the initial epidemic doubling time t d is 5 days. The epidemic origin has a population size of 7 million and is seeded with 10 infections at time 0. The mobility rate w ij is 5 × 10 -6 , 5 × 10 -5 , and 5 × 10 -4 per day, which span the realistic range for populations with 1-10 million people in the WAN ( Supplementary Fig. 1 ). d-f Quantile-quantile (Q-Q) plots for the analytical and simulated quantiles of T 1 ij , T 5 ij , and T 10 ij across 100 epidemic scenarios randomly generated from the following parameter space using Latin-hypercube sampling: doubling time t d and generation time T g both between 3 and 30 days, seed size s i between 1 and 100. Each epidemic scenario is coupled with a set of network parameters randomly generated with mobility rate w ij between 10 -6 and 10 -3 and population size N i between 0.1 and 10 million. Simulated quantiles in each scenario are compiled using 10,000 stochastic realizations. In the Q-Q plots, deviations from the diagonal indicate discrepancies between the analytical and simulated quantiles. Data points are colored in blue if the number of exportations is n or above with probability 1, and yellow otherwise. Insets show the corresponding histograms of percent error in E½T n ij origin to any population k ∈ D i,2 can be regarded as a twopopulation model, in which the adjusted mobility rate is α ψ and the epidemic in the origin grows exponentially at rate λ ψ . We term this procedure path reduction. By induction, we can recursively apply path reduction to characterize the EATs with comparable accuracy for all populations in D i,3 , D i,4 , etc. Figure 2f and Supplementary Fig. 5 verify this claim (e.g., the percent error in expected EAT is uniformly below 4%).

WAN. The accuracy of our framework for the WAN-SPT implies that for each (acyclic) path ψ connecting an arbitrary population k to the epidemic origin, the epidemic arrival process for population k along this path can be approximated as an NPP with intensity function α ψ exp λ ψ t À Á . In the entire WAN, each population may be connected to the epidemic origin via multiple paths (hence the difference in EATs between WAN-SPT and WAN, as shown in Supplementary Fig. 6 ), some of which may intersect and are therefore dependent. We conjecture that the dependence among such paths is sufficiently weak such that the overall epidemic arrival process for any population k is well approximated by the superposition of the NPPs 22 that correspond to these pseudo-independent paths. That is, if Ψ ik is the set of all acyclic paths connecting population k to the epidemic origin, the epidemic arrival process for population k can be well approximated by an NPP with intensity function P ψ2Ψ ik α ψ exp λ ψ t À Á . Figure 3 and Supplementary Fig. 7 show that our framework is accurate for all populations and epidemic scenarios.

Public health applications. Our framework provides both analytical and computational advancements for studying global spread of epidemics. First, not only can our framework be easily used to forecast EATs for all populations in the WAN, but it also analytically elucidates the dependence of EATs on the epidemiologic parameters (growth rate and seed size) and the network properties of the WAN (air-traffic volume and connectivity). Second, our framework provides closed-form probability distributions (Eq. 1) to support likelihood-based inference of key epidemiologic parameters from surveillance data on global and local spread. We exemplify the public health applications of our framework by retrospectively applying it to the 2009 influenza pandemic and the 2014 Ebola epidemic as follows.

In our first case study, we infer the transmissibility of the 2009 pandemic influenza A/H1N1 virus in Greater Mexico City following the formulation in Balcan et al. 26 Shortly after the pandemic influenza A/H1N1 virus was first detected in the USA and Mexico in April 2009, many countries enhanced their surveillance to monitor importations of pandemic infections. As such, data on EATs for these countries were deemed more reliable than epidemic curve data, which are typically confounded by reporting behavior and surveillance capacity [26] [27] [28] . Using GLEAM simulations powered by supercomputers to perform maximumlikelihood analyses of EATs for 12 countries seeded by Mexico, Balcan et al. 26 29 and other studies 27, 28, 30, 31 . The reduction in computational complexity and requirement provided by our framework translates into substantial improvement for timeliness and efficiency in situational awareness.

In our second case study, we analyze the 2014 West African Ebola epidemic in Montserrado and Margibi, Liberia (Montserrado henceforth for brevity). Specifically, we apply our framework to retrospectively nowcast the reporting proportion of Ebola cases (and hence the total number of cases) in Montserrado, and forecast the time to the next international case exportation from Montserrado assuming that the local epidemic would continue to grow exponentially at the nowcasted growth rate and the forward air-traffic would remain constant.

The Montserrado Ebola epidemic started in May 2014 32, 33 . By September 2014, two indigenous Ebola cases had been exported from Montserrado to other nations via commercial air travel: The first to Lagos, Nigeria, on 20 July 2014 34 ; and the second to Dallas, USA, on 19 September 2014 35 . By combining these global spread data with the World Health Organization patient database 33 on the weekly number of confirmed and probable cases in Montserrado and accounting for the effect of the opening of new Ebola treatment units in August 36,37 , we can use our framework to express the likelihood function in simple analytical form (see case study on the 2014 Liberian Ebola outbreak in Methods). We estimate that the reporting proportion (and hence the total number of cases) would have been statistically identifiable starting from 6 July 2014 onwards. We estimate that by 6 July 2014, the confirmed and probable cases only accounted for 18% (95% credible interval 7−33%) of all Ebola cases in Montserrado. The opening of new treatment units during August increased the reporting proportion to 30% (15−48%) by 17 August 2014, which is congruent with an independent estimate 38 based on capture-recapture sampling of raw patient records over a similar time horizon (34%; 95% confidence interval 26−50%). Retrospective real-time forecasts of the time to next exportation are consistent with the observed exportation times (namely 20 July and 19 September) except on 21−28 July during which the next exportation occurred later than predicted. The prediction errors on 21−28 July could be attributed to travel restrictions started in August 8 , the effect of which could not be included in the forecasts until they have actually occurred during August. If travel restrictions could have been foreseen on 21−28 July and incorporated into the forecasts (as a counterfactual scenario for illustration), the observed case exportation on 19 September 2014 would be consistent with the forecast range (Fig. 4b) . These conclusions are robust against temporal variations in epidemic growth rate (see case study on the 2014 Liberian Ebola outbreak in Methods).

In summary, our framework for characterizing the dynamics underlying global spread of epidemics comprises five approximations: (i) a closed-form pdf for EAT for any two directly connected populations (Eq. 1); (ii) adjustment for hub-effects; (iii) adjustment for continuous seeding; (iv) path reduction; and (v) path superposition. Approximation (i) is the indispensable centerpiece of our framework, whereas the necessity of approximations (ii)-(v) would depend on the specific application. Hubeffect adjustment is necessary when estimating the times of case exportation for populations that are directly connected to multiple populations and have relatively high outbound mobility rates. Continuous-seeding adjustment is necessary when estimating the times of case exportation for all populations except the epidemic origin (for which seeding is assumed to occur only at time 0). Path reduction and superposition are developed for simplifying computation as well as generating insights regarding global spread dynamics. In terms of computation, path reduction is required for populations that are three or more degrees of separation from the epidemic origin in a given acyclic path (which typically account for <20% of all populations), whereas path superposition is used for all populations (however, path superposition may not be necessary for populations that are directly connected to the epidemic origin with high mobility rates because the indirect paths have only minor impact on their EATs; see Supplementary Fig. 6 ). In terms of insights, the accuracy of path reduction implies that the epidemic arrival process from the epidemic origin to any given population along any given acyclic path ψ can be accurately approximated as an NPP with intensity function α ψ exp λ ψ t À Á , whereas the accuracy of path superposition implies that the dependence of multiple paths connecting a given population to the epidemic origin is relatively weak for the purpose of estimating EAT.

Although approximation (ii)-(v) are all necessary for estimating EAT in the WAN (Figs. 1-3) , they are not needed in our case studies on inference of transmission parameters: In the 2009 pandemic influenza A/H1N1 case study, we follow the inference formulation in Balcan et al. 26 which included only populations that are directly connected to Mexico City in the WAN. In the 2014 Ebola case study, the inference formulation tracks the timing of only two case exportations without the need to stratify them by outbound populations (see Case study on the 2014 Liberian Ebola outbreak in Methods).

Our study has several limitations. First, we did not consider age structure because the OAG air-traffic data do not have age information. If data are available for stratifying mobility rates and incidence by age, our framework should remain valid if the mobility rates w ij are calculated as the cross-product of agespecific mobility rates and age distribution of the disease. Second, we have assumed that each imported case spawns an exponentially growing infection tree with probability 1, whereas if we account for stochasticity in transmission dynamics, each imported case will fail to spawn an exponentially growing infection tree with probability p = 1/R 0 16 . Because this effect is similar to that of border control 19 , we conjecture that our framework can be extended to account for such stochasticity in transmission dynamics by discounting w ij with 1 − p. Third, we present our framework in the context of the classic SIR model. Nonetheless, our results can be generalized to all SE m I n R models 39 (see Generalizing to SE m I n R models in Methods). Fourth, we have not accounted for seasonality effects which may be strong and geographically heterogeneous for diseases such as seasonal influenza 40, 41 . Although the epidemic dynamics will certainly be less analytically tractable in the presence of seasonality (e.g., the pdf of the EAT can no longer be well approximated by simple closedform expressions as Eq. 1), we conjecture that the new analytics introduced here, namely adjustments for the hub-effect and continuous seeding as well as path reduction and superposition, will be useful for building a more general framework for global spread of epidemics. Finally, in our case studies, we have implicitly assumed that surveillance data were available in near realtime for nowcasting and forecasting, whereas in reality the availability of reliable data would likely incur longer lead times, and hence the timeliness of situational awareness implied here should be interpreted within such context.

In summary, we have developed a novel framework that can accurately characterize how global spread of epidemics depends on the infectious disease epidemiology and network properties of the WAN. Together with state-of-the-art global epidemic simulators such as GLEAM, our framework advances the frontiers of the next-generation informatics for pandemic preparedness and responses. is the basic reproductive number and T g,j is the mean generation time in population j. Let β j = R 0,j /T g,j be the disease transmission rate in population j, and w jk be the mobility rate from population j to k. The stochastic metapopulation model with G populations is specified by the following equations where Δt is a very small time interval: 

where X jk (t),Y jk (t), and Z jk (t) are the kth component of X j (t), Y j (t), and Z j (t), respectively. Multinomial(n, p 1 ,...,p G ) denotes a multinomial random variable with n trials and probabilities p 1 ,...,p G . We use Δt = 0.05 days in all of our simulations.

The global epidemic simulator. We build the global simulator using 2015 worldwide flight booking data from the Official Airline Guide (OAG, https://www. oag.com/) and the Gridded Population of the World Version 4 (GPWv4, http:// sedac.ciesin.columbia.edu/data/collection/gpw-v4/) data set from the NASA Socioeconomic Data and Applications Center (SEDAC) at Columbia University. Worldwide air-transportation data: Our OAG worldwide flight booking data set contains all air bookings that have taken place in all commercial airports worldwide during 2015. Each data record contains the following information for a flight route: (i) origin airport, (ii) destination airport, (iii) connecting airports (if any), and (iv) passenger bookings for each month. The city and country served by each airport and the coordinates of each airport are known. The raw data comprises 0.947 million records. Parameterizing the WAN using these raw data would therefore generate 0.947 million connections in the network, which is beyond our computational capacity and unnecessary for an accurate description of global spread (because the WAN is densely connected) 3-6,12,42-44 . As such, we perform the following steps to exclude flight routes with weak traffic from the WAN without compromising the realism of the global epidemic simulator:

1. We exclude all routes with no bookings for one or more months during 2015. 2. We exclude all routes in which the origin or destination is a remote area with very small population size (e.g., hamlets, settlements, or communities in Alaska and Northern Canada). 3. We exclude all routes with strong seasonality as measured by normalized information entropy 11 : H ij ¼ À 1 logð12Þ P 12 m¼1 ρ ijm log ρ ijm , where ρ ijm ¼ F ijm = P 12 m¼1 F ijm and F ijm denotes the number of air bookings from origin airport i to destination airport j in month m. The measure H ij ranges between 0 and 1, and decreases as temporal variation in air-traffic increases (e.g., if air-traffic is the same across all months, then H ij = 1). On the basis of the distribution of H ij in our OAG raw data, we exclude all routes with H ij < 0.8 ( Supplementary Fig. 8a) . Global population data: The GPWv4 data set integrates the highest resolution census data from the 2010 round of Population and Housing Censuses collected from hundreds of national statistics departments and organizations 45, 46 . GPWv4 provides eight different data sets, most of which are specialized geospatial metadata that partition the global population into a grid of cells with resolution of 30 arcsecond (~1 km at the equator). We use the vector data set "Administrative Unit Center Points with Population Estimates, v4 (2000, 2005, 2010, 2015, 2020)" 47 , because it provides all the information that we need to build the global epidemic model, e.g., the coordinates of centroid are available for each of the~12.5 million administrative census units (ACUs).

The WAN model: We combine our OAG data with the GPWv4 data to calculate the population size of the catchment area of each airport as follows:

1. We use the coordinates of the centroids of all ACUs and airports to calculate the great circle distance for all possible combinations of ACUs and airports within the same country. We use a Voronoi-like tessellation algorithm proposed by Balcan et al. 23 to link each ACU to its serving airport (i.e., the closest airport in its country). In this algorithm, we impose the constraint that the great circle distance between any pair of ACU and airport cannot exceed 200 km, according to the distribution of great circle distance for all combinations of ACUs and airports ( Supplementary Fig. 8b ). This reflects a reasonable upper bound on the distance of land transportation for reaching an airport 23 . Without this constraint, the algorithm may generate unreasonably large catchment areas for airports located in sparsely populated regions. Among the 7,995,985 ACUs with human habitats, only 45,692 are excluded from our model because of this constraint. The total population size served by an airport is the sum of populations for all ACUs assigned to that airport.

and realism of our global epidemic simulator, we exclude all routes having less than 3000 air passengers throughout the year (Fig. S8c-f ). This simplification is in line with the passenger threshold reported by Khan et al. 48 and hence has little impact on the accuracy of global spread dynamics. 3. In our OAG data set, some metropolitans (e.g., London, New York City, and Shanghai) and tourist locations (e.g., Hawaii and Canary Islands) have multiple airports. We model each of these locations as a single population by merging its serving airports and the corresponding catchment areas. 4. The daily air-traffic of each connection F ij is the average number of air passengers per day for that connection during the year of 2015. The ensemble of all connections shows a high degree of statistical symmetry, F ij ≈ F ji (R 2 = 0.9981), as in refs. [3] [4] [5] [6] [7] 11, 18, 20, [42] [43] [44] . As such, we symmetrize the air-traffic between each pair of populations by setting F ij = F ji = (F ij + F ji )/2. In summary, the WAN in our global metapopulation epidemic model comprises 54,106 connections and 2309 populations and preserves more than 92% of the global air bookings.

Details on assumption 1. Assumption 1 is stated as follows: suppose populations j and k are directly connected in the WAN and only population j is infected. Exportation of infections from population j to k is an NPP 22 with intensity function w jk I j (t) where I j (t) is the disease prevalence in population j at time t.

Previous studies 19,21 on global spread have made similar assumptions. For populations j and k mentioned above, the exportation process of infections from population j to population k clearly satisfies conditions 1 with intensity function w jk I j (t). If the mobility rate w jk is sufficiently small, the number of exportations is only a very small proportion of the disease prevalence in population j, and hence conditions 2 and 3 are also satisfied.

The two-population model analysis. Population i is the epidemic origin and only connected to population j. Let s i and λ i be the seed size and the initial epidemic growth rate. Let X ij be the total number of infections imported by population j over 1. X ij is Poisson distributed with mean A i T g F ij , where T g is the mean generation time, A i is the final attack rate in population i and F ij is the daily average number of passengers traveling from population i to j. That is, P(X ij = n) = f Poisson (n,A i T g F ij ). 2. Applying the framework of NPP 22 , we express the pdf of T n ij conditional on X ij ≥ n as

Supplementary Figure 2 shows that the pdf in Eq. S1 is very accurate for all realistic epidemic scenarios. If assumption 2 is also valid, i.e., I i ðtÞ ¼ s i exp λ i t ð Þ, then P (X ij ≥ n) = 1 and Eq. S1 can be simplified to

which is Eq. 1 in the main text with α ij = s i w ij . The corresponding cumulative distribution function (cdf) is given by

where Γ is the lower incomplete gamma function. The expected EAT is given by

u m du is the exponential integral. If α ij ≪ λ i and γ denotes the Euler constant, we obtain the following approximation

which is congruent with the EAT statistic in Gautreau et al. 18 for estimating the order of epidemic arrival across different populations. The expected time of the nth exportation is given by

For any positive integers m and n such that m < n, the pdf of T n ij À T m ij conditional on T m ij is simply f nÀm tjλ i ; α ij exp λ i T m ij which corresponds to the time of the (n − m)th exportation for an epidemic with seed size s i exp λ i T m ij . Using this relation recursively, we deduce that the joint pdf of T 1

which is the basis that supports our likelihood-based inference framework. By the same token,

The WAN-SPT analysis. Hub-effect: Suppose the epidemic origin (population i) is directly connected to one or more populations, one of which is population j (as illustrated in Fig. 2a) . In the deterministic version of our metapopulation epidemic model (see WAN metapopulation epidemic model in Methods), the disease prevalence in population i during the exponential growth phase is well approximated by the differential equation

where the actual growth rate of the disease prevalence in population i is λ i À P k w ik . This differential equation leads us to make the following conjecture: In our original stochastic model, in which the epidemic arrival process for population j is essentially an NPP with intensity function being the second term of the above equation (i.e., w ij I i ), we can estimate the EAT for population j using the results from the two-population model (The two-population model analysis in Methods) in which population i is exporting cases to population j at mobility rate w ij (viewed as a stochastic process) and the disease prevalence in population i is growing exponentially at rate λ ij = λ i − ∑ k ≠ j w ik (viewed as a deterministic process). The hub-adjusted growth rate λ ij can be interpreted as the rate at which disease prevalence in population i is growing exponentially before population j imports its first case from population i. Note that the hub-adjusted rate λ ij = λ i − ∑ k ≠ j w ik is not the same as the actual growth rate, namely λ i À P k w ik . To see this, consider the two-population model in which population i is only connected to population j. In this case, the EAT distribution is given by Eq. 1, which requires λ ij to be the hubadjusted rate λ i − ∑ k ≠ j w ik = λ i but not the actual growth rate λ i À P k w ik ¼ λ i À w ij . Continuous seeding: Consider the path connecting the epidemic origin to population k via population j, i.e., ψ : i → j → k. Let λ ij and λ jk be the hub-adjusted growth rate in populations i and j for this path. Under assumption 2, the prevalence in population j at time t that are spawned by the mth infection imported

where I{ ⋅ } is the indicator function. Therefore, the total prevalence in population j at time t is

The NPP intensity function for the exportation of infections from population j to population k is w jk I j (t).

Conditional on I j and hence T 1 ij , T 2 ij ,…, the pdf of T n ik is g n tjw jk I j

where the joint pdf of T 1 ij ¼ t 1 ,T 2 ij ¼ t 2 ,…, is simply the product of f 1 t m jλ ij ; w ij s i exp λ ij t mÀ1 À Á À Á for m = 1, 2, … (see Eq. S2). As described in the main text, we make the certainty equivalent assumption (CEA) that conditional

As such, conditional on T 1 ij , we approximate I j with

Eq. 2 and the previous section).

The resulting unconditional pdf of T n ik is E T 1 ij g n tjw jk I CEA

Path reduction: Consider the path ψ: i → j → k in the previous section. We can approximate the pdf E T 1 ij g n tjw jk I CEA j h i for T n ik with f n (t|λ ψ ,α ψ ), where λ ψ and α ψ are obtained by minimizing the relative entropy 25 for n = 1 (the first exportation)

This is a simple two-dimensional optimization problem. The accuracy of such path reduction ( Fig. 2f and Supplementary Fig. 5 ) implies that the spread of epidemics from the origin to any population k ∈ D i,2 can be regarded as a two-population model, in which (i) the adjusted mobility rate is α ψ and (ii) the epidemic in the origin grows exponentially at rate λ ψ . Next, consider the path ϕ: i → j → k → m, i.e. m ∈ D i,3 . Using path reduction, we can approximate ϕ with ϕ': i → k → m where the adjusted mobility rate and epidemic growth rate in the origin for the i → k leg are α ψ and λ ψ , respectively. The arrival times of imported cases in population m ∈ D i,3 (i.e., T n im , n = 1, 2, …) can then be estimated using the tools (i.e., adjustments for hub-effect and continuous seeding) that we have developed for D i,2 populations. The arrival times of imported cases for population D i,c , c = 4, 5, …, can be estimated analogously.

The WAN analysis. Superposition of paths: Let population i be the epidemic origin and consider population k ∈ D i,c , i.e., population k is c degrees of separation from the epidemic origin 24 . Superposition of NPPs for paths connecting population i to k is implemented as follows. As in the main text, let Ψ ik be the set of all acyclic paths connecting the epidemic origin to population k. Enumeration of all paths in Ψ ik for every population in the WAN is computationally prohibitive 49 (and unnecessary). Instead, we approximate Ψ ik with the 25 "fastest" paths from population i to k that are identified using the following algorithm:

1. Use the depth-first search algorithm 49 to identify the set of acyclic paths from the epidemic origin to population k that have at most c + 2 connections. We denote this set by Ω ik and assume that all the paths not in Ω ik have negligible contribution to the EAT for population k. 2. Define the distance between any two directly connected populations a and b as À ln w ab ð Þ, which is analogous to the distance metric in Brockmann and Helbing 20 , namely 1 À ln w ab = P b w ab À Á . We choose to use this distance metric because (as described in The two-population model analysis in Methods) if population j is directly connected to population i, then E T 1

given α ij ≪ λ i , where γ denotes the Euler constant and α ij = s i w ij . This indicates that the expected EAT is proportional to À ln w ij À Á . 3. Based on our distance metric in step 2, identify the 100 shortest paths in Ω ik by sorting in an ascending order. Denote the resulting set by Ω S ik . 4. For each path ψ 2 Ω S ik , use hub-effect adjustment, continuous-seeding adjustment and path reduction developed in the WAN-SPT analysis to calculate λ ψ and α ψ and the corresponding expected EAT, namely 1 λψ exp αψ λψ E 1 αψ λψ .

5. Approximate Ψ ik with the 25 paths in Ω S ik that have the smallest expected EATs computed in step 4 (i.e. the 25 "fastest" paths). We choose to use the 25 fastest paths in Ω S ik to approximate Ψ ik because Supplementary Fig. 9 shows that the accuracy of EAT estimates would slightly worsen if we use only the 10 fastest paths in Ω S ik while there is little improvement in performance if we use the 50 fastest or all paths in Ω S ik .

Generalizing to SE m I n R models. In the main text, our framework is built using the SIR model within each population. In this section, we describe how to generalize our framework to SE m I n R models 39 in which:

1. The duration of latency is gamma distributed with mean D E and m subclasses (i.e., with shape m and rate b E = m/D E ); 2. The duration of infectiousness is gamma distributed with mean D I and n subclasses (i.e., with shape n and rate b I = n/D I ). For any given population, let S(t) be the number of susceptible individuals, E i (t) the number of individuals in the ith latent subclass, and I j (t) the number of individuals in the jth infectious subclass. The SE m I n R system is described by the following differential equations: During the early stage of the epidemic (such that S(t) ≈ N), the prevalence of latent and infectious individuals both grows exponentially at rate λ which is the solution to the following equation 39 :

That is, the prevalence of latent and infectious individuals are well approximated by E exp λt ð Þ and I exp λt ð Þ, respectively, where E and I depend on the initial conditions and parameters of the differential equation systems (the analytical expressions of E and I can be obtained by solving the linearized system with S(t) = N). As such, if a proportion 1 − p E and 1 − p l of the latent and infectious individuals refrain from air travel because of their infections, then the seed size s 0 in the main text is simply p E E þ p I I.

Case study on the 2009 influenza A/H1N1 pandemic. As described in the main text, by integrating our framework into the inference formulation in Balcan et al. 26 , we express the likelihood function for the EATs for the 12 countries seeded by Mexico (see Supplementary Table 1) as

in which population i (the epidemic origin) is Greater Mexico City 50 where the epidemic began in mid-February to early-March 2009 27, 30 , t j is the observed EAT for population j which can be exact (A) or left-censored (B), λ ij = λ i − ∑ k ≠ j w ik is the hub-adjusted growth rate, α ij is the adjusted mobility rate. Because the air travel data were not reported in Balcan et al. 26 , we use the air travel data published in Fraser et al. 27 in which the basic reproductive number R 0 was estimated from the number of confirmed cases in different countries seeded by Mexico during March-April 2009. Supplementary Table 1 shows the EAT data from Balcan et al. 26 and the air-passenger data from Fraser et al. 27 The population size of Greater Mexico City in 2009 was 17.6 million 27 . We assume that the epidemic started with a single infected individual (i.e., s i = 1) in Greater Mexico City between 18 February and 14 March 2009 based on the documentation in surveillance reports 29 and other studies 27, 28, 30, 31 (Fig. 4a) . We adopt the natural history model described in Balcan et al. 26 : (i) the mean generation time is T g = 3.6 days with mean latent duration of 1.1 days; (ii) the latent and infectious duration are exponentially distributed (regardless of symptoms). Under these assumptions, the basic reproductive number is R 0 = (1 + λ i × mean latent period)(1 + λ i × mean infectious period) 51 . Balcan et al.

assumed that 67% of infections are symptomatic and 50% of symptomatic infections refrained from traveling by air. As such, we discount the mobility rates by multiplying w ij with 0.5 × 0.67 = 0.335.

In this case study, R 0 is the only parameter subject to inference. We assume non-informative flat prior and use the Metropolis-Hasting algorithm 52, 53 to estimate the posterior distributions of R 0 . We use five MCMC chains and initialize each chain with an R 0 value randomly chosen between 1 and 10. The trace plot and Geweke diagnostic indicate that each MCMC chain converges within 5000 iterations and the autocorrelation of the samples in the MCMC chain is essentially 0 when the lag is larger than 10 steps. As such, we estimate the posterior distribution of R 0 by running the Metropolis-Hasting algorithm for 110,000 iterations with a burn-in of 10,000 iterations and a thinning interval of 10. The Gelman-Rubin diagnostic indicates that all five chains converge to the same posterior distribution.

Case study on the 2014 Liberian Ebola outbreak. In 2014, the first laboratory confirmed Ebola case in Montserrado, Liberia, developed symptoms during the week of 5 May 2014 32, 33 . During this Ebola epidemic, two Ebola cases were exported from Montserrado to the following populations via international commercial air travel 34,35 :

1. Lagos, Nigeria on 20 July 2014 (t 1 ); 2. Dallas, USA on 19 September 2014 (t 2 ).

Montserrado and Margibi were the major epicenter in Liberia during the 2014 West African Ebola epidemic 32, 33, 54 , and they are served by the two contiguous Liberian commercial airports that have international flights (i.e., Roberts International Airport and Spriggs Payne Airport). In this case study, we apply our framework to estimate the reporting proportion and the total number of Ebola cases in Montserrado and Margibi (Montserrado hererafter for brevity) between 5 May 2014 (the approximate start time of this epidemic) and 21 September 2014 (the last day of the week during which the last exportation occurred). Based on ref. 54 , we assume that the latent period and the incubation period were the same. We assume that infectious cases did not travel by air (due to their symptoms), and exportations comprised only air travel of latent individuals (who had not yet developed symptoms). We note that there was some evidence 55 that the case exported to Lagos had already developed symptoms when he boarded the flight. Therefore, we include this case in our main analysis but exclude him in the sensitivity analysis. Results from both analyses are essentially the same (Supplementary Fig. 10) .

Let time 0 be 5 May 2014 and T be 21 September 2014. We denote May, June, July, August and September 2014 by months 1 to 5, respectively. Denote the last day of month k since time 0 by τ k , and the two observed times of case exportations since time 0 by t 1 and t 2 , respectively. We assume that the incidence rate was (i) 0 before 5 May 2014, (ii) i 0 on 5 May 2014, and (iii) i 0 exp λt ð Þ thereafter, i.e., this epidemic grew exponentially at rate λ between 5 May and 21 September 2014. The incubation period has been estimated to be gamma distributed with shape m = 1.41 and rate b E = 0.154 (which correspond to mean 9.2 days and standard deviation 7.7 days) 54 . Hence, symptomatic cases occurred at rate inc sym ðtÞ ¼

where g is the pdf of the incubation period, and Γ is the lower incomplete gamma function. Accordingly, the number of new symptomatic Ebola cases in the kth week since time 0 was

Let θ k be the probability that a true case with onset in week k was reported as confirmed or probable cases. New Ebola treatment units were established in Montserrado in early August 2014 36, 37 . As such, we assume that θ k = θ before if week k ended before 4 August 2014, and θ k = θ after otherwise. The likelihood for the observed number of confirmed and probable cases is

where y k is the observed number of confirmed and probable cases with onset in week k and f binomial is the binomial pdf. The observed weekly number of confirmed and probable Ebola cases in Montserrado is obtained from the World Health Organization (WHO) patient database 33 .

Our OAG data set also contains the monthly number of flight bookings in 2014. Supplementary Table 2 shows the monthly outbound mobility rates from Montserrado during May-September 2014 in this OAG database. We denote the outbound mobility rate from population i during month k by W k (i.e., W k = ∑ j w ijk , where w ijk is the daily mobility rate from population i to j during month k). Air travel restrictions were implemented starting in August 8 , which presumably resulted in a substantial proportion of canceled flight bookings (in particular for August 2014, see Supplementary Table 2 ). These abnormal cancellations were not registered in the OAG database. Therefore, as an approximation, we assume that the actual mobility rate in August 2014 was the same as that in September 2014.

According to our framework, if population i has seed size s, epidemic growth rate λ, and outbound mobility rate w, the probability that population i has no exportation up to time t is 1 À F 1 tjλ; sw ð Þ¼exp À sw λ exp λt ð Þ À 1 ð Þ and the probability density that population i has its first exportation at time t is

Given the incidence rate i 0 exp λ i t ð Þ, the seed size of latent infections was effectively

Þ . To see this, consider an SE m I n R system during the exponential growth phase with incidence rate i 0 exp λ i t ð Þ. The prevalence of latent individuals during this phase is well approximated by the following system:

; m :

Solving these differential equations gives In summary, we infer (λ,i 0 ,θ before ,θ after ) using the likelihood L λ; i 0 ; θ before ; θ after ð Þ ¼ L inc λ; i 0 ; θ before ; θ after ð Þ L export λ; i 0 ð Þ:

Note that θ after is defined only after 3 August 2014 and hence not inferred until then. We assume non-informative flat priors for all parameters and use Gibbs sampling 52 to estimate the posterior distributions of (λ, i 0 , θ before , θ after ). We use five MCMC chains and initialize each chain with a starting point randomly generated from the following ranges: lnð2Þ=λ (i.e., the doubling time) between 1 and 100 days, i 0 between 1 and 100, θ before between 0 and 1, and θ after between 0 and 1. The trace plot and Geweke diagnostic indicate that each MCMC chain converges within 100,000 iterations and the autocorrelation of the samples in the MCMC chain drops below 0.05 when the lag is larger than 2000 steps. As such, we estimate the posterior distribution of (λ, i 0 , θ before , θ after ) by running the Gibbs sampling for 5.5 million iterations with a burn-in of 0.5 million iterations and a thinning interval of 5000. The Gelman-Rubin diagnostic indicates that all five chains converge to the same posterior distribution. Given an estimate of (λ, i 0 , θ before , θ after ), the cumulative number of symptomatic Ebola cases up to time t was:

and the reporting proportion up to the end of week K was X K k¼1 y k =Cð7KÞ:

The nowcasted posterior estimates of the epidemic doubling time (which is simply ln(2)/λ) and initial incidence rate (i 0 ) are temporally consistent until mid-August after which both began to increase significantly. This suggests that the epidemic growth rate might have dropped since mid-August, which is plausible in view of substantial increase in mitigation efforts and resources starting in early August 36, 37 . As such, we perform a sensitivity analysis by assuming that the epidemic doubling time changed from D 1 to D 2 starting on 4 August 2014. Supplementary Figure 10 shows that our main result, namely the estimates of reporting proportion, remain essentially the same. For the scenario unadjusted for travel restrictions (Fig. 4b, bottom panel) , the retrospective real-time forecasts of the time to next international exportation are obtained by (i) assuming that mobility rates during the forecasted time period were the same as the most current mobility rates and (ii) sampling (λ, i 0 , θ before , θ after ) from their posterior distributions.

Code availability. Code is available on request from the authors.

Data availability. Global population data (raw data) that support the findings of this study are available from the Gridded Population of the World Version 4 (GPWv4) database at http://sedac.ciesin.columbia.edu/data/collection/gpw-v4. Restrictions apply to the availability of the worldwide air-traffic data set from the Official Airline Guide (https://www.oag.com/), which were used under license for the current study. Source data for case studies ( Fig. 4; Supplementary Fig. 10 

A mathematical model for the global spread of influenza

Strategies for mitigating an influenza pandemic

Delaying the international spread of pandemic influenza

Unifying viral genetics and human transportation data to predict the global transmission dynamics of human influenza H3N2

Forecast and control of epidemics in a globalized world

Will travel restrictions control the international spread of pandemic influenza?

Assessment of the Middle East respiratory syndrome coronavirus (MERS-CoV) epidemic in the Middle East and risk of international spread using a novel maximum likelihood analysis approach

Assessing the impact of travel restrictions on international spread of the 2014 West African Ebola epidemic

Potential for Zika virus introduction and transmission in resource-limited countries in Africa and the Asia-Pacific region: a modelling study

Spread of Zika virus in the Americas

The role of the airline transportation network in the prediction and predictability of global epidemics

Real-time numerical forecast of global epidemic spreading: case study of 2009 A/H1N1pdm

Meta)population dynamics of infectious diseases

Metapopulation dynamics

Large-scale spatial-transmission models of infectious disease

Modeling Infectious Diseases in Humans and Animals

Modelling dynamical processes in complex socio-technical systems

Global disease spread: Statistics and estimation of arrival times

A simple explanation for the low impact of border control as a countermeasure to the spread of an infectious disease

The hidden geometry of complex, networkdriven contagion phenomena

Fluctuation effects in metapopulation models Percolation and pandemic threshold

Stochastic Processes 2nd edn

Multiscale mobility networks and the spatial spreading of infectious diseases

Collective dynamics of 'small-world' networks

Elements of Information Theory 2nd edn

Seasonal transmission potential and activity peaks of the new influenza A(H1N1): a Monte Carlo likelihood analysis based on human mobility

Pandemic potential of a strain of Influenza A (H1N1): early findings

Use of cumulative incidence of novel influenza A/H1N1 in foreign travelers to estimate lower bounds on cumulative incidence in Mexico

Outbreak of swine-origin influenza A (H1N1) virus infection-Mexico

The transmissibility and control of pandemic influenza A (H1N1) virus

Initial human transmission dynamics of the pandemic (H1N1) 2009 virus in North America

Evolution and spread of Ebola virus in Liberia

Ebola data and statistics

Transmission dynamics and control of Ebola virus disease outbreak in Nigeria

Ebola virus disease cluster in the United States

Strategies for containing Ebola in West Africa

Impact of interventions and the incidence of Ebola virus disease in Liberia-implications for future epidemics

Use of capture-recapture to estimate underreporting of Ebola virus disease

Appropriate models for the management of infectious diseases

Environmental predictors of seasonal influenza epidemics across temperate and tropical climates

Inference of seasonal and pandemic influenza transmission dynamics

Sampling for global epidemic models and the topology of an international airport network

The cost of simplifying air travel when modeling disease spread

Hedging against antiviral resistance during the next influenza pandemic using small stockpiles of an alternative chemotherapy

Documentation for the Gridded Population of the World, Version

Taking advantage of the improved availability of Census Data: A first look at the Gridded Population of the World, Version 4

Administrative Unit Center Points with Population Estimates (NASA Socioeconomic Data and Applications Center (SEDAC)

Spread of a novel influenza A (H1N1) virus via global airline transportation

An Introduction

Characterizing the epidemiology of the 2009 Influenza A/H1N1 Pandemic in Mexico

How generation intervals shape the relationship between growth rates and reproductive numbers

Markov Chain Monte Carlo in Practice

Fractional dosing of yellow fever vaccine to extend supply: a modelling study

WHO Ebola Response Team. Ebola virus disease in West Africa-The first 9 months of the epidemic and forward projections

Ebola virus disease outbreak-Nigeria

We thank M. Lipsitch, B. J. Cowling, J. M. Read, K. Leung, H. Choi and Y. Zhang for helpful discussions. We thank C. K. Lam for assistance in data processing and technical support. We thank the Official Airline Guide, Center for International Earth Science Information Network at Columbia University, and World Health Organization (WHO) for providing their databases. This research was conducted in part using the research computing facilities and advisory services offered by 

Supplementary Information accompanies this paper at https://doi.org/10.1038/s41467-017-02344-z.Competing interests: The authors declare no competing financial interests.Reprints and permission information is available online at http://npg.nature.com/ reprintsandpermissions/ Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/.