The impact of multiple non-pharmaceutical interventions for China-bound travel on domestic COVID-19 outbreaks

Objectives Non-pharmaceutical interventions (NPIs) implemented on China-bound travel have successfully mitigated cross-regional transmission of COVID-19 but made the country face ripple effects. Thus, adjusting these interventions to reduce interruptions to individuals’ daily life while minimizing transmission risk was urgent. Methods An improved Susceptible-Infected-Recovered (SIR) model was built to evaluate the Delta variant’s epidemiological characteristics and the impact of NPIs. To explore the risk associated with inbound travelers and the occurrence of domestic traceable outbreaks, we developed an association parameter that combined inbound traveler counts with a time-varying initial value. In addition, multiple time-varying functions were used to model changes in the implementation of NPIs. Related parameters of functions were run by the MCSS method with 1,000 iterations to derive the probability distribution. Initial values, estimated parameters, and corresponding 95% CI were obtained. Reported existing symptomatic, suspected, and asymptomatic case counts were used as the training datasets. Reported cumulative recovered individual data were used to verify the reliability of relevant parameters. Lastly, we used the value of the ratio (Bias2/Variance) to verify the stability of the mathematical model, and the effects of the NPIs on the infected cases to analyze the sensitivity of input parameters. Results The quantitative findings indicated that this improved model was highly compatible with publicly reported data collected from July 21 to August 30, 2021. The number of inbound travelers was associated with the occurrence of domestic outbreaks. A proportional relationship between the Delta variant incubation period and PCR test validity period was found. The model also predicted that restoration of pre-pandemic travel schedules while adhering to NPIs requirements would cause shortages in health resources. The maximum demand for hospital beds would reach 25,000/day, the volume of PCR tests would be 8,000/day, and the number of isolation rooms would reach 800,000/day within 30 days. Conclusion With the pandemic approaching the end, reexamining it carefully helps better address future outbreaks. This predictive model has provided scientific evidence for NPIs’ effectiveness and quantifiable evidence of health resource allocation. It could guide the design of future epidemic prevention and control policies, and provide strategic recommendations on scarce health resource allocation.


Introduction
COVID-19 has created a global challenge that demands researchers, policymakers, and governments address multiple dimensions which go far beyond the implications of human health and well-being (1)(2)(3)(4). Scientific evidence has indicated that non-pharmaceutical interventions (NPIs) are effective measures to contain a pandemic and ease pressures on healthcare systems (5)(6)(7). NPIs are actions, apart from getting vaccinated and taking medicine, that people can take to help slow the spread of illnesses, also known as mitigation strategies (8-10). It includes travel restrictions, contact tracing, PCR tests, measures in social distancing, personal protection, and quarantines (6,11,12). The implementation of such interventions while maintaining social stability is a challenge to all countries. As a country consisting of more than 1.4 billion or 18% of the world's population, China's high population density, high volume, speed, and non-locality of human mobility would provide perfect conditions for the virus to spread (13,14). When highly transmissible Delta and Omicron variants resulted in massive surges in COVID-19 cases from December 2021 (15,16), China saw the largest spike for the past 2 years, despite determinedly pursuing one of the world's strictest virus elimination policies. When a local COVID-19 case occurred, mandatory interventions would be taken to cut off the transmission chain and terminate the outbreak in time to achieve maximum effectiveness with minimum cost. After years of exploration, such strategies' implementation received remarkable results in containing regional cases (17,18). However, it required extensive community involvement, government funding guarantees, application of new technology, motivation, and constraint mechanisms. Such a strategy created indefinable impacts on regional social development (19,20). Thus, knowing how to maximize the advantages of strategy in outbreak control while avoiding damaging the development of the country was critically important. Due to the combined use of NPIs in the strategy, we decided to quantify the impact of different NPIs. Extensive research was conducted by using a time-varying modelinginformed approach and focusing on the following three interventions in this paper: inbound flight restrictions, PCR tests, and centralized quarantine measures.
Severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) viral spread patterns were shaped by the high volume of cross-country mobility networks (21). In response to the pandemic, China reduced inbound flight schedules from 10,000 per week in 2019 to 500 per week in recent years (22), and international arrivals were reduced from approximately 162.5 million in 2019 to 30.4 million in 2020 (23). In July 2021, the aviation authority updated requirements-passengers were required to complete a PCR test within 5 days of embarkation and provide negative test results before boarding, as the government tried to further reduce the risk of imported cases (24). However, from July 1 to July 31, 2021, 1,213 confirmed COVID-19 cases were reported across the country, compared with 1,893 cases in August and 1,264 cases in September (25). Although travel restrictions and PCR tests were proven as useful practices (26), the theoretical basis of those strategies and how to strategically align them with a country's development was not studied.
There was a high level of agreement that the adoption of travel measures led to important changes in the dynamics of the early phases of the COVID-19 pandemic (27). Flight restrictions may have led to additional reductions in the number of exported and imported cases on the international scale, but such limitations (up to 90% of traffic) had only a modest effect unless combined with a 50% or higher reduction of transmission in the community (28). With the occurrence of domestic COVID-19 outbreaks, the association between international travel and the implementation of NPIs has not been identified. NPIs such as centralized hospitalization for mild and moderate patients could reduce disease transmissions and enhance protection for healthy and unhealthy individuals (29, 30). Nevertheless, the efficacy of mandatory isolation for international travelers at a designated place in a given period was not discussed. Research on Hong Kong-bound air passengers indicated that home quarantine was less effective than a centralized quarantine strategy initially but showed similar efficacy in the later phase (31). However, the effectiveness of self-isolation, transmission rate within the family cluster, related disease burden, and consumption of public health resources were not mentioned. According to a study published by United States Centers for Disease Control and Prevention, the transmission of SARS-CoV-2 among household members was common, and secondary infection rates were higher and occurred rapidly, with approximately 75% of infections identified within 5 days of the index patient's illness onset (32). Substantial transmission occurred whether the index patient was an adult or a child, leaving no one healthy enough to help other family members.
Mathematical models and time series analyses have been widely used to study the pandemic and predict the trend. Researchers used a time-dependent SIR model to track the transmission and recovery rate at time t and presented less than 3% of one-day prediction errors (33). But the effects of NPIs were not discussed in the research. Another time-varying SIRD model was also developed to capture possible changes in the epidemic behavior, due for example to containment measures enforced by authorities or modifications of the epidemic characteristics and to the effect of advanced antiviral treatments in Italy (34). However, the research team did not take the interaction effects between containment measures and international travel bans into consideration. To infer more accurate parameter estimates and reduce uncertainties, scholars used real datasets of COVID-19 cases via an SEIR model with time-varying transmission and reporting rates to perform 1-week ahead predictions and generated more realistic interpretations (35). Despite that, this model was designed to predict the number of under-reported active cases not for NPIs evaluation, strategic planning, and resource allocation.
Thus, we would develop epidemiological models to simulate the domestic spread of SARS-CoV-2 sparked by passengers who had followed NPIs, such as inbound travel restrictions, quarantine measures, and PCR tests. However, the traditional epidemiological models fail to show the real-time implications of NPIs implementations, delayed symptoms, and test results. To present the time-varying effects, we developed a homogenous hybrid dynamic Susceptible-Infectious-Recovered (SIR) model to quantify such implications. The model can capture multiple data resources rather than a single dataset and generate a more robust estimation of the underlying dynamics of transmission from noisy data. Furthermore, it clearly described the synergistic effects of multiple interventions, such as face masks and social distancing. By combining an improved SIR model with four datasets collected from July 21 to August 30, 2021, we explored the sustained humanto-human transmission relationship between the inbound travelers and the domestic outbreaks under effective NPIs. Based on the Frontiers in Public Health 03 frontiersin.org simulation results, we formed a comprehensive model to quantify the impact of each NPIs and predicted the trend of future outbreaks based on the implementation of these NPIs. The goal was to ① explore the relationship between the imported cases and the development of the domestic epidemic, ② discuss how to adjust existing prevention and control strategies based on our findings, and ③ prepare sufficient health resources in advance while preventing health systems become overwhelmed. Moving forward, we would like to explore the balance point in epidemic prevention and international travel restrictions that could minimize the disruptions to social development.

Model assumptions for consideration
The total population was 1,411,478,724 except for Hong Kong, Macau, Taiwan, and about 300,000 who are naturally immune (36).
• Assuming that the population is closed, meaning that there are no births and deaths. Population migration status change is considered during the study period, but they are dynamically stable, then • Assuming the population is homogeneously distributed and individuals mix uniformly. • Assuming that the infectiousness of symptomatic and asymptomatic individuals is the same in a real-world scenario (37). • Assuming that the recovered patients are negligible during the early stage of the pandemic and their presence will likely not affect the disease transmission (38)(39)(40). • Assuming that symptomatic and asymptomatic cases will be moved into convalescence after rehabilitation due to COVID-19 immunity after infection. • Assuming the effect of vaccines, average delays between symptom onset and test results are constant. • Assuming all inboard and abroad travelers have performed the PCR tests, centralized quarantine, and completed treatments at designated hospitals.

A homogenous hybrid network-based model of SARS-CoV-2 transmission
The SIR model was used to model the spread of infectious diseases among a fixed population. This classic compartment model divided the population into susceptible (S), infected (I), and recovered (R) individuals and track the transitions of individuals among these states. It is a deterministic model of a homogeneous population with wellmixed interactions. Since China is continually updating its prevention and control measures, we extend the SIR modeling framework to nine classes: susceptible (S), carried (C ), asymptomatic infected (I a ), symptomatic infected (I s ), recovered (R), quarantined (Q), dead (D), immigrated (L i ), and emigrated (L e ) to study the SARS-CoV-2 transmission on dynamic networks. Especially asymptomatic infected (I a ) are individuals who show no symptoms but PCR test positive, and virus-carrier compartment (C) represents individuals who show no symptoms and PCR test negative but infectivity. Furthermore, quarantined (Q), immigrated (L i ), and emigrated (L e ) compartments are designed to analyze the effectiveness of NPIs, such as the inbound flight restriction, PCR test, and centralized quarantine.
In the system of improved SIR model (Figure 1), α 0 t ( ) represents the percentage of inbound passengers. They are required to stay in a designated place for X days upon arrival and receive closed-loop care. A portion λ of Q will move to S, a portion δ i of Q will move to I s , and a portion δ q of Q will move to I a . Once they entered into the susceptible group S, there is a risk ratio β of S to move some of them into C and diagnosed as I s or I a by the transfer rate of ε and e q respectively. In addition, a portion p of S determined by close contact and sub-close contact tracing will move to quarantined Q. In the meantime, a portion of q i and q r represent the I a will move to I s and R. With the above, since population fraction in compartments S C Q I I L L R D a s i e , , , , , , , , varies with time t (in days), we assume S(t) + C(t) + Q(t) + I a (t) + I s (t) + L i (t) + L e (t) + R(t) + D(t) = N(t)==N, the following kinetic equation is obtained. Initial values, conditions, and descriptions are presented in Table 1. to inbound flight restrictions. The parameters λ δ , i , and δ q are dependent on centralized quarantine measures. Those four dependent variables are mainly changed by the independent variables, i.e., x 0 , x 1 , x 2 , and γ . x 0 represents the validity period of the PCR test, x 1 is the number of international flights, x 2 is the strength of the centralized quarantine measure, γ is the weight parameter related to the incubation period of SARS-CoV-2.
α 0 t ( ) as the main explanatory variable, signifies the proportion of the population migrating to China from other countries. We have modeled the population entry rate via the contribution of the validity period of the PCR test and the restrictions on international flights according to the characteristic of immigration by actual data tracing, shown as the formula (1): To simulate the number of international flights, we set parameters γ is the weight parameter only affected by the effective duration of the PCR test x 0 . After we draw a curve of best fit, we find the effects of PCR test validity period setting are in line with the logarithmic function. This means the virus incubation period could influence the test validity period (46). When the test validity period is shorter than the incubation period, the effect of the validity period of the PCR test conforms to the significant variation part of the logarithmic function, so set γ = 1. If the test validity period is longer than the incubation period, the effect of the validity period of the PCR test conforms to the gently part of the logarithmic function, so set γ = 100000 0 x .
The parameter λ is the release ratio at the end of the quarantine, which follows an exponential distribution with parameters c1 and c (47): The parameter δ i is the probability of the quarantine measure to the symptomatic infectious individuals, and the parameter δ q is the probability of the quarantine measure to the asymptomatic infectious individuals (48). Additionally, 0 ∆ , ρ 0 , η 0 , 1 ∆ , ρ 1 , and η 1 are all the Improved SIR model on SARS-CoV-2 transmission. Dashed lines are influence parameters refer to a real-world scenario where the untraceable infections were reported. For example, untraceable infections that caused by contaminated cold-chain products θ 2 and infections rate θ 1 (t) that trigger local outbreaks. Solid lines are transition probability of compartments. Parameters α 0 (t), λ, δ i and δ q are related to the NPIs implemented for Chinabound travelers and α 1 , α 2 , and α 3 are the outbound parameters; r i ,d i are the recovery rate and q r is the death rate; p, β, ε, e q , q i are the transition probability. Furthermore, the arrows represent the direction of transition/influence between compartments. With above, the initial values and detailed values are presented in Tables 1, 2.
Frontiers in Public Health 05 frontiersin.org fitting parameters, we also derive the 95% confidence interval (CI), which is shown in Table 2: In the context of infectious disease control, curtailing interactions between infected and susceptible populations, reducing the infectiousness of symptomatic patients, reducing the susceptibility of susceptible individuals, and scaling up such intervention coverage to accommodate rapid increases in the number of suspected cases are well-known strategies for minimizing pandemic spread (49). China has adopted measures conforming to China's conditions based on the strategic theory, i.e., local management. When an outbreak occurs, a local management strategy will be implemented in that particular city. To model the local management policy concretely, dynamic parameter θ 1 t ( ) varying with time is introduced to the improved model. The parameter is determined by the number of cities with infected cases and the population of each city. To enhance the generation ability of the model, we set the city size equal to 4,000,000 residents (50). Since the centralized quarantine strategy of inbound flights is managed in a closed loop, and researches show the majority of domestic outbreaks were caused by contaminated imported cold-chain food (51,52) which was less traceable, we set θ 2 as the probability of infection caused by cold-chain propagation.

Parameters setting and parameters estimation
According to VariFlight, there were an average of 16,707 inbound immigrants and 12,310 outbound emigrates per day. Deidentified aggregated data collected from July 21, 2021, to August 30, 2021, was used to fit the inbound parameter α 0 t ( ), the outbound parameters α α 1 2 , , and α 3 (54). To study the impact of the scenario with the normal inbound flights on the domestic outbreaks and economic , , , , , , β ε ). Parameters p and β were defined via reference (40). However, as the Delta variant continued to mutate, the early transmission rate β was lower than the current variation (Table 3). In addition, the average incubation period of the Delta virus was about 4.4 days and about two-thirds of those infectious cases were symptomatic (55), corresponding to ε + = e q 1, as a result, we set the transfer rate ε as 1/4.4*2/3. Furthermore, according to the study report (53), the average recovery period was between 14 and 28 days, thus we set r i and q r equal to (1/28, 1/14). Lastly, historical data has shown zero deaths during the selected period, so d i was set as zero.
To investigate how NPIs implementation impacts the outbreak duration or the turning point, the logic parameters ( 0 ∆ , ρ 0 , η 0 , 1 ∆ , ρ 1 , η 1 , e, l, e 1 , l 1 , q i , c, c 1 ) associated with fitting functions were estimated by Monte Carlo Stochastic Simulation (MCSS) approach. To get the probability distribution for variables related to population behaviors, a large number of simulation repetitions were needed to stabilize the frequency distributions. Parameters were randomly generated within a range equal to their best fit to the observed data or literature via efficient Python software, then we ran the MCSS method with 1,000 iterations to derive the probability distribution of those variables. Finally, we obtained the initial values and estimated parameters of the model, and listed parameters, initial values, as well as corresponding 95% CI in Table 2.
We further compared the prediction results with three training datasets to determine their final parameters solution aiming to minimize RMSE. To verify the validation of the SIR model and estimated parameters, we compared the model with the testing dataset. Predictive results indicated that the estimated values were in very good agreement with real reported data and that the estimated parameter values can be used to predict the future development trend of COVID-19 in mainland China.  (Figures 2A-C). Reported cumulative recovered individuals were used as a testing dataset to generate ( Figure 2D). To verify the model's reliability, root mean square error (RMSE) was adopted to cross-validate the predicted results and the realworld results. Since a smaller RMSE result refers to a better fitting result, by putting the weight vector quantity (1,0.1,1) to training datasets to reach a goal of minimum RMSE, we obtained the optimal parameters solution. Finally, for reliability verification, the optimal parameters were assigned to the target model to obtain the predicted results and compare it with the trend of cumulative recovered individuals.

Model verification of reliability, stability, and sensitivity
To verify the model's stability while generating the best model fitting result, we identified the equilibrium point between variance and bias, and set the value of ratio (Bias 2 /Variance) in the interval [0.5,1.3], based on bias-variance dilemma theory ( Table 3).
The sensitivity of NPIs on infected cases was tested in this section. Since the amount of three intervention combinations was 2,197, it was unrealistic to observe the effect of simultaneous changes on infected cases. In this paper, the changing influence of each NPIs on infected cases was observed while the other two NPIs maintain normal. , We completed sensitivity analysis on each travel-related intervention with input parameters, for example, the parameters e，e 1 of validity period setting of PCR test, the parameters l, l 1 of the control of inbound flights, and the parameters c, 1 c , 0 ∆ , ρ 0 , η 0 , 1 ∆ , ρ 1 , and η 1 of the strength of centralized quarantine. To quantify the parameter sensitivity of each intervention, we set the number of infected cases caused by current travel interventions as N*. Then the intensity of each intervention was set to vary around its mean 20%, to derive N i . We calculated the relative error of input parameters of each intervention according to the formula [abs (N * −N i )/N * ], as listed here (Table 4). We could observe that the input parameters sensitivity of the validity period setting of the PCR test was the highest, and the sensitivity of the control of inbound flights was the lowest. Thus, the results showed the input parameters of the PCR test were more stable than the other two types of input parameters.

Demand for health resources
The prevalence of COVID-19 worldwide will increase the risk of local transmission. Our model has described a scenario on how to allocate health resources in preparation for possible outbreaks when international flights have been reduced from 1,366 to 79. Figure 3 showed the predictive demand for hospital beds, PCR test volume, and centralized isolation rooms.
Firstly, Figures 3A,D,G showed how the number of international flights impacts hospital bed demand. We stipulated the number of beds in use was configured to be equal to the number of infected individuals to visualize hospital bed occupancy based on the China CDC's requirements (56).  Based on the analysis, we found that when the number of international flights was doubled (x 1 = 160), the number of hospital beds in use would increase by 83%, the PCR test volume would increase by 44%, and the number of isolation rooms in need was doubled. The results showed that the growth in the number of international flights had the greatest impact on isolation room demand. When the number of international flights increased from 79 to 1,366, the demand for hospital beds raised to 25,000/day, the PCR test volume was up to about 8,000/ day, and 800,000/day isolation rooms within 30 days were in need in preventing the spread of the epidemic. Our simulation results indicated that, under those epidemic prevention and control strategies, China was not ready to fully resume pre-pandemic international travels due to excessive demand for health resources. Additionally, the prevalence of COVID-19 in the surrounding countries would increase the probability of a sizable domestic outbreak. To prevent excess demand for health resources, the implementation of an aggressive disease prevention and control strategy was recommended.
As the virus continues to evolve, China is likely to readjust its preventive policies, we will discuss how these future modifications would impact the spread of disease and demand for health resources in the follow-up study.
Frontiers in Public Health 09 frontiersin.org

Effectiveness of NPIs and risk warning of domestic outbreaks
Our retrospective model has indicated that NPIs on travel requirements have successfully contained the spread of the virus. In this section, we will discuss the control of the number of inbound flights, the validity period setting of PCR tests, and the strength of centralized quarantine. By observing and analyzing changes in the number of infected cases and level of intervention implementation, the result will show the effectiveness of NPIs and risk warning of future domestic outbreaks. Figure 4 shows that with the increase of validity period setting of PCR test, the number of symptomatic and asymptomatic infected cases will continue to grow until converging to a stable state presenting no effect of the intervention. Figure 4 shows there are no major fluctuations in the number of infected cases when the validity period of the PCR test is in a 4-day window. However, when we adjust the validity period to 7 days and more, the number of infected cases will be in a stable state. Our results show it is necessary and urgent to set a PCR test time requirement before travelers' arrival. Secondly, the simulation shows that the validity period of the PCR test is closely related to the incubation period of the Delta variant, thus, the test validity period is suggested to be set within 4 days. To maximize impacts, the validity period should not exceed 7 days.  Figure 6 shows how centralized quarantine influences the number of infected cases. With the extension of the quarantine period, the number of infected cases will continue to grow. It can be observed that the impact of the intervention is still remarkable within threshold 35 on preventing the spread of the epidemic and the number of infected cases is converging to a stable situation when exceeding threshold 35. The model also indicates that 17 days of centralized quarantine would  The validity period setting of PCR test vs. infected cases count. The number of inbound flights vs. infected cases count.

The strength of centralized quarantine
effectively prevent disease spread. The quarantine benefit will diminish after 17 days benchmark and reach a stable state after 35 days.  Figure 7B. To quantify the difference in the infected cases between 2(79)-and 3(79)-day in the restricted scenario, we used RMSE to measure the gap, deriving about 53.431. For the small difference between 2(1366) and Studies performed in the United Kingdom and the United States indicated that the effectiveness of any single NPIs was likely to be limited, combining multiple interventions was worthy of further study (57). Scholars also indicated that the effectiveness of travel bans in reducing the spread of infectious diseases, and the relative effectiveness of NPIs for controlling the pandemic has gone largely unstudied (58,59). Therefore, our proposed model played a significant role in estimating the combined effects of NPIs implementation and predicting the demand for isolation rooms, PCR test volume, and hospital beds. The results could provide scientific guidance for nationwide strategic planning and policy implementation and also bridge the theoretical gaps between international travel controls and related effectiveness of the NPIs.

Comprehensive review of all interventions
On one hand, we analyzed how inbound flights would impact the distribution of health resources in response to a possible local outbreak. The model quantified the impact of local virus carriers The strength of centralized quarantine vs. infected cases count. Frontiers in Public Health 12 frontiersin.org that supported PCR testing arrangements for community screening. The number of infected cases and quarantined population could support the allocation of hospital beds and the configuration of isolation rooms. Thus, we recommended that the government should restore inbound flight numbers appropriately with sufficient medical supply in response to the increase in daily infected cases.
On the other hand, our model and related results provided scientific evidence that supported the design and implementation of existing interventions. The results indicated that comprehensive interventions of a two-day PCR test, 79 inbound flights per day, and 17 days of centralized quarantine were effective in stabilizing domestic disease transmission. In addition, the modeling effort also provided theoretical advice for future adjustment. When the epidemic prevention and control goal is to treat and monitor the health status of all infected individuals, limiting the inbound flight number to a small scale is recommended. When the priority is to treat severe and critical cases in hospitals and monitor the health status of individuals who have mild or no symptoms at home, resuming a regular inbound flight schedule is recommended.

Application of the improved SIR model from a micro perspective
The risk estimation of COVID-19 importation can be applied to identify the effectiveness of travel-related control measures (60). However, the connection between imported cases and local outbreaks was not studied. In our model, parameters θ 1 and θ 2 were key factors to understand and mitigate domestic outbreak risks, and also represented mathematical logic interaction associated with the domestic outbreak and global pandemic status. Going further, the current improved SIR model provided more heuristic thinking for constructing new models for domestic outbreaks affected by various factors.

Application of the improved SIR model at other variants of SARS-CoV-2, such as omicron
The new variant SARS-CoV-2 Omicron demonstrated partial vaccine escape and high transmissibility, with early reports indicating lower severity of infection (47) and reduced risk of hospitalization (61) than pre-existed variants. We would like to extend the deltafocused simulation model and related control strategy parameters to Omicron and discuss the applicability and sustainability of the continued implementation of such strategies in combating the new variants in our future research.

Limitations
Our study has several limitations. First, our model did not consider individuals' preventative behaviors. Secondly, we only considered the nationwide prevention strategies and did not dive into detailed strategies enacted at the province and city levels. To minimize such impacts, we adopted reasonable assumptions about epidemiological parameters and aspects of human behaviors that contributed to disease transmission. Although the results showed that our conclusions were remarkably robust, this model was highly sensitive to the quality of input parameters. Thus, we cautiously selected parameters and values based on literature research results and research data. In the proposed homogeneous hybrid model, the population and individuals were distributed and mixed homogeneously and uniformly. Disturbances, such as economic status, political environments, living environments, cultural influences, etc. remained the same. In the meantime, the transmission coefficient, and average delay between symptom onset and test results were constant, and the effect analysis of vaccines, the reporting delays, and testing delays were not captured, which would lead to the requiring hospitalization or developing severe COVID-19 stochastic by nature. Since the purposed model was set to make conservative predictions, when a new variant presents different severity, infectiousness, and immune escape features, we need to convert the purposed model with updated parameters and generate up-to-date predictions. Finally, the model neglected the stochastic effects at low-case numbers. When imported infections were reported, especially when testing was required, having or not a population-scale outbreak was a matter of probabilities; differential equation models cannot capture this accurately. Furthermore, for a disease like COVID-19 with such an overdispersed individual variation of infectiousness (62), outbreaks were likely to die out if very few cases were introduced (63).

Conclusion
Our finding indicated that restriction on inbound flight numbers played a key role in preventing and controlling the epidemic, but the combined use of other NPIs would maximize the effect in preventing additional transmission. Centralized quarantine days should be set in between 17 to 35 days for the Delta variant. The validity period of the PCR test was related to the disease incubation period, and the valid time should be less than 7 days. In addition, when the disease incubated, the PCR test period did not have a significant impact on epidemic control. More importantly, the model estimated that if recovering the pre-pandemic inbound travel strategy in 2019, the number of hospital beds would reach 25,000 per day, the volume of PCR tests would be 8,000 per day, and the isolation capacity would be nearly 800,000 per day within 30 days to maintain the same achievement of preventing outbreaks. All in all, our improved model, which can robustly generate scenarios, will help understand the tradeoffs between different strategies, and further guide the health resources preparation and allocation.

Data availability statement
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding authors.