A Fractional COVID-19 Model with Efficacy of Vaccination

Abstract

This paper develops a fractional-order model of COVID-19 with vaccination. The model is well designed by including both the efficacy and inefficacy of vaccinations in humans. Besides calculating the reproduction number, equilibrium points and the feasibility region are also determined. Stability analysis for the proposed model around equilibrium points is discussed. Fixed-point theory is employed to identify the singularity of the solution. Adomian decomposition and Laplace integral transformation are combined to obtain the solution. We present the solutions graphically to analyze the contributions of the disease dynamics based on different values of the fractional order. This study seeks an in-depth understanding of COVID-19 transmission to improve health outcomes.

1. Introduction

The World Health Organization issued a public health emergency of international concern on 30 January 2020 and a pandemic alert in March 2020 in response to the rapid and extensive spread of Coronavirus Disease 2019. Initially, the disease turned out to be an epidemic in mainland China first hitting Wuhan, in the province of Hubei [1]. India had the second-highest confirmed cases in the world. The direct transmission of SARS-CoV-2 was conclusively proven by clinical evidence in January 2020 [2]. Four consecutive lockdowns came into effect in India as a preventive measure. In the absence of vaccines, social distancing serves as the best control measure against COVID-19 [3].

As a part of the vaccination program, India has established a NEGVAC (National Expert Group on Vaccine Administration for COVID-19) committee to develop guidelines for the COVID-19 vaccine administration [4]. There have been vaccination campaigns conducted around the world to combat COVID-19. Several types of vaccines are available for COVID-19 [5]. Furthermore, India has exported 35,793,000 doses of vaccine internationally (commercial exports). Innovative technologies are tested during the development of vaccines to determine if they work. Vaccine development, which is dependent on public approval, determines the effectiveness of the vaccination program based on perception and beliefs surrounding the vaccine [6]. In addition, a significant portion of the population remains uncertain about the vaccine. The perception of the vaccine among people may change with the advent of awareness programs and the improvement of vaccine outcomes [7].

Fractional-order differential equations are the most effective tool for studying biological and engineering systems. Several fractional-order derivatives comprise numerical models of physical and biological systems [8,9]. There are numerous motivations to employ fractional order, but the main focus is on dealing with the memory dynamics evident in many biological systems. Differential equations derived using the fractional derivative have several applications to analyze various infectious disease transmission dynamics such as HIV/AIDS, TB, and others [10].

The results in [11] indicate that even in the absence of social distancing, only vaccination to people can significantly reduce the overall infected population. Other approaches related to vaccination strategies have been considered, for example, in [12]. Some fractional models have been considered [13] to derive the optimal control of COVID-19 transmission. In [14], the authors proposed the fraction model and dynamics of COVID-19 and numerical simulations are performed using Laguerre collocation technique. In [15], a generalized fractional-order SEIR model is proposed. Specifically, the paper shows that, according to real data from the USA, the considered fractional model has a good prediction ability for the epidemic trend in the next two weeks. In [16], numerical simulations of fractional-order modeling of COVID-19 in the case of Wuhan (China) were carried out. In particular, the Adams–Bashforth numerical scheme was used in simulations of the Caputo–Fabrizio fractional-order derivative.

The aim of this paper is to develop the fractional-order model to describe the dynamics of COVID-19 using the Caputo–Fabrizio fractional derivative. The main contribution of this paper is to prove the existence and uniqueness of the solution using fixed-point theory. In order to reduce infection in a community, this paper examines the dynamic behavior of the system with vaccination. A summary of the paper’s content is as follows. We formulate an SVEITR model for COVID-19 in Section 2. We calculate equilibrium points and the reproduction number in Section 3. Section 4 discusses the existence and uniqueness of the fractional system solution. The Adomian decomposition method with Laplace integral transform is employed to obtain the analytical solution whose graphical results appear in Section 5 and Section 6.

2. Mathematical Model

Let us consider the mathematical model for understanding the dynamics of COVID-19 spread with vaccination. The model is formulated based on assumptions involving six compartments: $S\left(t\right)$ and $V\left(t\right)$ represent the people who are likely to contract the infection and are vaccinated, $E\left(t\right)$ is the exposed population, $I\left(t\right)$ is the infected population, $T\left(t\right)$ is individuals who are receiving treatment after infection, $R\left(t\right)$ is the recovered population. The total population is $N\left(t\right)=S\left(t\right)+V\left(t\right)+E\left(t\right)+I\left(t\right)+T\left(t\right)+R\left(t\right)$. It is formulated as:

$$\begin{array}{ccc}\hfill \frac{dS}{dt}& =& \alpha -{\beta }_{1}SE-{\beta }_{2}SI-(\mu +\psi )S+\varsigma T\hfill \\ \hfill \frac{dV}{dt}& =& \psi S-(\iota +\mu )V+\varpi R\hfill \\ \hfill \frac{dE}{dt}& =& {\beta }_{1}SE+{\beta }_{2}SI-(\sigma +\mu )E\hfill \\ \hfill \frac{dI}{dt}& =& \sigma E-(\eta +d+\mu )I\hfill \\ \hfill \frac{dT}{dt}& =& \eta I-(\gamma +\mu +\varsigma )T\hfill \\ \hfill \frac{dR}{dt}& =& \gamma T-(\mu +\varpi )R\hfill \end{array}$$

(1)

where $\alpha$: the recruitment rate of individuals susceptible to the disease, ${\beta }_1$: the rate at which infectious agents transfer from exposed individuals to susceptible humans,

${\beta }_2$: the rate at which infection is transferred from infected persons to susceptible persons via objects or surfaces used by them,

$\mu$: natural death rate,

$\varsigma$: the removal rate of the treatment population to susceptible population,

$\sigma$: the rate of the exposed population getting infected,

$\eta$: the rate of individuals receiving treatment after knowing symptoms,

d: disease-induced mortality rate,

$\gamma$: recovery rate,

$\varpi$: the removal rate of the treatment population to vaccinated population,

$\psi$: vaccination rate,

$\iota$: rate of vaccine inefficacy in people.

With $S\left(0\right)=S_0,V\left(0\right)=V_0,E\left(0\right)=E_0,I\left(0\right)=I_0,T\left(0\right)=T_0,R\left(0\right)=R_0$ as the initial conditions.

In the model, we apply the definitions of fractional derivatives [17]. The fractional derivatives and integrals are incorporated, and the time derivative takes the form of a Caputo–Fabrizio (CF) fractional derivative [18]. For $g\in H^1(c,d)$ and $d>c$, the Caputo–Fabrizio derivative of fractional order $n\in (0,1)$ for g is given by

$$\begin{array}{c}\hfill {}^{CF}{D}^{n}g\left(t\right)=\frac{M\left(n\right)}{(1-n)}{\int }_c^{t}exp\left(\frac{-n}{1-n}(t-v)\right)g^{\prime }\left(v\right)dv\end{array}$$

where $t\ge 0,M\left(n\right)$ is a normalization function that depends on n and $M\left(0\right)=M\left(1\right)=1$. Hence, we introduce an auxiliary parameter called $\lambda$ for $n\in (0,1)$ [19,20,21]. The transmission model is as follows:

$$\begin{array}{ccc}\hfill \frac{1}{{\lambda }^{n-1}}{}^{CF}{D}_t^{n}S\left(t\right)& =& \alpha -{\beta }_{1}S\left(t\right)E\left(t\right)-{\beta }_{2}S\left(t\right)I\left(t\right)-(\mu +\psi )S\left(t\right)+\varsigma T\left(t\right)\hfill \\ \hfill \frac{1}{{\lambda }^{n-1}}{}^{CF}{D}_t^{n}V\left(t\right)& =& \psi S\left(t\right)-(\iota +\mu )V\left(t\right)-\varpi R\left(t\right)\hfill \\ \hfill \frac{1}{{\lambda }^{n-1}}{}^{CF}{D}_t^{n}E\left(t\right)& =& {\beta }_{1}S\left(t\right)E\left(t\right)+{\beta }_{2}S\left(t\right)I\left(t\right)-(\sigma +\mu )E\left(t\right)\hfill \\ \hfill \frac{1}{{\lambda }^{n-1}}{}^{CF}{D}_t^{n}I\left(t\right)& =& \sigma E\left(t\right)-(\eta +d+\mu )I\left(t\right)\hfill \\ \hfill \frac{1}{{\lambda }^{n-1}}{}^{CF}{D}_t^{n}T\left(t\right)& =& \eta I\left(t\right)-(\gamma +\mu +\varsigma )T\left(t\right)\hfill \\ \hfill \frac{1}{{\lambda }^{n-1}}{}^{CF}{D}_t^{n}R\left(t\right)& =& \gamma T\left(t\right)-(\mu +\varpi )R\left(t\right)\hfill \end{array}$$

(2)

with the same initial conditions.

3. Model Analysis

Theorem 1.

The closed set Ω is positively invariant under fractional system (2).

Proof. 

From (2), we obtain

$$\begin{array}{c}\hfill \frac{1}{{\lambda }^{n-1}}{}^C{D}_t^{n}N\left(t\right)=\alpha -\mu N\left(t\right)\end{array}$$

After some calculations, we have

$$\begin{array}{c}\hfill N\left(t\right)\le \frac{\alpha }{\mu }-\left(\frac{\alpha -\mu N_0}{\mu }\right)e^{-\mu t}\end{array}$$

Thus, $N\left(0\right)\le \frac{\alpha }{\mu }$ then $t>0$, $N\left(t\right)\le \frac{\alpha }{\mu }$. Therefore, the closed set $\Omega$ is positively invariant. □

3.1. Equilibria Points

By equating each relationship in the fractional system with zero, this section determines the equilibrium points of the fractional system. The disease-free equilibrium exists when there is no disease; it is denoted by $E_0$, where $E_0=\left(\frac{\alpha }{\mu +\psi },\frac{\alpha \psi }{\mu (\mu +\psi )},0,0,0,0\right)$.

The endemic equilibrium is denoted by $E^{*}$, $E^{*}=(S^{*},V^{*},E^{*},I^{*},T^{*},R^{*})$, where

$$\begin{array}{ccc}\hfill T^{*}& =& \frac{\eta I^{*}}{\gamma +\mu +\varsigma }\hfill \\ \hfill R^{*}& =& \frac{\gamma \eta I^{*}}{(\mu +\varpi )(\gamma +\mu +\varsigma )}\hfill \\ \hfill E^{*}& =& \frac{(\eta +d+\mu )I^{*}}{\sigma }\hfill \\ \hfill S^{*}& =& \frac{\alpha }{\mu +\psi }+A{I}^{*}\hfill \\ \hfill V^{*}& =& \frac{\alpha \psi (\gamma +\mu +\varsigma )+\left[\psi A(\mu +\varpi )(\gamma +\mu +\varsigma )+\varpi \gamma \eta \right]I^{*}}{(\iota +\mu )(\mu +\varpi )(\gamma +\mu +\varsigma )}\hfill \\ \hfill I^{*}& =& \frac{(\sigma +\mu )(\eta +d+\mu )(\mu +\psi )-{\beta }_1\alpha (\eta +d+\mu )-{\beta }_2\alpha \sigma }{{\beta }_{1}A(\eta +d+\mu )+{\beta }_{2}A\sigma }\hfill \end{array}$$

where $A=\frac{\varsigma \eta \sigma -(\sigma +\mu )(\eta +d+\mu )(\gamma +\mu +\varsigma )}{(\gamma +\mu +\varsigma )(\mu +\psi )\sigma }$.

3.2. Basic Reproduction Number

By considering the fractional system (2), the $R_0$ results from the subsequent matrices.

$$F=\frac{1}{{\lambda }^{n-1}}\left[\begin{array}{ccc}0& \frac{{\beta }_1\alpha }{\mu +\psi }& \frac{{\beta }_2\alpha }{\mu +\psi }\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]$$

$$V=\frac{1}{{\lambda }^{n-1}}\left[\begin{array}{ccc}\sigma +\mu & 0& 0\\ -\sigma & \eta +d+\mu & 0\\ 0& -\eta & \gamma +\mu +\varsigma \end{array}\right]$$

The reproduction number [22] for the system (2) is

$$\begin{array}{c}\hfill R_0=\frac{\alpha \sigma \left[{\beta }_1(\gamma +\mu +\varsigma )+{\beta }_2\eta \right]}{(\mu +\psi )(\sigma +\mu )(\eta +d+\mu )(\gamma +\mu +\varsigma )}\end{array}$$

(3)

4. Stability of the System

This section examines the stability of the system.

Theorem 2.

The equilibrium point $E_0$ is locally asymptotically stable in the system (2).

Proof. 

The reduced matrix of the system (2) is

$$F=\frac{1}{{\lambda }^{n-1}}\left[\begin{array}{cccccc}-\mu -\psi & 0& -\frac{{\beta }_1\alpha }{\mu +\psi }& -\frac{{\beta }_2\alpha }{\mu +\psi }& \varsigma & 0\\ \psi & -\iota -\mu & 0& 0& 0& \varpi \\ 0& 0& \frac{{\beta }_1\alpha }{\mu +\psi }-\sigma -\mu & \frac{{\beta }_2\alpha }{\mu +\psi }& 0& 0\\ 0& 0& \sigma & -\eta -d-\mu & 0& 0\\ 0& 0& 0& \eta & -\gamma -\mu -\varsigma & 0\\ 0& 0& 0& 0& \gamma & -\mu -\varpi \end{array}\right]$$

The eigenvalues of the above matrix are $-(\iota +\mu ),-(\mu +\varpi ),-(\mu +\psi ),-(\gamma +\mu +\varsigma )$ negative. The remaining eigenvalues of the matrix obtained by the equation $c_0{\lambda }^2+c_1\lambda +c_2=0$ where

$$\begin{array}{ccc}\hfill c_0& =& 1,\hfill \\ \hfill c_1& =& \sigma +\eta +d+2\mu -\frac{{\beta }_1\alpha }{\mu +\psi }\hfill \\ \hfill c_2& =& (1-R_0)(\sigma +\mu )(\eta +d+\mu )+\frac{{\beta }_1\alpha }{\mu +\psi }(\sigma -\eta -d-\mu )+\frac{{\beta }_2\alpha \sigma }{\mu +\psi }\frac{\eta -\gamma -\mu -\varsigma }{\gamma +\mu +\varsigma }\hfill \end{array}$$

Hence the model (2) is locally asymptotically stable if $R_0<1$. □

Theorem 3.

The equilibrium point $E_0$ is globally asymptotically stable in the system (2).

Proof. 

Consider the Lyapunov function,

$$\begin{array}{c}\hfill \overline{L}=x_1(S-S_0)+x_2(V-V_0)+x_{3}E+x_{4}I+x_{5}T+x_{6}R\end{array}$$

$$\begin{array}{c}\hfill \frac{d\overline{L}}{dt}=-[\mu (S+V+E+I+T+R)+\iota V+dI-\alpha ]\end{array}$$

Therefore, $\frac{d\overline{L}}{dt}<0$. Hence proved. □

5. A Fractional Approach

The unique solution of the fractional system (2) is investigated and becomes

$$\begin{array}{ccc}\hfill \frac{1}{{\lambda }^{n-1}}{}^{CF}{D}_t^{n}S\left(t\right)& =& G_1(t,S\left(t\right)),\hfill \\ \hfill \frac{1}{{\lambda }^{n-1}}{}^{CF}{D}_t^{n}V\left(t\right)& =& G_2(t,V\left(t\right)),\hfill \\ \hfill \frac{1}{{\lambda }^{n-1}}{}^{CF}{D}_t^{n}E\left(t\right)& =& G_3(t,E\left(t\right)),\hfill \\ \hfill \frac{1}{{\lambda }^{n-1}}{}^{CF}{D}_t^{n}I\left(t\right)& =& G_4(t,I\left(t\right)),\hfill \\ \hfill \frac{1}{{\lambda }^{n-1}}{}^{CF}{D}_t^{n}T\left(t\right)& =& G_5(t,T\left(t\right)),\hfill \\ \hfill \frac{1}{{\lambda }^{n-1}}{}^{CF}{D}_t^{n}R\left(t\right)& =& G_6(t,R\left(t\right)),\hfill \end{array}$$

(4)

Using the fractional derivative theorem, we obtain

$$\begin{array}{ccc}\hfill S\left(t\right)-S\left(0\right)& =& \frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{\int }_0^t{G}_1(\tau ,S){(t-\tau )}^{n-1}d\tau ,\hfill \\ \hfill V\left(t\right)-V\left(0\right)& =& \frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{\int }_0^t{G}_2(\tau ,V){(t-\tau )}^{n-1}d\tau \hfill \\ \hfill E\left(t\right)-E\left(0\right)& =& \frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{\int }_0^t{G}_3(\tau ,E){(t-\tau )}^{n-1}d\tau \hfill \\ \hfill I\left(t\right)-I\left(0\right)& =& \frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{\int }_0^t{G}_4(\tau ,I){(t-\tau )}^{n-1}d\tau \hfill \\ \hfill T\left(t\right)-T\left(0\right)& =& \frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{\int }_0^t{G}_5(\tau ,T){(t-\tau )}^{n-1}d\tau \hfill \\ \hfill R\left(t\right)-R\left(0\right)& =& \frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{\int }_0^t{G}_6(\tau ,R){(t-\tau )}^{n-1}d\tau \hfill \end{array}$$

(5)

Theorem 4.

The kernel $G_1$ satisfies the Lipschitz condition and contraction if the $0\le ({\beta }_1{h}_1+{\beta }_2{h}_2+(\mu +\psi ))<1$ holds.

Proof. 

Let us consider for S and $S_1$,

$$\begin{array}{cc}\hfill \left|\right|G_1(t,S)-G_1(t,S_1)\left|\right|=& \left|\right|-{\beta }_{1}E\left(t\right)(S\left(t\right)-S_1\left(t\right))-{\beta }_{2}I\left(t\right)(S\left(t\right)-S_1\left(t\right))\hfill \\ \hfill & -(\mu +\psi )(S\left(t\right)-S_1\left(t\right))\left|\right|\hfill \end{array}$$

$$\begin{array}{cc}\hfill & \le {\beta }_1\left|\right|E\left(t\right)\left|\right|\left|\right|S\left(t\right)-S_1\left(t\right)\left|\right|+{\beta }_2\left|\right|I\left(t\right)\left|\right|\left|\right|S\left(t\right)-S_1\left(t\right)\left|\right|+(\mu +\psi )\left|\right|S\left(t\right)-S_1\left(t\right)\left|\right|\hfill \\ \hfill & \le \left[{\beta }_1{k}_1+{\beta }_2{k}_2+\mu +\psi \right]\left|\right|S-S_1\left|\right|\hfill \end{array}$$

Hence, $b_1={\beta }_1{k}_1+{\beta }_2{k}_2+\mu +\psi ,$ where $\left|\right|E\left(t\right)\left|\right|=k_1,\left|\right|I\left(t\right)\left|\right|=k_2,$ is bounded function. Hence,

$$\begin{array}{c}\hfill \left|\right|G_1(t,S)-G_1(t,S_1)\left|\right|\le b_1\left|\right|S\left(t\right)-S_1\left(t\right)\left|\right|\end{array}$$

(6)

Therefore, if $0\le {\beta }_1{h}_1+{\beta }_2{h}_2+\mu +\psi <1$ the $G_1$ is a contraction and the Lipschitz condition is obtained. Similarly, we obtain the Lipschitz condition for all the relations.

$$\begin{array}{ccc}\hfill \left|\right|G_2(t,V)-G_2(t,V_1)\left|\right|& \le & b_2\left|\right|V\left(t\right)-V_1\left(t\right)\left|\right|\hfill \\ \hfill \left|\right|G_3(t,E)-G_3(t,E_1)\left|\right|& \le & b_3\left|\right|E\left(t\right)-E_1\left(t\right)\left|\right|\hfill \\ \hfill \left|\right|G_4(t,I)-G_4(t,I_1)\left|\right|& \le & b_4\left|\right|I\left(t\right)-I_1\left(t\right)\left|\right|\hfill \\ \hfill \left|\right|G_5(t,T)-G_5(t,T_1)\left|\right|& \le & b_5\left|\right|T\left(t\right)-T_1\left(t\right)\left|\right|\hfill \\ \hfill \left|\right|G_6(t,R)-G_6(t,R_1)\left|\right|& \le & b_6\left|\right|R\left(t\right)-R_1\left(t\right)\left|\right|\hfill \end{array}$$

where $b_2=(\iota +\mu ),b_3=(\sigma +\mu ),b_4=(\eta +d+\mu ),b_5=(\gamma +\mu +\varsigma ),b_6=(\mu +\varpi )$ are all bounded functions if $0\le b_i<1,$ then $G_i,i=2,3,4,5,6$ are contractions.

Consider the following recursive forms, according to the system (5),

$$\begin{array}{ccc}\hfill H_{1r}\left(t\right)& =& S_r\left(t\right)-S_{r-1}\left(t\right)=\frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{\int }_0^t(G_1(\tau ,S_{r-1})-G_1(\tau ,S_{r-2})){(t-\tau )}^{n-1}d\tau \hfill \\ \hfill H_{2r}\left(t\right)& =& V_r\left(t\right)-V_{r-1}\left(t\right)=\frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{\int }_0^t(G_2(\tau ,V_{r-1})-G_2(\tau ,V_{r-2})){(t-\tau )}^{n-1}d\tau \hfill \\ \hfill H_{3r}\left(t\right)& =& E_r\left(t\right)-E_{r-1}\left(t\right)=\frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{\int }_0^t(G_3(\tau ,E_{r-1})-G_3(\tau ,E_{r-2})){(t-\tau )}^{n-1}d\tau \hfill \\ \hfill H_{4r}\left(t\right)& =& I_r\left(t\right)-I_{r-1}\left(t\right)=\frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{\int }_0^t(G_4(\tau ,I_{r-1})-G_4(\tau ,I_{r-2})){(t-\tau )}^{n-1}d\tau \hfill \\ \hfill H_{5r}\left(t\right)& =& T_r\left(t\right)-T_{r-1}\left(t\right)=\frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{\int }_0^t(G_5(\tau ,T_{r-1})-G_5(\tau ,T_{r-2})){(t-\tau )}^{n-1}d\tau \hfill \\ \hfill H_{6r}\left(t\right)& =& R_r\left(t\right)-R_{r-1}\left(t\right)=\frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{\int }_0^t(G_6(\tau ,R_{r-1})-G_6(\tau ,R_{r-2})){(t-\tau )}^{n-1}d\tau \hfill \end{array}$$

with the initial conditions.

Consider the equation,

$$\begin{array}{cc}\hfill \left|\right|H_{1r}\left(t\right)\left|\right|=& \left|\right|S_r\left(t\right)-S_{r-1}\left(t\right)\left|\right|\hfill \\ \hfill =& \left|\right|\frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{\int }_0^t(G_1(\tau ,S_{r-1})-G_1(\tau ,S_{r-2})){(t-\tau )}^{n-1}d\tau \left|\right|\hfill \\ \hfill \le & \frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{\int }_0^t\left|\right|G_1(\tau ,S_{r-1})-G_1(\tau ,S_{r-2})){(t-\tau )}^{n-1}\left|\right|d\tau \hfill \end{array}$$

with the condition (6).

$$\begin{array}{c}\hfill \left|\right|H_{1r}\left(t\right)\left|\right|\le \frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{b}_1{\int }_0^t\left|\right|H_{1(r-1)}\left(\tau \right)\left|\right|d\tau \end{array}$$

(7)

Similarly, we obtain,

$$\begin{array}{c}\hfill \left|\right|H_{2r}\left(t\right)\left|\right|\le \frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{b}_2{\int }_0^t\left|\right|H_{2(r-1)}\left(\tau \right)\left|\right|d\tau \\ \hfill \left|\right|H_{3r}\left(t\right)\left|\right|\le \frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{b}_3{\int }_0^t\left|\right|H_{3(r-1)}\left(\tau \right)\left|\right|d\tau \\ \hfill \left|\right|H_{4r}\left(t\right)\left|\right|\le \frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{b}_4{\int }_0^t\left|\right|H_{4(r-1)}\left(\tau \right)\left|\right|d\tau \\ \hfill \left|\right|H_{5r}\left(t\right)\left|\right|\le \frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{b}_5{\int }_0^t\left|\right|H_{5(r-1)}\left(\tau \right)\left|\right|d\tau \\ \hfill \left|\right|H_{6r}\left(t\right)\left|\right|\le \frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{b}_6{\int }_0^t\left|\right|H_{6(r-1)}\left(\tau \right)\left|\right|d\tau \end{array}$$

(8)

This can be written as

$$\begin{array}{c}\hfill S_r\left(t\right)=\sum _{j=1}^n{H}_{1j}\left(t\right),V_r\left(t\right)=\sum _{j=1}^n{H}_{2j}\left(t\right),E_r\left(t\right)=\sum _{j=1}^n{H}_{3j}\left(t\right)\\ \hfill I_r\left(t\right)=\sum _{j=1}^n{H}_{4j}\left(t\right),T_r\left(t\right)=\sum _{j=1}^n{H}_{5j}\left(t\right),R_r\left(t\right)=\sum _{j=1}^n{H}_{6j}\left(t\right)\end{array}$$

Next, we must prove the existence of the solution. □

Theorem 5.

The fractional system (2) gives the system of solutions if there exist $t_1$ such that $\frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{t}_1{b}_i<1$i.

Proof. 

Using a recursive technique, (7) and (8) can be written as

$$\begin{array}{c}\hfill \left|\right|H_{1r}\left(t\right)\left|\right|\le \left|\right|S_r\left(0\right)\left|\right|{\left[\frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{b}_{1}t\right]}^r,\left|\right|H_{2r}\left(t\right)\left|\right|\le \left|\right|V_r\left(0\right)\left|\right|{\left[\frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{b}_{2}t\right]}^r\\ \hfill \left|\right|H_{3r}\left(t\right)\left|\right|\le \left|\right|E_r\left(0\right)\left|\right|{\left[\frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{b}_{3}t\right]}^r,\left|\right|H_{4r}\left(t\right)\left|\right|\le \left|\right|I_r\left(0\right)\left|\right|{\left[\frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{b}_{4}t\right]}^r\\ \hfill \left|\right|H_{5r}\left(t\right)\left|\right|\le \left|\right|T_r\left(0\right)\left|\right|{\left[\frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{b}_{5}t\right]}^r,\left|\right|H_{6r}\left(t\right)\left|\right|\le \left|\right|R_r\left(0\right)\left|\right|{\left[\frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{b}_{6}t\right]}^r\end{array}$$

We assume

$$\begin{array}{c}\hfill S\left(t\right)-S\left(0\right)=S_r\left(t\right)-B_{1r}\left(t\right),V\left(t\right)-V\left(0\right)=V_r\left(t\right)-B_{2r}\left(t\right),\\ \hfill E\left(t\right)-E\left(0\right)=E_r\left(t\right)-B_{3r}\left(t\right),I\left(t\right)-I\left(0\right)=I_r\left(t\right)-B_{4r}\left(t\right),\\ \hfill T\left(t\right)-T\left(0\right)=T_r\left(t\right)-B_{5r}\left(t\right),R\left(t\right)-R\left(0\right)=R_r\left(t\right)-B_{6r}\left(t\right),\end{array}$$

where

$$\begin{array}{cc}\hfill \left|\right|B_{1r}\left(t\right)\left|\right|=& \left|\right|\frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{\int }_0^t(G_1(\tau ,S)-G_1(\tau ,S_{r-1}))dd\tau \left|\right|\hfill \\ \hfill \le & \frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{b}_1\left|\right|S-S_{r-1}\left|\right|t\hfill \end{array}$$

(9)

$$\begin{array}{c}\hfill \left|\right|B_{1r}\left(t\right)\left|\right|\le {\left[\frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}t\right]}^{n+1}{b}_1^{n+1}h\end{array}$$

The above equation can be obtained by repeating the method.

$$\begin{array}{c}\hfill \left|\right|B_{1r}\left(t\right)\left|\right|\le {\left[\frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{t}_1\right]}^{n+1}{b}_1^{n+1}h\end{array}$$

We obtain $\left|\right|B_{1r}\left(t\right)\left|\right|\to 0$ as r tends to ∞. Hence, this can be obtained for all functions. This completes the proof. □

Assume the system (2) has another solution $S_1\left(t\right),V_1\left(t\right),E_1\left(t\right),I_1\left(t\right),T_1\left(t\right),R_1\left(t\right)$. We have

$$\begin{array}{c}\hfill S\left(t\right)-S_1\left(t\right)=\frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{\int }_0^t(G_1(\tau ,S)-G_1(\tau ,S_1))d\tau \end{array}$$

we have

$$\begin{array}{c}\hfill \left|\right|S\left(t\right)-S_1\left(t\right)\left|\right|=\frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{\int }_0^t\left|\right|(G_1(\tau ,S)-G_1(\tau ,S_1))\left|\right|d\tau \end{array}$$

$$\begin{array}{c}\hfill \left|\right|S\left(t\right)-S_1\left(t\right)\left|\right|\le \frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)}{b}_{1}t\left|\right|S\left(t\right)-S_1\left(t\right)\left|\right|\end{array}$$

Therefore,

$$\begin{array}{c}\hfill \left|\right|S\left(t\right)-S_1\left(t\right)\left|\right|\left(1-\frac{1}{{\lambda }^{1-n}\Gamma \left(n\right)b_{1}t}\right)\le 0\end{array}$$

(10)

Hence, $\left|\right|S\left(t\right)-S_1\left(t\right)\left|\right|=0$. Therefore, $S\left(t\right)=S_1\left(t\right)$. Similarly, we can obtain this result for all the equations in the model.

6. Numerical Simulation

An approximate method for the solution of the nonlinear problem is called the Adomian decomposition method [23] and provides solutions for both ordinary and fractional differential equations [24]. To obtain results, Laplace Transform is used for model (2)

$$\begin{array}{ccc}\hfill \mathit{L}\left[S\right(t\left)\right]-S\left(0\right)& =& \frac{s+p(1-s)}{s}\mathit{L}[\alpha -{\beta }_{1}SE-{\beta }_{2}SI-(\mu +\psi )S+\varsigma T]\hfill \\ \hfill \mathit{L}\left[V\right(t\left)\right]-V\left(0\right)& =& \frac{s+p(1-s)}{s}\mathit{L}[\psi S-(\iota +\mu )V+\varpi R]\hfill \\ \hfill \mathit{L}\left[E\right(t\left)\right]-E\left(0\right)& =& \frac{s+p(1-s)}{s}\mathit{L}[{\beta }_{1}SE+{\beta }_{2}SI-(\sigma +\mu )E]\hfill \\ \hfill \mathit{L}\left[I\right(t\left)\right]-I\left(0\right)& =& \frac{s+p(1-s)}{s}\mathit{L}[\sigma E-(\eta +d+\mu )I]\hfill \\ \hfill \mathit{L}\left[T\right(t\left)\right]-T\left(0\right)& =& \frac{s+p(1-s)}{s}\mathit{L}[\eta I-(\gamma +\mu +\varsigma )T]\hfill \\ \hfill \mathit{L}\left[R\right(t\left)\right]-R\left(0\right)& =& \frac{s+p(1-s)}{s}\mathit{L}[\gamma T-(\mu +\varpi )R]\hfill \end{array}$$

(11)

The series of the solution is considered to be [25],

$$\begin{array}{ccc}\hfill S\left(t\right)& =& \sum _{k=0}^{\infty }{S}_k\left(t\right),V\left(t\right)=\sum _{k=0}^{\infty }{V}_k\left(t\right),E\left(t\right)=\sum _{k=0}^{\infty }{E}_k\left(t\right)\hfill \\ \hfill I\left(t\right)& =& \sum _{k=0}^{\infty }{I}_k\left(t\right),T\left(t\right)=\sum _{k=0}^{\infty }{T}_k\left(t\right),R\left(t\right)=\sum _{k=0}^{\infty }{R}_k\left(t\right)\hfill \end{array}$$

$$\begin{array}{ccc}\hfill S\left(t\right)E\left(t\right)& =& \sum _{k=0}^{\infty }{A}_k(S,E),S\left(t\right)I\left(t\right)=\sum _{k=0}^{\infty }{B}_k(S,I)\hfill \end{array}$$

where $A_k(S,E),B_k(S,I)$ is used as Adomian polynomial. After some manipulation, the system (11) becomes

$$\begin{array}{cc}\hfill \mathit{L}\left[\sum _{k=0}^{\infty }{S}_k\left(t\right)\right]=S\left(0\right)+& \frac{s+p(1-s)}{s}\mathit{L}\left[\alpha -{\beta }_1\sum _{k=0}^{\infty }{A}_k(S,E)-{\beta }_2\sum _{k=0}^{\infty }{B}_k(S,I)\right.\hfill \\ \hfill & \left.-(\mu +\psi )\sum _{k=0}^{\infty }{S}_k\left(t\right)+\varsigma \sum _{k=0}^{\infty }{T}_k\left(t\right)\right]\hfill \\ \hfill \mathit{L}\left[\sum _{k=0}^{\infty }{V}_k\left(t\right)\right]=V\left(0\right)+& \frac{s+p(1-s)}{s}\mathit{L}\left[\psi \sum _{k=0}^{\infty }{S}_k\left(t\right)-(\iota +\mu )\sum _{k=0}^{\infty }{V}_k\left(t\right)+\varpi \sum _{k=0}^{\infty }{R}_k\left(t\right)\right]\hfill \\ \hfill \mathit{L}\left[\sum _{k=0}^{\infty }{E}_k\left(t\right)\right]=E\left(0\right)+& \frac{s+p(1-s)}{s}\mathit{L}\left[{\beta }_1\sum _{k=0}^{\infty }{A}_k(S,E)+{\beta }_2\sum _{k=0}^{\infty }{B}_k(S,I)\right.\hfill \\ \hfill & \left.(\sigma +\mu )\sum _{k=0}^{\infty }{E}_k\left(t\right)\right]\hfill \\ \hfill \mathit{L}\left[\sum _{k=0}^{\infty }{I}_k\left(t\right)\right]=I\left(0\right)+& \frac{s+p(1-s)}{s}\mathit{L}\left[\sigma \sum _{k=0}^{\infty }{E}_k\left(t\right)-(\eta +d+\mu )\sum _{k=0}^{\infty }{I}_k\left(t\right)\right]\hfill \\ \hfill \mathit{L}\left[\sum _{k=0}^{\infty }{T}_k\left(t\right)\right]=T\left(0\right)+& \frac{s+p(1-s)}{s}\mathit{L}\left[\eta \sum _{k=0}^{\infty }{I}_k\left(t\right)-(\gamma +\mu +\varsigma )\sum _{k=0}^{\infty }{T}_k\left(t\right)\right]\hfill \\ \hfill \mathit{L}\left[\sum _{k=0}^{\infty }{R}_k\left(t\right)\right]=R\left(0\right)+& \frac{s+p(1-s)}{s}\mathit{L}\left[\gamma \sum _{k=0}^{\infty }{T}_k\left(t\right)-(\mu +\varpi )\sum _{k=0}^{\infty }{R}_k\left(t\right)\right]\hfill \end{array}$$

After some manipulation, we obtain

$$\begin{array}{ccc}\hfill \mathit{L}\left[S_{k+1}\left(t\right)\right]& =& \frac{s+p(1-s)}{s}\mathit{L}[\alpha -{\beta }_1{A}_k(S,E)-{\beta }_2{B}_k(S,I)-(\mu +\psi )S_k\left(t\right)+\varsigma T_k\left(t\right)]\hfill \\ \hfill \mathit{L}\left[V_{k+1}\left(t\right)\right]& =& \frac{s+p(1-s)}{s}\mathit{L}[\psi S_k\left(t\right)-(\iota +\mu )V_k\left(t\right)+\varpi R_k\left(t\right)]\hfill \\ \hfill \mathit{L}\left[E_{k+1}\left(t\right)\right]& =& \frac{s+p(1-s)}{s}\mathit{L}[{\beta }_1{A}_k(S,E)+{\beta }_2{B}_k(S,I)-(\sigma +\mu )E_k\left(t\right)]\hfill \\ \hfill \mathit{L}\left[I_{k+1}\left(t\right)\right]& =& \frac{s+p(1-s)}{s}\mathit{L}[\sigma E_k\left(t\right)-(\eta +d+\mu )I_k\left(t\right)]\hfill \\ \hfill \mathit{L}\left[T_{k+1}\left(t\right)\right]& =& \frac{s+p(1-s)}{s}\mathit{L}[\eta I_k\left(t\right)-(\gamma +\mu +\varsigma )T_k\left(t\right)]\hfill \\ \hfill \mathit{L}\left[R_{k+1}\left(t\right)\right]& =& \frac{s+p(1-s)}{s}\mathit{L}[\gamma T_k\left(t\right)-(\mu +\varpi )R_k\left(t\right)]\hfill \end{array}$$

(12)

Exercising the transform on both sides in (12), we have

$$\begin{array}{c}\hfill S_0\left(t\right)=S_0,V_0\left(t\right)=V_0,E_0\left(t\right)=E_0,I_0\left(t\right)=I_0,T_0\left(t\right)=T_0,R_0\left(t\right)=R_0\end{array}$$

$$\begin{array}{ccc}\hfill S_1\left(t\right)& =& [\alpha -{\beta }_1{S}_0\left(t\right)E_0\left(t\right)-{\beta }_2{S}_0\left(t\right)I_0\left(t\right)-(\mu +\psi )S_0\left(t\right)+\varsigma T_0\left(t\right)](1+p(t-1)),\hfill \\ \hfill V_1\left(t\right)& =& [\psi S_0\left(t\right)-(\iota +\mu )V_0\left(t\right)+\varpi R_0\left(t\right)](1+p(t-1)),\hfill \\ \hfill E_1\left(t\right)& =& [{\beta }_1{S}_0\left(t\right)E_0\left(t\right)+{\beta }_2{S}_0\left(t\right)I_0\left(t\right)-(\sigma +\mu )E_0\left(t\right)](1+p(t-1)),\hfill \\ \hfill I_1\left(t\right)& =& [\sigma E_0\left(t\right)-(\eta +d+\mu )I_0\left(t\right)](1+p(t-1)),\hfill \\ \hfill T_1\left(t\right)& =& [\eta I_0\left(t\right)-(\gamma +\mu +\varsigma )T_0\left(t\right)](1+p(t-1)),\hfill \\ \hfill R_1\left(t\right)& =& [\gamma T_0\left(t\right)-(\mu +\varpi )R_0\left(t\right)](1+p(t-1))\hfill \end{array}$$

$$\begin{array}{ccc}\hfill S_2\left(t\right)& =& [\alpha -{\beta }_1{S}_1\left(t\right)E_1\left(t\right)-{\beta }_2{S}_1\left(t\right)I_1\left(t\right)-(\mu +\psi )S_1\left(t\right)+\varsigma T_1\left(t\right)](1+p(t-1)),\hfill \\ \hfill V_2\left(t\right)& =& [\psi S_1\left(t\right)-(\iota +\mu )V_1\left(t\right)+\varpi R_1\left(t\right)](1+p(t-1)),\hfill \\ \hfill E_2\left(t\right)& =& [{\beta }_1{S}_1\left(t\right)E_1\left(t\right)+{\beta }_2{S}_1\left(t\right)I_1\left(t\right)-(\sigma +\mu )E_1\left(t\right)](1+p(t-1)),\hfill \\ \hfill I_2\left(t\right)& =& [\sigma E_1\left(t\right)-(\eta +d+\mu )I_1\left(t\right)](1+p(t-1)),\hfill \\ \hfill T_2\left(t\right)& =& [\eta I_1\left(t\right)-(\gamma +\mu +\varsigma )T_1\left(t\right)](1+p(t-1)),\hfill \\ \hfill R_2\left(t\right)& =& [\gamma T_1\left(t\right)-(\mu +\varpi )R_1\left(t\right)](1+p(t-1))\hfill \end{array}$$

The solution can be expressed as

$$\begin{array}{ccc}\hfill S\left(t\right)& =& \sum _{j=0}^{\infty }{S}_j\left(t\right),V\left(t\right)=\sum _{j=0}^{\infty }{V}_j\left(t\right),E\left(t\right)=\sum _{j=0}^{\infty }{E}_j\left(t\right)\hfill \\ \hfill I\left(t\right)& =& \sum _{j=0}^{\infty }{I}_j\left(t\right),T\left(t\right)=\sum _{j=0}^{\infty }{T}_j\left(t\right),R\left(t\right)=\sum _{j=0}^{\infty }{R}_j\left(t\right)\hfill \end{array}$$

7. Graphical Discussion

The solution of the system (2) is calculated using the

Table 1. Values of the parameters for the fractional model (2).

Parameter	Values	Source
$\alpha$	0.5	Assumed
${\beta }_1$	2.55	[26]
${\beta }_2$	2.28	Assumed
$\varsigma$	0.86	[27]
$\psi$	0.036	[28]
$\iota$	0.006	[29]
$\varpi$	0.002	Assumed
$\sigma$	0.21	[30]
$\mu$	0.53	[31]
$\eta$	0.07	[32]
$\gamma$	0.91	[33]
d	0.012	[34]

.

Figure 1, Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6 show the approximate solutions of various population compartments corresponding to different fractional orders. The density of the susceptible population is decreasing according to the corresponding fractional-order depicted in Figure 1. In Figure 2, the number of the vaccinated population was rapidly increasing. According to Figure 3 and Figure 4, the infected population density increases following exposure to infection. Figure 5 and Figure 6 show that the total number of people in treatment is on the rise, as many people have been cured of the disease through proper treatment, resulting in an increasing number of people in recovery.

8. Conclusions

This paper presents an SVEITR COVID-19 model with a fractional-order derivative of Caputo–Fabrizio. We compute the reproduction number, equilibrium points, and feasible region. Fixed-point theory is applied to show that a unique solution exists for the fractional system. The Adomian method coupled with Laplace integral transform yields the approximate solution of the proposed model. The graphical presentation provides a better understanding of the dynamics, and the proposed technique compares favorably with the other method in terms of its speed of convergence. At the same time, COVID-19 is persistent and largely uncontrollable due to the continual change in the spread of information, a factor which will be explored in more depth in the future by bifurcation analysis, which combines fractional derivatives with optimal control. In future work, this model may be modified by introducing control variables to analyze optimal control strategies.