Stability Analysis of an Extended SEIR COVID-19 Fractional Model with Vaccination Efficiency.

This work is aimed at presenting a new numerical scheme for COVID-19 epidemic model based on Atangana-Baleanu fractional order derivative in Caputo sense (ABC) to investigate the vaccine efficiency. Our construction of the model is based on the classical SEIR, four compartmental models with an additional compartment V of vaccinated people extending it SEIRV model, for the transmission as well as an effort to cure this infectious disease. The point of disease-free equilibrium is calculated, and the stability analysis of the equilibrium point using the reproduction number is performed. The endemic equilibrium's existence and uniqueness are investigated. For the solution of the nonlinear system presented in the model at different fractional orders, a new numerical scheme based on modified Simpson's 1/3 method is developed. Convergence and stability of the numerical scheme are thoroughly analyzed. We attempted to develop an epidemiological model presenting the COVID-19 dynamics in Italy. The proposed model's dynamics are graphically interpreted to observe the effect of vaccination by altering the vaccination rate.

1. Introduction

Emerging in December 2019, COVID-19 became the biggest global threat. The whole world faced a critical situation, and a panic was created throughout the globe. The World Health Organization (WHO) declared it as an international health emergency. Not much was known about this virus as it was a new virus that has not been identified earlier in humans. So, relying on the observations and symptoms shown by the people who caught the virus, strategies of isolation, lockdown, and social distancing were adopted to prevent the disease from spreading. With the passage of time, the virus becomes more threatening because of its new variants. To investigate the dynamics and transmission of the virus, researchers utilized the concept of mathematical modeling [1-4]. These models actually do not provide the cure for the infectious diseases but can contribute in predicting the dynamics and behaviour of these infections through simulations.

Researchers rely more on fractional models because of their nonlocal behaviour, hereditary properties, and memory effects rather than integer order models. Besides, a number of experimental facts exhibit that the natural dynamics go along with fractional calculus. Ahmad et al. [5] utilized the modified Euler's method to study the dynamics of the model for COVID-19 based on fractional differential equations and simulate the results by using the available data of initial days of the spread in Wuhan City. Zakary et al. [6] investigated the situation of the spread of the Corona Virus under the effect of two classes, quarantined and other who did not respect the quarantine strategy using the mathematical model, and an optimal control strategy to reduce the infections in Morocco. Ahmad et al. [7] presented a brief analysis on the WHO reported data for Pakistan and investigated it using the SEIR fractional model. Askar et al. [8] considered the SITR fractional model to forecast the transmission of deadly virus along with imposed lockdown in India. Bushnaq et al. [9] investigated the effect of physical distance to control the spread of the virus through fractional order model. Abdulwasaa et al. [10] proposed a fractal fractional model for forecasting deaths and new cases of COVID-19 outbreak in India. Baba and Rihan [11] used the Caputo-Fabrizio fractional order model to analyze the dynamics of COVID-19 variants. A fractional order delay differential model for COVID-19 infection was provided by Rihan and Gandhi [12] to understand what triggers the severity of symptoms and disease of polluted lung and respiratory system. A stochastic SIAQR epidemic model was proposed by Rihan and Alsakaji [13] to understand the dynamics of COVID-19. Furati et al. [14] tested the data taken from China, Saudi Arabia, Brazil, and Italy by proposing a numerical scheme using fractional order model performing the two-step generalized exponential time-differencing method. Alqahtani [15] numerically analyzed and investigated the stability of the SIR epidemic model considering the health system. Recently another fractional order derivative with nonsingular kernel named as Generalized Hattaf Fractional (GHF) derivative [16] is introduced. We can see the use of this derivative for epidemiological models in recent literature. By creating a new numerical method, Hattaf [17] investigated the qualitative properties of solutions to fractional differential equations with the new GHFs derivative, including stability, asymptotic stability, and Mittag-Leffler stability and later applied the obtained analytical and numerical results to a biological nonlinear system derived from epidemiology. Cheneke et al. [18] developed and examined a coinfection of HIV and cholera model using GHF derivative and the solution's behaviour is deciphered. In present work we have used ABC derivative for our proposed model. However, we can extend our work using the GHF derivative for fractional order mathematical models in future.

Despite all the measures taken to control the virus, COVID-19 remained a threat because of its high propagation rate leaving the majority of the population susceptible to the disease and escalating the need to develop a suitable and effective vaccine. The process to gain immunity against COVID-19 is still under discussion. Most of the people infected with COVID-19 developed the immunity against the virus, but still it is unknown how long this immunity will last as there are some reported cases of getting the infection for a second time. In this situation, the only solution to gain herd immunity is mass vaccination. Hence, to achieve this milestone, vaccination is required for a huge portion of the global population. In this regard, WHO helped governments to launch and refine their vaccination programmes along with UNICEF, Gavi, and partners. According to reference [19], 43.5% of the world population has been given at least one shot of COVID-19 vaccine till the 21st of September 2021. Modified mathematical models with an additional compartment of the vaccinated class are discussed in recent literature to foresee the efficacy, effect, and dynamics of the vaccination. Meng et al. [20] established an SEIRV model to differentiate between the effects of mandatory and voluntary vaccination methods on heterogeneous networks. Ramos et al. [21] investigated the impact of variants and vaccines on the spread of COVID-19 applying the model to Italy as a test case. Tetteh et al. [22] analyzed the mass and ring vaccination strategies to control the transmission of SARS-CoV-2 by varying vaccine efficacy rate. Harizi et al. [23] used the SIRV model to check the future dynamics of the pandemic by varying vaccination rates in Canada. Karabay et al. [24] presented an SEIR simulator to model effective immunization and sterilizing outcomes in order to assess the vaccination effect. Tuteja [25] numerically solved the modified SEIR model including the vaccinated class to investigate the viral dynamics of the epidemic. Ghostine [26] proposed an SEIR model with an extension of the vaccination compartment to investigate the effect of vaccination in Saudi Arabia. Wintachaia and Prathomb [27] presented a detailed stability analysis of the COVID-19 epidemic model to investigate the vaccine efficiency. Alsakaji et al. [28] examined the dynamics of COVID-19 in the UAE using an exptended SEIR epidemic model with vaccination, temporal delays, and random noise. For transmission of the infectious disease, Kou et al. [29] established a multiscale agent-based model to investigate the impact of nonpharmaceutical interventions and vaccination. Acuña-Zegarra. [30] formulated a COVID-19 optimal control problem to study the efficacy and responses of the vaccine induced immunity. Albani et al. [31] estimated the parameters of SEIR like epidemic model using daily reports of COVID-19 of New York and Chicago to analyze the impact caused by vaccination delay.

Motivated from the above studies, we have presented in this work a modified SEIRV epidemic model with an additional compartment of vaccinated individuals to examine the effectiveness and efficacy of the vaccine by varying its rate considering the ABC derivative. We have used ABC fractional order derivative because a number of experimental facts exhibit that the natural dynamics go along with fractional calculus as the fractional order derivatives are nonlocal in behaviour and possess the hereditary properties and memory effects. Also, we have constructed a new, more efficient, and effective numerical scheme for fractional differential equation with Atangana-Baleanu derivative in Caputo sense using Simpson's 1/3 method for the approximation of the integral function involved in the problem.

The paper comprises different sections. In section 2, some basic definitions are recalled which have been used in the work later. SEIRV epidemic model for COVID-19 is presented in section 3. The equilibrium points and the basic reproduction number using next generation method are calculated, and the stability of these points is also analyzed in this section. For the solution of the proposed model modified Simpson's 1/3 method alongside complete stability and convergence analysis, it is presented in section 4. Conversion of mathematical model to fractional order model is given in section 5. Numerical simulations are also presented in the same section.

2. Preliminaries

This part will review some of the fundamental definitions of fractional calculus, which will be useful in the subsequent work.

Definition 1 .

For x ∈ ℍ1(0, T) and 0 < ϑ < 1, the ϑth–order ABC derivative and integral of x(t) are defined [32] as where J(ϑ) represents the normalization function given by with J(ϑ)| = 1 and E is the Mittag-Leffler function defined as E(z) = ∑∞(z/Γ(nϑ + 1)).

The lemma below presents the Gronwall inequality [33], which will be used for the proof of our main result.

Lemma 1 .

The Gronwall inequality for fractional order differential equation states

Let C1 > 0 independent of δt, C2 ≥ 0 and {u} satisfy the inequality with 0 < ϑ ≤ 1. Then where the parameter E is the Mittage-Leffler function. Particularly when ϑ = 1, inequality (5) becomes

In this study, we used an SEIR model with vaccination effectiveness to forecast the COVID-19 condition when a vaccine is released.

3. Mathematical Model

Mathematical models are of great importance in predicting the behaviour of the viral epidemic diseases. There are a number of models which can be used to describe and investigate the spread of the disease. Recently COVID-19 emerged as a viral epidemic, and we are not aware enough about the behaviour and spread of this virus. To investigate the COVID-19 situation under the effect of vaccination, we consider the following SEIRV epidemic disease model:

To define the model equations, the total population, represented by N, is separated into susceptible, exposed, infected, recovered, and vaccinated classes denoted by S(t), E(t), I(t), R(t), and V(t). The flow chart of the above system is shown in Figure 1.

Figure 1

SEIRV Model.

The above system of equation (7) is subject to the initial conditions

and the following are the parameters:

μ: per-capita natural death rate,

γ: the rate at which infectious people recover,

ρ: average fatality rate for virus-infected people,

β: disease transmission probability per contact times the amount of contacts made per unit of time (dimensionless),

Λ: the birth rate per-capita,

ε: progression rate from E to I,

1/γ: contagious period,

1/ε: incubation period,

δ: vaccination rate,

σ: vaccine ineffectiveness, 0 ≤ σ ≤ 1,

1 − σ: vaccine effectiveness,

having T units per unit of time.

3.1. Stability Analysis of SEIRV COVID-19 Model

In this section, stability analysis will be carried out to find the disease-free equilibrium and endemic equilibrium points. In order to get these points, each of the equations given in system of equation (7) will be equated to zero.

3.2. Disease-Free Equilibrium and the Basic Reproduction Number (R0)

Because there is no disease propagation in disease-free equilibrium, so I = 0, and hence we get so the disease-free equilibrium point is 𝔼0(S0, E0, I0, R0, V0) = (Λ/(μ + δ), 0, 0, 0, δΛ/μ(μ + δ)).

The average number of secondary infections induced by an infectious individual is given by the basic reproduction number R0. It estimates the growth of the virus outbreak. The basic reproduction number R0 is determined using the next generation method [34]. Using and from system of equation (7), the Jacobian matrices of the disease-free equilibrium are as follows:

Here FU−1 is the next generation matrix and we can obtain the basic reproduction number as R0 = ρ(FU−1), so

Substituting the values of S0 and V0, the following value of the reproduction number is obtained for the system of equation (7):

Theorem 1 .

If R0 < 1 then the COVID-19 model of equation (7) disease-free equilibrium point is locally asymptotically stable.

Proof

The Jacobian of system of equation (7) is given as

Jacobian around disease-free equilibrium point 𝔼0((Λ/μ + δ), 0, 0, 0, (δΛ/μ(μ + δ))) is

We can observe from the preceding matrix that three of its eigenvalues λ1,2 = −μ, λ3 = −μ − δ are negative and we are left with

Now from above matrix, the last two eigenvalues are where ϖ1 = (μ + ε) and ϖ2 = (γ + μ + ρ), and for R0 < 1, we have λ4 < 0 and λ5 < 0.

Hence for R0 < 1 the disease-free equilibrium is locally asymptotically stable.

3.3. Existence and Uniqueness of Endemic Equilibrium

In order to investigate the existence and uniqueness of the endemic equilibrium 𝔼∗ for R0 > 1, we will set all the derivatives of model of equation (7) equals to zero. Let 𝔼∗ = (S∗, E∗, I∗, R∗, V∗) be the endemic equilibrium of model of equation (7). Then by solving the following equations we can obtain the endemic equilibrium.

Now, solving first, third, and fourth equations of equation (17), we get

Also equation (17) gives

The following quadratic equation for I∗ is obtained: where

It is evident that A is always positive while C is negative for R0 > 1. Hence, there exists a positive and unique value of I∗, which results in a unique endemic equilibrium for R0 > 1.

4. Numerical Approximation

Consider the nonlinear fractional order differential equation where 0 < ϑ < 1. The vector function f in the differential equation is real valued and continuous and satisfies the Lipschitz condition.

To derive the numerical scheme for the above mentioned nonlinear fractional order differential equation, the given interval [0, T] will be subdivided into equal intervals with time step δτ = T/N, where δτ − δτ = δτ for n = 0, 1, 2, ⋯, N − 1. For convenience, the approximate solution of x(τ) will be denoted by x with δτ = nδτ, n = 0, 1, 2, ⋯, N.

Applying the integral operator of equation (2) on equation (22), we get where the integral I in the above equation is the ABC integral. Using equation (2) to discretize the integral I, we get

f(x(s)) is approximated on each subinterval [τ, τ], by the following piecewise quadratic interpolation polynomial. where τ = (τ + τ)/2 is mid point of the interval [τ, τ].

Also, the value of f(x(τ), τ) will be approximated using the interpolation

Putting equation (27) in equation (26), we get

Using above approximation, equation (25) takes the form

Evaluating all the integrals given in equation (30), we obtain where

Thus we can obtain the following modified fractional Simpson's 1/3 scheme for equation (22) in ABC sense:

4.1. Convergence and Stability Analysis

The following lemmas are very important to examine the proposed scheme's convergence and stability.

Lemma 2 .

Assume that g(t) ∈ C3([0, T]), then there exist a constant C such that

Using the Lagrange interpolation error formula, we have

From above lemma it is clear that the truncation error of the scheme of equation (33) is O((δτ)2).

Lemma 3 .

For 0 < ϑ < 1, the weights b, i = 0, 1, 2, ⋯, n − 2 given in equation (33) have the following order of magnitude:

According to the definition of b(i = 1, 2, ⋯, n − 2), it deduces that implies

Note that 0 < ϑ < 1, so the coefficients of power of 1/(n − i) are less than 1. Thus, the series given above is convergent, and hence the proof is completed for b(i = 1, 2, ⋯, n − 2). Similar strategy can be adopted for b, so we neglect the proof.

Theorem 2 .

The numerical scheme (33) is conditionally convergent. That is, for 0 < ϑ ≤ 1 and a sufficiently small δτ, there exist a constant C1 such that

To prove the theorem, considering the difference between actual solution and approximate solution, we get

Using Lipschitz condition of equation (23) after applying triangle inequality and using Lemma 4.1, we get

Now using Lemma 4.2 under the condition (L/J(ϑ))(1 − ϑ + (ϑδτ/Γ(ϑ + 3))) < 1, we get where g(δτ, ϑ) = (J(ϑ)Γ(ϑ + 3))/(J(ϑ)Γ(ϑ + 3) − L((Γ(ϑ + 3)(1 − ϑ) + ϑδτ)) For any 0 < ϑ < 1 and a sufficiently small δτ, there exist a constant C such that

Hence by applying Gronwall inequality given in Lemma 2.1 and using the property (43), we conclude

Theorem 3 .

For 0 < ϑ ≤ 1 and a sufficiently small δτ, there exist a constant C1 such that the numerical scheme (33) is stable.

To prove the theorem, let the perturbation of x and x0 be denoted by and , respectively, for n = 0, 1, ⋯, N − 1. Adding this perturbation to the numerical scheme (33) we get

Subtracting equation (33) from equation (45), we get

Using the Lipschitz condition followed by triangle inequality, we get

Now if (L/J(ϑ))(1 − ϑ + (ϑδτ/Γ(ϑ + 3))) < 1, then using Lemma 4.2 we get where g(δτ, ϑ) = J(ϑ)Γ(ϑ + 3)/(J(ϑ)Γ(ϑ + 3) − L(Γ(ϑ + 3)(1 − ϑ) + ϑδτ).

Using property (43), we have

Hence by applying the Gronwall inequality given in Lemma 2.1, we conclude

5. Fractional Model

Employing Atangana-Baleanu fractional order derivative in Caputo sense of order ϑ on system of equation (7), we have the following fractional model, along with initial conditions

Adopting the previously defined procedure, model of equation (51) takes the form

Furthermore applying the numerical scheme(33), system (53) becomes where a are already defined.

5.1. The Positivity and Boundness of the Solution

In this subsection we will discuss the positive nature of the epidemiological model of equation (51). We will recall the following lemma for the proof of the main theorem on the positivity of the solution of the fractional model (51).

Lemma 4 .

Generalized Mean Value Theorem: let ϕ(τ) ∈ C[a, b] and Atangana-Baleanu fractional order derivative, ABCDϕ(τ) ∈ C(a, b]for0 < ϑ ≤ 1, then with 0 ≤ u ≤ τ, ∀τ ∈ (a, b].

Remark 1 .

Let ϕ(τ) ∈ C[0, b] and Atangana-Baleanu fractional order derivative, ABCD ∈ (0, b) for 0 < ϑ ≤ 1. Lemma 5.1 shows that if  ABCDϕ(τ) > 0∀ τ ∈ (0, b], then the function ϕ(τ) is nondecreasing and if  ABCDϕ(τ) ≤ 0∀ τ ∈ (0, b] then ϕ(τ) is nonincreasing for all τ ∈ [0, b].

Theorem 4 .

The solution to the nonnegative initial conditions of equation (52) of the proposed model of equation (51) is unique and limited in R+5.

The existence and uniqueness of the model solution of equations (51) and (52) can be achieved in the time interval (0, ∞) by [35]. The nonnegative region R+5 would seem to be a positive invariant region that must be represented. We observe from the model of equation (51)

From Remark 1 and system of equation (56), the solution (S(0), E(0), I(0), R(0), V(0)) ∈ R+5. Also, in each line bounded the nonnegative octant, the vector field points will remain in R+5. Therefore, the fractional model of equation (51) of a solution ((S(τ), E(τ), I(τ), R(τ), V(τ)) is nonnegative if the initial condition is set positively to invariant.

6. Numerical Results and Discussion

Considering the case of Italy, we applied the proposed model to analyse the impact of vaccination. The vaccination campaign for COVID-19 started on 27th December 2020 in Italy managed by the Ministry of Health. In the early months of the vaccination campaign, the Italian government targeted health doctors and administrative staff as well as guests and nursing home personnel followed the vaccination of elderly people and public service personnel in second phase of the campaign. Italian government bought four out of seven, WHO approved vaccines to facilitate its people against the virus. So far, 73.3% of the total population has received at least one dose and 65.5% of the total population is fully vaccinated. As of 17 September 2021, top ten regions of Italy for the total number of administered COVID-19 vaccine doses can be seen in Figure 2. Moreover the number of people who received at least one dose and fully vaccinated according to the data [19] are presented in Figure 3.

Figure 2

Number of COVID-19 vaccine doses administered in Italy (top ten regions).

Figure 3

Number of Vaccinated People.

Further, graphical presentation of the numerical simulation for the fractional model (51) are given in this section. Table 1 shows the estimated model parameters and initial values used to evaluate the impact of vaccination on the COVID-19 epidemic.

Table 1

Parameter and initial values for the numerical simulations.

Parameter	Value	Parameter	Value	Initial time model variables	Value
N	60000000	ε	1/4.02	S(0)	59988539
μ	0.0000373	Λ	μ × N	E(0)	11460
γ	1/5.25	ϑ	0.00144	I(0)	1
β	0.75	σ	0.25	R(0)	0

Figure 4 shows the impact by varying the vaccination rates on different classes of the proposed model with δ = 0 representing the model output without vaccination. Although there is no remarkable change in the graph of susceptible class, we can observe a significant decrease in infected and exposed classes by increasing the vaccination rate. We can see the peak of the infected as well as that of exposed individuals gradually decreases with the increase in effective rate of vaccination. From Figure 4, we can observe that for exposed class there is a difference of almost 7 × 106 individuals in peak values for no vaccination to a vaccination rate δ = 0.012. Similarly, the peak values for no vaccination compared to a vaccination rate of δ = 0.012 for the infected class differ by almost 7 × 106 individuals. So, interpreting graphically the dynamics of the proposed model to observe the effect of vaccination by altering the vaccination rate suggest that, mass vaccination along with the implementation of the control measures until the achievement of complete herd immunity is needed to control the spread of the disease. Hence, majority of susceptible population should be vaccinated to prevent the new waves of COVID-19. The vaccination of the susceptible class will control the spread of the disease resulting in reduction of the mortality rate, which can clearly be seen by recovered class presented in Figure 4. Simulation results show that immunization is an essential strategy for lowering the incidence and severity of COVID-19 infections in the workplace, in communities, and across the globe. Therefore, it is advocated for policymakers to increase public trust and confidence in the efficiency and security of vaccinations as well as the expertise and dependability of the institutions that provide them. Figures 5 and 6 show a comparison of integer order model, fractional model for fractional order ϑ = 0.975,0.95,0.925, 0.9, and real data of total number of deaths D = ϑI owing to the COVID-19 outbreak in Italy. Figure 5 presents the comparison for first 80 days starting from 23 February 2020 to 13 May 2020 of the pandemic before vaccination. We can see that after 40 days, real data is better explained by fractional order model in comparison to integer order. Figure 6 is assigned to display the comparison from 1 July 2021 to 19 September 2021 along with administered vaccination effect considering the initial values of the compartments as of 1 July 2021 . The relevant COVID-19 data can be found at https://github.com/pcmdpc/COVID-19. Also we can see from our proposed model, fractional order models are better at fitting the real data than classical models because fractional order derivative and integral operators, due to their nonlocal nature, depend not only on their current state but also on all of their previous states, making them more effective than other classical deterministic operators in predicting the biological model's future state.

Figure 4

Impact of vaccination rate on all the four susceptible, exposed, infected, and recovered classes.

Figure 5

Comparison of the number of people who died for various fractional orders vs. real data (before vaccination).

Figure 6

Comparison of the number of people who died for various fractional orders vs. real data (after vaccination).

7. Concluding Remarks

Within the Atangana-Baleanu Caputo fractional operator, we studied alternative conditions for qualitative characterization and dynamics of the suggested model. Calculating disease-free and endemic equilibrium points, as well as the reproduction number, are used to examine the prerequisites for local stability. A newly developed numerical scheme is used to obtain the solution of the concerned fractional order model. The proposed scheme's stability and convergence are briefly examined. To depict the behaviour of the COVID-19 virus, a graphical description of the dynamics of the relevant system is developed within different fractional orders under the effect of vaccination suggesting the need of mass vaccination to gain the herd immunity in order to prevent the further spread of the disease. Due to the limited access of the available vaccine, other preventive measures should also be continued.