The total population was 1,411,478,724 except for Hong Kong, Macau, Taiwan, and about 300,000 who are naturally immune (36).

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 well-mixed 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 a , I s , L i , L e , R , D 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.

Multipronged interventions have considerable positive effects on minimizing the spread of outbreaks, decreasing the reproduction number, and reducing total infections. To further clarify the mechanism of interventions and additive effect on epidemic, parameters α 0 ( t ) , λ , δ i , and δ q related to the NPIs implemented for China travelers are constructed in the improved SIR model. Especially α 0 ( t ) is a comprehensive parameter determined by the parameter τ ( t ) related to interventions PCR test and the parameter ϑ ( t ) related 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 τ ( t ) and ϑ ( t ) varying with time t . τ ( t ) represents the contribution of the effective duration of the PCR test and ϑ ( t ) represents the contribution of the number of inbound flights on population entry rate at time t . Then we find τ ( t ) is linear to the weight parameter e 1 (43), and the weight portion is e times the reciprocal relationship with the number of inbound flights x 1 and is logarithmic with the effective duration of nucleic acid testing x 0 (44). ϑ ( t ) is linear to the weight parameter l 1 , and the weight portion is l times reciprocal relationship with the effective duration of the PCR test by fitting to the data (45):

γ 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 x 0 .

The parameter λ is the release ratio at the end of the quarantine, which follows an exponential distribution with parameters c 1 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 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.

July 21, 2021, was set as the starting date of this study. The initial value of S ( 0 ) was collected from the Seventh National Population Census. The initial values of existing symptomatic cases I s ( 0 ) , existing asymptomatic cases I a ( 0 ) , existing suspected cases Q ( 0 ) , cumulative recovered individuals R ( 0 ) , and cumulative deaths D ( 0 ) were captured from July 21, 2021, based on the National Health Commission of China reports. L e ( 0 ) and L i ( 0 ) were collected from VariFlight since July 21, 2021. Since the incubation period is around 4 days, the existing virus-carried cases C ( 0 ) were set to equal to the new domestic case count after (0 + 4) days, i.e., July 25, 2021. Based on VariFlight data and travel requirements, all international flights’ capacity were set to equal to 50% of the original capacity. For better versatility, the average population for medium-sized cities in China was set as 4,000,000 (53).

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 development, we collected the historical data from July 1 to 14, 2019, to simulate the future flow trend of inbound travelers, observe the development trend of COVID-19 and summarize recommendations.

Theoretically, without considering the epidemiological characteristics of SARS-CoV-2, this generic improved SIR model could provide estimation with the above parameters ( p , β , ε , r i , d i , e q , q r ). 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–D were simulated based on existing symptomatic, suspected, and asymptomatic cases and cumulative recovered individual datasets, reported by the National Health Commission of China from July 21, 2021, to August 30, 2021. The reported existing cases were set as training datasets to generate (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 real-world 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.

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²/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. Especially, Based on July 21, 2021, to August 30, 2021, NPIs requirements ( x 0 = 2 , x 1 = 79 , x 2 = 17 ), 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 , c 1 , Δ 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.

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).

Figures 3B,E,H simulated the demand for PCR tests, which was achieved by the product of the obtained number of virus carriers and their highest transmission coefficient. Lastly, Figures 3C,F,I indicated the isolation rooms demand varies with the number of inbound flights. Based on the current quarantine requirements, one person per private hotel room, we could configure the isolation rooms in unit proportion with the isolated population.

Based on the analysis, we found that when the number of international flights was doubled (x₁ = 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.

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.

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.

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 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.

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.

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 delta-focused 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.

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 over-dispersed individual variation of infectiousness (62), outbreaks were likely to die out if very few cases were introduced (63).

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.