A Fractional Ordered COVID-19 Model Incorporating Comorbidity and Vaccination

Abstract

The primary goal of this research is to investigate COVID-19 transmission patterns in West Bengal, India in 2021; the first Coronavirus illness (COVID-19) in West Bengal was revealed on 17 March 2020. We employed the modified Susceptible-Asymptomatic-Vaccinated-Comorbidity-Infectious-Recovered (SAVICR) compartmental model as part of fractional orders because of the uncertainty created by the limited Coronavirus (COVID-19) information. In this article, two sub-compartments (Normal Infected and Infected with Co-morbidity) has been considered with vaccinated class, which is relevant in the present situation. We have studied the dynamical analysis of the system and also studied sensitivity of the parameters for West Bengal framework. We have also considered an optimal control problem taking social distancing (non-pharmaceutical treatments) as a control parameter along with vaccination.

1. Introduction

Coronavirus disease (COVID-19) is a condition caused by the newly found Coronavirus SARS-COV-2. The majority of patients infected with COVID-19 will have mild to moderate symptoms and will recover without therapy. When an infected individual coughs, sneezes, or exhales, the virus that causes COVID-19 is primarily transferred by droplets. These droplets are too heavy to float in the air and fall to the ground or other surfaces swiftly. If you are in close contact to someone who has COVID-19, you can be infected by breathing in the virus or touching a contaminated surface and then touching your eyes, nose, or mouth. West Bengal, India was first affected by Coronavirus (COVID-19) on 17 March 2020; in Kolkata a sum of 1,343,442 COVID-19 positive cases have been confirmed by West Bengal’s Health and Family Welfare Department, until 28 May 2021 [1].

Various countries have used non-pharmaceutical treatments (NPIs) such as masking, social distancing, and good hygiene practice to combat the COVID-19 pandemic. This procedure aids in slowing, but not stopping, the progress of the disease. In order to eradicate COVID-19 as a pandemic, effective immunization tactics based on NPIs are required. To overcome COVID-19, a combination of fundamental techniques is unavoidably necessary under the optimal control problem. The actual number of infected patients in West Bengal is anticipated to be higher than the official count due to the restricted number of tests undertaken [1]. According to the WHO, the genuine asymptomatic transmission rates are unknown [2]. As a result, we have taken into account the fact that the asymptomatic class does not transmit sickness and that the Coronavirus is mostly distributed by infected individuals (symptomatically infected and infected with co-morbidity). In phases of co-morbidity (such as heart disease, diabetes, pulmonary disease, and so on) within infected humans, several studies may be available. Co-morbidity or medical disease puts everyone at a higher risk of infection than healthy persons [3]. From statistical databases, it is evident that co-morbidity risk condition in India is pronounced. Patients suffering from hypertension or diabetes mellitus are more likely to acquire a severe course and illness development [4,5,6]. The immunization program has a significant impact on pandemic eradication [7].

Calculus of fractional order can be thought of as an abstraction of the order of differentiation in which the fractional order is substituted by the integer order. However, some properties of classical integer order systems, such as Leibnitz’s rule and chain rules, are not preserved in fractional calculus [8]. An integer system is sometimes unable to express memory-based and hereditary property phenomena [9,10]. There are two main approaches in fractional calculus, namely, the continuous and discrete approaches. The discrete approach is based on Grunwald–Letnikov functional derivative whereas the continuous approaches are based on the Riemann–Liouville, Caputo derivative [11]. It has been observed that the data gathered from real-world events fit better with fractional-order systems. Diethelm [12] has compared the numerical solutions of the fractional-order system and integer-order system, and concluded that the fractional-order system gives more relevant interpretation than integer-order system. Recently, we have contributed to fractional-order dynamical research in the epidemiological field [13,14,15,16]. Many researchers have made significant contributions to various COVID-19 models [14,17,18,19,20,21,22,23,24]. Because the fractional derivative is a generalization of the integer-order derivative, fractional-order modeling has been used to investigate the disease transmission dynamics. In addition, the integer-order differentiation is local, whereas the fractional differentiation is not so. This behavior helps in the simulation of epidemic situations. Furthermore, the fractional derivative has the capability to improve the system’s stability zone. The calculus of fractional order system adds an additional parameter to the modeling framework, which helps in numerical simulations. The prior models are highly useful for analyzing COVID-19 transmission; however, they ignore co-morbidity, vaccinated classes, and the West Bengal pandemic situation. These facts, as well as the benefits of calculus of fraction order, compel us to build the proposed model on COVID-19 in the Caputo fractional framework.

In this study, a modified-fractional-order SAVICR model for two sub-compartments of infected patients (with or without co-morbidity) and vaccinated patients was developed (Section 2). Next, we estimated the model’s fundamental reproduction number. The solutions’ uniqueness, non-negativity, and boundedness were confirmed for the system’s well-posedness (Section 4). We also studied the local stability of disease-free equilibrium point and endemic equilibrium (Section 4, Section 4.4). We considered a control problem by modifying the previous system with controlling vaccination and social distancing (Section 5). Furthermore, we numerically investigated the dynamical system in relation to the parameter values associated with the West Bengal situation in 2021 (Section 6). Section 7 concludes with some key points.

2. Model Formulation

Fractional differentiation operator of the Caputo type was first introduced in the year 1967 by Michele Caputo [25,26].

Definition 1 

(Refs. [25,27]). The Caputo derivative with order $\epsilon \in (0,1)$ for a absolute continuous and differentiable function g on $[0,\infty +)$ is defined as:

$${}_0^c{D}_t^{\epsilon }g\left(t\right)=\left\{\begin{array}{c}{\displaystyle \frac{1}{\Gamma (1-\epsilon )}{\int }_{t_0}^t\frac{g^{{}^{\prime }}\left(s\right)}{{(t-s)}^{{}^{\epsilon }}}ds}\hfill \\ \\ {\displaystyle g^{{}^{\prime }}\left(t\right),\epsilon =1.}\hfill \end{array}\right.$$

In our whole context, we have used ${}_0{D}^{\epsilon }$ instead of ${}_0^c{D}_t^{\epsilon }$ and $0<\epsilon <1$. Caputo derivative is defined only for differentiable functions while the functions that have no first-order derivative might have no Caputo derivative.

Applying Caputo fractional differential equations, the following six compartmental models have been created:

$${\displaystyle \begin{array}{c}{\displaystyle {}_{t_0}^c{D}_t^{\epsilon }S\left(t\right)={\Lambda }^{\epsilon }-{\delta }^{\epsilon }S\left(t\right)-{\alpha }^{\epsilon }S\left(t\right)-{\beta }^{\epsilon }S\left(t\right)I\left(t\right),S\left(0\right)>0,}\hfill \\ \\ {\displaystyle {}_{t_0}^c{D}_t^{\epsilon }A\left(t\right)={\beta }^{\epsilon }S\left(t\right)I\left(t\right)-{\varphi }^{\epsilon }A\left(t\right)+{\sigma }^{\epsilon }{\beta }^{\epsilon }V\left(t\right)I\left(t\right)-{\delta }^{\epsilon }A\left(t\right),A\left(0\right)\ge 0,}\hfill \\ \\ {\displaystyle {}_{t_0}^c{D}_t^{\epsilon }V\left(t\right)={\alpha }^{\epsilon }S\left(t\right)-{\sigma }^{\epsilon }{\beta }^{\epsilon }V\left(t\right)I\left(t\right)-{\delta }^{\epsilon }V\left(t\right),V\left(0\right)>0,}\hfill \\ \\ {\displaystyle {}_{t_0}^c{D}_t^{\epsilon }I\left(t\right)=\xi {\varphi }^{\epsilon }A\left(t\right)-({\delta }_2^{\epsilon }+{\rho }^{\epsilon })I\left(t\right),I\left(0\right)\ge 0,}\hfill \\ \\ {\displaystyle {}_{t_0}^c{D}_t^{\epsilon }C\left(t\right)=(1-\xi ){\varphi }^{\epsilon }A\left(t\right)-({\delta }_1^{\epsilon }+{\gamma }^{\epsilon })C\left(t\right),C\left(0\right)\ge 0,}\hfill \\ \\ {\displaystyle {}_{t_0}^c{D}_t^{\epsilon }R\left(t\right)={\rho }^{\epsilon }I\left(t\right)+{\gamma }^{\epsilon }C\left(t\right)-{\delta }^{\epsilon }R\left(t\right),R\left(0\right)\ge 0.}\hfill \\ \end{array}}$$

(1)

Here ${}_{t_0}^c{D}_t^{\epsilon }$ is the notation of fractional derivative of order $\epsilon \in (0,1)$ in Caputo sense with initial time $t_0\ge 0$ (assuming that $t_0=0$). For simplicity we use ${}_0{D}^{\epsilon }$ instead of ${}_{t_0}^c{D}_t^{\epsilon }$. System 1 is dimensionally correct as both sides have same time dimension $tim{e}^{-\epsilon }$. It is also noticed that the equilibrium points and reproduction number contain $\epsilon$ as all the parameters contain $\epsilon$ in power form. We have omitted the power $\epsilon$ of all parameters (for simplicity) in theoretical analysis. In numerical simulations, we have taken into account the power $\epsilon$ of all parameters. Now, system 1 transforms to:

$${\displaystyle \begin{array}{c}{\displaystyle {}_0{D}^{\epsilon }S\left(t\right)=\Lambda -\delta S\left(t\right)-\alpha S\left(t\right)-\beta S\left(t\right)I\left(t\right),S\left(0\right)>0,}\hfill \\ \\ {}_0{D}^{\epsilon }A\left(t\right)=\beta S\left(t\right)I\left(t\right)-\varphi A\left(t\right)+\sigma \beta V\left(t\right)I\left(t\right)-\delta A\left(t\right),A\left(0\right)\ge 0,\hfill \\ \\ {}_0{D}^{\epsilon }V\left(t\right)=\alpha S\left(t\right)-\sigma \beta V\left(t\right)I\left(t\right)-\delta V\left(t\right),V\left(0\right)>0,\hfill \\ \\ {}_0{D}^{\epsilon }I(t)=\xi \varphi A(t)-\left({\delta }_2+\rho \right)I(t),I(0)\ge 0\hfill \\ \\ {}_0{D}^{\epsilon }C(t)=(1-\xi )\varphi A(t)-\left({\delta }_1+\gamma \right)C(t),C(0)\ge 0\hfill \\ \\ {}_0{D}^{\epsilon }R(t)=\rho I(t)+\gamma C(t)-\delta R(t),R(0)\ge 0\hfill \end{array}}$$

(2)

It is presumed $\delta \le {\delta }_2\le {\delta }_1$. Schematic diagram of system 2 is shown in Figure 1.

The susceptible, asymptomatically infected, Vaccinated class, symptomatically infected without co-morbidity, symptomatically infected with co-morbidity and recovered or removed population at time t are represented by $S\left(t\right),A\left(t\right),V\left(t\right),I\left(t\right),C\left(t\right)$, and $R\left(t\right)$, respectively.

Table 1. Description of parameters used in the system 2.

$\Lambda$	Recruitment rate of S
$\sigma$	Transmission rate at which vaccinated people become asymptomatic
$\beta$	Transmission rate at which susceptible people become asymptomatic
$\delta$	Morbidity rate
${\delta }_1$	Disease-induced death rate of symptomatically infected with co-morbidity class
${\delta }_2$	Disease induced death rate of symptomatically infected without co-morbidity class
$\rho$	Recovery rate of symptomatically infected without co-morbidity class $\left(I\right)$
$\varphi$	Conversion rate of asymptomatic infected people to symptomatic
$\xi$	Fraction of asymptomatic population enters to infected class without co-morbidity $\left(I\right)$
$\gamma$	Recovery rate of symptomatically infected with co-morbidity class (C)
$\alpha$	Rate of vaccination

provides a brief description of all parameters.

3. Equilibrium Points and Basic Reproduction Number

The disease-free equilibrium point $E_0$ and endemic equilibrium point $E_1$ of system (2) are given below.

• $E_0=\left({\displaystyle \frac{\Lambda }{\delta +\alpha },0,\frac{\alpha \Lambda }{\delta (\delta +\alpha )},0,0,0}\right)$
• $E_1=(S^{*},A^{*},V^{*},I^{*},C^{*},R^{*}).$

Here

$$\begin{array}{c}{\displaystyle S^{*}=\frac{\Lambda }{\alpha +\delta +\beta I^{*}}}\hfill \\ \\ {\displaystyle A^{*}=\frac{{\delta }_2+\rho }{\xi \varphi }{I}^{*}}\hfill \\ \\ {\displaystyle V^{*}=\frac{\alpha \Lambda }{[\alpha +\delta +\beta I^{*}][\sigma \beta I^{*}+\delta ]}}\hfill \\ \\ {\displaystyle I^{*}=\frac{\Lambda \xi \varphi }{(\varphi +\delta )({\delta }_2+\rho )}-\frac{\delta }{\sigma \beta }}\hfill \\ \\ {\displaystyle C^{*}=\frac{(1-\xi )({\delta }_2+\rho )}{\xi (\gamma +{\delta }_1)}{I}^{*}}\hfill \\ \\ {\displaystyle R^{*}=\left[{\displaystyle \frac{\rho }{\delta }+\frac{\gamma (1-\xi )({\delta }_2+\rho )}{\delta \xi ({\delta }_1+\gamma )}}\right]I^{*}}\hfill \end{array}$$

(3)

$E_1$ remains in ${\mathbb{R}}_{+}^6$ if

$${\displaystyle \Lambda \xi \varphi \sigma \beta \ge \delta (\varphi +\delta )({\delta }_2+\rho )}$$

The basic reproduction number, indicated by $R_0$, is the estimated number of secondary cases generated by the infection of a single individual. The basic reproduction number $R_0$ is calculated at disease-free equilibrium point $E_0$ by the help of the NGM (Next Generation Matrix) whose dominant eigenvalue is $R_0$ [28,29]. The components A, I, and C in system (2) are explicitly appended with a disease transmission. The matrices $\overline{F},\overline{V}$ represent respective new infection and transition matrices.

Consider $u={(A,I,C)}^T$, the subsystem of system (2) can be represented as:

$${\displaystyle {}_0{D}_t^{\epsilon }u=\overline{F}\left(u\right)-\overline{V}\left(u\right),}$$

$$\overline{F}\left(u\right)=\left[\begin{array}{c}\beta SI+\sigma \beta VI\\ \\ 0\\ \\ 0\\ \end{array}\right]$$

and

$$\overline{V}\left(u\right)=\left[\begin{array}{c}\varphi A+\delta A\\ \\ ({\delta }_2+\rho )I-\xi \varphi A\\ \\ ({\delta }_1+\gamma )C-(1-\xi )\varphi A\\ \end{array}\right]$$

$$D\overline{F}\left(E_0\right)=F$$

$$D\overline{V}\left(E_0\right)=V$$

$D\overline{F}\left(u\right)$ and $D\overline{V}\left(u\right)$ are the Jacobian of $\overline{F},\overline{V}$ at disease-free equilibrium $E_0$ where

$$F=\left[\begin{array}{ccccc}{\displaystyle 0}& & {\displaystyle \frac{\Lambda \beta }{\delta +\alpha }+\sigma \beta \frac{\alpha \Lambda }{\delta (\delta +\alpha )}}& & 0\\ \\ {\displaystyle 0}& & 0& & 0\\ \\ {\displaystyle 0}& & 0& & 0\\ \end{array}\right]$$

$$V=\left[\begin{array}{ccccc}{\displaystyle (\delta +\varphi )}& & 0& & 0\\ \\ {\displaystyle -\xi \varphi }& & {\displaystyle {\delta }_2+\rho }& & 0\\ \\ {\displaystyle -(1-\xi )\varphi }& & 0& & {\displaystyle {\delta }_1+\gamma }\end{array}\right]$$

Thus, we obtain

$${\displaystyle R_0=\frac{\xi \beta \Lambda \varphi }{(\varphi +\delta )(\delta +\alpha )({\delta }_2+\rho )}+\frac{\alpha \xi \varphi \beta \sigma \Lambda }{\delta (\varphi +\delta )(\delta +\alpha )({\delta }_2+\rho )}.}$$

(4)

Reproduction number $R_0$ is a function of $\epsilon$ because each parameter is a function of $\epsilon$. The value of $\epsilon$ has been fixed for analysis reasons. If we modify the value of $\epsilon$, it will affect all other parametric values, including the value of $R_0$.

4. Basic Analysis of the System (2)

This section has carefully accomplished certain basic results such as existence, non-negativity, boundedness, and stability of system (2).

4.1. Preliminaries of Caputo Fractional Calculus

The following definitions and theorems are essential for further theoretical study.

Lemma 1 

(Ref. [30]). Consider a continuous function on $[a,b]$ where $a>0$ and ${}_0{D}^{\epsilon }\psi \left(t\right)$ is continuous on $\left(a,b\right]$, then

$$\psi \left(x\right)=\psi \left(a\right)+\frac{1}{\Gamma \left(\epsilon \right)}{(x-a)}^{\epsilon }.{}_0{D}^{\epsilon }\psi \left(\eta \right),$$

where $0<a\le \eta \le x$, $\forall x\in \left(a,b\right]$.

Remark 1. 

If ${}_0{D}^{\epsilon }\psi \left(t\right)\ge 0\left({}_0{D}^{\epsilon }\psi \left(t\right)\le 0\right)$ for all $t\in (a,b)$, then $\psi \left(t\right)$ is a non-decreasing (non-increasing) function for $t\in \left[a,b\right]$.

Definition 2 

(Ref. [8]). Mittag–Leffler functions of single and double parameters are described below:

$$E_{\epsilon }\left(w\right)=\sum _{k=0}^{\infty }\frac{w^k}{\Gamma (\epsilon k+1)}\mathrm{and}{E}_{{\epsilon }_1,{\epsilon }_2}\left(w\right)=\sum _{k=0}^{\infty }\frac{w^k}{\Gamma ({\epsilon }_{1}k+{\epsilon }_2)},\mathrm{where}\epsilon ,{\epsilon }_1,{\epsilon }_2\in {\mathbb{R}}_{+}.$$

Theorem 1 

(Ref. [31]). Suppose $g\left(t\right)$ is differentiable with exponential order and $0<\epsilon <1$. Further, if ${}_0{D}^{\epsilon }g\left(t\right)$ is piece-wise continuous on $\left[0,\infty \right)$, then

$$\mathcal{L}\left\{{}_0{D}^{\epsilon }g\left(t\right)\right\}=s^{\epsilon }F\left(s\right)-\sum _{j=0}^{n-1}{s}^{\epsilon -j-1}{g}^j\left(0\right),$$

where the Laplace transform of $g\left(t\right)$ is represented by $F\left(s\right)=\mathcal{L}\left\{g\left(t\right)\right\}$.

Theorem 2 

(Ref. [32]). For any complex number M and ${\epsilon }_1,{\epsilon }_2\in (0,1)$, we have

$$\mathcal{L}\left\{t^{{\epsilon }_2-1}{E}_{{\epsilon }_1,{\epsilon }_2}\left(M{t}^{{\epsilon }_1}\right)\right\}=\frac{s^{{\epsilon }_1-{\epsilon }_2}}{(s^{{\epsilon }_1}-M)},$$

where the real part of $s>{∥M∥}^{\frac{1}{{\epsilon }_1}}$ and $E_{{\epsilon }_1,{\epsilon }_2}$ is the two parametric Mittag–Leffler function.

Theorem 3 

(Ref. [27]). Consider:

$${}_0{D}_t^{\epsilon }x\left(t\right)=\Phi \left(x\right),$$

with $0<\epsilon <1$, $x\in {\mathbb{R}}^n)$. The equilibrium points of the above system are solutions to the equation $\Phi \left(x\right)=0$. An equilibrium is locally asymptotically stable if all eigenvalues (${\lambda }_i$) of the Jacobian matrix ${\displaystyle J=\frac{\partial \Phi }{\partial x}}$ calculated at the equilibrium points satisfy ${\displaystyle \left|arg\left({\lambda }_i\right)\right|>\frac{\epsilon \pi }{2}}$.

4.2. Existence and Uniqueness

Lemma 2 

(Ref. [33]). Consider the system:

$${}_0{D}^{\epsilon }x\left(t\right)=g(t,x)$$

(5)

with initial condition $x\left(0\right)=x_0$, where $\epsilon \in \left(0,1\right)$, $g:\left[0,\infty \right)\times \Omega \to {\mathrm{I}\mathrm{R}}^n,\Omega \subseteq {\mathrm{I}\mathrm{R}}^n$, if $g(t,x)$ satisfies local Lipschitz condition with respect to x, then there exists a unique solution of (5) on $\left[0,\infty \right)\times \Omega$.

In our context $\Omega =\{(S,A,V,I,C,R)\in {\mathbb{R}}^6:max(\left|S\right|,\left|A\right|,\left|V\right|,\left|I\right|,\left|C\right|,\left|R\right|)\le M<\infty \}$.

Now, we have the following theorem.

Theorem 4.

The system (2) with initial condition $X\left(0\right)=\left(S\right(0),A(0),V(0),I(0),C(0),R(0\left)\right)\in \Omega$, possess a unique solution $X\left(t\right)\in \Omega$, $\forall t\ge 0$.

Proof of this Theorem is given in Appendix A.

4.3. Non-Negativity and Boundedness

Here, we first studied the non-negativity and then we demonstrated that all solutions of system (2) are bounded.

Theorem 5.

The solutions of $SAVICR$ model (2) with initial values lie in ${\Gamma }_{+}={\mathbb{R}}_{+}^6$.

Proof. 

$${\displaystyle {}_0{D}^{\epsilon }S\left(t\right){|}_{S\left(t\right)=0}=\Lambda >0}$$

(6)

$${\displaystyle {}_0{D}^{\epsilon }A\left(t\right){|}_{A\left(t\right)=0}=\beta SI+\sigma \beta VI}$$

(7)

$${\displaystyle {}_0{D}^{\epsilon }V\left(t\right){|}_{V\left(t\right)=0}=\alpha S}$$

(8)

$${\displaystyle {}_0{D}^{\epsilon }I\left(t\right){|}_{I\left(t\right)=0}=\xi \varphi A}$$

(9)

$${\displaystyle {}_0{D}^{\epsilon }C\left(t\right){|}_{C\left(t\right)=0}=(1-\xi )\varphi A}$$

(10)

$${\displaystyle {}_0{D}^{\epsilon }R\left(t\right){|}_{R\left(t\right)=0}=\rho I+\gamma C}$$

(11)

From (6), we have

$${}_0{D}^{\epsilon }{S\left(t\right)|}_{S\left(t\right)=0}=\Lambda >0.$$

We can deduce from Lemma 1 that $S\left(t\right)$ is increasing in near $S\left(t\right)=0$, and $S\left(t\right)$ is unable to go across the axis $S\left(t\right)=0$. As a result, S is non-negative and similarly we can show $V\left(t\right)>0$ for all $t\ge 0$ (from Equation (8)). We now assert that $A\left(t\right)$ commences in ${\Gamma }_{+}$ cannot be negative. Contrarily, there exists ${\tau }_1$ such that $A\left(t\right)$ traverses the $A\left(t\right)=0$ axis for the first time at $t={\tau }_1$ with the following possibilities.

$$\left\{\begin{array}{c}{\displaystyle A\left(t\right)>0,\mathrm{for}0\le t<{\tau }_1,}\hfill \\ A\left({\tau }_1\right)=0,\hfill \\ A\left({\tau }_1^{+}\right)<0.\hfill \end{array}\right.$$

Relation (7) revels that ${\displaystyle {}_0{D}^{\epsilon }A\left(t\right){|}_{A\left({\tau }_1\right)=0}=\beta S\left({\tau }_1\right)I\left({\tau }_1\right)+\sigma \beta V\left({\tau }_1\right)I\left({\tau }_1\right)}$. and this leads the following cases.

• Case1: Suppose $\beta S\left({\tau }_1\right)I\left({\tau }_1\right)+\sigma \beta V\left({\tau }_1\right)I\left({\tau }_1\right)\ge 0$, then from Lemma 1, $A\left(t\right)$ is not decreasing in the neighborhood of $t={\tau }_1$, and $A\left({\tau }_1^{+}\right)=0$. As a result, a contradiction occurs.
• Case2: On the other hand $\beta S\left({\tau }_1\right)I\left({\tau }_1\right)+\sigma \beta V\left({\tau }_1\right)I\left({\tau }_1\right)<0$, which discloses $I\left({\tau }_1\right)$ must be negative. If $I\left({\tau }_1\right)<0$, then a ${\tau }_2$$(0<{\tau }_2<{\tau }_1)$ exists with

  $$\left\{\begin{array}{c}{\displaystyle I\left(t\right)>0,\mathrm{for}0\le t<{\tau }_2,}\hfill \\ I\left({\tau }_2\right)=0,\hfill \\ I\left({\tau }_2^{+}\right)<0.\hfill \end{array}\right.$$

From (9), we obtain

$${\displaystyle {}_0{D}_t^{\epsilon }I\left(t\right){|}_{I\left({\tau }_2\right)=0}=\xi \varphi A\left({\tau }_2\right)>0,}$$

which refutes $I\left({\tau }_2^{+}\right)<0$. Therefore, $I\left(t\right)\ge 0$, $\forall t\in [0,\infty )$ and also $A\left(t\right)\ge 0$, $\forall t\in [0,\infty )$.

Once again from (10) and (11), we can conclude $C\left(t\right),R\left(t\right)$ are non-negative.

As a result, all solutions of system (2) that begin with ${\Gamma }_{+}$ are constricted to the region ${\Gamma }_{+}$. □

Theorem 6 

(Boundedness). System (1) has uniformly bounded solutions $X\left(t\right)=(S,A,V,I,C,R)$.

Proof. 

It is clear from the first equation of (2) that

$${\displaystyle {}_0{D}^{\epsilon }S\left(t\right)\le \Lambda -\delta S}$$

Using Laplace transforms $\left(\mathcal{L}\left\{·\right\}\right)$, we have

$$\begin{array}{c}{\displaystyle s^{\epsilon }\mathcal{L}\left\{S\left(t\right)\right\}-s^{\epsilon -1}S\left(0\right)+\delta \mathcal{L}\left\{S\left(t\right)\right\}\le \frac{\Lambda }{s},}\hfill \\ \\ {\displaystyle \Rightarrow \mathcal{L}\left\{S\left(t\right)\right\}\le A\frac{s^{\epsilon -(1+\epsilon )}}{s^{\epsilon }+\delta }+S\left(0\right)\frac{s^{\epsilon -1}}{s^{\epsilon }+\delta }}\hfill \end{array}$$

Using inverse Laplace transforms (using Theorem 2):

$${\displaystyle S\left(t\right)\le S\left(0\right)E_{\epsilon ,1}(-\delta t^{\epsilon })+\Lambda t^{\epsilon }{E}_{\epsilon ,\epsilon +1}(-\delta t^{\epsilon })}$$

Therefore, ${\displaystyle S\left(t\right)\le M_1[E_{\epsilon ,1}(-\delta t^{\epsilon })+\delta t^{\epsilon }{E}_{\epsilon ,\epsilon +1}(-\delta t^{\epsilon })]=\frac{M_1}{\Gamma \left(1\right)}=M_1,M_1=max\left\{\frac{\Lambda }{\delta },S\left(0\right)\right\}.}$

From the properties of Mittag–Leffler function: [34]:

$${\displaystyle E_{{\alpha }_1,{\alpha }_2}\left(z\right)=z{E}_{{\alpha }_1,{\alpha }_1+{\alpha }_2}\left(z\right)+\frac{1}{\Gamma \left({\alpha }_2\right)}}$$

In this circumstance,

$${\displaystyle E_{\epsilon ,1}(-\delta t^{\epsilon })=(-\delta t^{\epsilon })E_{\epsilon ,\epsilon +1}(-\delta t^{\epsilon })+\frac{1}{\Gamma \left(1\right)}}$$

(12)

Let $N\left(t\right)=S\left(t\right)+A\left(t\right)+V\left(t\right)+I\left(t\right)+C\left(t\right)+R\left(t\right)$.

Now,

$$\begin{array}{ccc}{\displaystyle {}_0{D}_t^{\epsilon }N\left(t\right)}\hfill & & {\displaystyle {=}_0{D}_t^{\epsilon }S\left(t\right){+}_0{D}_t^{\epsilon }A\left(t\right){+}_0{D}_t^{\epsilon }V\left(t\right){+}_0{D}_t^{\epsilon }I\left(t\right){+}_0{D}_t^{\epsilon }C\left(t\right){+}_0{D}_t^{\epsilon }R\left(t\right)}\hfill \\ \\ {\displaystyle }\hfill & & =\Lambda -\{\delta S\left(t\right)+\delta A\left(t\right)+\delta V\left(t\right)+{\delta }_{2}I\left(t\right)+{\delta }_{1}C\left(t\right)+\delta R\left(t\right)\}\hfill \\ \\ {\displaystyle }\hfill & & \le \Lambda -{\delta }_{m}N\left(t\right),\mathrm{where}{\delta }_m=min\{{\delta }_1,{\delta }_2,\delta \}\hfill \end{array}$$

Hence,

$${\displaystyle {}_0{D}_t^{\epsilon }N\left(t\right)+{\delta }_{m}N\left(t\right)\le \Lambda }$$

Using Laplace transformation, we have (using Theorem 1):

$$\begin{array}{c}{\displaystyle s^{\epsilon }F\left(s\right)-s^{\epsilon -1}N\left(0\right)+{\delta }_{m}F\left(s\right)\le \frac{\Lambda }{s},\mathrm{where}F\left(s\right)=\mathcal{L}\left\{N\left(t\right)\right\}}\hfill \\ \\ {\displaystyle \Rightarrow F\left(s\right)\le \Lambda \frac{s^{-1}}{s^{\epsilon }+{\delta }_m}+\frac{N\left(0\right)s^{\epsilon -1}}{s^{\epsilon }+{\delta }_m}=\frac{s^{\epsilon -1}N\left(0\right)}{s^{\epsilon }+{\delta }_m}+\frac{\Lambda s^{\epsilon -(1+\epsilon )}}{s^{\epsilon }+{\delta }_m}}\hfill \end{array}$$

Using inverse Laplace transforms (using Theorem 2):

$${\displaystyle N\left(t\right)=N\left(0\right)E_{\epsilon ,1}(-{\delta }_m{t}^{\epsilon })+\Lambda t^{\epsilon }{E}_{\epsilon ,\epsilon +1}(-{\delta }_m{t}^{\epsilon })}$$

(13)

From (12) and (13), we obtain

$${\displaystyle N\left(t\right)\le M_2[E_{\epsilon ,1}(-{\delta }_m{t}^{\epsilon })+{\delta }_m{t}^{\epsilon }{E}_{\epsilon ,\epsilon +1}(-{\delta }_m{t}^{\epsilon })]=\frac{M_2}{\Gamma \left(1\right)}=M_2,M_2=max\left\{\frac{\Lambda }{{\delta }_m},N\left(0\right)\right\}}$$

Thus $S\left(t\right),N\left(t\right)$ are bounded and hence the solution $\left(S\right(t),A(t),V(t),I(t),C(t),R(t\left)\right)$ is bounded in $\left\{(S,A,V,I,C,R)\right|S+A+V+I+C+R\le M_2;S\le M_1\}$ for $t\in [0,\infty )$. □

4.4. Local Stability

For simplicity of analytical studies, we have reduced the system (2) by omitting the sixth equation as R does not present in top five equations of system (2). The dynamics of R also be obtained from the dynamics of $S,A,V,I,C$. The following is the reduced system:

$$\begin{array}{c}{\displaystyle {}_0{D}^{\epsilon }S\left(t\right)=\Lambda -\delta S\left(t\right)-\alpha S\left(t\right)-\beta S\left(t\right)I\left(t\right),S\left(0\right)>0,}\hfill \\ \\ {}_0{D}^{\epsilon }A\left(t\right)=\beta S\left(t\right)I\left(t\right)-\varphi A\left(t\right)+\sigma \beta V\left(t\right)I\left(t\right)-\delta A\left(t\right),A\left(0\right)\ge 0,\hfill \\ \\ {}_0{D}^{\epsilon }V\left(t\right)=\alpha S\left(t\right)-\sigma \beta V\left(t\right)I\left(t\right)-\delta V\left(t\right),V\left(0\right)>0,\hfill \\ \\ {}_0{D}^{\epsilon }I(t)=\xi \varphi A(t)-\left({\delta }_2+\rho \right)I(t),I(0)\ge 0\hfill \\ \\ {}_0{D}^{\epsilon }C(t)=(1-\xi )\varphi A(t)-\left({\delta }_1+\gamma \right)C(t),C(0)\ge 0\hfill \end{array}$$

(14)

In this casewe have to analyze the stability of equilibriumpoints $E_0^{*}=\left({\displaystyle \frac{\Lambda }{\delta +\alpha },0,\frac{\alpha \Lambda }{\delta (\delta +\alpha )},0,0}\right)$ and $E_1^{*}=\left(S^{*},A^{*},V^{*},I^{*},C^{*}\right)$, where

$$\begin{array}{c}{\displaystyle S^{*}=\frac{\Lambda }{\alpha +\delta +\beta I^{*}}}\hfill \\ \\ {\displaystyle A^{*}=\frac{{\delta }_2+\rho }{\xi \varphi }{I}^{*}}\hfill \\ \\ {\displaystyle V^{*}=\frac{\alpha \Lambda }{[\alpha +\delta +\beta I^{*}][\sigma \beta I^{*}+\delta ]}}\hfill \\ \\ {\displaystyle I^{*}=\frac{\Lambda \xi \varphi }{(\varphi +\delta )({\delta }_2+\rho )}-\frac{\delta }{\sigma \beta }}\hfill \\ \\ {\displaystyle =\frac{\delta }{\sigma \beta }\left({\displaystyle \frac{\Lambda \xi \varphi \sigma \beta }{\delta (\varphi +\delta )({\delta }_2+\rho )}-1}\right)}\hfill \\ \\ {\displaystyle =\left(\frac{R_0}{K_0}-1\right)\frac{\delta }{\sigma \beta }}\hfill \\ \\ {\displaystyle C^{*}=\frac{(1-\xi )({\delta }_2+\rho )}{\xi (\gamma +{\delta }_1)}{I}^{*}}\hfill \\ \end{array}$$

and $R_0=\frac{\xi \beta \mathrm{\Lambda }\varphi }{(\varphi +\delta )(\delta +\alpha )\left({\delta }_2+\rho \right)}+\frac{\alpha \xi \varphi \beta \sigma \mathrm{\Lambda }}{\delta (\varphi +\delta )(\delta +\alpha )\left({\delta }_2+\rho \right)},K_0=\left(\frac{\delta }{(\delta +\alpha )\sigma }+\frac{\alpha }{\alpha +\delta }\right)$

For $E_1^{*}$ to exist in feasible region ${\mathbb{R}}_{+}^5$, it is necessary and sufficient that $R_0>K_0$. From simple calculation it is clear that $K_0\ge 1$ if $\sigma \le 1$. If $R_0=K_0$, then the disease-free equilibrium and endemic equilibrium points will be confluent.

The following lemma is essential for analysis of system (14).

Lemma 3 

(Refs. [35,36]). Suppose $P\left(\lambda \right)={\lambda }^n+a_1{\lambda }^{n-1}+a_2{\lambda }^{n-2}+...+a_n$ be the characteristic equation of system $\left(14\right)$ evaluated at equilibrium point. For $n=4$, $\nabla \left(P\right)$ is the discriminant of the characteristic equation $P\left(\lambda \right)={\lambda }^4+a_1{\lambda }^3+a_2{\lambda }^2+a_3\lambda +a_4$ which is represented as $\nabla \left(P\right)=$$\left|\begin{array}{ccccccc}1& a_1& a_2& a_3& a_4& 0& 0\\ 0& 1& a_1& a_2& a_3& a_4& 0\\ 0& 0& 1& a_1& a_2& a_3& a_4\\ 4& 3a_1& 2a_2& a_3& 0& 0& 0\\ 0& 4& 3a_1& 2a_2& a_3& 0& 0\\ 0& 0& 4& 3a_1& 2a_2& a_3& 0\\ 0& 0& 0& 4& 3a_1& 2a_2& a_3\end{array}\right|$

Lemma 4 

(Ref. [37]). Suppose $\nabla \left(P\right)<0,a_1>0,a_2>0,a_3>0,a_4>0$ hold, then the equilibrium point is asymptotically stable if ${\displaystyle a_2=\frac{a_1{a}_4}{a_3}+\frac{a_3}{a_1}}$. If ${\displaystyle a_2\ne \frac{a_1{a}_4}{a_3}+\frac{a_3}{a_1}}$ then ${\displaystyle \epsilon <\frac{1}{3}}$ must hold.

Theorem 7.

The disease-free equilibrium point $E_0^{*}$ of system (14) is locally asymptotically stable if and only if $R_0<1$.

Proof. 

Proof is given in Appendix A. □

Theorem 8.

The endemic equilibrium $E_1$ of system (14) is unstable.

Proof. 

Proof is given in Appendix A. □

4.5. Sensitivity Analysis

It is clear that $R_0$ relies on different parameters and value of $R_0$ is $0.1146$ for

Table 2. The values of the parameters used in system (1) corresponding to the situation in West Bengal (India).

Parameters	Values	Source
$\Lambda$	$0.00004$	[1]
$\beta$	$0.3945$	Estimated
$\delta$	$1/(70.4\times 365)$	[1]
${\delta }_1$	$1/(60\times 365)$	Assumed
${\delta }_2$	$1/(65\times 365)$	Assumed
$\varphi$	$1/5.2$	[39]
$\xi$	$0.3$	Assumed
$\rho$	$0.1245$	Assumed
$\sigma$	$0.08$	Estimated
$\gamma$	$0.1240$	Assumed
$\alpha$	$3\times {10}^{-4}$	Estimated
$\epsilon$	0.7	Estimated

. The normalized forward sensitivity index with regard to each parameter has been derived as follows to analyze the sensitivity of $R_0$ to any parameter (say, $\theta$) [29,38]:

$${\displaystyle {\Omega }_{\theta }^{R_0}=\frac{\partial R_0}{\partial \theta }\frac{\theta }{R_0}}$$

The sensitivity index can be affected by different factors, but it may also be constant or unrelated to them. These values are critical for estimating the sensitivity of parameters, which should be performed with caution because even minor changes in a parameter can produce significant quantitative changes. It is not necessary to be cautious while estimating a parameter with a lower sensitivity index because little changes in that parameter create small changes. We focused on the parameters $\Lambda ,\alpha ,\beta ,\varphi ,\rho ,and\sigma$ because we could not manage $\delta ,\xi$.

Table 3. Normalized forward sensitivity index for different parameters correspondence to Table 2.

Parameters	Sensitivity Index
$\beta$	+1
$\varphi$	0.0104
$\rho$	−0.9964
$\alpha$	−0.8704
$\Lambda$	+1
$\sigma$	+0.0152

shows the sensitivity index values for the parameters $\Lambda ,\alpha ,\beta ,\varphi ,\rho ,\sigma$ that match Table 2. From Table 3, it is revealed that infection transmission rate $\beta$, recovery rate of without co-morbidity class, and vaccination rate are highly sensitive. According to the sensitivity index (Table 3), the transmission rate from symptomatically infected people is the most sensitive of all the criteria for lowering disease prevalence, and vaccination lowers the number of symptomatically infected people. As a result, reducing virus transmission through social distancing and vaccination, as well as other preventative measures, would aid in the timely management of this pandemic situation.

5. Fractional-Order Optimal Control

One of the most important preventative precautions that everyone should do is to keep social distances and maintain adequate cleanliness [40]. Aside from that, other vaccines are now accessible, and people are being advised to take the appropriate amount to avoid future infection. As a result, these tactics have been implemented into this system in order to limit the quick transfer of information. Our goal is to reduce the number of susceptible and contaminated people together with the cost of applying combined control measures $u,v$ (social distance and vaccination). Though Ding et al. [41] and Agarwal et al. [42] have made early contributions to optimal control theory in fractional calculus, recently many researchers have significantly contributed in the field of fractional-order control [43,44].

One of the most important results in classical and fractional optimal control is the Pontryagin Principle [43,45]. The method is similar to what was used to solve the traditional integer-order optimal control problem. We have already studied some control problems using the Pontryagin Principle for fractional-order systems [15,16]. The main distinction is that in a fractional-order optimal control problem, the adjoint equations are in the Right-Riemann–Liouville derivative (${}_t^{RL}{D}_T^{\epsilon }$) of order $\epsilon$ but the co-state equations are in Caputo differential equations. System (2) is reintroduced by implementing some control interventions that can reduce the disease load. We have introduced two controls, $u,v$, where v is considered as intensity of vaccination and u can be considered as non-pharmaceutical interventions (complete or partial lockdown, social distancing, etc.) with $0\le u,v\le 1$. Here, $u=0$ portrays no such non-pharmaceutical interventions maintained in definite time zone ([0,$T_f$]), while $u=1$ represents the scenario where full $(100\%)$ interventions is maintained; whereas $v=0$ denotes no vaccination and $v=1$ stands for fully vaccinated.

Consider $U=\left\{(u,v)\in L_1^{\infty }(0,T_f),0\le u,v\le 1,t\in (0,T_f)\right\}$ be the control space with measurable control function $u,v$. Our main objective is to minimize the cost function J in $[0,T_f]$ by finding optimal solution as follows:

$$J(u^{*},v^{*})=J_{min\left(u\right(t),v(t\left)\right)\in U}$$

(15)

$${\displaystyle J(u,v)={\int }_0^{T_f}\left[{\displaystyle m_{1}I\left(t\right)+m_{2}C\left(t\right)+\frac{m_3}{2}{u}^2\left(t\right)+\frac{m_4}{2}{v}^2\left(t\right)}\right]dt}$$

(16)

subject to

$$\begin{array}{c}{\displaystyle {}_0{D}^{\epsilon }S\left(t\right)=\Lambda -(\delta +\alpha v)S-\beta SI(1-u),S\left(0\right)>0,}\hfill \\ \\ {\displaystyle {}_0{D}^{\epsilon }A\left(t\right)=(1-u)\beta SI+\sigma \beta (1-u)VI-(\delta +\varphi ),A\left(0\right)>0,}\hfill \\ \\ {\displaystyle {}_0{D}^{\epsilon }V\left(t\right)=v\alpha S-(1-u)\sigma \beta VI-\delta V,V\left(0\right)>0,}\hfill \\ \\ {\displaystyle {}_0{D}^{\epsilon }I\left(t\right)=\xi \varphi A-({\delta }_2+\rho )I,I\left(0\right)>0}\hfill \\ \\ {\displaystyle {}_0{D}^{\epsilon }C\left(t\right)=(1-\xi )\varphi A-({\delta }_1+\gamma )C,C\left(0\right)>0,}\hfill \\ \\ {\displaystyle {}_0{D}^{\epsilon }R\left(t\right)=\rho I+\gamma C-\delta R,R\left(0\right)>0.}\hfill \end{array}$$

(17)

The next theorem shows the existence of optimal solution of control problem.

Theorem 9.

Suppose the control function $x=(u,v)\in U$ be measurable on $[0,T_f]$ and $0\le u,v\le 1$. Then there subsists an optimal control $x^{*}=(u^{*},v^{*})$, which minimizes $J(u,v)$ of (20) with

$$\begin{array}{c}{}_t^{RL}{D}_{T_f}^{\epsilon }{\lambda }_1\left(t\right)=(1-u)({\lambda }_2-{\lambda }_1)+{\lambda }_1(\delta +\alpha v)\hfill \\ \\ -{\lambda }_3\alpha v\hfill \\ \\ {}_t^{RL}{D}_{T_f}^{\epsilon }{\lambda }_2\left(t\right)={\lambda }_2(\delta +\varphi )-{\lambda }_4\left(\xi \varphi \right)-{\lambda }_5\left\{(1-\xi )\varphi \right\}\hfill \\ \\ {}_t^{RL}{D}_{T_f}^{\epsilon }{\lambda }_3\left(t\right)={\lambda }_3(\delta +\sigma \beta (1-u)I)-{\lambda }_2\sigma \beta (1-u)I,\hfill \\ \\ {}_t^{RL}{D}_{T_f}^{\epsilon }{\lambda }_4\left(t\right)={\lambda }_3\sigma \beta (1-u)V-{\lambda }_2\{\sigma \beta (1-u)V+(1-u)\beta S\}\hfill \\ \\ +(1-u)\beta S{\lambda }_1+{\lambda }_4(\delta +\rho )-{\lambda }_6\rho -m_1,\hfill \\ \\ {}_t^{RL}{D}_{T_f}^{\epsilon }{\lambda }_5\left(t\right)={\lambda }_5({\delta }_1+\gamma )-{\lambda }_6\gamma -m_2,\hfill \\ \\ {}_t^{RL}{D}_{T_f}^{\epsilon }{\lambda }_6\left(t\right)=\delta {\lambda }_6.\hfill \\ \end{array}$$

where ${\lambda }_i\left(T_f\right)=0,i=1,2,...,6$ and

$$\begin{array}{c}{u}^{*}=max\{min\{\overline{u},1\},0\};v^{*}=max\{min\{\overline{v},1\},0\}\hfill \\ \\ {\displaystyle \overline{u}=\frac{({\lambda }_2-{\lambda }_1)S^{*}{I}^{*}\beta +\sigma \beta V^{*}{I}^{*}}{m_3},\overline{v}=\frac{({\lambda }_1-{\lambda }_3)\alpha I^{*}{S}^{*}}{m_4}}\hfill \end{array}$$

(18)

and $S^{*},A^{*},V^{*},I^{*},C^{*},R^{*}$ are the corresponding optimal state variables of (17).

Proof. 

Consider the following Hamiltonian function:

$$\begin{array}{c}{\displaystyle H=m_{1}I\left(t\right)+m_{2}C\left(t\right)+\frac{m_3}{2}{u}^2\left(t\right)+\frac{m_4}{2}{v}^2\left(t\right)+{\lambda }_1{(}_0{D}^{\epsilon }S\left(t\right))+{\lambda }_2{(}_0{D}^{\epsilon }A\left(t\right))+{\lambda }_3{(}_0{D}^{\epsilon }V\left(t\right))}\hfill \\ \\ {\displaystyle +{\lambda }_4{(}_0{D}^{\epsilon }I\left(t\right))+{\lambda }_5{(}_0{D}^{\epsilon }C\left(t\right))+{\lambda }_6{(}_0{D}^{\epsilon }R\left(t\right))}\hfill \end{array}$$

(19)

where adjoint variables are denoted by ${\lambda }_i,i=1,2,3,4,5,6$ with ${\lambda }_i\left(T_f\right)=0(i=1,2,3,4,5,6$).

For control system (17), let us consider that $u^{*},v^{*}$ are applied optimal control interventions with optimal state variables $S^{*},A^{*},V^{*},I^{*},C^{*},R^{*}$:

$$\begin{array}{c}{\displaystyle {}_t^{RL}{D}_{T_f}^{\epsilon }{\lambda }_1\left(t\right)=-\frac{\partial H}{\partial S}=(1-u)({\lambda }_2-{\lambda }_1)+{\lambda }_1(\delta +\alpha v)}\hfill \\ \\ -{\lambda }_3\alpha v\hfill \\ \\ {\displaystyle {}_t^{RL}{D}_{T_f}^{\epsilon }{\lambda }_2\left(t\right)=-\frac{\partial H}{\partial A}={\lambda }_2(\delta +\varphi )-{\lambda }_4\left(\xi \varphi \right)-{\lambda }_5\left\{(1-\xi )\varphi \right\}}\hfill \\ \\ {\displaystyle {}_t^{RL}{D}_{T_f}^{\epsilon }{\lambda }_3\left(t\right)=-\frac{\partial H}{\partial V}={\lambda }_3(\delta +\sigma \beta (1-u)I)-{\lambda }_2\sigma \beta (1-u)I,}\hfill \\ \\ {\displaystyle {}_t^{RL}{D}_{T_f}^{\epsilon }{\lambda }_4\left(t\right)=-\frac{\partial H}{\partial I}={\lambda }_3\sigma \beta (1-u)V-{\lambda }_2\{\sigma \beta (1-u)V+(1-u)\beta S\}}\hfill \\ \\ +(1-u)\beta S{\lambda }_1+{\lambda }_4(\delta +\rho )-{\lambda }_6\rho -m_1,\hfill \\ \\ {\displaystyle {}_t^{RL}{D}_{T_f}^{\epsilon }{\lambda }_5\left(t\right)=-\frac{\partial H}{\partial C}={\lambda }_5({\delta }_1+\gamma )-{\lambda }_6\gamma -m_2,}\hfill \\ \\ {\displaystyle {}_t^{RL}{D}_{T_f}^{\epsilon }{\lambda }_6\left(t\right)=-\frac{\partial H}{\partial R}=\delta {\lambda }_6.}\hfill \end{array}$$

(20)

From optimality conditions [45], we have obtained the optimal conditions:

$$\begin{array}{c}{\displaystyle \frac{\partial H}{\partial u}{|}_{u=\overline{u}}=m_3\overline{u}+{\lambda }_1\beta S^{*}{I}^{*}-{\lambda }_2(\beta S^{*}{I}^{*}+\sigma \beta V^{*}{I}^{*})+{\lambda }_3\sigma \beta V^{*}{I}^{*}=0}\hfill \\ \\ {\displaystyle \frac{\partial H}{\partial v}{|}_{v=\overline{v}}=m_4\overline{v}-{\lambda }_1{S}^{*}\alpha +{\lambda }_3\alpha S^{*}=0}\hfill \end{array}$$

(21)

So, we have

$${\displaystyle \overline{u}=\frac{({\lambda }_2-{\lambda }_1)S^{*}{I}^{*}\beta +\sigma \beta V^{*}{I}^{*}}{m_3},\overline{v}=\frac{({\lambda }_1-{\lambda }_3)\alpha I^{*}{S}^{*}}{m_4}}$$

(22)

So, in U, we have

$$u^{*}=\left\{\begin{array}{c}{\displaystyle 0\mathrm{if}\frac{\partial H}{\partial u}<0}\hfill \\ \\ {\displaystyle \overline{u}\mathrm{if}\frac{\partial H}{\partial u}=0}\hfill \\ \\ {\displaystyle 1\mathrm{if}\frac{\partial H}{\partial u}>0}\hfill \end{array}\right.$$

(23)

$$v^{*}=\left\{\begin{array}{c}{\displaystyle 0\mathrm{if}\frac{\partial H}{\partial v}<0}\hfill \\ \\ {\displaystyle \overline{v}\mathrm{if}\frac{\partial H}{\partial v}=0}\hfill \\ \\ {\displaystyle 1\mathrm{if}\frac{\partial H}{\partial v}>0}\hfill \end{array}\right.$$

(24)

and

$$u^{*}=max\{min\{\overline{u},1\},0\};v^{*}=max\{min\{\overline{v},1\},0\}$$

(25)

where ${\displaystyle \overline{u}=\frac{({\lambda }_2-{\lambda }_1)S^{*}{I}^{*}\beta +\sigma \beta V^{*}{I}^{*}}{m_3},\overline{v}=\frac{({\lambda }_1-{\lambda }_3)\alpha I^{*}{S}^{*}}{m_4}.}$ □

6. Numerical Simulations

For numerical purposes, we have used the MatLab software in combination with Roberto Garrappa’s fd12 Matlab function for fractional differential equations [46]. For the optimal control problem, the iterative scheme (Euler’s forward and backward) has been used [15].

Table 4. Day by day case report between 1 June 2021 and 20 October 2021 of West Bengal (India).

Day	Total Infected Cases	Model Value
1/6/2021	9424	9000
6/6/2021	7682	3385
11/6/201	5274	2978
16/6/2021	3268	1925
21/6/2021	2184	1651
26/6/2021	1933	1444
1/07/2021	1478	1286
6/07/2021	1367	1163
11/07/2021	997	1054
16/07/2021	831	962
21/07/2021	752	893
26/07/2021	657	832
31/07/2021	769	777
05/08/2021	732	728
10/08/2021	639	683
15/08/2021	663	648
20/08/2021	673	758
25/08/2021	758	685
30/08/2021	652	557
04/09/2021	677	590
09/09/2021	680	560
14/09/2021	702	540
19/09/2021	619	500
24/09/2021	707	470
30/09/2021	651	452
05/10/2021	618	420
10/10/2021	619	405
15/10/2021	530	393
20/10/2021	857	380

depicts the pandemic scenario of West Bengal in the period 1 June 2021 to 20 October 2021. We have also considered system (1) where all parameters contain $\epsilon$. The value of reproduction number $R_0$ is 0.0152 according to Table 2, which suggests the non-existence of an endemic equilibrium point and the existence of asymptotically stable disease free equilibrium. Numerical simulations have been compared to the model’s conclusions to real-world data from WHO [47] and worldometers [48] reports. The total population of West Bengal is around ${\displaystyle 9.2\times {10}^7}$ [1]. We considered $t=1$ day as the minimum time unit and end time is $T_f=60$ days. We have also considered $\Phi =I+C$ as the total infected population due to COVID-19 in West Bengal. The symptomatically infected cases (with or without co-morbidity), from 1 June 2021 to 20 October 2021, are listed in Table 4 with 5 days interval [1]. Initially, we have simulated our system for 60 days, and then we extended our limit to 150 days. We saw that in the beginning, the time series of total infected population ($\Phi =I+C$) does not quite fit with actual data; however, around t = 30 days, the model data fit well with real data when the order of the derivative was 0.7 (Figure 2). Next, we performed numerical simulation using Table 2 from 1 June to 30 October 2021 (Figure 3) and real data are taken up to 20 October 2021.

For the optimal control problem, we have assumed $m_1=100,m_2=100,m_3=10,m_4=1000$, all are in INR (Indian Rupee). It is observed from Figure 4 and Figure 5 that infected population with co-morbidity is decreasing if we control both factors (vaccination and social distancing) simultaneously. After effectively executing control policy, the number of infected individuals with co-morbidity will be lowered from 15,860 to 14,970 (Figure 4 and Figure 5). Figure 6 depicts the time series of the optimal control variables and the optimal cost of implementation of the control policy. It has been noted that the effect of controls $u^{*},v^{*}$ is minimal in the time span of 60 days. Figure 7 suggests the variation of total infected population ($\Phi$) when the combined control strategy is applied. For $\epsilon =0.7$, we saw that the population of the total infected class ($\Phi$) is not far from the realistic scenario for West Bengal (Figure 8 and Figure 9). We have estimated $\epsilon$ by graphical method using MatLab (Figure 9). Other parameters have been estimated in similar graphical method using MatLab as mentioned in Table 2 (Figure 10). In Figure 10, the variation of time series of different state variables have been produced for different values of $\epsilon$ for the long time period (t = 300). Figure 11, Figure 12 and Figure 13 depict the time series of state variables, optimal control variables, and optimal cost for longer time periods (t = 250 days).

7. Conclusions

In real-world dynamical processes, such as epidemic spread, fractional calculus plays a critical role. We looked into the emergence of a modified SAVICR epidemic model that included memory effects. The co-morbidity and vaccinated classes have been added. Our model accurately simulates the reality of the West Bengal (India) outbreak and predicts the daily number of confirmed COVID-19 cases. The size of the infected population with co-morbidity is important in estimating the number of Intensive Care Units (ICU) needed. Our proposed model is compatible with other places/countries. The order of derivative can differ from region to region. If we vary the order of derivatives while keeping other parametric values fixed, the results will be different (Figure 3). This demonstrates that the order of derivative is important in system simulation (1). We can see that the number of people with a co-morbidity and asymptomatic population decreases significantly at 200 days by conducting long-term simulation (Figure 10). Finally, it has shown that the optimal control strategy is implementing social distancing policies and lowering economic costs by combining non-pharmaceutical interventions with vaccination.

Our work reveals the validity of the proposed model in West Bengal (India) scenario. The main advantage of this model is that it is simple and fits with real data for long time periods, as mentioned in Figure 3 (1 June 2021 to 20 October 2021). Theoretically, it is established that the endemic equilibrium point is not asymptotically stable, which guarantees the extinction of the current pandemic (Theorem 6). One disadvantage of this model is that the model values and real data of the total infected population differ too much in the initial stage; further, vaccinated people may be affected by COVID-19, which is a factor not included in our model. It is observed that the number of infected individuals can be controlled by proper vaccination and NPI, and the strategies proposed to control COVID-19 through proper implementation of NPI and vaccination events can contribute to our society. Further, the study of the pulsed immunization strategy remains within the scope of our future studies. There are several fractional-order systems, namely Caputo–Fabrizio, Caputo–Hadamard, Atangana–Baleanu–Caputo derivative, and Reisz derivative, which are extensively used for studying epidemiological and biological systems. We wish to conduct a comparative study of the proposed model in different fractional-order frameworks.