Policies for Easing COVID-19 Pandemic Travel Restrictions

During the COVID-19 pandemic, many countries implemented international travel restrictions that aimed to contain viral spread while still allowing necessary cross-border travel for social and economic reasons. The relative effectiveness of these approaches for controlling the pandemic has gone largely unstudied. Here we developed a flexible network meta-population model to compare the effectiveness of international travel policies, with a focus on evaluating the benefit of policy coordination. Because country-level epidemiological parameters are unknown, they need to be estimated from data; we accomplished this using approximate Bayesian computation, given the nature of our complex stochastic disease transmission model. Based on simulation and theoretical insights we find that, under our proposed policy, international airline travel may resume up to 58% of the pre-pandemic level with pandemic control comparable to that of a complete shutdown of all airline travel. Our results demonstrate that global coordination is necessary to allow for maximum travel with minimum effect on viral spread.

1 Models and methods 1

.1 Local travel model
We first consider a local epidemiological model as in Warne et al. (2020) (1). In this local model, for each country, at a given time, its population's status is divided into 6 mutually exclusive compartments: susceptible (S), undetected infected (I), active confirmed (A), confirmed recovered (R), confirmed deceased (D) and unconfirmed recovered (R u ). Its dynamic states evolve as: where α is the transmission rate, γ is the identification rate, β is the recovery rate, and δ is the death rate. Suppose that for a given country the status of its population at time t is X(t) = [S(t), I(t), A(t), R(t), D(t), R u (t)], and θ = (α, β, δ, γ) represents the parameter of the statistical model for the country. Using the tau leaping method by Gillespie (2001) (2), the status of its population at time (t + τ ) evolves as X(t + τ ) = X(t) + 5 i=1 Y j h j (X(t))τ ν i .

Global travel model
Our global epidemiological model model is built based on the local model by utilizing travel flow data as follows. For a given country i, suppose the status of its population at the end of and the parameter of the statistical model for this country is θ i = (α i , β i , δ i , γ i ). On day t, the epidemic state in country i is updated via two steps. First, the state evolves based on country i's internal population. Second, the state evolves based on external factors, here the inflow of airline travelers from other countries and the outflow of airline travelers to other countries.
We consider changes due to internal effects first. For country i, the transition from t − 1 to t is characterized by the shift from X i (t − 1) to X i (t) = [S i (t), where and Y j,i (t − 1), j = 1, · · · , 5, are Poisson distributed with rates and P i (t − 1) is the size of the population in country i on day (t − 1).
The travel data specify how many new individuals enter the country on day t from each of the disease states. The current state is updated as represents the six compartments of people entering the country on day t and f out i (t) represents the six compartments of people leaving the country on day t. Due to temperature checks and other approaches for screening travelers, we assume that all active confirmed cases are unable to travel. We also assume that deceased individuals do not travel between countries. Consequently, the compartments in f in i (t) and f out i (t) only include four of the six disease states: susceptible (S), undetected infected (I), recovered confirmed (R), and recovered unconfirmed (R u ). Travelers in the recovered confirmed (R) and recovered unconfirmed (R u ) states do not impact the epidemiological state of destination population. However, data on all four categories is not readily available. While each country keeps track of the total number of confirmed recovered each day, they do not necessarily keep track of how many of them leave the country. Therefore, we take a conservative approach and assume that each traveler either belongs to the S category or the I category, meaning travelers bring some potential risk when they arrive in a new country as undetected infected will likely spread the disease and susceptible individuals reduce population immunity and can proliferate disease spread. In where I in i (t) and I out i (t) are the number of undetected infected that enter and leave country countries. We model the number of undetected infected people who are leaving country i for country j at day t using a multinomial distribution with probabilities based on travel network data. In other words, {I out ij (t)} 1≤j =i≤n ∼ M (I out i (t), {p ij (t)} 1≤j =i≤n ), where p ij (t) = T out ij (t) T out i (t) and T out ij (t) is the total number of travelers leaving country i for j at day t. Let us denote S out ij (t) as the number of susceptible people who travel from country i to country j at day t. Then . Therefore, at the end of day t, the six states for country i are updated

Travel regulation policies
Our goal is to find the value p so that the number of undetected infected I + (1), I + (2), · · · , I + (T ) stay below a given threshold c. More specifically, we consider two types of regulation, an average control policy and a probability control policy described below: 1. Regulation in terms of average control, where we find a proportion p such that the average number of undetected cases each day in the next T days stays below a threshold c, i.e., E(I + (1)), E(I + (2)), · · · , E(I + (T )) < c, where E denotes the expectation.
2. Regulation in terms of probability control, where we find a proportion p such that the probability of undetected cases each day in the next T days being lower than a threshold c is at least at π, i.e. P (I + (1) < c, I + (2) < c, · · · , I + (T ) + < c) ≥ π.
The following lemmas gives us the proportion p that satisfies the above requirements.
In other words, we need a p that satisfies: p < (c/m) 1/k ψ − 1, ∀k = 1, · · · , T , or p < (c(1−π)/I(0)) 1/k ψ − 1, ∀k = 1, · · · , T . In conclusion, for a probability control level π and a threshold c in the next T days, the required p is Remark: The probability control policy is more conservative than the average control policy.
Under the same threshold c, the difference between the two policies is the factor (1 − π) in the 8 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted April 20, 2021. ; https://doi.org/10.1101/2021.04.14.21255465 doi: medRxiv preprint numerator of the probability control strategy. If we want the probability control to have at least 0.9 of the threshold c of average control, this factor becomes 0.1. As a result, the proportion of p in probability control is much smaller than the proportion of p in average control. If we want to use the probability control with a probability of at least 0.9, we need to set up the threshold c in probability control higher than the threshold c in average control to make sure that our p is non-negative.
Example 1. Here we give one example of using the average control policy to regulate the travel.
For simplicity, we consider a small world with only three countries, with the following initial states and true parameter values: and {p 31 (t)} t=1,··· ,7 that can regulate airline travel from country 2 to country 1 and from country 3 country 1 such that for the next T = 7 days, the number of undetected infected cases in the arriving country will not exceed c = 70 cases on average per day. Applying Lemma 1a, we can find parameter p for country 1 as p = min (c/I 0 ) 1/k /ψ − 1 k=1,··· ,7 = 0.035.
So the sequence of the number of undetected infected imported cases that country 1 can 9 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted April 20, 2021. ; https://doi.org/10.1101/2021.04.14.21255465 doi: medRxiv preprint accept each day as (1 + p)ψ i I(0) i=1,··· ,7 = (2, 2, 2, 2, 2, 2, 2). Because country 1 has two "neighbors", so country 1 can accept about 1 undetected infected from each neighbor each day during the regulation period. The next step is to predict, for each day, the number of undetected infected travelers from countries 2 and 3 that can enter country 1 for the next 7 days when full travel is allowed. We get these numbers by simulating data given the true parameters under the fully open scenario. We full travel is allowed. The final step is obtaining the regulation sequence that country 1 can allow country 2 and country 3 to enter its border. The regulation sequence that country 1 can allow for country 2 to enter its border during the regulation period is obtained by dividing the number of daily undetected infected cases that country 1 can tolerate from country 2 by the daily estimated number of imported undetected infected cases from country 2 if full travel is allowed. Notice that if the daily proportion is greater or equal to 1, we set it to 1. Repeat the same procedure, we can also find the regulation sequence that country 1 can allow for country 3 to enter its border.
Following the above steps, we can find that in the next 7 days, the regulation sequence of proportions of people who can move from country 2 to country 1 is (0.103, 0.076, 0.051, 0.035, 0.022, 0.014, 0.009), and the sequence of proportions of people who can move from country 3 to country 1 is (1, 1, 1, 1, 1, 1, 1). Compared to the fully open scenario, using the average control approach with the threshold of 70 cases during the 7-day regulation period, about 6.03%

10
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted April 20, 2021. ; of travelers from country 2 are allowed enter the country 1, and all travelers from country 3 are allowed to enter country 1. Overall, the volume of inbound travelers in country 1 is about 88.64% of the normal levels.
In practice, the value of the transmission rate α varies over time, and we therefore provide an additional lemma that generalizes Lemma 1 to address the aspect of varying α.

The average control requirement is satisfied if
2. The probability control requirement is satisfied if Proof: 1. Average control. Follow the same argument as in the proof of Lemma 1, we have Repeating the argument until we reach day T yields : E(I(T ) + ) < (1 + p) T ψ T max I(0). Therefore, we want to find a p such that (1 + p)ψ max I(0), (1 + p) 2 ψ 2 max I(0), · · · , (1 + p) T ψ T max I(0) < c. Hence the value p that 11 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted April 20, 2021. ; https://doi.org/10.1101/2021.04.14.21255465 doi: medRxiv preprint satisfies the requirements for average control is

Let us denote
. Similar to Lemma 1, the sequence of I * (1), I * (2), · · · , I * (T ) forms a non-negative supermartingale sequence. Hence, following the same arguments as in Lemma 1, we have for a probability control level π and a threshold c in the next T days, the required p is p = min (c(1 − π)/I(0)) 1/k ψ max − 1 k=1,··· ,T

Computational
There are many variants of ABC, but they are all based on a comparison of observed and simulated data, which in most cases requires specification of data summary statistics, a distance measure, and a scalar distance threshold . The most basic ABC algorithm, the so-called acceptreject method, starts by simulating a parameter value from a prior distribution and then uses the model, given this parameter value, to generate one realization of data. If the distance between the summary statistics for the observed data and the summary statistics for the simulated data is less than or equal to , the sampled parameter value is retained; otherwise, it is discarded.
The collection of accepted parameter values constitutes a sample from an approximation of the 12 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.  (3)(4)(5)(6). In this paper, we use the variant from Drovandi and Pettitt (2011) (6), called replenishment ABC (RABC). For its implementation, we use the R package protoABC from Ebert (2020) (7). This package is very flexible as the users can employ any priors, generative models, and distance functions.
A commonly used distance is the Euclidean distance due to its simple form. In our problem setting, this Euclidean distance L S(Data (i) ), S(Data) can be written as:  (8)(9)(10)). This is why we also consider the following weighted Euclidean distance, with weights given by the inverse variances:

13
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 20, 2021. ; where {sd(A(t)), sd(R(t)), sd(D(t))} t=1,··· ,T are the prior predictive standard deviations of A(t), R(t), D(t) at each time step. These standardization values are obtained by first generating N (N being large) parameters values θ = (α, β, δ, γ) from the prior π(θ), and then generating one realization of simulated data for each of them. The standard deviation at each time step is then calculated based on these N simulated data, giving {sd (A(t)), sd (R(t)), sd (D(t))} t=1,··· ,T .
In our simulation study, we choose N = 5000.
To improve the predictive quality of ABC algorithms, we can also use additional information from parameter estimates to make new summary statistics. Under our model assumptions, we can learn additional knowledge about how parameters can be estimated. We are thus going to include as summary statistics estimates of our epidemiological parameters. For simplicity, we first limit our discussion to the local model, where each country is considered separately. The choice of the distance for the global model will be discussed in Section 1.5.
Under our model assumptions, we have: So for a given sequence of known {A(t)} t=1,··· ,T , the sequence of } t=1,··· ,T has β as common mean value. Therefore, we can use its median value to estimate β. If our algorithm generated a reasonable θ (i) , then the data generated by θ (i) should also gives us a sequence } t=1,··· ,T can be used to estimate the death rate δ, and adding the term |median{ D(t)−D(t−1) improve the estimation of δ.
We now try to learn the transmission rate α under our model assumptions. We have S(t) = 14 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 20, 2021. ; https://doi.org/10.1101/2021.04.14.21255465 doi: medRxiv preprint S(t − 1) + Poisson α S(t−1)I(t−1) P , ∀t = 1, · · · , T , where P is the total population size of the country. So for a given S(t − 1), Unfortunately, S(t − 1) and I(t − 1) are hidden states and not available in our data. To use the above strategy to improve the estimation of α, we need to reconstruct these hidden states based on the available data {(A(t), R(t), D(t))} t=1,··· ,T . Because our model is stochastic, all values would change each time we rerun the model. However, based on the available data {(A(t), R(t), D(t))} t=1,··· ,T we can reconstruct the mean realization that adopts these three states. Let us denote U (t) the total number of confirmed cases at time t, ∆U (t − 1) the number of new confirmed cases at time t as, and ∆R u (t − 1) the number of new undocumented recover cases at day t. Note that From the local model we have: Equation (5) . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
Using (5), we obtain ∆R u (i) and (6), the average value of R u (t) can be reconstructed as Equations (5) and (7) tell us that based on the available data {(A(t), R(t), D(t)} t=1,··· ,T if the identification rate γ and the recovered rate β are available to us, then we can reconstruct the average realization of I(t) and R u (t), ∀t = 1, · · · , T − 1. As a result the average realization category at time t can also be reconstruct as where P is the country population. Overall, based on the observed data {(A(t), R(t), D(t))} t=1,··· ,T , suppose that the identification rate γ and the recover rate β are available to us. The average realization that adopts 16 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 20, 2021. ; https://doi.org/10.1101/2021.04.14.21255465 doi: medRxiv preprint {(A(t), R(t), D(t))} t=1,··· ,T can be reconstructed up to time T − 1 as Therefore α can be estimated as the median of the sequence Unfortunately, β and γ are not known in advance and need to be estimated.
We therefore use the testing argument to recover the average realization and estimate α. This argument is based on the following observation. In step 1 of the ABC algorithm, a parameter value θ (i) = α (i) , β (i) , δ (i) , γ (i) is generated and available to us. If γ (i) , β (i) are correctly specified as γ, β of the underlying true parameter θ = (α, β, δ, γ). We would expect that the median ,··· ,T should also give us a good estimator for the underlying true α value. Therefore the distance |median{ S(t−1)−S(t) , S(Data) should be close to 0. On the other hand, if the generated parameters γ (i) , β (i) are far away from the underlying true param- Based on this observation, we then add the term in L S(Data (i) ), S(Data) to improve the estimation for α.

17
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 20, 2021. ; https://doi.org/10.1101/2021.04.14.21255465 doi: medRxiv preprint Finally, our proposed distance for the model is designed as follows:

Marginal approach to parameter estimation
In the following, we will discuss how to use ABC to estimate parameters in each country for our global model. The challenge of using ABC to estimate the global model parameters is that many parameters need to be estimated. Therefore directly using ABC to estimate all the parameters for all countries at once may result in very unstable parameter estimations and will be extremely computationally intensive. We propose a marginal estimating approach to estimate each country parameter for the global model separately while still taking into account the travel data.
For simplicity, let us first consider a given country m with the represented parameter as the number of active confirm cases, accumulated recover confirmed, and accumulated death confirmed at country k on day t, respectively. We denote T (t) = [T ij (t)] i,j=1,··· ,n is the global traveling matrix at day t, where T ij (t) is the number of travelers from country i to country j at day t. Notice that T ij (t) = 0 if i = j, ∀t ∈ 1, · · · , T . With the available data from the global model , D k (t)} k=1,··· ,n;t=1,··· ,T and the travel data {T (t)} t=1,··· ,T , we need to estimate θ i .
Before introducing our estimation procedure, we rewrite how our global model evolved for 18 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 20, 2021. ; https://doi.org/10.1101/2021.04.14.21255465 doi: medRxiv preprint the country m at day t and (t + 1): Step 1a The internal pandemic evolves at day t: The internal pandemic situation in the country Step 1b The external pandemic added at day t: From the travel data, Step 2a The internal pandemic evolves at day (t + 1): The internal pandemic situation in the country evolves from X + m (t) to X m (t + 1) = [S m (t + 1), and h 5 (X m (t)) = β m I + m (t).

19
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 20, 2021. ; https://doi.org/10.1101/2021.04.14.21255465 doi: medRxiv preprint As shown in the model's evolving process, step 1b and step 2b make the global model behave differently from the local model. Therefore, if we can estimate quantities S in for each day t, then we can use the marginal approach to estimate each country's parameters separately. The two quantities S out m (t), I out m (t) come from inside the considered country m. Therefore the amount can be calculated during the data generation process of the ABC algorithm. Our main task now is estimating S in m (t), I in m (t).
We have S in We discuss the procedure for estimating I out jm (t). To estimate I out jm (t) we need to estimate the pandemic situation in country j at day t, i.e. we need to estimate: . Then based on these compartments and the travel data T jm (t), we can estimate I out . From the global model we have:

20
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 20, 2021. ; Since Y 2,j (t) ∼ Poisson γ j I + j (t − 1) , therefore, we have Similarly we also have From (8) and (9), we have Applying (10) for t = 1, · · · , T − 1, we have the sequence of relationships: ) . Therefore, based on the available data at country j as {A j (t), R j (t), D j (t)} t=1,··· ,T , we can approximate the average realization of the sequence of undetected infected people in country j up to time (T − 1) by I j (0), I + j (1) = ∆U j (T −1) . In addition, we also have R u j (t) = R u j (t−1)+Y 5,j (t−1), where Y 5,j (t−1) ∼ Poisson βI + j (t− 1) . Therefore, } t=1,··· ,T can be used to approximate β. We denote this median value asβ.
∆R u j (i) and (11) tells us that the average value of R u (t) at 21 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
So the average realization of the pandemic at country j which adopts {(A j (t), R j (t), D j (t)} t=1,··· ,T as its active confirmed, recover confirmed and death confirmed can be reconstructed up to This procedure of estimating the average realization of a given country j based on the available data {(A j (t), R j (t), D j (t))} t=1,··· ,T is completed. As a result, this gives us the estimation of S out jm (t) and I out jm (t). This means we can estimate S in m (t), and I in m (t). Therefore, the underlying true parameter θ m in a given country m can be approximated marginally by using the above estimating procedure.
We now discuss the proposed distance when estimating θ m in a given country m by ABC marginally. Following the same argument as in Section 1.4, instead of using the commonly used Euclidean distance to estimate θ m we first need to standardize each sequence, then we try to learn each parameter under our model assumptions. From the available data    } t=1,··· ,T | to improve the estimation for the death rate δ m . We discuss the transmission rate α m , at a time step t + 1, we have:

22
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 20, 2021. ; . Notice that S m (t + 1) = S + m (t + 1) − S in m (t + 1) + S out m (t + 1). The hidden states S + m (t) and I + m (t) can be reconstructed as above or by using the testing argument as above. We discuss the approach by using the testing argument.
From (8) we have: . Therefore the average realization of The average realization of the pandemic at country m which adopts {(A m (t), R m (t), D m (t)} t=1,··· ,T as its active confirmed, recover confirmed and death confirmed can be reconstructed up to time Similarly as above, in step 1 of the ABC algorithm, the parameter θ is generated and available to us. If γ } 1,··· ,T −1 would help to improve the estimation for α m . Finally, the proposed global distance in the calibrating step of the ABC algorithm is designed as follow:

23
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 20, 2021. ; where 2 Simulation studies 2.1 Simulation 1: Performance of different RABC distance metrics For our first simulation study, we limit our model to the analysis of only one country, i.e., we only use the internal model. We here demonstrate the impact of the choice of the distance in ABC algorithms and which one to choose in our epidemiological framework.
We simulate N = 200 sets of parameters and data, in an ABC fashion, by first simulating a parameter value from the prior and using it to generate data according to the model. We treat these N simulations as our test data set to assess how accurately the true parameters are recovered by ABC using various distance functions. The simulation proceeds as follows.
Step 1. Generating data and parameters: For i ∈ {1, · · · , N } (N large), we generate U (0, 1), and γ (i) ∼ U (0, 1). Based on the parameters and the stochastic model, we generate a data set Data (i) corresponding to θ (i) . If the generated data set Data (i) satisfies certain conditions making it sufficiently real-world like, such as having the number of confirmed accumulated deaths greater than 1% and lower than 30% of total confirmed cases and having the number of accumulated recovered cases at least twice the number of accumulated deaths, then we keep θ (i) as a true parameter value to be estimates and treat the generated data {A CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 20, 2021. ; https://doi.org/10.1101/2021.04.14.21255465 doi: medRxiv preprint real observed data. We repeat the process until we obtain 200 underlying true parameter values θ (i) and the corresponding 200 datasets {A ) is used as the best candidate for estimating the underlying true θ (i) .
Step 3. Evaluating parameter estimates: For each iteration i, i ∈ {1, · · · , 200}, we evaluate estimation accuracy in terms of the absolute bias, absolute relative bias, interquartile range, and coverage rate of the interquartile for each parameter α (i) , β (i) , δ (i) , γ (i) and its average. These accuracy measurements are defined as follows. For a given parameter α (i) the absolute bias is defined as |α (i) − α (i) |, and the absolute relative bias is defined as |α Similarly for β (i) , γ (i) , and δ (i) . Average absolute bias for all parameters is defined as

)/4 and average absolute relative bias for all parameters is defined as
For each parameter, we also calculate the interquartile range (IQR) of the posterior, denoted IQR (i) , which is the difference between the third and the first quartile of the resulting ABC posterior distribution. Furthermore, we check whether the IQR of the posterior covers the underlying true parameter value, which we use to calculate the coverage rate CR (i) . Then the average of the IQR for all four parameters and the average of the coverage rate is calculated to characterize the overall performance of the two distance metrics. Finally, the average over the 200 iterations of these accuracy metrics is calculated, which we use as our overall accuracy metrics for comparing the performance of the 25 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 20, 2021. ; https://doi.org/10.1101/2021.04.14.21255465 doi: medRxiv preprint two RABC distance metrics. Table S1 summarizes the different prediction accuracy measures for the two distances. This table shows that the distance we proposed increases the estimation accuracy in terms of relative bias. The two types of bias are much smaller compared to using Euclidean distance. We also observe that the IQR for the proposed distance is higher than the IQR for the Euclidean distance, but the proposed distance also yields more narrow IQR. This means that our proposed distance metric more frequently correctly bounds the true parameter values. The IQR is about 2.5 times smaller when using the proposed distance metric instead of Euclidean distance. Figure S1 shows boxplots of the IQR for the two distance metrics. is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

27
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

Simulation 2: Performance of local and global estimation
In this simulation study, we investigate the accuracy of three different estimation procedures for the global travel model consisting of three countries. (We limit this investigation to three countries for simplicity.) The first procedure uses a local approach with the Euclidean distance to estimate each country's parameters independently and ignores the travel between the countries.
We call this estimation procedure Euclidean local, and we use it as a benchmark to be compared with the other two approaches. Then we consider a global estimation procedure as discussed in Section 1.5 to estimate each country's parameters. Here we use two distance metrics, Euclidean distance and the distance proposed in Section 1.5. We call these estimation procedures Euclidean global and Proposed global, respectively. The simulation is set up as follows.
Step 1. Generating data and parameters: For i ∈ {1, · · · , N } (N large), we generate the j ∼ U (0, 1), and γ (i) j ∼ U (0, 1). Based on the parameters and the stochastic model, we generate a data set Data (i) corresponding to θ (i) . If the generated data set Data (i) satisfies the conditions described above for Simulation 1, we retain θ (i) and treat it as the underlying true parameter value; we also retain the data (t,j) }, for j = 1, 2, 3, and treat them as the observed data from these three countries. We repeat the procedure until we have 500 parameter values and their corresponding data sets {A For simplicity, we fix the initial condition of the six compartments in the model as and set the simulation period T = 84 days for all i. Each day, the number of outbound travelers from country j, j = 1, 2, 3, is drawn from a normal distribution with mean µ j = P j * 0.0003 and standard deviation sd j = 0.05 * µ j , where P j is the size of the population of country j.
Those outbound travelers will enter one of the neighboring countries with proportions that are proportional to the sizes (populations) of the target countries. For example, if there are n 1 people leaving country 1, the number of them entering country 2 is n 1 P 2 P 2 +P 3 and the number of them entering country 3 is n 1 P 3 P 2 +P 3 .
Step 2. Estimating parameters: For each iteration i, i ∈ {1, · · · , 500}, based on the sequence (t,j) } j=1,2,3 , we first naively use RABC with the local estimation approach and Euclidean distance to estimate θ (i) . Then we use RABC with the global estimation approach with the two distance metrics to estimate the underlying true parameters θ (i) . Thenθ is obtained as the median of the RABC posterior samples and is used to estimate the underlying true θ (i) j , for j = 1, 2, 3.
Step 3. Evaluating parameter estimates: For each iteration i, i ∈ {1, · · · , 500}, for each country, we evaluate the accuracy of our parameter estimates based on the absolute bias, absolute relative bias, interquartile range (IQR), and coverage rate of IQR for each parameter α (i) , β (i) , δ (i) , γ (i) and its average as in simulation 1. The final accuracy measurements are calculated by averaging the accuracy measurements across all three countries. When averaging accuracy measures over multiple countries, we consider two weighted averages, one having equal weights for all countries regardless of their population sizes and the other weighted based on relative population sizes. In the latter, the weights are P 1 P 1 +P 2 +P 3 for country 1, P 2 P 1 +P 2 +P 3 for country 2, and P 3 P 1 +P 2 +P 3 for country 3. Tables S2 and S3 show the overall accuracy of different estimation procedures using equal weights for each country (Table S2) and using population-based weights for each country (Table   S3). The two tables convey the same message: using a local approach to estimate the parameters 29 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 20, 2021. ; Table S2: Accuracy of Euclidean distance and our proposed distance when using RABC to estimate the four parameters of the global model. For simplicity, we consider a small world of just three countries with different population sizes; here each country has the same weight when computing the overall accuracy.

Accuracy
Estimation is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 20, 2021. ; Table S3: Accuracy of Euclidean distance and our proposed distance when using RABC to estimate the four parameters of the global model. For simplicity, we consider a small world of just three countries with different population sizes; here each country's contribution to the overall accuracy is weighted based on the size of its population.

Accuracy
Estimation is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 20, 2021. ; https://doi.org/10.1101/2021.04.14.21255465 doi: medRxiv preprint in the travel model is not appropriate. As shown in these tables, the Euclidean local estimation procedure yields the highest bias, largest interquartile range, largest 95 % percentile range, and lowest coverage. The performance is better for the Euclidean global procedure. As expected, the proposed distance, which takes into account the travel model, performs best of the three.

Simulation 3: Effectiveness of travel regulation
In this simulation study, we study the effectiveness of different travel regulation policies. We compare the percentages of people allowed to travel under each policy and the pandemic situation in the country adopting the policy. The simulation is set up as follows.
Step 1. Generating data and parameters: For i ∈ {1, · · · , N } (N large), we generate the and γ (i) j ∼ U ( , 1 − ). We chose = 0.001 to make sure that the generated parameters do not fall at the boundaries of the priors and cause the generation of atypical data. We also added some constraints to ensure the parameter values are reasonable by only keeping between 0.47 and 6.47 as reported for different regions around the world (11). To investigate the effectiveness of travel regulations, we use one more constraint to set the reproduction number R 0 in these 4 countries in 4 different zones, where country 1 has R 0 between 0.47 and 0.9 , country 2 has R 0 between 0.9 and 1, country 3 has R 0 between 1 and 1.1, and country 4 has R 0 between 1.1 and 6.47. The initial conditions of each country are generated randomly as (S (i) 10). Based on the parameters and the stochastic model, we generate a data set Data (i) corresponding to θ (i) . If the generated data set Data (i) satisfies the conditions described for Simulation 1 above, we keep θ (i) and treat it as the underlying true value of the 32 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

33
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 20, 2021. ; https://doi.org/10.1101/2021.04.14.21255465 doi: medRxiv preprint P-4 The country requires a 14-day quarantine for travelers from high-risk countries only, i.e., countries with the average number of active confirmed daily cases greater than 20 in 100000 people during the last 2 weeks, and no quarantine for arrivals from other countries. Policy effectiveness is evaluated based on two factors: the percentage of people allowed to travel and the pandemic situation in the country once the policy is adopted.
1. The percentage of people allowed to enter the country under each policy is denoted Tc.
This number is calculated using the number of people allowed to travel inbound to the country divided by the total number of people willing to enter the country.

34
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 20, 2021. ; https://doi.org/10.1101/2021.04.14.21255465 doi: medRxiv preprint 2. The percentage of people that will travel due to each policy is denoted Te. This number is an adjusted version Tc. If a 14-day quarantine is applied to a country, we assume that only 5% of the normal number of travelers from this country are willing to travel under this policy.
The choice of 5% is based on the data provided by Korea Tourism Organization. (Korea is one of the countries that require a 14-day quarantine for all arrivals.) This adjustment gives us more insights concerning the effect of the 14-day quarantine requirement. After this adjustment, the percentage of expected inbound travelers is obtained by using the number of expected inbound travelers divided by the normal number of inbound travelers.
The effectiveness of policies on the epidemic in the considered country is evaluated based on 7 factors.
1. Percentage of active confirmed imported cases that enter the country due to each policy. This number is calculated using the total number of inbound traveling active confirmed cases that eventually become active confirmed cases, divided by the total number of inbound travelers during the regulation period. We denote this category as IA.
2. Percentage of undetected infected imported cases entering the country due to each policy. This number is obtained using the total number of undetected infected cases traveling inbound divided by the total number of inbound travelers during the regulation period. We denote this category as II.
3. Percentage of undetected infected imported cases when quarantining after entering the country. A policy that does not require quarantine is equivalent to a 0 -day quarantine. This number is obtained by taking the total number of undetected infected inbound travelers after quarantine divided by the total number of inbound travelers during the regulation period. We denote this category as IIQ.

35
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 20, 2021. ; 4. Relative change in total new cases (detected and undetected), denoted as RU. This number is calculated as the difference in the total number of cases at the end of the regulation period and the beginning of the regulation period, divided by the total number of cases at the beginning of the regulation period.
5. Relative change in total new active confirmed cases, denoted as RA. This number is calculated similarly to RU above but instead of using the number of cases, all counts are based on the number of active confirmed cases.
6. Percent change in total new cases, denoted as PU. This number is calculated as the difference in the total number of cases at the end of the regulation period and the beginning of the regulation period, divided by the population of the country. 7. Percent change in total confirmed cases, denoted as PA. This number is calculated as the difference in the total number of confirmed cases at the end of the regulation period and the beginning of the regulation period, divided by the population of the country.
We generate 1000 stochastic realizations conditional on the estimated parameters and the estimated initial conditions at the beginning of the regulation period. For each realization, we calculate the above metrics, and we report the 0.025 and 0.975 percentile values of each based on the 1000 realizations. To give a fair judgment on the effectiveness of travel regulation on the pandemic, we stratify the metrics by dividing countries to three different groups, where Group 1 corresponds to countries with an effective reproduction number R(t) lower than 0.9, group 2 corresponds to countries with R(t) between 0.9 and 1.1, and Group 3 for countries with an R 0 greater than 1.1. Notice that for our model, following Diekmann et al. (2009) (12), we can show that the effective reproduction number R i (t) of a given country i is R i (t) = S (t)+ i α i P i (t)(β i +γ i ) . The overall average across these 200 iterations of the above metrics for each group of countries is calculated and used as the final measurement to compare the effectiveness of different policies.

36
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted April 20, 2021. ; https://doi.org/10.1101/2021.04.14.21255465 doi: medRxiv preprint In Table S5, we see that under P-4 the number of expected inbound travelers, Te, is higher than P-5, the simplified version of our proposed average control policy. However, under P-4, the percent of undetected infected after done with quarantine, IIQ enter countries of Group 1 is about (0.03%, 0.03%) and Group 2 is about (0.03%, 0.04%). These values quite high compare to (0.00%, 0.00%) for both groups as in P-5. This is because the average of increased cases each day in the last 14 days was used to decide which countries belong to a green zone or red zone. However, the number of undetected infectious cases may grow very fast in the green zone countries, and in the absence of quarantine, undetected infectious cases from green zone countries may spread the disease fast in the arrival country.

Simulation 4. Effectiveness of policy coordination
In this simulation study, we study the effectiveness of a global response on the pandemic in terms of the percentage of people allowed to travel and the overall worldwide pandemic situation under a coordinated policy. Simulations are set up similarly to those in Section 2.3 but here we consider a world of 8 countries, where countries 1 and 2 with R(0) greater between 1.1 and 6.47, countries 3 and 4 with R(0) between 1 and 1.1, countries 5 and 6 with R(0) between 0.9 and 1, and countries 7 and 8 with R(0) from 0.47 to 0.9.
We consider 8 different policy coordination scenarios: S-1 All countries are fully open and allow all travel.
S-2 All countries are fully closed and do not allow any airline travel .
S-3 All countries require a 14-day quarantine for all arrivals.
S-4 All countries use the simplified version of the proposed average control policy.

37
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted April 20, 2021. ; https://doi.org/10.1101/2021.04.14.21255465 doi: medRxiv preprint S-5 Countries 1, 3, 5, 7 require a 14-day quarantine for all arrivals, and countries 2, 4, 6, 8 allow no inbound travel.