Next Article in Journal
Specialised Knowledge for Teaching Geometry in a Primary Education Class: Analysis from the Knowledge Mobilized by a Teacher and the Knowledge Evoked in the Researcher
Next Article in Special Issue
Fractional Growth Model with Delay for Recurrent Outbreaks Applied to COVID-19 Data
Previous Article in Journal
Deep Learning Approach to Mechanical Property Prediction of Single-Network Hydrogel
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Fractional Ordered COVID-19 Model Incorporating Comorbidity and Vaccination

1
Department of Mathematics, Indian Institute of Engineering Science and Technology, Shibpur, Howrah 711103, India
2
Institute of Research and Development of Processes, University of the Basque Country, 48940 Leioa, Spain
*
Author to whom correspondence should be addressed.
Mathematics 2021, 9(21), 2806; https://doi.org/10.3390/math9212806
Submission received: 29 September 2021 / Revised: 30 October 2021 / Accepted: 1 November 2021 / Published: 4 November 2021
(This article belongs to the Special Issue Mathematical Methods and Models in Epidemiology)

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 ε ( 0 , 1 ) for a absolute continuous and differentiable function g on [ 0 , + ) is defined as:
0 c D t ε g ( t ) = 1 Γ ( 1 ε ) t 0 t g ( s ) ( t s ) ε d s g ( t ) , ε = 1 .
In our whole context, we have used 0 D ε instead of 0 c D t ε and 0 < ε < 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:
t 0 c D t ε S ( t ) = Λ ε δ ε S ( t ) α ε S ( t ) β ε S ( t ) I ( t ) , S ( 0 ) > 0 , t 0 c D t ε A ( t ) = β ε S ( t ) I ( t ) ϕ ε A ( t ) + σ ε β ε V ( t ) I ( t ) δ ε A ( t ) , A ( 0 ) 0 , t 0 c D t ε V ( t ) = α ε S ( t ) σ ε β ε V ( t ) I ( t ) δ ε V ( t ) , V ( 0 ) > 0 , t 0 c D t ε I ( t ) = ξ ϕ ε A ( t ) ( δ 2 ε + ρ ε ) I ( t ) , I ( 0 ) 0 , t 0 c D t ε C ( t ) = ( 1 ξ ) ϕ ε A ( t ) ( δ 1 ε + γ ε ) C ( t ) , C ( 0 ) 0 , t 0 c D t ε R ( t ) = ρ ε I ( t ) + γ ε C ( t ) δ ε R ( t ) , R ( 0 ) 0 .
Here t 0 c D t ε is the notation of fractional derivative of order ε ( 0 , 1 ) in Caputo sense with initial time t 0 0 (assuming that t 0 = 0 ). For simplicity we use 0 D ε instead of t 0 c D t ε . System 1 is dimensionally correct as both sides have same time dimension t i m e ε . It is also noticed that the equilibrium points and reproduction number contain ε as all the parameters contain ε in power form. We have omitted the power ε of all parameters (for simplicity) in theoretical analysis. In numerical simulations, we have taken into account the power ε of all parameters. Now, system 1 transforms to:
0 D ε S ( t ) = Λ δ S ( t ) α S ( t ) β S ( t ) I ( t ) , S ( 0 ) > 0 , 0 D ε A ( t ) = β S ( t ) I ( t ) ϕ A ( t ) + σ β V ( t ) I ( t ) δ A ( t ) , A ( 0 ) 0 , 0 D ε V ( t ) = α S ( t ) σ β V ( t ) I ( t ) δ V ( t ) , V ( 0 ) > 0 , 0 D ε I ( t ) = ξ ϕ A ( t ) δ 2 + ρ I ( t ) , I ( 0 ) 0 0 D ε C ( t ) = ( 1 ξ ) ϕ A ( t ) δ 1 + γ C ( t ) , C ( 0 ) 0 0 D ε R ( t ) = ρ I ( t ) + γ C ( t ) δ R ( t ) , R ( 0 ) 0
It is presumed δ δ 2 δ 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 ( t ) , A ( t ) , V ( t ) , I ( t ) , C ( t ) , and R ( t ) , respectively. Table 1 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 = Λ δ + α , 0 , α Λ δ ( δ + α ) , 0 , 0 , 0
  • E 1 = ( S * , A * , V * , I * , C * , R * ) .
Here
S * = Λ α + δ + β I * A * = δ 2 + ρ ξ ϕ I * V * = α Λ [ α + δ + β I * ] [ σ β I * + δ ] I * = Λ ξ ϕ ( ϕ + δ ) ( δ 2 + ρ ) δ σ β C * = ( 1 ξ ) ( δ 2 + ρ ) ξ ( γ + δ 1 ) I * R * = ρ δ + γ ( 1 ξ ) ( δ 2 + ρ ) δ ξ ( δ 1 + γ ) I *
E 1 remains in R + 6 if
Λ ξ ϕ σ β δ ( ϕ + δ ) ( δ 2 + ρ )
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 F ¯ , V ¯ represent respective new infection and transition matrices.
Consider u = ( A , I , C ) T , the subsystem of system (2) can be represented as:
0 D t ε u = F ¯ ( u ) V ¯ ( u ) ,
F ¯ ( u ) = β S I + σ β V I 0 0
and
V ¯ ( u ) = ϕ A + δ A ( δ 2 + ρ ) I ξ ϕ A ( δ 1 + γ ) C ( 1 ξ ) ϕ A
D F ¯ ( E 0 ) = F
D V ¯ ( E 0 ) = V
D F ¯ ( u ) and D V ¯ ( u ) are the Jacobian of F ¯ , V ¯ at disease-free equilibrium E 0 where
F = 0 Λ β δ + α + σ β α Λ δ ( δ + α ) 0 0 0 0 0 0 0
V = ( δ + ϕ ) 0 0 ξ ϕ δ 2 + ρ 0 ( 1 ξ ) ϕ 0 δ 1 + γ
Thus, we obtain
R 0 = ξ β Λ ϕ ( ϕ + δ ) ( δ + α ) ( δ 2 + ρ ) + α ξ ϕ β σ Λ δ ( ϕ + δ ) ( δ + α ) ( δ 2 + ρ ) .
Reproduction number R 0 is a function of ε because each parameter is a function of ε . The value of ε has been fixed for analysis reasons. If we modify the value of ε , 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 ε ψ ( t ) is continuous on a , b , then
ψ ( x ) = ψ ( a ) + 1 Γ ( ε ) ( x a ) ε . 0 D ε ψ ( η ) ,
where 0 < a η x , x a , b .
Remark 1. 
If 0 D ε ψ ( t ) 0 0 D ε ψ ( t ) 0 for all t ( a , b ) , then ψ ( t ) is a non-decreasing (non-increasing) function for t a , b .
Definition 2 
(Ref. [8]). Mittag–Leffler functions of single and double parameters are described below:
E ε ( w ) = k = 0 w k Γ ( ε k + 1 ) and E ε 1 , ε 2 ( w ) = k = 0 w k Γ ( ε 1 k + ε 2 ) , where ε , ε 1 , ε 2 R + .
Theorem 1 
(Ref. [31]). Suppose g ( t ) is differentiable with exponential order and 0 < ε < 1 . Further, if 0 D ε g ( t ) is piece-wise continuous on 0 , , then
L 0 D ε g ( t ) = s ε F ( s ) j = 0 n 1 s ε j 1 g j ( 0 ) ,
where the Laplace transform of g ( t ) is represented by F ( s ) = L g ( t ) .
Theorem 2 
(Ref. [32]). For any complex number M and ε 1 , ε 2 ( 0 , 1 ) , we have
L t ε 2 1 E ε 1 , ε 2 ( M t ε 1 ) = s ε 1 ε 2 ( s ε 1 M ) ,
where the real part of s > M 1 ε 1 and E ε 1 , ε 2 is the two parametric Mittag–Leffler function.
Theorem 3 
(Ref. [27]). Consider:
0 D t ε x ( t ) = Φ x ,
with 0 < ε < 1 , x R n ) . The equilibrium points of the above system are solutions to the equation Φ ( x ) = 0 . An equilibrium is locally asymptotically stable if all eigenvalues ( λ i ) of the Jacobian matrix J = Φ x calculated at the equilibrium points satisfy arg ( λ i ) > ε π 2 .

4.2. Existence and Uniqueness

Lemma 2 
(Ref. [33]). Consider the system:
0 D ε x ( t ) = g ( t , x )
with initial condition x ( 0 ) = x 0 , where ε 0 , 1 , g : 0 , × Ω I R n , Ω I R n , if g ( t , x ) satisfies local Lipschitz condition with respect to x, then there exists a unique solution of (5) on 0 , × Ω .
In our context Ω = { ( S , A , V , I , C , R ) R 6 : max ( S , A , V , I , C , R ) M < } .
Now, we have the following theorem.
Theorem 4.
The system (2) with initial condition X ( 0 ) = ( S ( 0 ) , A ( 0 ) , V ( 0 ) , I ( 0 ) , C ( 0 ) , R ( 0 ) ) Ω , possess a unique solution X ( t ) Ω , t 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 S A V I C R model (2) with initial values lie in Γ + = R + 6 .
Proof. 
0 D ε S ( t ) | S ( t ) = 0 = Λ > 0
0 D ε A ( t ) | A ( t ) = 0 = β S I + σ β V I
0 D ε V ( t ) | V ( t ) = 0 = α S
0 D ε I ( t ) | I ( t ) = 0 = ξ ϕ A
0 D ε C ( t ) | C ( t ) = 0 = ( 1 ξ ) ϕ A
0 D ε R ( t ) | R ( t ) = 0 = ρ I + γ C
From (6), we have
0 D ε S ( t ) | S ( t ) = 0 = Λ > 0 .
We can deduce from Lemma 1 that S ( t ) is increasing in near S ( t ) = 0 , and S ( t ) is unable to go across the axis S ( t ) = 0 . As a result, S is non-negative and similarly we can show V ( t ) > 0 for all t 0 (from Equation (8)). We now assert that A ( t ) commences in Γ + cannot be negative. Contrarily, there exists τ 1 such that A ( t ) traverses the A ( t ) = 0 axis for the first time at t = τ 1 with the following possibilities.
A ( t ) > 0 , for 0 t < τ 1 , A ( τ 1 ) = 0 , A ( τ 1 + ) < 0 .
Relation (7) revels that 0 D ε A ( t ) | A ( τ 1 ) = 0 = β S ( τ 1 ) I ( τ 1 ) + σ β V ( τ 1 ) I ( τ 1 ) . and this leads the following cases.
  • Case1: Suppose β S ( τ 1 ) I ( τ 1 ) + σ β V ( τ 1 ) I ( τ 1 ) 0 , then from Lemma 1, A ( t ) is not decreasing in the neighborhood of t = τ 1 , and A ( τ 1 + ) = 0 . As a result, a contradiction occurs.
  • Case2: On the other hand β S ( τ 1 ) I ( τ 1 ) + σ β V ( τ 1 ) I ( τ 1 ) < 0 , which discloses I ( τ 1 ) must be negative. If I ( τ 1 ) < 0 , then a τ 2 ( 0 < τ 2 < τ 1 ) exists with
    I ( t ) > 0 , for 0 t < τ 2 , I ( τ 2 ) = 0 , I ( τ 2 + ) < 0 .
From (9), we obtain
0 D t ε I ( t ) | I ( τ 2 ) = 0 = ξ ϕ A ( τ 2 ) > 0 ,
which refutes I ( τ 2 + ) < 0 . Therefore, I ( t ) 0 , t [ 0 , ) and also A ( t ) 0 , t [ 0 , ) .
Once again from (10) and (11), we can conclude C ( t ) , R ( t ) are non-negative.
As a result, all solutions of system (2) that begin with Γ + are constricted to the region Γ + . □
Theorem 6 
(Boundedness). System (1) has uniformly bounded solutions X ( t ) = ( S , A , V , I , C , R ) .
Proof. 
It is clear from the first equation of (2) that
0 D ε S ( t ) Λ δ S
Using Laplace transforms ( L · ) , we have
s ε L S ( t ) s ε 1 S ( 0 ) + δ L S ( t ) Λ s , L S ( t ) A s ε ( 1 + ε ) s ε + δ + S ( 0 ) s ε 1 s ε + δ
Using inverse Laplace transforms (using Theorem 2):
S ( t ) S ( 0 ) E ε , 1 ( δ t ε ) + Λ t ε E ε , ε + 1 ( δ t ε )
Therefore, S ( t ) M 1 [ E ε , 1 ( δ t ε ) + δ t ε E ε , ε + 1 ( δ t ε ) ] = M 1 Γ ( 1 ) = M 1 , M 1 = max Λ δ , S ( 0 ) .
From the properties of Mittag–Leffler function: [34]:
E α 1 , α 2 ( z ) = z E α 1 , α 1 + α 2 ( z ) + 1 Γ ( α 2 )
In this circumstance,
E ε , 1 ( δ t ε ) = ( δ t ε ) E ε , ε + 1 ( δ t ε ) + 1 Γ ( 1 )
Let N ( t ) = S ( t ) + A ( t ) + V ( t ) + I ( t ) + C ( t ) + R ( t ) .
Now,
0 D t ε N ( t ) = 0 D t ε S ( t ) + 0 D t ε A ( t ) + 0 D t ε V ( t ) + 0 D t ε I ( t ) + 0 D t ε C ( t ) + 0 D t ε R ( t ) = Λ { δ S ( t ) + δ A ( t ) + δ V ( t ) + δ 2 I ( t ) + δ 1 C ( t ) + δ R ( t ) } Λ δ m N ( t ) , where δ m = min { δ 1 , δ 2 , δ }
Hence,
0 D t ε N ( t ) + δ m N ( t ) Λ
Using Laplace transformation, we have (using Theorem 1):
s ε F ( s ) s ε 1 N ( 0 ) + δ m F ( s ) Λ s , where F ( s ) = L N ( t ) F ( s ) Λ s 1 s ε + δ m + N ( 0 ) s ε 1 s ε + δ m = s ε 1 N ( 0 ) s ε + δ m + Λ s ε ( 1 + ε ) s ε + δ m
Using inverse Laplace transforms (using Theorem 2):
N ( t ) = N ( 0 ) E ε , 1 ( δ m t ε ) + Λ t ε E ε , ε + 1 ( δ m t ε )
From (12) and (13), we obtain
N ( t ) M 2 [ E ε , 1 ( δ m t ε ) + δ m t ε E ε , ε + 1 ( δ m t ε ) ] = M 2 Γ ( 1 ) = M 2 , M 2 = max Λ δ m , N ( 0 )
Thus S ( t ) , N ( t ) are bounded and hence the solution ( S ( t ) , A ( t ) , V ( t ) , I ( t ) , C ( t ) , R ( t ) ) is bounded in { ( S , A , V , I , C , R ) | S + A + V + I + C + R M 2 ; S M 1 } for t [ 0 , ) . □

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:
0 D ε S ( t ) = Λ δ S ( t ) α S ( t ) β S ( t ) I ( t ) , S ( 0 ) > 0 , 0 D ε A ( t ) = β S ( t ) I ( t ) ϕ A ( t ) + σ β V ( t ) I ( t ) δ A ( t ) , A ( 0 ) 0 , 0 D ε V ( t ) = α S ( t ) σ β V ( t ) I ( t ) δ V ( t ) , V ( 0 ) > 0 , 0 D ε I ( t ) = ξ ϕ A ( t ) δ 2 + ρ I ( t ) , I ( 0 ) 0 0 D ε C ( t ) = ( 1 ξ ) ϕ A ( t ) δ 1 + γ C ( t ) , C ( 0 ) 0
In this casewe have to analyze the stability of equilibriumpoints E 0 * = Λ δ + α , 0 , α Λ δ ( δ + α ) , 0 , 0 and E 1 * = S * , A * , V * , I * , C * , where
S * = Λ α + δ + β I * A * = δ 2 + ρ ξ ϕ I * V * = α Λ [ α + δ + β I * ] [ σ β I * + δ ] I * = Λ ξ ϕ ( ϕ + δ ) ( δ 2 + ρ ) δ σ β = δ σ β Λ ξ ϕ σ β δ ( ϕ + δ ) ( δ 2 + ρ ) 1 = R 0 K 0 1 δ σ β C * = ( 1 ξ ) ( δ 2 + ρ ) ξ ( γ + δ 1 ) I *
and R 0 = ξ β Λ ϕ ( ϕ + δ ) ( δ + α ) δ 2 + ρ + α ξ ϕ β σ Λ δ ( ϕ + δ ) ( δ + α ) δ 2 + ρ , K 0 = δ ( δ + α ) σ + α α + δ
For E 1 * to exist in feasible region R + 5 , it is necessary and sufficient that R 0 > K 0 . From simple calculation it is clear that K 0 1 if σ 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 ( λ ) = λ n + a 1 λ n 1 + a 2 λ n 2 + . . . + a n be the characteristic equation of system ( 14 ) evaluated at equilibrium point. For n = 4 , ( P ) is the discriminant of the characteristic equation P ( λ ) = λ 4 + a 1 λ 3 + a 2 λ 2 + a 3 λ + a 4 which is represented as ( P ) = 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 3 a 1 2 a 2 a 3 0 0 0 0 4 3 a 1 2 a 2 a 3 0 0 0 0 4 3 a 1 2 a 2 a 3 0 0 0 0 4 3 a 1 2 a 2 a 3
Lemma 4 
(Ref. [37]). Suppose ( P ) < 0 , a 1 > 0 , a 2 > 0 , a 3 > 0 , a 4 > 0 hold, then the equilibrium point is asymptotically stable if a 2 = a 1 a 4 a 3 + a 3 a 1 . If a 2 a 1 a 4 a 3 + a 3 a 1 then ε < 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 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, θ ) [29,38]:
Ω θ R 0 = R 0 θ θ 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 Λ , α , β , ϕ , ρ , a n d σ because we could not manage δ , ξ . Table 3 shows the sensitivity index values for the parameters Λ , α , β , ϕ , ρ , σ that match Table 2. From Table 3, it is revealed that infection transmission rate β , 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 R L D T ε ) of order ε 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 u , v 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 = ( u , v ) L 1 ( 0 , T f ) , 0 u , v 1 , t ( 0 , T f ) 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 ( u ( t ) , v ( t ) ) U
J ( u , v ) = 0 T f m 1 I ( t ) + m 2 C ( t ) + m 3 2 u 2 ( t ) + m 4 2 v 2 ( t ) d t
subject to
0 D ε S ( t ) = Λ ( δ + α v ) S β S I ( 1 u ) , S ( 0 ) > 0 , 0 D ε A ( t ) = ( 1 u ) β S I + σ β ( 1 u ) V I ( δ + ϕ ) , A ( 0 ) > 0 , 0 D ε V ( t ) = v α S ( 1 u ) σ β V I δ V , V ( 0 ) > 0 , 0 D ε I ( t ) = ξ ϕ A ( δ 2 + ρ ) I , I ( 0 ) > 0 0 D ε C ( t ) = ( 1 ξ ) ϕ A ( δ 1 + γ ) C , C ( 0 ) > 0 , 0 D ε R ( t ) = ρ I + γ C δ R , R ( 0 ) > 0 .
The next theorem shows the existence of optimal solution of control problem.
Theorem 9.
Suppose the control function x = ( u , v ) U be measurable on [ 0 , T f ] and 0 u , v 1 . Then there subsists an optimal control x * = ( u * , v * ) , which minimizes J ( u , v ) of (20) with
t R L D T f ε λ 1 ( t ) = ( 1 u ) ( λ 2 λ 1 ) + λ 1 ( δ + α v ) λ 3 α v t R L D T f ε λ 2 ( t ) = λ 2 ( δ + ϕ ) λ 4 ( ξ ϕ ) λ 5 { ( 1 ξ ) ϕ } t R L D T f ε λ 3 ( t ) = λ 3 ( δ + σ β ( 1 u ) I ) λ 2 σ β ( 1 u ) I , t R L D T f ε λ 4 ( t ) = λ 3 σ β ( 1 u ) V λ 2 { σ β ( 1 u ) V + ( 1 u ) β S } + ( 1 u ) β S λ 1 + λ 4 ( δ + ρ ) λ 6 ρ m 1 , t R L D T f ε λ 5 ( t ) = λ 5 ( δ 1 + γ ) λ 6 γ m 2 , t R L D T f ε λ 6 ( t ) = δ λ 6 .
where λ i ( T f ) = 0 , i = 1 , 2 , . . . , 6 and
u * = max { min { u ¯ , 1 } , 0 } ; v * = max { min { v ¯ , 1 } , 0 } u ¯ = ( λ 2 λ 1 ) S * I * β + σ β V * I * m 3 , v ¯ = ( λ 1 λ 3 ) α I * S * m 4
and S * , A * , V * , I * , C * , R * are the corresponding optimal state variables of (17).
Proof. 
Consider the following Hamiltonian function:
H = m 1 I ( t ) + m 2 C ( t ) + m 3 2 u 2 ( t ) + m 4 2 v 2 ( t ) + λ 1 ( 0 D ε S ( t ) ) + λ 2 ( 0 D ε A ( t ) ) + λ 3 ( 0 D ε V ( t ) ) + λ 4 ( 0 D ε I ( t ) ) + λ 5 ( 0 D ε C ( t ) ) + λ 6 ( 0 D ε R ( t ) )
where adjoint variables are denoted by λ i , i = 1 , 2 , 3 , 4 , 5 , 6 with λ i ( T f ) = 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 * :
t R L D T f ε λ 1 ( t ) = H S = ( 1 u ) ( λ 2 λ 1 ) + λ 1 ( δ + α v ) λ 3 α v t R L D T f ε λ 2 ( t ) = H A = λ 2 ( δ + ϕ ) λ 4 ( ξ ϕ ) λ 5 { ( 1 ξ ) ϕ } t R L D T f ε λ 3 ( t ) = H V = λ 3 ( δ + σ β ( 1 u ) I ) λ 2 σ β ( 1 u ) I , t R L D T f ε λ 4 ( t ) = H I = λ 3 σ β ( 1 u ) V λ 2 { σ β ( 1 u ) V + ( 1 u ) β S } + ( 1 u ) β S λ 1 + λ 4 ( δ + ρ ) λ 6 ρ m 1 , t R L D T f ε λ 5 ( t ) = H C = λ 5 ( δ 1 + γ ) λ 6 γ m 2 , t R L D T f ε λ 6 ( t ) = H R = δ λ 6 .
From optimality conditions [45], we have obtained the optimal conditions:
H u | u = u ¯ = m 3 u ¯ + λ 1 β S * I * λ 2 ( β S * I * + σ β V * I * ) + λ 3 σ β V * I * = 0 H v | v = v ¯ = m 4 v ¯ λ 1 S * α + λ 3 α S * = 0
So, we have
u ¯ = ( λ 2 λ 1 ) S * I * β + σ β V * I * m 3 , v ¯ = ( λ 1 λ 3 ) α I * S * m 4
So, in U, we have
u * = 0 if H u < 0 u ¯ if H u = 0 1 if H u > 0
v * = 0 if H v < 0 v ¯ if H v = 0 1 if H v > 0
and
u * = max { min { u ¯ , 1 } , 0 } ; v * = max { min { v ¯ , 1 } , 0 }
where u ¯ = ( λ 2 λ 1 ) S * I * β + σ β V * I * m 3 , v ¯ = ( λ 1 λ 3 ) α 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 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 ε . 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 9.2 × 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 Φ = 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 ( Φ = 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 ( Φ ) when the combined control strategy is applied. For ε = 0.7 , we saw that the population of the total infected class ( Φ ) is not far from the realistic scenario for West Bengal (Figure 8 and Figure 9). We have estimated ε 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 ε 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.

Author Contributions

Conceptualization, M.D., G.S. and M.D.l.S.; Methodology, M.D., G.S. and M.D.l.S.; Investigation, M.D., G.S. and M.D.l.S.; Formal analysis, M.D., G.S. and M.D.l.S.; Writing—original draft preparation, M.D., G.S. and M.D.l.S.; Writing—review and editing, M.D., G.S. and M.D.l.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Spanish Government for its support through grant RTI2018-094336-B-I00 (MCIU/AEI/FEDER, UE) and to the Basque Government for its support through grant IT1207-19.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data used to support the findings of this study are included in the references within the article.

Acknowledgments

The authors are grateful to the anonymous referees and Hannah Zheng, Assistant Editor for their careful reading, valuable comments, and helpful suggestions, which have helped them to improve the presentation of this work significantly. The third author (Manuel De la Sen) is grateful to the Spanish Government for its support through grant RTI2018-094336-B-I00 (MCIU/AEI/FEDER, UE) and to the Basque Government for its support through grant IT1207-19.

Conflicts of Interest

The authors declare that they have no conflict of interest regarding this work.

Appendix A

Proof of Theorem 4. 
Let us denote X = ( S , A , V , I , C , R ) and X ¯ = ( S ¯ , A ¯ , V ¯ , I ¯ , C ¯ , R ¯ ) . Consider a map L ( X ) = ( L 1 ( X ) , L 2 ( X ) , L 3 ( X ) , L 4 ( X ) , L 5 ( X ) , L 6 ( X ) ) , where
L 1 ( X ) = Λ δ S α S β S I L 2 ( X ) = β S I ϕ A + σ β V I δ A , L 3 ( X ) = α S σ β V I δ V L 4 ( X ) = ξ ϕ A ( δ 2 + ρ ) I L 5 ( X ) = ( 1 ξ ) ϕ A ( δ 1 + γ ) C L 6 ( X ) = ρ I + γ C δ R
For any X , X ¯ Ω :
L ( X ) L ( X ¯ ) = L 1 ( X ) L 1 ( X ¯ ) + L 2 ( X ) L 2 ( X ¯ ) + L 3 ( X ) L 3 ( X ¯ ) + L 4 ( X ) L 4 ( X ¯ ) + L 5 ( X ) L 5 ( X ¯ ) + L 6 ( X ) L 6 ( X ¯ ) = ( δ + α ) ( S S ¯ ) β ( I S I ¯ S ¯ ) + β ( I S I ¯ S ¯ ) ( ϕ + δ ) ( A A ¯ ) + σ β ( V I V ¯ I ¯ ) + | ξ ϕ ( A A ¯ ) + ( δ 2 + ρ ) ( I I ¯ ) | + | ( 1 ξ ) ϕ ( A A ¯ ) ( γ + δ 1 ) ( C C ¯ ) | + ρ ( I I ¯ ) + γ ( C C ¯ ) δ ( R R ¯ ) + α ( S S ¯ ) σ β ( V I V ¯ I ¯ ) δ ( V V ¯ ) ( δ + 2 α ) S S ¯ + 2 β I S S ¯ I ¯ + 2 σ β V I I ¯ V ¯ + ( δ + 2 ϕ ) A A ¯ + ( 2 ρ + δ ) I I ¯ + ( δ 1 + 2 γ ) C C ¯ + δ V V ¯ + δ R R ¯ ( δ + 2 α ) | S S ¯ | + 2 β | I S S ¯ I + S ¯ I S ¯ I ¯ | + 2 σ β | V I V ¯ I + V ¯ I I ¯ V ¯ | + ( δ + 2 ϕ ) | A A ¯ | + ( 2 ρ + δ ) | I I ¯ | + δ 1 + 2 γ | C C ¯ | + δ | V V ¯ | + δ | R R ¯ | ( δ + 2 α ) | S S ¯ | + 2 β M ( | I I ¯ | + | S S ¯ | ) + 2 δ β ( | V V ¯ | + | I I ¯ | ) + ( δ + 2 ϕ ) | A A ¯ | + ( 2 ρ + δ ) | I I ¯ | + δ 1 + 2 γ | C C ¯ | + δ | V V ¯ | + δ | R R ¯ | F 1 | S S ¯ | + F 2 | A A ¯ | + F 3 | V V ¯ | + F 4 | I I ¯ | + F 5 | C C ¯ | + F 6 | R R ¯ | F || X X ¯ || ,  where  F = max F 1 , F 2 , F 3 , F 4 , F 5 , F 6 ,
and
F 1 = ( δ + 2 β M + 2 α ) F 2 = ( δ + 2 ϕ ) F 3 = 2 β M + 2 σ β M + δ 2 + 2 ρ F 4 = δ 1 + 2 γ F 5 = { 2 σ β M + δ } F 6 = δ
Hence Lipschitz’s condition is satisfied. Therefore, Lemma 2 confirms that the system (2) posses a unique solution of system (2). □
Proof of Theorem 5. 
For disease-free equilibrium,
J ( E 0 * ) = ( δ + α ) 0 0 β Λ δ + α 0 0 ( δ + ϕ ) 0 β Λ δ + α + σ β α Λ δ ( δ + α ) 0 α 0 δ σ β α Λ δ ( δ + α ) 0 0 ξ ϕ 0 ( ρ + δ 2 ) 0 0 ( 1 ξ ) ϕ 0 0 ( δ 1 + γ )
Characteristic equation is:
( λ + δ 2 + ρ ) ( λ + δ + α ) ( λ + δ ) P ( λ ) = 0
where
P ( λ ) = λ 2 + λ ( ϕ + δ + δ 2 + ρ ) + ( ϕ + δ ) ( δ 2 + ρ ) ξ ϕ Λ β δ + α + σ β α Λ δ ( δ + α ) .
The three eigenvalues are ( ρ + δ 2 ) , ( δ + α ) , δ . For these eigenvalues, we have arg ( λ ) = π > ε π 2 . The rest two eigenvalues can be found from:
λ 2 + λ ( ϕ + δ + δ 2 + ρ ) + ( ϕ + δ ) ( δ 2 + ρ ) ξ ϕ Λ β δ + α + σ β α Λ δ ( δ + α ) = 0 .
Discriminant of the quadratic Equation (A1) is
( ϕ + δ + δ 2 + ρ ) 2 4 ( ϕ + δ ) ( δ 2 + ρ ) ξ ϕ Λ β δ + α + σ β α Λ δ ( δ + α ) = ( ϕ + δ δ 2 ρ ) 2 + 4 ξ ϕ Λ β δ + α + σ β α Λ δ ( δ + α ) 0
So, all the roots of the quadratic Equation (15) are real. The both roots will be negative if
( ϕ + δ ) ( δ 2 + ρ ) ξ ϕ Λ β δ + α + σ β α Λ δ ( δ + α ) > 0 ξ β Λ ϕ ( δ + α ) ( δ 2 + ρ ) ( δ + ϕ ) α ξ ϕ β σ Λ δ ( ϕ + δ ) ( δ + α ) ( δ 2 + ρ ) + 1 > 0 R 0 + 1 > 0 R 0 < 1
Using Theorem 3, the asymptotic stability of E 0 can be reached. □
Proof of Theorem 6. 
For endemic equilibrium point,
J ( S * , A * , V * , I * , C * ) = m 11 m 12 m 13 m 14 m 15 m 21 m 22 m 23 m 24 m 25 m 31 m 32 m 33 b 34 b 35 m 41 m 42 m 43 m 44 m 45 m 51 m 52 m 53 m 54 m 55
where
m 11 = ( δ + α ) β I * m 12 = 0 m 13 = 0 m 14 = β S * m 15 = 0 m 21 = β I * m 22 = ( δ + ϕ ) m 23 = σ β I * m 24 = β S * + σ β V * m 25 = 0 m 31 = α m 32 = 0 m 33 = ( σ β I * + δ ) m 34 = σ β V * m 35 = 0 m 41 = 0 m 42 = ξ ϕ m 43 = 0 m 44 = ( δ 2 + ρ ) m 45 = 0 m 51 = 0 m 52 = ( 1 ξ ) ϕ m 53 = 0 m 54 = 0 m 55 = ( δ 1 + γ )
Characteristic equation is:
( λ + δ 1 + γ ) Q ( λ ) = 0 ,
where Q ( λ ) = λ 4 + d 1 λ 3 + d 2 λ 2 + d 3 λ + d 4 and
d 1 = ( m 11 + m 12 + m 13 + m 14 ) d 2 = m 12 m 21 m 13 m 31 m 23 + m 32 + m 22 m 33 m 14 m 41 m 24 m 42 m 34 m 43 + ( m 22 + m 33 ) m 44 + m 11 ( m 22 + m 33 + m 44 ) d 3 = m 11 ( m 23 m 32 m 22 m 33 + m 24 m 42 + m 34 m 43 ) + m 14 ( m 22 m 41 + m 33 m 41 m 21 m 42 m 31 m 43 ) + m 24 ( m 33 m 42 m 32 m 43 ) + m 34 ( m 22 m 43 m 23 m 42 ) m 44 ( m 11 m 22 m 23 m 32 + ( m 11 + m 22 ) m 33 ) + m 13 ( m 21 m 32 m 34 m 41 + m 31 ( m 22 + m 44 ) ) + m 12 ( m 23 m 31 m 24 m 41 + m 21 ( m 33 + m 44 ) ) d 4 = m 12 m 41 ( m 24 m 33 m 23 m 34 ) + m 11 m 42 ( m 23 m 34 m 24 m 33 ) + m 24 m 43 ( m 32 m 11 m 12 m 31 ) + m 34 m 43 ( m 12 m 21 m 11 m 22 ) + m 14 { m 23 ( m 32 m 41 m 31 m 42 ) + m 21 ( m 33 m 42 m 32 m 43 ) + m 22 ( m 31 m 43 m 33 m 41 ) } + { m 12 ( m 23 m 31 m 21 m 33 ) + m 11 ( m 22 m 33 m 32 m 23 ) } + m 13 { m 22 ( m 34 m 41 m 31 m 44 ) + m 24 ( m 31 m 42 m 32 m 41 ) + m 21 ( m 32 m 44 m 34 m 42 ) }
One root of the characteristic Equation (A3) is ( δ 1 + γ ) and for this eigenvalue we have arg ( λ ) = π > ε π 2 . Other four roots can be found from Q ( λ ) = 0 .
The system (14) is asymptotically stable around endemic equilibrium E 1 * if the following conditions fulfilled:
  • If ( Q ) < 0 , d 1 > 0 , d 2 > 0 , d 3 > 0 , d 4 > 0 and ε < 1 3
  • If ( Q ) < 0 , d 1 > 0 , d 2 > 0 , d 3 > 0 , d 4 > 0 and d 2 = d 1 d 4 d 3 + d 3 d 1 , for ε ( 0 , 1 ) ,
where Q ( λ ) = λ 4 + d 1 λ 3 + d 2 λ 2 + d 3 λ + d 4 .
Now, d 4 > 0 is the essential condition for stability (asymptotically stable) of the endemic equilibrium point E 1 * (Lemma 4). From (29), we have
d 4 = β S * ( σ β I * α ξ ϕ ) β I * ( σ β I * + δ ) ξ ϕ ( δ + ϕ ) ( δ + α + β I * ) ( σ β I * + δ ) < 0 .
Hence, the endemic equilibrium is unstable. □

References

  1. Wikipedia Contributors. COVID-19 Pandemic in Westbengal. Wikipedia, The Free Encyclopedia. 3 August 2020, 15:22 UTC. Available online: https://https://en.wikipedia.org/wiki/COVID-19_pandemic_in_West_Bengal (accessed on 22 October 2020).
  2. ‘We Do Not Actually Have That Answer Yet’: WHO Clarifies Comments on Asymptomatic Spread of COVID-19. Available online: https://www.statnews.com/2020/06/09/who-comments-asymptomatic-spread-covid-19/ (accessed on 8 September 2021).
  3. Guan, W.; Liang, W. Comorbidity and its impact on 1590 patients with covid-19 in china: A nationwide analysis. Eur. Respir. J. 2020, 55, 2000547. [Google Scholar] [CrossRef] [Green Version]
  4. Gupta, R.; Hussain, A.; Misra, A. Diabetes and COVID-19: Evidence, current status and unanswered research questions. Eur. J. Clin. Nutr. 2020, 74, 864–870. [Google Scholar] [CrossRef]
  5. Lee, S.C.; Son, K.J.; Han, C.H.; Jung, J.Y.; Park, S.C. Impact of comorbid asthma on severity of coronavirus disease (COVID-19). Sci. Rep. 2020, 10, 21805. [Google Scholar] [CrossRef]
  6. Paramasivam, A.; Priyadharsini, J.V.; Raghun hakumar, S.; Elumalai, P. A novel COVID-19 and its effects on cardiovascular disease. Hypertens. Res. 2020, 43, 729–730. [Google Scholar] [CrossRef] [PubMed]
  7. Ghostine, R.; Gharamti, M.; Hassrouny, S.; Hoteit, I. An Extended SEIR Model with Vaccination for Forecasting the COVID-19 Pandemic in Saudi Arabia Using an Ensemble Kalman Filter. Mathematics 2021, 9, 636. [Google Scholar] [CrossRef]
  8. Podlubny, I. Fractional Differential Equations; Academic Press: San Diego, CA, USA, 1999. [Google Scholar]
  9. Du, M.; Wang, Z.; Hu, H. Measuring memory with the order of fractional derivative. Sci. Rep. 2013, 3, 3431. [Google Scholar] [CrossRef]
  10. Teodoro, G.S.; Machado, J.T.; de Oliveira, E.C. A review of definitions of fractional derivatives and other operators. J. Comput. Phys. 2019, 388, 195–208. [Google Scholar] [CrossRef]
  11. Podlubny, I. Geometric and Physical Interpretation of Fractional Integration and Fractional Differentiation. Fract. Calc. Appl. Anal. 2002, 5, 367–386. [Google Scholar]
  12. Diethelm, K. Efficient Solution of Multi-Term Fractional Differential Equations Using P(EC)mE Methods. Computing 2003, 71, 305–319. [Google Scholar] [CrossRef]
  13. Das, M.; Samanta, G.P. Optimal Control of Fractional Order COVID-19 Epidemic Spreading in Japan and India 2020. Biophys. Rev. Lett. 2020, 15, 207–236. [Google Scholar] [CrossRef]
  14. Das, M.; Samanta, G.P. Stability analysis of a fractional ordered COVID-19 model. Comput. Math. Biophys. 2021, 9, 22–45. [Google Scholar] [CrossRef]
  15. Das, M.; Samanta, G.P.; De la Sen, M. Stability Analysis and Optimal Control of a Fractional Order Synthetic Drugs Transmission Model. Mathematics 2021, 9, 703. [Google Scholar] [CrossRef]
  16. Das, M.; Samanta, G.P. Optimal control of a fractional order epidemic model with carriers. Int. J. Dynam. Control 2021, 1–22. [Google Scholar] [CrossRef]
  17. Ghosh, S.; Samanta, G.P.; Mubayi, A. Comparison of Regression Approaches for Analyzing Survival Data in the Presence of Competing Risks: An Application to COVID-19. Lett. Biomath. 2021, 8, 29–47. [Google Scholar]
  18. Ghosh, S.; Samanta, G.P.; Nieto, J.J. Application of non-parametric models for analyzing survival data of COVID-19 patients. J. Infect. Public Health 2021, 14, 1328–1333. [Google Scholar] [CrossRef]
  19. Khan, M.A.; Atangana, A. Modeling the dynamics of novel coronavirus (2019-nCov) with fractional derivative. Alex. Eng. J. 2020, 59, 2379–2389. [Google Scholar] [CrossRef]
  20. Liu, Z.; Magal, P.; Seyd, O.; Webb, G. A COVID-19 epidemic model with latency period. Infect. Dis. Model. 2020, 5, 323–337. [Google Scholar] [CrossRef]
  21. Rajagopal, K.; Hasanzadeh, N.; Parastesh, F.; Hamarash, I.I.; Jafari, S.; Hussain, I. A fractional-order model for the novel coronavirus (COVID-19) outbreak. Nonlinear Dyn. 2020, 101, 711–718. [Google Scholar] [CrossRef]
  22. Saha, S.; Samanta, G.P.; Nieto, J.J. Epidemic model of COVID-19 outbreak by inducing behavioural response in population. Nonlinear Dyn. 2020, 102, 455–487. [Google Scholar] [CrossRef]
  23. Saha, S.; Samanta, G.P. Modelling the role of optimal social distancing on disease prevalence of COVID-19 epidemic. Int. J. Dyn. Control 2021, 9, 1053–1077. [Google Scholar] [CrossRef]
  24. Shaikh, A.S.; Shaikh, I.N.; Nisar, K.S. A mathematical model of COVID-19 using fractional derivative: Outbreak in India with dynamics of transmission and control. Adv. Differ. Equ. 2020, 2020, 373. [Google Scholar] [CrossRef] [PubMed]
  25. Caputo, M. Linear models of dissipation whose Q is almost frequency independent. Ann. Geophys. 1966, 19, 383. [Google Scholar] [CrossRef]
  26. Caputo, M. Linear models of dissipation whose Q is almost frequency independent—II. Geophys. J. Int. 1967, 13, 529. [Google Scholar] [CrossRef]
  27. Petras, I. Fractional-Order Nonlinear Systems: Modeling Aanlysis and Simulation; Higher Education Press: Beijing, China, 2011. [Google Scholar]
  28. Diekmann, O.; Heesterbeek, J.A.P.; Roberts, M.G. The construction of next-generation matrices for compartmental epidemic models. J. R. Soc. Interface 2010, 7, 7873–7885. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  29. Van den Driessche, P.; Watmough, J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math. Biosci. 2002, 180, 29–48. [Google Scholar] [CrossRef]
  30. Odibat, Z.; Shawagfeh, N. Generalized Taylor’s formula. Appl. Math. Comput. 2007, 186, 286–293. [Google Scholar] [CrossRef]
  31. Liang, S.; Wu, R.; Chen, L. Laplace transform of fractional order differential equations. Electron. J. Differ. Equ. 2015, 2015, 139. [Google Scholar]
  32. Kexue, L.; Jigen, P. Laplace transform and fractional differential equations. Appl. Math. Lett. 2011, 24, 2019–2023. [Google Scholar] [CrossRef] [Green Version]
  33. Li, Y.; Chen, Y.; Podlubny, I. Stability of fractional-order nonlinear dynamic systems: Lyapunov direct method and generalized Mittag-Leffler stability. Comput. Math. Appl. 2010, 59, 1810–1821. [Google Scholar] [CrossRef] [Green Version]
  34. Haubold, H.J.; Mathai, A.M.; Saxena, R.K. Mittag-leffler functions and their applications. J. Appl. Math. 2011, 2011, 298628. [Google Scholar] [CrossRef] [Green Version]
  35. Ahmed, E.; El-Sayed, A.M.A.; El-Saka, H. On some Routh-Hurwitz conditions for fractional order differential equations and their applications in Lorenz, Rössler, Chua and Chen systems. Phys. Lett. A 2006, 358, 1–4. [Google Scholar] [CrossRef]
  36. Gelf, I.M.; Kapranov, M.M.; Zelevinsky, A.V. Discriminants, Resultants, and Multidimensional Determinants; Birkhäuser: Boston, MA, USA, 1994; ISBN 978-0-8176-3660-9. [Google Scholar]
  37. Matouk, A.E. Stability conditions, hyperchaos and control in a novel fractional order hyper-chaotic system. Phys. Lett. A 2009, 373, 2166–2173. [Google Scholar] [CrossRef]
  38. Zhao, H.; Mousseau, V.A. Extended Forward Sensitivity Analysis for Uncertainty Quantification. Nucl. Technol. 2013, 181, 184–195. [Google Scholar] [CrossRef] [Green Version]
  39. Li, Q.; Guan, X.; Wu, P.; Wang, X.; Zhou, L. Early transmission dynamics in wuhan, china, of novel coronavirus-infected pneumonia. N. Engl. J. Med. 2020, 382, 1199–1207. [Google Scholar] [CrossRef] [PubMed]
  40. Abel, T.; McQueen, D. The COVID-19 pandemic calls for spatial distancing and social closeness: Not for social distancing! Int. J. Public Health 2020, 65, 231. [Google Scholar] [CrossRef] [Green Version]
  41. Ding, Y.; Wang, Z.; Ye, H. Optimal Control of a Fractional-Order HIV Immune System with Memory. IEEE Trans. Control Syst. Technol. 2012, 20, 763–769. [Google Scholar] [CrossRef]
  42. Agarwal, O.P. A general formulation and solution scheme for fractional optimal control problems. Nonlinear Dyn. 2004, 38, 323–337. [Google Scholar] [CrossRef]
  43. Ndaïrou, F.; Torres, D.F.M. Distributed-Order Non-Local Optimal Control. Axioms 2020, 9, 124. [Google Scholar] [CrossRef]
  44. Kheiri, H.; Jafri, M. Optimal control of a fractional-order model for the HIV/AIDS epidemic. Int. J. Biomath. 2018, 11, 1850086. [Google Scholar] [CrossRef]
  45. Guo, T.L. The Necessary Conditions of Fractional Optimal Control in the Sense of Caputo. J. Opt. Theory Appl. 2013, 156, 115–126. [Google Scholar] [CrossRef]
  46. Garrappa, R. On linear stability of predictor-corrector algorithms for fractional differential equations. Int. J. Comput. Math. 2010, 87, 2281–2290. [Google Scholar] [CrossRef]
  47. WHO. Coronavirus Disease (COVID-19) Pandemic. Available online: https://www.who.int/emergencies/diseases/novel-coronavirus-2019 (accessed on 24 October 2021).
  48. COVID-19 Coronavirus Pandemic. 2021. Available online: https://www.worldometers.info/coronavirus/country/india/ (accessed on 24 October 2021).
Figure 1. The system 2 is represented schematically.
Figure 1. The system 2 is represented schematically.
Mathematics 09 02806 g001
Figure 2. Real Data vs. Model Data of total infected ( I + C ) cases between 1 June 2021 and 31 July 2021.
Figure 2. Real Data vs. Model Data of total infected ( I + C ) cases between 1 June 2021 and 31 July 2021.
Mathematics 09 02806 g002
Figure 3. Real Data vs. Model Data of total infected ( I + C ) cases between 1 June 2021 and 20 October 2021.
Figure 3. Real Data vs. Model Data of total infected ( I + C ) cases between 1 June 2021 and 20 October 2021.
Mathematics 09 02806 g003
Figure 4. Time series of state variables when no control is implemented.
Figure 4. Time series of state variables when no control is implemented.
Mathematics 09 02806 g004
Figure 5. Time series of state variables when combined control strategy is implemented.
Figure 5. Time series of state variables when combined control strategy is implemented.
Mathematics 09 02806 g005
Figure 6. Time series of optimal control variables and optimal cost when combined control strategy is implemented.
Figure 6. Time series of optimal control variables and optimal cost when combined control strategy is implemented.
Mathematics 09 02806 g006
Figure 7. Time series of total infected population ( I + C ) when control is applied.
Figure 7. Time series of total infected population ( I + C ) when control is applied.
Mathematics 09 02806 g007
Figure 8. Time series of total infected population ( I + C ) and infected cases with co-morbidity when ε = 0.6 , 0.7 , 0.8 , 0.9 .
Figure 8. Time series of total infected population ( I + C ) and infected cases with co-morbidity when ε = 0.6 , 0.7 , 0.8 , 0.9 .
Mathematics 09 02806 g008
Figure 9. Variation of total infected population ( I + C ) and infected cases with co-morbidity when ε = 0.6 , 0.7 , 0.8 , 0.9 .
Figure 9. Variation of total infected population ( I + C ) and infected cases with co-morbidity when ε = 0.6 , 0.7 , 0.8 , 0.9 .
Mathematics 09 02806 g009
Figure 10. Variation of time series of all state variables for t = 300 days when ε = 0.7 , 0.8 , 0.9 .
Figure 10. Variation of time series of all state variables for t = 300 days when ε = 0.7 , 0.8 , 0.9 .
Mathematics 09 02806 g010
Figure 11. Time series of state variables when combined control strategy is implemented T f = 250 , ε = 0.7 .
Figure 11. Time series of state variables when combined control strategy is implemented T f = 250 , ε = 0.7 .
Mathematics 09 02806 g011
Figure 12. Time series of optimal control variables when combined control strategy is implemented T f = 250 , ε = 0.7 .
Figure 12. Time series of optimal control variables when combined control strategy is implemented T f = 250 , ε = 0.7 .
Mathematics 09 02806 g012
Figure 13. Time series of optimal cost and total infected population ( Φ ) when combined control strategy is implemented T f = 250 , ε = 0.7 .
Figure 13. Time series of optimal cost and total infected population ( Φ ) when combined control strategy is implemented T f = 250 , ε = 0.7 .
Mathematics 09 02806 g013
Table 1. Description of parameters used in the system 2.
Table 1. Description of parameters used in the system 2.
Λ Recruitment rate of S
σ Transmission rate at which vaccinated people become asymptomatic
β Transmission rate at which susceptible people become asymptomatic
δ Morbidity rate
δ 1 Disease-induced death rate of symptomatically infected with co-morbidity class
δ 2 Disease induced death rate of symptomatically infected without co-morbidity class
ρ Recovery rate of symptomatically infected without co-morbidity class ( I )
ϕ Conversion rate of asymptomatic infected people to symptomatic
ξ Fraction of asymptomatic population enters to infected class without co-morbidity ( I )
γ Recovery rate of symptomatically infected with co-morbidity class (C)
α Rate of vaccination
Table 2. The values of the parameters used in system (1) corresponding to the situation in West Bengal (India).
Table 2. The values of the parameters used in system (1) corresponding to the situation in West Bengal (India).
ParametersValuesSource
Λ 0.00004 [1]
β 0.3945 Estimated
δ 1 / ( 70.4 × 365 ) [1]
δ 1 1 / ( 60 × 365 ) Assumed
δ 2 1 / ( 65 × 365 ) Assumed
ϕ 1 / 5.2 [39]
ξ 0.3 Assumed
ρ 0.1245 Assumed
σ 0.08 Estimated
γ 0.1240 Assumed
α 3 × 10 4 Estimated
ε 0.7Estimated
Table 3. Normalized forward sensitivity index for different parameters correspondence to Table 2.
Table 3. Normalized forward sensitivity index for different parameters correspondence to Table 2.
ParametersSensitivity Index
β +1
ϕ 0.0104
ρ −0.9964
α −0.8704
Λ +1
σ +0.0152
Table 4. Day by day case report between 1 June 2021 and 20 October 2021 of West Bengal (India).
Table 4. Day by day case report between 1 June 2021 and 20 October 2021 of West Bengal (India).
DayTotal Infected CasesModel Value
1/6/202194249000
6/6/202176823385
11/6/20152742978
16/6/202132681925
21/6/202121841651
26/6/202119331444
1/07/202114781286
6/07/202113671163
11/07/20219971054
16/07/2021831962
21/07/2021752893
26/07/2021657832
31/07/2021769777
05/08/2021732728
10/08/2021639683
15/08/2021663648
20/08/2021673758
25/08/2021758685
30/08/2021652557
04/09/2021677590
09/09/2021680560
14/09/2021702540
19/09/2021619500
24/09/2021707470
30/09/2021651452
05/10/2021618420
10/10/2021619405
15/10/2021530393
20/10/2021857380
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Das, M.; Samanta, G.; De la Sen, M. A Fractional Ordered COVID-19 Model Incorporating Comorbidity and Vaccination. Mathematics 2021, 9, 2806. https://doi.org/10.3390/math9212806

AMA Style

Das M, Samanta G, De la Sen M. A Fractional Ordered COVID-19 Model Incorporating Comorbidity and Vaccination. Mathematics. 2021; 9(21):2806. https://doi.org/10.3390/math9212806

Chicago/Turabian Style

Das, Meghadri, Guruprasad Samanta, and Manuel De la Sen. 2021. "A Fractional Ordered COVID-19 Model Incorporating Comorbidity and Vaccination" Mathematics 9, no. 21: 2806. https://doi.org/10.3390/math9212806

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop