Mathematical modeling of mutated COVID-19 transmission with quarantine, isolation and vaccination

1.   Introduction

Corona Virus Disease 2019 (COVID-19), caused by SARS-CoV-2, has been seriously impacting the world since the end of 2019 ^[1]. As of April 13, 2022, there have been over 499 million confirmed cases of COVID-19, including over 6 million deaths all over the world ^[2].

Since the outbreak of COVID-19, many scholars have used different methods to study its transmission ^[3], such as neural network ^[4], time series analysis ^[5] and kinetic modeling ^[6].

Dynamics models have played a key role in studying the development of things and assessing the impact of different intervention strategies ^[7,8]. In particular, during the COVID-19 outbreak, multiple dynamics models have been developed to investigate its spread and control.

Before the availability of effective vaccines, non-pharmaceutical interventions strategies (NPIs) were taken into account controlling the pandemic transmission, such as lockdowns, the use of face-masks in public, quarantine of close contacts, isolation of confirmed cases, among others ^[9]. Ngonghala et al. ^[10] proposed a detailed mathematical model of COVID-19 that combined quarantine and isolation to assess the impact of NPIs on combating and mitigating the burden of COVID-19. Their results indicated that the early implementation (and enhancement of effectiveness) of these intervention measures is obviously critically-important to curtail COVID-19 transmission. NPIs could mitigate the spread of the virus, but it is not the permanent solution to the problem ^[11].

Since January 2021, more than ten types of vaccines against SARS-CoV-2 have been fully or limitedly approved and used clinically ^[12]. And many countries are starting vaccination campaigns, which has improved the global epidemic to some extent ^[13]. Iboi et al. ^[14] developed a mathematical model for assessing the impact of an anti-COVID-19 vaccine on the spread of COVID-19 in the United States. This study showed that the prospect of eliminating COVID-19 in the US using the imperfect hypothetical vaccine is promising. Moreover, the integration of vaccination and physical distancing interventions, can result in better pandemic prevention ^[15]. Zou et al. ^[16] established a multiple patch coupled model to explore the impact of vaccination and NPIs on the control of COVID-19 and the results indicated that effective vaccination has positive impact on prevention of pandemic transmission and the joint implementation of vaccination and quarantine ensures the controllability of the epidemic.

However, new variants of COVID-19 with high transmission rates are spreading around the world and vaccines are less effective in preventing infection than they were for earlier virus variants, which have brought new challenges to the elimination of pandemic ^[17,18]. Recently, Li et al. ^[19] developed a mathematical model of COVID-19 transmission with imperfect vaccination to explore effective and reasonable plans to prevent the spread of Delta variant. Their results found that the optimal control measure is to dynamically adjust three control measures, namely, vaccination, isolation and nucleic acid testing, to achieve the lowest number of infections at the lowest cost. Truszkowska et al. ^[20] developed an extremely detailed mathematical model, to study the combined effect of booster shot administration and testing practices in this stage of the pandemic.

Motivated by the implement of NPIs, the promotion of vaccination strategy and the emerging of variants of COVID-19, we consider a mathematical model for mutated COVID-19 with quarantine, isolation and vaccination, make some theoretical analysis of the model, verify the applicability of the model by the COVID-19 pandemic mainly caused by Delta variant in South Korea, and finally, assess the impact of quarantine and vaccination on eliminating mutated COVID-19 transmission.

The organization of this paper is as follows. In Section 2, the basic model formulation is discussed. In Section 3 we derive the disease free equilibrium and basic reproduction number $\mathscr{R}_{0}$. The local asymptotic stability and global asymptotic stability of the model at the disease free equilibrium are proved. In Section 4, the numerical simulation is carried out to verify the correctness of the proof. In Section 5, according to the epidemic data of COVID-19 in South Korea from August 22 to October 20, 2021, unknown parameters are obtained by two-stage fitting, and then the epidemic is predicted, which is compared with the real data to verify the applicability of the model. In Section 6, the impact of quarantine and vaccination on the control of COVID-19 are assessed by numerical simulation. Finally, in Section 7, we give some discussions.

2.   Model formulation

We develop an epidemic model, which incorporates the main interventions being implemented to curtail COVID-19 transmission (such as quarantine, isolation and vaccination). We stratify the population as susceptible $(S_{u})$, vaccinated $(V_{u})$, exposed $(E_{u})$, asymptomatically-infectious $(I_{a})$, symptomatically-infectious $(I_{u})$ and recovered $(R)$ compartments, and further stratify the population to include quarantined susceptible $(S_{q})$, quarantined vaccinated $(V_{q})$, quarantined exposed $(E_{q})$ and hospitalized $(I_{h})$ compartments, so that the total population of the model $N(t)$ is divided into 10 compartments, namely

For developing the mathematical model, the basic assumptions are as follows:

(ⅰ) The recruitment rate (either by birth or by immigration) of the population is given by a constant rate $\varLambda$ and they are all susceptible.

(ⅱ) Susceptible individuals move to exposed compartment by effective contacts with exposed, symptomatically-infectious or asymptomatically-infectious individuals and after latency period, they become infectious and move to infectious compartment.

(ⅲ) After recovery the individuals have immunity to COVID-19 in a short time.

(ⅳ) A part of vaccinated individuals will lose their immunity and rejoin the susceptible compartment.

(ⅴ) After the immunity period, the vaccinated individuals become susceptible again.

(ⅵ) Hospitalization is completely effective.

(ⅶ) The exposed individuals and the asymptomatically-infectious individuals are also infectious, but the infectivity is weaker than that of symptomatically-infectious individuals, measuring by coefficients $\eta_{u}$ and $\eta_{a}$ respectively$(0 < \eta_{u}, \eta_{a} < 1)$.

(ⅷ) The contact rate $c$ is the same for the individuals in symptomatically-infectious $I_{u}$, non-quarantined exposed $E_{u}$ and the asymptomatically-infectious $I_{a}$ compartments.

(ⅸ) The latency period is defined as the days from exposure to the onset of illness.

(ⅹ) The asymptomatically-infectious compartment $I_{a}$ includes those with mild symptoms or no clinical symptoms of the COVID-19 at the end of the latency period.

(ⅹⅰ) Considering that the number of people that an infected person can contact per unit time is finite, a more realistic standard incidence is used in the model.

(ⅹⅱ) The vaccinated individuals will not lose immunity during the quarantine period.

(ⅹⅲ) If the quarantine individuals have symptoms, they will be sent to the hospital.

The assumptions above lead to the flow diagram of the model depicted in Figure 1.

It is worth noting that hospitalization indicates isolation at the hospital in this study. With contact tracing, a proportion, $q$, of individuals exposed to the virus is quarantined. The quarantined individuals, if infected, move to the compartment $E_{q}$ and and the remaining individuals can either stay in compartment $V_{q}$ or $S_{q}$, depending on whether they are vaccinated or not. For those individuals who are missed from contact tracing, the situation is similar.

The model is given by the following deterministic system of nonlinear differential equations.

where $\lambda(t) =$$\dfrac{c(I_{u}+\eta_{u}E_{u}+\theta\eta_{u}E_{q}+\eta_{a}I_{a})}{S_{u}+V_{u}+E_{u}+I_{u}+I_{a}+R+\theta(S_{q}+E_{q}+V_{q})}$.

The initial conditions for system (2.1) are as follows:

The state variables and parameters of the model are described in Tables 1 and 2, respectively.

Let $H(t) = (S_{u}(t), S_{q}(t), V_{u}(t), V_{q}(t), E_{u}(t), E_{q}(t), I_{u}(t), I_{h}(t), I_{a}(t)$ $, R(t)).$

Lemma 1. Consider the initial data $H_{i}(0)\geq0, i = 1, 2, ..., 10.$ Then, for any $t > 0$, we have the non-negative solution for the system (2.1).

Proof. Consider $t_{1} = sup\{t > 0:H_{i}(t) > 0, i = 1, 2, ..., 10\}$. So, $t_{1} > 0$. The following result is given using the first equation of the model (2.1),

Then, it can be written as,

Thus,

so that

For the remaining equations, we take the same steps to show $H_{i}(t) > 0$ $(i = 2, 3, ..., 10)$ for every $t > 0$.

Let the closed and biologically feasible region be $\Omega$, shown by $\Omega = \biggl\{(S_{u}, S_{q}, V_{u}, V_{q}, E_{u}, E_{q}, I_{u}, I_{h}, I_{a}, R)\in R^{10}_{+}, 0\leq S_{u}, S_{q}, V_{u}, V_{q}, E_{u}, E_{q}, I_{u}, I_{h}, I_{a}, R, N\leq\dfrac{\varLambda}{d}\biggr\}$.

Next, we present the following results for this feasible region.

Lemma 2. The region $\Omega$ is a positively invariant and attracting region for system $(2.1)$ with the non-negative initial conditions (2.2).

Proof. Summing up the ten equations in system (2.1) and considering that all parameters of the model are non-negative, we get

Now integrating both sides of the above inequality and using the comparison theorem ^[21], we obtain

Clearly, $0 < N(t)\leq\dfrac{\varLambda}{d}$, as $t\rightarrow \infty$. If $N(0)\leq\dfrac{\varLambda}{d}$, then $N(t)\leq\dfrac{\varLambda}{d}$. Thus, the region $\Omega$ is positive invariant and attracts all the possible solutions of the model (2.1). We will consider the dynamic behavior of model (2.1) on $\Omega$.

3.   Theoretical analysis

3.1.   basic reproduction number $\mathscr{R}_{0}$

Obviously, the model (2.1) always has a disease free equilibrium (DFE) $P_{0}(S_{0}, 0, V_{0}, 0, 0, 0, 0, 0, 0, 0)$, where

The basic reproduction number $\mathscr{R}_{0}$ will be calculated using the next generation matrix method ^[22,23]. It follows that the next generation operator matrices, $F$ and $V$ for the new infection terms and the transition terms, are given, respectively,

where, $M_{1} = c\beta(1-q)\dfrac{S_{0}+\varepsilon_{v}V_{0}}{\Lambda/d}$, $M_{2} = c\beta q\dfrac{S_{0}+\varepsilon_{v}V_{0}}{\Lambda/d}$, $K_{1} = \sigma_{u}+\alpha_{u}+d$, $K_{2} = \sigma_{q}+d$, $K_{3} = \phi_{u}+\gamma_{u}+\delta_{u}+d$, $K_{4} = \gamma_{h}+\delta_{h}+d$, $K_{5} = \sigma_{a}+\gamma_{a}+\delta_{a}+d$. By the definition of the basic reproduction number $\mathscr{R}_{0}$, we get

where,

The quantity $\mathscr{R}_{0}$ is the basic reproduction number of the model (2.1). It measures the average number of the new COVID-19 infections generated by a typical infective individual introduced into a population where basic public health interventions (quarantine, isolation, vaccination, etc.) are implemented.

3.2.   Stability analysis of the model

In this section, we analyse the stability of the model (2.1) at the DFE. First, we show the local stability of system (2.1) at the DFE.

Theorem 1. The DFE of the model $(2.1)$ is locally asymptoticallystable if $\mathscr{R}_{0} < 1$, and unstable if $\mathscr{R}_{0} > 1$.

Proof. Let $\bar{H} = (E_{u}(t), E_{q}(t), I_{u}(t), I_{a}(t), S_{u}(t), S_{q}(t), V_{u}(t), V_{q}(t), I_{h}(t), R(t))$, we can obtain the Jacobian matrix of the system corresponding to $\bar{H}$ at $P_{0}$, and partition the matrix.

It is given by

Using the Laplace theorem, it's easy to find that the eigenvalues of $J(P_{0})$ are determined by the eigenvalues of $J_{1}$ and $J_{4}$.

Consider the eigenvalues of $J_{4}$ first.

It is obvious that $J_{4}$ has always four negative eigenvalues $\lambda_{1} = -d$, $\lambda_{2} = -(\delta_{h}+\gamma_{h}+d)$, $\lambda_{3} = -(p_{1}\omega+d)$ and $\lambda_{4} = -(\omega+d)$, and the other eigenvalues of $J_{4}$ are determined by the equation

Clearly, $\lambda_{5}+\lambda_{6} = -(\xi_{v}+\eta_{v}+2d) < 0$, $\lambda_{5}\cdot\lambda_{6} = \xi_{v}d+\eta_{v}d+d^{2} > 0$, so $\lambda_{5}$ and $\lambda_{6}$ are both negative.

Then consider the eigenvalues of $J_{1}$. Through calculation and integration, we have

where,

If ${\mathscr{R}}_{0} < 1$, we have $\dfrac{\eta_{u}}{K_{1}}M_{1} < 1$ and $\dfrac{\theta\eta_{u}}{K_{2}}M_{2} < 1$, so it can be observed that $a_{1} > 0$.

Because of $\dfrac{a_{4}}{K_{1}K_{2}K_{3}K_{5}} = 1-\mathscr{R}_{0}$, we obtain $a_{4} > 0$ if $\mathscr{R}_{0} < 1$.

Due to

then $a_{3} > 0.$

Meanwhile, we can also obtain that

The first inequality is established due to $\left(\dfrac{\eta_{u}}{K_{1}}+\dfrac{f_{1}\sigma_{u}}{K_{1}K_{3}}\right)M_{1} < 1$, if $\mathscr{R}_{0} < 1$. It is apparent from the above analysis that it satisfies the Routh-Hurwitz stability criterion. Thus, it follows that the DFE $P_{0}$ is locally asymptotically stable if $\mathscr{R}_{0} < 1$. On the other hand, we have $\dfrac{a_{4}}{K_{1}K_{2}K_{3}K_{5}} = 1-\mathscr{R}_{0} < 0$ if $\mathscr{R}_{0} > 1$, then $a_{4} < 0$, which shows that at least one of the eigenvalues of $J_{1}$ is positive. Therefore, the DFE $P_{0}$ is unstable if $\mathscr{R}_{0} > 1$. This completes the proof.

Theorem 2. The DFE of the model $(2.1)$ is globallyasymptotically stable in $\Omega$ if $\mathscr{R}_{0} < 1$.

Proof. Let $\boldsymbol {X} = (E_{u}, E_{q}, I_{u}, I_{h}, I_{a})^{T}$. It can be stated that

with $F$ and $V$ are matrices on calculating $\mathscr{R}_{0}$. Let $\boldsymbol {u} = (\eta_{u}, \theta\eta_{u}, 1, 0, \eta_{a})$. It then follows from the fact $\mathscr{R}_{0} = \rho(FV^{-1}) = \rho(V^{-1}F)$ and direct calculation that $\boldsymbol {u}$ is a left eigenvector of the matrix $V^{-1}F$, i.e., $\boldsymbol {u}V^{-1}F = \mathscr{R}_{0}\boldsymbol {u}$.

Consider a Liapunov function

Differentiating the above equation, we have

If $\mathscr{R}_{0} < 1$, the equality $\dfrac{d\mathscr{L}_{0}}{dt} = 0$ implies that $\boldsymbol {uX} = 0$. This leads to $E_{u} = E_{q} = I_{u} = I_{a} = 0$ by noting that all components of $\boldsymbol {u}$ except for the fourth element are positive.

Clearly, let

then the largest invariant set of the system in $E$ is $M =$$\left\{ E_{u} = E_{q} = I_{u} = I_{a} = 0\right\}$. By LaSalle's Invariant Principle ^[24], we get $\lim\limits_{t\rightarrow \infty}E_{u}\left(t\right) = \lim\limits_{t\rightarrow \infty}E_{q}\left(t\right) = \lim\limits_{t\rightarrow \infty}I_{u}\left(t\right) = \lim\limits_{t\rightarrow \infty}I_{a}\left(t\right) = 0$. Then, it is easy to solve the 8th equation of the model (2.1), and we have

Using L'Hospital's rule,

In a similar way, for the 2nd, 4th and 10th equation of the model (2.1), we have $\lim\limits_{t\rightarrow \infty}S_{q}\left(t\right) = \lim\limits_{t\rightarrow \infty}V_{q}\left(t\right) = \lim\limits_{t\rightarrow \infty}R\left(t\right) = 0$. Hence, we can obtain the limit system of (2.1) as follows.

It is easy to solve that

Clearly, we can get $\lim\limits_{t\rightarrow \infty}S_{u}\left(t\right) = S_{0}$, $\lim\limits_{t\rightarrow \infty}V_{u}\left(t\right) = V_{0}$, so $\left(S_{0}, V_{0}\right)$ is globally attractive for the limit system (3.2). By discussion of linearization systems for (3.2), it can be known that $\left(S_{0}, V_{0}\right)$ is locally asymptotically stable if $\mathscr{R}_{0} < 1$. Combining global attraction of $\left(S_{0}, V_{0}\right)$, we can obtain that $\left(S_{0}, V_{0}\right)$ is globally asymptotically stable for the limit systems (3.2). It follows that $P_{0}(S_{0}, 0, V_{0}, 0, 0, 0, 0, 0, 0, 0)$ is globally attractive for the original system (2.1). In addition, combining Theorem 1, we eventually concluded that $P_{0}$ is globally asymptotically stable if $\mathscr{R}_{0} < 1$.

4.   Numerical simulation of the system (2.1)

In this section, we present the numerical solution of the system (2.1).

We can get $\mathscr{R}_{0} = 0.0147 < 1$ and $P_{0} = (48.7179, 0, 51.2821, 0, 0, 0, 0, 0, 0, 0, 0, 0)$ by using the parameter values given in Table 2. The dynamical behavior for different initial conditions of model (2.1) is presented in Figure 2. It's easy to find that, $S_{u}$ tends to 48.7179, $S_{q}$ tends to 0, $V_{u}$ tends to 51.2821, $V_{q}$ tends to 0, $E_{u}$ tends to 0, $E_{q}$ tends to 0, $I_{u}$ tends to 0, $I_{h}$ tends to 0, $I_{a}$ tends to 0, and $R$ tends to 0. It shows that system (2.1) has a DFE and it is globally asymptotically stable with different initial values when $\mathscr{R}_{0} < 1$. The numerical simulation results support the case stated in Theorem 2.

5.   Empirical analysis

5.1.   Data sources and methods

In order to show that the model is practical, the COVID-19 pandemic mainly caused by Delta variant in South Korea is analysed and verified in this section. The experimental data come from the Central Disaster Management Headquarters of Korea ^[25]. About 97% of confirmed cases were infected by the Delta variant on August 22, 2021 ^[26], so August 22, 2021 is chosen as the starting date of the data set. The data in this paper is the official data of COVID-19 in Korea from August 22, 2021 to October 20, 2021.

In order to verify the training effect, the data set is divided into training set (from August 22 to October 14, 2021, 54 days in total) and verification set (from October 15 to 20, 2021, 6 days in total) according to the ratio of 9:1. The unknown parameters of the model are obtained by fitting the real data of the training set, and the prediction effect of the model is tested by using the real data of the verification set. According to the Korea Disease Control and Prevention Agency (KDCA) ^[27], since September 24, 2021, the vaccinated individuals can be exempted from quarantine if they show no symptoms after close contact with the confirmed individuals. Therefore, this paper considers the training set in two stages. The first stage is from August 22, 2021 to September 23, 2021, and the second stage is from September 24, 2021 to October 14, 2021.

5.2.   Parameter fitting

Before fitting the parameters, we fix some parameters from previous literatures and studies to reduce the complexity.

Since the recommended duration in quarantine for people suspected of being exposed to COVID-19 in South Korea was 14 days ^[28], we set $\omega$ = 1/14 per day. Further, Zhang et al. estimated that the mean incubation period of Delta variant was 4.4 days ^[29], which was shorter than that of the original strain reported by Yin et al. (4.4 vs. 6.8) ^[30]. Thus, we consider $\sigma_{u}$ = $\sigma_{q}$ = 1/4.4 per day. Some studies have suggested that most COVID-19 infections (over 95.6%) show moderate, mild, or asymptomatic infections, about 3.08% show severe symptoms (but without requiring ICU admission), and 1.32% show critically-severe symptoms requiring ICU admission ^[31,32]. Consequently, we set $f_{2}$ = 0.044. According to ^[33], the proportion of asymptomatically-infectious individuals is 16%. Hence, we get $f_{1}$ = 0.80304 and $1-f_{1}-f_{2}$ = 0.15296 respectively. The modification parameter for the relative infectiousness of asymptomatically-infectious individuals ($\eta_{a}$) was estimated from ^[31,34] to be 0.5, so we set $\eta_{a} = 0.5$. Due to the strict quarantine policy in South Korea ^[28], we assume that the ineffectiveness of quarantine $\theta = 0$. The average lifespan of Korean is 83.3 years ^[35]. Therefore, we have the natural death rate $d = 1/(83.3\times365)$ per day. The latest total population of South Korea is about $N$(0) = 51,821,669 ^[36], hence, the recruitment rate $\Lambda$ is obtained from $\Lambda/d$ = $N(0)$, so $\Lambda = 1704.41$ per day. The World Health Organization (WHO) explicitly wrote in the document ^[37] that the critical characteristic or the minimal requirement is to confer protection for at least six months. Therefore, we assume that the average duration of vaccine immunity is 180 days, then $\eta_{v} = 1/180$ per day. It can be seen from ^[38] that COVID-19 vaccine was on average 86.6% effective in preventing infections, so we select $\varepsilon_{v}$ = 0.134. According to ^[39], the maximum and minimum values of disease-induced mortality rate from August 22 to October 14 were 0.93% and 0.82%, respectively, and we take the average value as 0.875%. Thus we set $\delta_{u} = \delta_{h} = 0.00875$ per day. To obtain estimation for $\delta_{a}$, we assume that $\delta_{a} = 0.5\delta_{u}$ (so that $\delta_{a} = 0.004375$ per day). According to ^[39] and calculation, vaccination rate $\xi_{v}$ is about $7.6\times10^{-3}$ per day. As can be seen from ^[27], the vaccinated individuals, like the unvaccinated individuals, needed to be quarantined, so $p_{1} = p_{2} = 1$ in the first stage. However, since the vaccinated individuals did not need to be quarantined since September 24 ^[27], $p_{1} = p_{2} = 0$ in the second stage. In addition, we fix the initial hospitalized population $I_{h}(0)$ and total population $N(0)$ as 27,959 and 51,821,669 respectively according to the data information. From ^[40], it can be known the number of newly confirmed cases per day confirmed by temporary screening. Let the vaccination proportion on August 22 be $\xi$. According to ^[41], $\xi$ = 22.5%. The default parameter values and the initial values of state variables are given in Tables 3 and 4, respectively.

The remaining parameters are obtained by fitting the data in the training set. The model is solved by using lsqcurvefit function in Matlab, which utilizes a nonlinear least squares method to estimate parameters. The fitting parameters in two stages are given in Table 5.

5.3.   Model fitting and prediction

In order to evaluate the fitting and prediction effect, the mean absolute percentage error (MAPE) is used in this paper. MAPE between the cumulative number of confirmed cases obtained by model fitting and the actual data in the training set is used to evaluate the fitting effect of the model. Similarly, MAPE between the cumulative number of confirmed cases obtained by model prediction and actual data in the verification set is used to evaluate the prediction effect of the model. MAPE is expressed as:

where, $C_{t}$ is the actual cumulative number of confirmed cases on the $t$ th day.

The best fits to the reported data of the two stages via our model are depicted in Figure 3. Further, using the above parameters, we make a prediction on the verification set, and the prediction result is shown in Figure 4. As can be seen, there is a good agreement between the prediction result and the real data. The fitting and prediction effect of the model are shown in Table 6.

6.   Numerical simulations of quarantine and vaccine impact

For simplicity of presentation, we make $\theta_{1} = 1-\theta$ and $\varepsilon = 1-\varepsilon_{v}$, i.e., $\theta_{1}$ is quarantine effectiveness and $\varepsilon$ is vaccine efficacy.

6.1.   Numerical simulations of quarantine impact

First of all, we simulate the impact of the quarantine rate ($q$) on the spread of COVID-19. As shown in the Figure 5(a), if the quarantine rate decreases, the cumulative number of deaths will sharply increase, which negatively affects prevention and control of COVID-19.

Next, we consider the impact of quarantine effectiveness ($\theta_{1}$) on COVID-19 transmission. The results are shown in Figure 5(b), which imply that implementing strict quarantine measures is an effective means to prevent the development of the COVID-19 pandemic.

6.2.   Numerical simulations of vaccine impact

We firstly assess the impact of the vaccination rate ($\xi_{v}$) on the prevention of COVID-19 transmission. The simulation results are presented in Figure 6(a), which show that a significantly large increasement in vaccination rate is necessarily needed for the control of COVID-19.

As can be seen from the comparison results in Figure 6(b), when we increase vaccine efficacy ($\varepsilon$), the cumulative number of deaths markedly decreases. Meanwhile, it is easy to find that the stronger the effectiveness of vaccine, the faster the elimination of pandemic.

Finally, we simulate the impact of rate of immunity loss ($\eta_{v}$) on dynamic transmission of the COVID-19. Depicted in Figure 6(c), the longer duration of vaccine immunity can have a better suppression effect on the spread of pandemic.

7.   Discussion

A mathematical model for mutated COVID-19 combining quarantine, isolation and vaccination is established in this paper, and it is proved that the model is globally asymptotically stable at the DFE if $\mathscr{R}_{0} < 1$, and the correctness of this result is verified by numerical simulation. The unknown parameters are obtained by fitting the COVID-19 epidemic data in South Korea from August 22 to October 14, 2021 in stages, and MAPE in two stages are 0.23 and 0.22%, respectively. The fitting effect of the model is satisfying. In the prediction of the data set from October 15 to October 20, 2021, MAPE is 0.047%, and the prediction effect is satisfactory, which also shows that the model is feasible. Numerical simulations indicate that it is essential to implement strict quarantine measures and strengthen contact tracing for combating the spread of COVID-19. Meanwhile, increasing vaccination rate and getting the vaccines with more effective and the longer duration of vaccine immunity have positive impact on the prevention of pandemic transmission.

However, there are still many issues worthy of further study. For instance, most parameters (such as quarantine rate and vaccination rate) are dynamically changing during the development of COVID-19, but these characteristics are not taken into account in the current model. Thus, dynamic model with stochastic parameters deserves further study. Besides, the severity of infection is different in reality, so we can consider further the confirmed individuals with common and severe symptoms. These problems are of great interest which can be left for studying in the near further.

Acknowledgements

The authors are grateful to the editors and the reviewers for their valuable comments and suggestions. This work was supported by the Fundamental Research Funds for the Central Universities (NO.2572018BC20).

Conflict of interest

The authors declare there is no conflict of interest.