Modeling-based design of adaptive control strategy for the effective preparation of ‘Disease X’

Abstract

This study aims at exploring a general and adaptive control strategy to confront the rapid evolution of an emerging infectious disease (‘Disease X’), drawing lessons from the management of COVID-19 in China. We employ a dynamic model incorporating age structures and vaccination statuses, which is calibrated using epidemic data. We therefore estimate the cumulative infection rate (CIR) during the first epidemic wave of Omicron variant after China relaxed its zero-COVID policy to be 82.9% (95% CI: 82.3%, 83.5%), with a case fatality rate (CFR) of 0.25% (95% CI: 0.248%, 0.253%). We further show that if the zero-COVID policy had been eased in January 2022, the CIR and CFR would have decreased to 81.64% and 0.205%, respectively, due to a higher level of immunity from vaccination. However, if we ease the zero-COVID policy during the circulation of Delta variant from June 2021, the CIR would decrease to 74.06% while the CFR would significantly increase to 1.065%. Therefore, in the face of a ‘Disease X’, the adaptive strategies should be guided by multiple factors, the ‘zero-COVID-like’ policy could be a feasible and effective way for the control of a variant with relative low transmissibility. However, we should ease the strategy as the virus matures into a new variant with much higher transmissibility, particularly when the population is at a high level of immunity.

Main

Throughout human history, there have been lots of alarming epidemics, such as the Plague of Justinian and the Black Death, all of which resulted in enormous losses and fatalities [1, 2]. While human public health standards advance with the development of society, the terrifying specter of devastating pandemics has never truly departed from humanity. The severe acute respiratory syndrome (SARS) that emerged in 2002 was caused by the SARS coronavirus, resulted in approximately 8,000 cases and 774 deaths [3, 4]. The H1N1 influenza pandemic in 2009, also known as swine flu, was caused by the H1N1 influenza virus, also caused approximately 200,000 deaths worldwide [5,6,7]. The Middle East Respiratory Syndrome (MERS), which emerged in 2012, is caused by the Middle East Respiratory Syndrome Coronavirus (MERS-CoV), first discovered in Saudi Arabia. It has continued to cause sporadic cases and outbreaks, and remains widely studied [8, 9]. The 2014 West African Ebola epidemic, while not a traditional pandemic, resulted in more than 28,000 cases and 11,000 deaths. It was classified by the World Health Organization as one of the most serious viruses in terms of its high fatality rate [10,11,12]. Lastly, there is COVID-19, which was discovered at the end of 2019 and has rapidly spread worldwide. Its impact on the world remains profound and enduring [13,14,15,16,17].

In recent times, an escalating awareness has emerged that a more perilous pandemic than COVID-19 may swiftly arise in the near future. The World Health Organization (WHO) refers to this potential threat as ‘Disease X’ [18]. ‘Disease X’ does not represent a specific existing ailment but symbolizes a novel, unknown disease to humanity. When it initially emerges, little is known about its characteristics. Whether it is lethal, highly contagious, or poses a threat to human lifestyles remains unclear. However, an increasing number of people are reaching a consensus: COVID-19 may be the first instance of ‘Disease X’, and the next ONE might already be circulating as a respiratory virus in animals [19,20,21,22]. It is both urgent and critical to prepare for the management of emerging infectious diseases. Developing a general and adaptive strategy based on our experience with existing epidemics would be a feasible and effective approach.

During the 4 years of the COVID-19 pandemic, we have developed a global toolkit of non-pharmaceutical interventions (NPIs), including wearing masks, social distancing, personal hygiene, testing, contact tracing, and isolating infected individuals [23]. Equipped with the NPI toolkit, the key question is how to design a proper or an improved strategy that controls the disease’s spread in the most cost-effective manner [24,25,26,27,28,29,30,31,32,33,34,35,36,37]. To be concluded, during the COVID-19 pandemic, there are two distinguished strategies adapted by different countries. Many countries implemented strict control measures during the early phase of the pandemic, and began to relax the strict policies to restore socio-economic activities [38,39,40]. And the focus of the strategy is to flatten the epidemic curve instead of clearing the case using NPIs. In contrast, China has implemented strict control policy (i.e., zero-COVID policy) for a long period till the end of 2022, where the original/Alpha variant and the Delta variant were circulating. The target of zero-COVID policy is to clear the case of each local outbreak in a short period, which has been demonstrated to be a feasible way [14, 27, 29, 33, 35, 41,42,43,44,45,46]. Some concise and useful mathematical models have been established to discuss China’s strict social restriction measures. For example, [43] estimated the 95% credible interval of epidemic parameters in China using data from January 21 to August 7, 2020. They also demonstrated the success and effectiveness of China’s strict measures. However, after the emergence of the Omicron variant, localized outbreaks in places like Shanghai, China, have indicated that even stringent intervention measures may reach their effectiveness limit [47,48,49,50,51,52,53,54,55,56,57,58]. Consequently, the Chinese government announced major changes to the country’s policies on COVID-19 (i.e., the ease of zero-COVID policy). However, there is still a lack of systematic quantitative evaluation of the COVID-19 control strategies used in China, including policy shifts and their optimization. Such evaluation is essential for providing robust qualitative support for the design of general and adaptive control strategies for an emerging infectious disease. Incorporating insights into this aspect falls within the scope of this study.

In conclusion, China, one of the world's most populous countries, followed a zero-COVID policy for three years until the end of 2022. With a national vaccination rate of 92.48% among its 1.4 billion people [59], China then shifted from strict control measures to more normalized interventions. Therefore, China's approach to controlling COVID-19 can serve as a foundational case study for developing a general and adaptive strategy for responding to rapidly evolving diseases like Disease X. The main purpose of this study is to generalize an improved and adaptive strategy from the control of the COVID-19 pandemic in China, by using the compartment model to conduct the retrospective analysis during the circulation of different variants of SARS-CoV-2. The key aspects for searching improved control strategies include if we can try to relax the strict control interventions at the much earlier phase (such as since the emerges of the Delta variant)? And which time would be the more proper for the ease of the zero-COVID policy when we face the Omicron variant with such high transmissibility?

Results

We established a mechanistic dynamic model for the transmission of COIVD-19 in the entire population of China. The model is a modified version of the classical SEIR model, designed to specifically reflect the entire process of epidemic transmission [60,61,62]. The proposed model includes age-based structures and immunity levels for different vaccine doses, covering susceptibility from the latent phase through mild to severe conditions, and including isolation or hospitalization (illustrated in Fig. 1(a, b), detailed in the Methods section).

Fig. 1

Model framework schematic. a Changes in population immunity levels. The immunity level ranks from low to high as low-level immunity, basic immunity and boosted immunity. Immunity enhancement through vaccination is considered from low immunity level to high immunity level and immunity wanning is considered from high immunity level to low immunity level. Specifically, for individuals who have never been infected, immunity can only be obtained through vaccination and for individuals who have been infected, natural immunity are acquired at the boosted immunity level. b Flowchart describing COVID-19 transmission, immunity waning considering primary and reinfection

Model calibration and epidemic projection

The model fitting results for the epidemic data during the time window before the announcement of the “20 Measures” are shown in Fig. 2(b). This announcement marked the end of the dynamic zero-COVID policy in China. It is important to note that the “20 Measures” relaxed policies on case isolation, mass PCR testing, contact tracing, social distancing, and travel restrictions. The “10 Measures” released on December 7, further relaxed these policies. Therefore, the control measures changed progressively rather than disappearing all at once. The relaxation of mass PCR testing policies directly led to inaccuracies in the reported number of new cases—these numbers were certainly lower than the actual number of new infections. Therefore, data after November 10 were not fitted. The estimated values of the key epidemic parameters are presented in Supplementary Information Table S1. Additionally, for some parameters directly assumed, we also conducted uncertainty and sensitivity analyses. Please refer to the supplementary information, Figure S2, for details. We then project the epidemic waves of the Omicron variant after the adjustment of the dynamic zero-COVID policy. Note that, as we only considered three age groups, we firstly adjusted the contact matrix obtained in [63] to our setting (see the details in method section). We assume that the easing of the strict control strategy results in the free transmission of the virus, hence deduce the transmission rate by setting the reproduction number as 8.2 for Omicron variant [64]. Consequently, we show the subsequent epidemic waves in Fig. 2(a) after the adjustment of the zero-COVID policy. The projection shows that the first epidemic waves reach the peak at December 24, 2022.

Fig. 2

Model calibration and simulations. a Daily new infections from November 1, 2022, to the end of June 2023. The lighter shaded portion represents 95% confidence interval. b Effective reproduction number ${R}_{e}\left(t\right)$. The shaded area represents a 10% perturbation. It drops below the threshold of 1 on December 29, 2022. c The early-stage epidemic curve of the initial impact wave. Three durations are divided according to the policies enforced in China. Data points represent actual data, the curve reflects the model’s fitting results, and the shaded area represents a 10% data perturbation. d, f, h. Daily new infection curves (d), daily new severe infection curves (f), and daily new death curves (h) in 0–2, 3–59 and 60 + age groups, respectively. The lighter shaded portion represents 95% confidence interval. e, g, i. The proportions of total infections (TI), total severe infections (TSI), and total deaths (TD) to the total population in age group 0–2 (e), 3–59 (g) and 60 + (i), respectively

Further, we show the daily new infections, daily new incidence of severe illness and death in the three age groups in Fig. 2(d, f, h), respectively. And the statistic results of the key epidemic indexes in the three age groups and the whole population are concluded in Table 1. Figure 2 and Table 1 show that the peak of daily new infections reached 46.63 million. By the end of January 2023, about $82.9\%$ ($95\%CI: 82.3\%, 83.5\%$) of China's population had been infected, representing the cumulative infection rate of the first epidemic wave. This is consistent with the results of cumulative infection rate (CIR) from reference [65]. In addition, we can find that the peak number of daily incidences of severe infections was projected as $\text{270,435}$ ($95\%CI: 248722,\, 294635$), which would result in high pressure on medical resources. To specific, incidence of the elderly ranked as the highest, aligning with the observations. This contrast becomes even more pronounced in the case of fatalities. Similarly, the overall CFR during the first epidemic wave in China was estimated at $0.25\% \left(95\%CI:0.248\%,0.253\%\right)$, with an estimated cumulative death toll of around $2.93$ million $\left(95\%CI:2.88million-2.98million\right)$. The elderly was the most affected group, with a CFR of 1%, which is 15 times higher than that of other age groups. This conclusion is also evident from Fig. 2 (e, g, i), which shows the proportions of total infections, severe infections, and deaths for specific age groups compared to the entire population during the first wave.

Table 1 Statistic results of the key epidemic indexes for the first epidemic wave

The first wave of shocks indeed caused a fairly large wave of infections, which far exceeded the casualties caused by all local epidemics in China in the past three years. From another standpoint, the sustained implementation of the zero-COVID policy delayed this impact by three years, providing three years of health protection for many individuals. This also allowed for better preparation of medical resources (such as ICU beds), thereby preventing medical overload to some extent. After deconstructing the real-world impact of policy relaxation, our subsequent task is to analyze an improved timing for easing the zero-COVID policy.

Improving the shifting of control policies

In this section, we focus on determining an improved timing for the adjustment of policy in balancing the control of epidemics and the economic development. To this end, we firstly search a better timing for easing the zero-COVID policy during the circulation of Omicron variant, given the shifting of the immunity in the population in China. In detail, Fig. 3 (a, b) shows the time series of daily new vaccination doses and the cumulative vaccination doses in China, indicating an increase in population immunity. However, a lot of clinical evidence has been well documented showing that individuals quickly lose their immunity against infection and severe illness after vaccination [66,67,68], indicating a decrease in immunity in the population. This motivates us to search an improved timing for the easing of zero-COVID policy by leveraging the protection of the mass vaccination program. To this end, we simulate the outcomes of the first large-scale epidemic wave of Omicron variant when we try to relax the zero-COVID policy in January 2022. The setting of vaccination-related parameters and the initial conditions of the vaccination are shown in Table 2. It should be mentioned that around 579 million doses have been additionally vaccinated to the population during January 2022 to November 2022. Specifically, 41,568,000 were first doses, 59,227,000 were second doses, and 478,299,000 were booster doses.

Fig. 3

a Data on new vaccination doses in China obtained from OWID. The start time of the simulation in the three scenarios is marked with a vertical dotted line. b Data on cumulative vaccination doses in China obtained from OWID. The start time of the simulation in the three scenarios is marked with a vertical dotted line. Histogram of c. CIR, e severe infection rate, and g CFR in three scenarios. d The collected data on China’s vaccination coverage at different time points played a crucial role in calibrating the model under baseline conditions and in setting the initial population distribution for various scenarios. f Simulated outbreak curves (daily infections) for China easing its zero-COVID policy under three scenarios

Table 2 Parameter settings for three scenarios

The simulation results are shown in Fig. 3 and Table 3. It follows from Fig. 3(f) that there will also a large-scale epidemic wave when we try to ease the strict control strategies under zero-COVID policy. It can be seen that the first epidemic wave peaks a little bit later (by 5 days) compared to the baseline scenario (the estimation of the real case considered in last section, that the zero-COVID policy is eased in November 2022). The CIR is estimated as 81.64%, representing an 1% decrease compared to the baseline scenario. The peak size is 52,208,764, which is 7% less than the baseline scenario in terms of population. Moreover, we find that the peak incidence number of daily severe infections significantly decreased 68,615, where the cumulative number of severe infections is estimated as 4,999,293 with a severe infection rate of 0.434% (compare to 0.484% in the baseline scenarios). The situation regarding deaths during the first large-scale epidemic wave has improved. Specifically, the peak daily deaths decreased to 80,886, cumulative deaths fell from 2.929 million to 2.336 million, and the CFR decreased from 0.25% to 0.203%.

Table 3 Comparison of main simulation results of the first wave infection in three scenarios

We further evaluate if the strict control interventions under the zero-COVID policy can be relaxed much earlier before the emerging of the variant with high transmissibility (i.e., the Omicron variant).

To this end, we simulated in Fig. 3(f) the outcomes of the epidemic wave considering that the zero-COVID policy was relaxed in June 2021, while the Delta variant was circulating in China. Figure 3(f) shows that, despite the Delta variant's lower transmissibility compared to the Omicron variant, there is still a potential for a large-scale epidemic wave with a cumulative infection rate of 74.06%. However, the epidemic characteristics changed significantly compared to the baseline scenario. For example, the first epidemic wave reached its peak with 29.71 million daily new infections approximately 44 days later than in the baseline scenario. Although with a much small CIR, the peak incidence of severe illness was approximately 6.98 million higher than that of the baseline case for the Omicron variant, while the severe infection rate increased from 0.484% to 1.32%. Regarding of the deaths, we find that the peak number of daily deaths can increase to around 289 thousand, while the cumulative death reaches 11.13 million, which is 280% higher than the baseline scenario. The CFR has increased from 0.25% to 1.065%, which is 326% higher than the baseline scenario.

Discussion

Since the beginning of the twenty-first century, emerging infectious diseases have surfaced intermittently, posing major challenges to society. It is anticipated that new infectious diseases, referred to as “Disease X”, will continue to arise in the near future. Common characteristics have been identified among various infectious disease epidemics, notably the rapid evolution of viruses and the successive emergence of new variants. Throughout the long history of infectious disease control, we have developed a suite of measures, including non-pharmaceutical interventions (NPIs) and vaccination, to curb the spread of emerging infections. However, vaccinations alone often fail to halt the spread of these diseases, either due to difficulties in vaccine development or because mass vaccination does not achieve herd immunity owing to rapid waning immunity. Similarly, NPIs cannot stop the global spread of emerging infectious diseases due to the prohibitive costs associated with implementing stringent control measures. Therefore, it is crucial to design a comprehensive and improved control policy that leverages the synergistic effects of NPIs and vaccination to combat Disease X to curb the epidemics at the most effective and cost-efficient way. It would be feasible to devise a general and adaptive control policy that draws on past experiences in tackling previous emerging infectious diseases. Here, adaptive strategy means an appropriate shifting of control policy following the quick evolving of the virus in terms of its transmissibility and the pathogenicity.

Considering that COVID-19 is regarded as the first instance of Disease X [18,19,20,21,22], and given the significant shifts in China's control policy for COVID-19, this study aims to develop a general and adaptive policy based on China’s experiences. We first projected the epidemic waves of the Omicron variant following the easing of the zero-COVID policy in November 2022 using a compartment model. This serves as the baseline for comparing whether there is a better timing for easing the policy. In the case of weak or even no control interventions, it can be seen that there is definitely a large-scale epidemic wave no matter when the zero-COVID policy is relaxed if the basic reproduction number exceeds 8. Note that, the large-scale wave occurs even when the vaccination coverage already exceeds 80% in China, which is owing to that a large ratio of population loss their immunity already before the vaccination coverage reaches a high level. Our results indicate that the easing of zero-COVID policy at January, 2022, can not only help to reduce the outbreak size and the incidence of severe infections but also save a significant cost of boosting dose of the mass vaccination comparing to the baseline case. Therefore, it is important to choose an improved timing for releasing the strict control intervention by leveraging the protection of mass vaccination. And our results shows that an improved timing for easing the strict control intervention can not only reduce the outbreak size of the epidemic wave but also greatly saving the control-cost of additional vaccination.

The natural question then arises: could the strict control interventions have been relaxed much earlier? To explore this, we simulated the outcomes of the epidemic if the zero-COVID policy had been lifted in January 2021, during the circulation of the Delta variant in China. The results indicate that there would still have been a large-scale epidemic wave, with an accumulative infection rate reaching 89% of the baseline scenario. However, the peak size of new infections would have decreased to 47% of the baseline, due to the Delta variant's lower transmissibility compared to that of the Omicron variant. Despite these findings, the situation would have worsened significantly, with severe infection rates and case fatality rates (CFR) increasing by 71% and 226%, respectively, compared to the baseline case. Given China's large population, this scenario would have led to a substantial increase in the number of deaths. Considering that China had 17 rounds of local outbreaks during 225 days of Delta variant circulation, with the highest peak at only 180 cases [58], it is suggested that quick case clearance should be pursued for variants with lower transmissibility. Each outbreak could potentially be controlled within a short period, such as about one month, based on China's experience with the Delta variant.

In conclusion, managing an emerging infectious disease (Disease X) requires comprehensive consideration and monitoring of several factors: the transmissibility of viral variants, the upper limit of NPI effectiveness, the speed of vaccine development and administration, and the waning efficacy of vaccines. These factors are crucial for formulating an adaptive control strategy. The proposed general and adaptive policy can be summarized as follows: initially, highly stringent control measures (such as a zero-COVID policy) may be employed to halt each outbreak when the virus's basic reproduction number is below 5. If stopping the spread quickly proves too costly, it may be necessary to relax these stringent measures and transition to normalized control strategies. This shift should occur at a more appropriate time, leveraging the protective effects of mass vaccination to maintain control over the disease spread.

There are still several limitations in this study. Firstly, since the model design did not specifically account for excess mortality among severe patients due to insufficient healthcare capacity, the estimated mortality figures in this paper are somewhat conservative. Secondly, the vaccine-related data we have is still incomplete, and missing vaccination data for each age group at certain time points limits our ability to perform a more detailed analysis. Lastly, although we performed three scenario analyses as preliminary experiments, designing dynamic and comprehensive strategies—such as NPIs and vaccination plans incorporating game theory, cost-effectiveness, and sensitive variant monitoring—remains a challenge. This will require further research.

Methods

Model formula

Based on a Susceptible- Exposed- Asymptomatic- Infectious- Quarantined- Hospitalized- Recovered- Death model structure (SEAIQHRD), we developed a dynamic model of COVID-19 transmission in the widespread outbreak caused by the SARS-CoV-2 Omicron variant in China, and the flow diagram is shown in Fig. 1. Considering the different levels of social activities of different age groups in the modeling, the population was divided into three age groups: 0–2 years old, 3–59 years old, and 60 years old and above. According to the dose of vaccine and individual immunity level, the population was further divided into four groups: unvaccinated group, one-dose vaccination group, two-dose vaccination group and booster vaccination group. Since almost everyone in China had never been infected with the COVID-19 before the easing of the ‘zero-COVID policy’, we divide the infection process into the first round of infection and reinfection. Therefore, each state variable or compartment variable has three superscripts $i, j,$ and $k$, which represent age group, immunity level group, and whether it is the first round of infection. For example, ${S}^{\text{2,1},0}$ represents people in the second age group and the second immunity group who have been infected and are susceptible again. It is worth noting that the division of three age groups can not only highlight the heterogeneity of the population, but also avoid the model structure from being too complex. And this is also related to the division of immunity level groups, because according to national policy, children aged 0–2 years are generally not vaccinated, while those over 60 years old are the national key protected groups. Meanwhile, the vaccination data available to us is categorized by age only for individuals aged 60 and above. The data for the 3–59 age group lacks further stratification. In fact, vaccination data largely determines how we set up the three age groups. In addition, although it appears that all three age groups fall into four immunity levels (a total of 12 different combinations), the first age group (0–2 years old) is always entirely in the first immunity level group (unvaccinated group) during the first round of infection. Immunity waning is considered to reflect the movement of people between different immunity level groups (e.g., ${S}^{i,\text{3,0}}$ to ${S}^{i,\text{1,0}}$). In addition, it is also used to reflect the movement of recovered people to susceptible people (e.g., ${R}^{i,\text{1,0}}$ to ${S}^{i,\text{1,1}}$). According to the individual state of infection, the population is divided into susceptible ($S$), exposed ($E$), infectious ($A, {I}_{m},{I}_{s}$), quarantine (${Q}_{A},{Q}_{m},H$), recovered ($R$) and dead ($D$) classes. We distinguish between asymptomatic infection ($A$), mild infection (${I}_{m}$), severe infection (${I}_{s}$), and also distinguish between the corresponding home quarantine (${Q}_{A}$, ${Q}_{m}$) and severe hospitalization ($H$).

We first sort out the dynamics of individual infections without superscripts. The transmission occurs when the susceptible population ($S$) closely contact with infectious individuals ($A$, ${I}_{m}$ and ${I}_{S}$) and become exposed class ($E$) with a transmission probability $\beta$ per contact. Two types of infectious individuals have calibration coefficients ${\varsigma }_{A}$ and ${\varsigma }_{{I}_{m}}$ for their infectivity relative to severe patients, respectively. Individuals in exposed ($E$) will enter into the asymptomatic population ($A$) with the transition rate $\left(1-\rho \right)\sigma$, where $\rho$ is the probability of developing symptoms and $1/\sigma$ is the incubation period. The other individuals in exposed will enter into the mild infection population (${I}_{m}$) and severe infection population (${I}_{S}$) with the transition rate $\left(1-\theta \right)\rho \sigma$ and $\theta \rho \sigma$, respectively. $\theta$ is the probability of severe illness among symptomatic patients. Further, these three types of infected individuals will enter into quarantined asymptomatic population (${Q}_{A}$), quarantined mild infection population (${Q}_{m}$) and hospitalized and isolated population ($H$) according to the home isolation rate (${\delta }_{A}$ and ${\delta }_{Im}$) and hospitalization rate (${\delta }_{Is}$) respectively. Among them, critically ill patients will be transferred to ${Q}_{m}$ with a discharge rate of ${\tau }_{H}$. All asymptomatic and mild symptom patients will move to the recovery class with recovery rates ${\gamma }_{A}$, ${\gamma }_{Im}$, ${\gamma }_{{Q}_{A}}$, and ${\gamma }_{{Q}_{m}}$ respectively. A disease-related death for the hospitalized population is considered with severe mortality rate $d$.

Finally, we assume that individuals with medium and high immunity levels (${S}^{i,\text{2,0}}$ and ${S}^{i,3,0}$) will move to the low immunity level group (${S}^{i,1,0}$) at rates ${\upomega }^{i,\text{2,0}}$ and ${\upomega }^{i,3,0}$ respectively. After experiencing the first round of infection, the recoveries (${R}^{i,j,0},j=\text{0,1},\text{2,3}$) will move to the susceptible class (${S}^{i,j,1}, j=\text{1,2},3$) at rates of immunity waning (${\omega }^{i,j,0}, j=\text{0,1},\text{2,3}$). It should be noted that since individuals in the first immunity level group acquire natural immunity after being infected and recovering, they will enter ${S}^{i,\text{1,1}}$ instead of ${S}^{i,\text{0,1}}$. The setting for immunity waning during the reinfection process is consistent with that for the first round of infection.

A more detailed description of the parameters and variables is given in Table S1. It is important to emphasize that, for the consistency of model representation, the superscripts of the compartment variables correspond one-to-one with the superscripts of each parameter. However, in the actual simulation, apart from the contact matrix $C$, which is only related to age groups, the other parameters are only related to immunity levels. Therefore, in Table S1, for clarity, we denote each parameter with a subscript and use the superscript ‘-’ to distinguish between reinfection and the first round of infection.

$$\begin{array}{l}\frac{d{S}^{i,\text{0,0}}}{dt}=-{\Lambda }^{i,\text{0,0}}\\ \frac{d{S}^{i,\text{1,0}}}{dt}=-{\Lambda }^{i,\text{1,0}}+{\upomega }^{i,\text{2,0}}{S}^{i,\text{2,0}}+{\upomega }^{i,\text{3,0}}{S}^{i,\text{3,0}}\\ \frac{d{S}^{i,\text{2,0}}}{dt}=-{\Lambda }^{i,\text{2,0}}-{\upomega }^{i,\text{2,0}}{S}^{i,\text{2,0}}\\ \frac{d{S}^{i,\text{3,0}}}{dt}=-{\Lambda }^{i,\text{3,0}}-{\upomega }^{i,\text{3,0}}{S}^{i,\text{3,0}}\\ \frac{d{S}^{i,\text{0,1}}}{dt}=-{\Lambda }^{i,\text{0,1}}\\ \frac{d{S}^{i,\text{1,1}}}{dt}=-{\Lambda }^{i,\text{1,1}}+{\upomega }^{i,\text{0,0}}{R}^{i,\text{0,0}}+{\upomega }^{i,\text{1,0}}{R}^{i,\text{1,0}}+{\upomega }^{i,\text{0,1}}{R}^{i,\text{0,1}}+{\upomega }^{i,\text{1,1}}{R}^{i,\text{1,1}}+{\upomega }^{i,\text{2,1}}{S}^{i,\text{2,1}}+{\upomega }^{i,\text{3,1}}{S}^{i,\text{3,1}}\\ \frac{d{S}^{i,\text{2,1}}}{dt}=-{\Lambda }^{i,\text{2,1}}+{\upomega }^{i,\text{2,0}}{R}^{i,\text{2,0}}+{\upomega }^{i,\text{2,1}}{R}^{i,\text{2,1}}-{\upomega }^{i,\text{2,1}}{S}^{i,\text{2,1}}\\ \frac{d{S}^{i,\text{3,1}}}{dt}=-{\Lambda }^{i,\text{3,1}}+{\upomega }^{i,\text{3,0}}{R}^{i,\text{3,0}}+{\upomega }^{i,\text{3,1}}{R}^{i,\text{3,1}}-{\upomega }^{i,\text{3,1}}{S}^{i,\text{3,1}}\\ \frac{d{E}^{i,j,k}}{dt}={\Lambda }^{i,j,k}-{\sigma }^{i,j,k}{E}^{i,j,k}\\ \frac{d{A}^{i,j,k}}{dt}=\left(1-{\rho }^{i,j,k}\right){\sigma }^{i,j,k}{E}^{i,j,k}-({\gamma }_{A}^{i,j,k}+{\delta }_{A}^{i,j,k}){A}^{i,j,k}\\ \frac{d{I}_{m}^{i,j,k}}{dt}=\left(1-{\theta }^{i,j,k}\right){\rho }^{i,j,k}{\sigma }^{i,j,k}{E}^{i,j,k}-({\gamma }_{Im}^{i,j,k}+{\delta }_{Im}^{i,j,k}){I}_{m}^{i,j,k}\\ \frac{d{I}_{s}^{i,j,k}}{dt}={\theta }^{i,j,k}{\rho }^{i,j,k}{\sigma }^{i,j,k}{E}^{i,j,k}-{\delta }_{Is}^{i,j,k}{I}_{s}^{i,j,k}\\ \frac{d{Q}_{A}^{i,j,k}}{dt}={\delta }_{A}^{i,j,k}{A}^{i,j,k}-{\gamma }_{{Q}_{A}}^{i,j,k}{Q}_{A}^{i,j,k}\\ \frac{d{Q}_{m}^{i,j,k}}{dt}={\delta }_{{I}_{m}}^{i,j,k}{I}_{m}^{i,j,k}+{\tau }_{H}^{i,j,k}{H}^{i,j,k}-{\gamma }_{{Q}_{m}}^{i,j,k}{Q}_{m}^{i,j,k}\\ \frac{d{H}^{i,j,k}}{dt}={\delta }_{Is}^{i,j,k}{I}_{s}^{i,j,k}-({\tau }_{H}^{i,j,k}+{d}^{i,j,k}){H}^{i,j,k}\\ \frac{d{R}^{i,j,k}}{dt}={\gamma }_{A}^{i,j,k}{A}^{i,j,k}+{\gamma }_{Im}^{i,j,k}{I}_{m}^{i,j,k}+{\gamma }_{{Q}_{A}}^{i,j,k}{Q}_{A}^{i,j,k}+{\gamma }_{{Q}_{m}}^{i,j,k}{Q}_{m}^{i,j,k}-{\upomega }^{i,j,k}{R}^{i,j,k}\\ \frac{d{D}^{i,j,k}}{dt}={d}^{i,j,k}{H}^{i,j,k}\end{array}$$

(1)

Particularly, the infection rate ${\Lambda }^{i,j,k}$ was given as follows.

$${\Lambda }^{i,j,k}=\sum\nolimits_{b=1}^{3}\frac{{c}_{ib}{\beta }^{i,j,k}{S}^{i,j,k}}{{N}^{b}}\sum\nolimits_{n=0}^{3}\sum\nolimits_{l=0}^{1}\left({\varsigma}_{A}^{b,n,l}{A}^{b,n,l}+{\varsigma }_{{I}_{m}}^{b,n,l}{I}_{m}^{b,n,l}+{I}_{s}^{b,n,l}\right)$$

(2)

Specifically, ${\Lambda }^{i,j,k}$ represents individuals who have moved out of the susceptible group (${\text{S}}^{i,j,k}$) due to being infected after coming into contact with an infected person. Since our contact matrix is derived based on three age groups, for a specific age group $\left(i\right)$, its contact rate with each age group is represented as ${c}_{ib}, (b=\text{1,2},3$). When contacting age group $b$, it also includes divisions based on four immunity levels and whether the individual is the first round of infection. Here, ${N}^{b}$​ represents the total population of age group $b$. It is placed in the denominator to represent the probability of encountering individuals from the aforementioned categories during each contact.

Basic reproduction number

The basic reproduction number ${R}_{0}$ plays a crucial role in disease control and the stability of the disease-free equilibrium (DFE) [69]. When ${R}_{0}>1$, the disease spreads within the population, rendering the DFE unstable as the rate of new infections surpasses the rate of recoveries. Conversely, when ${R}_{0}<1$, the disease cannot sustain itself, leading to a stable DFE and eventual elimination of the infection over time. This threshold concept is fundamental for designing control strategies, such as achieving sufficient vaccination coverage or implementing interventions to reduce ${R}_{0}$ below 1, thereby eradicating the disease. General epidemic control strategies aim to minimize ${R}_{0}$ to ensure the disease’s extinction.

The basic reproduction number of the above model is calculated using the method of the next-generation matrix [70]. First, we need to identify two vectors $\mathcal{F}$ and $\mathcal{V}$. $\mathcal{F}=\left[\begin{array}{c}{\Lambda }_{24\times 1}\\ {0}_{72\times 1}\end{array}\right]$, ${\Lambda }_{24\times 1}$ represents a column vector of 24 dimensions, whose elements are ${\Lambda }^{i,j,k}$,

$$\mathcal{V}=\left[\begin{array}{c}\begin{array}{c}{\sigma }_{24\times 24}{E}_{24\times 1}\\ \left({{\gamma }_{A}}_{24\times 24}+{{\delta }_{A}}_{24\times 24}\right){A}_{24\times 1}-\left({I}_{24\times 24}-{\rho }_{24\times 24}\right){\sigma }_{24\times 24}{E}_{24\times 1}\end{array}\\ \begin{array}{c}\left({{\gamma }_{{I}_{m}}}_{24\times 24}+{{\delta }_{{I}_{m}}}_{24\times 24}\right){{I}_{m}}_{24\times 1}-\left({I}_{24\times 24}-{\varrho }_{24\times 24}\right){\rho }_{24\times 24}{\sigma }_{24\times 24}{E}_{24\times 1}\\ \left({{\phi }_{{I}_{s}}}_{24\times 24}+{{\delta }_{{I}_{s}}}_{24\times 24}\right){{I}_{s}}_{24\times 1}-{\varrho }_{24\times 24}{\rho }_{24\times 24}{\sigma }_{24\times 24}{E}_{24\times 24}\end{array}\end{array}\right]$$

Here, ${E}_{24\times 1}$, ${A}_{24\times 1}$, ${{I}_{m}}_{24\times 1}$, ${{I}_{s}}_{24\times 1}$ represent column vectors of 24 dimensions respectively. Its elements are ${E}^{i,j,k}$, ${A}^{i,j,k}$, ${I}_{m}^{i,j,k}$, ${I}_{s}^{i,j,k}$. ${I}_{24\times 24}$ is a 24-dimensional identity matrix. ${\sigma }_{24\times 24}$, ${{\gamma }_{A}}_{24\times 24}$, ${{\gamma }_{{I}_{m}}}_{24\times 24}$, ${{\delta }_{A}}_{24\times 24}$, ${{\delta }_{{I}_{m}}}_{24\times 24}$, ${{\delta }_{{I}_{s}}}_{24\times 24}$, ${{\phi }_{{I}_{s}}}_{24\times 24}$, ${\rho }_{24\times 24}$, ${\varrho }_{24\times 24}$ represents a diagonal matrix of 24 dimensions, whose diagonal elements are ${\sigma }^{i,j,k}$, ${\gamma }_{A}^{i,j,k}$, ${\gamma }_{{I}_{m}}^{i,j,k}$, ${\delta }_{A}^{i,j,k}$, ${\delta }_{{I}_{m}}^{i,j,k}$, ${\delta }_{{I}_{s}}^{i,j,k}$, ${\phi }_{{I}_{s}}^{i,j,k}$, ${\rho }^{i,j,k}$, ${\varrho }^{i,j,k}$. For example,

$${\sigma }_{24\times 24}=\left[\begin{array}{cc}{\sigma }_{12\times 12}^{k=0}& {0}_{12\times 12}\\ {0}_{12\times 12}& {\sigma }_{12\times 12}^{k=1}\end{array}\right],$$

$${\sigma }_{12\times 12}^{k=0}=\left[\begin{array}{ccc}\begin{array}{ccc}{\sigma }^{\text{1,0},0}& 0& 0\\ 0& {\sigma }^{\text{2,0},0}& 0\\ 0& 0& {\sigma }^{\text{3,0},0}\end{array}& \cdots & \begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\\ \vdots & \ddots & \vdots \\ \begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}& \cdots & \begin{array}{ccc}{\sigma }^{\text{1,4},0}& 0& 0\\ 0& {\sigma }^{\text{2,4},0}& 0\\ 0& 0& {\sigma }^{\text{3,4},0}\end{array}\end{array}\right],$$

$${\sigma }_{12\times 12}^{k=1}=\left[\begin{array}{ccc}\begin{array}{ccc}{\sigma }^{\text{1,0},1}& 0& 0\\ 0& {\sigma }^{\text{2,0},1}& 0\\ 0& 0& {\sigma }^{\text{3,0},1}\end{array}& \cdots & \begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\\ \vdots & \ddots & \vdots \\ \begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}& \cdots & \begin{array}{ccc}{\sigma }^{\text{1,4},1}& 0& 0\\ 0& {\sigma }^{\text{2,4},1}& 0\\ 0& 0& {\sigma }^{\text{3,4},1}\end{array}\end{array}\right]$$

The construction of the remaining matrices follows a similar approach.

So, we have

$$F=\left[\begin{array}{ccc}\begin{array}{c}{0}_{24\times 24}\\ {0}_{24\times 24}\\ \begin{array}{c}{0}_{24\times 24}\\ {0}_{24\times 24}\end{array}\end{array}& \begin{array}{c}{F}_{\text{1,2}}\\ {0}_{24\times 24}\\ \begin{array}{c}{0}_{24\times 24}\\ {0}_{24\times 24}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}{F}_{\text{1,3}}\\ {0}_{24\times 24}\\ \begin{array}{c}{0}_{24\times 24}\\ {0}_{24\times 24}\end{array}\end{array}& \begin{array}{c}{F}_{\text{1,4}}\\ {0}_{24\times 24}\\ \begin{array}{c}{0}_{24\times 24}\\ {0}_{24\times 24}\end{array}\end{array}\end{array}\end{array}\right], V=\left[\begin{array}{ccc}\begin{array}{c}{V}_{\text{1,1}}\\ {V}_{\text{2,1}}\\ \begin{array}{c}{V}_{\text{3,1}}\\ {V}_{\text{4,1}}\end{array}\end{array}& \begin{array}{c}{0}_{24\times 24}\\ {V}_{\text{2,2}}\\ \begin{array}{c}{0}_{24\times 24}\\ {0}_{24\times 24}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}{0}_{24\times 24}\\ {0}_{24\times 24}\\ \begin{array}{c}{V}_{\text{3,3}}\\ {0}_{24\times 24}\end{array}\end{array}& \begin{array}{c}{0}_{24\times 24}\\ {0}_{24\times 24}\\ \begin{array}{c}{0}_{24\times 24}\\ {V}_{\text{4,4}}\end{array}\end{array}\end{array}\end{array}\right]$$

Here, ${F}_{\text{1,2}}$, ${F}_{\text{1,3}}$, and ${F}_{\text{1,4}}$ represent the Jacobian matrices composed of the partial derivatives of the ${\Lambda }_{24\times 1}$ vector with respect to the ${A}_{24\times 1}$, ${{I}_{m}}_{24\times 1}$, and ${{I}_{s}}_{24\times 1}$ vectors, respectively. Taking ${F}_{\text{1,2}}$ as an example,

$${F}_{\text{1,2}}=\frac{\partial {{\Lambda }_{24\times 1}}^{T}}{\partial {{A}_{24\times 1}}^{T}}=\frac{\partial \left({\Lambda }^{\text{1,0},0},{\Lambda }^{\text{2,0},0},{\Lambda }^{\text{3,0},0},\dots ,{\Lambda }^{\text{1,4},0},{\Lambda }^{\text{2,4},0},{\Lambda }^{\text{3,4},0},{\Lambda }^{\text{1,0},1},{\Lambda }^{\text{2,0},1},{\Lambda }^{\text{3,0},1},\dots ,{\Lambda }^{\text{1,4},1},{\Lambda }^{\text{2,4},1},{\Lambda }^{\text{3,4},1}\right)}{\partial \left({A}^{\text{1,0},0},{A}^{\text{2,0},0},{A}^{\text{3,0},0},\dots ,{A}^{\text{1,4},0},{A}^{\text{2,4},0},{A}^{\text{3,4},0},{A}^{\text{1,0},1},{A}^{\text{2,0},1},{A}^{\text{3,0},1},\dots ,{A}^{\text{1,4},1},{A}^{\text{2,4},1},{A}^{\text{3,4},1}\right)}{|}_{24\times 24}$$

Here, ${\left(\bullet \right)}^{T}$ denotes the transpose of the matrix. ${F}_{\text{1,3}}$ and ${F}_{\text{1,4}}$ follow a similar structure.

${V}_{\text{1,1}}$, ${V}_{\text{2,1}}$, ${V}_{\text{2,2}}$, ${V}_{\text{3,1}}$, ${V}_{\text{3,3}}$, ${V}_{\text{4,1}}$ and ${V}_{\text{4,4}}$ are all diagonal matrices of 24 × 24, and

$${V}_{\text{1,1}}={\sigma }_{24\times 24, } {V}_{\text{2,1}}=-\left({I}_{24\times 24}-{\rho }_{24\times 24}\right){\sigma }_{24\times 24}, {V}_{\text{2,2}}={{\gamma }_{A}}_{24\times 24}+{{\delta }_{A}}_{24\times 24},$$

$${V}_{\text{3,1}}=-\left({I}_{24\times 24}-{\varrho }_{24\times 24}\right){\rho }_{24\times 24}{\sigma }_{24\times 24}, {V}_{\text{3,3}}={{\gamma }_{{I}_{m}}}_{24\times 24}+{{\delta }_{{I}_{m}}}_{24\times 24},$$

$${V}_{\text{4,1}}={\varrho }_{24\times 24}{\rho }_{24\times 24}{\sigma }_{24\times 24}, {V}_{\text{4,4}}={{\phi }_{{I}_{s}}}_{24\times 24}+{{\delta }_{{I}_{s}}}_{24\times 24}$$

It can be calculated that,

$${V}^{-1}=\left[\begin{array}{ccc}\begin{array}{c}{V}_{\text{1,1}}^{-1}\\ {-{V}_{\text{1,1}}^{-1}V}_{\text{2,1}}{V}_{\text{2,2}}^{-1}\\ \begin{array}{c}{-{V}_{\text{1,1}}^{-1}V}_{\text{3,1}}{V}_{\text{3,3}}^{-1}\\ {-{V}_{\text{1,1}}^{-1}V}_{\text{4,1}}{V}_{\text{4,4}}^{-1}\end{array}\end{array}& \begin{array}{c}{0}_{24\times 24}\\ {V}_{\text{2,2}}^{-1}\\ \begin{array}{c}{0}_{24\times 24}\\ {0}_{24\times 24}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}{0}_{24\times 24}\\ {0}_{24\times 24}\\ \begin{array}{c}{V}_{\text{3,3}}^{-1}\\ {0}_{24\times 24}\end{array}\end{array}& \begin{array}{c}{0}_{24\times 24}\\ {0}_{24\times 24}\\ \begin{array}{c}{0}_{24\times 24}\\ {V}_{\text{4,4}}^{-1}\end{array}\end{array}\end{array}\end{array}\right],$$

$$F{V}^{-1}=\left[\begin{array}{ccc}\begin{array}{c}{G}_{\text{1,1}}\\ {0}_{24\times 24}\\ \begin{array}{c}{0}_{24\times 24}\\ {0}_{24\times 24}\end{array}\end{array}& \begin{array}{c}{G}_{\text{1,2}}\\ {0}_{24\times 24}\\ \begin{array}{c}{0}_{24\times 24}\\ {0}_{24\times 24}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}{G}_{\text{1,3}}\\ {0}_{24\times 24}\\ \begin{array}{c}{0}_{24\times 24}\\ {0}_{24\times 24}\end{array}\end{array}& \begin{array}{c}{G}_{\text{1,3}}\\ {0}_{24\times 24}\\ \begin{array}{c}{0}_{24\times 24}\\ {0}_{24\times 24}\end{array}\end{array}\end{array}\end{array}\right],$$

Then, the basic reproduction number is ${R}_{0}=\mathcal{R}\left(F{V}^{-1}{|}_{t=0}\right)=\mathcal{R}\left({G}_{\text{1,1}}{|}_{t=0}\right)$, ${G}_{\text{1,1}}=-{F}_{\text{1,2}}{V}_{\text{1,1}}^{-1}{V}_{\text{2,1}}{V}_{\text{2,2}}^{-1}-{F}_{\text{1,3}}{V}_{\text{1,1}}^{-1}{V}_{\text{3,1}}{V}_{\text{3,3}}^{-1}-{F}_{\text{1,4}}{V}_{\text{1,1}}^{-1}{V}_{\text{4,1}}{V}_{\text{4,4}}^{-1}$. And the effective reproduction number ${R}_{t}=\mathcal{R}\left(F{V}^{-1}{|}_{t}\right)$.

Contact matrix

Since we obtained an 84-age group contact matrix from the references, and what we need for this study is a contact matrix among three age groups, it is necessary to process the original contact matrix. ${c}_{ij}$ represents the per capita contact rate between age group $i$ and age group $j$ ($\text{i},\text{j}=\text{1,2},\dots ,84$), defined as.

$${c}_{ij}={f}_{ij}/{n}_{i},$$

Where ${f}_{ij}$ represents the total number of contacts between age group $i$ and age group $j$, where ${f}_{ij}={f}_{ji}$ is met, and ${n}_{i}$ represents the population of age group $i$. Below, we divide the 84 age groups into three large age groups, where ${F}_{uv}$ corresponds to the total number of contacts between age group $u$ and age group $v$ ($\text{u},\text{v}=\text{1,2},3$). The specific calculation process is as follows:

$${F}_{11}=\sum_{i=1}^{2}\sum_{j=1}^{2}{c}_{ij}{n}_{i}, {F}_{12}=\sum_{i=1}^{2}\sum_{j=3}^{59}{c}_{ij}{n}_{i}, {F}_{13}=\sum_{i=1}^{2}\sum_{j=60}^{84}{c}_{ij}{n}_{i},$$

$${F}_{21}=\sum_{i=3}^{59}\sum_{j=1}^{2}{c}_{ij}{n}_{i}, {F}_{22}=\sum_{i=3}^{59}\sum_{j=3}^{59}{c}_{ij}{n}_{i}, {F}_{23}=\sum_{i=3}^{59}\sum_{j=60}^{84}{c}_{ij}{n}_{i},$$

$${F}_{31}=\sum_{i=60}^{84}\sum_{j=1}^{2}{c}_{ij}{n}_{i}, {F}_{32}=\sum_{i=60}^{84}\sum_{j=3}^{59}{c}_{ij}{n}_{i}, {F}_{33}=\sum_{i=60}^{84}\sum_{j=60}^{84}{c}_{ij}{n}_{i},$$

$${N}_{1}=\sum_{i=1}^{2}{n}_{i}, {N}_{2}=\sum_{i=3}^{59}{n}_{i}, {N}_{3}=\sum_{i=60}^{84}{n}_{i},$$

$${C}_{uv}={F}_{uv}/{N}_{u}$$

Here, ${n}_{i}$ and ${c}_{ij}$ are available data [63]. The specific values of the contact matrix ${C}_{uv}$ is given in Table 1.

Data

The multi-source epidemiological data on the local COVID-19 outbreak in China caused by the SARS-CoV-2 Omicron variant were obtained from the Chinese Center for Disease Control and Prevention. Detailed information is provided in Supplementary information Figure S1 [71]. The Chinese government announced two press conferences on November 10 and December 7, 2022, respectively, addressing the new measures for COVID-19 prevention and control [72]. This marked a significant shift away from the zero-COVID policy. Consequently, we obtained the reported daily case numbers from November 1 to December 23. It is important to note that, Due to the substantial reduction in the scope of mass Polymerase Chain Reaction (PCR) testing mandated by the 20 measures, the reported data no longer align with the actual number of new cases since November 10, 2022.

Additionally, according to national COVID-19 infection statistics from the Chinese Center for Disease Control and Prevention, the first wave reached its peak between December 22 and December 25 [71]. Despite the absence of mass PCR testing in China, the absolute peak value of new infections is unreliable; however, the relative timing of reaching the peak is considered reliable. Therefore, we utilized data from November 1 to November 10, along with the peak timings of the first wave, for model fitting.

The vaccine administration data were sourced not only from the official website of the Chinese National Health Commission but also from Our World in Data (OWID) [58]; however, only partial datasets were available from both sources. The fitting of baseline conditions and the corresponding time-point vaccine data for the two additional scenarios are presented in Fig. 3(d) in the main text [71, 73,74,75,76,77]. It is important to note that due to the deletion or hiding of much information about cases and vaccine administration on the official website of the National Health Commission of China, we can only rely on citations from other websites that report the commission’s announcements as indirect references.

Model fitting process

We calibrated model (1) using the least squares method (LSM) with COVID-19 data from China, all implemented using MATLAB software.

First, we utilized the calculation formula for the basic reproduction number (${R}_{0}$), assuming that all susceptibles (i.e., the entire population) were unvaccinated ($\sum_{i}{S}^{i,\text{0,0}}\left(0\right)={N}_{all}$) and excluding the effects of all control measures. With all other parameters held constant, we estimated the original transmissibility ${\beta }_{0}$. Then we use LSM only in order to estimate the initial susceptible population ${S}^{i,j,0}\left(0\right)$, initial mild infection population of primary infection ${I}_{m}^{i,j,0}\left(0\right)$, the calibration coefficient ​${f}_{sd}$ for the transmissibility ${\beta }_{0}$, which reflects the reduction in transmissibility due to social distancing measures.

We fixed some parameters in model (1) by reviewing literature, which included the incubation period ($1/\sigma$), rate of immunity waning (${\omega }^{i,j,0}$), the vaccine effectiveness, contact matrix ($C$), the probability of developing symptoms in unvaccinated exposed individuals at three ages (${\rho }_{0}$), the probability of severe symptoms among unvaccinated people in three age groups (${\theta }_{0}$), the proportion of severe ill unvaccinated patients admitted to the ICU in three age groups (${\alpha }_{0}$), and severe mortality among unvaccinated people in three age groups (${d}_{0}$). Additionally, the correction coefficients for the force of infection of asymptomatic and mildly symptomatic individuals relative to severe cases,

${\varsigma }_{A}$​ and ${\varsigma }_{m}$​, are both assumed to be 1. The immunity waning rate (${\omega }^{i,j,1}$)​ after the first infection is assumed to be the same as (${\omega }^{i,j,0}$)​. Due to the relaxation of home quarantine policies, we set the home quarantine rate to 0. Additionally, based on empirical data, we made assumptions regarding the hospitalization rate, the average length of hospital stay for severe patients, and the recovery rates. Sensitivity analyses were conducted for all the assumed parameters mentioned above. Please refer to the Supplementary Information section on sensitivity analysis.

Except for the initial susceptible population and the mildly infected population, the initial values of all other compartments in the model are set to zero, as shown in Table S1. The initial distribution of susceptible individuals is set based on the vaccination data for the partially vaccinated population and the elderly that we obtained, as shown in Fig. 3(d). Additionally, supported by the estimation in reference [65], we assume that 15% of the population did not participate in the first wave of the pandemic.

$$Loss={\sum }_{t=Nov 1st}^{Nov 10th}{\left({Data}_{t}-{\sum }_{i,j,k}{\sigma }^{i,j,k}{E}_{t}^{i,j,k}\right)}^{2}+\psi {\left({Data}_{peak time}-{Model}_{peak time}\right)}^{2}$$

We defined the loss function for the LSM method as shown above. The loss function consists of two parts: the first part represents the cumulative residuals of daily new cases, and the second part represents the residuals for the peak time of the first wave. In MATLAB, we use the ‘fmincon’ function to solve constrained minimization problems. The Loss function it calls is set as described. The magnitude of the first part, ${Data}_{t}$ for the 10-day period, is different from the magnitude of the second part, ${Data}_{peak time}$. To help ‘fmincon’ converge more effectively to the optimal solution, we need to assign weights to the two parts of the Loss. Through numerical experiments, we set $\psi$ to 10, as detailed in the Data Availability section.

Scenario setting for the improved easing time

We constructed a model that takes into account different immunity levels and estimated the initial population sizes for various groups. Combining vaccination data curves obtained from official website of China National Center for Disease Control and Prevention and OWID [59], we can essentially provide the distribution of population sizes with differing immunity levels at certain time points over three years. Consequently, the model can obtain initial values at any given moment, enabling simulations of the epidemic scenario in China under different timeframes for easing the zero-COVID policy.

The new vaccination doses and the cumulative vaccination doses in China were shown in Fig. 3(a, b), it can be found that the vaccination rate exhibits a slow-fast-slow trend, resembling an S-shaped curve. Within the vaccination intervals observed over the last two years, it is clear that, from the standpoint of the entire population, the highest average immunity level is currently occurring during the most rapid phase of vaccine administration. These can provide us with some references for setting up scenarios.

Below, we illustrate how to set the relevant parameters for a scenario using an example. If we choose June 2021 as a scenario, according to the data we collected shown in Fig. 3(d), the proportion of people newly vaccinated is 43.62% [71, 73,74,75,76,77]. It is noteworthy that the transmissibility of the Delta variant was not as potent as that of the Omicron variant. In our simulations, once the policy is eased, the vaccination program also remains at this moment. We depict the initial distribution of susceptible individuals based on the age distribution of the Chinese population obtained from the National Bureau of Statistics of China [78]. At that point, China only had a population consisting of individuals without immunity and those with a single dose of vaccine. Therefore, the population newly vaccinated with one dose, 43.62%, is divided by age proportion to obtain the number of individuals with one-dose immunity. Then, subtracting this number from the total population of each age group gives the number of individuals with zero-dose immunity. Note that infants and toddlers aged 0–2 years are assumed not to be vaccinated. To better reflect reality, it is necessary to adjust the key parameters for simulating the Delta variant. According to the literature [79,80,81], we can assume that the regeneration number of the Delta variant ${R}_{\delta }$ is about 5.08 and the case fatality rate of it is 3 times that of Omicron, as reflected in the parameter for severe mortality $d$.

Similarly, if we choose January 2022 as another scenario, approximately 2.07 billion doses of vaccine are administered within six months. As of January 6, 2022, approximately 88.5% of the population had received at least one vaccine doses, and 85% of the population had received a second dose. Furthermore, the universal administration of the second dose commenced in June 2021 [73]. Thus, in this scenario, the proportion of the population at the baseline immunity level is 85%, while the population with zero-dose immunity is 11.5%. Hence, the count of individuals with single-dose immunity can be derived accordingly. Finally, specific initial population distributions are provided based on different age groups’ proportions.

In the scenario analysis, a single simulation is conducted using only the average values of the estimated parameters.