Country distancing increase reveals the effectiveness of travel restrictions in stopping COVID-19 transmission

Despite a number of successful approaches in predicting the spatiotemporal patterns of the novel coronavirus (COVID-19) pandemic and quantifying the effectiveness of non-pharmaceutical interventions starting from data about the initial outbreak location, we lack an intrinsic understanding as outbreak locations shift and evolve. Here, we fill this gap by developing a country distance approach to capture the pandemic’s propagation backbone tree from a complex airline network with multiple and evolving outbreak locations. We apply this approach, which is analogous to the effective resistance in series and parallel circuits, to examine countries’ closeness regarding disease spreading and evaluate the effectiveness of travel restrictions on delaying infections. In particular, we find that 63.2% of travel restrictions implemented as of 1 June 2020 are ineffective. The remaining percentage postponed the disease arrival time by 18.56 days per geographical area and resulted in a total reduction of 13,186,045 infected cases. Our approach enables us to design optimized and coordinated travel restrictions to extend the delay in arrival time and further reduce more infected cases while preserving air travel. Network approaches are key to understand epidemic spreading, inherently driven by human mobility patterns and constrained by transport systems. In this work, the authors develop a country distance framework to capture the spread of COVID-19 on top of the airline network, analyzing the effectiveness of mobility restrictions in the presence of multiple outbreaks and suggesting strategies for optimized coordinated travel restrictions.

G iven 41,570,883 confirmed cases of COVID-19 and 1,134,940 deaths worldwide as of 23 October 2020 1-3 , the need for world to deploy non-pharmaceutical interventions prior to comprehensive consideration is urgent [4][5][6] . Today's high population density and the high volume, speed, and nonlocality of human mobility provide perfect conditions for an epidemic to spread [7][8][9][10] and simultaneously raise the challenges related to the development of non-pharmaceutical intervention strategies on the timescale that modern diseases spread [11][12][13] . Although the practice of quarantine and social distancing protocols can drastically reduce the virus propagation locally 6,[14][15][16] , the global COVDI-19 pandemic patterns are shaped by the global mobility network (GMN), which determines when and where the disease arrives globally 17,18 . Consequently, the straightforward way to diminish the international importation of COVID-19 involves imposing radical travel restrictions (i.e., entry bans, global travel bans, and lockdowns) [19][20][21] , which reduce the entry of airline passengers into a country. According to available data as of 1 June 2020, 187 geographical areas imposed the entry bans, 111 geographical regions imposed global travel bans, and 38 geographical areas imposed the full lockdowns to prevent their citizens and tourists from traveling overseas 5,22 . However, researchers demonstrated that these travel restrictions were only effective at the beginning of the outbreak 23 . When such interventions fail to control the initial outbreak, they instead disrupt the healthcare aid and support, business, and cause extensive and profound social and economic damage 24,25 . Therefore, it is crucial to assess and impose effective travel restrictions to avoid unnecessary and costly responses to COVID-19 from the side of governments 26 .
Measuring the effectiveness of travel restrictions often relies on the specific epidemic models 4,[27][28][29][30] , which require accurate estimation of the disease's epidemiological parameters, such as the basic reproductive number (R 0 ). However, the parameter estimations are often unreliable due to the daily changes in underreported cases and various errors arising from the lack of diagnosis tests 11,31,32 . Furthermore, these models are difficult to calibrate due to incomplete information (i.e., partial network topology 33 or unknown dynamics 19,34 ). Overall, it is unclear how much detail is required to achieve a certain level of predictive accuracy. Human mobility plays a significant role in understanding hidden spatiotemporal spreading patterns 35 and enables us to predict the arrival time 17,36,37 and number of infected cases 38 . In particular, effective distance, a method introduced by Brockmann and Helbing 17 , measures the mobility from the initial outbreak location (OL) to the target geographical area by discarding other redundant connections. This method allows to predict the arrival time and epidemic wavefront without knowing the epidemiological parameters and has already been demonstrated useful in the 2009 H1N1 pandemic and the 2003 SARS epidemic. Additionally, the initial OL's mobility outflow is a vital predictor for the log-transformed cumulative infections in the destination locations 38,39 , and was validated by the Wuhan's outflow to each prefecture in mainland China in the early stage of the COVID-19 pandemic.
Here, by simulating the spread of disease with the metapopulation susceptible-infected-removed (SIR) model, we show that the effective distance method is successful in predicting the evolution of the pandemic when Mainland China is the only OL. As shown in Fig. 1a, b, the effective distance accurately estimates countries' arrival times with an R-squared value R 2 = 0.87 and log-transformed cumulative infected cases with R 2 = 0.88. However, when multiple OLs appeared, effective distance weakly correlated with disease dynamics with R 2 ≤ 0.4 (see in Fig. 1c, d). We argue that the observed underperformance is due to the effective distance's underlying assumption of a single OL, which neglects other OLs' potential in exporting viruses to other areas through air travel. The presence of varying OLs is essentially related to the variability in countries' health responses to COVID-19. In fact, the time-varying nature of OLs poses challenges in capturing the ongoing pandemic's spatiotemporal pattern and evaluating the effects of travel restrictions. This feature prompts new mathematical tools to advance the understanding of global disease dynamics under unprecedented non-pharmaceutical interventions. To address this gap, we propose a country distancing method that captures global disease dynamics when OLs evolve. Our approach takes inspiration from the connection laws in electric circuits with resistances in series and parallel configurations. We apply this approach to quantify the effectiveness and efficiency of existing travel restrictions and design optimized and coordinated travel restrictions to maximize their impact in slowing infections.

Results
Mathematical framework of country distancing. We tested the international spread of the disease on the GMN (see Supplementary Table 1). The GMN, which is provided by the Official Aviation Guide, is presented by a matrix of airline passenger influx F among 250 geographical areas, representing countries, dependent territories, and special areas of geographical interest. Here, F mn (F mn ∈ F) expresses the airline passenger influx from area n to the area m. For a single OL, the diseases may propagate to a geographical area through distinct paths, the shortest of which predicts the infection to the destination 17 . For area n and its connected area m, their effective distance is d mjn ¼ 1 À log P mn , indicating that a larger fraction of air travel P mn (P mn ¼ F mn ∑ k F kn ) means a smaller distance, and vice versa. Then, for an arbitrary area m that can be reached by n through a path τ, the effective distance is the sum of effective lengths along the links of the shortest path, d mjn ¼ min τ ∑ ði;jÞ2τ d jji . We call this series connection law (see Fig. 2a) because it is analogous to the effective resistance in series circuits defined as R = ∑ i R i , where R i is the resistance that is connected along a chain.
The number of OLs may grow or shrink, which escalates or diminishes the importation risk for the other geographical areas and simultaneously increases/reduces geographical areas' distance to the risky sources, respectively. Using the meta-population SIR model, the disease evolution on this model becomes segmented according to the presence of OLs (see "Methods" and Supplementary Note 2). Motivated by the effective distance 17 and propagation processes simulated by the meta-population SIR model with multiple OLs, we derive a formula that examines areas' closeness to all existing OLs through the parallel connection law. For example, the disease propagates from two OLs, n and c, to the destination geographical area, m, with effective distances d m|n , and d m|c , respectively (see Fig. 2b). The overall likelihood of transmitting from both OLs satisfies e d mj n;c f g / where N I is the OL set and M is the total number of geographical areas (M = 250). Note that a large set of OLs may lead to small country distancing, and a larger portion of passenger influx leads Global disease dynamics cannot be predicted by effective distance when the outbreak locations (OLs) are evolving. a, b depict the global disease dynamics, i.e., arrival times T m and log-transformed infected cases I m (t) (e.g., t = 50) as functions of effective distance from the initial OL, respectively. However, when the single OL evolves to multiple OLs (e.g., five OLs), see c and d, the effective distance from the initial OL fails to correlate with arrival times T m and log-transformed infected cases I m (t), with an R-squared value (R 2 ) < 0.4. Circles represent geographical areas, and those belonging to the same continent are presented in the same color. Both arrival time and infected cases are obtained from the meta-population susceptible-infected-removed model, with epidemiological parameters provided by the literature 55,56 . Global disease dynamics using country distancing when the outbreak locations (OLs) are evolving. By deriving the series and parallel connection law (a, b) for global disease transmission, we derive the country distancing method to capture the global disease spread dynamics, i.e., the arrival times and new infected cases. R n and R c represent the resistances in the circuit. Nodes n, c, and m represent the geographical areas and the directed links between the nodes represent the traffic flow. The disease propagates from one OL (e.g., n) to destination m, on par with series connection law; whereas the disease propagates from multiple OLs (e.g., n and c) to destination m, on par with parallel connection law. c, d separately show the R-squared values (R 2 ) for evaluating global disease dynamics versus country distancing and effective distance. Country distancing correlates strongly with remaining arrival times T m (t) (c) and log-transformed new infected cases I m (t) − I m (t − 1) (d) when multiple OLs are present. The R 2 values for country distancing keep above than 0.7 although OLs are shifting, whereas the R 2 values for effective distance decline sharply given that multiple OLs emerge. e, f are the slopes of linear correlations for c and d, representing the speed of arrival times and speed of infection. The shaded areas are the 95% confidence intervals for the speeds. Outbreak locations are defined as the geographic areas whose active confirmed case greater than 0.01% population. On 1 June 2020, 146 geographical areas are OLs. Similar to Fig. 1, both arrival time and infected cases are obtained from the meta-population susceptible-infected-removed model. to a smaller effective distance and may further lead to smaller country distancing.
Global disease dynamics for the evolving outbreak locations. Two fundamental properties describe the main spread dynamics of the COVID-19 pandemic: the arrival time (T m ), i.e., the date of the first confirmed case, and the cumulative infected cases [I m (t)] in a geographical area m. As numerous undetected, missing, undiagnosed, or unreported COVID-19 cases result in biased arrival times and infected cases in the collected real-world dataset 31,32,40 , we simulate the spread of COVID-19 by adopting the meta-population SIR model 41 in the segmented time interval according to the presence of different OL sets. The model is validated by (1) defining the OL N I (t) as the geographic areas having greater than 0.01% of the infected population at time t; (2) ratifying the strengths for entry bans, global travel ban, and lockdown in 90% effective in limiting travels according to the International Civil Aviation Organization 42 ; and (3) assuming domestic trips alter the local infection rate. The simulated data resemble the arrival times and infected cases within the given ranges as of 1 June 2020 (see Supplementary Table 3). We conducted extensive sensitivity analysis and model validation using different combinations of parameters, as shown in Supplementary Note 2.
Evidence shows that human mobility determines arrival times 12,17,43 and infected cases 38 when there is only one OL. However, these approaches are not suitable for the presence of multiple OLs because it is unclear how each OL contributes to the arrival time and infected cases in each geographical area. Our approach compresses multiple OLs to a single one and calculates the corresponding effective mobility using the series and parallel connection law for global disease transmission. As demonstrated in Fig. 2c, d, country distancing correlates with the simulated arrival times and the logarithm of the simulated infected cases with a coefficient of determination R 2 > 0.7 regardless of the number of OLs at time t. In addition, R 2 for country distancing is consistently far above its value for effective distance, which is ≤0.3. For a robustness check, we also test the predictability of country distancing for different criteria of selecting OLs in Supplementary Figs. 1-5. All the result indicates that instead of effective distance, country distancing is an excellent predictor of global disease dynamics, especially for the evolving OLs. Thus, for each time t, we formulate functional linear relationships between country distancing and remaining arrival times T m (t) [T m (t) = T m − t] and new infected cases I m (t + 1) − I m (t) as: where v(t) and u(t) are the slopes (speeds), as shown in Fig. 2e, f, representing the rates of change of arrival times/infected cases relative to country distancing. The different OLs and temporal infection rates exhibit variability in speeds and further estimating the effectiveness of travel restrictions (see Supplementary Notes 3 and 4).
Status quo of existing travel restrictions. As of 1 June 2020, 625 travel restrictions have been imposed worldwide. An overview of implemented travel restrictions in given geographical areas is shown in Fig. 3a. As shown in Supplementary Table 2 and Supplementary Fig. 6, 184 geographical areas imposed 476 entry bans from 21 January 2020 to 17 March 2020. These areas that imposed entry bans deny access to noncitizens at some specific geographical regions such as mainland China, South Korea, Japan, Iran, and the Schengen Area. From 11 March 2020, when approximately half of the countries were infected worldwide, entry bans were not sufficient to lower the risk of coronavirus importation from the infected regions 44 . Consequently, 111 geographical areas imposed the global travel ban to prevent overseas travelers from entering their areas except for their residents, and 38 geographical areas imposed full/national lockdowns to deter people entering and exiting their countries.
Measuring effectiveness and efficiency of travel restrictions using country distancing. According to the International Civil Aviation Organization 42 , 90% of passenger seats were lost in the second quarter of 2020. As illustrated in Fig. 3b, for the area n that implemented the travel restrictions (see Eqs. (5)- (7)), (1) the entry ban reduces the passenger influx from banned areas to n with strength α = 90%; (2) the global travel ban reduces the passenger influx from all neighboring areas to enter n with strength β = 90%; (3) the lockdown reduces the passenger influx entering/leaving n with strength γ = 90% for full lockdown. In other words, travel restrictions induce a reduction in passenger influx, leading to larger country distancing (named country distancing increase, see Eqs. (8)- (15)) in Fig. 3c. We map the country distancing increase at each geographic area to the changes in arrival times (named as arrival time delay (ATD)) and changes in infected cases (named as infected case reduction (ICR)) according to the linear relationships in Eq. (2). Given the slopes of the linear correlations, namely, v(t) and u(t) in Eq. (2), the effectiveness of travel restrictions in slowing disease dynamics could be quantified by the world's average ATD and total ICR with the loss of air travel caused by them. Through all existing travel restrictions, 36.3% of passenger influx is reduced from GMN, the arrival times of the disease in the world are delayed by 18.56 (95% CI, 15.84-21.28) days on average, and infected cases are reduced by 13,186,045 (95% CI, 3,992,055-40,148,831) cases. However,~63.2% of travel restrictions are ineffective, resulting in zero country distancing increase. Figure 4a, b shows the effectiveness and efficiency of the remaining 36.8% travel restrictions by representing their ATD/ ICR against the cost of losing passenger influx. The means of lost passenger influx (vertical line) and ATD/ICR (horizontal line) divide the travel restrictions into four blocks, i.e., effective and efficient block (top left), effective and inefficient block (top right), ineffective and efficient block (bottom left), and ineffective and inefficient block (bottom right). In addition to the lockdowns imposed by mainland China, the lockdowns imposed by Spain (@ES), New Zealand (@NZ), and South Africa (@ZA) produce crucial ATD for the world. Concurrently, lockdowns established by Italy (@IT) and Spain (@ES) as well as the global travel ban enforced by Turkey (*TR) carry significant ICR for the world. However, these extreme travel restrictions also aggravate the loss of passenger influx. Several travel restrictions, e.g., the entry ban imposed by the United States (US-CN), Hong Kong (HK-CN), and Italy (IT-CN) on mainland China produce comparable ATD or ICR but with considerably less loss of passenger influx. Moreover, 505 (80.8%) travel restrictions generate <0.01 days of ATD and <1000 cases of ICR for the world, suggesting the ineffectiveness of travel restrictions.
Quantification of the impact of travel restrictions through country distancing increase. Figure 4 indicates that entry bans to the OLs are most effective and efficient. The areas implementing such entry bans and their descendants in the shortest path tree are distant from the OLs. As shown in Fig. 5a, by implementing an entry ban on mainland China on 31 January 2020, the United States increased country distancing by 40.96 in total. This means that the entry ban produces 8.10 days of ATD for the 28 areas, which are descendants of the United States (shaded area). By imposing a global travel ban on 21 January 2020 (see Fig. 5b), North Korea's country distancing rises by 2.30 (13.53 days of ATD, 0 cases of ICR) but induces no country distancing change for other countries, because North Korea is a leaf node. The findings suggest that the travel restrictions imposed by the areas with more descendant areas in the shortest path tree are more effective in distancing the world from coronavirus importation risk. The OLs, which are the sources of the shortest path tree, have the most significant influence. By imposing the national lockdown on 8 February 2020 (see Fig. 5c Fig. 5d-f, the distance to OLs is reduced. Substantial intervention efforts cause relatively smaller country distancing increases and fewer health benefits in delaying infection. For example, see Fig. 5d, e, the entry ban imposed by the United States to Japan on 1 March 2020 and the global travel ban imposed by New Zealand on 19 March 2020 both leads to <0.07 ATD and <2000 ICR. The same phenomenon is observed in Italy, one of the 38 sources in the shortest path tree, which imposed a national lockdown on 23 March 2020 (see Fig. 5f). This lockdown led to a 3.75 increase in country distancing. The lockdown also produces 0.007 (95% CI, 0.006-0.008) days of ATD per geographical area and 13,590 (95% CI, 10,569-17,191) cases of ICR in total, accounting for nearly 0% of ATD and 0.1% of ICR generated by all travel restrictions. Our findings indicate that two factors determine the effectiveness of travel restrictions. One is the geographical areas' position that imposed travel restrictions on the shortest path tree; another is the implementation date. Insufficient consideration of both factors leads to ineffectiveness and further inefficiency of travel restrictions.
Geographic areas' contributions and gains. Geographical areas implemented different travel restrictions. By integrating the health benefits of travel restrictions enforced by the same area, we find that mainland China (CN) has a dominant contribution followed by Australia (AU), the United States (US), New Zealand (NZ), the Netherlands (NL), and Russia (RU), as shown in Fig. 6a (also see Supplementary Fig. 10). This finding suggests that mainland China is the most influential area in mitigating the spread of COVID-19. On the other hand, Italy (IT), Germany (DE), Hong Kong (HK), Turkey (TR), and Taiwan (TW) contribute a great number of ICRs by imposing travel restrictions. Optimal and coordinated travel restrictions. We find that most of the existing travel restrictions are ineffective for two reasons: (1) the travel restrictions are imposed by geographical areas in an uncoordinated manner out of self-interest, failing to contribute to  Fig. 3. The two-letter codes for geographical areas can be found in Supplementary Table 5. the global public interest; (2) the sole travel restriction is not enacted in an optimal time and at optimal locations for the most significant self-interest. Furthermore, these inefficient travel restrictions have created a substantial unnecessary loss of passenger influx, ultimately damaging the global economy and social stability 25 . These findings prompted us to design the strategic plans for when and where to impose each travel restriction tailored to the real-time national context. Specifically, we formulate a bi-objective optimization problem of maximizing travel restrictions' the effectiveness (country distancing increase) and efficiency (the loss of airline passenger influx) (see Eq. (16)). If governments worldwide coordinately implemented the optimal travel restrictions, the infection could be dramatically delayed and avoided.
Using the Non-dominated Sorting Genetic Algorithm (NSGA-II), we obtain non-dominated solutions for each entry ban and present the approximate optimal solution, which has the largest country distancing increase in Fig. 7. Our numerical results show that the optimal and coordinated travel restrictions significantly outperform the existing travel restrictions in the three features, i.e., lost passenger influx, the average ATD, and the total ICR. As shown in Fig. 7a-c, the coordinated travel restrictions reach of an The star area's travel restriction directly increases its descendants nodes' distance from the OL. As the number of OLs grows in d-f (multiple OLs merge to the root), the geographical areas' country distancing is reduced. The colors of the areas in the trees correspond to the colors on the corresponding maps. See Supplementary Fig. 12 for the shortest path tree before any travel restrictions are implemented. Unlike the existing travel restrictions suggesting that mainland China contributes most ATD and ICR for the world, the coordinated travel restrictions work as a whole-of-government and whole-of-society approach with many geographical areas contributing substantially to an increased ATD and ICR for the world 45 . As shown in Fig. 7d, mainland China occupies small portions of the ATD and ICR contributions, which decline to 1.75% and~0% from 58.26% and 49.43%, respectively, when a coordinated approach is adopted. Concurrently, other geographical areas contribute more when using a coordinated approach. For example, the United States' portion of ATD and ICR contributions could rise to 9.36% and 2.95% from 4.89% and 0.46%, respectively.

Discussion
We developed a country distancing method that captures the global disease dynamics when OL evolve. Our method enables us to systematically quantify the effectiveness and efficiency of travel restrictions (i.e., entry bans, global travel bans, and lockdowns), which delay or prevent infections at the cost of loss of air travels. Our analysis confirms findings from existing studies 4,35,39,46 , which conclude that travel restrictions are more effective at the early stage of the pandemic. Our study reveals that the coexistence of multiple OLs substantially weakens the effectiveness of travel restrictions. In addition, travel restrictions' effectiveness is related to the geographical areas' capability to cut possible disease propagation paths from existing OLs to other countries. By maximizing the health impact of travel restrictions, in other words, increasing the country distancing from OLs and minimizing the loss of air travel, we find that the well deployment of entry bans to OLs via global joint efforts as early as possible is sufficient to fight effectively against COVID-19. The optimized and globally coordinated travel restrictions enable the sustainable suppression of transmission at a low level 47,48 , without the use of radical approaches, such as the global travel bans.
Three main limitations can lead one to overestimate/underestimate travel restrictions' effectiveness: (1) incomplete and biased travel restrictions data; (2) homogeneous assumptions on the strengths of different travel restrictions; (3) ignorance that the combined effect of travel restrictions and local anti-contiguous policies 6 , such as social distancing policy, work from home, and school closure, has a substantial effect on reducing the rate of global spread 21,49 . Understandably, these three limitations could be overcome with more available data of the daily GMN. Furthermore, COVID-19 undertesting and underreporting 50 would overvalue the impact of travel restrictions on slowing infection. Nevertheless, the approach in this study and its implications may help to control the spread of COVID-19.
The key advantage of the country distancing approach is that it captures the global diffusion pattern despite the heterogeneous responses of governments to the pandemic and varying OLs. The idea of quantifying the health impact of travel restrictions is useful in the ongoing pandemic and future diseases. Although many geographical areas attempt to ease their travel restrictions by escalating the testing and mandatory quarantines for the incoming travelers, it is clear that addressing how people/infected people move is key to controlling the disease. Pinpointing the significance of controlling OLs provides insights to properly distribute and administer a COVID-19 vaccine when it becomes available. In summary, the use of country distancing approach complements existing studies highlighting the significance of mobility from both the OLs and hubs in GMN. Furthermore, this study also seeks to balance optimally positive health effects of travel restrictions with their negative impact on individuals' free movement. Knowing that the COVID-19 pandemic is more than a health crisis and may last until 2022 51 , geographical areas are continuously enduring the coronavirus importation risk from other infected areas and enduring social instability. We emphasize the necessity of a global joint effort to implement travel restrictions when it seems impossible to curb the spread of COVID-19 with isolated actions. Although limitations exist, this study indicates that whole-of-government and whole-of-society approaches are necessary to fight against coronavirus and strengthen pandemic preparedness in the future 45 .

Methods
Model with multiple outbreak locations. We adopt the meta-population SIR model 17,52,53 to simulate the spread of disease in GMN G when OLs are changing. In the meta-population SIR model, each geographical area n has a population size of Ω n . The evolving disease within the population is governed by three states (i.e., susceptible s n ¼ S n Ω n , infectious i n ¼ I n Ω n , removed r n ¼ R n Ω n ), and the disease evolution between populations is described by the travel influx coupling term p mn . For each set of OL N I (t): _ s n ¼ ÀaðtÞs n i n σði n =εÞ þ cðtÞ∑ m≠n p mn ðtÞðs m À s n Þ _ i n ¼ aðtÞs n i n σði n =εÞ À bðtÞi n þ cðtÞ∑ m≠n p mn ðtÞði m À i n Þ _ r n ¼ bðtÞi n þ cðtÞ∑ m≠n p mn ðtÞðr m À r n Þ with the initial condition given as: from t to t + 1. We show the definitions of the parameters as below:  Fig. 3 and the two-letter codes for geographical areas can be found in Supplementary Table 5. • N I ðtÞ ¼ fnj8n 2 N; I o n ðtÞ ≥ θg are the defined OL set, consisting of the geographical areas whose reported cases I o n ðtÞ are greater than the threshold θ. The time range for N I (t) is from 8 December 2019 to 1 June 2020. For the initial outbreak date t 0 = 2019 − 12 − 08, mainland China is the OL with I o n ðt 0 Þ ¼ 1.
• a(t) is the infection rate, and b(t) is the remove rate, which is defined as the sum of recovery rate and death rate.
• cðtÞ ¼ ∑ m;n F mn ðtÞ ∑ n Ω n is the mobility rate.
• p mn ðtÞ ¼ F mn ðtÞ ∑ k F kn ðt 0 Þ is the fraction of travel influx, where the influx matrix F mn (t) is altered by travel restrictions; see "Passenger influx reduction".
η (η = 4) is the gain parameter for the threshold ε. (1) An entry ban leads to a α decline in passenger influx from banned areas m to n s , i.e., F s n s m ¼ F sÀ1 n s m ð1 À αÞ: ð5Þ (2) A global travel ban results in a reduction in passenger influx from neighboring areas m to n s by β, i.e., F s n s m ¼ F sÀ1 n s m ð1 À βÞ: ð6Þ (3) A lockdown reduces passenger influx by γ from neighboring areas m to n s and passenger influx from n s to neighboring areas m, i.e., F s mn s ¼ F sÀ1 mn s ð1 À γÞ; F s n s m ¼ F sÀ1 n s m ð1 À γÞ: The other weighted links that are not influenced by s th travel restrictions remain the same, i.e., F s mn ¼ F sÀ1 mn . Notably, F 0 = F. Thus, ∑ m;n F s mn À F 0 mn is the lost passenger influx until the s th travel restriction is imposed.
Country distancing increase. Based on the influx matrix F s and the set of OLs N s I , we could measure the country distancing of geographical area m when s th travel restriction is imposed with: where d mjðn i ;F s Þ is the effective distance from source n i to destination location m based on the influx-fraction matrix P s . Here, P s mn ¼ F s mn ∑ k F kn , which is derived from F s and F. We know that travel restrictions increase the country distancing by decreasing passenger influx in F s , whereas the number of OLs N s I decreases country distancing by promoting importation risk from multiple OLs. To better understand the impact of travel restrictions on country distancing, we exclude the influence of OLs D mjN s I : where: Given that travel restrictions, which are collected until 1 June 2020, continuously reduce the airline passengers, the country distancing difference with and without the s th travel restriction is always non-negative. In detail: Because ΔD s mk ≥ 0, we call the non-negative difference ΔD s mk as the country distancing increase at geographical area m caused by s th travel restriction imposed by area k.
Arrival time delay and infected case reduction. Country distancing is a good predictor for the arrival times and cumulative infected cases, as illustrated in Eq. (2). Given the country distancing increase (ΔD s mk ) at area m caused by travel restrictions s implemented by area k, we could map the country distancing increase (ΔD s mk ) to the delayed infection-ATD and ICR. Specifically, ATD, i.e., the difference in arrival times with and without s th travel restriction is denoted as: where v(t s ) is the speed for arrival times at time t s when s th travel restriction is implemented. It should be noted that the ATD is only applicable to uninfected areas. See Algorithm 1 in Supplementary Note 3 for the details. On the other hand, ICR, i.e., the difference of newly infected cases with and without travel restrictions for area m at time t is: where I ? m ðtÞ ¼ I o m ðt þ 1Þ À I o m ðtÞ is the observed new infected increase in area m from day t + 1 to day t and D ◊ m ðtÞ ¼ ∑ s2fsjt s ≤ tg ΔD s mk is the accumulative country distancing increase due to the travel restrictions implemented by time t. Thus, as of the end date (t S = 1 June 2020), the accumulative ICR at area m is: by assuming that daily ICRs in geographical areas are independent. For s th travel restrictions implemented by area k, the ICR at area m is: Given that travel restrictions continuously reduce the airline passengers and continuously distance geographical areas, the country distancing increase accumulates over time. Simultaneously, the ICR accumulates over time and grows exponentially as the country distancing increase accrues. However, as time progresses, more or less OLs are present, an the intertwined effect of OLs disables existing travel restrictions. See Algorithm 2 in the Supplementary Note 3 for the details on obtaining ICR when multiple OLs are present. Substituting the upper and lower bounds of the 95% confidence interval for v(t) and u(t) in Eqs. (12) and (13), we obtain the upper and lower limits for ATD and ICR. Thus, for s th travel restriction implemented by area k, we could calculate its global average ATD as ∑ m ∑ k ΔT s mk M and total ICR as ∑ m ∑ k ΔI s mk . For each geographic area k, we could calculate its total contribution of ATD as ∑ s ∑ m ΔT s mk M and the total contribution of ICR as ∑ s ∑ m ΔI s mk . Furthermore, the gain of ATD in geographic area k is ∑ s ∑ m ΔT s km M and the gain of ICR in geographical area k is ∑ s ∑ m ΔI s km . Because the two matrices ΔT s mk and ΔI s mk are asymmetric, i.e., ΔT s mk ≠ΔT s km and ΔI s mk ≠ΔI s km , such that countries' contributions differ from their gains.
Bi-objective problem. To address the ineffectiveness and inefficiency of the existing travel restrictions, we formulate the coordinated and targeted travel restrictions as a bi-objective problem of maximizing their influence on delaying infection (maximizing the entry bans' country distancing increase) and minimizing the loss of airline passenger influx in GMN: For the s th entry ban implemented by area k, which reduces the passenger influx from arbitrary area n to area k, we find its corresponding optimized solution θ s 0 (θ s 0 2 Θ s 0 ) by selecting one airline link from the GMN. The set Θ s 0 is the solution set for the s 0th entry ban. Then each airline link ðn; kÞ 2 θ s 0 follows the formulation F s 0 kn ¼ F s 0 À1 kn ð1 À αÞ. The optimal solution should ensure that the lost passenger influx is minimal while maximizing the global country distancing increase. To solve the problem of minimizing the loss of passenger influx and maximizing country distancing increase, we adopt the NSGA-II 54 , a well-known fast sorting and elite multi-objective genetic algorithm. This algorithm can find the solutions that are not dominated by any other solutions and are closer to the true Pareto optimal front in the solution space. Given that we have found that effective travel restrictions prevent air travels from OLs, we prioritize the solution of entry bans to OLs. The procedure for generating non-dominated fronts follows the algorithm proposed in the literature 54 . Based on experiments, we finally choose a population of size 100, a crossover probability of 0.5, and a mutation probability of 0.5 to solve the problems of balanced travel restrictions. The algorithm terminates after it runs 1000 generations.

Data availability
All data needed to evaluate the paper's conclusions are presented in the Supplementary Note 1. Additional data related to this paper may be requested from the authors.