Modeling outbreaks of COVID-19 in China: The impact of vaccination and other control measures on curbing the epidemic

ABSTRACT

This study aims to examine the development trend of COVID-19 in China and propose a model to assess the impacts of various prevention and control measures in combating the COVID-19 pandemic. Using COVID-19 cases reported by the National Health Commission of China from January 2, 2020, to January 2, 2022, we established a Susceptible-Exposed-Infected-Asymptomatic-Quarantined-Vaccinated-Hospitalized-Removed (SEIAQVHR) model to calculate the COVID-19 transmission rate and R_t effective reproduction number, and assess prevention and control measures. Additionally, we built a stochastic model to explore the development of the COVID-19 epidemic. We modeled the incidence trends in five outbreaks between 2020 and 2022. Some important features of the COVID-19 epidemic are mirrored in the estimates based on our SEIAQVHR model. Our model indicates that an infected index case entering the community has a 50%–60% chance to cause a COVID-19 outbreak. Wearing masks and getting vaccinated were the most effective measures among all the prevention and control measures. Specifically targeting asymptomatic individuals had no significant impact on the spread of COVID-19. By adjusting prevention and control parameters, we suggest that increasing the rates of effective vaccination and mask-wearing can significantly reduce COVID-19 cases in China. Our stochastic model analysis provides a useful tool for understanding the COVID-19 epidemic in China.

KEYWORDS:

• COVID-19
• public health
• compartmental model
• stochastic

Background

The first case of COVID-19 was reported on December 1, 2019 in Wuhan city of Hubei province, China, after which the virus quickly spread to other parts. The novel coronavirus spread rapidly worldwide, leading to a pandemic. China implemented strict measures, including wearing masks and quarantining patients, resulting in successful control of the COVID-19 outbreak within its borders. However, since March 2020, China has experienced multiple COVID-19 outbreaks on several occasions.

Since the emergence of COVID-19, researchers from various fields have studied and analyzed its data. Infectious disease models, such as SIR, SEIR, and improved SEIR models, have been widely established, with the SEIR model being the most commonly used.^1–5 These models have been used to estimate the basic reproduction number of the virus and predict the epidemic trend in different regions.^6–9 However, deterministic models can no longer fully capture the complex dynamics of diseases like COVID-19 due to the constantly changing and unpredictable environment. Therefore, researchers have increasingly turned to stochastic infectious disease models, which take into account random interference in the transmission of the disease.^10–14 These studies have shown that COVID-19 has a high probability of spreading to most people, making it crucial to implement rapid and efficient contact detection, management, and timely case isolation to control outbreaks.

In the following section, we constructed a deterministic model and a stochastic model, which were fitted with actual data to analyze the development trend of COVID-19 in China. We evaluated the effects of various prevention and control measures and compared the results of the stochastic model with the deterministic model.

Methods

This study focuses on COVID-19 cases reported by the National Health Commission of the People’s Republic of China between March 4, 2020, and January 21, 2022. We will analyze the severe COVID-19 period between 2020 and 2022, calculate the effective reproduction number of different time periods, compare the COVID-19 outbreaks and their development trend in China, and evaluate the effectiveness of prevention and control measures. Furthermore, we will compare the differences between the stochastic and deterministic models.

The deterministic mode

A typical approach when modeling the spread of an infectious disease in a community is to formulate a compartmental model. We divide the total human population N, into nine compartments: susceptible individuals S, exposed individuals E, symptomatic individuals I, asymptomatic individuals A, quarantined susceptible Q_s, quarantined exposed Q_e, vaccinated individuals V, hospitalized individuals H and recovered individuals R. The dynamics explained in this subsection are displayed in Figure 1, also see Tables 1 and 2 for the meaning of the remaining parameters.

Figure 1. Flow diagram of the models. Each box represents an individual compartment and the arrows represent transitions between the compartments. Nine individual compartments are susceptible individuals S, exposed individuals E, symptomatic individuals I, asymptomatic individuals A, quarantined susceptible Q_s, quarantined exposed Q_e, vaccinated individuals V, hospitalized individuals H and recovered individuals R, respectively. The transition rates are shown next to the arrow.

Table 1. Description of variables in the model (1.1).

Variables	Descriptions
N	Number of total population
S	Number of susceptible individuals
E	Number of exposed individuals
I	Number of symptomatic individuals
A	Number of asymptomatic individuals
Qs	Number of susceptible individuals in quarantine
Qe	Number of exposed individuals in quarantine
V	Number of vaccinated individuals
H	Number of hospitalized individuals
R	Number of recovered individuals

Table 2. The meaning and value of each parameter in the model.

Parameter	Meaning	Units	Value
β1	Rate of transmission of susceptible with exposed human	Persons^−1day^−1	Estimated
β2	Rate of transmission of susceptible with symptomatic human	Persons^−1day^−1	Estimated
β3	Rate of transmission of susceptible with asymptomatic human	Persons^−1day^−1	Estimated
μ	Rate of effective vaccination	–	Estimated
φ	Rate of effective wearing masks	–	0.05
Os	Quarantine rate of susceptible individuals	–	Estimated
Oe	Quarantine rate of exposed individuals	–	Estimated
ω	Rate of incubation	Day^−1	[0.071,0.33]^Citation15,Citation16
p	Fraction of exposed individuals showing clinical symptoms	–	[0.3–0.85]^Citation16
λ	Rate of hospitalized of symptomatic individuals	Day^−1	[0.02,0.1]^Citation13
κ	Rate of hospitalized of asymptomatic individual	Day^−1	Estimated
γ1	Recovery rate of symptomatic	Day^−1	0.3^Citation13
γ2	Recovery rate of asymptomatic	Day^−1	0.2^Citation17
γ3	Recovery rate of hospitalized individuals	Day^−1	0.3^Citation18,Citation19

In order to build this model (1.1), we made the following assumptions:

1. Susceptible persons can be infected by exposed individuals, symptomatic individuals and asymptomatic individuals, and become exposed after the incubation period. The transmission rate of exposed individuals is ${\beta }_1$, symptomatic individuals is ${\beta }_2$ and asymptomatic individuals is ${\beta }_3$.

2. The individuals may become symptomatic or asymptomatic after being infected by SARS-CoV-2, and the probability of becoming symptomatic is p.

3. The probability that a susceptible individual will be quarantined is $o_s$, and that an exposed individual will be quarantined is $o_e$.

4. ${\gamma }_1$, ${\gamma }_2$, ${\gamma }_3$ are the recovery rates of symptomatic, asymptomatic, and hospitalized individuals, respectively.

5. Effective mask wearing rate in the population is $\mathit{\phi }$.

6. Susceptible people have $\mathit{\mu }$ probability of vaccination and become immune in a short time after vaccination.

7. The recovered individuals do not acquire the infection again and are not going to infect others.

Based on these assumptions, we proposed a schematic diagram (Figure 1) and provided the corresponding model Equations (1.1):

$\mathit{dS}/\mathit{dt}=-\left\{\mathit{S}\left[1-\left(\mathit{\varphi }+\mathit{\mu }+o_S\right)\right]\right\}\left({\beta }_1\mathit{E}+{\beta }_2\mathit{I}+{\beta }_3\mathit{A}\right)-\mathit{\mu S}-o_S\mathit{S}$

$\mathit{dE}/\mathit{dt}=\left\{\mathit{S}\left[1-\left(\mathit{\varphi }+\mathit{\mu }+o_S\right)\right]\right\}\left({\beta }_1\mathit{E}+{\beta }_2\mathit{I}+{\beta }_3\mathit{A}\right)-o_e\mathit{E}-\mathit{\omega E}$

$\mathit{dI}/\mathit{dt}=\mathit{\omega pE}-\mathit{\lambda I}-{\gamma }_1\mathit{I}$

$\mathit{dA}/\mathit{dt}=\mathit{\omega }\left(1-\mathit{p}\right)\mathit{E}-\mathit{\kappa A}-{\gamma }_2\mathit{A}$

$\mathit{d}{Q}_s/\mathit{dt}=o_S\mathit{S}$

$\mathit{d}{Q}_e/\mathit{dt}=o_e\mathit{E}$

$\mathit{dV}/\mathit{dt}=\mathit{\mu S}$

$\mathit{dH}/\mathit{dt}=\mathit{\lambda I}+\mathit{\kappa A}-{\gamma }_3\mathit{H}$

$\mathit{dR}/\mathit{dt}={\gamma }_1\mathit{I}+{\gamma }_2\mathit{A}+{\gamma }_3\mathit{H}$

Effective reproduction number

According to next-generation matrix method, we determine the expression for the effective reproduction rate R_t. In order to this derivation, we add birth rate($\mathrm{\Lambda }$) and death rate(c) into the model (1.1), get model (1.2):

$\mathit{dS}/\mathit{dt}=\mathrm{\Lambda }-\left\{\mathit{S}\left[1-\left(\mathit{\varphi }+\mathit{\mu }+o_S\right)\right]\right\}\left({\beta }_1\mathit{E}+{\beta }_2\mathit{I}+{\beta }_3\mathit{A}\right)-\mathit{\mu S}-o_S\mathit{S}-\mathit{cS}$

$\mathit{dE}/\mathit{dt}=\left\{\mathit{S}\left[1-\left(\mathit{\varphi }+\mathit{\mu }+o_S\right)\right]\right\}\left({\beta }_1\mathit{E}+{\beta }_2\mathit{I}+{\beta }_3\mathit{A}\right)-o_e\mathit{E}-\mathit{\omega E}-\mathit{cE}$

$\mathit{dI}/\mathit{dt}=\mathit{\omega pE}-\mathit{\lambda I}-{\gamma }_1\mathit{I}-\mathit{cI}$

$\mathit{dA}/\mathit{dt}=\mathit{\omega }\left(1-\mathit{p}\right)\mathit{E}-\mathit{\kappa A}-{\gamma }_2\mathit{A}-\mathit{cA}$

$\mathit{d}{Q}_s/\mathit{dt}=o_S\mathit{S}-\mathit{c}{Q}_s$

$\mathit{d}{Q}_e/\mathit{dt}=o_e\mathit{E}-\mathit{c}{Q}_e$

$\mathit{dV}/\mathit{dt}=\mathit{\mu S}-\mathit{cV}$

$\mathit{dH}/\mathit{dt}=\mathit{\lambda I}+\mathit{\kappa A}-{\gamma }_3\mathit{H}-\mathit{cH}$

$\mathit{dR}/\mathit{dt}={\gamma }_1\mathit{I}+{\gamma }_2\mathit{A}+{\gamma }_3\mathit{H}-\mathit{cR}$

By adding multiple equations in the model (1.2), we get:

$(\mathit{S}+\mathit{E}+\mathit{I}+\mathit{A}+Q_s+Q_e+\mathit{V}+\mathit{H}+\mathit{R}{)}^{\prime }\le \mathrm{\Lambda }-\mathit{c}\left(\mathit{S}+\mathit{E}+\mathit{I}+\mathit{A}+Q_s+Q_e+\mathit{V}+\mathit{H}+\mathit{R}\right)$

So, there is the following set:

$\underset{\mathit{t}\to +\mathrm{\infty }}{\mathbf{lim}}(\mathit{S}+\mathit{E}+\mathit{I}+\mathit{A}+Q_s+Q_e+\mathit{V}+\mathit{H}+\mathit{R})\le \frac{\mathrm{\Lambda }}{c}.$

$\mathrm{\Omega }=\left\{\left(\mathit{S},\mathit{E},\mathit{I},\mathit{A},Q_s,Q_e,\mathit{V},\mathit{H},\mathit{R}\right)\in R_{+}^9:\mathit{S}+\mathit{E}+\mathit{I}+\mathit{A}+Q_s+Q_e+\mathit{V}+\mathit{H}+\mathit{R}\le \frac{\mathrm{\Lambda }}{c}\right\}$

is the positively invariant set of the model (1.2).

According to next-generation matrix method, let $\mathit{x}={\left(\mathit{E},\mathit{I},\mathit{A}\right)}^T$.

So, there is the following set:

$\mathit{dx}/\mathit{dt}=\mathit{F}\left(\mathit{x}\right)-\mathit{V}\left(\mathit{x}\right).$

The new infectious terms and transition terms of the model (1.2) are given by

$\mathit{F}\left(\mathit{x}\right)=\left(\begin{array}{c}\left\{\mathit{S}\left[1-\left(\mathit{\varphi }+\mathit{\mu }+o_S\right)\right]\right\}\left({\beta }_1\mathit{E}+{\beta }_2\mathit{I}+{\beta }_3\mathit{A}\right)\\ 0\\ 0\end{array}\right),$

$\mathit{V}\left(\mathit{x}\right)=\left(\begin{array}{c}{o}_e\mathit{E}+\mathit{\omega E}+\mathit{cE}\\ -\mathit{\omega pE}+\mathit{\lambda I}+{\gamma }_1\mathit{I}+\mathit{cI}\\ -\mathit{\omega }\left(1-\mathit{p}\right)\mathit{E}+\mathit{\kappa A}+{\gamma }_2\mathit{A}+\mathit{cA}\end{array}\right).$

$\mathit{F}$ and $\mathit{V}$ are the Jacobi matrices of $\mathit{F}\left(\mathit{x}\right)$ and $\mathit{V}\left(\mathit{x}\right)$. It follows that

$\mathit{F}=\left(\begin{array}{ccc}\mathit{S}{\beta }_1\left(1-\left(\mathit{\varphi }+\mathit{\mu }+o_S\right)\right)\frac{\mathrm{\Lambda }}{\mu +o_S+\mathit{c}}& \mathit{S}{\beta }_2\left(1-\left(\mathit{\varphi }+\mathit{\mu }+o_S\right)\right)\frac{\mathrm{\Lambda }}{\mu +o_S+\mathit{c}}& \mathit{S}{\beta }_3\left(1-\left(\mathit{\varphi }+\mathit{\mu }+o_S\right)\right)\frac{\mathrm{\Lambda }}{\mu +o_S+\mathit{c}}\\ 0& 0& 0\\ 0& 0& 0\end{array}\right),$

$\mathit{V}=\left(\begin{array}{ccc}{o}_e+\mathit{\omega }+\mathit{c}& 0& 0\\ -\mathit{\omega p}& \mathit{\lambda }+{\gamma }_1+\mathit{c}& 0\\ -\mathit{\omega }\left(1-\mathit{p}\right)& 0& \mathit{\kappa }+{\gamma }_2+\mathit{c}\end{array}\right),$

Where

$\mathit{F}{V}^{-1}=\left(\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ 0& 0& 0\\ 0& 0& 0\end{array}\right),$

$a_{11}={\beta }_1\frac{[1-(\varphi +\mu +o_S)]\mathrm{\Lambda }}{(o_e+\mathit{\omega }+\mathit{c})(\mathit{\mu }+o_S+\mathit{c})}+{\beta }_2\frac{[1-(\varphi +\mu +o_S)]\mathit{\omega p}\Lambda }{(o_e+\mathit{\omega }+\mathit{c})(\mathit{\lambda }+{\gamma }_1+\mathit{c})(\mathit{\mu }+o_S+\mathit{c})}\mathit{\hspace{0.17em}}+{\beta }_3\frac{[1-(\varphi +\mu +o_S)]\mathit{\omega }(1-\mathit{p})\mathrm{\Lambda }}{(o_e+\mathit{\omega }+\mathit{c})(\mathit{\kappa }+{\gamma }_2+\mathit{c})(\mathit{\mu }+o_S+\mathit{c})},$

$a_{12}=\mathit{S}{\beta }_2\mathrm{\Lambda }\frac{1-\left(\varphi +\mu +o_S\right)}{\left(\mathit{\lambda }+{\gamma }_1+\mathit{c}\right)\left(\mathit{\mu }+o_S+\mathit{c}\right)},$

$a_{13}=\mathit{S}{\beta }_3\mathrm{\Lambda }\frac{1-\left(\varphi +\mu +o_S\right)}{\left(\mathit{\kappa }+{\gamma }_2+\mathit{c}\right)\left(\mathit{\mu }+o_S+\mathit{c}\right)}.$

Because the effective reproduction number is the same as the FV⁻¹ spectrum radius:

$R_t=\frac{\mathrm{\Lambda }\left[1-\left(\varphi +\mu +o_S\right)\right]}{\left(o_e+\mathit{\omega }+\mathit{c}\right)\left(\mathit{\mu }+o_S+\mathit{c}\right)}\left[{\beta }_1+{\beta }_2\frac{\mathit{\omega p}}{\lambda +{\gamma }_1+\mathit{c}}+{\beta }_3\frac{\mathit{\omega }\left(1-\mathit{p}\right)}{\kappa +{\gamma }_2+\mathit{c}}\right].$

In this study:

$R_t=\frac{\left[1-\left(\varphi +\mu +o_S\right)\right]}{\left(o_e+\mathit{\omega }\right)\left(\mathit{\mu }+o_S\right)}\left[{\beta }_1+{\beta }_2\frac{\mathit{\omega p}}{\lambda +{\gamma }_1}+{\beta }_3\frac{\mathit{\omega }\left(1-\mathit{p}\right)}{\kappa +{\gamma }_2}\right].$

The stochastic model

We have extended our deterministic model (1.1) to a corresponding stochastic model, recognizing that the spread of infectious diseases is subject to stochastic perturbations. The study of the stochastic model is motivated by the observation that dynamics can be significantly impacted by small variations in the parameters. In many instances, infections are confined to small localities or geographic regions, making the setup of a stochastic model appropriate. As a result, important and relevant information can be derived from the stochastic model.

Assuming that noise directly disturbs each population segments in proportion, the stochastic model (1.3) is described as follows:

$\begin{array}{rl}& \mathit{dS}(\mathit{t})=\left[-\left\{\mathit{S}(\mathit{t})[1-(\mathit{\varphi }+\mathit{\mu }+o_S)]\right\}({\beta }_1\mathit{E}(\mathit{t})+{\beta }_2\mathit{I}(\mathit{t})+{\beta }_3\mathit{A}(\mathit{t}))-\mathit{\mu S}(\mathit{t})-o_S\mathit{S}(\mathit{t})\right]\mathit{dt}+{\sigma }_1\mathit{S}(\mathit{t})\mathit{d}{B}_1(\mathit{t})\\ & \mathit{dE}(\mathit{t})=\left[\left\{\mathit{S}(\mathit{t})[1-(\mathit{\varphi }+\mathit{\mu }+o_S)]\right\}({\beta }_1\mathit{E}(\mathit{t})+{\beta }_2\mathit{I}(\mathit{t})+{\beta }_3\mathit{A}(\mathit{t}))-o_{\mathrm{e}}\mathit{E}(\mathit{t})-\mathit{\omega E}(\mathit{t})\right]\mathit{dt}+{\sigma }_2\mathit{E}(\mathit{t})\mathit{d}{B}_2(\mathit{t})\\ & \mathit{dI}(\mathit{t})=\left[\mathit{\omega p}(1-o_e)\mathit{E}(\mathit{t})-\mathit{\lambda I}(\mathit{t})-{\gamma }_1\mathit{I}(\mathit{t})\right]\mathit{dt}+{\sigma }_3\mathit{I}(\mathit{t})\mathit{d}{B}_3(\mathit{t})\\ & \mathit{dA}(\mathit{t})=\left[\mathit{\omega }(1-\mathit{p})(1-o_e)\mathit{E}(\mathit{t})-\mathit{\kappa A}(\mathit{t})-{\gamma }_2\mathit{A}(\mathit{t})\right]\mathit{dt}+{\sigma }_4\mathit{A}(\mathit{t})\mathit{d}{B}_4(\mathit{t})\\ & \mathit{d}{Q}_s(\mathit{t})=\left(o_S\mathit{S}(\mathit{t})\right)\mathit{dt}+{\sigma }_5{Q}_s(\mathit{t})\mathit{d}{B}_5(\mathit{t})\\ & \mathit{d}{Q}_e(\mathit{t})=\left(o_{\mathrm{e}}\mathit{E}(\mathit{t})\right)\mathit{dt}+{\sigma }_6{Q}_e(\mathit{t})\mathit{d}{B}_6(\mathit{t})\\ & \mathit{dV}(\mathit{t})=\left(\mathit{\mu S}(\mathit{t})\right)\mathit{dt}+{\sigma }_7\mathit{V}(\mathit{t})\mathit{d}{B}_7(\mathit{t})\\ & \mathit{dH}(\mathit{t})=\left(\mathit{\lambda I}(\mathit{t})+\mathit{\kappa A}(\mathit{t})-{\gamma }_3\mathit{H}(\mathit{t})\right)\mathit{dt}+{\sigma }_8\mathit{H}(\mathit{t})\mathit{d}{B}_8(\mathit{t})\\ & \mathit{dR}(\mathit{t})=\left({\gamma }_1\mathit{I}(\mathit{t})+{\gamma }_2\mathit{A}(\mathit{t})+{\gamma }_3\mathit{H}(\mathit{t})\right)\mathit{dt}+{\sigma }_9\mathit{R}(\mathit{t})\mathit{d}{B}_9(\mathit{t})\end{array}$

$B_i\left(\mathit{t}\right)\left(\mathit{i}=1,2,3,4,5,6,7,8,9\right)$ is an independent standard Wiener process with $B_i\left(0\right)=0$. And for any $\mathit{t}\ge 0,{\sigma }_i\left(\mathit{t}\right)\left(\mathit{i}=1,2,3,4,5,6,7,8,9\right)$ is a continuous bounded function, where ${\sigma }_i^2\left(\mathit{t}\right)\left(\mathit{i}=1,2,3,4,5,6,7,8,9\right)$ is the strength of the Wiener process.

Estimation

We use Matlab software to estimate the parameters of the model. The fourth-order Runge-Kutta-Fehlberg Method, with tolerance set at 0.001, was used to solve the differential equation, and curve fitting. The parameter estimates were obtained by fitting the data with Matlab software. The routines applied for estimation is as flows: First, the time period data were selected for analysis. Then we used Matlab software to write the program of the ordinary differential equation, and finally the parameters are estimated with the actual data to get the corresponding results. The corresponding stochastic model adds random terms in the deterministic model. For estimation of stochastic models, we write the program of the corresponding stochastic differential equation, and finally estimates the parameters with the actual data to get the corresponding results.

We selected five periods for analyzes and estimation. The periods selected are the severe COVID-19 periods between 2020 and 2022, which are January 2, 2021 to February 14, 2021, July 4, 2021 to September 8, 2021, September 8, 2021 to October 14, 2021, October 14, 2021 to November 24, 2021 and November 24, 2021 to January 21, 2022.

Sensitivity analysis

The forward sensitivity index of normalization of function f with variable x is ${\xi }_{\mathrm{x}}^{\mathrm{u}}=\frac{\mathrm{x}}{\mathrm{u}}\times \frac{\mathrm{\partial }\mathrm{u}}{\mathrm{\partial }\mathrm{x}}$. According to the formula R_t, R_t is affected by four parameters $\text{bbeta}$, $\mathit{\mu }$, $\mathrm{m}$ and $\mathrm{w}$, among which parameters $\mathit{\mu }$, $\mathrm{m}$ and $\mathrm{w}$ are uncontrollable variables with fixed parameter values in the process of disease transmission, so we only need to analyze the sensitivity changes of R_t with $\text{bbeta}$.

We find

${\mathit{\xi }}_{{\mathit{\beta }}_1}^{{\mathrm{R}}_{\mathrm{t}}}=\frac{{\text{bbeta}}_1}{{\mathrm{R}}_{\mathrm{t}}}\times \frac{\mathrm{\partial }{\mathrm{R}}_{\mathrm{t}}}{\mathrm{\partial }{\text{bbeta}}_1}=\frac{{\text{bbeta}}_1\left(\mathit{\lambda }+{\mathit{\gamma }}_1\right)\left(\mathit{\kappa }+{\mathit{\gamma }}_2\right)}{{\text{bbeta}}_1\left(\mathit{\lambda }+{\mathit{\gamma }}_1\right)\left(\mathit{\kappa }+{\mathit{\gamma }}_2\right)+{\text{bbeta}}_2\mathit{\omega }\mathrm{p}\left(\mathit{\kappa }+{\mathit{\gamma }}_2\right)+{\text{bbeta}}_3\mathit{\omega }\left(1-\mathrm{p}\right)\left(\mathit{\lambda }+{\mathit{\gamma }}_1\right)}$

${\mathit{\xi }}_{{\text{bbeta}}_2}^{{\mathrm{R}}_{\mathrm{t}}}=\frac{{\text{bbeta}}_2}{{\mathrm{R}}_{\mathrm{t}}}\times \frac{\mathrm{\partial }{\mathrm{R}}_{\mathrm{t}}}{\mathrm{\partial }{\text{bbeta}}_2}=\frac{{\text{bbeta}}_2\mathit{\omega }\mathrm{p}\left(\mathit{\kappa }+{\mathit{\gamma }}_2\right)}{{\text{bbeta}}_1\left(\mathit{\lambda }+{\mathit{\gamma }}_1\right)\left(\mathit{\kappa }+{\mathit{\gamma }}_2\right)+{\text{bbeta}}_2\mathit{\omega }\mathrm{p}\left(\mathit{\kappa }+{\mathit{\gamma }}_2\right)+{\text{bbeta}}_3\mathit{\omega }\left(1-\mathrm{p}\right)\left(\mathit{\lambda }+{\mathit{\gamma }}_1\right)}$

${\mathit{\xi }}_{{\text{bbeta}}_3}^{{\mathrm{R}}_{\mathrm{t}}}=\frac{{\text{bbeta}}_3}{{\mathrm{R}}_{\mathrm{t}}}\times \frac{\mathrm{\partial }{\mathrm{R}}_{\mathrm{t}}}{\mathrm{\partial }{\text{bbeta}}_3}=\frac{{\text{bbeta}}_3\mathit{\omega }\left(1-\mathrm{p}\right)\left(\mathit{\lambda }+{\mathit{\gamma }}_1\right)}{{\text{bbeta}}_1\left(\mathit{\lambda }+{\mathit{\gamma }}_1\right)\left(\mathit{\kappa }+{\mathit{\gamma }}_2\right)+{\text{bbeta}}_2\mathit{\omega }\mathrm{p}\left(\mathit{\kappa }+{\mathit{\gamma }}_2\right)+{\text{bbeta}}_3\mathit{\omega }\left(1-\mathrm{p}\right)\left(\mathit{\lambda }+{\mathit{\gamma }}_1\right)}$

${\mathit{\xi }}_{\mathit{\mu }}^{{\mathrm{R}}_{\mathrm{t}}}=\frac{\mathit{\mu }}{{\mathrm{R}}_{\mathrm{t}}}\times \frac{\mathrm{\partial }{\mathrm{R}}_{\mathrm{t}}}{\mathrm{\partial }\mathit{\mu }}=\frac{\mathit{\mu }\left(\mathit{\varphi }-1\right)}{\left[1-\left(\mathit{\varphi }+\mathit{\mu }+{\mathrm{o}}_{\mathrm{S}}\right)\right]\left(\mathit{\mu }+{\mathrm{o}}_{\mathrm{S}}\right)}$

${\mathit{\xi }}_{{\mathrm{o}}_{\mathrm{s}}}^{{\mathrm{R}}_{\mathrm{t}}}=\frac{{\mathrm{o}}_{\mathrm{s}}}{{\mathrm{R}}_{\mathrm{t}}}\times \frac{\mathrm{\partial }{\mathrm{R}}_{\mathrm{t}}}{\mathrm{\partial }{\mathrm{o}}_{\mathrm{s}}}=\frac{\mathit{\mu }\left(\mathit{\varphi }-1\right)}{\left[1-\left(\mathit{\varphi }+\mathit{\mu }+{\mathrm{o}}_{\mathrm{s}}\right)\right]\left(\mathit{\mu }+{\mathrm{o}}_{\mathrm{s}}\right)}$

${\mathit{\xi }}_{{\mathrm{o}}_{\mathrm{e}}}^{{\mathrm{R}}_{\mathrm{t}}}=\frac{{\mathrm{o}}_{\mathrm{e}}}{{\mathrm{R}}_{\mathrm{t}}}\times \frac{\mathrm{\partial }{\mathrm{R}}_{\mathrm{t}}}{\mathrm{\partial }{\mathrm{o}}_{\mathrm{e}}}=-\frac{{\mathrm{o}}_{\mathrm{e}}}{\left({\mathrm{o}}_{\mathrm{e}}+\mathit{\omega }\right)}$

${\mathit{\xi }}_{\mathit{\lambda }}^{{\mathrm{R}}_{\mathrm{t}}}=\frac{\mathit{\lambda }}{{\mathrm{R}}_{\mathrm{t}}}\times \frac{\mathrm{\partial }{\mathrm{R}}_{\mathrm{t}}}{\mathrm{\partial }\mathit{\lambda }}=-\frac{\mathrm{p\omega }{\text{bbeta}}_2}{\left[{\text{bbeta}}_1+{\text{bbeta}}_2\frac{\mathit{\omega }\mathrm{p}}{\mathit{\lambda }+{\mathit{\gamma }}_1}+{\text{bbeta}}_3\frac{\mathit{\omega }\left(1-\mathrm{p}\right)}{\mathit{\kappa }+{\mathit{\gamma }}_2}\right]{(\mathit{\lambda }+{\mathit{\gamma }}_1)}^2}$

${\mathit{\xi }}_{\mathit{\kappa }}^{{\mathrm{R}}_{\mathrm{t}}}=\frac{\mathit{\kappa }}{{\mathrm{R}}_{\mathrm{t}}}\times \frac{\mathrm{\partial }{\mathrm{R}}_{\mathrm{t}}}{\mathrm{\partial }\mathit{\kappa }}=-\frac{\left(1-\mathrm{p}\right)\mathit{\omega }{\text{bbeta}}_3}{\left[{\text{bbeta}}_1+{\text{bbeta}}_2\frac{\mathit{\omega }\mathrm{p}}{\mathit{\lambda }+{\mathit{\gamma }}_1}+{\text{bbeta}}_3\frac{\mathit{\omega }\left(1-\mathrm{p}\right)}{\mathit{\kappa }+{\mathit{\gamma }}_2}\right]{(\mathit{\kappa }+{\mathit{\gamma }}_2)}^2}$

Result

Since March 4, 2020, there has been a significant decrease in the number of active COVID-19 cases in China, with no significant new outbreaks reported until December 2020. Subsequently, there has been a gradual emergence of COVID-19 outbreaks in China. Serious outbreaks of COVID-19 were observed in January, July, September, November, and December of 2021 (refer to Figure 2).

Figure 2. Trend of active cases of COVID-19 in China from March 4, 2020 to January 21, 2022.

Deterministic model

The output of the fitted model and the number of active COVID-19 cases in five different major epidemics are shown in Figure 3. Tables 3 and 4 shows the parameter estimates and effective reproduction number of the five time periods.

Figure 3. Plots of the output of the fitted model and the observed active COVID-19 case. (a) January 2, 2021 to February 14, 2021; (b) July 4, 2021 to September 8, 2021; (c) September 8, 2021 to October 14, 2021; (d) October 14, 2021 to November 24, 2021; e: November 24, 2021 to January 21, 2022.

Table 3. Mean values of estimated parameters ${\beta }_1$, ${\beta }_2$, ${\beta }_3$, $\mathit{\mu }$, $o_s$, $o_e$ and $\mathit{\kappa }$. These parameters are explained in Table 2: ${\beta }_1$, ${\beta }_2$, and ${\beta }_3$ refer to transmission rate of susceptible with exposed, symptomatic and asymptomatic human, respectively. $\mathit{\mu }$ refers to rate of effective vaccination. ${\mathrm{o}}_{\mathrm{s}}$ and ${\mathrm{o}}_{\mathrm{e}}$refer to quarantine rate of susceptible and exposed individuals, respectively. $\mathit{\kappa }$ refers to rate of hospitalized of asymptomatic individual.

	${\beta }_1$	${\beta }_2$	${\beta }_3$	$\mathit{\mu }$	$o_s$	$o_e$	$\mathit{\kappa }$
January 2, 2021 to February 14, 2021	$4.572\times {10}^{-11}$	$7.4754\times {10}^{-10}$	$7.4736\times {10}^{-11}$	0.0327	0.00193	0.0024	0.00365
July 4, 2021 to September 8, 2021	$2.1611\times {10}^{-10}$	$3.8818\times {10}^{-11}$	$3.3442\times {10}^{-11}$	0.0191	$4.1\times {10}^{-4}$	0.00761	0.0157
September 8, 2021 to October 14, 2021	$2.5373\times {10}^{-10}$	$1.0211\times {10}^{-10}$	$1.0237\times {10}^{-11}$	0.00198	0.01945	$1.79\times {10}^{-4}$	0.00409
October 14, 2021 to November 24, 2021	$1.3825\times {10}^{-9}$	$8.0474\times {10}^{-10}$	$4.8803\times {10}^{-11}$	0.00795	0.00758	0.0027	0.06867
November 24, 2021 to January 21, 2022	$1.3313\times {10}^{-9}$	$4.0647\times {10}^{-11}$	$3.3456\times {10}^{-11}$	0.0234	0.0063	0.0203	0.00585

Table 4. Mean values of estimated effective reproduction number.

	Rt
January 2, 2021 to February 14, 2021	3.5
July 4, 2021 to September 8, 2021	2.4
September 8, 2021 to October 14, 2021	1.6
October 14, 2021 to November 24, 2021	1.9
November 24, 2021 to January 21, 2022	2.6

The highest effective reproduction number of the epidemic was observed from January 2, 2021 to February 14, 2021, followed by the period from November 24, 2021 to January 21, 2022, while the lowest was noted from September 8, 2021 to October 14, 2021. This observation may be related to the implementation of prevention and control measures, population mobility, and the scope of the epidemic. Patients in the incubation period are more infectious than asymptomatic patients and symptomatic patients. Furthermore, a higher vaccination rate is observed during severe outbreaks. Figure 3 illustrates the fitting of the model with real COVID-19 data for all datasets. It is important to emphasize that the results of the model are effective only for a short time and depend entirely on the current epidemic situation. Predictions for the long term are neither feasible nor appropriate, given the dynamic and evolving nature of disease and interventions.

Sensitivity analysis

According to the sensitivity index analysis, R_t was most sensitive to parameters ${\text{bbeta}}_1$, $\mathit{\mu },\mathit{\hspace{0.17em}}{\mathrm{o}}_{\mathrm{s}}$, and positively correlated with ${\text{bbeta}}_1$, ${\text{bbeta}}_2$, ${\text{bbeta}}_3$, negatively correlated with $\mathit{\mu },\mathit{\hspace{0.17em}}{\mathrm{o}}_{\mathrm{s}}$. In theory, in order to control the spread of COVID-19, ${\text{bbeta}}_1$ ${\text{bbeta}}_2,\mathit{\hspace{0.17em}}{\text{bbeta}}_3$ can be reduced, and the effect of reducing ${\text{bbeta}}_1$ may be better. At the same time, increasing the effective vaccination rate of susceptible individuals and improving the quarantine rate of susceptible individuals can better control the spread of COVID-19. The effect of quarantine rates of exposed individuals and hospitalization rates of symptomatic and asymptomatic individuals on the transmission of COVID-19 is not significant (Figure 4).

Figure 4. R_t sensitivity analysis. (a) January 2, 2021 to February 14, 2021; (b) July 4, 2021 to September 8, 2021; (c) September 8, 2021 to October 14, 2021; (d) October 14, 2021 to November 24, 2021; (e) November 24, 2021 to January 21, 2022. Sensitivity index of R_t with respect to parameters ${\text{bbeta}}_1$, ${\text{bbeta}}_2$, ${\text{bbeta}}_3$, $\mathit{\mu },\mathit{\hspace{0.17em}}{\mathrm{o}}_{\mathrm{s}}$, ${\mathrm{o}}_{\mathrm{e}}$ $\mathit{\lambda },$ and $\mathit{\kappa },$ respectively. These parameters are explained in Table 2: ${\text{bbeta}}_1$, ${\text{bbeta}}_2$, and ${\text{bbeta}}_3$ refer to transmission rate of susceptible with exposed, symptomatic and asymptomatic human, respectively. $\mathit{\mu }$ refers to rate of effective vaccination. ${\mathrm{o}}_{\mathrm{s}}$ and ${\mathrm{o}}_{\mathrm{e}}$ refer to quarantine rate of susceptible and exposed individuals, respectively. $\mathit{\lambda }$ refers to rate of hospitalized of symptomatic individuals. $\mathit{\kappa }$ refers to rate of hospitalized of asymptomatic individual.

Stochastic model

We applied the stochastic differential equation (SDE) to simulate each of the five outbreaks 100 times. We set the Wiener process intensity for each equation to ${\mathit{\sigma }}_1=$ 0.00025, ${\mathit{\sigma }}_2={\mathit{\sigma }}_4={\mathit{\sigma }}_5={\mathit{\sigma }}_7={\mathit{\sigma }}_8={\mathit{\sigma }}_9=$ 0.025, ${\mathit{\sigma }}_3={\mathit{\sigma }}_6=$ 0.035. Initially, we compared the average of the 100 random model fits to the actual data and the results of the deterministic model. The figure demonstrates that the mean of the stochastic model is similar to the fitting results of the deterministic model, except for the period from July 4, 2021, to September 8, 2021 (Figure 5). To further investigate the effect of randomness on the infected individuals, we graphed 100 simulation runs for both symptomatic (I) and asymptomatic (A) infected individuals. We observed that asymptomatic individuals might be more numerous than symptomatic patients due to random interference (Supplementary 1 and Supplementary 2).

Figure 5. Output of the stochastic model and deterministic model and observed active COVID-19 cases. (a) January 2, 2021 to February 14, 2021; (b) July 4, 2021 to September 8, 2021; (c) September 8, 2021 to October 14, 2021; (d) October 14, 2021 to November 24, 2021; (e) November 24, 2021 to January 21, 2022.

In order to enhance the comprehension of the outcomes of the stochastic model, we generated a distribution plot of the total number of individuals affected by 100 simulations of the five outbreaks. The order of the figures is consistent with Figure 3. Our findings reveal that the likelihood of a relatively minor outbreak is considerably low, leading to only a few individuals getting infected. Conversely, the probability of a significant epidemic is noticeably more than 60%. Furthermore, even two minor outbreaks in September and October can result in significant harm (Supplementary 3).

As an illustration, we examined the preventive and control measures for the epidemic that occurred from January 2, 2021 to February 14, 2021, by utilizing stochastic models. We evaluated the impact of measures such as wearing masks, quarantine, vaccination, and hospitalization, while varying intervention rates from 0 to 0.8. These images from left to right correspond to intervention measures of 0, 0.2, 0.4, 0.6, and 0.8 in size, respectively. Our analysis revealed that as the rates of quarantine for susceptible and exposed individuals, effective vaccination, effective mask-wearing, and hospitalization increase, there is a gradual decrease in the likelihood of a major outbreak. Notably, the hospitalization rate of asymptomatic individuals did not have a significant effect on the epidemic. We further observed that effective mask-wearing, quarantine of exposed individuals, and effective vaccination can shift the distribution of infected individuals toward a positive skewness. (See Supplementary 4–9 for detailed results.)

Discussion

Since the emergence of COVID-19 in Wuhan, China in December 2019, rigorous control measures have been implemented resulting in successful containment of the large-scale spread of the disease by April 2020. Despite this achievement, the intricate global landscape has imposed significant challenges on the efforts to prevent and manage the epidemic, leading to the reemergence of the COVID-19 epidemic in the country.

Modeling is an important method to describe the transmission of COVID-19 in China and assess the effectiveness of prevention and control measures. The SEIAQVHR model we constructed is methodologically robust and is the extensions the classical SEIR model. We proposed deterministic and stochastic models to model the incidence trends during five COVID-19 outbreaks in China. Unlike other models, in our deterministic model we considered the impact of asymptomatic infections, quarantine, vaccination and hospitalization on the spread of the epidemic. Our stochastic model is an extension of the deterministic model, where we considered random interference in the transmission of the disease to make the results more realistic. Our SEIAQVHR model was proposed to assessed the impact of various preventive and control measures on China’s efforts to combat the COVID-19 pandemic and provide a reference for the prevention and control of emerging infectious diseases that may occur in the future.

In this study, we analyzed the course of COVID-19 in China and developed a deterministic compartment model to describe disease transmission in the population. Our model included measures such as mask-wearing, hospitalization of symptomatic and asymptomatic patients, quarantine for susceptible and exposed individuals, and vaccination of susceptible individuals. The result of the fitted model indicate that the phase of the epidemic differed from the actual situation at the beginning. The starting point of the time period we selected is at the point where the previous epidemic tends to end and the next one does not start. At this point, the number of existing patients has a downward trend. Since the model fits the data according to the ordinary differential equation, the model will have a decrease at first, but the number of active cases is obviously increasing. It is inevitable that the model fitting results differ from the actual data. The parameters included in the model, such as recovery rate, isolation rate, etc., are fixed, but the actual situation is constantly changing, and the corresponding parameters are also in the process of changing. Therefore, the initial phase of the epidemic differed from the actual situation and the model wasn’t able to fit it well. After that, the actual situation tends to stabilize and the model fitted the actual data better. We will further improve this in subsequent studies. We will consider setting the parameters as changing parameters, which are changing over time.

Our estimated parameters showed that the transmission rate was highest among exposed individuals, followed by symptomatic individuals, and lowest among asymptomatic individuals. This suggests that exposed individuals, who may not display clinical symptoms, could spread the disease to susceptible individuals for an extended period of time. Therefore, we recommend focusing on this population, as ignoring them could lead to a severe COVID-19 epidemic. Our R_t calculation indicated that the outbreak was most severe between 2 January 2021 and 4 February 2021, similar to the actual situation. During this period, the effective vaccination rate was also found to be higher than other periods, possibly due to the promotion of vaccines during the severe epidemic. These results highlight the importance of controlling disease transmission rates to prevent rapid spread among the population.

The results of the stochastic model indicate that it is a useful tool for analyzing disease spread and providing additional information beyond what is obtained from deterministic models. Our analysis of five outbreaks revealed that most of the outbreaks had a high probability (>60%) of developing into serious outbreaks, with a potential total patient count exceeding 20,000, highlighting the highly infectious nature of COVID-19. This finding warrants attention, as it suggests that even with strict prevention and control measures, large-scale outbreaks are still possible.

We conducted a comprehensive analysis of parameters associated with interventions such as face masks, quarantine, vaccination, and hospitalization, using a stochastic model. Our results suggest that higher intervention efforts are required to control disease outbreaks within a shorter period in China. The reduction in the probability of a major outbreak mainly depends on the rate of effective vaccination, which directly reduces the number of susceptible individuals and thus reduces the spread of the disease in the population. Quarantine of susceptible and exposed individuals is a crucial measure to prevent the spread of the disease. We found that the higher the quarantine rate, the lower the probability of susceptible individuals coming into contact with infectious individuals (exposed, symptomatic, and asymptomatic individuals), thus reducing the spread of the disease. Hospitalization of symptomatic individuals has a significant impact on controlling the spread of the disease, but the hospitalization rate of asymptomatic individuals has little effect. Wearing masks is also an effective measure to control the spread of COVID-19. Our findings suggest that increasing the effective wearing rate of masks can be as effective in reducing the spread of COVID-19 as vaccination. However, no single intervention strategy can completely eliminate COVID-19, as even an increase of 80% in any one intervention would still result in a large number of patients. Nonetheless, a combination of multiple interventions can theoretically eliminate COVID-19.²⁰

While there are signs that COVID-19 may transit into an endemic phase, it is still unclear when this stage will be reached. However, our modeling approach provides some evidence how to manage future outbreaks of COVID-19 or other infectious diseases with similar characteristics. We found vaccination and non-pharmacological interventions like mask-wearing as most effective in curbing an epidemic outbreak. While vaccination may not completely stop transmission it can reduce the probability of transmission and wearing masks adds to this effect thereby stopping further spread of the infection in the community. These findings can not only help governments prepare for future outbreaks, but also provide an important reference for optimizing intervention plans for this and similar diseases.

Authors’ contribution

Conceptualization and Writing – original draft: Wenting Zha and Han Ni; Formal analysis and Software: Wenting Zha, Wentao Kuang and Yuxi He; Investigation and Data curation: Xuewen Yang and Jin Zhao; Validation and Methodology: Nan Ni, Nan Zhou and Yuxi He; Visualization and Project administration: Wentao Kuang, Liuyi Fu and Haoyun Dai; Writing – review & editing: Yuan Lv, Nan Zhou and Xuewen Yang; Funding acquisition: Yuan Lv and Wenting Zha. All authors have read and agreed to the published version of the manuscript.

Code

The research involves a significant amount of computational code, and we will only provide some of the key code snippets here:

s=randn(1,44);

for idx=1:length(T)-1

S(idx+1)=S(idx)-S(idx)*(1-(b+c+ds))*(a1*E(idx)+a2*I(idx)

 +a3*A(idx))-b*S(idx)-ds*S(idx)+0.00025*s(idx)*S(idx);

E(idx+1)=E(idx)+S(idx)*(1-(b+c+ds))*(a1*E(idx)+a2*I(idx)

 +a3*A(idx))-de*E(idx)-f*E(idx)+0.025*s(idx)*E(idx);

I(idx +1)=I(idx)+f*g*E(idx)-h*I(idx)-j*I(idx)+0.035*s(idx)*I(idx);

A(idx+1)=A(idx)+f*(1-g)*E(idx)-i*A(idx)-k*A(idx)+0.025*s(idx)

 *A(idx);

Qs(idx +1)=Qs(idx)+ds*S(idx)+0.025*s(idx)*Qs(idx);

Qe(idx +1)=Qe(idx)+de*E(idx)+0.035*s(idx)*Qe(idx);

V(idx +1)=V(idx)+b*S(idx)+0.025*s(idx)*V(idx);

H(idx +1)=H(idx)+h*I(idx)+i*A(idx)-l*H(idx)+0.025*s(idx)*H(idx);

R(idx +1)=R(idx)+j*I(idx)+k*A(idx)+l*H(idx)+0.025*s(idx)*R(idx);

Q(idx +1)=Qs(idx)+Qe(idx);

Disclosure statement

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

Data availability statement

Data sources for this study: https://ourworldindata.org/; https://coronavirus.jhu.edu/data.

Supplementary material

Supplemental data for this article can be accessed on the publisher’s website at https://doi.org/10.1080/21645515.2024.2338953.

Additional information

Funding

This research was funded by [Hunan Provincial Department of Education] grant number [JG2018B041], [Hunan Natural Fund] grant number [2020JJ4059], [Changsha City Science and Technology Plan Project] grant number [kq2001025] and [Research project of Hunan Provincial Health Commission] grant number [202212054651].