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

Abstract: Multiple variants of SARS-CoV-2 have emerged but the effectiveness of existing COVID19 vaccines against variants has been reduced, which bring new challenges to the control and mitigation of the COVID-19 pandemic. In this paper, a mathematical model for mutated COVID-19 with quarantine, isolation and vaccination is developed for studying current pandemic transmission. The basic reproduction number R0 is obtained. It is proved that the disease free equilibrium is globally asymptotically stable if R0 < 1 and unstable if R0 > 1. And numerical simulations are carried out to illustrate our main results. The COVID-19 pandemic mainly caused by Delta variant in South Korea is analyzed by using this model and the unknown parameters are estimated by fitting to real data. The epidemic situation is predicted, and the prediction result is basically consistent with the actual data. Finally, we investigate several critical model parameters to access the impact of quarantine and vaccination on the control of COVID-19, including quarantine rate, quarantine effectiveness, vaccination rate, vaccine efficacy and rate of immunity loss.


Introduction
Corona Virus Disease 2019 , 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 facemasks 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 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.

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: (i) The recruitment rate (either by birth or by immigration) of the population is given by a constant rate Λ and they are all susceptible.
(ii) 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.
(iii) After recovery the individuals have immunity to COVID-19 in a short time.
(iv) A part of vaccinated individuals will lose their immunity and rejoin the susceptible compartment.
(v) After the immunity period, the vaccinated individuals become susceptible again.
(vi) Hospitalization is completely effective.
(vii) 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 η u and η a respectively(0 < η u , η a < 1).
(viii) The contact rate c is the same for the individuals in symptomatically-infectious I u , nonquarantined exposed E u and the asymptomatically-infectious I a compartments.
(ix) The latency period is defined as the days from exposure to the onset of illness. (x) 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.
(xi) 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.
(xii) The vaccinated individuals will not lose immunity during the quarantine period.
(xiii) 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.
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)). Population of non-quarantined exposed individuals E q Population of quarantined exposed individuals I u Population of non-hospitalized symptomatically-infectious individuals I h Population of hospitalized symptomatically-infectious individuals I a Population of asymptomatically-infectious individuals with mild or no clinical symptoms of the disease R Population of recovered individuals Then, for any t > 0, we have the non-negative solution for the system (2.1).
The following result is given using the first equation of the model (2.1), Then, it can be written as, Thus, For the remaining equations, we take the same steps to show H i (t) > 0 (i = 2, 3, ..., 10) for every t > 0. □ Latency period for non-quarantined (quarantined) exposed individuals 5.1  Let the closed and biologically feasible region be Ω, shown by Ω = Next, we present the following results for this feasible region.
The region Ω is a positively invariant and attracting region for system (2.1) with the nonnegative 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 Thus, the region Ω 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 Ω.
The basic reproduction number 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, By the definition of the basic reproduction number R 0 , we get where, The quantity 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.

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.
, we can obtain the Jacobian matrix of the system corresponding toH 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 , and the other eigenvalues of J 4 are determined by the equation > 0, then a 3 > 0. Meanwhile, we can also obtain that The first inequality is established due to 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 R 0 < 1. On the other hand, we have Proof. Let 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 R 0 . Let u = (η u , θη u , 1, 0, η a ). It then follows from the fact R 0 = ρ(FV −1 ) = ρ(V −1 F) and direct calculation that u is a left eigenvector of the matrix V −1 F, i.e., uV −1 F = R 0 u.

Consider a Liapunov function
Differentiating the above equation, we have If R 0 < 1, the equality dL 0 dt = 0 implies that uX = 0. This leads to E u = E q = I u = I a = 0 by noting that all components of u except for the fourth element are positive. Clearly, let f 2 σ u E u (τ) + rσ q E q (τ) + ϕ u I u (τ) + σ a I a (τ) e K 4 τ dτ e −K 4 t .
Using L'Hospital's rule, In a similar way, for the 2nd, 4th and 10th equation of the model (2.1), we have lim t→∞ S q (t) = lim t→∞ V q (t) = lim t→∞ R (t) = 0. Hence, we can obtain the limit system of (2.1) as follows.
It is easy to solve that Clearly, we can get lim t→∞ S u (t) = S 0 , lim t→∞ V u (t) = V 0 , so (S 0 , V 0 ) is globally attractive for the limit system (3.2). By discussion of linearization systems for (3.2), it can be known that (S 0 , V 0 ) is locally asymptotically stable if R 0 < 1. Combining global attraction of (S 0 , V 0 ), we can obtain that (S 0 , V 0 ) 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 R 0 < 1. In this section, we present the numerical solution of the system (2.1). We can get 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 R 0 < 1. The numerical simulation results support the case stated in Theorem 2.

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  Figure 2. The dynamical behavior for different initial conditions of non-quarantined susceptible individuals S u , quarantined susceptible individuals S q , non-quarantined vaccinated individuals V u , quarantined vaccinated individuals V q , non-quarantined exposed individuals E u , quarantined exposed individuals E q , non-hospitalized symptomatically-infectious individuals I u , hospitalized infectious individuals I h , asymptomatically-infectious individuals I a and recovered individuals R, subfigures (a)-(j) represent them, respectively.
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.

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 ω = 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 σ u = σ 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 (η a ) was estimated from [31,34] to be 0.5, so we set η a = 0.5. Due to the strict quarantine policy in South Korea [28], we assume that the ineffectiveness of quarantine θ = 0. The average lifespan of Korean is 83.3 years [35]. Therefore, we have the natural death rate d = 1/(83.3 × 365) per day. The latest total population of South Korea is about N(0)=51,821,669 [36], hence, the recruitment rate Λ is obtained from Λ/d = N(0), so Λ = 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 η 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 ε v =0.134. According to [39], the maximum and minimum values of diseaseinduced 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 δ u = δ h = 0.00875 per day. To obtain estimation for δ a , we assume that δ a = 0.5δ u (so that δ a = 0.004375 per day). According to [39] and calculation, vaccination rate ξ v is about 7.6 × 10 −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 ξ. According to [41], ξ=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.

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.
Derived from the number of newly confirmed cases per day confirmed by temporary screening, ω and ξ V u 11,457,611  The cumulative number of confirmed cases

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 (θ 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.

Numerical simulations of vaccine impact
We firstly assess the impact of the vaccination rate (ξ 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 (ε), 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 (η 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.

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 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.