A Detailed Mathematical Analysis of the Vaccination Model for COVID-19

1  Introduction

Coronaviruses have three branches known as alpha, beta, and gamma. Recently, various strains of SARS-CoV-2 have emerged, including the most destructive and most dangerous delta variant, SARS-CoV-2 [1]. Human coronaviruses have been first distinguished during the 1960s. The first case reported of COVID-19 in Wuhan City in China’s Hubei Province on December 31, 2019, has been found to be contaminated with a new COVID-19 that has never been seen. As reported [2], this disease is believed to be transmitted from animal creatures to people, like SARS-CoV and MERS-CoV. Now the disease is sent from one individual to another. On March 23, 2022, according to the global statistics of this epidemic, the number of infected cases was estimated at 480,165,010, while the number of people recovered reached 414,597,310. By this date, there were 6,144,249 deaths related to the disease worldwide [3]. In this regard, an investigation of the patients who died found that most of them were elderly or patients who had been diagnosed with chronic diseases such as heart disease, diabetes, lung and kidney disease, etc.

Despite the different non-pharmaceutical control systems against COVID-19, a portion of the vaccines that have acquired crisis use approval (EUA) by the United States Centers for Disease Control and Prevention (CDC) are the Pfizer-BioNTech with 95% efficacy, the Moderna immunization with 94.5% viability, and the Janssen vaccine with 67% adequacy, and many others, see [4]. The vaccines have been successful against SARS-CoV-2 infections, including asymptomatic infections and symptomatic cases, severe COVID-19 disease, and deaths.

KSA began vaccinating people against the Coronavirus on December 17, 2020. On November 03, 2021, the use of the Pfizer vaccine was approved for the age groups 5–11 years. In the early stages of the vaccine, citizens and residents over 65 were targeted, which ended on February 18, 2021, whereas the second phase was launched on February 18, 2021, it targeted citizens and residents over 50 years old. The third phase targeted all citizens and residents wishing to receive the vaccine.

On the other hand, time series models are employed in the estimate process of variable behavior and other phenomena, as well as their future trends in diseases, epidemiology, climate sciences, economics, management, and other sciences. Lately, mathematical and statistical modeling have been used to forecast the behavior and patterns of some epidemics and diseases. The prediction method for these models includes four fundamental procedures, such as determining the model, estimation of unknown parameters, diagnostic process, and prediction process. The models (MA, AR, ARCH, ARMA, GARCH, ARIMA) are among the most commonly used time series models for forecasting [5–8].

In this context, the processes of making a mathematical model of a problem, interpreting the solution, validating the model, then making the model ready for utilization are not processes that can be conquered directly. Particularly these days, large numbers of the problems are complicated, nonlinear, have memory impact, or possess a stochastic construction; therefore, modeling methods and solutions specific to these problems must be developed.

Fractional calculus has shown gigantic improvement in applications to various real-world problems in different fields such as continuum mechanics, electromagnetic theory, and biological mathematics. Fractional calculus has become a substantial mathematical tool for the investigation of nonlinear derivative problems, see [9–13]. In modeling biological systems, fractional derivatives play a significant role in considering the nonlocality and memory impact properties that perfectly fit the test data of memory phenomena in various disciplines, such as mechanics, epidemiology, and psychology. The memory effect clarifies that the future condition of the fractional operator of a certain function relies upon its recorded behavior and present status. By using fractional derivatives, especially Caputo fractional derivatives, several real-world systems have been studied successfully in biomathematics and engineering [14–17]. The fractional derivative’s definition has different methodologies with various kernels, so the researchers are keen on picking the best one. Specialists are drawn in by the nonsingular kernels as the classical singular kernels experienced difficulties in modeling some physical phenomena. Caputo et al. [18] overcome this trouble by expanding the Caputo fractional derivative to a nonsingular kernel. There are some interesting properties of the Caputo-Fabrizio operator that have been discussed by Losada et al. in [19]. ABC fractional derivative was given by Atangana et al. in [20], they considered the generalized Mittag-Leffler function as a nonsingular and nonlocal kernel. Some of the generalized characteristics of the ABC operator are introduced in [21,22]. This new operator commanded the attention of investigators because of its brilliant memory description [23,24]. Recent mathematical models on COVID-19 under the ABC operator show different approaches to manifest the transmission of the disease, see [25–32]. In this regard, other papers dealt extensively with modeling and analysis of the transmission of the COVID-19 pandemic in order to limit the spread of this epidemic, see [33–38]. Diagne et al. [39] formulated a mathematical model of the COVID-19 transmission mechanism with Vaccination and Treatment. A mathematical model of COVID-19 with the effect of vaccines was constructed by Gokbulut et al. [40]. They made the analysis of the model, then the idea of vaccination was used with a change in the vaccination rate among the population. Yavuz et al. [41] developed a mathematical model to reveal the effects of vaccine treatment on COVID-19 based on a system of ordinary differential equations.

The literature shows that very few works available on COVID-19 models are in the form of ABC fractional derivatives, even though they describe the dynamic behavior accurately compared to the classical derivatives. To our knowledge, there is no available literature on the COVID-19 fractional model with an expectation of the vaccination effect on KSA. Therefore, this is the main motivation behind our work. Also, the comparison of the real data with different fractional order simulations is one of our unique aims in this work. In addition, we extend and generalize the model studied by Yavuz et al. [41], based on a system of fractional differential equations involving ABC fractional derivatives.

This work is organized as follows. Section 2 provides a statistical analysis of COVID-19 vaccination in KSA. Some background material and auxiliary results of ABC fractional calculus are given in Section 3. The fractional mathematical model and its fundamental properties are described in Section 4. In Section 5, we prove the theoretical results of the given model. Sections 6 and 7 give numerical results and simulation results by estimating the parameters. In the last section, we provide brief conclusions of the study.

2  Statistical Analysis

In this section, we discuss some empirical results. First, we analyze and present the descriptive data of the vaccinated persons in KSA. Then we predict the number of people who will be fully vaccinated during the next period.

2.1 Descriptive Data Analysis

The latest statistics about the vaccination of people against the spread of the Coronavirus epidemic are collected from the official data published by the Saudi Ministry of Health as the fully vaccinated people reached 144,057 during the period January 06, 2021, to June 27, 2021, with an average of 8,401 per day. Then, data is collected on the number of fully vaccinated people daily during the period June 28, 2021, to February 15, 2022. Table 1 indicates that the number of vaccinated people in July 2021 reached 6,493,739, with an average of 209,476 people per day. Then it increased in August 2021 to reach 6,833,828 with an average of 220,446 people per day. Hence, we observe that the number of vaccinated people decreased continuously until the number of vaccinated people during the fifteen days of February 2022 reached 264,633. In general, the number of fully vaccinated people during the studied period reached 23,964,566, which represents 67.81% of the total population in KSA, which is 35,340,680 people. For more details, see Figs. 1 and 2.

Figure 1: Number of people fully vaccinated in the KSA during the period June 01, 2021 to February 15, 2022

Figure 2: Daily numbers of people completed vaccinations in KSA from June 28, 2021 to February 15, 2022

Table 2 summarizes the number of daily doses taken by people against coronavirus disease in KSA during the period June 28, 2021, to February 15, 2022. We find that the doses increased continuously from January 2021 to April 2021, then decreased in May 2021 and June 2021. After that, it increased dramatically in July 2021 and August 2021 to reach 19,252,318 doses with 33% of the total doses given during the study period. Then it decreased during the period from September 2021 to December 2021. The total doses given amounted to 57,964,264 during the studied period, with an average of 248,773 doses per day.

2.2 Prediction Strategies

In this portion, we are collecting data on the fully vaccinated people in KSA for the period June 01, 2021 to February 15, 2022. Moreover, the stability tests of these data are examined to use the prediction process as unit root tests and estimation of coefficients (ACF & APCF). Phillips-Peron’s and Dickey-Fuller’s tests show that the data series is unstable at the level, which means that there is a general trend in the series. As shown in Fig. 2. In this regard, we process the data by converting it to the logarithm and taking the first differences to remove the effect of the general trend from the series to be used in estimation and prediction, see Fig. 3. Furthermore, we repeat the unit root tests as in Table 3. Hence, we find that the calculated values are more significant than the critical values and Prob. = 0.0000, which is less than 5%, this indicates that the data series does not have a unit root and has become stable.

Figure 3: Transformation of the data of completely vaccinated people in KSA to the first difference

We apply ARIMA models for forecasting through the Eviews program on the data series about the number of vaccinated people completely in KSA. As well, we estimate ARIMA models by using the Ordinary Least Squares method, as shown in Table 4. On the other hand, the models that have no statistical significance are excluded and the comparison only in the statistically significant models, and we are choosing a proper model that achieves the best models in which the coefficient of determination is greater, less variance, less volatility, and less value to AIC indicator.

Moreover, we check the ARIMA (1, 1, 2) model by testing residuals and the shape of the autocorrelation and partial autocorrelation coefficients. It follows from Fig. 4 that the actual values match the estimated values, and this indicates that the differences were small, as well as the efficiency of the model and its suitability in the forecasting process.

Figure 4: Evaluation of the ARIMA model for actual and fitted values and residual limits

We see through the chosen model a decrease in the number of people who will be fully dosed to (927164) during the predictive period from February 16, 2022, to April 14, 2022. For more details, see Fig. 5 and Table 6.

Figure 5: Prediction of the number of fully vaccinated people in KSA with confidence intervals (95%) for the period February 16, 2022 to April 15, 2022

Finally, we estimate a linear model using the least-squares method to test the predictive ability of the model. Indeed, we take the actual values as a dependent variable and the estimated values as an independent variable. We conclude that the closer the estimated parameter to one, the more closely the estimated values are to the actual values. Through the results of Table 5, it is clear that the predicted parameter is close to one (0.955688), and this indicates the quality of the model in prediction and that there is a convergence of the predicted values from the actual values and is statistically significant. Additionally, we find that the value of F = 1287.919, whereas Prob. = 0.0000 is less than 0.05, which indicates the model is significant and good.

3  Auxiliary Results

Definition 3.1. [20] Let ρ∈[0,1], and ω∈H1(a,b), b > a. Then, the ABC fractional derivative of order ρ is defined by

ABCDa+ρω(t)=ℵ(ρ)1−ρ∫atEρ(ρρ−1(t−θ)ρ)ω′(θ)dθ,t>0,(1)

The normalization function ℵ(ρ) satisfies ℵ(0)=ℵ(1)=1, where Eρ is the Mittag-Leffler function defined by

Eρ(z)=∑k=0∞zkΓ(kρ+1), Re(ρ)>0,z∈C.

Definition 3.2. [20] The correspondent AB fractional integral is given by

ABℑa+ρω(t)=1−ρℵ(ρ)ω(t)+ρℵ(ρ)Γ(ρ)∫at(t−θ)ρ−1ω(θ)dθ,ρ∈[0,1].

Lemma 3.1. [22] Let ω(t) be a function defined on [a,b] and 0<ρ≤1. Then

ABABCℑa+ρABCDa+ρω(t)=ω(t)−ω(a).

Lemma 3.2. [20] The Laplace transform of the ABC fractional derivative is defined by

LABCD0+ρω(t)=ℵ(ρ)1−ρsρL{ω(t)}(s)−sρ−1ω(0)sρ+ρ1−ρ.

Lemma 3.3. [20] The fractional differential system ABCD0+ρω(t)=f(t), ω(0)=ω0, gives the unique solution

ω(t)=ω0+1−ρℵ(ρ)f(t)+ρℵ(ρ)Γ(ρ)∫at(t−θ)ρ−1f(θ)dθ,ρ∈[0,1].

4  Model Derivation in ABC Operator

Yavuz et al. [41] developed a mathematical model to reveal the effects of vaccine treatment on COVID-19 described by a system of ODEs as follows:

{S⋅(ϰ)=Λ−(αE(ϰ)+m+μ)S(ϰ),E⋅(ϰ)=αE(ϰ)S(ϰ)+pV(ϰ)E(ϰ)−(κI(ϰ)+c+μ)E(ϰ),I⋅(ϰ)=κI(ϰ)E(ϰ)−(z+μ+σ)I(ϰ),V⋅(ϰ)=mS(ϰ)−(pE(ϰ)+μ)V(ϰ),R⋅(ϰ)=zI(ϰ)+cE(ϰ)−μR(ϰ).(2)

With initial conditions

S(0)=S0,E(0)=E0,I(0)=I0,V(0)=V0,R(0)=R0.

In the present work, we make extend the model (2) by replacing the time derivative with the ABC fractional derivative. With this alteration, the right-and left-hand sides will not have the same dimensions. To conquer this affair, we add an auxiliary parameter K with the dimension of s, to change the nonlocal fractional operator so that the sides possess the same dimension [43]. Herewith, we reformulate the following fractional system:

{1K1−ΘABCDΘS(ϰ)=Λ−(αE(ϰ)+m+μ)S(ϰ),1K1−ΘABCDΘE(ϰ)=αE(ϰ)S(ϰ)+pV(ϰ)E(ϰ)−(κI(ϰ)+c+μ)E(ϰ),1K1−ΘABCDΘI(ϰ)=κI(ϰ)E(ϰ)−(z+μ+σ)I(ϰ),1K1−ΘABCDΘV(ϰ)=mS(ϰ)−(pE(ϰ)+μ)V(ϰ),1K1−ΘABCDΘR(ϰ)=zI(ϰ)+cE(ϰ)−μR(ϰ).(3)

With initial conditions

S(0)=S0≥0,E(0)=E0≥0,I(0)=I0≥0,V(0)=V0≥0,R(0)=R0≥0,(4)

where ABCDΘ is the ABC fractional derivative of order Θ∈(0,1], the functions S(ϰ), E(ϰ), I(ϰ), V(ϰ), R(ϰ) and their fractional derivatives are continuous at ϰ ≥ 0, and (S,E,I,V,R)∈Λ5, where Λ:=C([0,T],R) is a Banach space endowed with the norm ‖X‖:=supt∈[0,T]{|X(t)|}. Furthermore, we set Ω:=(Λ5,‖X‖) with the norm ||X||=||S||+||E||+||I||+||V||+||R||.

Here S(ϰ) is the class of Susceptible Individuals, E(ϰ) is the class of Exposed Individuals, I(ϰ) is the class of Infected individuals, V(ϰ) is the class of Individuals Vaccinated, R(ϰ) is the class of Recovered Individuals. The model parameter values and source are given in Table 7.

Total population at time ϰ, denoted by N(ϰ) and given by N(ϰ) = S(ϰ) + E(ϰ) + I(ϰ) + V(ϰ) + R(ϰ).

4.1 Positivity and Boundedness

In this portion, we show that the non-negative domain R5 is positively invariant region and the solutions of model are bounded, where R+5={Φ∈R5:Φ(ϰ)≥0}.

Theorem 4.1. There is a unique solution Φ(ϰ)=(S(ϰ),E(ϰ),I(ϰ),V(ϰ),R(ϰ))T for the models (3) and (4) on ϰ ≥ 0 in [0, T] (0<T<∞) and the solution will remain in R5. Furthermore, the solution is stated in the region Δ⊂R+5, described by

Δ={(S(ϰ),E(ϰ),I(ϰ),V(ϰ),R(ϰ))∈R5:N(ϰ)≤Λμ}.

Proof. From the model (3), we have

ABCDΘS(ϰ)=K1−Θ(Λ−(αE(ϰ)+m+μ)S(ϰ)),ABCDΘE(ϰ)=K1−Θ(αE(ϰ)S(ϰ)+pV(ϰ)E(ϰ)−(κI(ϰ)+c+μ)E(ϰ)),ABCDΘI(ϰ)=K1−Θ(κI(ϰ)E(ϰ)−(z+μ+σ)I(ϰ)).ABCDΘV(ϰ)=K1−Θ(mS(ϰ)−(pE(ϰ)+μ)V(ϰ)).ABCDΘR(ϰ)=K1−Θ(zI(ϰ)+cE(ϰ)−μR(ϰ)).

The norm and all assumptions of the classical results are valid. It follows that:

{ABCDΘS(ϰ)≥−K1−Θ(αE(ϰ)+m+μ)S(ϰ),ABCDΘE(ϰ)≥−K1−Θ(κI(ϰ)+c+μ)E(ϰ),ABCDΘI(ϰ)≥−K1−Θ(z+μ+σ)I(ϰ)ABCDΘV(ϰ)≥−K1−Θ(pE(ϰ)+μ)V(ϰ)ABCDΘR(ϰ)≥−K1−ΘμR(ϰ)

For K>0 and ϰ > 0, we have

{S(ϰ)≥S0EΘ(−ΘK1−Θ(α||E||+m+μ)ϰΘℵ(Θ)−(1−Θ)K1−Θ(α||E||+m+μ)),E(ϰ)≥E0EΘ(−ΘK1−Θ(κ||I||+c+μ)ϰΘℵ(Θ)−(1−Θ)K1−Θ(κ||I||+c+μ)),I(ϰ)≥I0EΘ(−ΘK1−Θ(z+μ+σ)ϰΘℵ(Θ)−(1−Θ)K1−Θ(z+μ+σ)),V(ϰ)≥V0EΘ(−ΘK1−Θ(p||E||+μ)ϰΘℵ(Θ)−(1−Θ)K1−Θ(p||E||+μ))R(ϰ)≥R0EΘ(−ΘK1−ΘμϰΘℵ(Θ)−(1−Θ)K1−Θμ).

This shows that if S0,E0,I0,V0, and R0 are positive, then S(ϰ), E(ϰ), I(ϰ), V(ϰ), and R(ϰ) are also positive. Now, we give a biologically feasible region for the suggested models (3) and (4). Consider summing all the equations in the model as follows

1K1−ΘABCDΘN(ϰ)=1K1−ΘABCDΘ(S(ϰ)+E(ϰ)+I(ϰ)+V(ϰ)+R(ϰ)).

It follows for the whole population that

ABCDΘN(ϰ)=K1−Θ(Λ−μN(ϰ)−σI(ϰ))≤K1−Θ(Λ−μN(ϰ)),

which implies that

ABCDΘN(ϰ)≤K1−Θ(Λ−μN(ϰ)).

After applying the Laplace transform, we have

N(s)≤(K1−Θ(1−Θ)Λϑ+ℵ(Θ)ϑN(0))sΘ−1(sΘ+λ)+K1−ΘΘΛϑs−1(sΘ+λ),

where ϑ:=[ℵ(Θ)+K1−Θ(1−Θ)μ], and λ=K1−ΘΘμ[ℵ(Θ)+K1−Θ(1−Θ)μ]. Using the inverse Laplace, we get

N(ϰ)=(K1−Θ(1−Θ)Λϑ+ℵ(Θ)N(0)ϑ)EΘ,1(−λϰΘ)+K1−ΘΘΛϑϰΘEΘ,Θ+1(−λϰΘ),

where Eα,β(z)=∑k=0∞zkΓ(kα+β), and L[ϰβ−1Eα,β(±aϰα)]=sα−β(sα∓a). Subsequently, the solution of the model with the nonnegative conditions in Δ stays in Δ. In this way, the region Δ is positively invariant and attracts all the solutions in R5. This shows that the models (3) and (4) is epidemiologically well posed and the solutions.

4.2 Equilibrium Points and Reproduction Number

In this part, we will find the equilibrium points of the COVID-19 model. By equating each equation of model (3) to the zero, we can write

0=K1−Θ(Λ−(αE+m+μ)S),0=K1−Θ(αES+pVE−(κI+c+μ)E),0=K1−Θ(κIE−(z+μ+σ)I).0=K1−Θ(mS−(pE+μ)V).0=K1−Θ(zI+cE−μR).

For K>0, we have

0=Λ−(αE+m+μ)S,0=αES+pVE−(κI+c+μ)E,0=κIE−(z+μ+σ)I.0=mS−(pE+μ)V.0=zI+cE−μR.

Hence, the disease free equilibrium (DFE) and endemic equilibrium (EE) are given by the following theorem.

Theorem 4.2. We have the next affirmations:

i. The model (3) has always DFE, E0=(S0,E0,I0,V0,R0), where S0=Λm+μ,E0=0,I0=0,V0=Λm(m+μ)μ, and R0=0.

ii. The model (3) has EE, E∗=(S∗,E∗,I∗,V∗,R∗), where

S∗=κΛαμ+ασ+κμ+κm+αzE∗=σ+μ+zκI∗=κ(Λ(αμ+mp)+αp(σ+μ+z))(κμ+pσ+pμ+pz)(κm+κμ+ασ+αμ+αz)−c+μκV∗=κ2Λm(κμ+pσ+pμ+pz)(κm+κμ+ασ+αμ+αz)R∗=κ2(m(μ(c(μ+σ)−μz)+Λpz)+μ(cμ(μ+σ)+z(αΛ−μ2)))κμ(κμ+pσ+pμ+pz)(κm+κμ+ασ+αμ+αz)+κ(μ+σ+z)+p(c(μ+σ)(μ+m)+z(αΛ−μ(μ+m)))κμ(κμ+pσ+pμ+pz)(κm+κμ+ασ+αμ+αz)+αμ(c(μ+σ)−μz)+αp(σ+μ+z)2(c(μ+σ)−μz)κμ(κμ+pσ+pμ+pz)(κm+κμ+ασ+αμ+αz).

By using the model (3) and next generation matrix method, the basic reproduction number R0 calculated as

R0=Λ(μα+mp)μ(μ+c)(μ+m).

For details, see [41].

4.3 Stability Analysis

Stability analysis of the equilibrium points is clarified in [41], exhaustively. In order of stability, the authors dealt with the next two theorems.

Theorem 4.3. The DFE point E0 of the epidemic model is locally asymptotically stable if R0<1, otherwise unstable.

Theorem 4.4. The EE point E∗ of the epidemic model is locally asymptotically stable if R0>1, otherwise unstable.

Here we comment that R0 gives us data about the transmission of infectious sickness. For example in the normal utilization of contamination models, when R0>1, the disease infection will actually be able to begin transmitting in society, yet if R0<1 it implies that the disease infection continues declining in the general public. Overall, the bigger the worth of R0 the troublesome it is to control the pandemic from spreading in the public arena.

5  Existence Criteria

In this section, we apply the Picard-Lindel method and the Laplace transform to investigate the existence and uniqueness of solution for preventive and curative to fractional COVID-19 disease model.

5.1 Iterative Scheme with Laplace Transform

Theorem 5.1. For 0<Θ≤1, the following ABC-type FDE

1K1−αABCDΘu(ϰ)=f(ϰ)(5)

has a unique solution, which is

u(ϰ)=u(0)+(1−Θ)κ1−Θℵ(Θ)f(ϰ)+Θκ1−Θℵ(Θ)Γ(Θ)∫0ϰ(ϰ−ν)Θ−1f(ν)dν.(6)

Proof. Applying the Laplace transform on both sides of Eq. (5) we obtain

L{1K1−ΘABCDΘu(ϰ)}(ϖ)=L{f(ϰ)}(ϖ),ϖ>0.

It follows from Theorem 3 in [20] that

ℵ(Θ)1−ΘϖΘL(u(ϰ))(ϖ)−ϖΘ−1u(0)ϖΘ+Θ1−Θ=K1−ΘL{f(ϰ)}(ϖ),

which is equivalent to

L(u(ϰ))(ϖ)=1ϖu(0)+(1−Θ)K1−Θℵ(Θ)L{f(ϰ)}(ϖ)+ΘK1−ΘϖΘℵ(Θ)L{f(ϰ)}(ϖ).

Applying the inverse Laplace transform give us

u(ϰ)=u(0)+(1−Θ)K1−Θℵ(Θ)f(ϰ)+L−1{ΘK1−ΘϖΘℵ(Θ)L{f(ϰ)}(ϖ)}(ϰ).

Now, we have

ΘK1−Θℵ(Θ)(1ϖΘ)=ΘK1−Θℵ(Θ)L{ϰΘ−1Γ(Θ)}(ϖ).

Let F(ϖ)=ΘK1−Θℵ(Θ)L{ϰΘ−1Γ(Θ)}(ϖ), and G(ϖ)=L{f(ϰ)}(ϖ). It follows form convolution theorem that

L−1{ΘK1−ΘϖΘℵ(Θ)L{f(ϰ)}(ϖ)}(ϰ)=L−1{F(ϖ)×G(ϖ)}(ϰ)=ΘK1−Θℵ(Θ){ϰΘ−1Γ(Θ)∗f(ϰ)}(ϖ)=ΘK1−Θℵ(Θ)Γ(Θ)∫0ϰ(ϰ−ϖ)Θ−1f(ϖ)dϖ.

Hence

u(ϰ)=u(0)+(1−Θ)K1−Θℵ(Θ)f(ϰ)+ΘK1−Θℵ(Θ)Γ(Θ)∫0ϰ(ϰ−ν)Θ−1f(ν)dν.

By using Theorem 5.1, our model is equivalent to

{S(ϰ)−S(0)=(1−Θ)K1−Θℵ(Θ){Λ−(αE(ϰ)+m+μ)S(ϰ)}+ΘK1−Θℵ(Θ)Γ(Θ)∫0ϰ(ϰ−r)Θ−1{Λ−(αE(r)+m+μ)S(r)}drE(ϰ)−E(0)=(1−Θ)K1−Θℵ(Θ){αE(ϰ)S(ϰ)+pV(ϰ)E(ϰ)−(κI(ϰ)+c+μ)E(ϰ)}+ΘK1−Θℵ(Θ)Γ(Θ)∫0ϰ(ϰ−r)Θ−1{αE(r)S(r)+pV(r)E(r)−(κI(r)+c+μ)E(r)}drI(ϰ)−I(0)=(1−Θ)K1−Θℵ(Θ){κI(ϰ)E(ϰ)−(z+μ+σ)I(ϰ)}+ΘK1−Θℵ(Θ)Γ(Θ)∫0ϰ(ϰ−r)Θ−1{κI(r)E(r)−(z+μ+σ)I(r)}drV(ϰ)−V(0)=(1−Θ)K1−Θℵ(Θ){mS(ϰ)−(pE(ϰ)+μ)V(ϰ)}+ΘK1−Θℵ(Θ)Γ(Θ)∫0ϰ(ϰ−r)Θ−1{mS(r)−(pE(r)+μ)V(r)}drR(ϰ)−R(0)=(1−Θ)K1−Θℵ(Θ){zI(ϰ)+cE(ϰ)−μR(ϰ)}+ΘK1−Θℵ(Θ)Γ(Θ)∫0ϰ(ϰ−r)Θ−1{zI(r)+cE(r)−μR(r)}dr(7)

The iterative scheme of the model (7) is given by

S0(ϰ)=S(0),E0(ϰ)=E(0),I0(ϰ)=I(0),V0(ϰ)=V(0),R0(ϰ)=R(0).

{Sn+1(ϰ)=(1−Θ)K1−Θℵ(Θ){Λ−(αEn(ϰ)+m+μ)Sn(ϰ)}+ΘK1−Θℵ(Θ)Γ(Θ)∫0ϰ(ϰ−r)Θ−1{Λ−(αEn(r)+m+μ)Sn(r)}drEn+1(ϰ)=(1−Θ)K1−Θℵ(Θ){αEn(ϰ)Sn(ϰ)+pVn(ϰ)En(ϰ)−(κIn(ϰ)+c+μ)En(ϰ)}+ΘK1−Θℵ(Θ)Γ(Θ)∫0ϰ(ϰ−r)Θ−1{αEn(r)Sn(r)+pVn(r)En(r)−(κIn(r)+c+μ)En(r)}drIn+1(ϰ)=(1−Θ)K1−Θℵ(Θ){κIn(ϰ)En(ϰ)−(z+μ+σ)In(ϰ)}+ΘK1−Θℵ(Θ)Γ(Θ)∫0ϰ(ϰ−r)Θ−1{κIn(r)En(r)−(z+μ+σ)In(r)}drVn+1(ϰ)=(1−Θ)K1−Θℵ(Θ){mSn(ϰ)−(pEn(ϰ)+μ)Vn(ϰ)}+ΘK1−Θℵ(Θ)Γ(Θ)∫0ϰ(ϰ−r)Θ−1{mSn(r)−(pEn(r)+μ)Vn(r)}drRn+1(ϰ)=(1−Θ)K1−Θℵ(Θ){zIn(ϰ)+cEn(ϰ)−μRn(ϰ)}+ΘK1−Θℵ(Θ)Γ(Θ)∫0ϰ(ϰ−r)Θ−1{zIn(r)+cEn(r)−μRn(r)}dr(8)

Taking the limit as n→∞, we anticipate getting the exact solution.

5.2 Existence of a Unique Solution

In this portion, we apply the Picard-Lindel method to investigate the existence of solution for the fractional COVID-19 disease model, which its mathematical represente are presented by:

{1K1−ΘABCDΘS(ϰ)=Λ−(αE(ϰ)+m+μ)S(ϰ),1K1−ΘABCDΘE(ϰ)=αE(ϰ)S(ϰ)+pV(ϰ)E(ϰ)−(κI(ϰ)+c+μ)E(ϰ),1K1−ΘABCDΘI(ϰ)=κI(ϰ)E(ϰ)−(z+μ+σ)I(ϰ).1K1−ΘABCDΘV(ϰ)=mS(ϰ)−(pE(ϰ)+μ)V(ϰ).1K1−ΘABCDΘR(ϰ)=zI(ϰ)+cE(ϰ)−μR(ϰ).(9)

For sake of simpilicity, we define the functions Ψi, i = 1, 2,…, 5, as follows:

{Ψ1(ϰ,S(ϰ))=Λ−(αE(ϰ)+m+μ)S(ϰ),Ψ2(ϰ,E(ϰ))=αE(ϰ)S(ϰ)+pV(ϰ)E(ϰ)−(κI(ϰ)+c+μ)E(ϰ),Ψ3(ϰ,I(ϰ))=κI(ϰ)E(ϰ)−(z+μ+σ)I(ϰ),Ψ4(ϰ,V(ϰ))=mS(ϰ)−(pE(ϰ)+μ)V(ϰ),Ψ5(ϰ,R(ϰ))=zI(ϰ)+cE(ϰ)−μR(ϰ).(10)

The model described in Eq. (9) becomes in the following form:

{1K1−ΘABCDΘS(ϰ)=Ψ1(ϰ,S(ϰ)),1K1−ΘABCDΘE(ϰ)=Ψ2(ϰ,E(ϰ)),1K1−ΘABCDΘI(ϰ)=Ψ3(ϰ,I(ϰ)),1K1−ΘABCDΘV(ϰ)=Ψ4(ϰ,V(ϰ)),1K1−ΘABCDΘR(ϰ)=Ψ5(ϰ,R(ϰ)),(11)