The modeling and analysis of the COVID-19 pandemic with vaccination and isolation: a case study of Italy

1.   Introduction

COVID-19, caused by the SARS-CoV-2 coronavirus, has become a widespread epidemic, the most dangerous infectious disease to people's health in the recent decade. The first documented instance of the sickness was in Wuhan, China, in December 2019, and subsequent cases have been recorded internationally ^[1]. On January 30, 2020, the World Health Organization (WHO) announced an international public health emergency, and on March 11, 2020, it declared a pandemic. According to the WHO, more than 536 million confirmed cases and more than 6.3 million deaths had been documented worldwide as of June 19, 2022.

In comparison to other viral infections such as influenza, COVID-19 has the following characteristics ^[2]: it has a long incubation period; the spread of this disease is difficult to control because it can be transmitted by asymptomatic patients ^[3], and it can be detected by polymerase chain reaction (PCR) tests. This virus can be fatal, especially in the elderly or those suffering from underlying conditions such as cancer, diabetes, or arterial hypertension. The primary known transmission mechanism is through direct social interaction between a vulnerable individual and an infected person ^[4]. Susceptible individuals can also develop the virus by direct contact with contaminated surfaces or objects, inhaling droplets ejected from the noses of symptomatic or asymptomatic infected individuals while spitting, coughing, or sneezing ^[5]. An infected person develops symptoms 5–6 days after the first infection, but the incubation period is usually 2–14 days ^[6,7]. A large proportion of patients with COVID-19 disease have no clinical symptoms or negligible disease symptoms. However, common symptoms in infected persons include a dry cough, high fever, discomfort, headache, severe weariness, and shortness of breath ^[8]. Patients with the most severe symptoms may develop pneumonia or possibly organ failure.

This pandemic has tested the public health systems of all countries. In order to control the spread of the disease, various countries have taken measures, such as quarantines, lockdowns, and social evacuation measures. COVID-19 has wreaked havoc in many countries. In this paper, we will apply our model to the case of COVID-19 transmission in Italy. Italy is one of the European countries most severely affected by COVID-19. The first cases in Italy were reported at the end of January 2020. In early March 2020, Italy was forced to carry out extensive containment and control to stabilize the epidemic. On March 19, 2020, 41,035 detected cases were reported in Italy, more than China and Iran combined, making it the country with the worst epidemic worldwide. As of March 19, 2020, the number of deaths in Italy reached 3405, making it the country with the highest number of deaths due to COVID-19 in the world. The Italian government has taken strong measures to avoid the continued spread of the virus among the population. These actions include social evacuation, mandatory use of masks, and temporary cancellation of cultural, sports, and educational activities.

Large-scale vaccination campaigns have been initiated in some countries, such as Italy, since December 2020. However, the global distribution of the COVID-19 vaccine is a complicated and lengthy process. Therefore, the effects of this fledgling vaccination campaign are not yet reflected in the number of infections reported daily.

Mathematical models have an important role in the study of the spread of diseases. They can help us to simulate, analyze and predict the spread of epidemics ^{[9,10,11,12,13,14]}. Many researchers have now used mathematical models to study the transmission of COVID-19, such as the models in ^{[15,16,17,18,19,20,21,22,23]}. ^[24] developed a stochastic discrete-time infectious disease model to evaluate the risk of COVID-19 resurgence in the second wave. ^[25] compared the situation and control measures during the outbreak phase in Italy with Guangdong Province, China, and they concluded that the control measures in Italy were not timely and effective at the beginning. ^[26] developed a SIFR model to investigate the spread of COVID-19 in Harbin, China. ^[27] analyzed the effect of mask-wearing on the spread of COVID-19 in New York City. ^[28] studied the optimal control of COVID-19 in Saudi Arabia. ^[29] conducted an in-depth study of nosocomial infections by developing a SEIHR mathematical model of COVID-19 transmission dynamics. ^[30] studied the spread of COVID-19 in Spain by a new extended SEIR model. There are also several studies on the effect of vaccination on the spread of COVID-19. ^[31] studied a SE (Is) (Ih) AR model combining vaccine and antiviral control. ^[32] studied the long-term dynamics of the IVRD model using efficient numerical methods. ^[33] proposed a mathematical model to study the effect of two doses of vaccine and non-pharmaceutical prophylaxis on the spread of COVID-19 in Bangladesh. In ^[34], the impact of vaccination on the control of COVID-19 in the United States in 2021 is studied using a mathematical model. ^[35] evaluated the effect of age-specific vaccination strategies on the number of infections and deaths using a SEIR model. ^[36] developed a patch model with the vaccine to study the spread of COVID-19 in China based on the mechanism of the spread of the COVID-19 epidemic. Several other studies have discussed vaccination prioritization strategies, for example ^[37,38].

The motivation for the study in this paper is ongoing vaccination to control the spread of novel coronaviruses in humans. Given the same rate of vaccine protection, how effective is vaccination in controlling symptomatic and asymptomatic infected individuals, respectively? What is the relationship between the rate of vaccine protection and the number of symptomatic infected persons, and the number of dead cases? Do we need to take other measures besides vaccination to control the spread of the epidemic? Here, we developed a SEIAR-based compartment model to simulate the spread of the virus after vaccination.

In the study of this paper, we propose a SEIAHRV model containing a vaccination compartment. The developed model is analyzed and numerically simulated to describe the propagation dynamics of COVID-19. The next-generation matrix approach is applied to derive the control reproduction number. We determine the uniqueness of the existence of disease-free equilibrium and endemic equilibrium and also analyze the stability of equilibria. In numerical simulations, we analyze the control effect of the vaccine on the number of asymptomatic infected, symptomatic infected, and fatal cases. We perform sensitivity analyses on different parameters to identify the model parameters that significantly impact the control reproduction number of COVID-19. The results of this paper provide recommendations for policymakers on the effective mitigation of the COVID-19 epidemic.

The results of the analysis of the role of vaccines show that vaccination is more effective in controlling the increase in the number of symptomatic infected individuals in the population; in reducing the likelihood of people becoming symptomatic and dying after infection, we should choose vaccines with a higher rate of vaccine protection. The results of the sensitivity analysis suggest that if vaccination alone is not enough to stifle the spread of the epidemic in Italy in susceptible individuals, manual interventions (e.g., wearing masks, maintaining social distance, strengthening quarantine of the population, etc.) must be continued after vaccination for the spread of the virus to be significantly limited. With no overload of medical resources, increasing the quarantine isolation rate for the population to the original 150% would reduce the number of detected cases by 81.25% after 150 days; keeping the current quarantine isolation rate unchanged would reduce the number of detected cases by 53.23% after 150 days; reducing the quarantine isolation rate to the original 50% would only reduce the number of detected cases by 4.83% after 150 days. After the relaxation of the quarantine of the population, the number of detected cases seemed to "decrease" in the short term because the infected people in the population were not detected in time, so the disease would spread rapidly after a while and eventually, the spread of the epidemic could not be controlled.

The rest of the paper is organized as follows. In Section 2, we present the established disease transmission model. Section 3 presents the equilibria of the model and the control reproduction number. Section 4 analyzes the stability of the equilibria. Section 5 fits the model parameters to simulate the epidemic's spread and analyzes the parameters' sensitivity. Finally, a brief discussion and model conclusions are given in Section 6.

2.   Spread model

According to the transmission characteristics of COVID-19 and preventive measures, we divided the population into eight groups, as shown in Table 1. Individuals in group $A$ or group $I$ are infected and infectious. The transformation relationship between compartments is shown in model (1), the state transfer schematic is shown in Figure 1, and the parameters are shown in Table 2. The range of parameters in Table 2 is $\left[\mathrm{0, 1}\right]$. $M$ is the number of new populations in the area per day, including the number of births and population in-migration.

where $N = S+E+I+A+H+R+V$.

For model (1), we have the following assumptions

(a) Individuals in compartment $H$ are entirely isolated and therefore have no opportunity to infect susceptible individuals. Asymptomatic individuals will be isolated at home immediately after being detected. Symptomatic individuals will be isolated in the hospital immediately after being detected.

(b) Individuals in compartment $E$, although infected with the virus, are still in the latent phase and are not infectious. Individuals in compartments $I$ and $A$ are both infectious.

(c) Individuals in compartment $V$ are vaccinated, but the vaccine may not protect everyone. Individuals for whom the vaccine is effective will be protected from infection, while unprotected individuals will return from compartment $V$ to compartment $S$, i.e., become susceptible again.

Let

be the solution of model (1) through any ${\rm{ \mathsf{ φ} }}: = {\left({\rm{ \mathsf{ φ} }}_{1}, {\rm{ \mathsf{ φ} }}_{2}, {\rm{ \mathsf{ φ} }}_{3}, {\rm{ \mathsf{ φ} }}_{4}, {\rm{ \mathsf{ φ} }}_{5}, {\rm{ \mathsf{ φ} }}_{6}, {\rm{ \mathsf{ φ} }}_{7}\right)}^{T}\in {\Omega }$, where

with ${\mathbb{R}}_{+} = [0, \infty).$

For the well-posedness and dissipativeness of model (1), we have the following proposition

Proposition 1. Model (1) with any ${\rm{ \mathsf{ φ} }}\in {\Omega }$ has the unique non-negative bounded solution $U\left(t\right)$ on $\left[0, \mathrm{\infty }\right)$.

Proof. From existence and uniqueness theorem of solutions of ordinary differential equations (ODEs) ^[39], the solution $U\left(t\right) = {\left({U}_{1}\left(t\right), {U}_{2}\left(t\right), {U}_{3}\left(t\right), {U}_{4}\left(t\right), {U}_{5}\left(t\right), {U}_{6}\left(t\right), {U}_{7}\left(t\right)\right)}^{T}$ with any ${\rm{ \mathsf{ φ} }}\in {\Omega }$ is unique on its maximum interval $[0, {T}_{\varphi })$ of existence.

First, let us show that the solution $U\left(t\right)$ is non-negative on $[0, {T}_{\varphi })$. For any given $b\in (0, {T}_{\varphi })$, and any given $\varepsilon > 0$, $U\left(t, \varepsilon \right) = {\left({U}_{1}\left(t, \varepsilon \right), {U}_{2}\left(t, \varepsilon \right), {U}_{3}\left(t, \varepsilon \right), {U}_{4}\left(t, \varepsilon \right), {U}_{5}\left(t, \varepsilon \right), {U}_{6}\left(t, \varepsilon \right), {U}_{7}\left(t, \varepsilon \right)\right)}^{T}$ is the solution of the following model

with ${\rm{ \mathsf{ φ} }}\in {\Omega }$ for $t\in \left[0, b\right]$. According to the dependence of solutions of ODEs on parameters ^[39], the solution $U\left(t, \varepsilon \right)$ is uniformly existent on $[0, b]$ for sufficiently small $\varepsilon > 0$. Since $\frac{d{U}_{i}\left(0, \varepsilon \right)}{dt} > 0$ if ${U}_{i}\left(0, \varepsilon \right) = 0$ for any $i\in {I}_{7} = \{\mathrm{1, 2}, \mathrm{3, 4}, \mathrm{5, 6}, 7\}$, we can claim that $U\left(t, \varepsilon \right)\gg {\bf{0}}$ on $(0, b)$. We prove the claim by contradiction. In fact, suppose that there exists a $\stackrel{-}{t}\in \left(0, b\right)$ such that ${U}_{i}\left(\stackrel{-}{t}, \varepsilon \right) = 0$ for some $i\in {I}_{7}$ and $U(t, \varepsilon)\gg {\bf{0}}$ for $t\in \left(0, \stackrel{-}{t}\right)$, where $\stackrel{-}{t} = \underset{1\le i\le 7}{\mathrm{min}}\left\{{t}_{i}\right\}, {t}_{i} = \mathrm{s}\mathrm{u}\mathrm{p}\left\{t\in \left(0, b\right)|{U}_{i}\left(x, \varepsilon \right) > 0, x\in (0, t]\right\}$. Then, we have

We denote $P = N-H-V = S+E+I+A+R$. For $t\in \left[0, \stackrel{-}{t}\right]$, we can obtain

where $\xi = \mathrm{m}\mathrm{a}\mathrm{x}\left\{\eta, \delta, \theta, \mu \right\}$. Hence, we can obtain $P\left(t, \varepsilon \right) > P\left(0, \varepsilon \right){e}^{-2\xi t} > 0$ for $t\in \left[0, \stackrel{-}{t}\right]$. Since ${U}_{i}\left(\stackrel{-}{t}, \varepsilon \right) = 0$ and $U\left(\stackrel{-}{t}, \varepsilon \right)\ge {\bf{0}}$, it follows from model (2) that $\frac{\mathrm{d}{U}_{i}(\stackrel{-}{t}, \varepsilon)}{\mathrm{d}t} > 0$. Clearly, this contradicts Eq (3). Therefore, we have $U\left(t, \varepsilon \right)\gg {\bf{0}}, t\in (0, b)$.

Now, let $\varepsilon \to {0}^{+}$, then it follows that $U\left(t, 0\right) = U\left(t\right)\ge {\bf{0}}$ for any $t\in [0, b)$. Consider any arbitrary $b\in (0, {T}_{\varphi })$, then it follows that $U\left(t\right)\ge {\bf{0}}$ for any $t\in [0, {T}_{\varphi })$. By model (1), we obtain $\frac{dN\left(t\right)}{dt} = M-\kappa H\left(t\right)-\mu N\left(t\right)\le M-\mu N\left(t\right).$ From comparison principle, $U\left(t\right)$ is bounded on ${\mathbb{R}}_{+}^{7}$. According to the extension theorem of the solution of ODEs ^[39], ${T}_{\varphi } = \infty$.

By the same argument as above, we can get Proposition 1.

Remark 1. From the proof of Proposition 1, it follows that the set ${\Omega }$ is positively invariant for model (1).

3.   Equilibria and control reproduction number

The equilibria and the control reproduction number of model (1) are discussed in this section.

It is clear that model (1) always has a disease-free equilibrium ${P}_{0} = {({S}_{0}, \mathrm{0, 0}, \mathrm{0, 0}, 0, {V}_{0})}^{T}$, where ${S}_{0} = \frac{M(1-\alpha +\mu)}{\mu (1-\alpha +\mu +\eta)}, {V}_{0} = \frac{M\eta }{\mu (1-\alpha +\mu +\eta)}$.

Next, we use the next-generation matrix method ^[40] to calculate the control reproduction number. Let

The Jacobi matrices of $F$ and $V$ at the disease-free equilibrium ${P}_{0}$ are as follows

Thus, we can get the control reproduction number ${R}_{c}$ as follows

Let the right-hand sides of the equations in model (1) equal to 0, we can get

Adding the first two equations in model (1) and then put Eq (4) into it, we can obtain

From Eqs (4) and (5), model (1) admits a unique endemic equilibrium ${P}_{e} = {({S}_{e}, {E}_{e}, {I}_{e}, {A}_{e}, {H}_{e}, {R}_{e}, {V}_{e})}^{T}\gg {\bf{0}}$ when and only when $0 < {E}_{e} < \frac{M}{\gamma +\mu }$. Since $P = S+E+I+A+R$, putting it into the first equation of model (1), we get

where $q = 1+\frac{\phi \gamma }{\delta +\mu }+\frac{(1-\phi)\gamma }{\theta +\mu }+\frac{\lambda }{\mu \left(\lambda +\kappa +\mu \right)}\left(\delta \frac{\phi \gamma }{\delta +\mu }+\theta \frac{\left(1-\phi \right)\gamma }{\theta +\mu }\right)$. Putting Eq (5) into Eq (6) and simplifying it, we have

According to Eqs (5) and (7), there holds

We can see that Eq (8) has a positive root

when and only when ${R}_{c} > 1$. Therefore, model (1) has a unique endemic equilibrium ${P}_{e}$ if and only if ${R}_{c} > 1$, where

By Eqs (4), (5) and (8), we found that the disease-free equilibrium ${P}_{0}$ is unique.

Based on the above analysis, we have the following theorem.

Theorem 1. Model (1) has a unique disease-free equilibrium ${P}_{0}$. And model (1) has a unique endemic equilibrium ${P}_{e}$ if and only if ${R}_{c} > 1$.

4.   Stability analysis of equilibria

In this section, the stability of the disease-free equilibrium and the endemic equilibrium is analyzed, with respect to ${R}_{c}$.

Theorem 2. The disease-free equilibrium ${P}_{0}$ is locally asymptotically stable if ${R}_{c} < 1$ and unstable if ${R}_{c} > 1$.

Proof. The characteristic equation of the linearized system of model (1) at ${P}_{0}$ is

It can be seen that Eq (9) has at least four roots with negative real parts ${X}_{1} = {X}_{2} = -\mu$, ${X}_{3} = -(\kappa +\lambda +\mu)$, ${X}_{4} = -(1-\alpha +\mu)$. For simplicity, let ${k}_{1} = \gamma +\mu, {k}_{2} = \delta +\mu, {k}_{3} = \theta +\mu$. Therefore, we only need to discuss the roots of the following equation

Equation (10) can be organized to obtain

Simplifying Eq (11), we have ${X}^{3}+{a}_{1}{X}^{2}+{a}_{2}X+{a}_{3} = 0$, where ${a}_{1} = {k}_{1}+{k}_{2}+{k}_{3}$, ${a}_{2} = {k}_{2}{k}_{3}+{k}_{1}{k}_{2}(1-{R}_{1})+{k}_{1}{k}_{3}(1-{R}_{2})$, ${a}_{3} = {k}_{1}{k}_{2}{k}_{3}(1-{R}_{c})$, ${R}_{1} = \frac{\beta \gamma \phi }{{k}_{1}{k}_{2}}$, ${R}_{2} = \frac{\beta \gamma \left(1-\phi \right)\tau }{{k}_{1}{k}_{3}}$. Obviously if ${R}_{c} < 1$, ${a}_{1} > 0$ and ${a}_{3} > 0$, and thus easy to determine ${a}_{1}{a}_{2} > {a}_{3}$. Therefore, according to the Routh-Hurwitz criterion it is known that the disease-free equilibrium is locally asymptotically stable if ${R}_{c} < 1$.

Then consider that if ${R}_{c} > 1$, Eq (11) can be written as

Let ${X}_{5}$, ${X}_{6}$ and ${X}_{7}$ be the three roots of Eq (11), and according to Veda's theorem we have

There are two cases about ${X}_{5}$, ${X}_{6}$ and ${X}_{7}$ implied by Eq (12)

(ⅰ) All three roots are real numbers, two of them are negative, and the other is positive.

(ⅱ) One of the three roots is real and positive, and the other two roots are conjugate complexes.

In summary, Eq (9) has a positive root if ${R}_{c} > 1$, so the disease-free equilibrium is unstable.

Theorem 3. The disease-free equilibrium ${P}_{0}$ is globally asymptotically stable if ${R}_{c} < 1$ in ${\Omega }$.

Proof. By Theorem 2, we only need to prove that the disease-free equilibrium ${P}_{0}$ is globally attractive for ${\mathrm{R}}_{\mathrm{c}} < 1$. From the above we can get

is the solution of model (1) with any ${\rm{ \mathsf{ φ} }}\in {\Omega }$. To show that ${P}_{0}$ is globally attractive, we just need to verify that $\omega \left(\varphi \right) = \left\{{P}_{0}\right\}$, where $\omega \left(\varphi \right)$ is the ${\rm{ \mathsf{ ω} }}$ limit set of $\varphi$. Let we define a function $\mathcal{L}$ on ${\Omega }$ as follows,

The derivative of $\mathcal{L}$ along $U\left(t\right)$ is given by

With some computations, we get

We can find that when ${R}_{c} < 1$, $\dot{\mathcal{L}}\left(U\left(t\right)\right)\le 0$. Hence if ${R}_{c} < 1$, $\mathcal{L}$ is a Lyapunov function on ${\Omega }$. By Corollary 2.1 in ^[41], we have $\dot{\mathcal{L}}\left(\phi \right) = 0$ for any $\phi \in \omega \left(\varphi \right)$.

For any $\psi \in {\rm{ \mathsf{ ω} }}\left({\rm{ \mathsf{ φ} }}\right)$, assume that $U\left(t\right) = {\left(S\left(t\right), E\left(t\right), I\left(t\right), A\left(t\right), H\left(t\right), R\left(t\right), V\left(t\right)\right)}^{T}$ is the solution of model (1) through $\psi$. Then by the invariance of ${\rm{ \mathsf{ ω} }}\left({\rm{ \mathsf{ φ} }}\right)$, we have that $U\left(t\right)\in {\rm{ \mathsf{ ω} }}\left({\rm{ \mathsf{ φ} }}\right)$ for all $t\in \mathbb{R}$. Next, by $\dot{\mathcal{L}}\left(U\left(t\right)\right) = 0$ and Eq (13), we have $E\left(t\right) = 0$ for $\forall t\in \mathbb{R}$. From the 3rd, 4th, 5th, and the 6th equations of model (1) and the invariance of ${\rm{ \mathsf{ ω} }}\left({\rm{ \mathsf{ φ} }}\right)$, it holds $I\left(t\right) = A\left(t\right) = H\left(t\right) = R\left(t\right) = 0$ for $t\in \mathbb{R}$. Therefore, model (1) can be simplified as

According to system (14), we can see that

From the invariance of ${\rm{ \mathsf{ ω} }}\left({\rm{ \mathsf{ φ} }}\right)$, it follows that $S+V = \frac{M}{\mu }$ for $t\in \mathbb{R}$. Hence, system (14) and the invariance of ${\rm{ \mathsf{ ω} }}\left({\rm{ \mathsf{ φ} }}\right)$ imply that $S = \frac{M(1-\alpha +\mu)}{\mu (1-\alpha +\mu +\eta)} = {S}_{0}$ and $V = \frac{M\eta }{\mu \left(1-\alpha +\mu +\eta \right)} = {V}_{0}$ for $t\in \mathbb{R}.$ Therefore, $\omega \left(\phi \right) = \left\{{P}_{0}\right\}$. Thus, the disease-free equilibrium ${P}_{0}$ is globally asymptotically stable.

By the proof approach of Theorem 2.1 in ^[42], we can obtain the local stability of ${P}_{e}$.

Theorem 4. The endemic equilibrium ${P}_{e}$ is locally asymptotically stable if ${R}_{c} > 1$.

Proof. The characteristic equation of the linearized system at ${P}_{e}$ is given by

where

Obviously, Eq (15) has negative roots $-\mu$ and $-\left(\kappa +\lambda +\mu \right)$. The other eigenvalues satisfy ${J}_{{P}_{e}} = 0$. Next, we show that any root $X$ of ${J}_{{P}_{e}} = 0$ has negative real part by contradiction. Assume that $X$ has a nonnegative real part. Then $\frac{X+b}{X+a}$ and $\frac{1}{X+a}$ also have a non-negative real part for $a, b > 0$. It follows from ${J}_{{P}_{e}} = 0$ that

By the properties of complex modulus and ${R}_{c} > 1$, it holds that

and

Clearly, this is a contradiction, and hence any root of Eq (15) has a negative real part. Therefore, ${P}_{e}$ is locally asymptotically stable for ${R}_{c} > 1$.

Theorems 2 and 3 show that if ${R}_{c} < 1$, the number of people infected with COVID-19, $I\left(t\right)$ and $A\left(t\right)$, will eventually be zero, i.e., the outbreak will be controlled. Theorem 4 shows that when ${R}_{c} > 1$, $I\left(t\right)$ and $A\left(t\right)$ will be stable around positive constants, respectively. This implies that COVID-19 will form an endemic disease. One certainly wants COVID-19 to disappear as soon as possible, which means that the control reproduction number ${R}_{c}$ should be controlled below 1.

According to the expression of ${R}_{c}$, the following measures can help control the spread of COVID-19:

(Ⅰ) Strengthen the quarantine of the population so that the number of infected individuals ($A$ and $I$) mixed in the population will be reduced (i.e., increase $\delta$ and $\theta$).

(Ⅱ) Implement control measures such as maintaining social distancing and wearing masks to reduce the rate of disease transmission (i.e., reduce $\beta$).

5.   Numerical studies and analysis

In this section, parameters are fitted and simulated to the model using real case data from Italy, and the role of vaccination is analyzed based on the simulation results. Then a sensitivity analysis of the parameters of the control reproduction number is performed.

5.1.   Parameter fitting

In this section, model (1) is used to simulate the spread of COVID-19 in Italy. Vaccination in Italy officially starts on December 27, 2020. The brand of vaccination is mainly the Pfizer vaccine ^[43]. The model simulation will use the protection rate of the Pfizer vaccine ^[44] as a parameter of the model. More than 3000 fully vaccinated (two doses) people have been reported since January 17, 2021. Therefore, we selected data from Italy from January 20, 2021, to June 20, 2021, for the simulation. Daily data on detected cases of COVID-19 in Italy from 2020 to 2021 are shown in Figure 2(b). Italy released an official ban ^[45] on non-essential population movements between regions, which came into effect on December 21, 2020. The point in time when the ban became effective, the point in time when the ban was lifted ^[46], and the point in time when the emergency started ^[47] have been marked in Figure 2(b). The case data used in the model include the number of cases that detected positive for the virus (they are hospitalized or in-home isolation), the number of deaths due to disease, and the number of recovered cases from the cure, all from actual data published in Italy ^[48]. We use $D\left(t\right)$ to denote the number of deaths due to disease, and the number of cases of death due to disease satisfies the equation $\frac{dD\left(t\right)}{dt} = \kappa H\left(t\right)$. With the increase in the number of medical personnel and supplies, as well as the accumulation of experience in treatment, the rate of death from disease gradually decreases, and the rate of cure increases. Therefore, it is assumed that both the mortality rate $\kappa = \kappa \left(t\right)$ and the cure rate $\lambda = \lambda \left(t\right)$ are functions of time, which is more in line with the actual situation.

Most cases that are detected as positive for the virus are asymptomatic patients, and home isolation will begin once they are detected as positive for the virus. Therefore, we assume that the number of positive virus detections is the number of isolated cases. We used the least squares method to fit the parameters of model (1). According to ^[49], the average life expectancy in Italy is 82, so we assume that $\mu$ is $\frac{1}{82\times 365}\approx 3.34\times {10}^{-5}$. The total population of Italy is estimated to be 59,325,000. The population mobility decreases during the epidemic, while the birth rate in Italy ^[50] is about 7 newborns per 1000 people, so $M$ is assumed to be 1170. The fitting results are shown in Figure 2(a), where the dots indicate the actual data, and the curves indicate the fitting results of the model. The fitted values of the parameters are shown in Table 3. The parameters of $\kappa \left(t\right)$ and $\lambda \left(t\right)$ are obtained by least squares fitting based on the official statistics of the number of death cases and the number of recovered cases. ${\rm{ \mathsf{ κ} }}\left(t\right)$ and ${\rm{ \mathsf{ λ} }}\left(t\right)$ are as follows

The initial values of the states of each compartment are shown in Table 4. The mean relative error can be used to evaluate the fitting effect, which is defined as follows

where $X\left(i\right)$ is the true value, $\widehat{X}\left(i\right)$ is the fitted value, and $n$ is the amount of data. The mean relative error of the model fit results for daily positive virus-detected cases can be calculated as ${E}_{MRE} = 4.21\%$.

The simulation results and the mean relative error show that the model can reflect the spread of the epidemic well. Based on the results in Table 3, it can be calculated that ${R}_{c}\approx 1.7417 > 1$, indicating that the epidemic will continue to spread in Italy. This conclusion is consistent with the realistic rebound of the epidemic in Italy in August 2021.

5.2.   Analysis of vaccine effects

Vaccination is an essential tool for pharmacological interventions to reduce the number of susceptible individuals and thus control the spread of diseases in the population. We introduced a vaccination compartment in the proposed model. We hypothesized that completing two doses of the vaccine implies only temporary immunity and that the individual may revert to the susceptible class due to the loss of immunity. Currently, the Italian government uses several ^[43] vaccines (Pfizer, Moderna, Johnson & AstraZeneca, etc.) against COVID-19. We used the symbols $\alpha$ and $\eta$ to model the effect of vaccination on the population, representing the vaccine protection rate for the susceptible population and the vaccination rate for the susceptible population. We considered seven different scenarios for the vaccine protection rate $\alpha$, which were unvaccinated, with protection rates ${\rm{ \mathsf{ α} }}$ of 50%, 60%, 70%, 80%, 90%, and 95%. We considered six different scenarios for the vaccination rate $\eta$, $50\%\eta$, $75\%\eta$, $100\%\eta$, $125\%\eta$, $150\%\eta$ and unvaccinated to simulate the spread of the epidemic in the population. Figures 3 and 4 show the effect of vaccination on the spread of COVID-19 transmission in Italy, respectively, indicating that the vaccine protection rate and vaccination rate have an important role in controlling the COVID-19 epidemic.

The time starting point in Figure 3(a)–3(h) is January 20, 2021, and the change in the number of cases relative to January 20 is calculated every 25 days. The values of all parameters are shown in Table 3, except for the different values of vaccine protection rate $\alpha$. The values in Figure 3(i) are the values in the first column of Figure 3(d) minus the values in the other six columns, and the values in Figure 3(j) are the values in the first column of Figure 3(f) minus the values in the other six columns.

As can be seen in Figure 3(a), (b), the higher the vaccine protection rate ${\rm{ \mathsf{ α} }}$, the earlier the "turning point" when the number of detected cases starts to decline, and the lower the peak number of detected cases. In the case of ${\rm{ \mathsf{ α} }}$ 70% in Figure 3(d), (f), for example, the number of asymptomatic infections increased by 37.90% after 25 days relative to January 20, while the number of asymptomatic infections without vaccination increased by 50.06% after 25 days relative to January 20, so the growth of asymptomatic infections decreased by $50.06\%\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }-\mathrm{ }37.90\%\mathrm{ } = \mathrm{ }12.16\%$; the number of symptomatic infected individuals increased by 40.93% after 25 days relative to January 20, while the number of symptomatic infected individuals at the time of unvaccinated increased by 53.66% after 25 days relative to January 20. Hence, the growth of symptomatic infected individuals decreased by $53.66\%-40.93\% = 12.73\%$. Similarly, in the case of the row with an interval of 25 days in Figure 3(d), (f), the increase in the number of symptomatic infected relative to the unvaccinated group decreased by 7.89%, 9.57%, 12.16%, 16.74%, 27.12%, and 39.76%, respectively, when the vaccine protection rate ${\rm{ \mathsf{ α} }}$ gradually increased; while the increase in the number of asymptomatic infected relative to the unvaccinated group decreased by 8.26%, 10.01%, 12.73%, 17.53%, 28.39%, and 41.61%, respectively, using the first column in Figure 3(d) minus the other six columns to obtain Figure 3(i), and using the first column in Figure 3(f) minus the other six columns to obtain Figure 3(j). By comparing the magnitude of the values, it can be found that the reduction in the increase in the number of symptomatic infections is greater than the reduction in the increase in the number of asymptomatic infections. Therefore, it can be seen that the vaccine is more effective in controlling the growth of symptomatic infections in the population when the vaccine protection rate $\alpha$ is the same. As can be seen in Figure 3(d), (f), the higher the vaccine protection rate $\alpha$, the lower the growth rate of asymptomatic and symptomatic cases. From Figure 3(g), (h), it can be seen that the higher the vaccine protection rate ${\rm{ \mathsf{ α} }}$, the lower the growth rate of the number of fatal cases. Vaccines make it possible to reduce the probability of a susceptible person becoming symptomatic after being infected and can better reduce the harm people suffer.

The time starting point in Figure 4 is January 20, 2021, and the change in cases relative to January 20 is calculated every 25 days. The values of all parameters are shown in Table 3, except for the different values of vaccination rate $\eta$.

From Figure 4(a), (b), it can be seen that the spread of the epidemic is controlled to some extent as soon as vaccination is started, regardless of the vaccination rate. From Figure 4(a)–(h), it can be seen that the higher the vaccination rate $\eta$, the lower the peak number of detected cases, and the earlier the "turning point" when the number of detected cases starts to decline, and the lower the number of symptomatic patients, asymptomatic patients, and new deaths.

The vaccine worked as expected to stifle the spread of the epidemic and reduce the harm people suffered. We must realize, however, that relying on vaccination alone will not completely control the spread of the outbreak. As seen in Figure 2(b), the number of detected cases decreased after the vaccination started. But after the ban on non-essential population movements between regions was lifted, the number of detected cases began to rise rapidly. The detected cases decreased until March 27, 2021, when the Italian Ministry of Health issued emergency measures. The above analysis shows that non-pharmaceutical interventions are also necessary to control the outbreak.

5.3.   Sensitivity analyses

In addition to vaccination, to reduce the spread of COVID-19, it is necessary to consider which control measures are most effective; in other words, which parameters should be analyzed to have the greatest impact on the control reproduction number ${R}_{c}$. This is the sensitivity of the parameters, which refers to the relative changes in the variables of interest caused by changes in some parameters. Assuming that the function $f$ is differentiable concerning the parameter $x$, the sensitivity index of $f$ to $x$ is defined as ^[53]:

The sensitivity index ${{\rm Y}}_{f}^{x}$ reflects the robustness of the function f to the variable $x$. Specifically, with the values of other variables (or parameters) held constant, if ${\mathrm{{\rm Y}}}_{f}^{x} > 0$, when $x$ increases (or decreases) by 1%, f increases (or decreases) by ${\mathrm{{\rm Y}}}_{f}^{x}\%$; if ${\mathrm{{\rm Y}}}_{f}^{x} < 0$, when $x$ increases (or decreases) by 1%, f decreases (or increases) by ${-{\rm Y}}_{f}^{x}\%$. The sensitivity indicators of the control reproduction number ${R}_{c}$ for each parameter are shown in Table 3. For example, the sensitivity index of ${R}_{c}$ to $\beta$ is 1 (${\mathrm{{\rm Y}}}_{{\mathrm{R}}_{\mathrm{c}}}^{{\rm{ \mathsf{ β} }}} = 1$), which means that if $\beta$ increases by 10% and the other parameters remain unchanged, ${R}_{c}$ will increase by 10%. Similarly, since ${\mathrm{{\rm Y}}}_{{R}_{c}}^{{\rm{ \mathsf{ δ} }}} = -0.0497$, $\delta$ increases by 10%, then ${R}_{c}$ decreases by 0.4%. As can be seen from Table 3, the parameters with larger absolute values of the susceptibility index are $\beta$, $\theta$, $\delta$, $\tau,$ $\mu,$ $\phi$ and $\gamma$, in that order. This indicates that effective measures to reduce ${R}_{c}$ include frequent hand washing, wearing masks, keeping social distance (reducing $\beta$); strengthening the detection of infected patients, and isolating as many detected infected individuals as possible, and ensuring that detected individuals are strictly isolated from susceptible individuals (increasing $\delta$ and $\theta$). This is the same result as the analysis in Section 4.

Some parameters can be adjusted for simulation and comparison to visualize the effect of corresponding prevention and control measures on COVID-19 transmission. According to the sensitivity analysis, $\beta$ has the greatest effect on the control reproduction number ${R}_{c}$, followed by $\theta$. we first focus on the effect of different $\beta$ on the spread of COVID-19. Keeping the other parameters unchanged, let $\beta$ be 100%, 80% and 50% of the original taken values, respectively. The simulation results are shown in Figure 5. The corresponding control reproduction numbers are 1.7417, 1.3934 and 0.8709, respectively. It can be seen that the smaller the $\beta$, the smaller the control reproduction number ${R}_{c}$ and the smaller the number of detected cases. This indicates that preventive and control measures for the population are indeed an effective way to control the spread of COVID-19.

Similarly, keeping the other parameters constant (where $\beta = 0.0507$), the $\theta$ are 50%, 100%, and 150% of the original values, respectively. The corresponding control reproduction numbers are calculated as 3.3929, 1.7417, and 1.1904, respectively. The simulation results are shown in Figures 6 and 7. The time starting point in Figure 7 is January 20, 2021, and the change in cases relative to January 20 is calculated every 25 days. The values of all parameters are shown in Table 3, except for the different values of quarantine isolation rate $\theta$ for asymptomatic infected persons.

It can be seen that the larger $\theta$ is, the smaller the control reproduction number ${R}_{c}$ is. It can be found in Figure 6 that the larger $\theta$ is in the early period, the more the number of detected cases; the larger $\theta$ is in the later period, the less the number of detected cases.

It can be found in Figure 7 that at the time interval of 25 days, with the increase of $\theta$, the changes in the number of detected cases relative to January 20 are a decrease of 18.54%, an increase of 2.22%, and an increase of 15.79%, respectively; at the time interval of 150 days, with the increase of $\theta$, the changes in the number of detected cases relative to January 20 are decreased by 4.83%, by 53.23%, and by 81.25%.

Combining the above analysis, we can get that the larger $\theta$ is, the more asymptomatic cases are detected and isolated in the early stage (the source of infection in the population becomes less), and the fewer cases are detected in the later stage; the smaller $\theta$ is, it seems that the number of detected cases decreases in the short term, but the number of detected cases will rise again in the later stage, and the control effect brought by quarantine and isolation is poor. The above analysis suggests that enhancing quarantine and isolation measures for the population is also an effective way to control the spread of COVID-19. Our advice to policymakers is that the government should select a more efficient vaccine and increase the intensity of quarantine ^[54] as much as possible while urging people to go out less and take personal protective measures, provided that medical resources are not overloaded.

6.   Conclusions

It has been more than two years since COVID-19 first appeared. Countless efforts and sacrifices have been made to control COVID-19. However, COVID-19 is still in a global pandemic phase. In particular, several variants of novel coronaviruses have been identified in several countries in the last year or so. These variants are more infectious, posing an even greater epidemic prevention and control challenge. Although universal vaccination against new coronaviruses has been implemented in some countries worldwide, vaccine capacity does not meet the current epidemic prevention needs. More vaccines and preventive measures with strict non-pharmacological therapeutic interventions are needed to prevent the further spread of the epidemic.

Based on this background, this paper develops a mathematical model to simulate the spread of COVID-19 in some countries that have been vaccinated. This paper analyzes the uniqueness of the existence of model disease-free and endemic equilibria and the stability of disease-free and endemic equilibria. The model parameters were estimated using actual data from Italy, and the disease transmission in Italy was simulated. The experimental results show that the model can simulate well the transmission dynamics of COVID-19 within some countries that have been vaccinated. Analysis of vaccine protection rates and vaccination rates revealed that vaccination better suppressed the increase in the number of symptomatic patients in the population and that the higher the vaccine protection rate, the greater the suppression effect. Increasing the vaccination rate was an effective preventive and control measure to control COVID-19. Through parametric sensitivity analysis, reducing population-to-population contact and increasing the rate of quarantine isolation of infected individuals are effective non-pharmaceutical interventions to control COVID-19. Increasing the rate of quarantine isolation of infected individuals will result in a significant increase in detected cases in the short term, but in the long term will significantly reduce the number of infections and lower the control reproduction number. If the rate of quarantine isolation of infected individuals is reduced, the control of COVID-19 will be poor in the long term.

Through the modeling and analysis in this paper, we hope to provide some useful insights for developing epidemic prevention measures to curb the spread of COVID-19. We must say that there are still many undetected cases in the real situation and our model can only estimate to some extent the number of infected people in the population. To build more accurate models to simulate and predict the transmission dynamics of COVID-19, we will fully consider the population and virus heterogeneity ^[55,56] in future studies and collect relevant real-world data.

Acknowledgments

This work was supported in part by the National Natural Science Foundation of China (Nos. 12171003 and 11901027), the China Postdoctoral Science Foundation (No. 2021M703426), the Pyramid Talent Training Project of BUCEA (JDYC20200327), the Bill & Melinda Gates Foundation (INV-005834) and the BUCEA Post Graduate Innovation Project (PG2022137).

Conflict of interest

The authors declare there is no conflict of interest.