Dynamics of Caputo Fractional Order SEIRV Epidemic Model with Optimal Control and Stability Analysis.

In mid-March 2020, the World Health Organization declared COVID-19, a worldwide public health emergency. This paper presents a study of an SEIRV epidemic model with optimal control in the context of the Caputo fractional derivative of order 0 < ν ≤ 1 . The stability analysis of the model is performed. We also present an optimum control scheme for an SEIRV model. The real time data for India COVID-19 cases have been used to determine the parameters of the fractional order SEIRV model. The Adam-Bashforth-Moulton predictor-corrector method is implemented to solve the SEIRV model numerically. For analyzing COVID-19 transmission dynamics, the fractional order of the SEIRV model is found to be better than the integral order. Graphical demonstration and numerical simulations are presented using MATLAB (2018a) software.

Introduction

The first instances of corona virus infection in humans in 1965, with symptoms that were comparable to the common cold were first reported by Tyrrell and Bynoe [1]. From December 2019, a new coronavirus known as SARS-COV-2 has been identified in a number of nations, infecting thousands of individuals worldwide with a higher mortality rate. The virus, however, proved deadly in Wuhan, Hubei Province, China, in 2019, after multiple modifications [2]. Henceforth, the virus was named SARS-CoV-2 and the disease was declared to be COVID-19. The spread of the disease with fatal consequences soon accelerated in several countries and was declared a pandemic by WHO on 11th March 2020. When an infected individual coughs or sneezes, the virus is spread through contact with them or by respiratory droplets. It can harm the respiratory tract in a variety of ways, from a typical cold to more serious diseases like SARS [3]. COVID-19 has recently engulfed most of the world’s countries, infecting a huge number of individuals. COVID-19 has had a significant impact on people’s lives and has resulted in significant economic losses. Coronavirus (formally known as COVID-19) has killed 4,336,797 individuals and infected 205,468,881 individuals throughout the world since it first emerged in China (as on Aug 12, 2021) [4].

Epidemic mathematical models may help in various aspects, including analyzing the role of epidemic propagation and suggesting effective control measures. The idea that epidemic disease transmission patterns may be expressed mathematically dates back to 1766, when Daniel Bernoulli published an article describing the impact of smallpox variolization on lifespan [5]. The epidemiological model created by Kermack and Mckendrick during 1927 is an example of a mathematical model that can be useful in studying the development and prevention of infectious illness [6]. Differential equations of integral order are commonly used in epidemiology modelling of biological processes [7-11]. Many studies have been conducted in this field of mathematics, and it has been shown that differential equations using fractional operators are effective in demonstrating epidemic models linked to many infectious illnesses [12-16].Tang, Wang, Li, and Bragazzi [17] proposed a compartmental mathematical approach that would combine the disease’s clinical progression, the patient’s epidemiological state, and intervention measures. There are several concerns, however, about the integral order of the differential equations in such models. Fractional order differential equations, a relatively new and growing topic of mathematical calculus, takes into account such constraints. Various approaches for constructing actual and approximate solutions to fractional order differential equations have emerged throughout time as a result of extensive research [18-23]. Because of its memory effects, fractional-order models have been investigated for mimicking actual occurrences in recent decades [24]. In [25], the SEIRA mathematical model is analyzed using the Atangana–Baleanu fractional derivative method with the Mittag–Leffler kernel. Furthermore, because fractional models have been shown to be capable of accurately representing chaotic systems, they have surfaced in a variety of domains dealing with chaos, including physics, biology, and finances [26-28].

Research Background and Motivation

Fractional derivatives are a powerful tool for describing memory and heredity features in a wide range of systems and phenomena. Fractional-order differential equations preserve the function’s fundamental information in stacked form. We have such an additional variable (order of the derivative) in fractional-order modeling that is beneficial for numerical methods. 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. When talking about real problems, the Caputo derivative is highly useful since it allows traditional starting and boundary conditions be included in the derivation, and the derivative of a constant is zero, that is not the case with the Riemann–Liouville fractional derivative.

In this study, a fractional SEIRV model with optimum control has been created using Caputo fractional-order differential equations, motivated by the aforementioned studies and the benefits of Caputo fractional-order differential equations.

The objectives of this work are:

Investigate the SEIRV model’s dynamical behavior and stability.

Determine the Basic Reproduction number and Equilibrium points.

The model system is subjected to an optimal control analysis by controlling ‘vaccination rate’ parameter.

Application of the Adam-Bashforth-Moulton predictor–corrector technique to obtain numerical solution.

The paper is organized as follows: The fractional operator is defined and results for the fractional operatorand the Laplace transform are provided in see section ‘Research Background and Motivation’. SEIRV model with fractional derivative in Caputo sense and the existence and uniqueness of the model solution, including positivity and boundedness are established in see section ‘Preliminaries of Fractional Calculus’. The stability analysis and stability criterion of the model system are discussed in see section ‘Stability Analysis’. We also provide an optimal control strategy for an SEIRV model using the control parameter “vaccination rate” in see section ‘SEIRV Model with Optimal Control’. In section ‘Adam-Bashforth-Moulton Predictor–Corrector Scheme for the SEIRV Model’, we perform Adam-Bashforth-Moulton predictor–corrector scheme for the SEIRV model. In section ‘Numerical Simulation and Discussion’, numerical simulation and discussion are presented via MATLAB. Finally, section ‘Conclusion’ includes conclusion of the paper.

Preliminaries of Fractional Calculus

In this section we use some fundamental definitions of fractional differential and integral operators.

The Caputo fractional derivative [29-35] of order is defined as

Definition 1

A function with fractional order , is defined as

Here the Gamma function is defined by .

Definition 2

The Caputo fractional derivative operator of order , is defined as

Definition 3

The Riemann–Liouville fractional derivative of order , is defined as

Definition 4

Let be a function, and The fractional derivative in Caputo sense is defined as.where denotes the normalization function with .

Definition 5

The Laplace transform for the fractional operator of order is defined as.

Definition 6

For and where denote complex plane, then.

Lemma 1

(Generalized Mean Value Theorem)

Let and if is continuous in then .where

Lemma 2

Let us the fractional order system as.

with and . For calculate the equilibrium points, we have . These equilibrium points are locally asymptotically stable if and only if each eigen value of the Jacobian matrix calculated at the equilibrium points satisfy .

Lemma 3

Assume that is a differentiable function. Then, for any

Model Formulation

The entire population is divided into five categories, namely, the susceptible individuals , the exposed individuals , the infected individuals , the recovered individuals and the vaccinated individuals at any time . Thus

The proposed SEIRV model with vaccination is depicted in Fig. 1 as a flow diagram.

Fig. 1

The SEIRV model is represented schematically

The SEIRV model [36, 37] with vaccination in the sense of integral order is defined as follows on the basis of the flow diagram:where

: birth rate of susceptible individuals,

: the rate of infection of susceptible individuals,

: the rate of mortality of all individuals,

: the rate of vaccination,

: the rate of progression from exposed to infected individuals,

: the recovery rate of infected individuals.

In this communication, we consider the SEIRV model using fractional order derivatives with Caputo operator of order .

It is found that the model system’s time dimension (3.3) is correct as both sides of the model system’s equations have dimension [38]. Next, let us consider and disregarded the super script of all parameters and the system becomes:

The initial conditions are

Positivity and boundedness of Solutions

Proposition

For all the variables are non-negative.

The closed region is positive invariant for the system (3.4).

Proof

From the model (3.4) we have.

Using Laplace transform, we obtain

Taking inverse Laplace transform, we have

According to Mittag–Leffler function,

Hence, .

And hence the model (3.4) is bounded above by .

Thus and are all positive functions, and the system (3.4) is positive invariant.

Existence and uniqueness of solution

This section demonstrates that the system (3.4) has a unique solution. To begin with, we rewrite system (3.4) in the following manner:where

Taking integral transform on both sides of the Eqs. (3.8), we obtain

The kernels , i = 1, 2, 3, 4,5, fulfil the Lipschitz condition and contraction, as demonstrated.

Theorem 1

satisfies the Lipschitz condition and contraction if the following condition holds:

For and ,

Suppose , where is a bounded function.

For , the Lipschitz condition is obtained, and if 1, then is a contraction.

In the same manner, satisfy the Lipschitz condition as follows:where

For , we obtain , then are contractions. Consider the following recursive patterns, as suggested by system (3.9):with and .

Throughout the above system, we compute the norm of its first equation, and then

.

Possessing Lipschitz’s condition (3.10), we have

In similar aspect, we obtain

As a result, we may write

Theorem 2

A system of solutions described by the SEIRV model (3.4) exists if there exists such that

From (3.11) and (3.12), we have.

Thus, the system is continuous and has a solution. Now we’ll explain how the functions listed above may be used to construct a model solution (3.9). We make the assumption that

So

We get the result by repeating the process.

as .

Similarly, we may establish that as .

To examine the uniqueness of the solution, we assume that there is another solution of the system, such as .

Then

Taking norm we have

From Lipschitz condition (3.10),

Thus

Theorem 3

The model system ( 3.4 ) has a unique solution, provided that

Assuming that condition (3.14) is valid,

Then . So, we have . Similarly, we can prove that , , , .

Equilibrium Points and Basic Reproduction Number

The disease-free equilibrium points and the epidemic equilibrium point of the system (3.4) are obtained from

We have  = () and ,

where  = , ,  = ,  = ,  = .

The basic reproduction number, indicated by , is the estimated number of secondary cases generated by infection of a single susceptible individual.

Using next generation matrix method [39, 40], the reproduction number () can be obtained from the leading eigen value of the matrix where,

Each parameter is obviously dependent on So is a function of For analysis purpose, we have fixed the value of . If we change the value of , then all other parametric values will be changed and this will change the value of .

Significance of Sensitivity Parameters

This section shows the impact of altering parameter values on the , reproduction number’s perceived usefulness. The crucial parameter, which might be a critical threshold for illness management, must be identified.

The following are the mathematical representations of ’s sensitivity index towards the parameters , :

It is inferred that some derivatives seem positive, and that the basic reproductive number increases as any of the aforementioned positive value parameters is increased. The proportionate reaction to the proportion stimulation is used to calculate the elasticity.

We have

As a result, we observed that and are positive. This means that increasing the values of the parameters and raises the value of the fundamental reproduction number . The fundamental reproduction number might vary a lot depending on how these factors are changed. A highly sensitive component should be computed with care, since even slight variations might result in significant quantitative systemic changes.

Stability Analysis

For stability analysis, the Jacobianmatrix of the system (3.4) at disease-free equilibrium point is given by wherewith ,,,,.

Theorem 4

When When  < 1, the disease-free equilibrium point of the system ( 3.4 ) is locally asymptotically stable, and when , it is unstable.

Proof

The characteristic equation of is given by determinant

The roots of the characteristic equation are , , and .

By Routh-Hurwitz Criterion, the roots are negative if and

Now

Since the first three roots are negative and other roots will be negative if and positive if .

Therefore the equilibrium point is locally asymptotically stable or unstable according as or .

Theorem 5

When  < 1, the disease-free equilibrium point of the system (3.4) is globally asymptotically stable, and unstable when .

Taking the appropriate Lyapunov function into consideration.

The time fractional derivative of the above function is

From (3.4) we get,

Now,

Since , it follows that

Hence if then . Hence, by LaSalle's extension to Lyapunov’s principle [41, 42], the disease-free equilibrium point is globally asymptotically stable and unstable if .

Theorem 6

If , the epidemic equilibrium is locally asymptotically stable.

The characteristic equation is (. where

with

Since , by Routh-Hurwitz Criterion, the system (3.4) is locally asymptotically stable at .

Theorem 7

The epidemic equilibrium is globally asymptotically stable if .

The Goh-Volterra form's non-linear Lyapunov function is defined as.

Using Lemma 3 and then the above function’s time fractional derivative is,

Using system (3.4) we get,

We have Eq. (3.4) in steady state,

Substituting Eq. (4.3) into (4.2) we have

Further simplification gives,

Adding all infected classes without a single star (*) from (4.4) to zero:

The steady state was somewhat perturbed between (3.4) and (4.5), yielding:

Substituting the expression from (4.6) into (4.4) gives:

.

We have an arithmetic mean that is greater than the geometric mean.

Hecne,

for .

Hence is a Lyapunov function. If , the epidemic equilibrium is globally asymptotically stable, according to LaSalle's Invariance Principle [42].

SEIRV Model with Optimal Control

Vaccination is an important tool in the fight against infectious illnesses. The vaccine against Covid-19 has recently been proven to be an effective means of stopping the disease’s transmission. Ding et al. [43] and Agarwal et al. [44] have contributed on optimum control theory in fractional calculus. Pontryagain’s maximal principle [45] strikes at the core of the concept of optimal control in fractional calculus. Our goal is to incorporate the effectiveness of vaccination through a control measure namely and to identify the best control to minimize the cost function of the control strategy.

subject to

where .

Theorem 8

Let be a measurable control function on , with having a value in [0,1]. Then an optimal control minimizing the objective function of (5.1) with.where

is the corresponding solution of the system (5.2).

The Hamiltonian has been analyzed in the following manner:with

, i = 1, 2, 3, 4, 5 are the adjoint variables with expressed in terms canonical equations:

As a result, the issue of determining that minimizes in the presence of (5.2) is recast as minimizing the Hamiltonian with regard to the control. We then establish the following optimum condition using the Pontryagin principle:

which may be solved using state and adjoint variables to yield

For the best control , take into account the control constraints as well as the sign of .

As a result, we haveand where .

By substituting to the equation, the optimal condition may be determined for the system (5.2).

Adam-Bashforth-Moulton Predictor–Corrector Scheme for the SEIRV Model

The Adams–Bashforth-Moulton strategy is the most widely used numerical technique for addressing any fractional order initial value problems.

Let’s look at the fractional differential equation below.where is the arbitrary real number, and the fractional differential operator is identical to the well-known Volterra integral equation in the Caputo sense.

Using Adam’s-Bashforth-Moulton predictor–corrector scheme, we explore the numerical solution of a fractional order SEIRV model with vaccination. The algorithm is described in the following manner.

Let

Corrector formulae:

Predictor formulae:

where

and

Numerical Simulation and Discussion

In this section, we perform rigorous numerical simulations to evaluate and verify the analytical results of our model system (3.4). Using mathematical software MATLAB (2018a version), we have employed Adam’s-Bashforth-Moulton predictor–corrector scheme to obtain numerical solution to the system (3.4).

We investigate numerical simulations of the model system (3.4) for India in the Caputo sense, using the parameters listed in Table 1. We estimate and anticipate the progression of the COVID-19 pandemic using recent Indian data up to the 10th of August 2021 [4]. In the India scenario, Table 1 is utilized for simulation. The following figures were produced to examine the behavior of the model (3.4) under various initial conditions.

Table 1

Estimated values of parameters for India:

Parameter	Value	References
Λ	0.0182	[46]
β	0.476	[47]
δ_1	0.071	[48]
δ_2	0.286	[47]
δ_0	0.0073	[46]
δ	0.01	Model to fit

Figure 2 demonstrates the dynamical behavior of all individuals for fractional order The comparison of the number of susceptible, infected, exposed and recovered individuals in case of and is quite obvious.The number of susceptible individualsis more for than Similar is the case with exposed individuals and infected individuals. However, in case of recoveredindividuals, it is just the opposite, due to obvious reasons. Now, the recovered individuals will be more in case of than in case of .

Fig. 2

Comparison of dynamical behaviour of all individuals with respect to time for fractional order , and

Figure 3 depicts the dynamical behavior of all individuals with a vaccination rate of 0.01 at fractional orders of and the value of is . The purpose of this study is to demonstrate the importance of the COVID-19 vaccination rate. When we enforced a vaccination rate, the basic reproduction number decreases.

Fig. 3

Dynamical behaviour of all individuals with respect to time with a vaccination rate, and fractional order

Figure 4a-e shows the time series of susceptible individuals, exposed individuals, infected individuals, recovered individuals and vaccinated individuals across a time period of [0,100] with optimal control taking fractional order ν = 1.

Fig. 4

With respect to time, the time series of the model system (3.3) corresponds to Table 1

Vaccination is a critical component in preventing people from COVID-19, and various ideas have been proposed in which vaccination rates are viewed as quite beneficial. As a result, the addition of the vaccination parameter decreases the reproduction number . For the simulation of the optimal control problem subject to the model (3.3) corresponding to Table 1 in the India scenario, we used a final time of . Figure 5 depicts the time series of optimal control variable and optimal cost .

Fig. 5

Time series of optimal control variable and optimal cost with parameter values corresponding to Table 1

Data Fitting and Model Validation

The data fitting and model validation of the system (3.4) for Infected population in Brazil are described in this section. From the 10th of April to the 19th of July, 2021, we compared the model values with the real scenario for Brazil. The total initial population of Brazil is around [4]. The parametric values are given in Table 2. We have taken  day as time unit and as final time. Table 3 recommends day wise Infected population from 10th April, 2021 to 19th July, 2021. Figure 6 depict time series solution of Infected population of the system (3.4) for Table 2 taking

Table 2

The estimated parametric values are as follows in Brazil

Parameter	Value	References
Λ	0.0187	[46]
β	0.32	[47]
δ_1	0.344	[48]
δ_2	0.041	[47]
δ_0	0.0063	[46]
δ	0.01	Model to fit

Table 3

Day wise Infected population of Brazil from 10th April, 2021 to 19th July, 2021

Day	Infected population
10/04/2021	1,269,000
20/04/2021	1,285,000
30/04/2021	1,270,000
10/05/2021	1,111,000
20/05/2021	1,068,000
30/05/2021	1,108,000
09/06/2021	1,128,000
19/06/2021	1,257,000
29/06/2021	1,227,000
09/07/2021	813,700
19/07/2021	825,000

Fig. 6

Time series solution of Infected population of the system (3.4) for Table 2 taking

Conclusion

In this paper we have discussed the optimal control of fractional order SEIRV model with vaccination as the control parameter . Based on the COVID-19 cases data in India, collected upto 10th August, 2021, we estimated the basic reproduction number to be 3.67 and with vaccination to be 1.55. The fractional-order derivatives are usually more suitable in modeling since the choice of the derivative order provides one more degree of freedom and this leads to better fit to the real time data with less error than the integer-order model. A comparison of the number of individuals in different compartments for has been presented in Fig. 2 in case of and . The stability analysis of the model shows that the system is locally as well as globally asymptotically stable at disease-free equilibrium point when  < 1 and at epidemic equilibrium when  > 1. Sensitivity analysis shows that is directly proportional to thebirth rate of susceptible individuals the rate of infection of susceptible individuals and the rate of progression from exposed to infected individuals , all of which may be controlled with the effective execution of vaccination drives. We have used the Pontryagin’s Maximum Principle to provide the necessary conditions needed for the existence of the optimal solution to the optimal control problem. Adam-Bashforth-Moulton predictor–corrector technique has been used to obtain numerical solutions to the system. Numerical simulations are presented using MATLAB to validate the efficacy and impact of the control parameter. It is evident that if the control measure is employed then the transmission of the disease may be checked and eradicated. Additionally, the optimal control value has been determined in Theorem 8 to minimize the cost of vaccination given by . We have assumed a final time for optimal control. 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 (Fig. 6). This demonstrates that the order of derivative is important in system simulation. We have comparatively studied the model values and real scenario of Brazil starts from 10th April 2021 and continues up to 100 days. It has been observed that our model fits with with realistic data.