We analyzed epidemiological data on 30,413,435 reported cases of COVID-19 in the Republic of Korea from February 26, 2021 to March 6, 2023, provided by the Korea Disease Control and Prevention Agency (KDCA) (18). The proportions of the Delta and Omicron variants among all COVID-19 reported cases were obtained from the covariance data (35). The time intervals used for our analysis were categorized into three periods, based on the globally dominant variants, which are Delta and Omicron, observed (18, 19), summarized in Supplementary Table S1. In the Republic of Korea, the emergence of COVID-19 cases with the Delta variant began in April 2021, and by 11 July, 2021, it accounted for more than 50% of the total cases (36, 37). Subsequently, from January 11, 2022, the Omicron variant was the predominant strain, representing more than 50% of the cases, resulting in a rapid increase in the number of COVID-19 cases (38). Therefore, we designated the three periods as “Pre-Delta” (February 26, 2021–July 10, 2021), “Delta” (July 11, 2021–January 10, 2022), and “Omicron” (January 11, 2022–March 6, 2023), and the time intervals were labeled as T_Pre-Delta, T_Delta, and T_Omicron, respectively.

We divided the total population into 17 age groups with interval of 5 years, ranging from 0–4 years to 80 years and older. Weekly vaccination data were extracted for the first, second, and third dose vaccinations administered to different age groups, as provided by the KDCA (35). Details on the Korean population size, segregated by age group and region for the year 2021, were obtained from Statistics Korea (39).

For estimating the probability distributions of the SIs, we first calculated the number of transmission pairs in the data based on infector onset dates. To account for the data reported daily, the discretized probability density function f t θ was defined at time t for the parameter of the distribution θ . For example, in gamma distribution, the parameter θ represents a vector of mean (μ) and standard deviation (SD) (σ) of the probability distribution, such that θ = μ σ . Then, the likelihood function for an SI is defined as

where m is the total number of pairs and d indicates the serial intervals for time period (i.e., d = d 1 , … , d m ). We used 70,414 infector–infectee pairs to estimate the SI distribution, while four commonly used distributions for epidemiological periods: gamma, log-normal, normal, and Weibull (40, 41) were employed to estimate the time period. The performance of each statistical model was compared by calculating the Akaike information criterion (AIC).

where K is the number of parameters used. Among the four commonly used statistical models, the best-fitted distribution was selected based on the minimum AIC values. h τ indicates the probability distributions of the SIs estimated from February 18, 2020 to March 6, 2023 during the total period. We defined the best-fitted distribution as h P r e − Delta and h Delta for the Pre-Delta period and for the Delta period, respectively.

To consider non-positive values in the SI data, the analysis involved two approaches: fitting the distributions to positive values only (truncated), and fitting the distributions to shifted data with 11-day delays added to each observation (shifted) (42). Thus, we assumed that pre-symptomatic transmissions could be accounted by adding delays to each observation. Hence, a more accurate representation of the underlying distribution was captured and meaningful insights from the data were derived.

We assessed R t to quantify the time-dependent variations in the average number of secondary cases generated per case during the course of the outbreak due to intrinsic (decline in susceptible individuals) and extrinsic factors, such as behavioral changes and implementation of public health measures (43–45) In Korea, the reported cases vary throughout the week, with notably lower counts observed on Saturday and Sunday. We applied a moving window using a 21-day window to address this variability. By using the smoothed data on COVID-19 cases, we estimated the evolution of R t for COVID-19 in the Republic of Korea.

Several studies estimated the most recent R t by simulating the progression of incident cases and applying the discretized probability distribution of the generated interval using renewal equations (13, 26, 43, 46). In a model study conducted across 131 countries, the impact of implementing and easing eight different NPIs on R t was examined (47). The reopening of schools; lifting of bans on public events, gatherings of 10 or more people, and stay-at-home orders; and easing of internal movement restrictions were found to increase R t . However, the effects of NPI implementation and easing were not immediate. Using maximum likelihood estimation (MLE) and sequential Bayesian methods, R 0 and R t were estimated (48).

R t was defined as the total number of incident cases I t arising at time t, divided by the discretized probability function h t , which was defined at time t with the lowest value of AIC for truncated distributions, as shown in Table 1.

To compute 95% credible intervals (95% CrIs) of R t , the bootstrapping method was applied to generate 100 samples from the Gamma distributions (26).

Data on the number of COVID-19 cases, proportion of variants, and vaccination coverage are presented in Figure 1. Figures 1A,B present the number of confirmed cases caused by Delta and Omicron variants. The number of cases increased with the emergence of new mutations, with the Delta variant causing the highest number of cases, reaching over 8,000. During Omicron, the number of cases peaked at approximately 600,000 before decreasing. The proportion of the variants are shown in Figure 1C, where the dashed lines indicate the start points of the time periods, T_Delta and T_Omicron during which the proportion of each variant exceeded 50%. The vaccination coverage is shown in Figure 1D. The first and second vaccine doses were administered during Pre-Delta, with the majority receiving the second dose during Delta. The third doses of booster shots were administered 2 months prior to the emergence of the Omicron variant. The vaccination coverage for each period is shown in Figure 1E. During Pre-Delta, the first dose accounted for about 40% of the total population, whereas the second dose accounted for only 10%. During Delta, the first and second doses were administered to 80% of the population. The third dose was administered to 40% of the population during Delta, while 70% were covered during Omicron.

Out of 30,413,435 COVID-19 cases, we reconstructed 70,414 transmission pairs from the known onset dates for the infectors and infected population. The SIs ranged from −11 to 17 days, and were estimated using truncated and shifted distributions for the entire period, as shown in Figure 2. Based on the AIC values, the truncated Weibull distribution provided the best fit for all three periods. The estimated mean SI for the total period was 3.85 days, with an SD of 3.41 days, as shown in Figure 2A. To compute R t , we employed the estimated SI using the truncated Weibull distribution for each period (Pre-Delta, Delta), as presented in Table 1. For Pre-Delta and Delta, the estimated mean SIs were 4.29 and 3.56 days, with SDs of 3.84 and 3.08 days, respectively. The estimated SIs for Pre-Delta and Delta from the truncated and shifted distributions are shown in Supplementary Figure S2. The different values of R t obtained from truncated and shifted distributions using the Gamma, Weibull, Lognormal methods are presented in Supplementary Figure S3.

The values of R t calculated using Eqs. (2) and (3), distinguished as R t and R x t , respectively, were compared, as presented in Figure 3. When R t was calculated based on the total number of cases, significant differences in values for each variant were observed. Calculated R t values using data from February 26, 2021 to January 10, 2022 and from January 11, 2022 to March 6, 2023 are presented in Figures 3A,B, respectively. Supplementary Table S2 summarized the NPI levels implemented in Korea. It is shown in Figures 3A,B, along with the COVID-19 confirmed cases. During the Pre-Delta period with NPI level 2, R t was consistently around 1. Subsequently, a point where R t exceed 1, we could interpret that it is closely related to the emergence of the Delta variant. This relationship is evident from the discrepancies between R t and R x t shown in Figure 3C. Furthermore, during this period, the reduction of the NPI level to 1 led to an increase in R t to 1.2. From July 11, 2021 onwards, R t decreased and remained at approximately 1. It believed that this was a consequence of the NPI level being intensified to 4. After the spread of the Omicron variant in 2022, R t increased to 1.4. On July 11, 2021, R x t for Pre-Delta was lower than R t calculated using the total number of cases, while R x t for Delta was higher. Starting from January 11, 2022, when the Omicron variant accounted for more than 50% of the total cases, R x t for Delta showed a significant decrease compared to R t calculated based on the total number of cases using Eq. (2), as shown in Figure 3D. Although the overall value of R t increased, analysis of R x t , specifically for Delta, revealed a decrease. Thus, despite the overall increase in transmission of the virus in the population, the measures implemented to control the Delta variant were effective in reducing its spread. The NPI intensity was gradually reduced in a phased restoration of daily life, leading to the lifting of social distancing after April 18, 2022. During the Omicron period, R t , exhibited higher volatility than other periods, shown in Figure 3E.

The values of R t and R v calculated for with and without immunity from variant x using Eqs. (3) and (5), distinguished as R x t and R v , x t , respectively, are shown in Figure 4. The monthly mean, SD, and CV of R x t and R v , x t are presented in Supplementary Table S3, while the statistics of their estimated values are summarized in Supplementary Table S4. The difference between R x t and R v , x t was not large because the vaccination coverage was not high during Pre-Delta period, as shown in Figure 4A. R x t considering immunity remained approximately 1, as shown in Figure 4B, and after December 31, 2021, it decreased to a value below 1. However, R v , x t without immunity, always remained greater than 1 during Delta. As shown in Figure 4C, R x t was less than 1 for some data points, whereas R v , x t always remained greater than 1, during Omicron. The boxplot of R x t and R v , x t for each variant is shown in Figure 4D. During Pre-Delta, a small difference existed between R x t and R v , x t , while a significant difference was observed during Delta and Omicron. Thus, a large variability existed for each variant, such that a high variability was observed during Delta, which subsequently decreased during Omicron.

CV was calculated from the monthly mean and SD of R x t and R v , x t and presented in Figure 4E. A small difference was observed between the CV calculated during Pre-Delta and Omicron, but a notable difference was observed during Delta. We defined 2 months with significant differences as “a” and “b,” which were calculated from the shaded area in Figure 4B. During “a,” R v , x t exhibited higher variability compared to R x t , accompanied by a substantial increase in R v , x t . Conversely, during “b,” R x t showed higher variability compared to R v , x t . Although the magnitudes of the changes were similar, smaller values of R x t resulted in larger variability. This period corresponded to the initiation of a third-dose vaccination campaign. Overall, the variability between R x t and R v , x t differed across the variants, with Delta characterized by significant differences and the start of the third-dose vaccination campaign.

We conducted a comparison of R t using Eq. (2) and R x , g t using Eq. (6) with the results illustrated in Figure 5. The location of each region in Korea is illustrated in Supplementary Figure S1. Figures 5A–G presents the results for Seoul Metropolitan Area, Gangwon, Chungcheong, Honam, Gyeongbuk, Gyeongnam, and Jeju. Figure 5H shows a regional box plot comparing R t and R x , g t . The discrepancy between R t and R x , g t is the largest in Jeju, while it is almost negligible in Seoul Metropolitan Area. In all regions, a difference between R t and R x , g t is observed around July 11, 2021, coinciding with the transition from the Pre-Delta to the Delta variant. During this period, Seoul Metropolitan Area, Chungcheong, Honam, and Gyeongnam exhibited similar patterns in R t and R x , g t , whereas Gangwon, Gyeongbuk, and Jeju displayed divergent patterns. Particularly in Jeju, R x , g t exhibits significant volatility during the Pre-Delta period, which is likely attributed to the small number of cases, and Gangwon and Gyeongbuk require additional analysis.

We computed R x , g t with immunity using Eq. (6) and R v , x , g t without immunity using Eq. (7) in seven regions of Korea, as aforementioned for Pre-Delta, Delta, and Omicron is shown in Supplementary Figures S4–S6, respectively. A comparison between the calculated maximum value of R x , g t and R v , x , g t by region g is presented in Figure 6. While comparing the R x , g t results for Pre-Delta, Delta, and Omicron, the severity of the variant viruses was difficult to determine. However, while examining the R v , x , g t results, the overall R v , x , g t values were higher during Delta, particularly in Gyeongnam and Gangwon. During Omicron, the R v , x , g t values were higher than those during Pre-Delta. The mean values of R x , g t and R v , x , g t are shown in Supplementary Figure S7.

The results for three regions, Seoul Metropolitan Area (Region A), Gangwon (Region B), and Gyeongbuk (Region E), were compared, as shown in Figure 7. Variations in the magnitudes of R v , x , g t across different variant periods were observed. Without considering immunity, the mean of R v , x , g t consistently increased over time in all regions. However, a significantly higher value of maximum R v , x , g t compared to an average value of R v , x , g t in each variant period suggested the occurrence of an event that led to a spike in cases in that region. During Pre-Delta, when the increase in the number of cases was relatively smaller compared to that during Delta and Omicron, the difference between the average and maximum values of R v , x , g t was relatively small. Supplementary Table S5 summarizes the maximum, mean, and minimum reproduction numbers for each period by region.

The values of monthly CV of R x , g t and R v , x , g t and vaccination coverage for the three regions are shown in Supplementary Figure S9. Supplementary Table S6 provides a summary of the monthly CV for the two indicators of reproduction number. Supplementary Table S7, on the other hand, outlines Distribution of age population size and vaccination coverage by region. During Pre-Delta, R x , g t and R v , x , g t exhibited similar patterns. However, during the second phase, when the Delta variant spread, a significant variability in R v , x , g t was observed. This period coincided with a rapid increase in the administration of the second and third vaccine doses, indicating an increase in immunity. In the third phase, a consistent pattern in CV across regions was observed, suggesting a similar trend. The number of individuals receiving the third vaccine dose significantly decreased, and owing to an already high number of infected individuals, the difference between R x , g t and R v , x , g t was reduced. In Seoul Metropolitan Area and Gyeongbuk, shown in Supplementary Figure S9, a rapid increase in the number of COVID infected individuals was observed in mid-September 2021. Thus, the variability in CV increased to approximately 0.2.