Impact of dynamic self-protection intensity on the COVID-19 pandemic: a case study in Shenzhen based on medical resources

ABSTRACT

Owing to the high uptake of vaccines and the low mortality of the Omicron variants, the Chinese government has shifted away from the ‘dynamic zero-COVID’ policy, which used to be effective for controlling the spread of the COVID-19 pandemic but with a high economic cost. This study proposes a Dynamic Self-Protection Intensity (DSPI) strategy, a liberalization policy can be used in the later stages of a pandemic. We investigate its impact on society particularly on the burden of medical needs. If citizens in Shenzhen have adopted the weak and strong DSPI strategies proposed in this paper at the beginning of the waiver of ‘dynamic zero-COVID’, it would effectively alleviate the run on medical resources by postponing the peak backwards to two or three lower peaks. This study also provides some insight into how to prevent and control future pandemic outbreaks.

KEYWORDS:

• COVID-19 pandemic
• dynamic self-protection intensity
• medical resource
• sensitivity analysis

1. Introduction

COVID-19 was reported in Wuhan, China in December 2019. The World Health Organization (WHO) announced that COVID-19 became a pandemic in March 2020. Based on established scientific knowledge and international consensus on the characteristics of the novel coronavirus and its associated disease, it was known for its large size and impressive growth rate, which forced many governments to adopt strong control measures such as the national lockdown and social distancing measures [1–3]. Tang et al. [4] studied how well different regions managed to control the spread of SARS-CoV-2 over the whole world. They revealed that East and Southeast Asian populations tended to be more collectivist and self-sacrificing, responding quickly to early signs of the pandemic and readily complying with most restrictions to control its spread. Hensel et al. [5] studied citizens' and governments' responses to the COVID-19 pandemic. They revealed that strong government reactions could correct misperceptions and reduce worries and depression. Uddin et al. [6] studied the seven most commonly used measures, varied notably by concerning their actions in relation to infection rate, disease rate, and timing of measures. The effectiveness of the public measures also depends on whether they are conducted at the right time and according to the right scheme. This is a decision-learning problem, attracting a lot of researchers' attention, see [4,7,8].

In the context of the COVID-19 pandemic, a decisive and swift response to the outbreak guided by scientific evidence may have lowered the number of infections considerably, thereby saving thousands of lives [9,10]. Kusumasari et al. [11] further discussed on how different governments played pivotal roles in responding to the pandemic with evidences, and concluded that they will continue to be crucial both to recovery and to the building of normal activities. To this end, the Chinese government used to be adhered to the policy of ‘dynamic zero-COVID’ by breaking the transmission chain and eliminating the transmission source through extending the scope of close-contact tracing, health-code usage and mass testing, which has effectively controlled the spread of COVID-19 [12].

However, the ‘dynamic zero-COVID’ policy can prevent the spread of the coronavirus with lower infection and mortality rates but may create burden on economy recovery [13]. The study in [14] also finds that stronger and faster anti-pandemic measures are positively associated with trust in the government of handling of the pandemic, but too strict policies usually have a negative impact on the economy and financial markets, resulting in significant income reductions, an increase in unemployment, and disruptions in the transportation, service, and manufacturing industries, among many other consequences [15]. China has decided to abandon the ‘dynamic zero-COVID’ policy after a three-year period of adherence. The center for disease control and prevention reported the numbers of hospitalized patients described as being in severe across the country peaked at the beginning of 2023. Stepping out the wave of the epidemic does not mean that the COVID-19 is over, since the possible emergence of new variants of the virus occurs among the population. It reached its second peak a month later in the private sector clinics. This paper aims to study the impact of people's autonomous awareness of protection on the development of the epidemic.

Due to the constant mutation of the virus, even those who have been infected will inevitably be re-infected by other strains. Due to the ongoing lower risk of novel coronavirus importation and the superior transmissibility of the Omicron variant, the coexistence of COVID-19 coronavirus and human beings is a common cognition of human beings. The most urgent and important thing people need to build, in the long term perspective of liberalization policy, is people's self-protection consciousness. In order to study the intuitive grasp, it is of great significance to study the spread of COVID-19 coronavirus in China. The accurate prediction of the COVID-19 pandemic is foundational to guiding public health policy-making and alleviating socio-economic consequences [16]. Modeling and forecasting the spread of epidemics remain a challenge, depending on social distancing, shelter-in-place orders, disease surveillance, contact tracing, isolation, and quarantine [17]. How to react after a novel infection has a big social and economic impact [18].

This paper will provide new insights into how the pandemic might be mitigated and why the intensity of interventions may need to be changed dynamically. We explore a better understanding of the interaction between dynamic self-protection strategy and the epidemic. Besides, three important factors (age, vaccination and drug) are considered in the modeling of the COVID-19 pandemic. These factors have been extensively discussed in numerous studies; for example, [19,20] discovered that the COVID-19 pandemic caused significantly higher morbidity and mortality among older adults than among the younger people. Xu et al. [21] and Douglas et al. [22] studied the effects of vaccinations and drug. Those are also the main factors when we build our dynamic self-protection intensity model. Another important part of modeling of the COVID-19 pandemic is to estimate the time-varying, instantaneous reproduction numbers, which quantify the viral transmission in real time. They are often defined through a mathematical compartmental model of the epidemic, like a stochastic SEIR model, whose parameters must be estimated from multiple time series of epidemiological data [23]. A general Bayesian method was proposed to study the effects of non-pharmaceutical interventions (NPIs) in reducing COVID-19 transmission in 11 European countries [24]. The lag of infections and the missing data generation mechanism of COVID-19 reported cases were studied in [25] where a regression model is constructed to learn the time-varying reproduction numbers. The prediction of the COVID-19 outbreak in China has attracted the attention of many scholars, such as [26–28].

The main purpose of this paper is to propose an adaptive epidemic self-protection grasp for the post-pandemic era in China. Our strategy aims to provide an approach for ordinary people to contend with the threat posed by the pandemic under the current state of ‘relaxation’ policy and anticipate its consequences in the months ahead. While this paper focuses on research conducted in Shenzhen, the approach can be applied to any other region in China.

As a key city in China's economic development, Shenzhen is a young and dynamic city that is also a super-large port city close to Hong Kong, making it particularly vulnerable to the pandemic. By examining the role of age structure in pandemics and improving real-time targeted forecasting of hospitalization and intensive care needs through dynamic self-protection awareness, we studied a multi-level heterogeneous SEIR compartment model through simulations under different scenarios of self-protection intensity, vaccines, and drug assignment.

In summary, the main contribution of this paper is threefold:

• This paper introduces a new strategy called ‘dynamic self-protection intensity’ (DSPI), corresponding to 5 Levels of self-protection strategies, which enables people to dynamically protect themselves from COVID-19 based on the existing medical treatment resources in Shenzhen;

• Dynamically combining the DSPI and the drug usage to balance the infected cases and the medical treatment resources in Shenzhen;

• Propose a principle of priority use of drugs in the elderly to significantly reduce the risk of exhausting medical resources.

Our simulations indicate that people in countries with aging populations must take aggressive self-protective measures to protect themselves. The specific mechanisms that have been researched will enable an understanding of how the DSPI strategy significantly mitigates the burden on the country and the risks faced by people in the event of a nationwide epidemic. Moreover, the statistical method employed in the Shenzhen case study can be extended to other cities or the entire country.

2. Methods

2.1. Definition of DSPI

During the initial phase of transitioning from a ‘dynamic zero-COVID’ policy, particularly following the first wave of the initial epidemic peak, it becomes imperative to adopt a certain level of self-protection strategies characterized by a reduced crowd contact rate, with the objective of curtailing the speed of epidemic transmission, alleviating the burden on the medical service system arising from the rapid surge in the number of patients, and avert the saturation of medical service resources and capacity. This intensified self-protection awareness is referred to as the spontaneous protective behavior against the epidemic, we call it dynamic self-protection intensity (DSPI) in this paper.

We will define five levels of DSPI and then design a strategy by taking those different measures of DSPI dynamically to balance the number of infections and the medical resources. Specifically,

• Level 1: Masks required in public.

• Level 2: In addition to Level 1, less often go to exhibits, events and large-scale meetings.

• Level 3: In addition to Level 2, less often visit entertainment places and indoor leisure places.

• Level 4: In addition to Level 3, less often go to supermarkets, less dine-in.

• Level 5: In addition to Level 4, less use public transportation, work at home if necessary.

The difference between the ‘dynamic zero-COVID’ policy and the DSPI strategy is that the ‘dynamic zero-COVID’ policy is a government-led epidemic prevention strategy that aims to avoid having confirmed cases in society by dynamically implementing different levels of epidemic prevention policies based on reported confirmed cases. In contrast, the DSPI strategy is a people-led epidemic prevention strategy that is based on citizens' self-protection measures. It dynamically implements DSPI based on the existing medical resources, with the number of hospital bed occupancies for COVID-19 treatment serving as the most important index of medical resources. Unlike the ‘dynamic zero-COVID’ policy, the DSPI is not intended to reduce the number of real-time infections, but rather to reduce the real-time demand for medical resources and prevent a medical run when people get sick. The DSPI allows for a lower degree of transmission of COVID-19 in society, as long as the burden on the healthcare system does not exceed its capacity.

2.2. DSPI strategy

The study conducted an analysis of the medical treatment resources available in Shenzhen for COVID-19 patients. The total number of hospitalizations in general wards was estimated to be 5,069, with 3,839 in two designated hospitals, 150 patients in five professional hospitals, 1,080 patients in 54 fever departments, and approximately 655 in the comprehensive ICU. It should be noted that the number of general wards considered in this study only pertains to those for COVID-19 treatment, excluding treatment beds for other diseases. The total population of Shenzhen, including the floating population, was set to be 21.82 million, which was used as the estimator of the actual administered population in Shenzhen. The proposed DSPI strategy was found to be insensitive to this parameter, and the results can be easily adapted to an adjusted value of the total population in Shenzhen.

The dynamic self-protection intensity (DSPI) strategy can alleviate the pressure on hospitalization treatment by allowing people to adjust their self-protection intensity based on the number of inpatients in hospital wards. We propose two DSPI strategies, namely Weak DSPI and Strong DSPI. In Weak DSPI, the level of self-protection intensity increases as hospital admissions reach 10%, 30%, 40%, and 50% of the capacity, corresponding to DSPI Levels 1–5, respectively. For instance, when the hospital admissions are under 10%, people take Level 1 DSPI. In contrast, in Strong DSPI, the level of DSPI will increase as hospital admissions reach only 5%, 10%, 15%, and 20%. Accordingly, people lower their DSPI level if the hospital admission rate is reduced. The dynamic scheme for the epidemic autonomous awareness based on the hospital beds occupancy for COVID-19 treatment is presented in Table 1, as shown in Figure 1.

Figure 1. Illustration of the weak DSPI strategy.

Table 1. Dynamic self-protection intensity based on the rate of hospital beds occupancy for COVID-19 treatment.

DSPI level	Weak DSPI	Strong DSPI
Level 1	0%–10%	0%–5%
Level 2	10%–30%	5%–10%
Level 3	30%–40%	10%–15%
Level 4	40%–50%	15%–20%
Level 5	50%+	20%+

3. Results

3.1. Prediction of COVID-19 infections

This subsection presents simulation results on COVID-19 infections under different strategies of DSPI: No DSPI, Weak DSPI, and Strong DSPI. The city of Shenzhen, with an estimated population of 21.82 million, is taken as an example. We consider a scenario where the number of infections evolves over a period of 180 days, based on the following parameters: 10 initial imported infections, 40% drug coverage, 80% drug effectiveness, and an optimistic vaccine escape rate. To model the dynamics of the epidemic, a heterogeneous multi-level SEIR compartment model is employed. Details on the model fitting and parameter specification can be found in the Supplementary Materials.

The DSPI strategy could lead to a multi-modal distribution of infections, with multiple spikes occurring in the periods after DSPI is adjusted according to the dynamic autonomous awareness. Figure 2 illustrates the number of total confirmed cases, the number of current confirmed cases, and the number of deaths under three different strategies. If only Level 1 DSPI is taken, almost the entire population of Shenzhen will be infected within 180 days. However, with the Weak DSPI and the Strong DSPI, around 17 million and 12 million cumulative confirmed cases would occur during the same period, respectively. Specifically:

1. In the absence of any DSPI, which means that people's self-protection intensity will always remain at Level 1, the peak of existing infections would be reached around the 68th day, with approximately 5.6 million people infected, and the final death toll during this period would be around 1,950 people.

2. When the Weak DSPI strategy was adopted, the first peak of current diagnoses would be brought forward to about 4.15 million on the 60th day, and the second peak would be about 800,000 on the 128th day. The total number of deaths during this time period would be around 1,650.

3. When the Strong DSPI strategy was adopted, the first peak of existing confirmed cases would reach about 2.1 million on the 57th day, the second peak would reach about 1.2 million on the 108th day, and the third peak would be approximately 800,000 on the 167th day. The total number of deaths during this time period would be around 1,100.

Figure 2. Number of total confirmed cases, current confirmed cases and total deaths in the coming 180 days under three types of DSPIs. The blue, brown, and yellow lines represent no DSPI, weak DSPI, strong DSPI respectively. (a) Total confirmed cases (b) Current confirmed cases (c) Total deaths.

Note that the increase in current confirmed cases under strong DSPI and weak DSPI seems to be faster than the situation with no DSPI. This is possibly because when the DSPI level is at Level 4 or 5, people will increase their testing frequency, allowing infected people to be confirmed more promptly. For a detailed discussion, please refer to the Appendix.

3.2. Prediction of medical resource needs

The core concept of DSPI involves the dynamic adaptation of self-protection intensity by individuals, in response to hospital ward occupancy rates exceeding a certain threshold. This dynamic strategy enables the prevention of medical resource depletion and overreaction, balancing social and economic impacts. Figure 3 demonstrates that:

1. The peak number of general ward and ICU inpatients, using only Level 1 DSPI strategy, will occur on the 71st day with 7,700 and 1,320 patients, respectively.

2. When individuals adopt the weak DSPI strategy, the first peak of general ward inpatients will occur on the 63rd day with 5,500 inpatients, followed by a second peak of 1,050 inpatients on the 130th day. Similarly, the first peak of ICU inpatients will be reached on the 65th day with 910 patients, with a second peak of 190 inpatients on the 128th day.

3. When individuals adopt the strong DSPI strategy, the first peak of general ward inpatients will occur on the 60th day with 2,900 inpatients, comprising 53% of total general inpatient wards. The second peak will occur on the 113rd day with approximately 1,700 inpatients, followed by a third peak of about 1,050 inpatients on the 168th day. The first peak of ICU inpatients will be reached on the 60th day with 470 patients, comprising 72% of total ICU patients. The second and third peaks of ICU inpatients will occur on the 112th day with 300 patients and the 168th day with approximately 190 patients, respectively.

Figure 3. Demands of general wards and ICUs in the coming 180 days under three types of DSPIs. The blue, brown and yellow lines represent no DSPI, weak DSPI, strong DSPI respectively. The horizontal dotted line is the maximum value of general wards and ICUs. (a) Demands of general wards (b) Demands of ICUs.

Figure 4 shows the predicted needs for hospitalization, ICU care, and deaths for different age groups and vaccination status over the next 180 days under three types of DSPI strategies: no DSPI, weak DSPI, and strong DSPI. The results show that as self-protection intensity of people increases, the number of hospitalizations, ICU care, and deaths decrease for all age groups and vaccination status. However, the proportion of hospitalizations, ICU care, and deaths in different age groups and vaccination status remains similar to each other across the three types of DSPI interventions. This suggests that vaccination status is an important factor in predicting the severity of COVID-19 outcomes, regardless of the level of self-protection intensity.

Figure 4. Number of general wards, ICUs, deaths in the coming 180 days for different age groups and vaccination groups under three types of DSPIs. The blue, green, orange, and brown bar graphs represent inpatients who are unvaccinated, first vaccinated, second vaccinated, and booster vaccinated respectively. (a) Demands of general wards under no DSPI (b) Demands of ICUs under no DSPI (c) Total deaths under no DSPI (d) Demands of general wards under weak DSPI (e) Demands of ICUs under weak DSPI (f) Total deaths under weak DSPI (g) Demands of general wards under strong DSPI (h) Demands of ICUs under strong DSPI (i) Total deaths under strong DSPI.

Two significant findings can be extracted from the analysis. Firstly, vaccines have played a pivotal role in preventing severe illnesses and fatalities caused by COVID-19. Booster immunization is effective in combating the virus, and the administration of three doses of the vaccine result in a substantial reduction in severe inpatients. The proportion of intensive care unit (ICU) patients who have received at least two vaccine doses is 9.6%, 8.9%, and 8.8% under no DSPI, weak DSPI, and strong DSPI, respectively. Furthermore, the proportion of people who complete booster immunization is 3.2%, 3.3%, and 2.7% in general wards, and 0.050%, 0.031%, and 0.087% in ICU when self-protection intensity varies from lax to strict. Note that the data for both general wards and ICUs exhibit fluctuations in line with the variation in the DSPI from lax to strict. This is due to the fact that the primary function of booster shots is to prevent fatalities, not infections. The shift from lax to strict in people's self-protection awareness, mainly attributed to a reduction in contact rates, leads to a decrease in infection rates. This indirectly results in the protective effects of vaccines against disease transmission becoming less apparent.

Secondly, the number of severe illnesses and deaths due to COVID-19 increases with age, with the highest risk of death being among those over 60 years old, while the lowest risk is among adolescents aged 0–17. Elderly individuals over 60 years of age accounted for over 75.9%, 76.5%, and 76.3% of the total COVID-19 deaths under no DSPI, weak DSPI and strong DSPI, although the corresponding proportions are 64.16%, 66.57%, and 64.43% in ICU. Inpatients over 60 years old had a case fatality rate of 15.96%, 16.54%, and 15.21% when self-protection intensity varies from loose to strict. The above data are relative numbers reflecting the proportion of deaths among the elderly. The absolute number of deaths decreases for the elderly under stronger DSPI conditions. Therefore, more attention must be given to the elderly population to minimize the loss of life caused by COVID-19.

The demand for medical resources presents multiple peaks when people adopt weak or strong DSPI strategies. Currently, if no DSPI or weak DSPI measures are implemented, there will be a shortage of medical resources. On the other hand, if strong DSPI measures are implemented, the number of infected people and medical resources can be adapted, indicating that the proposed DSPI measures will not overwhelm the medical system as long as elderly immunization coverage is improved to a sufficiently high level.

3.3. Sensitivity analysis

It is important to note that the results presented above are based on specific parameters, including 10 initial imported infections, 40% drug coverage, 80% drug effectiveness, and an optimistic vaccine escape rate. To ensure the robustness of the results, a sensitivity analysis is necessary to examine how the results change when the parameters vary. Certain parameters may have a significant impact on the model and therefore may strongly influence the final results. In this section, we describe and analyze influential parameters under three dynamic DSPI strategies. The sensitivity analysis of parameters forms the foundation for a monitoring or control procedure for process variability.

In this paper, the focus is on six parameters: the number of initially imported cases, transmission rate, drug coverage, drug effectiveness, drug allocation, and levels of immune escape. Among these parameters, drug coverage and drug allocation have a significant impact on the demand for medical resources. Therefore, the sensitivity analysis results of these two parameters are presented in this section, while the results of the other parameters can be found in the Appendix.

3.3.1. Sensitivity analysis for the drug coverage

The impact of varying drug coverage on the demand for medical resources is shown in Figure 5. If there is no drug available and no DSPI is implemented, the peak number of inpatients in general wards and ICUs will occur around the 70th day with 11,800 and 2,000 people respectively, and the total number of deaths will be approximately 2,900. When no DSPI is implemented, only if the drug coverage rate exceeds 80%, can the general wards match the needs of existing medical resources in Shenzhen. However, even then, the ICUs will not be sufficient unless the drug coverage rate approaches 100% if the drug effectiveness parameter is set at 80%. Increasing the drug coverage rate generally slows the rate of epidemic scale growth.

Figure 5. Sensitivity experiments of drug coverage by medical demand (general wards, ICUs, deaths). The blue, orange, green, red, purple and brown lines represent 0, 0.2, 0.4, 0.6, 0.8, 1 respectively. The horizontal dotted line is the maximum value of general wars and ICUs. (a) Demands of general wards under no DSPI (b) Demands of ICUs under no DSPI (c) Total deaths under no DSPI (d) Demands of general wards under weak DSPI (e) Demands of ICUs under weak DSPI (f) Total deaths under weak DSPI (g) Demands of general wards under strong DSPI (h) Demands of ICUs under strong DSPI (i) Total deaths under strong DSPI.

When people adopt the DSPI strategy, the demand for general wards and ICUs will still have multiple but lower peaks due to dynamic control of the pandemic. In the case of weak DSPI without drugs, the first peak in the number of inpatients in general wards would occur around the 60th day with 6,200 people, and the second peak would occur around the 127th day with 2,600 people. If the drug coverage is greater than 60%, the needs of general wards at the two peaks can be met. The first peak in the number of inpatients in ICUs would occur around the 62nd day with 1,000 people and the second peak would occur around the 126th day with 450 people, even without drugs. In conclusion, in the case of weak DSPI, the needs of general wards can be met in line with medical resources in Shenzhen as long as the drug coverage is more than 60%, while ICUs can be adequately resourced only when the drug coverage is more than 80%.

If people adopt the Strong DSPI strategy, the demand for general wards and ICUs will exhibit three distinct peaks. Specifically, the number of general ward inpatients would reach a peak of 3,000 individuals around the 58th day, followed by peaks of 2,000 and 1,600 individuals around the 105th and 155th days, respectively, in the absence of drug intervention. It is worth noting that all of these peaks fall below the existing medical resources in Shenzhen. Similarly, the number of inpatients in ICUs would exhibit three peaks, all of which would remain below the capacity of medical resources in Shenzhen. In the strong DSPI scenario, the adequacy of medical resources would be ensured, regardless of the drug coverage rate, provided that the DSPI strategy is widely adopted.

3.3.2. Sensitivity analysis for the drug allocation

As can be seen from the age-structured infection figures above, there has been an increase in the number of infections among older age groups. It is noteworthy that inpatients in general wards, ICUs, and deaths have responded to this upward shift with a significant spike in the over 60 age group. This may be attributed to the fact that elderly people are more susceptible to coronavirus due to their weaker immune systems. Additionally, they are more likely to suffer from underlying conditions such as heart disease, lung disease or diabetes, which can further weaken their overall health condition, as described in [28].

To reduce hospitalization demand and prevent medical runs, it is recommended to maintain drug coverage at 40% and prioritize drug treatment distribution to the older patient group. Assuming drug effectiveness is 80% with an optimistic level of vaccine escape ratio, and people take the Weak DSPI with convex continuity of transmission rate (see details in the Appendix), the results are shown in Figure 6. Allocating the drug to the 60+ age group with 100% drug coverage, 80% drug effectiveness, and an optimistic level of vaccine escape ratio under Weak DSPI could reduce demand in general wards by 50%, ICUs by 65%, and deaths by nearly 62.5%. Therefore, we suggest that the drug for the over 60 age group be fully covered, while the coverage rate for the other age groups could be lower (over 37%). This suggestion is that due to the fact that the elderly population in Shenzhen is relatively small,the priority of drugs for the elderly population has only a minor impact on the other groups.

Figure 6. Hospitalization impact on priority of drug assignment for elder. The blue and orange lines represent the same drug coverage as 40% among all people, priority of drug treatment for people over 60 years old, respectively. The green horizontal dotted line represent the maximum value of general wards and ICUs. (a) Demand of genaral wards (b) Demand of ICUs (c) Total deaths.

4. Conclusion and discussion

Based on our study and sensitivity analysis, it is evident that maintaining an adequate level of self-protection awareness for citizens and utilizing drugs and vaccines are critical for effective pandemic control. Among these factors, regulating the use of drugs and vaccines is relatively easier for the government to implement, and enhancing self-protection awareness is a basic requirement for citizens.

The use of Weak DSPI can effectively delay the onset of infection in some individuals, leading to a distinctive double-peak pattern in the curve of confirmed cases. The resulting curve of hospitalizations in both general and ICU wards also exhibits a double-peak structure, with each peak exhibiting a lower demand for medical resources. However, even with a drug coverage rate of 40% across all age groups, a shortage of medical resources may still occur during the initial wave of the epidemic.

The adoption of the Strong DSPI results in a three-peak structure for the curves of confirmed cases, hospitalizations in general wards, and ICU admissions. Notably, the peaks of confirmed cases are lower when the DSPI is at the strong level, and the demand for medical resources during these peaks is less than 70% of the total medical resources when drug coverage for all ages is 40%. Our study underscores the crucial role of drugs in pandemic control, but their efficacy is contingent on achieving certain thresholds, namely, drug coverage of at least 40% and drug effectiveness of at least 80%.

Furthermore, our analysis suggests that the allocation of the drug can have a significant impact on pandemic control. Based on our analysis of the composition of hospitalized and deceased patients, we found that individuals aged 60 and above accounted for more than 64% and 75% of the ICU needs and deaths, respectively. Therefore, if the total amount of drugs remains unchanged, prioritizing drug allocation to the over 60 age group can result in a 50% reduction in the number of hospitalizations in general wards, as well as a 65% reduction in the number of hospitalizations and deaths in the ICU.

In conclusion, the adoption of DSPI by the people in Shenzhen can effectively control the pandemic and meet the medical needs of critically ill patients. This approach can be extended to other regions in China with corresponding medical resources and can potentially be implemented on a larger scale across the country. The heterogeneous multi-level SEIR compartment model used in this study can be further improved to simulate pandemic control in different regions and evaluate the effectiveness of various government-led or people-led anti-epidemic programs.

However, it is worth noting that while this study provides valuable insights into the potential effectiveness of a DSPI approach in mitigating the impact of the COVID-19 pandemic on hospitalization demand, it is limited in its focus on only one important factor: the occupancy rate of general wards which may not fully capture the complexity of the healthcare system and resource allocation. As such, it is important to acknowledge that there may be other critical variables to consider in decision-making around pandemic control, such as healthcare worker resources, testing resources, and drug resources, among others. It is therefore essential to approach pandemic control measures with a multidisciplinary perspective and remain adaptable to changing circumstances and emerging data, while continually striving to improve and optimize response strategies according to the development of the pandemic. This study also provides some insight into how to prevent and control future pandemic outbreaks.

Acknowledgments

Author Contributions Conceptualization, Chao Liu, Zhen Zhang and Jian Qing Shi; Data curation, Xinxin Han; Formal analysis, Bin Zhu; Funding acquisition, Dongfeng Gu; Methodology, Chao Liu, Guanpeng Li and Jian Qing Shi; Project administration, Jian Qing Shi; Resources, Jie Huang; Software, Haonan Long and Pengzhen Chen; Supervision, Zhen Zhang and Yanqing Hu; Validation, Yanqing Hu; Writing – original draft, Chao Liu and Guanpeng Li; Writing – review & editing, Jie Huang.

Availability of Data and Materials

All data generated or analyzed during this study are included in this published article and its Appendix. All data used in the study is publicly available.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

Chao Liu was supported by the Shenzhen University Research Start-up Fund for Young Teachers (No. 868-000001032037); Jian Qing Shi was supported by Shenzhen Fundamental Research Program (No. JCYJ20220818100602005); Dongfeng Gu was supported by the Shenzhen Science and Technology Innovation Committee (No. JSGG20220301090202005).

Notes on contributors

Chao Liu

Chao Liu is an assistant professor in the College of Economics, Shenzhen University, Postal address: No. 1066, Xueyuan Avenue, Nanshan District, Shenzhen, Guangdong Province, China.

Haonan Long

Haonan Long is undergraduate student in the department of statistics and data science, Southern University of Science and Technology. Postal address: No. 1088, Xueyuan Avenue, Nanshan District, Shenzhen, Guangdong Province, China.

Guanpeng Li

Guanpeng Li is a PhD student in the School of Public Health and Emergency Management, Southern University of Science and Technology. Postal address: No. 1088, Xueyuan Avenue, Nanshan District, Shenzhen, Guangdong Province, China.

Pengzhen Chen

Pengzhen Chen is undergraduate student in the department of statistics and data science, Southern University of Science and Technology. Postal address: No. 1088, Xueyuan Avenue, Nanshan District, Shenzhen, Guangdong Province, China.

Zhen Zhang

Zhen Zhang is a professor in the department of mathematics, Southern University of Science and Technology. Postal address: No. 1088, Xueyuan Avenue, Nanshan District, Shenzhen, Guangdong Province, China.

Zhen Zhang belong to the National Center for Applied Mathematics in Shenzhen, Shenzhen, China.

Jie Huang

Jie Huang is professor in the School of Public Health and Emergency Management, Southern University of Science and Technology. Postal address: No. 1088, Xueyuan Avenue, Nanshan District, Shenzhen, Guangdong Province, China.

Bin Zhu

Bin Zhu is professor in the School of Public Health and Emergency Management, Southern University of Science and Technology. Postal address: No. 1088, Xueyuan Avenue, Nanshan District, Shenzhen, Guangdong Province, China.

Xinxin Han

Xinxin Han is professor in the School of Public Health and Emergency Management, Southern University of Science and Technology. Postal address: No. 1088, Xueyuan Avenue, Nanshan District, Shenzhen, Guangdong Province, China.

Yanqing Hu

Yanqing Hu is professor in the department of statistics and data science, Southern University of Science and Technology. Postal address: No. 1088, Xueyuan Avenue, Nanshan District, Shenzhen, Guangdong Province, China.

Jian Qing Shi

Jian Qing Shi is professor in the department of statistics and data science, Southern University of Science and Technology. Postal address: No. 1088, Xueyuan Avenue, Nanshan District, Shenzhen, Guangdong Province, China.

Jian Qing Shi belong to the National Center for Applied Mathematics in Shenzhen, Shenzhen, China.

Dongfeng Gu

Dongfeng Gu is professor in the School of Public Health and Emergency Management, Southern University of Science and Technology. Postal address: No. 1088, Xueyuan Avenue, Nanshan District, Shenzhen, Guangdong Province, China.

Appendices

Appendix 1. Mathematical Materials

A.1. SEIQRD model

To characterize the epidemic of COVID-19, we generalize the classical SEIR model by introducing six differentstates: S, E, I, Q, R, D, denoting the number of susceptible cases, exposed cases (infected but not yet infectious, in a latent period), infectious cases (with infectious capacity but not yet quarantined), quarantined cases (confirmed and infected), recovered cases and closed cases (or deaths), respectively. The susceptible population S was divided into four stages based on the vaccine level. $S0$ corresponds to the susceptible population that has not been vaccinated. $S1$ corresponds to the susceptible population that has received only one dose of the vaccine. $S2$ corresponds to the susceptible population that has completed basic immunization. $S3$ corresponds to the susceptible population that has completed enhanced immunization. Note that states Q and D were isolated, so they were not transmissible, while the remaining four states were free to move and spread socially; see the flow chart in Figure A1.

Figure A1. SEIQRD model.

In the SEIQRD model, we would not implement contact tracing of people infected with COVID-19. That is, we only quarantine those who have been diagnosed with COVID-19; people in states S and E would not be quarantined. This is different from the early work on the prevention and control of COVID-19 in China; see [13,27]. In this case, the SEIQRD model can be described as follows: $\begin{array}{rl}& \frac{\mathrm{d}{S}_i}{\mathrm{d}t}=-\frac{\text{β}(1-{\nu }_i)I{S}_i}{N},i=0,1,2,3\\ & \frac{\mathrm{d}E}{\mathrm{d}t}=\frac{\sum \text{β}I(1-{\nu }_i)S_i}{N}-\mathrm{\gamma E}\\ & \frac{\mathrm{d}I}{\mathrm{d}t}=\mathrm{\gamma E}-\mathrm{\lambda I}\\ & \frac{\mathrm{d}Q}{\mathrm{d}t}=\mathrm{\lambda I}-\mathrm{\delta Q}-\mathrm{kQ}\\ & \frac{\mathrm{d}R}{\mathrm{d}t}=\mathrm{\delta Q}\\ & \frac{\mathrm{d}D}{\mathrm{d}t}=\mathrm{kQ}\\ & N=\sum _{i=0}^3{S}_i+E+I+R\end{array}$ The target is to estimate the transmission parameters of COVID-19 without any dynamic self-protection awareness(DSPI) among individuals. For model fitting, data from March 1 to March 27, 2022 were chosen. It was noted that Shanghai did not adopt the national nucleic acid policy from March 1 to March 27. The lack of nucleic acid testing in this period results in a difference between the reported confirmed cases and the actual infected cases, which will be one of the major challenges for future research.

This paper will fit the model parameters based on the reported confirmed cases. To reduce the fitting complexity, we fix some parameters in the SEIQRD model; the meaning and value of these parameters can be seen in Table A2. The transition rate from the exposed group (E) to the infected group (I) is chosen to be $\gamma =0.2777$, which is derived from the study in the literature [29], which showed that the average incubation period of the Omicron variant was 3.6 days. The transition rate from state Q to state R is chosen as $\delta =0.1195$, which is derived from the hospitalization probability, the recovery rate of general wards, the hospital conversion rate, the recovery rate of ICUs, and the average length of hospital stay for each period, referring to [28]. The transition rate from state Q to state D is estimated at k = 0.0011, the value is derived from the mortality rate of general wards and ICUs in the reference [1], the mortality data for the fifth wave of the Hong Kong epidemic in March 2022 validated the rationality of this parameter. According to the actual conditions in China, the vast majority of the vaccines vaccinated by Chinese people are inactivated vaccines made in China. Therefore, this article does not consider other mRNA or hybrid immunity for the time being. The vaccine protective parameters of the inactivated vaccine are chosen to be the same as [28], and the values of the three vaccine effectiveness parameters are shown in Table A2, the vaccination data for Shanghai and Shenzhen up to October, 2022, are shown in Table A1.

Table C1. Vaccination status as of October 28, 2022.

Location	Total population	Unvaccinated	First dose	Second dose	Booster
Shanghai	24871100	618900	24252200	23093900	11932100
Shenzhen	21819000	1146500	20672500	19651000	13309900

Table C2. Parameters used in SEIQRD model.

Parameters	Definition	Value or distribution
β	Transition rate from susceptible individuals to exposed individuals	$N(1.0342,0.3261)$
γ	Transition rate from exposed individuals to infected individuals	0.2777
λ	Transition rate from infected individuals to confirmed individuals	$N(0.2985,0.1600)$
δ	Recovery rate of the confirmed individuals	0.1195
κ	mortality of the confirmed individuals	0.0011
${\nu }_1$	Rate of protection given only one dose of vaccine against infection in susceptible class	5.6%
${\nu }_2$	Rate of protection given two doses against infection in susceptible class	9.1%
${\nu }_3$	Rate of protection given booster immunization against infection in susceptible class	17%

The final two parameters, β and λ, must be estimated using reported and confirmed data. We will give the detailed procedure for the estimation of these two parameters. We aimed to provide a distribution of these two parameters rather than estimating a fixed value of these two parameters, such as [13,28]. The rolling sampling method was an appropriate statistical method to estimate the probability distribution of unknown parameters β and λ. The sampled data from March 1 to March 27 is a continuous time point with length $h=(10,27;1)$. The notation $(a,b;c)$ represents a vector with values ranging from a to b with interval c. For example, when h = 10, a total of 17 sets of data could be obtained by rolling sampling from March 1 to March 27, 2022, in which the first set of data is from March 1 to 10 with length h = 10. From this viewpoint, 132 sampling groups were obtained by rolling sampling to fit the parameters, so as to learn the probability distributions of the two parameters. The 132 estimators for the parameters were obtained and presented in Figure A2. The normal distribution $N(1.0342,0.3261)$ and $N(0.2985,0.1600)$ can be used to approximate the distribution of these two parameters as shown in Table A1.

Figure A2. The Distribution of estimation of parameter β and λ. (a) The Distribution of β (b) The Distribution of λ.

Figure A3 shows the scatter diagram of the 132 sampled estimations of parameters β and λ, which represents a linear regression model as $\lambda =-0.1435+0.3795\cdot \text{β}$, $R^2=0.7783$ with the P-value smaller than 0.05. Therefore, there is a significant linear positive correlation between the probability of successful transmission from exposed individuals to infected individuals and the transition rate of infected individuals to confirmed individuals; that is, as the probability of successful transmission increases, the transition rate of infected individuals to confirmed individuals also increases. It is tempting to conclude that the relationship is one based on cause and effect: by increasing the probability of successful transmission, we could always ensure a high transition rate of infected individuals to confirmed cases. This result is consistent with the real world situation; that is, when the epidemic spreads faster, people's awareness of self-protection about the epidemic will be enhanced which will make the people increase their frequency of nucleic acid detection accordingly, so the transition rate from infected individuals to confirmed cases will also increases.

Figure A3. Scatter plot of parameters β and λ for the 132 sampled estimations.

A.2. SEIR-HML model

In order to estimate the demand for medical resources under the impact of the epidemic in detail, we divided the population by age group and vaccination status. The more detailed values and parameters are provided in Table A3. Specifically, the SEIR-HML model divides the population into six age groups, four vaccination levels, and three different DSPI ideology, including no DSPI, weak DSPI, and strong DSPI. The flow chart of the proposed SEIR-HML model is shown in Figure A4.

Figure A4. Illustration of the SEIR-HML model.

Table C3. Notations in SEIR-HML model.

Parameters	Description
$V\left(t\right)$	Number of susceptible populations at time t
$E\left(t\right)$	Number of latent populations at time t
$I\left(t\right)$	Number of infected populations at time t
$Q\left(t\right)$	Number of asymptomatic populations at time t
$Q_A\left(t\right)$	Number of asymptomatic populations at time t
$Q_S^{\mathrm{preO}}\left(t\right)$	Number of patients use treatment at time t
$Q_S^{\mathrm{preR}}\left(t\right)$	Number of symptomatic populations will recover by treatment at time t
$Q\left(t\right)$	Number of asymptomatic populations at time t
$U^{\mathrm{preD}}\left(t\right)$	Number of patients in ICU die at time t
$U^{\mathrm{preR}}\left(t\right)$	Number of patients in ICU recover at time t
$H^{\mathrm{preD}}\left(t\right)$	Number of patients in general ward die at time t
$H^{\mathrm{preR}}\left(t\right)$	Number of patients in general ward recover at time t
$D\left(t\right)$	Total number of deaths at time t
$R\left(t\right)$	Total number of people recovered at time t
$T_{11}$	Number of people switching from susceptible group to latent group at time t
$T_{21}$	Number of people switching from latent group to infection group at time t
$T_{31}$	Number of infected people to confirmed people at time t
$T_{41}$	Number of infected people to confirmed symptomatic people at time t
$T_I^{\mathrm{Qs}}$	Number of infected people to confirmed symptomatic people at time t
$T_S^{\mathrm{Drug}}$	Number of confirmed symptomatic people to medicated group at time t
$T_S^{\mathrm{nonDrug}}$	Number of confirmed symptomatic people to non-drug group at time t
$T_{42}$	Number of infected people to confirmed asymptomatic people at time t
$T_s^{\mathrm{Drug}}$	Number of confirmed symptomatic people to medicated group at time t
$T_{51}$	Number of infected people to confirmed asymptomatic people at time t
$T_{52}$	Number of confirmed symptomatic people to non-drug group will be admitted to hospital at time t
$T_{53}$	Number of symptomatic people to be hospitalized population at time t
$T_{61}$	Number of confirmed symptomatic people to medicated group recover at time t
$T_{62}$	Number of confirmed symptomatic people to non-drug group recover at time t
$T_{63}$	Number of symptomatic people to recovered group at time t
$T_Q^H$	Number of symptomatic people to general ward at time t
$T_Q^U$	Number of symptomatic people to ICU at time t
$T_{71}$	Number of ICU patients died at time t
$T_{72}$	Number of $U_{a,t}^{\mathrm{preD}}$ to death
$T_{81}$	Number of ICU patients recovered at time t
$T_{82}$	Number of $U_{a,t}^{\mathrm{preR}}$ to recovered group
$T_{91}$	Number of $T_Q^H$ to death
$T_{92}$	Number of general patients died at time t
$T_{x1}$	$T_Q^H-T_{91}$
$T_{x2}$	Number of general patients recovered at time t

The SEIR-HML model is a generalization of the above SEIQRD model, which models the transmission of COVID-19 for each age group separately in a more detailed version. Figure A4 shows the transmission process for age group a, in which the first three parameters $\text{β},\gamma ,\lambda$ are all the same for people of all ages, while the other parameters will change with age group. In particular, for the age group a, confirmed individuals have a $p_a^s$ chance of being symptomatic, whereas others are asymptomatic. The model assumes that asymptomatic infected individuals would isolate themselves at home without drug treatment. For symptomatic infected individuals, drug treatment would be carried out at the current level. It means that when drug coverage reaches at $p_a^{\mathrm{drug}}$, there would be $p_a^{\mathrm{drug}}$ proportion of confirmed symptomatic patients who would receive the drug treatment. In the absence of drug therapy, each diagnosed person would be hospitalized with the probability of $p_a^h$. Assuming that the effective rate of drug protection from hospitalization is ${ϵ}_{\mathrm{drug}}$, the population receiving drug treatment would be hospitalized with the probability of ${ϵ}_{\mathrm{drug}}\cdot p_a^h$. Every inpatient admitted to the general ward has a probability of $1-p_a^{\mathrm{icu}}$ and every inpatient admitted to the ICU has a probability of $p_a^{\mathrm{icu}}$. In particular, patients in the general wards have a certain probability of being transferred to the ICU. The average time of transfer is assumed to be the Poisson distribution with a mean of 6. The mortality rates in general wards and ICUs are $p_a^{\mathrm{HD}}$ and $p_a^{\mathrm{UD}}$, respectively. The time if took for patients to recovery from general ward and ICUs was assumed to follow a Poisson distribution with mean value $\frac{1}{{\gamma }^{\mathrm{HR}}}$ and $\frac{1}{{\gamma }^{\mathrm{UR}}}$, respectively. The time of patients from general wards and ICUs to death was assumed to be a Poisson distribution with a mean value of $\frac{1}{{\gamma }^{\mathrm{HD}}}$ and $\frac{1}{{\gamma }^{\mathrm{UD}}}$, respectively.

The discrete dynamic equation model (SEIR-HML) for age group a can be specified as follows, where the subscript of age group notation a was omitted for convenience. $\begin{array}{rl}& V\left(t\right)=V(t-1)-T_{11}\\ & E\left(t\right)=E(t-1)+T_{11}-T_{21}\\ & I\left(t\right)=I(t-1)+T_{21}-T_{31}\\ & Q_A\left(t\right)=Q_A(t-1)+T_{41}-T_{42}\\ & Q_S^{\mathrm{preO}}\left(t\right)=Q_S^{\mathrm{preO}}(t-1)+T_{51}+T_{52}-T_{53}\\ & Q_S^{\mathrm{preR}}\left(t\right)=Q_S^{\mathrm{preR}}(t-1)+T_{61}+T_{62}-T_{63}\\ & U^{\mathrm{preD}}\left(t\right)=U^{\mathrm{preD}}(t-1)+T_{71}-T_{72}\\ & U^{\mathrm{preR}}\left(t\right)=U^{\mathrm{preR}}(t-1)+T_{81}-T_{82}\\ & H^{\mathrm{preD}}\left(t\right)=H^{\mathrm{preD}}(t-1)+T_{91}-T_{92}\\ & H^{\mathrm{preR}}\left(t\right)=H^{\mathrm{preR}}(t-1)+T_{x1}-T_{x2}\\ & D\left(t\right)=D(t-1)+T_{72}+T_{92}\\ & R\left(t\right)=R(t-1)+T_{42}+T_{63}+T_{82}+T_{x2}\\ & T_{41}=T_{51}+T_{52}+T_{61}+T_{62}\\ & T_{53}=T_{71}+T_{81}+T_{91}+T_{x1}\\ & Q\left(t\right)=Q_A\left(t\right)+Q_S^{\mathrm{preO}}\left(t\right)+Q_s^{\mathrm{preR}}\left(t\right)\end{array}$ The meanings of the notations above can be seen in Table A3. The random distribution of each conversion volume between compartments (T variables) was shown in the following. Denote $I_{\mathrm{all}}$ is the total number of all age groups in the compartment of the confirmed population. $\begin{array}{rl}& T_{11}\sim \mathrm{Binomial}(V,(1-\nu )\cdot \text{β}\cdot I_{\mathrm{all}})\\ & T_{21}\sim \mathrm{Binomial}(E,\gamma )\\ & T_{31}\sim \mathrm{Binomial}(I,\lambda )\\ & T_{41}=T_{31}-T_I^{Q_S}\\ & T_I^{Q_S}\sim \mathrm{Binomial}(T_{31},p^S(1-V{E}^S\left)\right)\\ & T_S^{\mathrm{Drug}}\sim \mathrm{Binomial}(T_I^{Q_S},p^{\mathrm{drug}})\\ & T_S^{\mathrm{nonDrug}}=T_I^{Q_S}-T_S^{\mathrm{Drug}}\\ & T_{51}\sim \mathrm{Binomial}(T_S^{\mathrm{Drug}},p^h\cdot (1-V{E}^h)\cdot (1-{ϵ}_{\mathrm{drug}}\left)\right)\\ & T_{61}=T_S^{\mathrm{Drug}}-T_{51}\\ & T_{52}\sim \mathrm{Binomial}(T_S^{\mathrm{nonDrug}},p^h\cdot (1-V{E}^h\left)\right)\end{array}$ $\begin{array}{rl}& T_{62}=T_S^{\mathrm{nonDrug}}-T_{52}\\ & T_{53}\sim \mathrm{Binomial}(Q_S^{\mathrm{preO}},1-e^{(-{\gamma }^{\mathrm{SH}})})\\ & T_{63}\sim \mathrm{Binomial}(Q_S^{\mathrm{preR}},1-e^{(-{\gamma }^I)})\\ & T_Q^U\sim \mathrm{Binomial}(T_{53},p^{\mathrm{icu}}\cdot (1-V{E}^{\mathrm{ICU}}\left)\right)\\ & T_Q^H=T_{53}-T_Q^U\\ & T_{42}\sim \mathrm{Binomial}(Q_A,1-e^{(-{\gamma }^I)})\\ & T_{71}\sim \mathrm{Binomial}(T_Q^U,p^{\mathrm{UD}}\cdot (1-V{E}^{\mathrm{death}}\left)\right)\\ & T_{81}=T_Q^U-T_{71}\\ & T_{72}\sim \mathrm{Binomial}(U^{\mathrm{preD}},1-e^{(-{\gamma }^{\mathrm{UD}})})\\ & T_{82}\sim \mathrm{Binomial}(U^{\mathrm{preR}},1-e^{(-{\gamma }^{\mathrm{UR}})})\\ & T_{91}\sim \mathrm{Binomial}(T_Q^H,p^{\mathrm{HD}}\cdot (1-V{E}^{\mathrm{death}}\left)\right)\\ & T_{x1}=T_Q^H-T_{91}\\ & T_{92}\sim \mathrm{Binomial}(H^{\mathrm{preD}},1-e^{(-{\gamma }^{\mathrm{HD}})})\\ & T_{x2}\sim \mathrm{Binomial}(H^{\mathrm{preR}},1-e^{(-{\gamma }^{\mathrm{HR}})})\end{array}$

A.3. Dynamic self-protection intensity

The dynamic self-protection intensity (DSPI) strategy is presented in Figure 1 in the main paper. Given the above mathematical models and parameters, the main focus is on how to determine the values of the parameters for the five levels of the dynamic self-protection measures taken by the people independently. This subsection will give a detailed illustration.

Recall that the parameters obtained in Section A.1 consist of two parts. One is for the fixed parameters, which are assigned values based on previous research [13,28]; the other is for the two important parameters estimated through rolling sampling. That is, the distributions of transmission rates β and λ are obtained, which correspond to the parameters under normal control (i.e. the Level 1 DSPI). The main concern now is how to take values for the other four levels of DSPI for parameters β and λ.

A direct idea is to control the frequency of nucleic acid detection to be ${\lambda }_1={\lambda }_2={\lambda }_3=\lambda$, where $\lambda =0.2985$ is the estimator from the Level 1 DSPI. According to Shenzhen's previous successful epidemic prevention and control strategy, ${\lambda }_4=0.5$ corresponds to nucleic acid testing once every two days, and ${\lambda }_5=1$ corresponds to nucleic acid testing once a day. As for the transmission rate β, note that the main objective of the 5 different levels of DSPI interventions was to control the scale of outbreak spread by controlling social distancing, i.e. the contact rate, so that the transmission rate β for different levels of DSPI depends mainly on the contact rate. According to the literature [26], under normal circumstances, the population contact rate in Shenzhen without any social distancing measures (corresponding to Level 1 SPA) is 7.9, and the population contact rate in lockdown scenarios (corresponding to level 5 DSPI) is 2.2. If the transmission rate of the Level 1 DSPI is ${\text{β}}_1=\text{β}=1.0342$, the transmission rate of the Level 5 DSPI is ${\text{β}}_5=\frac{2.2}{7.9}\text{β}$. The parameter estimates under the DSPI measures of Level 2-Level 4 are obtained in the form of equal difference series, that is, ${\text{β}}_2=0.8196\text{β}$, ${\text{β}}_3=0.6393\text{β}$, ${\text{β}}_4=0.4589\text{β}$, respectively. These parameters are shown in Table A5.

Table C4. Parameters used in SEIR-HML model.

Parameters	Definition	Value or distribution
β	Probability of successful transmission	$N(1.0342,0.3261)$
γ	Transition rate of latent individuals to the infected class	0.2777
λ	Transition rate of infected individuals to confirmed individuals	$N(0.2985,0.1600)$
ν	To simplify the formula, ν is used to replace ${\nu }_0$, ${\nu }_1$, ${\nu }_2$ or ${\nu }_3$, corresponding to the vaccine protection for unvaccinated, first dose, second dose and booster, respectively	see the following four lines
$V{E}^S$	Rate of vaccine ν protection against infection	${\nu }_0^S=0\mathrm{\%}$, ${\nu }_1^S=5.6\mathrm{\%}$, ${\nu }_2^S=9.1\mathrm{\%}$, ${\nu }_3^S=17.0\mathrm{\%}$
$V{E}^h$	Rate of vaccine ν protection against entry into general wards	${\nu }_0^h=0$, ${\nu }_1^h=46.2\mathrm{\%}$, ${\nu }_2^h=78.8\mathrm{\%}$, ${\nu }_3^h=98.1\mathrm{\%}$
$V{E}^{\mathrm{ICU}}$	Rate of vaccine ν protection against entry into ICU	${\nu }_0^{\mathrm{ICU}}=0$, ${\nu }_1^{\mathrm{ICU}}=46.2\mathrm{\%}$, ${\nu }_2^{\mathrm{ICU}}=78.8\mathrm{\%}$, ${\nu }_3^{\mathrm{ICU}}=98.1\mathrm{\%}$
$V{E}^{\mathrm{death}}$	Rate of vaccine ν protection against death	${\nu }_0^{\mathrm{death}}=0$, ${\nu }_1^{\mathrm{death}}=56.3\mathrm{\%}$, ${\nu }_2^{\mathrm{death}}=83.2\mathrm{\%}$, ${\nu }_3^{\mathrm{death}}=98.4\mathrm{\%}$
$p^s$	Probability of unvaccinated people have symptoms	Age 0–17: 7.16%; Age 18–39: 8.87%; Age 40–59: 12.09%; Age 60+: 15.89%
$p^{\mathrm{drug}}$	Drug coverage among patients	0/0.2/0.4/0.6/0.8/1
$p^h$	Probability of hospitalization	Age 0–2: 11.09%; Age 3–11: 6.37%; Age 12–17: 15.01%; Age 18–29: 7.41%; Age 30–39: 9.78%; Age 40–49: 8.02%; Age 50–59: 17.06%; Age 60–69: 20.45%; Age 70+: 65.60%
${ϵ}_{\mathrm{drug}}$	Effectiveness of drug treatment	0.5/0.8
$p^{\mathrm{UD}}$	Probability of ICU inpatients	Age 0–2: 3.40%; Age 3–11: 6.74%; Age 12–17: 9.81%; Age 18–29: 12.48%; Age 30–39: 10.45%; Age 40–49: 20.66%; Age 50–59: 33.32%; Age 60–69: 36.45%; Age 70+: 37.47%
$p^{\mathrm{HD}}$	Probability of general inpatients	Age 0–2: 0.26%; Age 3–11: 0.52%; Age 12–17: 0.75%; Age 18–29: 0.96%; Age 30–39: 0.80%; Age 40–49: 3.23%; Age 50–59: 5.75%; Age 60–69: 8.68%; Age ≥ 70: 24.98%
$1/{\gamma }_I$	Average time from symptoms to recovery	5.64
$1/{\gamma }_{\mathrm{SH}}$	Average time from symptoms to general wards	2.2
$1/{\gamma }_{\mathrm{HD}}$	Average time from general wards to death	8
$1/{\gamma }_{\mathrm{UD}}$	Average time from ICU to death	8
$1/{\gamma }_{\mathrm{HR}}$	Average time from general wards to death	6
$1/{\gamma }_{\mathrm{UR}}$	Average time from ICU to recovery	8

Table C5. Estimates of parameter based on usage proportional of general wards under different DSPI scenario.

SPA Level	Weak DSPI	Strong DSPI	Estimates of parameters
1	0%–10%	0%–5%	$\text{β}={\text{β}}_1,\lambda =0.2985$
2	10%–30%	5%–10%	$\text{β}={\text{β}}_2,\lambda =0.2985$
3	30%–40%	10%–15%	$\text{β}={\text{β}}_3,\lambda =0.2985$
4	40%–50%	15%–20%	$\text{β}={\text{β}}_4,\lambda =0.5$
5	50%+	20%+	$\text{β}={\text{β}}_5,\lambda =1$

Table C6. Vaccination coverage by age in Shenzhen.

Age group	Unvaccinated	First dose	Second dose	Booster
0–17	0.26061	0.26061	0.12404	0
18–29	0.31093	0.31093	0.27903	0.26172
30–39	0.20041	0.20041	0.27961	0.32348
40–49	0.10335	0.10335	0.16635	0.22059
50–59	0.06989	0.06989	0.10978	0.14679
60+	0.05481	0.05481	0.04119	0.04742

The transmission rate values listed above are a reasonable choice for various levels of DSPI. Furthermore, we perform sensitivity analysis for the transmission rate β to vary in a continuous manner from Level 1-Level 5 with the above choice of β being only one special case; see the details in Section A.5, the value of ${\text{β}}_1-{\text{β}}_5$ can refer to Tables A7 and A8 under weak and strong DSPI.

Table C7. The curve $\text{β}\left(t\right)$ for three continuous scenario under weak DSPI $(\text{β}=1.0342)$.

Intensity	$\text{β}\left(t\right)\left(x\right(t)<0.5)$
Convex	$\text{β}\left(t\right)=\text{β}\cdot \left(\right(e^4-0.2785\left)/\right(e^4-1)-0.7215\cdot e^{8\cdot x\left(\mathrm{tc}\right)}/(e^4-1\left)\right)$
Lower Convex	$\text{β}\left(t\right)=\text{β}\cdot \left(\right(0.2785-e^{-4}\left)/\right(1-e^{-4})+0.7215\cdot e^{-8\cdot x\left(t\right)}/(1-e^{-4}\left)\right)$
Linear	$\text{β}\left(t\right)=\text{β}\cdot (1-1.425\cdot x(t\left)\right)$

Table C8. The curve $\text{β}\left(t\right)$ for three continuous scenario under strong DSPI $(\text{β}=1.0342)$.

Intensity	$\text{β}\left(t\right)\left(x\right(t)<0.2)$
Convex	$\text{β}\left(t\right)=\text{β}\cdot \left(\right(e^4-0.2785\left)/\right(e^4-1)-0.7215\cdot e^{20\cdot x\left(t\right)}/(e^4-1\left)\right)$
Lower Convex	$\text{β}\left(t\right)=\text{β}\cdot \left(\right(0.2785-e^{-4}\left)/\right(1-e^{-4})+0.7215\cdot e^{-20\cdot x\left(t\right)}/(1-e^{-4}\left)\right)$
Linear	$\text{β}\left(t\right)=\text{β}\cdot (1-3.6075\cdot x(t\left)\right)$

A.4. Parameters and data summary

This subsection will give an illustration of the notations and parameters used in the proposed SEIR-HML model, and the summary of data used in this paper is also presented, see Table A1–A8 at the end of this manuscript.

Appendix 2. Additional Simulation results

In the main paper, we presented a sensitivity analysis for parameters related to drug coverage and drug allocation. In this section, we will describe and illustrate the sensitivity analysis for several influential parameters under three dynamic DSPI intervention actions. Specifically, we will focus on four parameters: the number of initial imported cases, transmission rate, drug effectiveness, and level of immune escape.

A.5. Sensitivity analysis for transmission rate

A reasonable choice for the transmission rate parameter is described in Section A.3, where the value is discretely taken for the five DSPI levels. In this subsection, we consider a continuous varying value for the transmission rate parameter $\text{β}\left(t\right)$, which is varied with the occupancy of the general wards $x\left(t\right)$. To describe the viewpoints of the government to the economy and epidemic, we divide the interventions into three types: prioritized economic protection, which corresponds to the ‘Loose first, then strict’ scenario; prioritize control of the spread of the epidemic, which corresponds to the ‘Strict first, then loose’ scenario; and epidemic prevention and control policies without priority (theoretical ideal state). The corresponding transmission rate value is shown in Figure A5 for the dynamic peak suppression strategy under weak DSPI scenario. While analyzing the effectiveness of transmission rate, other parameters are fixed at given values with 40$\mathrm{\%}$ drug coverage rate, 80% drug efficiency, 10 initial infusion infections and a low virus vaccine escape rate.

1. Epidemic prevention and control policies that prioritize economic protection (convex continuity) ${\text{β}}^{\left(1\right)}\left(t\right)=\text{β}\cdot (a_1+b_1\cdot e^{m\cdot x\left(t\right)}),$ when $x\left(t\right)\le 0.5$ and ${\text{β}}^{\left(1\right)}\left(t\right)=0.2785\text{β}$ when $x\left(t\right)>0.5$ in the weak DSPI scenario.

2. Epidemic prevention and control policies that prioritize control of the spread of the epidemic (lower convex continuity) ${\text{β}}^{\left(2\right)}\left(t\right)=\text{β}\cdot (a_2+b_2\cdot e^{-m\cdot x\left(t\right)}),$ when $x\left(t\right)\le 0.5$ and ${\text{β}}^{\left(2\right)}\left(t\right)=0.2785\text{β}$ when $x\left(t\right)>0.5$ in the weak DSPI scenario.

3. Epidemic prevention and control policies without priority (linear continuity) ${\text{β}}^{\left(3\right)}\left(t\right)=\text{β}\cdot (a_3+b_3\cdot x(t\left)\right),$ $x\left(t\right)\le 0.5$ and ${\text{β}}^{\left(3\right)}\left(t\right)=0.2785\text{β}$ when $x\left(t\right)>0.5$ in the weak DSPI scenario. The strong DSPI scenario is similar. Here, β is the estimated transmission rate parameter and the other values of $a_1,b_1,a_2,b_2,a_3,b_3$ for the weak and strong DSPI scenarios are shown in Tables A7 and A8, respectively.

Figure A5. Three scenario of transmission rate parameter varied with the occupancy of the general wards under weak DSPI which can achieve dynamic peak suppression.

When individuals adopt the weak DSPI, compared with linear continuity self-protection measure which can successfully limit the demand of the general wards but fails in the demand of ICU, the convex continuity self-protection measure will increase the values of the need for general wards and ICU by 28$\mathrm{\%}$, making medical resources insufficient during the first outbreak peak of the epidemic but delaying the occurrence time of the second peaks by about 18 days with a 20$\mathrm{\%}$ reduction in medical resources. However, the lower convex continuity self-protection measure will decrease the general wards' demands in the first outbreak by 50$\mathrm{\%}$ which only equals 46$\mathrm{\%}$ of the total number of wards, and delay the occurrence of the first outbreak by 10 days. Its second outbreak will occur at the same time as the linear one, but the peak values will be raised by about 45$\mathrm{\%}$. The total deaths would be lower from the convex continuity scenario (‘Loose first, then strict’ scenario) to the linear continuity scenario (without priority) and the lower convex continuity scenario (‘Strict first, then loose’ scenario). The corresponding DSPI scenario would be implemented as the Figure A6 shows.

Figure A6. Demands for general wards and ICUs, deaths and PHSM levels implemented under weak PHSM strategy for three types of transmission rates. (a) Demands of general wards (b) Demands of ICUs (c) Number of Deaths (d) PHSM intervention levels.

In the case of strong DSPI, similar conclusions can be drawn and the results are shown in Figure A7. To conclude, as the transmission rate varies from ‘Loose first, then strict’ policy to ‘Strict first, then loose’ scenario, the demands on medical resources and deaths will be lower and lower, which shows the strict first scenario is more friendly for individuals to reduce the risk of infection of the COVID-19 epidemic.

Figure A7. Demands of general wards and ICUs, deaths and SPA levels implemented under strong DSPI strategy for three types of transmission rate. (a) Demands of general wards (b) Demands of ICUs (c) Number of Deaths (d) PHSM intervention levels.

A.6. Sensitivity analysis for drug effectiveness

In this subsection, we will test the impact of drug effectiveness under the conditions of 10 initial input infected cases, 40$\mathrm{\%}$ drug coverage rate, and a low virus vaccine escape rate. Figure A8 shows the sensitivity analysis results for all three DSPIs with drug effectiveness ranging between 50$\mathrm{\%}$ and 80$\mathrm{\%}$.

Figure A8. Sensitivity analysis for drug effectiveness under three DSPI scenarios, where the blue and orange curves correspond to the 50% and 80% effectiveness of the drug, respectively. The green horizontal dotted line represents the maximum value of general wards and ICUs. (a) Demand of general ward without DSPI (b) Demand of ICU without DSPI (c) Death without DSPI (d) Demand of general ward under weak DSPI (e) Demand of ICU under weak DSPI (f) Death under weak DSPI (g) Demand of general ward under strong DSPI (h) Demand of ICU under strong DSPI (i) Death under strong DSPI.

Without any DSPI, compared to 50$\mathrm{\%}$ drug effectiveness, the peak demand of general ward, the peak demand of ICU, and the cumulative deaths of 80$\mathrm{\%}$ drug effectiveness decreased respectively by 66$\mathrm{\%}$, 55$\mathrm{\%}$, and 57$\mathrm{\%}$, and the times of these peaks were the same. When weak DSPI is adopted, the reduction ratios of the above three values are 50$\mathrm{\%}$, 44$\mathrm{\%}$, 55$\mathrm{\%}$ respectively, but the occurrence time of the general ward's demand peak for the 80$\mathrm{\%}$ drug effectiveness will be a week later than the one for 50$\mathrm{\%}$ drug effectiveness. When strong DSPI is implemented, there will be two demand peaks. As we vary the drug effectiveness from 50% to 80%, the occurrence time of the first peak and the second peak will be 5 or 35 days later, while the values of the first demand peaks will be reduced by about 25$\mathrm{\%}$ and the reduction rate of the death is over 50$\mathrm{\%}$. In summary, the improvement of drug effectiveness will help alleviate the medical resource shortage. Therefore, developing more potent drugs for the treatment of COVID-19 is still of great importance for fighting the epidemic.

A.7. Sensitivity analysis for immune escape

This subsection will examine the effect of immune escape on the infection scale under conditions of initial infusion infection in 10 people, 40$\mathrm{\%}$ drug coverage, 80$\mathrm{\%}$ drug efficiency under 3 different DSPIs. Low and high immune escape scenarios are considered; similar analysis and parameters can be found in reference [28].

It can be seen from Figure A9 that in the case of no DSPI, the values of the demand peaks of general wards and ICUs will decrease 50$\mathrm{\%}$ and 30$\mathrm{\%}$ when the level of immune escape is low, but the largest demands are still above the maximum capacity of the existing resources. As for the total death, the low immune escape will bring a 22$\mathrm{\%}$ reduction compared to the high immune escape scenario. In the case of weak DSPI, the values of the first demand peaks of general wards and ICUs will decrease 10$\mathrm{\%}$ and 5$\mathrm{\%}$ when the level of immune escape is low and the numbers in the second peaks are 80$\mathrm{\%}$ and 70$\mathrm{\%}$ and the largest demands that occur in the first outbreak are still over the maximum capacity of the existing resources. The demands on general wards, ICUs, and total death in the case of strong DSPI show no significant differences between low immune escape and high immune escape. To sum up, changes in the immune escape only have a significant impact under a more relaxed prevention strategy. When the general public's autonomous consciousness of epidemic prevention becomes more stringent, there is no obvious difference in the demand for medical resources and deaths between the high and low escape scenarios.

Figure A9. Sensitivity analysis for levels of immune escape under three DSPI scenarios, where the blue and orange curves correspond to the low and high immune escape of the vaccine protectiveness. The green horizontal dotted line represents the maximum value of general wards and ICUs. (a) Demand of general ward without DSPI (b) Demand of ICU without DSPI (c) Death without DSPI (d) Demand of general ward under weak DSPI (e) Demand of ICU under weak DSPI (f) Death under weak DSPI (g) Demand of general ward under strong DSPI (h) Demand of ICU under strong DSPI (i) Death under strong DSPI.

A.8. Sensitivity analysis for initial imported cases

This subsection shows the sensitivity analysis of the initial infected population under 3 different DSPIs with the conditions of 40$\mathrm{\%}$ drug coverage, 80$\mathrm{\%}$ drug effectiveness, and a low viral vaccine escape rate. Three initial imported cases are considered: 10, 30, and 50.

It can be seen from Figure A10 that the change in the initial number of people infected will only shift the occurrence time of each peak of the epidemic without changing the number of peaks or the value of these peaks in all three DSPI cases. When the initial number of infections rises from 10 to 30, each peak comes about 6 days earlier, and the number will be 10 if it comes to 50 initial infected cases. Therefore, the number of initial input cases only has a little impact on the spread of the COVID-19 epidemic, which is consistent with the conclusions in [28].

Figure A10. Sensitivity analysis for initial input cases under three DSPI scenarios, where the blue, orange and green curves correspond to the initial input cases 10, 30, and 50, respectively. The red horizontal dotted line represents the maximum value of general wards and ICUs. (a) Demand of general ward without DSPI (b) Demand of ICU without DSPI (c) Death without DSPI (d) Demand of general ward under weak DSPI (e) Demand of ICU under weak DSPI (f) Death under weak DSPI (g) Demand of general ward under strong DSPI (h) Demand of ICU under strong DSPI (i) Death under strong DSPI .