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

Authors: Mubashara Wali; Sadia Arshad; Jianfei Huang
Journal: Computational and Mathematical Methods in Medicine
 (2022)

---

## Abstract

This work is aimed at presenting a new numerical scheme forCOVID-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.


---

## Body

## 1. Introduction

Emerging in December2019, 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 section2, 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.
Forx∈ℍ10,T and 0<ϑ<1, the ϑth–order ABC derivative and integral of xt are defined [32] as
(1)ABCDϑxt≔Jϑ1−ϑ∫0tEϑ−ϑ1−ϑt−ξϑẋξdξ,(2)ABCIϑxt≔1−ϑJϑxt+ϑJϑΓϑ∫0tt−ξϑ−1xξdξ,where Jϑ represents the normalization function given by
(3)Jϑ=1−ϑ+ϑΓϑ,with Jϑϑ=0,1=1 and Eϑ is the Mittag-Leffler function defined as Eϑz=∑n=0∞zn/Γnϑ+1.The lemma below presents the Gronwall inequality [33], which will be used for the proof of our main result.Lemma 2.1.
The Gronwall inequality for fractional order differential equation states
LetC1>0 independent of δt, C2≥0 and {un} satisfy the inequality
(4)un≤δtϑC1∑l=0n−1n−lϑ−1uj+C2,l=0,1,⋯,n−1,nδt≤T,with 0<ϑ≤1. Then
(5)un≤C2EϑC1ΓϑTϑ,nδt≤T,where the parameter Eϑ is the Mittage-Leffler function. Particularly when ϑ=1, inequality (5) becomes
(6)un≤C2eC1T,nδt≤T.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. RecentlyCOVID-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:
(7)Ṡ=Λ−μS−δS−βSIN,Ė=βSIN−μ+εE+σβVIN,İ=εE−γ+μ+ρI,Ṙ=γI−μR,V̇=δS−μV−σβVIN.To define the model equations, the total population, represented byN, is separated into susceptible, exposed, infected, recovered, and vaccinated classes denoted by St, Et, It, Rt, and Vt. 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
(8)S0=S0,E0=E0,I0=I0,R0=R0,V0=V0,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 SEIRVCOVID-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 NumberR0

Because there is no disease propagation in disease-free equilibrium, soI=0, and hence we get
(9)S0=Λμ+δ,E0=0,R0=0,V0=δΛμμ+δ,so the disease-free equilibrium point is E0S0,E0,I0,R0,V0=Λ/μ+δ,0,0,0,δΛ/μμ+δ.The average number of secondary infections induced by an infectious individual is given by the basic reproduction numberR0. It estimates the growth of the virus outbreak. The basic reproduction number R0 is determined using the next generation method [34]. Using Ė and İ from system of equation (7), the Jacobian matrices of the disease-free equilibrium are as follows:
(10)F=0βS0+σβV000,U=μ+ε0−εγ+μ+ρ,U−1=1μ+εγ+μ+ργ+μ+ρ0εμ+ε,FU−1=1μ+εγ+μ+ρ0βS0+σβV000γ+μ+ρ0εμ+ε,FU−1=1μ+εγ+μ+ρβεS0+σβεV0βS0+σβV0μ+ε00.HereFU−1 is the next generation matrix and we can obtain the basic reproduction number as R0=ρFU−1, so
(11)R0=βεS0+βεσV0μ+εγ+μ+ρ.Substituting the values ofS0 and V0, the following value of the reproduction number is obtained for the system of equation (7):
(12)R0=βεΛμ+δσμμ+δμ+εγ+μ+ρ.Theorem 3.1.
IfR0<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
(13)J=−μ−δ−βIN0−βSN00βIN−μ+εβSN+σβVN0σβIN0ε−γ+μ+ρ0000γ−μ0δ0−σβVN0−μ.
Jacobian around disease-free equilibrium pointE0Λ/μ+δ,0,0,0,δΛ/μμ+δ is
(14)J0=−μ−δ0−βΛμ+δ000−μ+εβΛμ+δ+σβδΛμμ+δ000ε−γ+μ+ρ0000γ−μ0δ0−σβδΛμμ+δ0−μ.
We can observe from the preceding matrix that three of its eigenvaluesλ1,2=−μ, λ3=−μ−δ are negative and we are left with
(15)J0=−μ+εβΛμ+δ+σβδΛμμ+δε−γ+μ+ρ.
Now from above matrix, the last two eigenvalues are(16)λ4=−12ϖ1+ϖ2+ϖ1−ϖ22+4ϖ1ϖ2R0,λ5=−12ϖ1+ϖ2−ϖ1−ϖ22+4ϖ1ϖ2R0,where ϖ1=μ+ε and ϖ2=γ+μ+ρ, and for R0<1, we have λ4<0 and λ5<0.
Hence forR0<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 equilibriumE∗ for R0>1, we will set all the derivatives of model of equation (7) equals to zero. Let E∗=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.
(17)Λ−μS∗−δS∗−βS∗I∗N=0,βS∗I∗N−μ+εE∗+σβV∗I∗N=0,εE∗−γ+μ+ρI∗=0,γI∗−μR∗=0,δS∗−μV∗−σβV∗I∗N=0.Now, solving first, third, and fourth equations of equation (17), we get
(18)S∗=Λμ+δ+βI∗/N,E∗=γ+μ+ρεI∗,R∗=γμI∗.Also equation (17) gives
(19)V∗=βS∗I∗μN−μ+εμE∗+δμS∗.The following quadratic equation forI∗ is obtained:
(20)AI∗2+BI∗+C=0,where
(21)A=σβ2ρμ+εγ+μ+ρ,B=μNμ+εγ+μ+ρβ+σβμ+εγ+μ+ρN−NσΛβ2ρε,C=μN2μ+δμ+εγ+μ+ρ1−R0.It is evident thatA 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.


## 3.1. Stability Analysis of SEIRVCOVID-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 NumberR0

Because there is no disease propagation in disease-free equilibrium, soI=0, and hence we get
(9)S0=Λμ+δ,E0=0,R0=0,V0=δΛμμ+δ,so the disease-free equilibrium point is E0S0,E0,I0,R0,V0=Λ/μ+δ,0,0,0,δΛ/μμ+δ.The average number of secondary infections induced by an infectious individual is given by the basic reproduction numberR0. It estimates the growth of the virus outbreak. The basic reproduction number R0 is determined using the next generation method [34]. Using Ė and İ from system of equation (7), the Jacobian matrices of the disease-free equilibrium are as follows:
(10)F=0βS0+σβV000,U=μ+ε0−εγ+μ+ρ,U−1=1μ+εγ+μ+ργ+μ+ρ0εμ+ε,FU−1=1μ+εγ+μ+ρ0βS0+σβV000γ+μ+ρ0εμ+ε,FU−1=1μ+εγ+μ+ρβεS0+σβεV0βS0+σβV0μ+ε00.HereFU−1 is the next generation matrix and we can obtain the basic reproduction number as R0=ρFU−1, so
(11)R0=βεS0+βεσV0μ+εγ+μ+ρ.Substituting the values ofS0 and V0, the following value of the reproduction number is obtained for the system of equation (7):
(12)R0=βεΛμ+δσμμ+δμ+εγ+μ+ρ.Theorem 3.1.
IfR0<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
(13)J=−μ−δ−βIN0−βSN00βIN−μ+εβSN+σβVN0σβIN0ε−γ+μ+ρ0000γ−μ0δ0−σβVN0−μ.
Jacobian around disease-free equilibrium pointE0Λ/μ+δ,0,0,0,δΛ/μμ+δ is
(14)J0=−μ−δ0−βΛμ+δ000−μ+εβΛμ+δ+σβδΛμμ+δ000ε−γ+μ+ρ0000γ−μ0δ0−σβδΛμμ+δ0−μ.
We can observe from the preceding matrix that three of its eigenvaluesλ1,2=−μ, λ3=−μ−δ are negative and we are left with
(15)J0=−μ+εβΛμ+δ+σβδΛμμ+δε−γ+μ+ρ.
Now from above matrix, the last two eigenvalues are(16)λ4=−12ϖ1+ϖ2+ϖ1−ϖ22+4ϖ1ϖ2R0,λ5=−12ϖ1+ϖ2−ϖ1−ϖ22+4ϖ1ϖ2R0,where ϖ1=μ+ε and ϖ2=γ+μ+ρ, and for R0<1, we have λ4<0 and λ5<0.
Hence forR0<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 equilibriumE∗ for R0>1, we will set all the derivatives of model of equation (7) equals to zero. Let E∗=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.
(17)Λ−μS∗−δS∗−βS∗I∗N=0,βS∗I∗N−μ+εE∗+σβV∗I∗N=0,εE∗−γ+μ+ρI∗=0,γI∗−μR∗=0,δS∗−μV∗−σβV∗I∗N=0.Now, solving first, third, and fourth equations of equation (17), we get
(18)S∗=Λμ+δ+βI∗/N,E∗=γ+μ+ρεI∗,R∗=γμI∗.Also equation (17) gives
(19)V∗=βS∗I∗μN−μ+εμE∗+δμS∗.The following quadratic equation forI∗ is obtained:
(20)AI∗2+BI∗+C=0,where
(21)A=σβ2ρμ+εγ+μ+ρ,B=μNμ+εγ+μ+ρβ+σβμ+εγ+μ+ρN−NσΛβ2ρε,C=μN2μ+δμ+εγ+μ+ρ1−R0.It is evident thatA 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(22)ABCDϑxτ=fxτ,τ,0<τ<T<∞,x=0=x0,where 0<ϑ<1. The vector function f in the differential equation is real valued and continuous and satisfies the Lipschitz condition.
(23)fx1τ−fx2τ≤Lx1τ−x2τ,L>0.To derive the numerical scheme for the above mentioned nonlinear fractional order differential equation, the given interval0,T will be subdivided into equal intervals with time step δτ=T/N, where δτn+1−δτn=δτ for n=0,1,2,⋯,N−1. For convenience, the approximate solution of xτn will be denoted by xn with δτn=nδτ, n=0,1,2,⋯,N.Applying the integral operator of equation (2) on equation (22), we get
(24)xτ=x0+Iϑfxτ,τ,0<τ<T<∞,where the integral Iϑ in the above equation is the ABC integral. Using equation (2) to discretize the integral Iϑ, we get
(25)ABCIϑfxτn+1,τn+1=1−ϑJϑfxτn+1,τn+1+ϑJϑΓϑ∫0τn+1τn+1−sϑ−1fxs,sds,n=0,1,2,⋯,N−1.fxs is approximated on each subinterval τi,τi+1, by the following piecewise quadratic interpolation polynomial.
(26)fxs,sτi,τi+1≈τi+1−sτi+1/2−sτi−τi+1τi−τi+1/2fxτi,τi+τi+1−sτi−sτi+1/2−τi+1τi+1/2−τifxτi+1/2,τi+1/2+τi+1/2−sτi−sτi+1−τiτi+1−τi+1/2fxτi+1,τi+1,where τi+1/2=τi+τi+1/2 is mid point of the interval τi,τi+1.Also, the value offxτi+1/2,τi+1/2 will be approximated using the interpolation
(27)fxτi+1/2,τi+1/2=12fxτi,τi+12fxτi+1,τi+1.Putting equation (27) in equation (26), we get
(28)fxs,sτi,τi+1≈τi+1−sτi+1/2−sτi−τi+1τi−τi+1/2fxτi,τi+τi+1−sτi−sτi+1/2−τi+1τi+1/2−τi12fxτi,τi+12fxτi+1,τi+1+τi+1/2−sτi−sτi+1−τiτi+1−τi+1/2fxτi+1,τi+1=2δτ2τi+1−sτi+1/2−s−τi+1−sτi−sfxτi,τi+τi+1/2−sτi−s−τi+1−sτi−sfxτi+1,τi+1.Using above approximation, equation (25) takes the form
(29)ABCIϑfxτn+1,τn+1=1−ϑJϑfxτn+1,τn+1+ϑJϑ2δτ2Γϑ∑i=0n∫τiτi+1τn+1−sϑ−1τi+1−sτi+1/2−s−τi+1−sτi−sfxτi,τi+τi+1/2−sτi−s−τi+1−sτi−sfxτi+1,τi+1(30)=1−ϑJϑfxτn+1,τn+1+ϑJϑ2δτ2Γϑ∫τ0τ1τn+1−sϑ−1τ1−sτ1/2−s−τ1−sτ0−sfxτ0,τ0ds+∫τ0τ1τn+1−sϑ−1τ0−sτ1/2−s−τ1−sτ0−sfxτ1,τ1ds+∑k=1n−1∫τiτi+1τn+1−sϑ−1τi+1−sτi+1/2−s−τi+1−sτi−sfxτi,τi+τi−sτi+1/2−s−τi+1−sτi−sfxτi+1,τi+1ds+∫tnτn+1τn+1−sϑ−1τn+1−sτn+1/2−s−τn+1−stn−sfxτn,τnds+∫tnτn+1τn+1−sϑ−1τn−sτn+1/2−s−τn+1−sτn−sfxτn+1,τn+1ds.Evaluating all the integrals given in equation (30), we obtain
(31)ABCIϑfxτn+1,τn+1=1−ϑJϑfxτn+1,τn+1+ϑδτϑJϑ∑i=0n+1bn+1,ifxτi,τi,n=0,1,2,⋯,N−1,where
(32)bn+1,i=1Γϑ+3ϑ+1ϑ+2n+1ϑ+ϑ+2nϑ+1−ϑ+2n+1ϑ+1i=0,−2ϑ+2n+1−iϑ+1+ϑ+2n+2−iϑ+1+ϑ+2n−iϑ+11≤i≤n,2+ϑi=n+1.Thus we can obtain the following modified fractional Simpson’s 1/3 scheme for equation (22) in ABC sense:
(33)xn+1=x0+1−ϑJϑfxτn+1,τn+1+ϑδτϑJϑ∑i=0n+1bn+1,ifxτi,τi,n=0,1,2,⋯,N−1.


### 4.1. Convergence and Stability Analysis

The following lemmas are very important to examine the proposed scheme’s convergence and stability.Lemma 4.1.
Assume thatgt∈C30,T, then there exist a constant C such that
(34)1Γϑ∫τoτn+1τn+1−sϑ−1gsds−δτϑ∑i=0n+1bigτi≤TϑCδτ22!Γϑ+1.Proof.
Using the Lagrange interpolation error formula, we have(35)1Γϑ∫τoτn+1τn+1−sϑ−1gsds−δτϑ∑i=0n+1bigτi≤12!Γϑ∑i=0n∫τiτi+1τn+1−sϑ−1g′′ξis−τis−τi+1ds≤Cδτ22!Γϑ∫τ0τn+1τn+1−sϑ−1ds≤TϑCδτ22!Γϑ+1.
From above lemma it is clear that the truncation error of the scheme of equation (33) is Oδτ2.Lemma 4.2.
For0<ϑ<1, the weights bi, i=0,1,2,⋯,n−2 given in equation (33) have the following order of magnitude:
(36)bi=On−iϑ−1.Proof.
According to the definition ofbii=1,2,⋯,n−2, it deduces that
(37)bi=n−iϑ−1Γϑ+3−2ϑ+2n−i21+1n−iϑ+1+ϑ+2n−i21+2n−iϑ+1+ϑ+2n−i2,bi=ϑ+2n−iϑ−1Γϑ+3−2n−i2+ϑ+1n−i+ϑ+1ϑ2!+ϑ+1ϑϑ−13!1n−i+ϑ+1ϑϑ−1ϑ−24!1n−i2+⋯+n−i2+2ϑ+1n−i2+ϑϑ+1+ϑ+1ϑϑ−13!2n−i+ϑ+1ϑϑ−1ϑ−24!2n−i2+⋯+n−i2,implies
(38)bi=ϑ+2n−iϑ−1Γϑ+32ϑ+1ϑϑ−1ϑ−24!1n−i2+6ϑ+1ϑϑ−1ϑ−2ϑ−35!1n−i3+14ϑ+1ϑϑ−1ϑ−2ϑ−3ϑ−46!1n−i4+⋯.
Note that0<ϑ<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 bii=1,2,⋯,n−2. Similar strategy can be adopted for bn, so we neglect the proof.Theorem 4.1.
The numerical scheme (33) is conditionally convergent. That is, for 0<ϑ≤1 and a sufficiently small δτ, there exist a constant C1 such that
(39)xτn+1−xn+1≤C1δτ2.Proof.
To prove the theorem, considering the difference between actual solution and approximate solution, we get(40)xτn+1−xn+1=1−ϑJϑfxτn+1,τn+1−fxn+1,τn+1+ϑJϑΓϑ∫0τn+1τn+1−sϑ−1fxs,sds−ϑδτϑJϑ∑i=0n+1bn+1,ifxi,τi=1−ϑJϑfxτn+1,τn+1−fxn+1,τn+1+ϑδτϑJϑ∑i=0n+1bn+1,ifxτi,τi−fxi,τi+ϑJϑ1Γϑ∫0τn+1τn+1−sϑ−1fxs,sds−δτϑ∑i=0n+1bn+1,ifxτi,τi,xτn+1−xn+1=1−ϑJϑfxτn+1,τn+1−fxn+1,τn+1+ϑδτϑJϑ∑i=0n+1fxτi,τi−fxi,τi+ϑJϑ1Γϑ∫0τn+1τn+1−sϑ−1fxs,sds−δτϑ∑i=0n+1bn+1,ifxτi,τi.
Using Lipschitz condition of equation (23) after applying triangle inequality and using Lemma 4.1, we get
(41)xτn+1−xn+1≤1−ϑJϑfxτn+1,τn+1−fxn+1,τn+1+ϑδτϑJϑ∑i=0n+1bn+1,ifxτi,τi−fxi,τi+ϑJϑ1Γϑ∫0τn+1τn+1−sϑ−1fxs,sds−δτϑ∑i=0n+1bn+1,ifxτi,τi≤1−ϑLJϑxτn+1−xn+1+ϑTϑJϑCδτ22Γϑ+1+ϑδτϑLJϑ∑i=0n+1bn+1,ixτi−xi≤LJϑ1−ϑ+ϑδτϑΓϑ+3xτn+1−xn+1+ϑTϑJϑCδτ22Γϑ+1+ϑδτϑLJϑ∑i=0nbn+1,ixτi−xi.
Now using Lemma4.2 under the condition L/Jϑ1−ϑ+ϑδτϑ/Γϑ+3<1, we get
(42)xτn+1−xn+1≤gδτ,ϑϑTϑJϑCδτ22Γϑ+1+gδτ,ϑϑδτϑLJϑ∑i=0nn−iϑ−1xτi−xi,where gδτ,ϑ=JϑΓϑ+3/JϑΓϑ+3−LΓϑ+31−ϑ+ϑδτϑ For any 0<ϑ<1 and a sufficiently small δτ, there exist a constant Cg such that
(43)1<gδτ,ϑ<Cg.
Hence by applying Gronwall inequality given in Lemma2.1 and using the property (43), we conclude
(44)xτn+1−xn+1≤C1δτ2.Theorem 4.2.
For0<ϑ≤1 and a sufficiently small δτ, there exist a constant C1 such that the numerical scheme (33) is stable.Proof.
To prove the theorem, let the perturbation ofxn and x0 be denoted by x~n and x~0, respectively, for n=0,1,⋯,N−1. Adding this perturbation to the numerical scheme (33) we get
(45)xn+1+x~n+1=x0+x~0+1−ϑJϑfxn+1+x~n+1+ϑδτϑJϑ∑i=0n+1bn+1,ifxi+x~i.
Subtracting equation (33) from equation (45), we get
(46)x~n+1=x~0+1−ϑJϑfxn+1+x~n+1−fxn+1+ϑδτϑJϑ∑i=0n+1bn+1,ifxi+x~i−fxi.
Using the Lipschitz condition followed by triangle inequality, we get(47)x~n+1=x~0+1−ϑJϑfxn+1+x~n+1−fxn+1+ϑδτϑJϑ∑i=0n+1bn+1,ifxi+x~i−fxi,≤x~0+1−ϑJϑfxn+1+x~n+1−fxn+1+ϑδτϑJϑ∑i=0n+1bn+1,ifxi+x~i−fxi,≤x~0+1−ϑLJϑx~n+1+ϑδτϑLJϑ∑i=0n+1bn+1,ix~i,≤x~0+LJϑ1−ϑ+ϑδτϑΓϑ+3x~n+1+ϑδτϑLJϑ∑i=0nbn+1,ix~i.
Now ifL/Jϑ1−ϑ+ϑδτϑ/Γϑ+3<1, then using Lemma 4.2 we get
(48)x~n+1≤gδτ,ϑx~0+gδτ,ϑϑδτϑLJϑ∑i=0nn−iϑ−1x~i,where gδτ,ϑ=JϑΓϑ+3/JϑΓϑ+3−LΓϑ+31−ϑ+ϑδτϑ.
Using property (43), we have
(49)x~n+1≤Cgx~0+CgϑδτϑLJϑ∑i=0nn−iϑ−1x~i.
Hence by applying the Gronwall inequality given in Lemma2.1, we conclude
(50)x~n+1≤C1x~0.


## 4.1. Convergence and Stability Analysis

The following lemmas are very important to examine the proposed scheme’s convergence and stability.Lemma 4.1.
Assume thatgt∈C30,T, then there exist a constant C such that
(34)1Γϑ∫τoτn+1τn+1−sϑ−1gsds−δτϑ∑i=0n+1bigτi≤TϑCδτ22!Γϑ+1.Proof.
Using the Lagrange interpolation error formula, we have(35)1Γϑ∫τoτn+1τn+1−sϑ−1gsds−δτϑ∑i=0n+1bigτi≤12!Γϑ∑i=0n∫τiτi+1τn+1−sϑ−1g′′ξis−τis−τi+1ds≤Cδτ22!Γϑ∫τ0τn+1τn+1−sϑ−1ds≤TϑCδτ22!Γϑ+1.
From above lemma it is clear that the truncation error of the scheme of equation (33) is Oδτ2.Lemma 4.2.
For0<ϑ<1, the weights bi, i=0,1,2,⋯,n−2 given in equation (33) have the following order of magnitude:
(36)bi=On−iϑ−1.Proof.
According to the definition ofbii=1,2,⋯,n−2, it deduces that
(37)bi=n−iϑ−1Γϑ+3−2ϑ+2n−i21+1n−iϑ+1+ϑ+2n−i21+2n−iϑ+1+ϑ+2n−i2,bi=ϑ+2n−iϑ−1Γϑ+3−2n−i2+ϑ+1n−i+ϑ+1ϑ2!+ϑ+1ϑϑ−13!1n−i+ϑ+1ϑϑ−1ϑ−24!1n−i2+⋯+n−i2+2ϑ+1n−i2+ϑϑ+1+ϑ+1ϑϑ−13!2n−i+ϑ+1ϑϑ−1ϑ−24!2n−i2+⋯+n−i2,implies
(38)bi=ϑ+2n−iϑ−1Γϑ+32ϑ+1ϑϑ−1ϑ−24!1n−i2+6ϑ+1ϑϑ−1ϑ−2ϑ−35!1n−i3+14ϑ+1ϑϑ−1ϑ−2ϑ−3ϑ−46!1n−i4+⋯.
Note that0<ϑ<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 bii=1,2,⋯,n−2. Similar strategy can be adopted for bn, so we neglect the proof.Theorem 4.1.
The numerical scheme (33) is conditionally convergent. That is, for 0<ϑ≤1 and a sufficiently small δτ, there exist a constant C1 such that
(39)xτn+1−xn+1≤C1δτ2.Proof.
To prove the theorem, considering the difference between actual solution and approximate solution, we get(40)xτn+1−xn+1=1−ϑJϑfxτn+1,τn+1−fxn+1,τn+1+ϑJϑΓϑ∫0τn+1τn+1−sϑ−1fxs,sds−ϑδτϑJϑ∑i=0n+1bn+1,ifxi,τi=1−ϑJϑfxτn+1,τn+1−fxn+1,τn+1+ϑδτϑJϑ∑i=0n+1bn+1,ifxτi,τi−fxi,τi+ϑJϑ1Γϑ∫0τn+1τn+1−sϑ−1fxs,sds−δτϑ∑i=0n+1bn+1,ifxτi,τi,xτn+1−xn+1=1−ϑJϑfxτn+1,τn+1−fxn+1,τn+1+ϑδτϑJϑ∑i=0n+1fxτi,τi−fxi,τi+ϑJϑ1Γϑ∫0τn+1τn+1−sϑ−1fxs,sds−δτϑ∑i=0n+1bn+1,ifxτi,τi.
Using Lipschitz condition of equation (23) after applying triangle inequality and using Lemma 4.1, we get
(41)xτn+1−xn+1≤1−ϑJϑfxτn+1,τn+1−fxn+1,τn+1+ϑδτϑJϑ∑i=0n+1bn+1,ifxτi,τi−fxi,τi+ϑJϑ1Γϑ∫0τn+1τn+1−sϑ−1fxs,sds−δτϑ∑i=0n+1bn+1,ifxτi,τi≤1−ϑLJϑxτn+1−xn+1+ϑTϑJϑCδτ22Γϑ+1+ϑδτϑLJϑ∑i=0n+1bn+1,ixτi−xi≤LJϑ1−ϑ+ϑδτϑΓϑ+3xτn+1−xn+1+ϑTϑJϑCδτ22Γϑ+1+ϑδτϑLJϑ∑i=0nbn+1,ixτi−xi.
Now using Lemma4.2 under the condition L/Jϑ1−ϑ+ϑδτϑ/Γϑ+3<1, we get
(42)xτn+1−xn+1≤gδτ,ϑϑTϑJϑCδτ22Γϑ+1+gδτ,ϑϑδτϑLJϑ∑i=0nn−iϑ−1xτi−xi,where gδτ,ϑ=JϑΓϑ+3/JϑΓϑ+3−LΓϑ+31−ϑ+ϑδτϑ For any 0<ϑ<1 and a sufficiently small δτ, there exist a constant Cg such that
(43)1<gδτ,ϑ<Cg.
Hence by applying Gronwall inequality given in Lemma2.1 and using the property (43), we conclude
(44)xτn+1−xn+1≤C1δτ2.Theorem 4.2.
For0<ϑ≤1 and a sufficiently small δτ, there exist a constant C1 such that the numerical scheme (33) is stable.Proof.
To prove the theorem, let the perturbation ofxn and x0 be denoted by x~n and x~0, respectively, for n=0,1,⋯,N−1. Adding this perturbation to the numerical scheme (33) we get
(45)xn+1+x~n+1=x0+x~0+1−ϑJϑfxn+1+x~n+1+ϑδτϑJϑ∑i=0n+1bn+1,ifxi+x~i.
Subtracting equation (33) from equation (45), we get
(46)x~n+1=x~0+1−ϑJϑfxn+1+x~n+1−fxn+1+ϑδτϑJϑ∑i=0n+1bn+1,ifxi+x~i−fxi.
Using the Lipschitz condition followed by triangle inequality, we get(47)x~n+1=x~0+1−ϑJϑfxn+1+x~n+1−fxn+1+ϑδτϑJϑ∑i=0n+1bn+1,ifxi+x~i−fxi,≤x~0+1−ϑJϑfxn+1+x~n+1−fxn+1+ϑδτϑJϑ∑i=0n+1bn+1,ifxi+x~i−fxi,≤x~0+1−ϑLJϑx~n+1+ϑδτϑLJϑ∑i=0n+1bn+1,ix~i,≤x~0+LJϑ1−ϑ+ϑδτϑΓϑ+3x~n+1+ϑδτϑLJϑ∑i=0nbn+1,ix~i.
Now ifL/Jϑ1−ϑ+ϑδτϑ/Γϑ+3<1, then using Lemma 4.2 we get
(48)x~n+1≤gδτ,ϑx~0+gδτ,ϑϑδτϑLJϑ∑i=0nn−iϑ−1x~i,where gδτ,ϑ=JϑΓϑ+3/JϑΓϑ+3−LΓϑ+31−ϑ+ϑδτϑ.
Using property (43), we have
(49)x~n+1≤Cgx~0+CgϑδτϑLJϑ∑i=0nn−iϑ−1x~i.
Hence by applying the Gronwall inequality given in Lemma2.1, we conclude
(50)x~n+1≤C1x~0.


## 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,
(51)DABCτϑSτ=Λ−μS−δS−βSIN,ABCDτϑEτ=βSIN−μ+εE+σβVIN,ABCDτϑIτ=εE−γ+μ+ϑI,ABCDτϑRτ=γI−μR,ABCDτϑVτ=δS−μV−σβVIN,along with initial conditions
(52)S0=S0,E0=E0,I0=I0,R0=R0,V0=V0.Adopting the previously defined procedure, model of equation (51) takes the form
(53)DABCτϑSτ=X1t,S,E,I,R,V,ABCDτϑEτ=X2t,S,E,I,R,V,DABCτϑIτ=X3t,S,E,I,R,V,DABCτϑRτ=X4t,S,E,I,R,V,DABCτϑVτ=X5t,S,E,I,R,V.Furthermore applying the numerical scheme(33), system (53) becomes
(54)Sτn+1=S0+1−ϑJϑX1τ,S,E,I,R,V+ϑδτϑJϑ∑i=0n+1bn+1,iX1τ,S,E,I,R,V,Eτn+1=E0+1−ϑJϑX2τ,S,E,I,R,V+ϑδτϑJϑ∑i=0n+1bn+1,iX2τ,S,E,I,R,V,Iτn+1=I0+1−ϑJϑX3τ,S,E,I,R,V+ϑδτϑJϑ∑i=0n+1bn+1,iX3τ,S,E,I,R,V,Rτn+1=R0+1−ϑJϑX4τ,S,E,I,R,V+ϑδτϑJϑ∑i=0n+1bn+1,iX4τ,S,E,I,R,V,V0+1−ϑJϑX5τ,S,E,I,R,V+ϑδτϑJϑ∑i=0n+1bn+1,iX5τ,S,E,I,R,V,where an+1,k 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 5.1.
Generalized Mean Value Theorem: letϕτ∈Ca,b and Atangana-Baleanu fractional order derivative,ABCDτϑϕτ∈Ca,bfor0<ϑ≤1, then
(55)ϕτ=ϕs+1−ϑJϑABCDϑϕτ+1JϑΓϑABCDϑϕuτ−sϑwith 0≤u≤τ,∀τ∈a,b.Remark 1.
Letϕτ∈C0,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 5.1.
The solution to the nonnegative initial conditions of equation (52) of the proposed model of equation (51) is unique and limited in R+5.Proof.
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)
(56)ABCDτϑSτS=0=Λ,ABCDτϑEτE=0=βSIN+σβVIN,ABCDτϑIτI=0=∈E,ABCDτϑRτR=0=γI,ABCDτϑVτV=0=δS.
From Remark1 and system of equation (56), the solution S0,E0,I0,R0,V0∈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.


## 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 5.1.
Generalized Mean Value Theorem: letϕτ∈Ca,b and Atangana-Baleanu fractional order derivative,ABCDτϑϕτ∈Ca,bfor0<ϑ≤1, then
(55)ϕτ=ϕs+1−ϑJϑABCDϑϕτ+1JϑΓϑABCDϑϕuτ−sϑwith 0≤u≤τ,∀τ∈a,b.Remark 1.
Letϕτ∈C0,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 5.1.
The solution to the nonnegative initial conditions of equation (52) of the proposed model of equation (51) is unique and limited in R+5.Proof.
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)
(56)ABCDτϑSτS=0=Λ,ABCDτϑEτE=0=βSIN+σβVIN,ABCDτϑIτI=0=∈E,ABCDτϑRτR=0=γI,ABCDτϑVτV=0=δS.
From Remark1 and system of equation (56), the solution S0,E0,I0,R0,V0∈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 forCOVID-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 17th 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 ofCOVID-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.
ParameterValueParameterValueInitial time model variablesValueN60000000ε1/4.02S059988539μ0.0000373Λμ×NE011460γ1/5.25ϑ0.00144I01β0.75σ0.25R00Figure4 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 23rd February 2020 to 13th 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 1st July 2021 to 19th September 2021 along with administered vaccination effect considering the initial values of the compartments as of 1st 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 theCOVID-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.



---

*Source: 3754051-2022-09-20.xml*