key: cord-1023717-3l3sy457
authors: Mahata, Animesh; Paul, Subrata; Mukherjee, Supriya; Das, Meghadri; Roy, Banamali
title: Dynamics of Caputo Fractional Order SEIRV Epidemic Model with Optimal Control and Stability Analysis
date: 2022-01-17
journal: Int J Appl Comput Math
DOI: 10.1007/s40819-021-01224-x
sha: 62f1d4764bc5570ffa416031f5f12748ce2c05ae
doc_id: 1023717
cord_uid: 3l3sy457

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 [Formula: see text] . 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.

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] . 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 andMotivation'. 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 andDiscussion', numerical simulation and discussion are presented via MATLAB. Finally, section 'Conclusion' includes conclusion of the paper.

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

The Caputo fractional derivative [29] [30] [31] [32] [33] [34] [35] of order 0 < ν ≤ 1 is defined as Definition 1 A function g : R + → R with fractional order 0 < ν ≤ 1, is defined as

Here the Gamma function is defined by ().

The Caputo fractional derivative operator of order 0 < ν ≤ 1, is defined as

The Riemann-Liouville fractional derivative of order 0 < ν ≤ 1, is defined as

Let g ∈ C[a, b], a < b, be a function, and 0 < ν ≤ 1. The fractional derivative in Caputo sense is defined as. where M(ν) denotes the normalization function with M(0) M(1) 1.

The Laplace transform for the fractional operator of order 0 < ν ≤ 1 is defined as.

For a 1 , a 2 ∈ R + and A ∈ C n×n where C denote complex plane, then.

Lemma 2 Let us the fractional order system as. 

calculated at the equilibrium points satisfy arg λ j > νπ 2 .

∈ R + is a differentiable function. Then, for any t > 0,

), g * ∈ R + , ∀ν ∈ (0, 1).

The entire population (N ) is divided into five categories, namely, the susceptible individuals (S), the exposed individuals (E), the infected individuals (I ), the recovered individuals (R) and the vaccinated individuals (V ) at any time t ≥ 0. Thus

The proposed SEIRV model with vaccination is depicted in Fig. 1 as a flow diagram. 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, δ 0 : the rate of mortality of all individuals, δ: the rate of vaccination, δ 1 : the rate of progression from exposed to infected individuals, δ 2 : the recovery rate of infected individuals. In this communication, we consider the SEIRV model using fractional order derivatives with Caputo operator of order 0 < ν < 1.

, 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 (time) −ν [38] . Next, let us consider t 0 0 and disregarded the super script ν of all parameters and the system becomes: The initial conditions are

Proposition For all t ≥ 0,the variables are non-negative.

The closed region {(S, E, I , R, V ) ∈ R 5 : 0 < N ≤ λ δ 0 }is positive invariant for the system (3.4) . Proof From the model (3.4) we have.

(3.5)

Using Laplace transform, we obtain

Taking inverse Laplace transform, we have

According to Mittag-Leffler function,

Hence, N (t)

And hence the model (3.4) is bounded above by λ δ 0 . Thus S, E, I , R, and V are all positive functions, and the system (3.4) is positive invariant.

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

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

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

Theorem 1 G 1 satisfies the Lipschitz condition and contraction if the following condition holds: 0 ≤ βa 1 + δ 0 + δ < 1.

Proof For S and S 1 ,

.

For G 1 , the Lipschitz condition is obtained, and if 0 ≤ βa 1 + δ 0 + δ < 1, then G 1 is a contraction.

In the same manner, G j , j 2, 3, 4, 5 satisfy the Lipschitz condition as follows:

, 

For j 2, 3, 4, 5, we obtain 0 ≤ K j < 1, then G j are contractions. Consider the following recursive patterns, as suggested by system (3.9):

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

In similar aspect, we obtain

As a result, we may write

.

Proof 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

We get the result by repeating the process.

Similarly, we may establish that P jn (t) → 0, j 2, 3, 4, 5 as n → ∞.

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

Then

Taking norm we have

From Lipschitz condition (3.10),

Thus

Proof Assuming that condition (3.14) is valid,

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

We have E 0 ( δ 0 +δ , 0, 0, 0, δ δ 0 (δ 0 +δ) ) and

The basic reproduction number, indicated by 0 , 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 ( 0 ) can be obtained from the leading eigen value of the matrix F V −1 where,

Therefore, the reproduction number 0 β δ 1 (δ 0 + δ)(δ 0 + δ 1 )(δ 0 + δ 2 )

. (3.16) Each parameter is obviously dependent on ν. So 0 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 0 .

This section shows the impact of altering parameter values on the 0 , 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 0 's sensitivity index towards the parameters , β, δ, δ 0 , δ 1 , δ 2 :

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

We have

As a result, we observed that E , E β and E δ 1 are positive. This means that increasing the values of the parameters , β and δ 1 raises the value of the fundamental reproduction number 0 . 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.

For stability analysis, the Jacobianmatrix of the system (3.4) at disease-free equilibrium point with P 11 −(δ + δ 0 ),P 13 −β S 0 ,P 22 −(δ 1 + δ 0 ),P 23 β S 0 , P 32 δ 1 ,P 33 −(δ 0 + δ 2 ), P 44 −δ 0 , P 55 −δ 0 , P 43 δ 2 .

Proof The characteristic equation of J 0 is given by determinant (P − λI 5 ) 0.

The roots of the characteristic equation are −δ 0 , −δ 0 , −(δ + δ 0 ) and λ 2 − λ(P 22 + P 33 ) + (P 22 P 33 − P 23 P 32 ) 0.

By Routh-Hurwitz Criterion, the roots are negative if (P 22 + P 33 ) < 0 and (P 22 P 33 − P 23 P 32 ) > 0. Now (P 22 P 33 − P 23 P 32 ) > 0

Since the first three roots are negative and other roots will be negative if 0 < 1 and positive if 0 > 1.

Therefore the equilibrium point E 0 is locally asymptotically stable or unstable according as 0 < 1 or 0 > 1.

Proof Taking the appropriate Lyapunov function into consideration.

The time fractional derivative of the above function is

From (3.4) we get,

Hence if 0 < 1, then C D ν t F(t) < 0. Hence, by LaSalle's extension to Lyapunov's principle [41, 42] , the disease-free equilibrium point E 0 is globally asymptotically stable and unstable if 0 > 1.

If 0 > 1, the epidemic equilibrium E 1 (S * , E * , I * , R * , V * )is locally asymptotically stable.

Since A > 0, C > 0, AB > C, by Routh-Hurwitz Criterion, the system (3.4) is locally asymptotically stable at E 1 .

The epidemic equilibrium E 1 is globally asymptotically stable if 0 > 1.

Proof 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,

(4.1)

Using system (3.4) we get,

We have Eq. (3.4) in steady state,

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, C D ν t W (t) ≤ 0 for 0 > 1. Hence W is a Lyapunov function. If 0 > 1, the epidemic equilibrium E 1 is globally asymptotically stable, according to LaSalle's Invariance Principle [42] .

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 w(0 ≤ w(t) ≤ 1) and to identify the best control w * to minimize the cost function J (w) of the control strategy.

The cost function J w * min(J (w(t)))with J (w)

Theorem 8 Let w(t)be a measurable control function on 0, t g , with w(t)having a value in [0, 1] . Then an optimal control w * minimizing the objective function J (w)of (5.1) with.

is the corresponding solution of the system (5.2).

Proof The Hamiltonian has been analyzed in the following manner:

with ε i (t), i 1, 2, 3, 4, 5 are the adjoint variables with ε i t g 0, expressed in terms canonical equations:

As a result, the issue of determining w * that minimizes H 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 w * , take into account the control constraints as well as the sign of ∂H ∂w .

Predictor formulae: 

. and j,n+1

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. 0.82. The comparison of the number of susceptible, infected, exposed and recovered individuals in case of δ 0 and δ 0.01 is quite obvious.The number of susceptible individualsis more for δ 0 than δ 0.01. 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 δ 0.01 than in case of δ 0. Figure 3 depicts the dynamical behavior of all individuals with a vaccination rate of 0.01 at fractional orders of ν 0.8, 0.9, 1 and the value of 0 is 1.55. 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. 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. 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 0 . 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 t g 100. Figure 5 depicts the time series of optimal control variable w * and optimal cost J * .

The data fitting and model validation of the system Table 1 209500000 [4]. The parametric values are given in Table 2 . We have taken t 1 day as time unit and t 100 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 ν 0.8, 0.85, 0.9.

In this paper we have discussed the optimal control of fractional order SEIRV model with vaccination as the control parameter w. Based on the COVID-19 cases data in India, collected upto 10th August, 2021, we estimated the basic reproduction number Time series of optimal control variable w * and optimal cost J * with parameter values corresponding to Table 1 0 without vaccination 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 ν 0.82 has been presented in Fig. 2 in case of δ 0 and δ 0.01. The stability analysis of the model shows that the system is locally as well as globally asymptotically stable at disease-free equilibrium point E 0 when 0 < 1 and at epidemic equilibrium E 1 when 0 > 1. Sensitivity analysis shows that 0 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 δ 1 , 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 w is employed then the transmission of the disease may be checked and eradicated. Additionally, the optimal control value w * has been determined in Theorem 8 to minimize the cost of vaccination given by J (w) t g 0 E + I + A 1 w 2 dt. We have assumed a final time t g 100 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 ν 0.85 with realistic data.

Funding None.

Data Availability There is nothing left after reporting all of the data needed for numerical simulations and comparisons in the tables and depicting it in the graphical representations.

Conflict of interest There are no conflicts of interest declared by the authors.

Ethical approval All of the authors participated directly and actively involved in the considerable effort that led to the publication of the article, and they will be held accountable for its content. This manuscript is the authors' original work, and it has never been published earlier. The manuscript is not being considered for publication anywhere at this time. The work accurately and completely represents the authors' own research and analysis.

The authors state that they have no known conflicting financial interests or personal connections that might have influenced the research presented in this publication.

Cultivation of viruses from a high proportion of patients with colds

A novel coronavirus outbreak of global health concern

Bernoulli was ahead of modern epidemiology

Contribution to mathematical theory of epidemics

An introduction to compartmental modeling for the budding infectious disease modeler

Solution of the epidemic model by Adomian decomposition method

COVID-19 infection: origin, transmission, and characteristics of human coronaviruses

Mathematical model of COVID-19 with comorbitity and controlling using non-pharmaceutical interventions and vaccination

Dynamics of COVID-19 transmission with comorbidity:a data driven modelling based approach

A study of behaviour for immune and tumor cells in immune genetic tumour model with non-singular fractional derivative

A fractional model for the dynamics of tuberculosis infection using Caputo-Fabrizio derivative

Transmission dynamics of fractional order typhoid fever model using Caputo-Fabrizio operator

A numeric analytic method for approximating a giving up smoking model containing fractional derivatives

Mathematical analysis of HIV/AIDS infection model with Caputo-Fabrizio fractional derivative

Estimation of the transmission risk of the 2019-nCoV and its implication for public health interventions

Approximate solutions of a generalized Hirota-Satsuma coupled KdV and a coupled mKdV systems with time fractional derivatives

Two analytical methods for time-fractional nonlinear coupled Boussinesq-Burger's equations arise in propagation of shallow water waves

Impulsive initial value problems for a class of implicit fractional differential equations

An iterative method for solving nonlinear functional equations

On linear stability of predictor-corrector algorithms for fractional differential equations

A prey-predator fractional order model with fear effect and group defense

The fractional dynamics of a linear triatomic molecule

Analysis of a basic SEIRA model with Atangana-Baleanu derivative

On a nonlinear dynamical system with both chaotic and non-chaotic behaviours: a new fractional analysis and control

A nonstandard finite difference scheme for the modelling and nonidentical synchronization of a novel fractional chaotic system

Hyperchaotic behaviours, optimal control, and synchronization of a nonautonomous cardiac conduction system

Theory and Applications of Fractional Differential Equations

Properties of the new fractional derivative without singular kernel

A new definition of fractional derivative without singular kernel

Laplace transform and fractional differential equations

Fractional-Order Nonlinear Systems: Modeling Aanlysis and Simulation

Laplace transform of fractional order differential equations. Electron

Stability of fractional-order nonlinear dynamic systems: Lyapunov direct method and generalized Mittag-Leffler stability

Stability and Hopf bifurcation analysis of an SVEIR epidemic model with vaccination and multiple time delay

Modeling the virus dynamics in computer network with SVEIR model and nonlinear incident rate

A commentary on fractionalization ofmulti-compartmental models

Dynamics analysis and optimal control strategy for a SIRS epidemic model with two discrete time delays

Estimating the size of COVID-19 epidemic outbreak

Global dynamics of an SEIR epidemic model with vertical transmission

Differential Equations and Dynamical Systems

Optimal control of a fractional-order hiv immune system with memory

A general formulation and solution scheme for fractional optimal control problems. Nonlinear Dyn

Pontryagin maximum principle for fractional ordinary optimal control problems

SEIR epidemic model and scenario analysis of COVID-19 pandemic

Optimal control applied to vaccination and treatment strategies for various epidemiological models

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations

The authors are grateful to the reviewers for their valuable comments and suggestions. 

As a result, we haveand w * max{min{w, 1}, 0} where w (ε 1 −ε 5 )S 2 A 1 . By substituting w * to the equation, the optimal condition may be determined for the system (5.2).

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 F r j0 is the arbitrary real number, ν > 0 and the fractional differential operator D ν t 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 h T m , t n nh, n 0, 1, 2, . . . ,m. Corrector formulae: