A New COVID-19 Pandemic Model including the Compartment of Vaccinated Individuals: Global Stability of the Disease-Free Fixed Point

Abstract

Owing to the COVID-19 pandemic, which broke out in December 2019 and is still disrupting human life across the world, attention has been recently focused on the study of epidemic mathematical models able to describe the spread of the disease. The number of people who have received vaccinations is a new state variable in the COVID-19 model that this paper introduces to further the discussion of the subject. The study demonstrates that the proposed compartment model, which is described by differential equations of integer order, has two fixed points, a disease-free fixed point and an endemic fixed point. The global stability of the disease-free fixed point is guaranteed by a new theorem that is proven. This implies the disappearance of the pandemic, provided that an inequality involving the vaccination rate is satisfied. Finally, simulation results are carried out, with the aim of highlighting the usefulness of the conceived COVID-19 compartment model.

1. Introduction

The first study of epidemics can be found in an Egyptian medical papyrus of about 3500 years ago [1]. However, mathematical models for precisely analyzing the spread of epidemics were introduced only at the beginning of the twentieth century [2]. In [3], a numerical investigation of discrete oscillating epidemic models was investigated. Recently, attention has been focused on the COVID-19 pandemic, which broke out in December 2019, but is still affecting social and economic life across the world. Because of this, several epidemic models have been proposed, described by integer-order differential/difference equations [4], as well as fractional differential/difference equations [5,6,7]. For example, some works have presented dynamic compartment models (described by integer-order operators) that analyze the evolution of the disease over time by dividing the communities into some classes (i.e., susceptible, exposed, infected, etc). In particular, ref. [8] propose a model that considers eight stages of infection for the COVID-19 disease: susceptible (S), infected (I), diagnosed (D), ailing (A), recognized (R), threatened (T), healed (H) and extinct (E). The dynamic model, called SIDARTHE, introduces the distinction between diagnosed and non-diagnosed individuals, because the former are typically isolated and hence less likely to spread the infection. In [9], the dynamic behavior of the classical susceptible-infectious-removed (SIR) model when applied to the transmission of COVID-19 disease has been studied. The model includes both a nonlinear removal rate (related to the hospital bed population ratio) and the effects of media on public awareness. In [10], the inverse problem in epidemiology is exploited via a SIR model. At first, the method is used for estimating the infectivity and recovery rates from real COVID-19 data. Then, the estimated rates have been used to compute the evolution of the COVID-19 disease. In [11], a novel SIR model for simulating the COVID-19 epidemic in Wuhan is presented. The transmission dynamic model in [11] is updated with real-time input data and enriched with additional data sources, in order to infer a preliminary set of clinical parameters that could guide public health decision-making. In [12], a nonlinear SIR epidemic dynamic system is introduced to model the spread of COVID-19 under the effect of social distancing imposed by the government measures. In [13], a novel deterministic mathematical model of the COVID-19 pandemic is developed. In particular, by using time-series plots and phase portraits, references [13,14] show that the conceived COVID-19 model exhibits chaotic behaviors. Referring to systems described by fractional operators, in [15], a non-integer order COVID-19 dynamic model is proposed, which incorporates the reinfection rate in the individuals recovered from the disease. In [16], a fractional-order COVID-19 model involving the Caputo derivative is introduced. Moreover, the stability of the steady-states and the existence of non-negative solutions are investigated. In [17], a fractional compartmental model for predicting the spread of the COVID-19 pandemic is developed. In particular, reference [17] discusses the conditions under which the disease-free and the endemic equilibrium points become asymptotically stable. In [18], the spread of the COVID-19 pandemic is modeled via the Caputo–Fabrizio derivative, as well as via the Atangana–Baleanu–Caputo derivative. Then, the local stability of the equilibrium points is investigated for both modeling approaches. In [19], a fractional dynamic model for evaluating the consequences of adaptive immune responses to the COVID-19 viral mutation is presented. In [14], the dynamics of a novel COVID-19 pandemic model described by fractional derivatives are investigated. Several chaotic behaviours, obtained by varying the order of the derivative, are proven to exist in the conceived model.

Some papers, starting in the beginning of the year 2021, began to take into account, in different ways, the role of vaccination in reducing the spread of the COVID-19 pandemic. A dynamic model based on eight state variables is proposed in [20], for instance, where vaccination is viewed as a preventive action (rather than as a system variable). In [21], a new discrete susceptible-exposed-infectious-recovered (SEIR) epidemic model is illustrated. In particular, a feedback vaccination control law on the susceptible population is incorporated to stabilize the system dynamics. In [22], novel SEIR dynamic models (both integer-order and fractional-order) that include vaccine rate are proposed. In [23], a new SEIR epidemic model (the so-called SE(Is)(Ih)AR epidemic model) is presented. The conceived system incorporates two control laws, a feedback vaccination law and an antiviral treatment control law. In [24], some basic properties of a SEIR COVID-19 epidemic model subject to vaccination and treatment controls are studied. In particular, stability, boundedness, and non-negativity of the state trajectories are investigated in detail. In [25], the impact of multiple vaccination strategies on the dynamics of a fractional COVID-19 model is analyzed. The existence and uniqueness of the system solution is proven using Banach’s fixed point theorem. In [26], a set of ordinary differential equations that use vaccinations as the control input signals describe how COVID-19 behaves in Iranian and Russian societies. A Lyapunov-based methodology is used in [27] to examine how vaccination affects COVID-19 inhibition in Canada.

However, many of the proposed systems describe the virus accurately but use compartments such that it is difficult to find the initial conditions and some parameters. For example, the presence of a cabin for people who carry the disease but do not show symptoms increases the accuracy of the system, but after applying it, we find it is impossible to find the initial condition for this cabin, as well as the rate of settlement in it or the rate of exit from it. However, some systems are easy to apply, but they neglect many compartments. Therefore, in this work we formulated a system that took into account all possible compartments and whose initial conditions can be found to make the system applicable, as well as without neglecting any possible compartment that may reduce the accuracy of the system. The result was a system characterized by the maximum possible accuracy and applicability.

This paper introduces a new COVID-19 model, which incorporates the number of immunized individuals as an additional state variable describing the system dynamics, in an effort to further the discussion of the mathematical modeling of epidemics. The study demonstrates the existence of two fixed points, a disease-free fixed point and an endemic fixed point. The proposed compartment model is described by integer-order differential equations. When an inequality involving the vaccination rate is satisfied, the pandemic vanishes, according to the results of a stability analysis of the disease-free fixed point. This finding of the proposed approach is noteworthy because it is rigorous (proven by a theorem), which may aid decision-makers in better understanding the epidemiological behavior of the disease over time. The arrangement of the manuscript is as follows. A brand new compartment model for explaining the spread of the COVID-19 pandemic is presented in Section 2. Five state variables, namely the Susceptible class S, the Recovered class R, the Infection class I, the Infection dangerous class Id and the Vaccinated class V, which denotes the number of vaccinated people, are used to describe the dynamics of the system. The existence of the solution for the analyzed COVID-19 model is demonstrated in Section 3. It is demonstrated in Section 4 that the system has two fixed points, namely the epidemic fixed point and the disease-free fixed point. Additionally, the fundamental reproduction number is calculated. The stability of the disease-free fixed point is thoroughly examined in Section 5. The pandemic vanishes provided that an inequality involving the vaccinated rate and some system parameters are satisfied, according to a new theorem that ensures the global stability of the disease-free fixed point. Finally, simulation results are presented in Section 6 with the goal of demonstrating the value of the novel COVID-19 compartment model.

2. A New Compartment Model including Vaccinated Individuals

We divided the study population $\left(N\right)$ into three main classes in order to better understand the behavior of the epidemic’s spread: a class for those who were exposed to the infection, a class for those who became infected, a class for those who died as a result of the disease. Each of these classes are also divided into secondary classes, so that the class of people who are exposed to infection is divided into three sub-classes: people who are exposed to infection but were not previously infected and did not receive a vaccine S, people who were previously infected but recovered from the disease and are at risk of being infected again R, and people who were vaccinated against the epidemic V. As for the class of infected persons, it is divided into two secondary classes: for persons with good immunity and for whom infection does not pose a great risk, suppose that their ratio in society is $\lambda ,$$\left(\lambda \le 1\right),$ such that the class of infected people from this group are assigned to sub-class I. The infected persons for whom the infection is dangerous consists of the elderly, pregnant women and people with chronic diseases (whose ratios in society are $\left(1-\lambda \right)$) $I_d$. And finally, the class of deaths due to the epidemic D. Now we will explain the immigration in each class and compare it to the others.

Susceptible classS: This class acquires $\Omega$ persons in the unit of time. In the event that the studied area is isolated, then $\Omega$ represents the birth number. This class loses persons who are exposed to infection at a rate $r_1(I+I_d)$, where $r_1=\frac{p_{1}k}{N}$, k is the average numbers of contacts per person (per unit of time), $p_1$ is the probability of contagion and N is the total population (it can be considered as the maximum value of the population). Additionally, this class loses patients who have received a vaccination against the disease at a rate of $\upsilon$ and at a natural death rate of $\mu$.

Recovered classR: This class acquires newly recovered individuals from classes I and $I_d$ at a rate of $\rho$, and loses persons who are exposed to infection at a rate of $r_2(I+I_d)$, where $r_2=\frac{p_{2}k}{N}$, which is less than $r_1$. Similar to S, this class loses patients who have received a vaccination against the disease at a rate of $\upsilon$ and at a natural death rate of $\mu$.

Vaccinated classV: This class represents people who have been vaccinated and may have been infected and cured of the virus (they are people from the class R) or people who have never been infected (people from the class S). This category is not immune to the virus, but the infection rate is lower than in other classes. This class acquires persons at the vaccination rate $\upsilon$, representing newly vaccinated persons from classes S and R. This class loses persons who are infected at a rate of $r_3(I+I_d)$, where $r_3=\frac{p_{3}k}{N}$, which is less than $r_1$ and $r_2$ (the probability of infection $p_3$ in this class is less due to vaccination). This class also loses deceased persons (natural deaths) at a rate of $\mu$.

Infection classI: This class acquires newly infected persons from classes S, R and V at a rate of $\left(\lambda \left(r_1+r_2\right)+r_3\right)\left(I+I_d\right)$. This class loses recovered persons at a rate of $\rho$, as well as due to natural deaths at a rate of $\mu$ (we assume that in this class there are no deaths due to the epidemic).

The class of dangerous infection$I_d$: This class acquires newly infected persons from classes S and R (the class V does not go into this class because vaccination protects against dangerous infection) with a rate of $\left(1-\lambda \right)\left(r_1+r_2\right)\left(I+I_d\right)$. This class loses recovered persons at a rate of $\rho$, a rate of $\mu$ (natural deaths) and due to death related to infection at a rate of $\delta$.

Death classD: This class acquires persons who have died due to the epidemic at rate of $\delta .$

From these, we can now make the following model, which gives the mathematical explanation for all of the above:

$$\left\{\begin{array}{c}\frac{dS}{dt}=\Omega -r_1\left(I\left(t\right)+I_d\left(t\right)\right)S\left(t\right)-\left(\mu +\upsilon \right)S\left(t\right),\hfill \\ \frac{dR}{dt}=\rho \left(I\left(t\right)+I_d\left(t\right)\right)-r_2\left(I\left(t\right)+I_d\left(t\right)\right)R\left(t\right)-\left(\upsilon +\mu \right)R\left(t\right),\hfill \\ \frac{dV}{dt}=\upsilon \left(S\left(t\right)+R\left(t\right)\right)+{\rho }^{\prime }I\left(t\right)-r_3\left(I\left(t\right)+I_d\left(t\right)\right)V\left(t\right)-\mu V\left(t\right),\hfill \\ \frac{dI}{dt}=\left(\lambda \left(r_{1}S\left(t\right)+r_{2}R\left(t\right)\right)+r_{3}V\left(t\right)\right)\left(I\left(t\right)+I_d\left(t\right)\right)-\left(\mu +\rho +{\rho }^{\prime }\right)I\left(t\right),\hfill \\ \frac{d{I}_d}{dt}=\left(1-\lambda \right)\left(r_{1}S\left(t\right)+r_{2}R\left(t\right)\right)\left(I\left(t\right)+I_d\left(t\right)\right)-\left(\mu +\delta +\rho \right)I_d\left(t\right),\hfill \\ \frac{dD}{dt}=\delta I_d\left(t\right).\hfill \end{array}\right.t\in {\mathbb{R}}^{+}.$$

(1)

We note that the last equation in the system is isolated and can be neglected, so we will neglect it in future studies to simplify the system slightly. Thus, we find the following system:

$$\left\{\begin{array}{c}\frac{dS}{dt}=\Omega -r_1\left(I\left(t\right)+I_d\left(t\right)\right)S\left(t\right)-\left(\mu +\upsilon \right)S\left(t\right),\hfill \\ \frac{dR}{dt}=\rho \left(I\left(t\right)+I_d\left(t\right)\right)-r_2\left(I\left(t\right)+I_d\left(t\right)\right)R\left(t\right)-\left(\upsilon +\mu \right)R\left(t\right),\hfill \\ \frac{dV}{dt}=\upsilon \left(S\left(t\right)+R\left(t\right)\right)-r_3\left(I+I_d\right)V\left(t\right)-\mu V\left(t\right),\hfill \\ \frac{dI}{dt}=\left(\lambda \left(r_{1}S\left(t\right)+r_{2}R\left(t\right)\right)+r_{3}V\left(t\right)\right)\left(I\left(t\right)+I_d\left(t\right)\right)-\left(\mu +\rho \right)I\left(t\right),\hfill \\ \frac{d{I}_d}{dt}=\left(1-\lambda \right)\left(r_{1}S\left(t\right)+r_{2}R\left(t\right)\right)\left(I\left(t\right)+I_d\left(t\right)\right)-\left(\mu +\delta +\rho \right)I_d\left(t\right).\hfill \end{array}\right.t\in {\mathbb{R}}^{+}.$$

(2)

With the following initial conditions:

$$S\left(0\right),R\left(0\right),V\left(0\right),I\left(0\right),I_d\left(0\right)\ge 0.$$

(3)

The sum of all the compounds represents the total number of living persons in the studied population $\left(N\right)$, thus:

$$N\left(t\right)=S\left(t\right)+R\left(t\right)+V\left(t\right)+I\left(t\right)+I_d\left(t\right).$$

(4)

This system must be subject to the following assumptions:

 Assumption A1. 

The number of new infections is directly proportional to infection rates $r_1$, $r_2$ and $r_3$.

 Assumption A2. 

The number of new infections is inversely proportional to the cure rate ρ and the vaccination rate υ.

3. Existence, Positivity and Invariant Region

3.1. Existence and Uniqueness

To prove the existence and uniqueness of the solution, we use the Cauchy–Lipshitz theorem. First, let the problem

$$\frac{dy}{dt}=f(t,y\left(t\right)),$$

(5)

where $f:U\to {\mathbb{R}}^n,$ is continuous and U is an open set of $\mathbb{R}\times {\mathbb{R}}^n.$

 Theorem 1 

(Cauchy–Lipschitz [28]). Suppose that f is locally Lipschitz in y, then problem (5) with the initial conditions $(t_0,y_0)\in \mathbb{R}\times {\mathbb{R}}^n$ implies a unique maximal solution.

Based on the above theorem, we find the following result regarding the existence and uniqueness of the solution to problem (2).

 Theorem 2. 

System (2) with the initial conditions (3) has a unique maximal solution.

 Proof. 

The second part of system (2) is continuous and also belongs to class $C^{\infty }$, with respect to $(S,R,V,I,I_d)$, hence it is a locally Lipschitzien. From it, according to the Cauchy–Lipschitz theorem system (2) has a unique maximal solution. □

3.2. Positivity and Invariant Region

We are interested here in studying whether the solution is positive. Because it must be, we will also find an invariant region such that when we take the initial conditions from it, the solution remains bonded and belongs to this region. First, we start with the positivity of the solution through the following Lemma.

 Lemma 1. 

The solution of (2) is positive when the initial conditions are positives.

 Proof. 

First, we note that the initial condition is positive and the solution is continuous. Therefore, in order for one of the components of the solution to become negative, it must first be null and the derivative at the zero point be negative. But we note from:

$$\left\{\begin{array}{c}{\left.\frac{dS}{dt}\right|}_{S=0}=\Omega ,\hfill \\ {\left.\frac{dR}{dt}\right|}_{R=0}=\rho \left(I+I_d\right),\hfill \\ {\left.\frac{dV}{dt}\right|}_{V=0}=\upsilon \left(S+R\right),\hfill \\ {\left.\frac{dI}{dt}\right|}_{I=0}=\left(\lambda \left(r_{1}S+r_{2}R\right)+r_{3}V\right)\left(I_d\right),\hfill \\ {\left.\frac{d{I}_d}{dt}\right|}_{I_d=0}=\left(1-\lambda \right)\left(r_{1}S+r_{2}R\right)I.\hfill \end{array}\right.$$

(6)

Hence, the first element that vanishes from $\left(S,R,V,I,I_d\right)$, its derivative at the point of vanishing is non-negative and therefore increases again and remains positive. Therefore, all elements of $\left(S,R,V,I,I_d\right)$ remain positive and never become negative as long as the initial condition is positive, i.e., $\left(S,R,V,I,I_d\right)\in {\mathbb{R}}_{+}^5$, where ${\mathbb{R}}_{+}^5=\{\left(x_1,x_2,x_3,x_4,x_5\right)\in {\mathbb{R}}^5$ and $x_i\ge 0,\mathrm{for}i=1.5\}$. □

Based on the following theorem:

 Theorem 3 

([29] Comparison Theorem). Let $f;g:\mathbb{R}\to \mathbb{R}$ two Lipschitz function. We consider the solutions $x\left(t\right)$ and $y\left(t\right)$ of the Cauchy problems:

$$\begin{array}{cc}\left\{\begin{array}{c}{x}^{\prime }\left(t\right)=f(t,x\left(t\right),\hfill \\ x\left(0\right)=x_0,\hfill \end{array}\right.\hfill & \left\{\begin{array}{c}{y}^{\prime }\left(t\right)=g(t,y\left(t\right),\hfill \\ y\left(0\right)=y_0.\hfill \end{array}\right.\end{array}$$

Suppose that $f(t;x)\le g(t;x)$ for all $(t;x)\in \mathbb{R}\times \mathbb{R}$ and that $x_0\le y_0$. Then $x\left(t\right)\le y\left(t\right)$ for all t.

This immediately follows for the invariant region, and we give the following result.

 Theorem 4. 

System (2) has

$$\Psi =\left\{\left(S,R,V,I,I_d\right)\in {\mathbb{R}}_{+}^{5}andS+R+V+I+I_d\le \frac{\Omega }{\mu }\right\},$$

(7)

as invariant region.

 Proof. 

Adding the equations of the system (2), we get:

$$\frac{dS\left(t\right)}{dt}+\frac{dR\left(t\right)}{dt}+\frac{dV\left(t\right)}{dt}+\frac{dI\left(t\right)}{dt}+\frac{d{I}_d\left(t\right)}{dt}=\Omega -\mu N\left(t\right)-\delta I_d\left(t\right),$$

which means

$$\frac{dN\left(t\right)}{dt}=\Omega -\mu N\left(t\right)-\delta I_d\left(t\right),$$

Since $I_d$ is positive, we get:

$$\frac{dN\left(t\right)}{dt}\le \Omega -\mu N\left(t\right),$$

where, according to Equation (4)

$$N\left(0\right)=S\left(0\right)+R\left(0\right)+V\left(0\right)+I\left(0\right)+I_d\left(0\right).$$

Applying the comparison theorem by replacing $x=N,x\left(0\right)=y\left(0\right)=N\left(0\right),f=\frac{dN}{dt}$ and $g=\Omega -\mu N\left(t\right)$, we find that

$$N\left(t\right)\le \frac{\Omega }{\mu }-\left(\frac{\Omega }{\mu }-N\left(0\right)\right)e^{-\mu t}.$$

So the following is achieved

$$0\le N\left(t\right)\le \frac{\Omega }{\mu },$$

where $N\left(0\right)\le \frac{\Omega }{\mu }.$

Therefore, the solution to the equation belongs to $\Psi$. □

4. Fixed Points and Basic Reproduction Number

In this section, we will study the fixed points and the basic reproduction number, which are important for studying the stability of epidemic systems.

4.1. Fixed Points

Finding the fixed points is necessary before studying the dynamics of the model (2), and finding the fixed points requires solving the equation:

$$\left\{\begin{array}{c}\Omega -r_1\left(I^{*}+I_d^{*}\right)S-\left(\mu +\upsilon \right)S^{*}=0,\hfill \\ \rho \left(I^{*}+I_d^{*}\right)-r_2\left(I^{*}+I_d^{*}\right)R^{*}-\left(\upsilon +\mu \right)R^{*}=0,\hfill \\ \upsilon \left(S^{*}+R^{*}\right)-r_3\left(I^{*}+I_d^{*}\right)V^{*}-\mu V^{*}=0,\hfill \\ \left(\lambda \left(r_1{S}^{*}+r_2{R}^{*}\right)+r_3{V}^{*}\right)\left(I^{*}+I_d^{*}\right)-\left(\mu +\rho \right)I^{*}=0\hfill \\ \left(1-\lambda \right)\left(r_1{S}^{*}+r_2{R}^{*}\right)\left(I^{*}+I_d^{*}\right)-\left(\mu +\delta +\rho \right)I_d^{*}=0.\hfill \end{array}\right.$$

(8)

The previous equation has the point $E_0=\left(\frac{\Omega }{\left(\mu +\upsilon \right)},0,\frac{\upsilon \Omega }{\mu \left(\mu +\upsilon \right)},0,0\right)$ as a solution. As can be seen, there is no disease at this point, so it is referred to as the disease-free fixed point, and we will study its stability later.

If we suppose that $\left(I^{*}+I_d^{*}\right)\ne 0$, we will get:

$$\begin{array}{c}\frac{\Omega }{r_1\left(I^{*}+I_d^{*}\right)+\left(\mu +\upsilon \right)}=S^{*},\\ \frac{\rho \left(I^{*}+I_d^{*}\right)}{r_2\left(I^{*}+I_d^{*}\right)+\left(\upsilon +\mu \right)}=R^{*},\\ \frac{\upsilon \left(S^{*}+R^{*}\right)}{r_3\left(I^{*}+I_d^{*}\right)+\mu }=V^{*},\\ \frac{\left(\lambda \left(r_1{S}^{*}+r_3{R}^{*}\right)+r_2{V}^{*}\right)}{\left(\mu +\rho \right)}=\frac{I^{*}}{\left(I^{*}+I_d^{*}\right)},\\ \frac{\left(1-\lambda \right)\left(r_1{S}^{*}+r_2{R}^{*}\right)}{\left(\mu +\delta +\rho \right)}=\frac{I_d^{*}}{\left(I^{*}+I_d^{*}\right)}.\end{array}$$

(9)

This system is a classical non-linear system, which can be solved numerically and then studied for its stability. Overall, this point is called the endemic equilibrium point $E^{*}=\left(S^{*},R^{*},V^{*},I^{*},I_d^{*}\right).$

4.2. Basic Reproduction Number

We will now calculate a key number called the basic reproduction number $R_0$. In the study of stability for the disease-free fixed point, this is crucial. We will follow the steps described in [30] to calculate this number, which represents the rate of new people being infected by one sick person until their recovery. We determine the class expressing the new infection, which corresponds to the final two equations in system (2), and hence it can be calculated using:

$$\begin{array}{c}\frac{dI}{dt}=\left(\lambda \left(r_{1}S\left(t\right)+r_{2}R\left(t\right)\right)+r_{3}V\left(t\right)\right)\left(I\left(t\right)+I_d\left(t\right)\right)-\left(\mu +\rho \right)I\left(t\right),\hfill \\ \frac{d{I}_d}{dt}=\left(1-\lambda \right)\left(r_{1}S\left(t\right)+r_{2}R\left(t\right)\right)\left(I\left(t\right)+I_d\left(t\right)\right)-\left(\mu +\delta +\rho \right)I_d\left(t\right).\hfill \end{array}$$

(10)

We reformulate this system as follows

$$\frac{d}{dt}\left(\begin{array}{c}I\\ I_d\end{array}\right)=\mathcal{F}-\mathcal{V},$$

where $\mathcal{F}$ is the rate of the appearance of new infections,

$$\mathcal{F}=\left(\begin{array}{c}\left(\lambda \left(r_{1}S\left(t\right)+r_{2}R\left(t\right)\right)+r_{3}V\left(t\right)\right)\left(I\left(t\right)+I_d\left(t\right)\right)\\ \left(1-\lambda \right)\left(r_{1}S\left(t\right)+r_{2}R\left(t\right)\right)\left(I\left(t\right)+I_d\left(t\right)\right)\end{array}\right),$$

and $\mathcal{V}$ is the rate of the transfer of individuals into other compartments,

$$\mathcal{V}=\left(\begin{array}{c}\left(\mu +\rho \right)I\left(t\right)\\ \left(\mu +\delta +\rho \right)I_d\left(t\right)\end{array}\right).$$

Calculating the Jacobian matrix F and V for $\mathcal{F}$ and $\mathcal{V}$ respectively at the disease-free fixed point, i.e., at $\left(S=\frac{\Omega }{\left(\mu +\upsilon \right)},R=0,V=\frac{\upsilon \Omega }{\mu \left(\mu +\upsilon \right)},I=0,I_d=0\right)$, we get:

$$F=\left(\begin{array}{cc}\frac{\Omega }{\left(\mu +\upsilon \right)}\left(\lambda r_1+\frac{\upsilon r_3}{\mu }\right)& \frac{\Omega }{\left(\mu +\upsilon \right)}\left(\lambda r_1+\frac{\upsilon r_3}{\mu }\right)\\ \left(1-\lambda \right)\left(\frac{\Omega r_1}{\left(\mu +\upsilon \right)}\right)& \left(1-\lambda \right)\left(\frac{\Omega r_1}{\left(\mu +\upsilon \right)}\right)\end{array}\right),$$

and

$$V=\left(\begin{array}{cc}\left(\mu +\rho \right)& 0\\ 0& \left(\mu +\delta +\rho \right)\end{array}\right).$$

The next generation matrix is:

$$F{V}^{-1}=\left(\begin{array}{cc}\frac{\Omega }{\mu }\frac{\upsilon r_3+\lambda \mu r_1}{\left(\mu +\upsilon \right)\left(\mu +\rho \right)}& \frac{\Omega }{\mu }\frac{\upsilon r_3+\lambda \mu r_1}{\left(\mu +\upsilon \right)\left(\mu +\delta +\rho \right)}\\ -\Omega r_1\frac{\lambda -1}{\left(\mu +\upsilon \right)\left(\mu +\rho \right)}& -\Omega r_1\frac{\lambda -1}{\left(\mu +\upsilon \right)\left(\mu +\delta +\rho \right)}\end{array}\right).$$

(11)

The basic reproductive number is given as the spectral radius (the greater length of the eigenvalues) of $F{V}^{-1}$. When calculating the characteristic polynomial $P\left(F{V}^{-1}\right),$ for $F{V}^{-1}$, we find

$$P\left(F{V}^{-1}\right)=X\left(X-\frac{\Omega }{\left(\mu +\upsilon \right)}\left(\frac{\left(1-\lambda \right)r_1}{\left(\mu +\delta +\rho \right)}+\frac{\upsilon r_3+\lambda \mu r_1}{\mu \left(\mu +\rho \right)}\right)\right)$$

We notice that $F{V}^{-1}$ has 0 and $\frac{\Omega }{\left(\mu +\upsilon \right)}\left(\frac{\left(1-\lambda \right)r_1}{\left(\mu +\delta +\rho \right)}+\frac{\upsilon r_3+\lambda \mu r_1}{\mu \left(\mu +\rho \right)}\right)$ as eigenvalues, thus

$$R_0=\frac{\Omega }{\left(\mu +\upsilon \right)}\left(\frac{\left(1-\lambda \right)r_1}{\left(\mu +\delta +\rho \right)}+\frac{\upsilon r_3+\lambda \mu r_1}{\mu \left(\mu +\rho \right)}\right).$$

(12)

5. Stability Analysis of the Disease-Free Fixed Point

What matters to us is the disappearance of the disease. In this section, we will formulate the conditions to ensure the stability of the disease-free fixed point, that is, formulate the conditions to ensure the disappearance of the disease. When applying the model (2) in a specific region, we need the initial conditions and the parameters of that region. The parameters are divided into two parts. Firstly, the fixed part, which cannot be modified, is represented by: the birth rate $\Omega$, natural death rate $\mu$, recovered rate $\rho$ and death due to infection rate $\delta$. Although it is possible to modify $\rho$ and $\delta$ by improving the health conditions, we assume that the authorities are doing everything possible, meaning that these rates are the best they can be. The second kind of parameters can be modified, which are: infection rates $r_1$, $r_2$ and $r_3$ and the vaccinated rate $\upsilon$. Authorities can impose measures such as closing some facilities and imposing quarantine to reduce the rate of infection that reduces rates $r_1$, $r_2$ and $r_3$. By increasing the vaccination rate, we can decrease classes S and R and increase class V, in which the probability of infection is less than the previous two classes.

5.1. Local Stability

The stability of the disease-free fixed point will be examined in this subsection using an imposed condition on $R_0.$

 Theorem 5. 

Assume that $R_0<1$. Therefore, the disease-free fixed point $E_0$ of system (2) is locally asymptotically stable.

 Proof. 

The system at $E_0$ has the Jacobian matrix:

$$J=\left(\begin{array}{ccccc}-\left(\mu +\upsilon \right)& 0& 0& -\frac{\Omega r_1}{\left(\mu +\upsilon \right)}& -\frac{\Omega r_1}{\left(\mu +\upsilon \right)}\\ 0& -\left(\upsilon +\mu \right)& 0& \rho & \rho \\ \upsilon & \upsilon & -\mu & -\frac{\Omega \upsilon r_3}{\mu \left(\mu +\upsilon \right)}& -\frac{\Omega \upsilon r_3}{\mu \left(\mu +\upsilon \right)}\\ 0& 0& 0& \left(\frac{\lambda r_1\Omega }{\left(\mu +\upsilon \right)}+\frac{\upsilon r_3\Omega }{\mu \left(\mu +\upsilon \right)}\right)-\left(\mu +\rho \right)& \left(\frac{\lambda r_1\Omega }{\left(\mu +\upsilon \right)}+\frac{\upsilon r_3\Omega }{\mu \left(\mu +\upsilon \right)}\right)\\ 0& 0& 0& \frac{\left(1-\lambda \right)r_1\Omega }{\left(\mu +\upsilon \right)}& \frac{\left(1-\lambda \right)r_1\Omega }{\left(\mu +\upsilon \right)}-\left(\mu +\delta +\rho \right)\end{array}\right)$$

(13)

Then, the characteristic polynomial:

$$P\left(J\right)=\left(X+\mu \right){\left(X+\mu +\upsilon \right)}^2\left(X^2+AX+B\right)$$

where

$$\begin{array}{c}A=\left(\mu +\rho \right)+\left(\mu +\delta +\rho \right)-\frac{\Omega \left(\mu r_1+\upsilon r_3\right)}{\mu \left(\mu +\upsilon \right)},\\ B=\left(\mu +\rho \right)\left(\mu +\delta +\rho \right)-\frac{\Omega \left(\mu r_1\left(\mu +\rho +\lambda \delta \right)+\upsilon r_3\left(\mu +\delta +\rho \right)\right)}{\mu \left(\mu +\upsilon \right)},\end{array}$$

As a result, the matrix J has two eigenvalues that are both less than zero: $-\mu$ as a normal eigenvalue and $-\left(\mu +\upsilon \right)$ as a double eigenvalue. The rest of the roots of the polynomial:

$$X^2+AX+B,$$

(14)

A and B can be written as follows:

$$\begin{array}{c}A=\left(\mu +\rho \right)+\left(\mu +\delta +\rho \right)\left(1-R_0\right)+\frac{\Omega \delta \left(\upsilon r_3+\lambda \mu r_1\right)}{\mu \left(\mu +\rho \right)\left(\mu +\upsilon \right)},\\ B=\left(\mu +\rho \right)\left(\mu +\delta +\rho \right)\left(1-R_0\right).\end{array}$$

If $R_0<1,$ then according to the Routh–Hurwitz criterion, the roots of the polynomial (14) have a negative real part. Therefore, all eigenvalues of matrix J have a negative real part, so the disease-free fixed point is locally asymptotically stable. □

This theorem can be applied easily, but it remains only local, that is, when the initial conditions are far from the fixed point. Hence, we need to study its global stability.

5.2. Global Stability

The last two equations of the system added together describe infection, and yield the following equation:

$$\frac{d}{dt}\left(I+I_d\right)=\left(r_{1}S\left(t\right)+r_{2}R\left(t\right)+r_{3}V\left(t\right)-\left(\mu +\rho \right)\right)\left(I\left(t\right)+I_d\left(t\right)\right)-\delta I_d\left(t\right).$$

Since I is positive:

$$\frac{d}{dt}\left(I+I_d\right)\le \left(r_{1}S\left(t\right)+r_{2}R\left(t\right)+r_{3}V\left(t\right)-\left(\mu +\rho \right)\right)\left(I\left(t\right)+I_d\left(t\right)\right).$$

Since $r_i=\frac{p_{i}k}{N},i=1,2,3$:

$$\begin{array}{cc}\frac{d}{dt}\left(I+I_d\right)\hfill & \le \left(\frac{p_{1}k}{N}S\left(t\right)+\frac{p_{2}k}{N}R\left(t\right)+\frac{p_{3}k}{N}V\left(t\right)-\left(\mu +\rho \right)\right)\left(I\left(t\right)+I_d\left(t\right)\right),\hfill \\ & \le \left(\frac{p_{1}k}{N}\left(S\left(t\right)+R\left(t\right)+V\left(t\right)\right)-\left(\mu +\rho \right)\right)\left(I\left(t\right)+I_d\left(t\right)\right),\hfill \\ & \le \left(p_{1}k-\left(\mu +\rho \right)\right)\left(I\left(t\right)+I_d\left(t\right)\right).\hfill \end{array}$$

Using the comparison theorem, we find that

$$\left(I+I_d\right)\left(t\right)\le \left(I+I_d\right)\left(0\right)e^{\left(p_{1}k-\left(\mu +\rho \right)\right)t}.$$

We note that if $\left(p_{1}k-\left(\mu +\rho \right)\right)<0,$ then $\left(I+I_d\right)\left(t\right)\to 0$ when $t\to \infty .$ Thus, we get the following result:

 Theorem 6. 

If

$$k<\frac{\mu +\rho }{p_1},$$

(15)

then the disease will disappear.

 Remark 1. 

We note that this condition is very logical. Factors contributing to the disappearance of the disease are when the rate of infection is small, or the recovery rate is increased, or the probability of injury is small.

This condition is sufficient for the disappearance of the disease, but it may not be achievable and may be very expensive. Sometimes it is not possible to reduce the rate of infection more than a certain limit. In the following, we will study the disappearance of the disease by studying the stability of the disease-free fixed point.

To prove global stability we use the method of Carlos Castillo-Chavez described in [31]. First, system (2) must be written in the form:

$$\begin{array}{c}\frac{dX}{dt}=F(X,Y),\\ \frac{dY}{dt}=G(X,Y),\end{array}$$

(16)

where $X=\left(S,R,V\right)$ denotes the number of uninfected individuals and $Y=\left(I,I_d\right)$ denotes the number of infected, infectious. $E_0=(X_0,0)$ denotes the disease-free equilibrium of this system.

$\left(H_1\right)$ In $\frac{dX}{dt}=F(X,0),$$X_0$ is globally asymptotically stable.

$\left(H_2\right)$$G(X,Y)=AY-\widehat{G}(X,Y),\widehat{G}(X,Y)\ge 0$ for $(X,Y)\in \Psi ,$ where A is an M-matrix (the off diagonal elements of A are non negative) and $\Psi$ is the region where the model makes biological sense.

If system (2) satisfies the above two conditions, then the following theorem holds:

 Theorem 7 

([31]). $E_0=(X_0,0),$ is a globally asymptotically stable equilibrium of (2) provided that $R_0<1$ and that assumptions $\left(H_1\right)$ and $\left(H_2\right)$ are satisfied.

First, we will prove that the two conditions $\left(H_1\right)$ and $\left(H_2\right)$ are satisfied for system (2): the system $\frac{dX}{dt}=F(X,0)$ is written as follows:

$$\begin{array}{c}\frac{dS}{dt}=\Omega -\left(\mu +\upsilon \right)S\left(t\right),\hfill \\ \frac{dR}{dt}=-\left(\mu +\upsilon \right)R\left(t\right),\hfill \\ \frac{dV}{dt}=\upsilon \left(S\left(t\right)+R\left(t\right)\right)-\mu V\left(t\right),\hfill \end{array}$$

and $X_0=\left(\frac{\Omega }{\left(\mu +\upsilon \right)},0,\frac{\upsilon \Omega }{\mu \left(\mu +\upsilon \right)}\right)$, and from it

$$S\left(t\right)=\frac{\Omega }{\left(\mu +\upsilon \right)}-\left(\frac{\Omega }{\left(\mu +\upsilon \right)}-S\left(0\right)\right)e^{-\left(\mu +\upsilon \right)t},$$

and

$$R\left(t\right)=R\left(0\right)e^{-\left(\mu +\upsilon \right)t}.$$

Additionally, from

$$\frac{d\left(S+R+V\right)}{dt}=\Omega -\mu \left(S+R+V\right)\left(t\right),$$

we get

$$\left(S+R+V\right)\left(t\right)=\frac{\Omega }{\mu }-\left(\frac{\Omega }{\mu }-\left(S+R+V\right)\left(0\right)\right)e^{-\upsilon t},$$

when $t\to \infty :$$S\left(t\right)\to \frac{\Omega }{\left(\mu +\upsilon \right)},R\left(t\right)\to 0,V\left(t\right)\to \frac{\upsilon \Omega }{\mu \left(\mu +\upsilon \right)}.$ Hence, the global stability of $X_0.$

We have, on the other hand:

$$G(X,Y)=\left(\begin{array}{c}\left(\lambda \left(r_{1}S\left(t\right)+r_{2}R\left(t\right)\right)+r_{3}V\left(t\right)\right)\left(I\left(t\right)+I_d\left(t\right)\right)-\left(\mu +\rho \right)I\left(t\right)\hfill \\ \left(1-\lambda \right)\left(r_{1}S\left(t\right)+r_{2}R\left(t\right)\right)\left(I\left(t\right)+I_d\left(t\right)\right)-\left(\mu +\delta +\rho \right)I_d\left(t\right)\hfill \end{array}\right),$$

thus

$$G(X,Y)=\widehat{A}Y-\widehat{G}(X,Y),$$

where

$$\widehat{A}=\left(\begin{array}{cc}\left(\frac{\lambda r_1\Omega }{\left(\mu +\upsilon \right)}+\frac{\upsilon r_3\Omega }{\mu \left(\mu +\upsilon \right)}\right)-\left(\mu +\rho \right)& \left(\frac{\lambda r_1\Omega }{\left(\mu +\upsilon \right)}+\frac{\upsilon r_3\Omega }{\mu \left(\mu +\upsilon \right)}\right)\\ \frac{\left(1-\lambda \right)r_1\Omega }{\left(\mu +\upsilon \right)}& \frac{\left(1-\lambda \right)r_1\Omega }{\left(\mu +\upsilon \right)}-\left(\mu +\delta +\rho \right)\end{array}\right)$$

and

$$\widehat{G}(X,Y)=\left(\begin{array}{c}\left(\left(\frac{\lambda r_1\Omega }{\left(\mu +\upsilon \right)}+\frac{\upsilon r_3\Omega }{\mu \left(\mu +\upsilon \right)}\right)-\lambda \left(\left(r_{1}S\left(t\right)+r_{2}R\left(t\right)\right)+r_{3}V\left(t\right)\right)\right)\left(I+I_d\right)\\ \left(1-\lambda \right)\left(\frac{r_1\Omega }{\left(\mu +\upsilon \right)}-\left(r_{1}S\left(t\right)+r_{2}R\left(t\right)\right)\right)\left(I+I_d\right)\end{array}\right)$$

$\widehat{A}$ is an M-matrix, thus we get the following result:

 Theorem 8. 

Suppose that $R_0<1$. If

$$\rho -\frac{r_2\Omega }{\left(\mu +\upsilon \right)}\le 0,$$

(17)

then $E_0$ is globally asymptotically stable.

 Proof. 

Suppose the initial conditions are as follows: $S\left(0\right)+R\left(0\right)\le \frac{\Omega }{\left(\mu +\upsilon \right)},V\left(0\right)\le \frac{\upsilon \Omega }{\mu \left(\mu +\upsilon \right)},$ such that

$$\left(\frac{r_1\Omega }{\left(\mu +\upsilon \right)}-\left(r_{1}S\left(t\right)+r_{2}R\left(t\right)\right)\right)\ge r_1\left(\frac{\Omega }{\left(\mu +\upsilon \right)}-\left(S\left(t\right)+R\left(t\right)\right)\right).$$

On the other hand:

$$\begin{array}{cc}\frac{d\left(S\left(t\right)+R\left(t\right)\right)}{dt}=\hfill & \Omega +\rho \left(I\left(t\right)+I_d\left(t\right)\right)-r_1\left(I\left(t\right)+I_d\left(t\right)\right)S\left(t\right)-r_2\left(I\left(t\right)+I_d\left(t\right)\right)R\left(t\right)\hfill \\ & -\left(\upsilon +\mu \right)\left(R\left(t\right)+S\left(t\right)\right),\hfill \end{array}$$

which gives us

$$\begin{array}{ccc}{\left.\frac{d\left(S\left(t\right)+R\left(t\right)\right)}{dt}\right|}_{S+R=\frac{\Omega }{\left(\mu +\upsilon \right)}}\hfill & =\hfill & \Omega +\rho \left(I\left(t\right)+I_d\left(t\right)\right)-r_1\left(I\left(t\right)+I_d\left(t\right)\right)S\left(t\right)-r_2\left(I\left(t\right)+I_d\left(t\right)\right)R\left(t\right)\hfill \\ & & -\left(\upsilon +\mu \right)\left(R\left(t\right)+S\left(t\right)\right)\hfill \\ & \le \hfill & \Omega +\rho \left(I\left(t\right)+I_d\left(t\right)\right)-\left(\left(\upsilon +\mu \right)+r_2\left(I\left(t\right)+I_d\left(t\right)\right)\right)\left(R\left(t\right)+S\left(t\right)\right)\hfill \\ & \le \hfill & \left(\rho -\frac{r_2\Omega }{\left(\mu +\upsilon \right)}\right)\left(I\left(t\right)+I_d\left(t\right)\right),\hfill \end{array}$$

and according to (15):

$${\left.\frac{d\left(S\left(t\right)+R\left(t\right)\right)}{dt}\right|}_{S+R=\frac{\Omega }{\left(\mu +\upsilon \right)}}\le 0,$$

thus

$$\left(S\left(t\right)+R\left(t\right)\right)\le \frac{\Omega }{\left(\mu +\upsilon \right)},$$

(18)

and

$$\left(1-\lambda \right)\left(\frac{r_1\Omega }{\left(\mu +\upsilon \right)}-\left(r_{1}S\left(t\right)+r_{2}R\left(t\right)\right)\right)\left(I+I_d\right)\ge 0.$$

On the other hand:

$${\left.\frac{dV}{dt}\right|}_{V=\frac{\upsilon \Omega }{\mu \left(\mu +\upsilon \right)}}=\upsilon \left(S\left(t\right)+R\left(t\right)\right)-r_3\left(I+I_d\right)\frac{\upsilon \Omega }{\mu \left(\mu +\upsilon \right)}-\mu \frac{\upsilon \Omega }{\mu \left(\mu +\upsilon \right)},$$

and from (18)

$${\left.\frac{dV}{dt}\right|}_{V=\frac{\upsilon \Omega }{\mu \left(\mu +\upsilon \right)}}\le -r_3\left(I+I_d\right)\frac{\upsilon \Omega }{\mu \left(\mu +\upsilon \right)}\le 0,$$

then $V\left(t\right)\le \frac{\upsilon \Omega }{\mu \left(\mu +\upsilon \right)},$ and from it

$$\left(\left(\frac{\lambda r_1\Omega }{\left(\mu +\upsilon \right)}+\frac{\upsilon r_3\Omega }{\mu \left(\mu +\upsilon \right)}\right)-\lambda \left(\left(r_{1}S\left(t\right)+r_{2}R\left(t\right)\right)+r_{3}V\left(t\right)\right)\right)\left(I+I_d\right)\ge 0.$$

Finally, $\widehat{G}(X,Y)\ge 0.$ According to Carlos Castillo-Chavez [31], $E_0$ is globally asymptotically stable. □

6. Numerical Simulations

In this section, we’ll apply the system under study to the country of Brazil and contrast the results with real data to determine how effective it is. Based on [32], we can divide the initial population as follows:

$$\begin{array}{ccccc}S\left(0\right)=17300532,& R\left(0\right)=30921318,& V\left(0\right)=169017000,& I\left(0\right)=616214,& I_d\left(0\right)=410809.\end{array}$$

(19)

The values of the system’s parameters can be calculated according to the same source [32], however, we discover the following:

$$\begin{array}{ccc}\Omega =287010;& \mu =0.00066;& \lambda =0.6;\\ r_1=4.4\times {10}^{-11};& r_2=3.215\times {10}^{-11};& r_3=1.4925\times {10}^{-11};\\ \rho =0.0352;& \upsilon =0.0166;& \delta =5.6259\times {10}^{-4}.\end{array}$$

(20)

We also take the real data of active cases in Brazil in the period from 23 July to 15 August 2022, shown in Figure 1.

We apply system (2) using the previous data, and we find the numerical simulation as shown in Figure 2.

Note that:

$$\rho -\frac{r_2\Omega }{\left(\mu +\upsilon \right)}=3.4665\times {10}^{-2}>0,\mathrm{and}{R}_0=0.19435.$$

(21)

Thus, the condition of Theorem 7 is not fulfilled, we can then say that $E_0$ is locally asymptotically stable (according to Theorem 7) but not necessarily globally asymptotically stable. We notice from the simulations in Figure 3 that the model gives a good result in terms of its prediction.

7. Conclusions

The number of immunized people is a new state variable in the COVID-19 compartment model that this paper has used to illustrate the dynamics of the system. The proposed model, which is described by integer-order differential equations and is demonstrated in this paper, has a disease-free fixed point and an endemic fixed point. We also proved a new theorem, which has assured the global stability of the disease-free fixed point. The theorem has highlighted that the pandemic can disappear, provided that an inequality involving the vaccination rate is satisfied. This is a noteworthy discovery of the suggested approach, as it may aid in the comprehension of decision-makers when it comes to the disease’s evolving epidemiological behavior. Last but not least, numerical simulations have been conducted to demonstrate the value of the COVID-19 epidemic model.