Both public and non-public data that were used in this study were from the Korea Disease Control and Prevention Agency (KDCA). Publicly available daily number of confirmed cases and vaccine administration were aggregated from daily presses and were used in the mathematical model simulation.37 Two types of information, symptom onset date and diagnosis date, were aggregated from the non-publicly available, individual based data and were used in the maximum likelihood estimation (MLE) process.

In this study, a deterministic mathematical model that includes age, vaccines, antiviral drugs, and influx of the omicron variant was developed. We consider eight age groups and two strains of COVID-19 virus, delta (δ) and omicron (o). These variants have different basic reproductive numbers, transmission rates, and severe rates. Fig. 1 illustrates the flow diagram of the mathematical model. Note that X indicates vaccine- or waning-related status of the host and i refers to age group. There are seven vaccine- or waning-related compartments (X); u (unvaccinated), w (unvaccinated, previously infected, but natural immunity has waned), v1 (two weeks before finishing primary doses), v2 (two weeks after finishing primary doses), wv (waned after primary doses), b (boostered), wb (waned after booster).

An unvaccinated host (u_i) moves to the v1ⁱ compartment after administration of the first dose and has partial vaccine effectiveness. Two weeks after receiving the second dose, the host has full vaccine effectiveness (v2_i). Later, vaccine-induced immunity wanes and so the host moves to the wv_i compartment. The host goes to b_i compartment after receiving a booster shot and later to wb_i, considering the waning of booster shots. In this study, we assume that the immunity against symptomatic infection wanes but the immunity against severe infection does not. This assumption is supported by studies, where a population-wide decline in effectiveness against infection was observed but effectiveness against hospitalization remained high and with no significant change over time.3839 For both variants, an exposed host E_X,i becomes infectious (I_X,i) and spreads the disease until case confirmation, and so the host moves to the Q_X,i compartment. After confirmation, the isolated host either develops mild symptoms M_X,i, including asymptomatic case, or severe/critical symptoms C_X,i. An isolated host with mild symptoms recovers (R_X,i), while those with severe symptoms may either recover or die. We assume that recovered individuals, whether vaccinated or not, develop natural immunity which wanes over time. Since recovered individuals retain protection against severe infection, they move to a different compartment w_i for unvaccinated, wv_i for primary-dosed, or wb_i for boostered, after the natural immunity has waned.39 It was demonstrated that in unvaccinated participants, the infection-acquired immunity waned after about 1 year but remained consistently high in previously vaccinated-participants, even for individuals who were infected 18 months prior.40 Hence, we use a different natural immunity waning rate for those who were unvaccinated (ζ) and vaccinated (ζ_v), with ζ_v < ζ. The parameter Γ_i represents the number of overseas entrant cases from age group i who are not screened and entered the local community. Its value is calculated using data on average daily number of overseas entrant cases across all ages from November 24 to December 31, 2021. The detailed description of the mathematical model, including the governing equations, can be found in the Supplementary Data 1.

Transmission rates among age groups were estimated using MLE. For MLE, we considered two events for a host at one unit time: not being infected and being infected. Individual based data provided by the KDCA were used for the MLE procedure to capture every infection event (also uninfected events) of age groups. Detailed formulation is described in the Supplementary Data 1.

To quantify the impact of NPIs, a time-dependent parameter μ is introduced to indicate the reduction in transmission caused by NPIs. For example, ignoring other factors, if the basic reproductive number is 2 and μ is 0.7, then the effective reproductive number becomes (1 − 0.7) × 2 = 0.6. We estimated the value of μ every week using least squares curve fitting method, by minimizing the difference between the cumulative incidence calculated using the model ($\int {\sum }_X{\sum }_i\alpha (I_{X,i}^{\delta }+I_{X,i}^o))dt$) and the available data. The model simulation time was done from August 1, 2021 to February 2, 2022, because the testing policy has been changed since February 3, 2022.41 Furthermore, we proceeded with parameter bootstrapping to examine the reliability of the estimation and data. Detailed description of the bootstrapping method and results are in the Supplementary Data 1. During the parameter estimation period, antiviral drugs were provided for certain age groups.2021 To apply impact of antiviral drugs, we simply set that severity rate of age over 60 and 40 has reduced since January 14 and February 21, 2022, by 80%, respectively. Note that the 80% severity reduction is between the lowest and highest values that were considered in a recent study.42 For example, delta variant’s severity rate for the unvaccinated individuals aged 60 to 69 reduces from 7.49% to 1.50% since January 14, 2022.

In this research, we performed a multifaceted approach to examine critical factors which affected the COVID-19 epidemic. First, we did a forecast considering different NPIs-related factors and vaccine hesitancy. Second, we examined the time-dependent sensitivity of screening measures to the disease spread since the omicron variant has arrived in Korea. Finally, we derived a manageable daily incidence number from the endemic equilibrium state of the mathematical model.

For the forecast, we extended the simulation time until the end of 2022 and varied the factors related to NPIs and booster shots. We set the range for the NPIs-related reduction factor (µ) from 0.4 to 0.65 in 0.05 increments (six scenarios), and the maximum number of daily booster shots as 300,000 or 100,000, to observe the influence of vaccine hesitancy. Furthermore, to consider the implementation of the antiviral drugs, we set that groups of age over 60 (age over 40) have reduced severity after February 10, 2022 (February 21, 2022).

We examined the impact of screening measures by varying the value of Γ_i by factors of 0.1 to 10, and the date of entry of omicron to the local community from November 24 to December 1, December 8, and December 22. The rest of the parameters are fixed to their values on Supplementary Tables 1 and 2.

As the number of cases rapidly increases, endemicity of COVID-19 becomes an issue. We performed an endemic equilibrium analysis to investigate how COVID-19 cases can be maintained on a manageable level. Ignoring age structure (i), vaccination-related history (X), and strains, and considering endemic equilibrium (assuming that there is natural death and birth in susceptible groups, therefore endemic equilibrium can exist), ordinary differential equations of confirmed (Q) and severe patients (C) are:

where p^* and γ^* are average severe rate and recovery rate, respectively. Combining the two results above, we get:

Because the number of severe patients should be below the severe bed capacity, C_max, the following inequality is formulated:

Considering that αI is the daily incidence and 1/γ^* is the average length of hospital stay, then the threshold value of the inequality, referred to as the manageable daily incidence, is given as follows:

The manageable daily incidence is a function of three input parameters, average length of hospital stay (t_H), average severe rate (p^*), and severe patient capacity (C_max).

The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Institutional Review Board of Konkuk University (7001355-202101-E-130). Informed consent was submitted by all subjects when they were enrolled.

Fig. 2 shows the transmission rate matrix, M₁, represented as a heatmap. The maximum value is 6.14, which is the value among age group 8. Estimated reproductive number from M₁ is 6.16, which is affected by NPIs but not by vaccine because reduced probability of being infected was considered in the MLE process. To exclude the effect of NPIs, the adjusted matrix M₂ was introduced using the basic reproductive number of the variant and the estimated effective reproductive number from the transmission rate matrix, M₁ (Supplementary Data 1). The adjusted matrix M₂ was applied into the mathematical model.

Fig. 3 shows that the daily and cumulative incidences from the model simulation fit the data well (Fig. 3A and B, respectively). Also, the number of administered severe patients captures the trend well (Fig. 3C), even if the parameters related to severe patients were not fitted but aggregated from references.

During social distancing level 4 (August to October 2021), the range of estimated µ was from 0.61 to 0.73, except near the national holiday season (Chuseok, September 20 to September 22, 2021) when µ dropped to 0.4. Since November, as the gradual recovery policy began, µ decreased and ranged from 0.31 to 0.41, and later becomes 0.52 as suspended gradual recovery was announced because the number of severe patients reached more than 1,100. The obtained µ estimates and the corresponding values of the effective reproductive number R_t are illustrated as horizontal lines and dashed curves, respectively, in Fig. 3A. Note that all the estimated values of µ were within the 95% confidence intervals of the parameter bootstrapping results and the details are displayed in Supplementary Fig. 1 and listed in the Supplementary Table 3. The proportion of omicron among new cases (magenta curve) increased from 7% to 71% from December 16, 2021 to January 16, 2022 and reached 97% by the end of the simulation period.

Age groups 1 and 2 (age under 30) showed the maximum and second maximum incidence among age groups in both phases, August 1 to December 31, 2021 and January 1 to February 3, 2022, respectively. Age group 6 (60 to 69) had the third highest incidence number in 2021 but third lowest in 2022. Age groups 8 and 7 (age over 70) had the minimum and second minimum incidence during the simulation period.

Forecast from February 3 to the end of 2022 considering various NPIs-related reduction factor (µ) and maximum number of daily booster shot administration, showing the range of confirmed cases and administered severe patients, are displayed in Fig. 4 If the maximum number of booster per day is 300,000 (100,000), the peak size of incidence and administered severe patients will range from 320,300 (420,900) to 1,409,200 (1,518,900) and 1,210 (1,530) to 5,120 (5,410) according to the value of µ which varies from 0.4 to 0.65, respectively. A secondary wave towards the end of 2022 is observed in each scenario, and the size of the secondary peak (incidence: less than 300,000, severe patient: less than 1,000) is smaller than the first peak. We display the data (red boxes) until March 13, 2022, before testing policy was changed to include positive rapid antigen test done in an accredited facility as a confirmed case.43

Fig. 5 shows the log-scaled simulation results using filled curves with different colors. Red, green, blue, and cyan areas indicate the ranges of daily incidence for the various numbers of overseas entrant cases (0.1Γ_i to 10Γ_i), initiated on November 24, December 1, December 8, and December 22, respectively. As the date is delayed, the ranges of incidence become narrow. The black curve indicates the incidence when Γ_i is set to its baseline value and initiated on November 24. The ratio of the maximum (minimum) to the baseline incidence value when Γ_i is initiated on November 24, December 1, December 8, and December 22 are 8.64 (0.13), 2.61 (0.84), 1.36 (0.96), and 1.04 (0.99), respectively. In particular, if the impact of screening (Γ_i) is varied since November 24, 2021, the range of values of the number of daily cases by February 3, 2022, is (2,730–190,420). On the contrary, if the screening is varied since December 22, 2021, the range of daily cases is (21,800–22,780).

Using Equation 1, if the severe rate, length of hospital stay, and severe patient capacity are 5%, 28 days, and 500, respectively, which might be similar to the early stage of COVID-19 in Korea, then the manageable number of incidence is approximately 360. Fig. 6 illustrates the manageable daily incidence if the severe patient capacity is fixed to 2,800 (Fig. 6A) or when the average length of hospital stay is set to 7 days (Fig. 6B). Fig. 6C shows the actual incidence data and theory-driven manageable incidence using aggregated data of each day, interpolated data of hospital stay, severe rate, and severe patient capacity.4445 Blue square indicates that the manageable daily incidence of Korea in mid-February 2022 is around 600,000 assuming that severe patient capacity is 2,800, average severe rate is 0.16%, and average length of hospital stay is 7 days. Moreover, it is visible that actual daily incidence exceeded the theory-driven manageable incidence in December 2021, when the government declared the suspended gradual recovery.