The heterogeneous mixing model of COVID-19 with interventions.

The emergence of mutant strains of COVID-19 reduces the effectiveness of vaccines in preventing infection, but remains effective in preventing severe illness and death. This paper established a heterogeneous mixing model of age groups with pharmaceutical and non-pharmaceutical interventions by analyzing the transmission mechanism of breakthrough infection caused by the heterogeneity of protection period under the action of vaccine-preventable infection with the original strain. The control reproduction number Rc of the system is analyzed, and the existence and stability of equilibrium are given by the comparison principle. Numerical simulation was conducted to evaluate the vaccination program and intervention measures in the customized scenario, demonstrating that the group-3 coverage rate p3 plays a key role in Rc. It is proposed that accelerating the rate of admission and testing is conducive to epidemic control by further fitting data of COVID-19 transmission in real scenarios. The findings provide a general modeling idea for the emergence of new vaccines to prevent infection by mutant strains, as well as a solid theoretical foundation for mainland China to formulate future vaccination strategies for new vaccines. This manuscript was submitted as part of a theme issue on "Modelling COVID-19 and Preparedness for Future Pandemics".

Introduction

The infectious disease COVID-19 caused by SARS-CoV-2 has triggered a global pandemic (WHO, 2020), posing unprecedented challenges to public health systems. To deal with the potential impact of the COVID-19 on social, the Chinese government and other governments have implemented a variety of intervention measures. The interventions are divided into non-pharmaceutical interventions (NPIs) and pharmaceutical interventions. The former includes mass nucleic acid testing, contact tracking, isolation, quarantine, and other measures, whereas the latter includes vaccines and targeted medications (Lai et al., 2020, Flaxman et al., 2020, Ngonghala et al., 2020). The most effective and convenient method of preventing and controlling COVID-19 is timely vaccination, and this strategy has gained global acceptance. On 31 December 2020, Pfizer/BioNTech’s BNT162b2 vaccine was the first vaccine certified by WHO for emergency use against the original strain. Since then, ten vaccines have been listed on the Emergency Use List and their clinical trial data show excellent efficacy and immunogenicity, with robustest vaccine efficacy of 95% (Thomas et al., 2021). Whether individually or as a group, vaccination can contribute to reducing infection, morbidity, severe illness and death. The original strain has mutated at important genetic loci, and five VOCs-Alpha (B.1.1.7), Beta (B.1.351), Gamma (p.1), Delta (B.1.617.2) and Omicron (B.1.1.529) - have been reported against which existing vaccines have limited effectiveness (Nasreen et al., 2022, Andrews et al., 2022, Bernal et al., 2021). The vaccine is less effective in preventing infection, but it does provide some protection against morbidity, severe illness, and death. In this paper, a mathematical model was established for the prevention of infection by vaccines against the original strains, and also provides theoretical guidance for scientific distribution when new vaccines against mutant strains become available in the future. This model also applies to the emergence of new vaccines that can prevent infection in the future, but not for scenarios where the ability of vaccines to prevent serious illness and death should be considered.

In the normalization stage, the scientific and active implementation of immunization considering population heterogeneity has more crucial practical significance for the herd immunity. The immune system of different age groups is heterogeneous, resulting in significant heterogeneity of the immune response after vaccination (Collier et al., 2021), and the immunogenicity of the single-dose vaccine decreases with age (Faro-Viana et al., 2022). Secondly, heterogeneity in prevalence was evident from the onset of the pandemic in the age distribution of cases (ODriscoll et al., 2021, Castro and Singer, 2021), a phenomenon that can be attributed to differences in population susceptibility or infectivity leading to inconsistent clinical symptom tendencies (Davies et al., 2020). In this paper, the relative susceptibility of the population is considered. It is likewise a result of the heterogeneity of the immune system in age groups. Finally, heterogeneity of contact between age groups is also a major reason for considering age groups in our model (Cui et al., 2019, Cui et al., 2022, Feng et al., 2015, Feng et al., 2020). The social stage is characterized by peer contact, while the private stage is dominated by the interaction between all age groups (Seno, 2020, Prem et al., 2021). Therefore, this paper establishes a mathematical model against the background of mutual transmission of age groups to discuss how to vaccinate can prevent the spread of COVID-19 and achieve herd immunity. Overall, vaccines provide extensive protection against SARS-CoV-2 infection, but the protection period of vaccines for individuals is limited and heterogeneous. Individual immunity weakens over time, especially after the protection period, and being vaccinated can still develop breakthrough infections and spread the virus (Keehner et al., 2021, Juthani et al., 2021). Therefore, breakthrough infection in modeling is not negligible. At the same time, with herd immunity not truly achieved globally and variants persisting, the risk of a resurgence of COVID-19 remains high if NPIs are abandoned. The inclusion of NPIs in the model is indisputable.

Understanding the role of age in the transmission of COVID-19 is critical for implementing interventions targeting the original strain. In the absence of vaccines and drugs, using age-structured subpopulation model (Babajanyan and Cheong, 2021, Cuevas-Maraver et al., 2021) or incorporating real-time contact data into COVID-19 dynamics model (Iyaniwura et al., 2022) can accurately capture differences in symptoms and behavior across age groups. Monod et al. show that targeted interventions to aged 20 to 49 is an important consideration in halting resurgent COVID-19 (Monod et al., 2021). Therefore, it is extremely important to consider age structure into a dynamic model. The compartment model with age structure was not only used by Gozzi et al. to study the interplay of vaccines rollout and behavioral dynamics (Gozzi et al., 2021), but also by Moore to predict the possible long-term dynamics of SARS-CoV-2 during the planned vaccine rollout (Moore et al., 2021). Jiménez et al. took into account the ways in which people of different ages respond to the virus and evaluated two vaccination programs (Jiménez-Rodríguez et al., 2021). At the same time, novel analytical methods have sprung up. Lovell-Read et al. considered the effects of heterogeneities age-related factors on the SARS-CoV-2 transmission for infected individuals of ages by branching process model (Lovell-Read et al., 2022). In addition, the optimization theory also provides a research idea for determining the vaccination program under the different target of determining the original strain (Han et al., 2021, Das et al., 2021, Acuña-Zegarra et al., 2021). Using an age-stratified model paired with optimization algorithms, Matrajt et al. determined optimal vaccine allocation for different metrics (Matrajt et al., 2021b), and strategies with one and two doses of vaccine under various degrees of viral transmission (Matrajt et al., 2021a). However, the mathematical model of vaccination alone is not appropriate for the normalized stage, and it is necessary to incorporate NPIs into the model. Zou et al. combined vaccination and quarantine strategies to derive the key quarantine rate for controlling transmission and the vaccination rate for achieving herd immunity (Zou et al., 2022). Using age-stratified compartmental model, Bauer quantified the rate of increase in NPIs relative to vaccination progress without overwhelming the healthcare system (Bauer et al., 2021), Choi estimated the infection probability for each age group under different levels of social distancing implemented in Korea and investigated the effective age-dependent vaccination strategies (Choi et al., 2021). Among of these references mentioned above with extensive content, breakthrough infections caused by heterogeneity in the protection period have not been well incorporated into the model, as well as representative normal measures such as mass testing and tracking are not well coupled with the immune system heterogeneity. Therefore, this paper proposes a heterogeneous mixing age-group model with NPIs and imperfect vaccines to capture age-specific characteristics, and attempts to solve the following questions: In the case of breakthrough infections, with specific vaccines for specific strains, what are the effects of the combination of different vaccination strategies and levels of NPIs on the spread of COVID-19? Are there differences for achieving herd immunity? What factors are extremely important in terms of reproduction number?

The paper is structured in the following way. Section 2 introduces the heterogeneous mixing age-group model with NPIs and imperfect vaccines. Section 3 derives the formal expression of reproduction number in single age group- and in multiple age group, and demonstrates the threshold dynamics determined by . In Section 4 presents numerical results of different vaccination programs in hypothetical scenarios and data fitting using actual COVID-19 data in Shanghai. Finally, we make some discussions in Section 5.

Model formulation

Considering the heterogeneity of contact patterns, transmission risk, infection risk and vaccination risk, the model comprises four age groups 0–2, 3–17, 18–59, 60. Group-1 called the infant group, has a lower risk of participating in transmission but a higher risk of vaccination, so vaccination is not considered. Group-2 called the adolescent group. Group-3 is the young and middle-aged group, which is characterized by high risk of transmission. Group-4 is the elderly group with many underlying diseases and a higher risk of death.

The population in group- are sorted into susceptible , being vaccinated but non-responsive , being vaccinated and responsive , immunity invalidation , latent , asymptomatic and not been tested , symptomatic and not been tested , being tested positive and isolated , hospitalized and recovered . Fig. 1 illustrates the basic modeling structure of the transmission mechanisms of COVID-19 between age group- and , with details described in Section 2.1.

Fig. 1

Illustration of the basic modeling structure. The light green plane and light blue plane are the maps of COVID-19 transmission in age group- and , respectively. The rose lines represent age group- becomes age group- due to aging.

Assumptions and descriptions

The establishment of the model mainly refers to the multiple outbreaks of COVID-19 in mainland China. The NPIs used are numerous and accurate, among which testing and contact tracking are special under the “dynamic zero-COVID” policy. Testing can quickly identify the infectious source, and tracking can help break the transmission chain. These methods can effectively curb the spread of COVID-19. The transmission mechanism in single age group- is detailed in Fig.S1 (see Supplementary Material), including the vaccination process in Fig.S2 and the COVID-19 transmission process. The interpretations of the variables and parameters are described in Table 1.

Table 1

Description of variables and parameters.

Variables	Description
S_k	Number of unvaccinated individuals who are fully susceptible in age group-k.
Sk0	Number of vaccinated non-responders in age group-k.
V_k	Number of successful vaccinator in age group-k.
Vk0	Number of people whose immunity fails after the immune period in age group-k.
E_k	Number of individuals who are in latent period without infectious in age group-k.
A_k	Number of asymptomatic infections who have not been tested in age group-k.
I_k	Number of symptomatic infections who have not been tested in age group-k.
Q_k	Number of people isolated after being tested positive for nucleic acid in age group-k.
H_k	Number of hospitalized patients in age group-k.
R_k	Number of individuals who recovered in age group-k.
N	Total population.
Parameters	Description
b	Natural birth rate/Natural death rate.
λ1k,λ2k,λ3k	Susceptibility of S_k, Sk0, Vk0.
C_jk	The number of contacts that a person in age group-j has with individuals in age group-k.
ɛ	Correction factor for the incidence of asymptomatic infection.
p_k	Vaccination coverage in age group-k, p_1=0.
α_k	Probability of successful vaccination in age group-k.
ξ_k	Probability of lifetime immunity.
β	Probability of E_k,A_k,I_k is found during mass nucleic acid testing period.
η	Probability of E_k,A_k,I_k is found during close contacts tracking period.
ω_k	Aging rate from group k to k+1, ω_0=ω_4=0.
1/θ_k	Protection period of vaccine in age group-k.
σ_k	The rate of E_k→A_k in age group-k.
ϕ_k	The rate of E_k→I_k in age group-k.
μ_k,δ_k,γ_k	The recover rate of A_k, I_k and Q_k in age group-k.
ρ_k	Proportion of hospitalizations in age group-k.
1/υ_k	Length of hospital stay before recovery in age group-k.
ψkeq	The rate of E_k→Q_k in age group-k.
ψkaq	The rate of A_k→Q_k in age group-k.
ψkiq	The rate of I_k→Q_k in age group-k.

The detail of vaccination under the heterogeneity of immune protection is shown in Fig.S2. The vaccination starts from the susceptible, and the immune protection period after being vaccinated is 0, fixed and infinite. If immune period is zero, that is, being vaccinated but non-responsive , which is characterized by being able to participate in transmission. Otherwise, the vaccination is considered successful. Such people belong to , in which a small proportion of the population gets lifelong immunity into . Others become after the immune period due to immune wane, and they can still participate in transmission, but whether the risk of prevent morbidity will be lower or higher than is unknown. Therefore, and were infected as breakthrough infection. Moreover, the doses is not considered in this model, i.e., the first vaccination is regarded as the vaccinated population.

To construct the corresponding model in Fig. 1, the following assumptions and descriptions are given:

(1) The paper does not take into account the fatality rate, which is currently extremely low in mainland China. The latent period is not considered infectious.

(2) For contact pattern of group-, is the average contact number, is the proportion of the group- individual contacts with members of specific group-. So the average number of group- contacted by a group- individual per day is , denoted . Therefore, incidence rate between infectious of group- and susceptible of group- is or , where is susceptibility of . It is worth noting that needs the balance condition to be established (Cui et al., 2019).

(3) The group- will naturally into the group over time, which is the aging rate . All newborns are susceptible of group-1.

(4) The class is nucleic acids tested positive and isolated, its input includes two aspects. The first is characteristic of infectious diseases, and follow the epidemiological pathology naturally into . The second is social behavior in disease prevention and control, that is, mass testing and close contact tracking. In the population within , is proportion of the hospitalized treatment, and remains in for medical observation.

(5) Close contact tracking is to track and isolate close contacts of a person who is positive as soon as being found. Suppose a symptomatic infected in group- is found. In that case, the symptomatic infected belonging to group- among its close contacts is , where the probability of close contacts being tracked is .

Therefore, the model can be modified to adapt to different stages of prevention and control. As with the Wuhan outbreak in 2020, there were no medicines or vaccines in the early stage of COVID-19, and the measures were mainly NPIs. The schematic diagram of transmission mechanism is shown in Fig.S3(a). Then vaccination was imminent, the transmission mechanism of spread and vaccination at the same time is shown in Fig.S1, just like the Nanjing outbreak in 2021. Now that the vaccination work has been basically completed, and the transmission mechanism shown in Fig.S3(b) has been formed, which is applicable to the outbreak situation of Shanghai in 2022.

Table 1 summarizes the notions of the main epidemiological parameters, which are determined by three aspects. First, the vaccine itself determines the variables: immunity wane rate , vaccination success probability , lifelong immunity probability . Particularly, vaccination success probability is not vaccine efficacy. Secondly, social behavior determines the variables: contact patterns , the rate of positive detected in mass testing is , the probability of close contacts being tracked is . Thirdly, virus determines variables: natural evolution rate of , , , , .

Population dynamics of age groups

Suppose that the natural birth rate or death rate is proportional to the population and the total population is at a steady state during the short period of COVID-19 prevalent. The transition between age groups is only with aging as shown in Fig.S4, corresponding to the system where is the natural birth or death rate, represents the aging rate of group-, is the population size of group-. denotes the total population, which remains constant and can be scaled to 1. Each group model with demography, varies with time . The system Eq. (2.1) is one-order linearly differential equations system. Let are the proportion of each group population to the total population respectively. Therefore, the normalized system as follows

The system Eq. (2.2) exists a unique positive equilibrium , and is globally asymptotically stable.

See Appendix A for proof. □

Transmission dynamics of COVID-19 in age groups

Considering the transmission mechanism of COVID-19 in age groups, the schematic diagram is shown in Fig.S1. Appendix B presents the corresponding non-autonomous system Eq. (B.1) established by population size, which is depends on time through the function . By using the limit as time goes to infinity and let be measured as the proportion of each compartment population in group- to the total population. Therefore we replace Eq. (B.1) with the following limiting system where is the Kronecker delta function, if , otherwise .

Note: The superscript of parameters also changes accordingly to maintain the consistency of system symbols, i.e. . Moreover, , for .

Model analysis

In this section, the control reproduction number of single age group and multi-group of system Eq. (2.3) is given, and it is further proved that the multi-group reproduction number is the threshold for the stability of disease-free equilibrium and the existence of endemic equilibrium.

Existence of disease-free equilibrium

Model Eq. (2.3) is well posed, i.e., nonnegative initial conditions lead to nonnegative solutions for .

From Lemma 3.1, the feasible region can be shown as nonnegative cone where is subset of hyperplane . To simplify analysis, equations , and can be ignored because they are decoupled. Therefore, only consider the initial conditions in the bounded area and the area is positive invariant.

This section explains that the system Eq. (2.3) has a disease-free equilibrium and is unique. Set each derivative equal to 0 with , then solve the algebraic equation. For , that is, linear inhomogeneous system It is not difficult to know that the above equation system has a unique solution , where Next, similar methods are used to solve algebraic equations of other variables, and the equilibrium are shown in Appendix C. Then system Eq. (2.3) has a unique disease-free equilibrium , where is the disease-free equilibrium of age group- as below (see Box I).

where each of these components are the coordinate components of .

Threshold analysis

The basic reproduction number is the most important index in epidemic model. It measures the internal transmission ability of infectious diseases without intervention. An infectious will lead to an outbreak if , and put pressure on health care system. Otherwise, there will be no outbreak. The State Council joint prevention and control mechanism against COVID-19 issues control schemes to make threshold less than 1, which is also called control reproduction number . In this part, the thresholds of multi-group and single group are given respectively.

Multi-group threshold analysis

Using the next-generation matrix method to calculate the threshold of system Eq. (2.3). Only the disease compartments , , can be considered. Let be the increasing rate of secondary infection, is evolution operator and represents the internal evolution law (e.g., natural birth and death, and movements among compartments). Calculate the Jacobian F and V at the disease-free equilibrium respectively, see Appendix D.

Hence, as follows Following Van den Driessche and Watmough (2002), is the next generation matrix and set where is the spectral radius of .

The effective matrix elements for calculating the threshold is . But this system is more difficult to calculate , the specific form is not calculated in this paper. Moreover, let has the Z sign pattern, is non-negative and is non-singular M-matrix, following Theorem 2 in Van den Driessche and Watmough (2002) is holds, where is the eigenvalue of }.

Single-group threshold analysis

Similar to the previous method, this part only considers the control reproduction number of a single age group, which further explains the biological significance of the results. It is worth noting that the aging rate still needs to be considered, because it is a natural phenomenon and is part of the population dynamics. Through analysis, the control reproduction number for single group- as

Taking group-2 as an example, only the intra-group transmission is considered, and contact matrix takes the diagonal element . indicates the total infection of three classes of susceptible , and infected by two classes of infectious and . In the system, all positive cases are admitted to isolation and are no longer involved in transmission. Therefore, the period from infectious to entering class is regarded as the infectious period, further explained as Fig.S5. The other groups are similar, the difference is that in because babies are not vaccinated. There is no aging rate in , i.e. .

Through numerical calculation, can be obtained for . In the biological sense, this conclusion can also be explained. includes not only the intra group transmission ability, i.e. , but also transmission due to natural contact between groups, as well as transmission within and between groups due to aging into the next age group.

Stability analysis of disease-free equilibrium

If, the disease-free equilibriumof system Eq. (2.3) is locally asymptotically stable.

Analysis of the stability of the disease-free equilibrium using linearization around the equilibrium. The Jacobian matrix at disease-free equilibrium is The block matrices and are both lower triangular matrices whose eigenvalues are less than 0. In other words, and . If , then . Therefore, all eigenvalues of the Jacobian have negative real parts, i.e., . From Theorem 2 in Van den Driessche and Watmough (2002), the disease-free equilibrium is locally asymptotically stable. □

If, the disease-free equilibriumof system Eq. (2.3) is global attractiveness.

To illustrate the global attractiveness, it is only necessary to prove that , , () is established when . First, consider auxiliary system all eigenvalues of the corresponding Jacobian have negative real parts, equilibrium , is locally asymptotically stable and equal to , . Furthermore, the auxiliary system Eq. (3.6) is linear system, equilibrium is globally asymptotically stable. From system Eq. (3.6) is quasi-monotone increasing system, , , , is holds by using comparison principle. Therefore, , , , .

In the following prove that , , tend to 0 if . Define , where of which It is well known that . And is continuous with small , there exits enough small such that . The following discusses the solution of Eq. (3.8) if . Taking , there exits some such that , , for . Set is quasi-monotone increasing system. Because , the solution of system Eq. (3.8) are all close to 0 if . By comparison theorem, On the other hand, from the theory of asymptotic autonomous systems (Thieme, 1992) can be The proof is complete. □

By Theorem 3.1, Theorem 3.2, the following results must be true.

If, the disease free equilibriumof system Eq. (2.3) is globally asymptotically stable.

Existence of endemic equilibrium

The purpose of subsection is to illustrate the existence of endemic equilibrium by using threshold condition and consistent persistence theory. To this end, the following definitions are made

If , there exits positive real number , such that for any initial with , the solution satisfies and .

On the contrary, assume that there exits positive real , such that for all initial values belong to with . Now consider perturbation system It is quasi-monotone increasing system linear system. Similar to system Eq. (3.6), system Eq. (3.10) has a unique equilibrium , which is globally asymptotically stable.

Now just consider equation, is a continuous function of . It is possible to find some positive real such that and . Furthermore, there is sufficiently large , such that for . So is holds if . Similarly, is holds for . Next, consider following For the sake of convenience, denote Then the right system is quasi-monotone increasing system. From , i.e. . The function is continues for any small value , so can take enough small such that . Additionally, which is a contradiction. Therefore, there is holds. □

The above Theorem illustrate that at least one infected compartment is not 0. Furthermore, Wang and Zhao (2004) uniform persistence Theorem is used to demonstrate that the original system has a positive equilibrium.

If, system Eq. (2.3) has at least a positive equilibrium. That is, for initial , there is a positive real number such that and k 1,2,3,4.

To prove that Eq. (3.12) is established, we will demonstrate that system Eq. (2.3) is uniform persistence with respect to . For simplicity, denote By definition, are positively invariant and is the bounded closed set of . The trajectory starting from eventually enters the set , and it can be seen that the system Eq. (2.3) is point-dissipative (i.e. the existence of a bounded globally attracting set).

Let is the solution semiflow of system Eq. (2.3) on , , i.e. Next we prove Assume ((0), ,(0), ,(0),(0),(0)) and the definition of , it suffices to know If not, for and there exits some and , such that So the assumption might be .

Set group-1, group-2, group-3 and group-4 divide into two sets and , where is the set of age group- satisfying , and is the set of age group- with at least one compartment greater than 0 in compartment . Neither set is empty, because set is nonempty by definition of , and as can be seen from , is also a non-empty set.

Without loss of generality, for any , assuming , then Then there is enough small , such that In addition, set any , Then there is enough small , such that Further, then there is enough small , such that To sum up, if , there is not belong to , which contradicts the assumption that . Hence, the expression Eq. (3.13) is holds.

The open set is positive invariant set, and is globally asymptotically stable. By using abstract persistence theory, is isolated invariant set in , such that the orbit from is almost close to equilibrium . Moreover, no subset of forms a cycle in and , where . Therefore, according to Theorem 4.6 in Thieme (1993), system Eq. (2.3) is uniformly persistence with respect to , i.e., Eq. (3.12) is satisfied. □

Numerical simulations

In this section, contact patterns appropriate for the paper are first discussed in Section 4.1. Secondly, the vaccination strategies were evaluated from different perspectives through custom scenarios in 4.2 and data fitting was performed on the number of hospitalized patients through the actual scenario of the Shanghai outbreak in 2022 see Section 4.3.

Contact pattern

From several outbreaks in mainland China (Liaoning Province, Hebei Province, Nanjing City and Zhengzhou City), epidemiological survey of cases revealed heterogeneity in the age structure, see Fig.S6. Many efforts have been made to find out contact patterns between age groups. Regrouping the research results of 2020 (Prem et al., 2021), that is, 16 groups of data (left panel) are merged into 4 groups (right panel) in Fig. 2. The principle of reorganizing the contact pattern is that the number of contacts before and after grouping is equal. The contact pattern corresponds to an irreducible matrix , in which each small cell quantitatively describes the average number of contacts between two groups. The -axis represents the contactee, and the -axis is contactor. The darker the color, the greater the average number of contacts. Taking 16 groups of contact matrices as an example, the color of the diagonal is darker than that of the same row and column, the contact between the same age group is the most, which is in line with the characteristics of crowd behavior. It also qualitatively depicts the social behavior between the same age group, such as work and study.

Fig. 2

Illustration of contact pattern. The left panel is contact pattern of 16 groups, and the right panel is 4 groups. The degree of color is proportional to the number of contacts.

Vaccination strategies

The part evaluates the impact of the combination of age group vaccinations on the prevention and control of COVID-19 in the custom scenarios. Taking the proportion of vaccinations required to achieve herd immunity, effective reproduction number and control reproduction number as the evaluation criteria, the numerical results are used to evaluate the optimal vaccination strategy under different demands. In the custom scenarios, with the population distribution of Taiyuan, Shanxi Province as the background provided by the Bureau of Statistics (NBS, 2021). The contact pattern adopts the results of Section 4.1. The vaccination rate was estimated from the vaccination data of Shanxi Provincial Health Commission, and baseline parameters are shown in Table 2. The initial value of infection are assumed according to the age distribution of cases in other provinces and cities. The initial values as and the initial value of other variables is 0. In general, the custom scenario is representative and conforms to the characteristics of the outbreak in mainland China.

Table 2

Baseline values.

Parameters	Values	Source
b	2.1356e−05	NBS (2021)
C_ij	1.742.662.750.30.6310.25.650.420.342.4912.710.80.362.285.811.66	Derived from Prem et al. (2021)
ϵ	0.05	He et al. (2020)
ω_k	[1/418.65 1/3510.06 1/17482.69 0]	ω_k=(1Mean age in group-k - Mean age in group-(k-1)−1Average life span)×1365
θ_k	[1/180 1/180 1/180 1/180]	Center for Drug Evaluation of NMPA (2020)
σ_k	[1/3 1/3 1/3 1/3]	Davies et al. (2020)
ϕ_k	[1/5.1 1/5.1 1/5.1 1/5.1]	Davies et al. (2020)
μ_k	[1/14 1/14 1/14 1/14]	Assumed
δ_k	[1/12.9 1/12.9 1/12.9 1/12.9]	Derived from Davies et al. (2020)
γ_k	[1/10 1/10 1/10 1/10]	Wang et al. (2020)
λ1/2/3k	[0.4 0.39 0.82 0.825]	Bubar et al. (2021)
ψkeq	[1/8 1/8 1/8 1/8]	Davies et al. (2020)
ψkaq	[1/5 1/5 1/5 1/5]	Davies et al. (2020)
ψkiq	[1/2.9 1/2.9 1/2.9 1/2.9]	Davies et al. (2020)
β	1/3	The population of custom scenario is about 5.3 million, so β=1/3 is adopted from the State Council statement.
η	0.85	Assumed
ξ_k	[0.001 0.001 0.001 0.001]	Assumed

Evaluation criteria 1: Total immunized ratio and effective reproduction number

The effective reproduction number is time-dependent and represents the expectation of secondary cases arising from a primary case infected at time . The value changes over the course of outbreak and is used to assess the effectiveness of the current control (Thompson et al., 2019). In this section, can be made by vaccination. In the absence of vaccination , which means that herd immunity can be achieved if the vaccine coverage is 69.92%. Similar to the calculation method of , can be solved numerically. In this paper, the reasons such as vaccine failure and non-response are considered, the overall vaccine coverage must be greater than 69.92%. In this part, we propose 11 vaccination strategies, whose effects and costs are shown in Table 3. The first column No. represents the coding of the vaccination strategies, and the second column is the specific description of the strategies. The coverage threshold set by the experiment and the actual total vaccinated ratio were denoted as TH.Cov and TIR respectively. () is the moment when for the first time. In addition, TIR and () were selected as evaluation indexes of strategy selection.

Table 3

Results of vaccination strategies.

		(A)    1/θ = 10  months			(B)    1/θ = 12  months
No.	Vaccination Strategy	TH.Cova	T(R_t<1)b	TIRc	TH.Cov	T(R_t<1)	TIR
M1	p_3 = 0.01	–d	–	–	–	–	–
		–	–	–	–	–	–
M2	p_2=p_4 = 0.01	–	–	–	–	–	–
		–	–	–	–	–	–
M3	p_3=p_4 = 0.01	–	–	–	–	–	–
		–	–	–	–	–	–
M4	p_2=p_3 = 0.01	–	–	–	–	329	78.06%
		–	–	–	–
M5	p_2=p_3=p_4 = 0.01	–	315	94.34%	–	228	88.01%
		–			–
M6-1	p_2=p_4 = 0.01 →ep_3 = 0.01	0.75	441	95.52%	0.75	368	88.57%
M6-2		0.1	325	94.45%	0.1	238	87.94%
M7-1	p_2=p_3 = 0.01 →p_4 = 0.01	0.75	472	95.96%	0.75	381	90.68%
M7-2		0.1	327	94.42%	0.1	239	88.05%
M8-1	p_3=p_4 = 0.01 →p_2 = 0.01	0.75	473	95.51%	0.75	382	90.58%
M8-2		0.1	327	94.42%	0.1	239	88.05%

TH.Cov means threshold of coverage.

T() is the time when for the first time.

TIR is the abbreviation of total immunized ratio, which is equal to total number of vaccinations divided by N.

– means no data results during the whole experiment period.

stands for switching vaccination groups.

The strategies also takes into account the extremely important factor of vaccine protection period . Once the protection period is over, the vaccinated individuals will still regard as “susceptible” and be infected. On 14 August 2020, the Center for Drug Evaluation of NMPA (2020) issued the Guidelines for Clinical Evaluation of Novel Coronavirus Preventive Vaccines (Trial), stating that “vaccines preferably provide protection for 1 year or more, with at least 6 months”. The experiments select 10 months and 12 months respectively, corresponding to situations A and B in Table 3. In addition, by the single day vaccination amount in Shanxi Province, the daily dose accounts for about 0.012 of the total population. Therefore, the vaccination proportion is selected to take . Fig.S7 and Fig.S8 show the curves of and immunization number under different vaccination methods. The right coordinate is the cumulative number of vaccinated, which is represented by solid lines of different colors. The left coordinate is , indicated by dark green dotted line that becomes solid if .

In Table 3, the No. M1 to M5 strategies from single-group single-stage (M1), double-group single-stage (M2-M4), and three-group single-stage (M5) to compare the effects of vaccination schemes for 20 months. Neither single-group single-phase nor double-group single-phase regimen achieved herd immunity throughout the trial period if months. But if three groups were vaccinated simultaneously, on the 315th day, see Fig.S7(a). In view of the vaccine supply, if doses are insufficient in the early stage and only some people can be vaccinated first, so which group is most effective in controlling an outbreak is worth investigating.

The No. M6-1 to M8-2 are multi-group multi-stage mixed vaccination methods. The vaccination coverage threshold was set to be 0.75 and 0.1, which were higher and lower than the herd immunity threshold 69%, respectively. If TH.Cov 0.75 and months, the numerical results show that M6-1 achieves herd immunity earliest, although the vaccination proportion is as high as 95.52%. The strategy “ 0.01 0.01” corresponds to Fig.S7(c), group-2 and group-4 were selected to be vaccinated at a rate of 0.01 in Phase I. When the number of vaccinations in the group reaches the established coverage threshold 0.75, stop (group-2) vaccination, and to vaccinate the other two groups (group-4 and group-3) at a rate of 0.01 in Phase II. If threshold 0.75 is reached in a certain group again, stop the vaccination (group-4) and the other group (group-3) is continued to be vaccinated until its coverage reaches 0.75, the vaccination is stopped. At this time, the three groups is vaccinated at the rate of 0.01 until the end. If TH.Cov 0.1, the effects of three strategies are almost the same. Compared with TH.Cov 0.75, both and TIR are decreased, which can be explained. In a short time, the immunization proportions of the three groups all reached 0.1, and thereafter it is actually M5. The indicator () in the table are extremely similar, with a maximum difference of 12 days compared with M5. However, is this time difference sufficient? It is not known to buy time for sufficient supply of vaccines. To be on the safe side, it may be safer to adopt the M6-1 scheme in the early stage for 10 months.

Next, we will discuss the situation where 12 months, and there is a significant improvement in each scheme compared to 10 months (Fig.S7, Fig.S8). The M4 scheme achieved herd immunity on the 329th day, when only 78.06% of the total population was vaccinated. Programs M7-1 and M8-1 completed herd immunization 91 days in advance, and the coverage of M6-1 is reduced by 6.95% compared with 10 months. For scheme M5, the three groups could reach after 228 days of continuous inoculation at the same time. The mechanism that causes these changes is already evident, and the immune protection period is prolonged, which is equivalent to the slowing the outflow of the effectively immunized population. Ideally, if the duration of protection approaches infinity or is longer than life expectancy, it is equivalent to lifelong immunity. In this way, herd immunity can be achieved if the vaccination coverage reaches 69.92%. It is worth noting that the curve drops sharply as soon as the “group-3” are vaccinated in Fig.S7.

In summary, if the immune protection period is short, M6-1 strategy can be adopted. Otherwise, M5 takes precedence, followed by M4.

Evaluation criteria 2: control reproduction number

The previous subsection evaluated the vaccination strategies in proposed scenarios with as the index. This subsection takes control reproduction number as the evaluation index and focuses on the influence of controllable factors on . It is not only related to the characteristics of the vaccine, but also to social behavior and the internal transmission mechanism of the disease. Within the scope of limited capabilities, the influence of controllable variables on the threshold and whether it has guiding significance for disease prevention and control should be considered in scientific research. For the baseline parameters at transmission in Table 2. The contribution of vaccination rate , , and social behavior variable to was assessed by univariate, bivariate, and multivariate methods.

Case 1. Univariate evaluation: continuous vaccination only for single age groups.

Ignoring the long-term adverse effects of vaccine on humans, only the effect of population heterogeneity of vaccination on was investigated. In general, vaccine coverage rate has positive effect on the reduction of , as shown in Fig.S9. The light pink is drop area, and the basic reproduction number is 3.327 in the absence of vaccination. When the coverage rate gradually increased to 0.01 in different age groups, decreased to different degrees, among which decreased by 59.2% if group-2 was inoculated alone. The change of is not obvious after the plunging area, but the partial enlarged view (i.e. ) shows that the three threshold curves still show a downward trend. What is more worth mentioning is that the effect of vaccination among adolescents and young adults is significant, while in elderly is slightly weaker.

Case 2. Bivariate evaluation: continuous vaccination for double age groups.

The contribution of different combinations to is assessed by using a dual vaccination method, that is, inoculating both age groups simultaneously. Using the characteristics of the 2D contour map is same line equal height, different combination methods can present a consistent effect, which is lower than that at the origin. The vertical level colorbar depicts , and the change of color level corresponds to the value change of . Fig. 3 panel (a) shows simultaneous inoculation of “group-3 and group-4”, the former plays a leading role, and the elderly group only slightly affects the threshold change when . Panel (c) is “group-4 and group-2”, the elderly has limited influence on the results, similar to (a). The most effective strategy is “group-3 and group-2”. Not only does decrease the most, but also the combination to reach a certain threshold is more flexible. Especially in the early stage of the vaccine, the phase III clinical trial of the young group is the earliest and the data is the most complete. Taking the young and middle-aged group as the leading group of vaccination coverage can reduce the coverage of other groups as much as possible.

Fig. 3

The 2D contour map of with . Panels (a)–(c) show the effect of simultaneous inoculation of “group-3 and group-4”, “group-3 and group-2”, and “group-4 and group-2” on respectively. The colorbar represents the value of .

Using gradient to further characterize the effect of coverage rate on . Assume that the coverage of two groups is () and the amount of change is (), the binary function linearly approximate and in vector form where is the gradient of at . changes the fastest along the gradient direction, the rate of change is the gradient modulus . In Fig. 4, changes little when and increases to a certain. In particular, panel (c) shows the direction of negative gradient along which changes the fastest and always has a positive effect on reducing the threshold. In Fig.S10 panel (c), the arrow of any initial value is always perpendicular to the axis. That is, in the “group-3 and group-4” combination mode, it is the most effective vaccination way to ensure that the young and middle-aged group should be vaccinated as much as possible, and the conclusion is consistent with panels (a) - (b).

Fig. 4

Panels (a)–(c) represent the , magnitudes of gradient , negative gradient direction represented by streamline of “group-2 and group-3” strategy respectively.

Case 3. Multivariate evaluation: continuous vaccination for three age groups.

The 3D slice map in Fig.S11 illustrates the influence of three groups of immunization coverage rates on . The three-dimensional coordinates are , and respectively. Taking slices of three dimensions respectively, the overall change of is obvious, with , playing a decisive role.

Case 4. Multivariate evaluation: influence of and on .

In the model, is the mass testing rate, i.e. the nucleic acid testing speed of all residents in a city with cases. The State Council has clarified that cities with a permanent population of less than 5 million have the ability to complete all nucleic acid tests within 2 days by coordinating the resources in the province. For those over 5 million have the ability to complete full testing within 3–5 days by coordinating provincial resources and national support. If the testing ability is increased to 2 days, i.e. , the effect of experimental simulation is extremely obvious (compare Fig. 3 and Fig.S12). The of dual-age combined vaccination mode decreased significantly, and the infectivity decreased by 30.2% compared with . The threshold of “group-3 and group-2” gradually tends to below 1 with the increase of coverage. Even if the vaccination is not carried out, reduced from 3.327 to 2.233. It can be seen that the controllable variable plays an extremely important role in the prevention and control of the epidemic. In particular, the transmission capacity is 71.2% lower than that only the whole people were tested but not vaccinated if and .

The relationship between the three variables and is shown in Fig. 5, where 0:0.2:1, 0:0.0001:0.001 and fixed other parameters. Light pink is the threshold surface, the upper part of the surface is . Rapid universal testing and small vaccinations for adolescents and young and middle-aged groups, both of the two measures are sufficient to against COVID-19 and to contain its spread. The color on each slice of panel (b) corresponds to the value of colorbar, which matches and remains consistent with the policy opinions reflected in panel (a). If , i.e., complete all testing within one day, . It is the most desirable scenario in the epidemic prevention and control, which has greatly suppresses the spread of the epidemic. However, this measure has also increased the pressure on medical staff and testing personnel. Subsequently, the mixed sampling detection technology issued by the State Council can be regarded as an important measure to increase .

Fig. 5

Left panel is the 3D curved surface of , the upper part of light pink surface is . Right panel is slice map correspond with left panel, X–Y plane is coverage rate combination of group-2 and 3, -axis represents the change of , the colorbar is the value of .

In the COVID-19 prevention and control measures, we should “run faster than the virus” - increase the to speed up the testing speed. On the other hand, we should “fight against the virus” - increase “”, encourages voluntary COVID-19 vaccinations, and works to ensure all people eligible for vaccination have access to it. To implement regular epidemic prevention and control measures, vaccination is the general trend. Only by establishing a universal immunity barrier can promote the comprehensive restoration of social economy and living order.

A case study: the Shanghai COVID-19 outbreak

On February 24, 2022, Shanghai reported its first asymptomatic local infection case. On April 11, differentiated prevention and control by region were carried out based on the actual situation, with a total of more than 600,000 local positive cases up to now. Therefore, in order to visually reproduce the situation of COVID-19 transmission in Shanghai, we eliminate the age heterogeneity within system Eq. (B.1) and obtain the adjusted model Eq. (E.1) in Appendix E. As of February 15, Shanghai has completed 55.49 million doses of COVID-19 vaccines, covering 95.1% of the city’s permanent population. In the process of data fitting, we did not consider continuing to vaccinate the population, so . The number of local hospitalized cases notified by Shanghai Municipal Health Commission from April 11 to May 16, 2022 was fitted based on the least square method.

Assume that is the fitting the number of local cases being treated in hospital at time , and its change with time is determined by the following ordinary differential equation The reported local hospitalized cases at time is . By using the least squares method to estimate the unknown parameter values to minimize the objective function The specific values of the initial variables and parameters used are shown in Table 4.

Table 4

Descriptions and estimations of initial variables and parameters of Shanghai outbreak.

Initial variables	Values	Sources	Initial variables	Values	Sources
S(0)	1 220 144	NBS (2021)	S^0(0)	0	Assumed
V(0)	7 367 347	Calculated	V^0(0)	16 313 409	Calculated
E(0)	200 000	LS	A(0)	28 500	Assumed
I(0)	2982	Assumed	Q(0)	205 617	Shanghai (2021)
H(0)	6921	Shanghai (2021)	R(0)	1786	Shanghai (2021)
Parameters	Values	Sources	Parameters	Values	Sources
C	2.3	Zhang et al. (2020)	b	1.9726e−05	NBS (2021)
ϵ	0.35	Cai et al. (2022)	θ	0.0081	Calculated
σ	1/1.2	Cai et al. (2022)	ϕ	1/1.2	Cai et al. (2022)
μ	1/5.64	Cai et al. (2022)	δ	1/8.2	Calculated
γ	1/10	Wang et al. (2020)	υ	1/6	Cai et al. (2022)
λ_1	1	Cai et al. (2022)	λ_2	1	Assumed
λ_3	0.6393	LS	ψ_E	1/4.5	LS
ψ_A	1/4	Assumed	ψ_I	1/2.2	Cai et al. (2022)
β	1/3	Assumed	η	0.85	Assumed
α	0	Assumed	ξ	0	Assumed
ρ	0.0137	LS	p	0	Assumed

As shown in Fig. 6(a), the fitted daily number of hospitalized cases did well with the reported data, suggesting that the adjusted model Eq. (E.1) is helpful to explain the transmission of Shanghai. More interestingly, the proportion of positive cases hospitalized in Table 4 is 0.01317, which shows that vaccination reduced the emergency hospitalizations after infection with Omicron. It can be seen from Fig. 6(b) that accelerating the admission rate of symptomatic and speeding up the testing rate have extremely important positive effects on controlling the spread of COVID-19 in Shanghai.

Fig. 6

The fitting results and PRCCs for Shanghai. Panel (a) is the fitting results of the reported local hospitalized cases in Shanghai from April 11 to May 16, 2022. Panel (b) is the value of partial rank correlation coefficient between control reproduction number and parameters , , , .

Discussion

The main purpose of this paper is to provide a general modeling framework for the spread of COVID-19 in China in the presence of NPIs and imperfect vaccines to prevent infection against the original strain, to preliminarily search for appropriate vaccination scheme. Using the numerical results of different vaccination programs, the effectiveness of active vaccination was fully demonstrated from the perspective of control reproduction number and effective reproduction number . However, the proposed model still has limitations in application scenarios. First, this work is not suitable to explain the role of vaccines in reducing infection for Omicron, but does not rule out the emergence of vaccines against mutant strains. The paper will prepare for the emergence of new vaccines in the future and provide modeling ideas. Secondly, although the study has given vaccination recommendations about different demands, there are still great challenges in the “last mile of battle against COVID-19”. The model cannot characterize the role of vaccines in reducing severe illness and death and fails to highlight the importance of vaccine doses in slowing the outbreak. Thirdly, NPIs are highly targeted and more suitable for mainland China. The fifth wave of the pandemic in Hong Kong showed significant age heterogeneity in severe illness and mortality. It will be future work to determine whether vaccination dose plays a significant role in preventing severe illness and death in an age stratification. Finally, the model takes the positive or not as the isolation standard, and there is no class corresponding to the actual data of confirmed or asymptomatic cases. In short, vaccines are not a panacea, but they are not impossible either. Implement normalized epidemic prevention and control measures, and comprehensively promote the construction of immune barrier. Every cloud has a silver lining, and one day we can win the battle against COVID-19.

CRediT authorship contribution statement

Moran Duan: Methodology, Software, Formal analysis, Writing – original draft. Zhen Jin: Conceptualization, Supervision, Funding acquisition, Writing – review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.