Fractional modeling of COVID-19 epidemic model with harmonic mean type incidence rate

Sowwanee Jitsinchayakul; Rahat Zarin; Amir Khan; Abdullahi Yusuf; Gul Zaman; Usa Wannasingha Humphries; Tukur A. Sulaiman

doi:10.1515/phys-2021-0062

Article Open Access

Fractional modeling of COVID-19 epidemic model with harmonic mean type incidence rate

Sowwanee Jitsinchayakul , Rahat Zarin , Amir Khan , Abdullahi Yusuf , Gul Zaman , Usa Wannasingha Humphries and Tukur A. Sulaiman

Published/Copyright: November 19, 2021

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Open Physics Volume 19 Issue 1

Abstract

Coronavirus disease 2019 (COVID-19) is a disease caused by severe acute respiratory syndrome coronavirus 2 (SARS CoV-2). It was declared on March 11, 2020, by the World Health Organization as a pandemic disease. Regrettably, the spread of the virus and mortality due to COVID-19 have continued to increase daily. The study is performed using the Atangana–Baleanu–Caputo operator with a harmonic mean type incidence rate. The existence and uniqueness of the solutions of the fractional COVID-19 epidemic model have been developed using the fixed point theory approach. Along with stability analysis, all the basic properties of the given model are studied. To highlight the most sensitive parameter corresponding to the basic reproductive number, sensitivity analysis is taken into account. Simulations are conducted using the first-order convergent numerical approach to determine how parameter changes influence the system’s dynamic behavior.

Keywords: epidemic model; stability analysis; numerical simulations; sensitivity analysis; COVID-19

1 Introduction

Bacteria, viruses, fungi, and parasites give rise to infectious diseases which are disorders caused by these organisms. As these organisms live in the human body and some of them live on the body and are the key factor for causing the disease, people are affected worldwide. China receives worldwide notoriety when an unknown virus attacks people and begins to destroy them. The agent of causation was later described as a new ailment and is known as coronavirus (severe acute respiratory syndrome coronavirus 2 [SARS CoV-2]). The authority in China has done the best to mitigate the outbreak, but sadly the virus kept spreading to the international domain. The dormancy period for coronavirus disease 2019 (COVID-19) is 11–14 days, while the common symptoms are shortness of breath, myalgia, coughing, and sneezing. Patients of comorbid conditions, and aged individuals are particularly susceptible to COVID-19. The COVID-19 transmission rate is very high and between 2.2 and 3.58 is the basic reproduction number. That is why it has spread across the globe and affected 213 nations. Therefore, on January 30, 2020, the WHO named COVID-19 as one of the global pandemic [1].

At first, the countries that are vulnerable to COVID-19 besides China were Iran and Italy. Iran shares border with Pakistan, and thousands of people come from Pakistan to Iran yearly for religious ceremonies. Some of these individuals were infected while coming back from Iran and they became the cause of the Pakistan spread of the ailment. In spite of the decision of the government to impose border closure, the case for the first time was announced by Karachi on February 26, 2020. The government employed quarantine measures to prevent the spread and casualties [1].

The introduction of lockdown to maintain social distance is to monitor COVID-19. In order to monitor the spread of the disease, this technique is an excellent measure. In the near future, a full lockdown could be the source of a major financial crisis. Specifically, lockdown in high-density nations may alleviate the rate of disease transmission, although it can be impossible to achieve maximum control. In order to keep a country’s economic status alive, there should be an acceptable balance between the two distinct features of safe free conditions and an absolute lockdown of government policies.

Mathematical modeling is regarded as an efficient method to explain the complex behavior of diseases [2,3, 4,5,6]. The mathematical modeling methods of the last century are also used by mathematicians. More recently, considerable attention has been given to the proposal of mathematical models in comprehending the ailment of infectious nature. Many researchers have devised models for the realization and regulation of the outbreak of transmissible diseases in a population. Infectious diseases are the second largest cause of death across the globe. For the study of transmissible ailment, mathematical models have been employed. Over the last few years, several researchers have been exploring infectious diseases and their mechanisms using different methods [7,8,9, 10,11,12]. This does not only help to control/spreading of infectious diseases but also aid in everyday life to prevent these diseases. Several researchers have researched epidemic models to examine and monitor various diseases, such as avian influenza, hepatitis B, tuberculosis, and leishmaniasis. The epidemic model is useful for both academia and everyday life [13,14, 15,16,17,18].

Because the existence and annihilation of COVID-19 are subject to numerous parameters of the affected system, we can not characterize the entire disease system throughout the globe using a single model. Models can be more detailed with the inclusion of detail and complexity, but this also makes mathematics very complex. The model is quite useful for appropriate care and formal vaccination. Motivated by this, we propose a mathematical model in this study that incorporates a quarantine class and government intervention steps such as lockdown, social distance coverage in the media, and enhancement of hygiene in public to reduce the transmission of the ailment.

Fractional calculus is the generalization of classical calculus. Mathematical models with noninteger order operators provide a better understanding of phenomena. Furthermore, models with fractional-order derivatives are capable capturing the fading memory and crossover behavior and provide a greater degree of accuracy. Mathematical models with fractional derivative give more insights about a disease under consideration [19,20,21, 22,23,24,25,26]. Different fractional operators with singular and nonsingular kernels were suggested in the literature [27,28, 29,30]. The applications of these fractional operators can be found in recent literature and references therein [31,32, 33,34]. Recently, very few COVID-19 models based on fractional-order operators are proposed. The authors in refs [35,36] considered a mathematical model for COVID-19 and cancer and hepatitis co-dynamics in fractional derivative and examined its results. Motivated by the above discussion, currently, in this work, we study the dynamical analysis of the COVID-19 model presented in ref. [37], considering the Atangana–Baleanu–Caputo (ABC) operator in order to gain more insights about the pandemic. The impact of important model parameters is shown for various values of the arbitrary fractional order of the ABC operator.

The rest of the article is organized as follows: The model formulation to the fractional-order derivative is presented in Section 2. Section 3 contains the basic reproduction number and local stability for the proposed model at corona-free equilibrium (CFE) points of fractional order derivative. Section 4 contains the existence and uniqueness of solutions for the ABC model. The Hyers–Ulam stability, determination of the important parameters, and their impact on the basic reproduction number through sensitivity analysis are presented in Section 5. The graphical results through biological discussion are depicted in Section 6. In Section 7, the work is summarized.

2 Model formulation

Mandal et al. [37] proposed the SEQIR epidemic model in which the total COVID-19 population N ( t ) is partitioned into the following sub-populations. At any time instant t , the human populations are subdivided into five time-dependent classes, namely, susceptible S ( t ) , exposed E ( t ) , hospitalized infected I ( t ) , quarantine Q ( t ) , and recovered or removed R ( t ) . The model is expressed as:

(2.1) S ( t ) = A − β ( 1 − ρ 1 ) ( 1 − ρ 2 ) W ( S , E ) + b 1 Q − d S − p S M , E ( t ) = β ( 1 − ρ 1 ) ( 1 − ρ 2 ) W ( S , E ) − b 2 E − α E − σ E − d E , Q ( t ) = b 2 E − b Q − c Q − d Q , I ( t ) = α E + c Q − ( η + d + δ ) I , R ( t ) = η I + σ E − d R + p S M .

Model parameters and their interpretation.

Parameters	Interpretation
A	The constant recruitment rate to the susceptible population
β	Disease transmission rate
ρ 1	The virus of COVID-19 is spreading when a vulnerable person comes into contact with an exposed person; therefore, we think that ρ 1 ( 0 < ρ 1 < 1 ) portion of susceptible human would maintain proper precaution measure
ρ 2	ρ 2 ( 0 < ρ 2 < 1 ) portion of the exposed class would take proper precaution measure for disease transmission (i.e., use of face mask, social distancing, and implementing hygiene) Therefore, the disease can only be transmitted to the ( 1 − ρ 1 ) S portion of susceptible individuals due to the contact of ( 1 − ρ 2 ) E portion of exposed individuals with a bi-linear disease transmission rate β
b 1	The individual becomes susceptible to the disease after the quarantine period
b 2	The portions of the exposed class go to the quarantine class
α	The portions of the exposed class go to the infected class
η	Recovery rate of the hospitalized infected populations I
σ	Recovery rate of the hospitalized exposed class E
d	Natural death rate
δ	COVID-19 induced death rate
M	The control policy parameter. To control the pandemic COVID-19, suitable governmental measures (like complete or semi lock-down, rationing system, continue media coverage on social isolation and improvement of public hygiene, home deliveries of essential commodities, to make an alternate source of income for job-losers during the lock-down, etc.) have been implemented by the various governmental and nongovernmental agencies. Thus, this policy may be considered as one of the effective control tools, and mainly the susceptible population of COVID-19 cases would be benefited due to this policy
p	Implementation rate of policy

In addition, W ( S , E ) is a feature known as the uptake function. The incidence of susceptible with exposed class is related to this feature. In this article, we provide the susceptible and exposed interaction, as in ref. [38], which has the form:

(2.2) W ( S , E ) ≔ 2 S E S + E .

In addition, the function described in (2.2) displays the susceptible and exposed harmonic mean. For several studies, the W ( S , E ) uptake function is considered to be the product of susceptible and exposed, known as the ref. [37] bi-linear incidence rate. The dynamics of the COVID-19 model were addressed by the authors of ref. [37] when considering the uptake feature to be the product of S and E . Apart from these, it is very interesting to consider W ( S , E ) ≔ 2 S E S + E as the uptake function. Since it is understood that if there are two quantities S and E in the sense that S , E ≥ 0 , then the following relation is true:

2 S E S + E ≤ S E ≤ S + E 2 .

In fact, the average of two values is a measure of the centrality of a data set. In addition, the geometric mean is mainly used to reach average data change ratios or rates. As far as harmonic means are concerned, as opposed to arithmetic or geometric means, it is less susceptible to a few broad values. It is used for highly skewed variables often. The readers are referred to refs [38,39, 40,41] for the biological interpretation of the harmonic mean.

2.1 Boundedness of the system

Here we examine the boundedness property of the system (2.1).

Theorem 1

The solutions of the system (2.1) are uniformly bounded.

Proof

We assume that X = S + E + Q + I + R . Therefore,

d X d t = d S d t + d E d t + d Q d t + d I d t + d R d t , d X d t + δ I = A − d X , i.e. , d X d t + d X ≤ A .

Integrating the above inequality and by applying the theorem of differential equation due to Birkhoff and Rota [42], we get

X ≤ A d + p M [ 1 − e − d t ] + X 0 e − d t .

Now for t → ∞

0 < X ≤ A d + p M .

Hence, all the solutions of (2.1) that are initiating in { R + 5 } are confined in the region

X ∈ R + 5 : 0 ≤ X ( S , E , Q , I , R ) < A d + p M + ε

for any ε > 0 and for t → ∞ . Hence the theorem is proved.

Remark 1

The conservative law as comes next satisfies (2.1):

(2.3) d N d t = A − d N ,

where N ( t ) represents the total population size with N ( t ) = S ( t ) + E ( t ) + Q ( t ) + I ( t ) + R ( t ) at any time with t . Equation (2.3) possesses the exact solution as comes next

N ( t ) = A d + N 0 − A d e − d t ,

with the initial condition

S ( 0 ) ≥ 0 , E ( 0 ) ≥ 0 , Q ( 0 ) ≥ 0 , I ( 0 ) ≥ 0 , R ( 0 ) ≥ 0 .

In addition to the above, we have

N ( 0 ) = S ( 0 ) + E ( 0 ) + Q ( 0 ) + I ( 0 ) + R ( 0 ) .

Thus,

S ( t ) ≥ 0 , E ( t ) ≥ 0 , Q ( t ) ≥ 0 , I ( t ) ≥ 0 , R ( t ) ≥ 0 .

It is to be noted that the solution is certainly positive.

Remark 2

V : [ 0 , ∞ ) → [ 0 , ∞ ) is the uptake function and fulfills some of the properties [43]:

V ( 0 ) = 0 , V ( S ) > 0 , for S > 0 .
lim P → ∞ S ( t ) S ( t ) + E ( t ) = L 1 , 0 < L 1 < ∞ ; where L 1 = 1 .
V ( S ) is continuously differentiable.
d V d t = E ( S + E ) 2 > 0 , which means that V is rising monotonically.

Definition 1

[44] In the context of ABC, an integral operator is provided as

(2.4) D 0 , t Ξ ABC [ f ( t ) ] = 1 − Ξ B ( Ξ ) f ( t ) + Ξ B ( Ξ ) Γ ( Ξ ) ∫ 0 t f ( s ) ( t − s ) Ξ − 1 d s .

2.2 Model in ABC sense

In this subsection, we reformulate the COVID-19 epidemic model (2.1) by incorporating the harmonic mean type incidence rate and definition (1). In order to observe the memory effects, the system (2.1), is fractionalized as:

(2.5) D 0 , t Ξ ABC [ S ( t ) ] = A Ξ − 2 β Ξ ( 1 − ρ 1 ) ( 1 − ρ 2 ) S E S + E + b 1 Ξ Q − d Ξ S − p Ξ S M Ξ , D 0 , t Ξ ABC [ E ( t ) ] = 2 β Ξ ( 1 − ρ 1 ) ( 1 − ρ 2 ) S E S + E − b 2 Ξ E − α Ξ E − σ Ξ E − d Ξ E , D 0 , t Ξ ABC [ Q ( t ) ] = b 2 Ξ E − b Ξ Q − c Ξ Q − d Ξ Q , D 0 , t Ξ ABC [ I ( t ) ] = α Ξ E + c Ξ Q − ( η Ξ + d Ξ + δ Ξ ) I , D 0 , t Ξ ABC [ R ( t ) ] = η Ξ I + σ Ξ E − d Ξ R + p Ξ S M Ξ .

it is subject to the following initial conditions:

S ( t ) , E ( t ) , Q ( t ) , I ( t ) , R ( t ) ≥ 0 .

2.3 Existence and positivity of the solution

We explore here the existence and positivity of the model (2.5).

Theorem 2

There exists a unique positive solution for the model (2.5) and remains in R + 5 .

Proof

To show that the solution of the system (2.5) is positive, we provide the result as follows:

(2.6) D 0 , t Ξ ABC [ S ( t ) ] ∣ S = 0 = A Ξ + b 1 Ξ Q ⩾ 0 , D 0 , t Ξ ABC [ E ( t ) ] ∣ E = 0 = 0 , D 0 , t Ξ ABC [ Q ( t ) ] ∣ Q = 0 = b 2 Ξ E ⩾ 0 , D 0 , t Ξ ABC [ I ( t ) ] ∣ I = 0 = α Ξ E + c Ξ Q ⩾ 0 , D 0 , t Ξ ABC [ R ( t ) ] ∣ R = 0 = η Ξ I + σ Ξ E + p Ξ S M Ξ ⩾ 0 ,

which shows that the model solution will stay in R + 5 for all time t ⩾ 0 . Furthermore, adding the equations in (2.5), we have

D 0 , t Ξ ABC N ( t ) = A − d N .

Furthermore, it follows

lim t → ∞ sup N ( t ) ⩽ A d

and thus, the biological feasible region for the model (2.5) can be shown by Θ ⊂ R + 5 , where

Θ = ( S ( t ) , E ( t ) , Q ( t ) , I ( t ) , R ( t ) ) ∈ R + 5 : N ( t ) ⩽ A d .

The model described above (2.5) for the COVID-19 in the fractional operator ABC is used further to obtain the associated results in the following section.□

3 Basic reproduction number R 0 and stability

E 0 represents free equilibrium point of the disease for (2.5) and is provided as

(3.1) E 0 = ( S 0 , E 0 , Q 0 , I 0 , R 0 ) = A Ξ d Ξ + p Ξ M Ξ , 0 , 0 , 0 , A Ξ p Ξ M Ξ d Ξ ( d Ξ + p Ξ M Ξ ) .

The infection’s spread and control are fundamentally linked to the R 0 threshold number. This total demonstrates the distribution and control of the infection. The infection disappears from the population on the off chance that this edge number is R 0 < 1 , and the infection-free equilibrium point exists that is asymptotically stable locally and globally. It helps with preventative measures to manage the spread of an epidemic. However, if, for example, R 0 > 1 the other way around, then under specific conditions, the endemic equilibria are locally and globally stable. The virus assumes the form of the plague and lives in society forever. According to this discussion, we are sure that R 0 measure the infection in the model which implies that R 0 only depend on the infected classes, i.e., ( E , Q , I ) . To find our model’s basic reproduction number, consider ( E , Q , I ) , the infected compartment in (2.5):

(3.2) D 0 , t Ξ ABC [ E ( t ) ] = 2 β Ξ ( 1 − ρ 1 ) ( 1 − ρ 2 ) S E S + E − b 2 Ξ E − α Ξ E − σ Ξ E − d Ξ E , D 0 , t Ξ ABC [ Q ( t ) ] = b 2 Ξ E − b Ξ Q − c Ξ Q − d Ξ Q , D 0 , t Ξ ABC [ I ( t ) ] = α Ξ E + c Ξ Q − ( η Ξ + d Ξ + δ Ξ ) I .

The Jacobian matrix J for the E 0 is provided by the next-generation matrix approach.

J = β Ξ ( 1 − ρ 1 ) ( 1 − ρ 2 ) − ( b 2 Ξ + α Ξ + σ Ξ + d Ξ ) 0 0 b 2 Ξ − ( b Ξ + c Ξ + d Ξ ) 0 α Ξ c Ξ − ( η Ξ + d Ξ + δ Ξ ) .

By separating J in the form of V and F , we attain

F = β Ξ ( 1 − ρ 1 ) ( 1 − ρ 2 ) 0 0 0 0 0 0 0 0

and

V = ( b 2 Ξ + α Ξ + σ Ξ + d Ξ ) 0 0 − b 2 Ξ ( b Ξ + c Ξ + d Ξ ) 0 − α Ξ − c Ξ + ( η Ξ + d Ξ + δ Ξ ) .

R 0 is provided by

R 0 = β Ξ ( 1 − ρ 1 ) ( 1 − ρ 2 ) ( b 2 Ξ + α Ξ + σ Ξ + d Ξ ) .

Theorem 3

The CFE point E 0 for (2.5) is asymptotically (locally) stable if R 0 < 1 and Ξ = 1 .

Proof

The Jacobian matrix for (2.5) at E 0 is

(3.3) J ( E 0 ) = − ( d Ξ + p Ξ M Ξ ) β Ξ ( 1 − ρ 1 ) ( 1 − ρ 2 ) b 1 Ξ 0 0 0 a 22 0 0 0 0 b 2 Ξ − ( b 1 Ξ + c Ξ + d Ξ ) 0 0 0 α Ξ c Ξ − ( η Ξ + d Ξ + δ Ξ ) 0 p Ξ M Ξ σ Ξ 0 η Ξ − d Ξ ,

where a 22 = β Ξ ( 1 − ρ 1 ) ( 1 − ρ 2 ) − ( b 2 Ξ + α Ξ + σ Ξ + d Ξ ) .

The characteristic equation of J ( E 0 ) looks:

[ ω 1 + d Ξ ] [ ω 2 + ( η Ξ + d Ξ + δ Ξ ) ] [ ω 3 + ( d Ξ + p Ξ M Ξ ) ] ( A ω 2 + B ω + C ) = 0 ,

where

(3.4) A = 1 , B = ( b 1 Ξ + c Ξ + d Ξ ) − β Ξ ( 1 − ρ 1 ) ( 1 − ρ 2 ) + ( b 2 Ξ + α Ξ + σ Ξ + d Ξ ) , = ( b 1 Ξ + c Ξ + d Ξ ) + ( b 2 Ξ + α Ξ + σ Ξ + d Ξ ) ( 1 − R 0 ) , C = − ( b 1 Ξ + c Ξ + d Ξ ) [ β Ξ ( 1 − ρ 1 ) ( 1 − ρ 2 ) − ( b 2 Ξ + α Ξ + σ Ξ + d Ξ ) ] , = ( b 1 Ξ + c Ξ + d Ξ ) ( b 2 Ξ + α Ξ + σ Ξ + d Ξ ) [ 1 − R 0 ] .

From the characteristic equation, it can be shown that the first three eigenvalues are negative, while we find a condition in which the real parts of the remaining two eigenvalues are negative.

We can see that B > 0 , A > 0 , and the C value can be negative or positive. By the instances described in ref. [45] for the value of C , it follows that R 0 < 1 , which completes the proof.□

4 Existence and uniqueness of the solutions

In this portion of the study, we will provide the existence of the unique solution for the proposed fractional-order model under the Atangana–Baleanu fractional operator in the Caputo sense by using the fixed point theory. The existence and uniqueness of the solution with regard to ABC will be revealed herein for (2.5). Suppose that a continuous real-value function denoted by the sup norm property containing a Banach space B ( J ) on J = [ 0 , b ] and P = B ( J ) × B ( J ) × B ( J ) × B ( J ) × B ( J ) with norm ∥ ( S , E , Q , I , R ) ∥ = ‖ S ‖ + ‖ E ‖ + ∥ Q ∥ + ∥ I ∥ + ‖ R ‖ , where ‖ S ‖ = sup t ∈ J ∣ S ( t ) ∣ , ‖ E ‖ = sup t ∈ j ∣ E ( t ) ∣ , ∥ Q ∥ = sup t ∈ j ∣ Q ( t ) ∣ , ∥ I ∥ = sup t ∈ j ∣ I ( t ) ∣ , ∥ R ∥ = sup t ∈ j ∣ R ( t ) ∣ . Utilizing the ABC fractional integral operator in both the sides of equation (2.5), we obtain

(4.1) S ( t ) − S ( 0 ) = D 0 , t Ξ ABC [ S ( t ) ] A Ξ − 2 β Ξ ( 1 − ρ 1 ) ( 1 − ρ 2 ) S E S + E + b 1 Ξ Q − d Ξ S − p Ξ S M Ξ } , E ( t ) − E ( 0 ) = D 0 , t Ξ ABC [ E ( t ) ] 2 β Ξ ( 1 − ρ 1 ) ( 1 − ρ 2 ) S E S + E − b 2 Ξ E − α Ξ E − σ Ξ E − d Ξ E , Q ( t ) − Q ( 0 ) = D 0 , t Ξ ABC [ Q ( t ) ] { b 2 Ξ E − b Ξ Q − c Ξ Q − d Ξ Q } , I ( t ) − I ( 0 ) = D 0 , t Ξ ABC [ I ( t ) ] { α Ξ E + c Ξ Q − ( η Ξ + d Ξ + δ Ξ ) I } , R ( t ) − R ( 0 ) = D 0 , t Ξ ABC [ R ( t ) ] { η Ξ I + σ Ξ E − d Ξ R + p Ξ S M Ξ } .

Now, on each of the above equations, utilizing definition (1), one reaches

(4.2) S ( t ) − S ( 0 ) = 1 − Ξ B ( Ξ ) K 1 ( Ξ , t , S ( t ) ) + Ξ B ( Ξ ) Γ ( Ξ ) × ∫ 0 t ( t − θ ) Ξ − 1 K 1 ( Ξ , θ , S ( θ ) ) d θ , E ( t ) − E ( 0 ) = 1 − Ξ B ( Ξ ) K 2 ( Ξ , t , E ( t ) ) + Ξ B ( Ξ ) Γ ( Ξ ) × ∫ 0 t ( t − θ ) Ξ − 1 K 2 ( Ξ , θ , E ( θ ) ) d θ , Q ( t ) − Q ( 0 ) = 1 − Ξ B ( Ξ ) K 3 ( Ξ , t , Q ( t ) ) + Ξ B ( Ξ ) Γ ( Ξ ) × ∫ 0 t ( t − θ ) Ξ − 1 K 3 ( Ξ , θ , Q ( θ ) ) d θ , I ( t ) − I ( 0 ) = 1 − Ξ B ( Ξ ) K 4 ( Ξ , t , I ( t ) ) + Ξ B ( Ξ ) Γ ( Ξ ) × ∫ 0 t ( t − θ ) Ξ − 1 K 4 ( Ξ , θ , I ( θ ) ) d θ , R ( t ) − R ( 0 ) = 1 − Ξ B ( Ξ ) K 5 ( Ξ , t , R ( t ) ) + Ξ B ( Ξ ) Γ ( Ξ ) × ∫ 0 t ( t − θ ) Ξ − 1 K 5 ( Ξ , θ , R ( θ ) ) d θ ,

where

(4.3) K 1 ( Ξ , t , S ( t ) ) = A Ξ − 2 β Ξ ( 1 − ρ 1 ) ( 1 − ρ 2 ) S E S + E + b 1 Ξ Q − d Ξ S − p Ξ S M Ξ , K 2 ( Ξ , t , E ( t ) ) = 2 β Ξ ( 1 − ρ 1 ) ( 1 − ρ 2 ) S E S + E − b 2 Ξ E − α Ξ E − σ Ξ E − d Ξ E , K 3 ( Ξ , t , Q ( t ) ) = b 2 Ξ E − b Ξ Q − c Ξ Q − d Ξ Q , K 4 ( Ξ , t , I ( t ) ) = α Ξ E + c Ξ Q − ( η Ξ + d Ξ + δ Ξ ) I , K 5 ( Ξ , t , R ( t ) ) = η Ξ I + σ Ξ E − d Ξ R + p Ξ S M Ξ .

The symbols K 1 , K 2 , K 3 , K 4 , and K 5 satisfy the Lipschitz condition only if S ( t ) , E ( t ) , Q ( t ) , I ( t ) , and R ( t ) possess an upper bound. Assuming that S ( t ) and S ∗ ( t ) are a few function, we reveal

(4.4) ∥ K 1 ( Ξ , t , S ( t ) ) − K 1 ( Ξ , t , S ∗ ( t ) ) ∥ = − 2 β Ξ ( 1 − ρ 1 ) ( 1 − ρ 2 ) E S + E + d Ξ + p Ξ M Ξ ( S ( t ) − S ∗ ( t ) ) ∥ .

Considering

η 1 ≔ − 2 β Ξ ( 1 − ρ 1 ) ( 1 − ρ 2 ) E S + E + d Ξ + p Ξ M Ξ

one reaches

(4.5) ∥ K 1 ( Ξ , t , S ( t ) ) − K 1 ( Ξ , t , S ∗ ( t ) ) ∥ ≤ η 1 ∥ S ( t ) − S ∗ ( t ) ∥ .

To proceed in the same manner, one reaches

(4.6) ∥ K 2 ( Ξ , t , E ( t ) ) − K 2 ( Ξ , t , E ∗ ( t ) ) ∥ ≤ η 2 ∥ E ( t ) − E ∗ ( t ) ∥ , ∥ K 3 ( Ξ , t , Q ( t ) ) − K 3 ( Ξ , t , Q ∗ ( t ) ) ∥ ≤ η 3 ∥ Q ( t ) − Q ∗ ( t ) ∥ , ∥ K 4 ( Ξ , t , I ( t ) ) − K 4 ( Ξ , t , I ∗ ( t ) ) ∥ ≤ η 4 ∥ I ( t ) − I ∗ ( t ) ∥ , ∥ K 5 ( Ξ , t , R ( t ) ) − K 5 ( Ξ , t , R ∗ ( t ) ) ∥ ≤ η 5 ∥ R ( t ) − R ∗ ( t ) ∥ ,

where

η 2 = − 2 β Ξ ( 1 − ρ 1 ) ( 1 − ρ 2 ) S S + E + ( b 2 Ξ + α Ξ + σ Ξ + d Ξ ) , η 3 = ∥ ( − ( b Ξ + c Ξ + d Ξ ) ) ∥ , η 4 = ∥ ( − ( η Ξ + d Ξ + δ Ξ ) ) ∥ , η 5 = ∥ ( − d Ξ ) ∥ .

This means that for all the five functions, the Lipschitz condition holds. In a recursive way, the expressions in (4.2) provide

(4.7) S n ( t ) − S ( 0 ) = 1 − Ξ B ( Ξ ) K 1 ( Ξ , t , S n − 1 ( t ) ) + Ξ B ( Ξ ) Γ ( Ξ ) × ∫ 0 t ( t − θ ) Ξ − 1 K 1 ( Ξ , θ , S n − 1 ( θ ) ) d θ , E n ( t ) − E ( 0 ) = 1 − Ξ B ( Ξ ) K 2 ( Ξ , t , E n − 1 ( t ) ) + Ξ B ( Ξ ) Γ ( Ξ ) × ∫ 0 t ( t − θ ) Ξ − 1 K 2 ( Ξ , θ , E n − 1 ( θ ) ) d θ , Q n ( t ) − Q ( 0 ) = 1 − Ξ B ( Ξ ) K 3 ( Ξ , t , Q n − 1 ( t ) ) + Ξ B ( Ξ ) Γ ( Ξ ) × ∫ 0 t ( t − θ ) Ξ − 1 K 3 ( Ξ , θ , Q n − 1 ( θ ) ) d θ , I n ( t ) − I ( 0 ) = 1 − Ξ B ( Ξ ) K 4 ( Ξ , t , I n − 1 ( t ) ) + Ξ B ( Ξ ) Γ ( Ξ ) × ∫ 0 t ( t − θ ) Ξ − 1 K 4 ( Ξ , θ , I n − 1 ( θ ) ) d θ , R n ( t ) − R ( 0 ) = 1 − Ξ B ( Ξ ) K 5 ( Ξ , t , R n − 1 ( t ) ) + Ξ B ( Ξ ) Γ ( Ξ ) × ∫ 0 t ( t − θ ) Ξ − 1 K 5 ( Ξ , θ , R n − 1 ( θ ) ) d θ ,

together with S 0 ( t ) = S ( 0 ) , E 0 ( t ) = E ( 0 ) , Q 0 ( t ) = Q ( 0 ) , I 0 ( t ) = I ( 0 ) , and R 0 ( t ) = R ( 0 ) . Now,

(4.8) χ S , n ( t ) = S n ( t ) − S n − 1 ( t ) = 1 − Ξ B ( Ξ ) ( K 1 ( Ξ , t , S n − 1 ( t ) ) − K 1 ( Ξ , t , S n − 2 ( t ) ) ) , + Ξ B ( Ξ ) Γ ( Ξ ) ∫ 0 t ( t − θ ) Ξ − 1 ( K 1 ( Ξ , θ , S n − 1 ( θ ) ) − K 1 ( Ξ , θ , S n − 2 ( θ ) ) ) d θ , χ E , n ( t ) = E n ( t ) − E n − 1 ( t ) = 1 − Ξ B ( Ξ ) ( K 2 ( Ξ , t , E n − 1 ( t ) ) − K 2 ( Ξ , t , E n − 2 ( t ) ) ) , + Ξ B ( Ξ ) Γ ( Ξ ) ∫ 0 t ( t − θ ) Ξ − 1 ( K 2 ( Ξ , θ , E n − 1 ( θ ) ) − K 2 ( Ξ , θ , E n − 2 ( θ ) ) ) d θ , χ Q , n ( t ) = I 1 n ( t ) − Q n − 1 ( t ) = 1 − Ξ B ( Ξ ) ( K 3 ( Ξ , t , Q n − 1 ( t ) ) − K 3 ( Ξ , t , Q n − 2 ( t ) ) ) , + Ξ B ( Ξ ) Γ ( Ξ ) ∫ 0 t ( t − θ ) Ξ − 1 ( K 3 ( Ξ , θ , Q n − 1 ( θ ) ) − K 3 ( Ξ , θ , Q n − 2 ( θ ) ) ) d θ , χ I , n ( t ) = I 2 n ( t ) − I n − 1 ( t ) = 1 − Ξ B ( Ξ ) ( K 4 ( Ξ , t , I n − 1 ( t ) ) − K 4 ( Ξ , t , I n − 2 ( t ) ) ) , + Ξ B ( Ξ ) Γ ( Ξ ) ∫ 0 t ( t − θ ) Ξ − 1 ( K 4 ( Ξ , θ , I n − 1 ( θ ) ) − K 4 ( Ξ , θ , I n − 2 ( θ ) ) ) d θ , Ξ R , n ( t ) = F n ( t ) − R n − 1 ( t ) = 1 − Ξ B ( Ξ ) ( K 5 ( Ξ , t , R n − 1 ( t ) ) − K 5 ( Ξ , t , R n − 2 ( t ) ) ) , + Ξ B ( Ξ ) Γ ( Ξ ) ∫ 0 t ( t − θ ) Ξ − 1 ( K 5 ( Ξ , θ , R n − 1 ( θ ) ) − K 5 ( Ξ , θ , R n − 2 ( θ ) ) ) d θ .

It is vital to observe that

S n ( t ) = ∑ i = 0 n Ξ S , i ( t ) , E n ( t ) = ∑ i = 0 n Ξ E , i ( t ) , Q n ( t ) = ∑ i = 0 n Ξ Q , i ( t ) , I n ( t ) = ∑ i = 0 n Ξ I , i ( t ) , R n ( t ) = ∑ i = 0 n Ξ R , i ( t ) .

Additionally, by using equations (4.5)–(4.6) and considering that

Ξ S , n − 1 ( t ) = S n − 1 ( t ) − S n − 2 ( t ) , Ξ E , n − 1 ( t ) = E n − 1 ( t ) − E n − 2 ( t ) , Ξ Q , n − 1 ( t ) = Q n − 1 ( t ) − Q n − 2 ( t ) , Ξ I , n − 1 ( t ) = I n − 1 ( t ) − I n − 2 ( t ) , Ξ R , n − 1 ( t ) = R n − 1 ( t ) − R n − 2 ( t ) ,

we reach

(4.9) ∥ χ S , n ( t ) ∥ ≤ 1 − Ξ B ( Ξ ) η 1 ∥ Ξ S , n − 1 ( t ) ∥ Ξ B ( Ξ ) Γ ( Ξ ) η 1 × ∫ 0 t ( t − θ ) Ξ − 1 ∥ Ξ S , n − 1 ( θ ) ∥ d θ , ∥ χ E , n ( t ) ∥ ≤ 1 − Ξ B ( Ξ ) η 2 ∥ Ξ E , n − 1 ( t ) ∥ Ξ B ( Ξ ) Γ ( Ξ ) η 2 × ∫ 0 t ( t − θ ) Ξ − 1 ∥ Ξ E , n − 1 ( θ ) ∥ d θ , ∥ χ Q , n ( t ) ∥ ≤ 1 − Ξ B ( Ξ ) η 3 ∥ Ξ Q , n − 1 ( t ) ∥ Ξ B ( Ξ ) Γ ( Ξ ) η 3 × ∫ 0 t ( t − θ ) Ξ − 1 ∥ Ξ Q , n − 1 ( θ ) ∥ d θ , ∥ χ I , n ( t ) ∥ ≤ 1 − Ξ B ( Ξ ) η 4 ∥ Ξ I , n − 1 ( t ) ∥ Ξ B ( Ξ ) Γ ( Ξ ) η 4 × ∫ 0 t ( t − θ ) Ξ − 1 ∥ Ξ I , n − 1 ( θ ) ∥ d θ , ∥ χ R , n ( t ) ∥ ≤ 1 − Ξ B ( Ξ ) η 5 ∥ Ξ R , n − 1 ( t ) ∥ Ξ B ( Ξ ) Γ ( Ξ ) η 5 × ∫ 0 t ( t − θ ) Ξ − 1 ∥ Ξ R , n − 1 ( θ ) ∥ d θ .

Theorem 4

Assuming that the following condition holds

(4.10) 1 − Ξ B ( Ξ ) η i + Ξ B ( Ξ ) Γ ( Ξ ) b Ξ η i < 1 , i = 1 , 2 , … , 5 .

Then, (2.5) has a unique solution for t ∈ [ 0 , b ] .

Proof

It has been seen that S ( t ) , E ( t ) , Q ( t ) , I ( t ) , and R ( t ) are bounded functions. In addition, as can be seen from equations (4.5) and (4.6), the symbols K 1 , K 2 , K 3 , K 4 , and K 5 hold for Lipchitz condition. Thus, considering equation (4.9) along with a recursive hypothesis, we reach

(4.11) ∥ χ S , n ( t ) ∥ ≤ ∥ S 0 ( t ) ∥ 1 − Ξ B ( Ξ ) η 1 + Ξ b Ξ B ( Ξ ) Γ ( Ξ ) η 1 n , ∥ χ E , n ( t ) ∥ ≤ ∥ E 0 ( t ) ∥ 1 − Ξ B ( Ξ ) η 3 + Ξ b Ξ B ( Ξ ) Γ ( Ξ ) η 2 n , ∥ χ Q , n ( t ) ∥ ≤ ∥ Q 0 ( t ) ∥ 1 − Ξ B ( Ξ ) η 3 + Ξ b Ξ B ( Ξ ) Γ ( Ξ ) η 3 n , ∥ χ I , n ( t ) ∥ ≤ ∥ I 0 ( t ) ∥ 1 − Ξ B ( Ξ ) η 4 + Ξ b Ξ B ( Ξ ) Γ ( Ξ ) η 4 n , ∥ χ R , n ( t ) ∥ ≤ ∥ R 0 ( t ) ∥ 1 − Ξ B ( Ξ ) η 5 + Ξ b Ξ B ( Ξ ) Γ ( Ξ ) η 5 n .

It can be notice that

(4.12) ∥ χ S , n ( t ) ∥ → 0 , ∥ χ E , n ( t ) ∥ → 0 , ∥ χ Q , n ( t ) ∥ → 0 , ∥ χ I , n ( t ) ∥ → 0 , ∥ χ R , n ( t ) ∥ → 0 ,

as n → ∞ . In addition, from equation (4.11) and applying the triangle inequality, for any k , we reveal

(4.13) ∥ S n + k ( t ) − S n ( t ) ∥ ≤ ∑ j = n + 1 n + k Z 1 j = Z 1 n + 1 − Z 1 n + k + 1 1 − Z 1 , ∥ E n + k ( t ) − E n ( t ) ∥ ≤ ∑ j = n + 1 n + k Z 2 j = Z 2 n + 1 − Z 2 n + k + 1 1 − Z 2 , ∥ Q n + k ( t ) − Q n ( t ) ∥ ≤ ∑ j = n + 1 n + k Z 3 j = Z 3 n + 1 − Z 3 n + k + 1 1 − Z 3 , ∥ I n + k ( t ) − I n ( t ) ∥ ≤ ∑ j = n + 1 n + k Z 4 j = Z 4 n + 1 − Z 4 n + k + 1 1 − Z 4 , ∥ R n + k ( t ) − R n ( t ) ∥ ≤ ∑ i = n + 1 n + k Z 5 j = Z 5 n + 1 − Z 5 n + k + 1 1 − Z 5 ,

with Z i = 1 − Ξ B ( Ξ ) η i + Ξ B ( Ξ ) Γ ( Ξ ) b Ξ η i < 1 by hypothesis. Now, in Banach B ( J ) space, S n , E n , Q n , I n , and R n can be seen as a Cauchy sequence. This showed that ref. [46] is a uniform convergent. By bearing the limits principle in equation (4.8) as n → ∞ confirms the uniqueness of the limit for the sequences of (2.5). This ensures that a unique solution for equation (2.5) exists under the condition (4.10).□

5 Hyers–Ulam stability

Definition 2

[47] The considered model (2.5) is Hyers–Ulam stable if there exists a real number C ≥ 0 , such that for every ε > 0 and for any solution Y ∈ C 1 ( G , R ) of the following inequality:

∣ D Ξ ABC Y ( t ) − Ψ ( t , Y ( t ) ) ∣ ≤ ε , t ∈ G

there is a unique solution W ∈ C 1 ( G , R ) of the considered model (2.5), such that

∣ Y ( t ) − W ( t ) ∣ ≤ C ε , t ∈ G .

Note: For the qualitative analysis, we define Banach space S = A × A × A × A × A , where A = C ( G ) under

the norm ‖ Y ‖ = ‖ ( S , E , Q , I , R ) ‖ = max t ∈ G [ ∣ S ( t ) ∣ + ∣ E ( t ) ∣ + ∣ Q ( t ) ∣ + ∣ I ( t ) ∣ + ∣ R ( t ) ∣ ] .

Definition 3

The A B C fractional integral system provided by equation (4.2) is called Hyers–Ulam stable if exist constants Δ i ≥ 0 , i ∈ N 5 satisfying: For every γ i > 0 , i ∈ N 5 , for

(5.1) S ( t ) − 1 − Ξ B ( Ξ ) K 1 ( Ξ , t , S ( t ) ) + Ξ B ( Ξ ) Γ ( Ξ ) × ∫ 0 t ( t − θ ) Ξ − 1 K 1 ( Ξ , θ , S ( θ ) ) d θ ≤ γ 1 , E ( t ) − 1 − Ξ B ( Ξ ) K 2 ( Ξ , t , E ( t ) ) + Ξ B ( Ξ ) Γ ( Ξ ) × ∫ 0 t ( t − θ ) Ξ − 1 K 2 ( Ξ , θ , E ( θ ) ) d θ ≤ γ 2 , Q ( t ) − 1 − Ξ B ( Ξ ) K 3 ( Ξ , t , Q ( t ) ) + Ξ B ( Ξ ) Γ ( Ξ ) × ∫ 0 t ( t − θ ) Ξ − 1 K 3 ( Ξ , θ , Q ( θ ) ) d θ ≤ γ 3 , I ( t ) − 1 − Ξ B ( Ξ ) K 4 ( Ξ , t , I ( t ) ) + Ξ B ( Ξ ) Γ ( Ξ ) × ∫ 0 t ( t − θ ) Ξ − 1 K 4 ( Ξ , θ , I ( θ ) ) d θ ≤ γ 4 , R ( t ) − 1 − Ξ B ( Ξ ) K 5 ( Ξ , t , R ( t ) ) + Ξ B ( Ξ ) Γ ( Ξ ) × ∫ 0 t ( t − θ ) Ξ − 1 K 5 ( Ξ , θ , R ( θ ) ) d θ ≤ γ 5 ,

there exist ( S ˙ ( t ) , E ˙ ( t ) , Q ˙ ( t ) , I ˙ ( t ) , R ˙ ( t ) ) , which are satisfying

(5.2) S ˙ ( t ) = 1 − Ξ B ( Ξ ) K 1 ( Ξ , t , S ( t ) ) + Ξ B ( Ξ ) Γ ( Ξ ) × ∫ 0 t ( t − θ ) Ξ − 1 K 1 ( Ξ , θ , S ˙ ( θ ) ) d θ , E ˙ ( t ) = 1 − Ξ B ( Ξ ) K 2 ( Ξ , t , E ( t ) ) + Ξ B ( Ξ ) Γ ( Ξ ) × ∫ 0 t ( t − θ ) Ξ − 1 K 2 ( Ξ , θ , E ˙ ( θ ) ) d θ , Q ˙ ( t ) = 1 − Ξ B ( Ξ ) K 3 ( Ξ , t , Q ( t ) ) + Ξ B ( Ξ ) Γ ( Ξ ) × ∫ 0 t ( t − θ ) Ξ − 1 K 3 ( Ξ , θ , Q ˙ ( θ ) ) d θ , I ˙ ( t ) = 1 − Ξ B ( Ξ ) K 4 ( Ξ , t , I ( t ) ) + Ξ B ( Ξ ) Γ ( Ξ ) × ∫ 0 t ( t − θ ) Ξ − 1 K 4 ( Ξ , θ , I ˙ ( θ ) ) d θ , R ˙ ( t ) = 1 − Ξ B ( Ξ ) K 5 ( Ξ , t , R ( t ) ) + Ξ B ( Ξ ) Γ ( Ξ ) × ∫ 0 t ( t − θ ) Ξ − 1 K 5 ( Ξ , θ , R ˙ ( θ ) ) d θ .

Such that

(5.3) ∣ S ( t ) − S ˙ ( t ) ∣ ≤ ζ 1 γ 1 , ∣ E ( t ) − E ˙ ( t ) ∣ ≤ ζ 2 γ 2 , ∣ Q ( t ) − Q ˙ ( t ) ∣ ≤ ζ 3 γ 3 , ∣ I ( t ) − I ˙ ( t ) ∣ ≤ ζ 4 γ 4 , ∣ R ( t ) − R ˙ ( t ) ∣ ≤ ζ 5 γ 5 .

Theorem 5

With Assumption J , the suggested model of fractional order (2.5) is Hyers–Ulam stable.

Proof

The proposed A B C fractional model (2.5) with the help of theorem (4) has a unique solution:

Let ( S ( t ) , E ( t ) , Q ( t ) , I ( t ) , R ( t ) ) be approximate solution of model (2.5) satisfying the system equations (4.2). Thus,

(5.4) ‖ S ( t ) − S ˙ ( t ) ‖ ≤ 1 − Ξ B ( Ξ ) ∥ K 1 ( Ξ , t , S ( t ) ) − K 1 ( Ξ , t , S ˙ ( t ) ) ∥ + Ξ B ( Ξ ) Γ ( Ξ ) ∫ 0 t ( t − θ ) Ξ − 1 × ∥ K 1 ( Ξ , t , S ( t ) ) − K 1 ( Ξ , t , S ˙ ( t ) ) ∥ d θ ≤ 1 − Ξ B ( Ξ ) + Ξ B ( Ξ ) Γ ( Ξ ) ϕ 1 ‖ S − S ˙ ‖ ,

(5.5) ‖ E ( t ) − E ˙ ( t ) ‖ ≤ 1 − Ξ B ( Ξ ) ∥ K 2 ( Ξ , t , E ( t ) ) − K 2 ( Ξ , t , E ˙ ( t ) ) ∥ + Ξ B ( Ξ ) Γ ( Ξ ) ∫ 0 t ( t − θ ) Ξ − 1 × ∥ K 2 ( Ξ , t , E ( t ) ) − K 1 ( Ξ , t , E ˙ ( t ) ) ∥ d θ ≤ 1 − Ξ B ( Ξ ) + Ξ B ( Ξ ) Γ ( Ξ ) ϕ 2 ‖ E − E ˙ ‖ ,

(5.6) ‖ Q ( t ) − Q ˙ ( t ) ‖ ≤ 1 − Ξ B ( Ξ ) ∥ K 3 ( Ξ , t , Q ( t ) ) − K 3 ( Ξ , t , Q ˙ ( t ) ) ∥ + Ξ B ( Ξ ) Γ ( Ξ ) ∫ 0 t ( t − θ ) Ξ − 1 × ∥ K 3 ( Ξ , t , Q ( t ) ) − K 3 ( Ξ , t , Q ˙ ( t ) ) ∥ d θ ≤ 1 − Ξ B ( Ξ ) + Ξ B ( Ξ ) Γ ( Ξ ) ϕ 3 ‖ Q − Q ˙ ‖ ,

(5.7) ‖ I ( t ) − I ˙ ( t ) ‖ ≤ 1 − Ξ B ( Ξ ) ∥ K 4 ( Ξ , t , I ( t ) ) − K 4 ( Ξ , t , Q ˙ ( t ) ) ∥ + Ξ B ( Ξ ) Γ ( Ξ ) ∫ 0 t ( t − θ ) Ξ − 1 × ∥ K 4 ( Ξ , t , I ( t ) ) − K 1 ( Ξ , t , I ˙ ( t ) ) ∥ d θ ≤ 1 − Ξ B ( Ξ ) + Ξ B ( Ξ ) Γ ( Ξ ) ϕ 4 ‖ I − I ˙ ‖ ,

(5.8) ‖ R ( t ) − R ˙ ( t ) ‖ ≤ 1 − Ξ B ( Ξ ) ∥ K 5 ( Ξ , t , R ( t ) ) − K 5 ( Ξ , t , R ˙ ( t ) ) ∥ + Ξ B ( Ξ ) Γ ( Ξ ) ∫ 0 t ( t − θ ) Ξ − 1 × ∥ K 5 ( Ξ , t , R ( t ) ) − K 5 ( Ξ , t , R ˙ ( t ) ) ∥ d θ ≤ 1 − Ξ B ( Ξ ) + Ξ B ( Ξ ) Γ ( Ξ ) ϕ 5 ‖ R − R ˙ ‖ .

Taking, γ i = ϕ i , Δ i = 1 − Ξ B ( Ξ ) + Ξ B ( Ξ ) Γ ( Ξ ) , this implies

(5.9) ‖ S ( t ) − S ˙ ( t ) ‖ ≤ γ 1 Δ 1 .

Likewise, in the same way, we provide

(5.10) ‖ E ( t ) − E ˙ ( t ) ‖ ≤ γ 2 Δ 2 ‖ Q ( t ) − Q ˙ ( t ) ‖ ≤ γ 3 Δ 3 ‖ I ( t ) − I ˙ ( t ) ‖ ≤ γ 4 Δ 4 ‖ R ( t ) − R ˙ ( t ) ‖ ≤ γ 5 Δ 5 .

A B C fractional integral method with the aid of equations (5.9), (5.10), and (4.2) is Hyers–Ulam and consequently the ABC-fractional order model (2.5) is Hyers–Ulam stable.□

5.1 Sensitivity analysis

In this subsection, we perform the sensitivity analysis for the model under consideration. The elasticity indices as well as the parameter values are depicted in Table 1. The sensitivity of various parameters with R 0 is demonstrated in Figures 1, 2, 3, 4, 5. Reaching the parameters that are useful in minimizing the infectious disease transmission is done by sensitivity analysis. Although its estimation is tedious for complex biological models, forward sensitivity analysis is considered a vital component of disease modeling. The ecologist and epidemiologist gained a lot of attention from the sensitivity study of R 0 .

Table 1

Sensitivity indices of the reproduction number R 0 against mentioned parameters

Parameters	S. index	Value	Parameters	S. index	Value
β	S β	1.00000000	d	S d	−0.4999885937
σ	S σ	−0.7142914285	b 2	S b 2	−0.3711929652
A	S A	1.00000000	ρ 1	S ρ 1	−0.7213098196
ρ 2	S ρ 1	−0.4285717142	α	S α	0.31632629344
M	S M	0.1245906249	p	S p	−0.2446571428

$Figure 1 Plot of R 0 {R}_{0} against M M and A A .$

Figure 1

Plot of R 0 against M and A .

$Figure 2 Plot of R 0 {R}_{0} against β \beta and d d .$

Figure 2

Plot of R 0 against β and d .

$Figure 3 Plot of R 0 {R}_{0} against ρ 2 {\rho }_{2} and d d .$

Figure 3

Plot of R 0 against ρ 2 and d .

$Figure 4 Plot of R 0 {R}_{0} against M M and d d .$

Figure 4

Plot of R 0 against M and d .

$Figure 5 Plot of R 0 {R}_{0} against β \beta and A A .$

Figure 5

Plot of R 0 against β and A .

Definition 1

The normalized forward sensitivity index of the R 0 that depends differentially on a parameter ω is defined as

(5.11) S ω = ω R 0 ∂ R 0 ∂ ω .

6 Numerical simulations and discussion

Parameters	Description	Values/ranges
A	Total recruitment	50
β	Disease transmission rate	[0.5,2.3]
ρ 1	Portion of S contact with E	(0,1)
ρ 2	Portion of E contact with S	(0,1)
d	Natural death rate	0.2
b 1	The rate that Q becomes S	0.25
b 2	The rate that E becomes quarantine	0.8
α	The rate that E becomes I	0.3
η	The rate that I becomes R naturally	0.25
σ	The rate that E becomes R naturally	0.2
c	The rate that Q becomes I	0.12
δ	The mortality rate for I	0.25
M	Policy parameter	0.8
p	Implementation rate of policy	0.78

The ABC version of (2.5) is numerically simulated here, as suggested in refs [48,49,50], using first-order convergent numerical methods. These numerical methods are reliable, conditionally stable, and convergent for solving fractional-order linear and nonlinear systems of ordinary differential equations.

Consider a general Cauchy problem of fractional order having autonomous nature

(6.1) D 0 , t Ξ ⋆ ( y ( t ) ) = g ( y ( t ) ) , Ξ ∈ ( 0 , 1 ] , t ∈ [ 0 , T ] , y ( 0 ) = y 0 ,

where y = ( a , b , c , w ) ∈ R + 4 is a real-valued continuous vector function which satisfies the Lipchitz criterion given as

(6.2) ‖ g ( y 1 ( t ) ) − g ( y 2 ( t ) ) ‖ ≤ M ∣ ∣ y 1 ( t ) − y 2 ( t ) ∣ ∣ ,

where M is a positive real Lipchitz constant.

Using the fractional-order integral operators, one obtains

(6.3) y ( t ) = y 0 + J 0 , t Ξ g ( y ( t ) ) , t ∈ [ 0 , T ] ,

where J 0 , t Ω is the fractional-order integral operator in Riemann–Liouville. Consider equispaced integration intervals over [ 0 , T ] with the fixed step size h ( = 1 0 − 2 for simulation ) = T n , n ∈ N . Suppose that y p be the approximation of y ( t ) at t = t p for p = 0 , 1 , … n . The numerical technique for the governing model under Caputo fractional derivative operator takes the form

(6.4) S p + 1 c = a 0 + h Ξ Γ ( Ξ + 1 ) ∑ k = 0 p ( ( p − k + 1 ) Ξ − ( p − k ) Ξ ) ( A Ξ − β Ξ ( 1 − ρ 1 ) ( 1 − ρ 2 ) S E + b 1 Ξ Q − d Ξ S − p Ξ S M Ξ ) , E p + 1 c = b 0 + h Ξ Γ ( Ξ + 1 ) ∑ k = 0 p ( ( p − k + 1 ) Ξ − ( p − k ) Ξ ) ( β Ξ ( 1 − ρ 1 ) ( 1 − ρ 2 ) S E − b 2 Ξ E − α Ξ E − σ Ξ E − d Ξ E ) , Q p + 1 c = d 0 + h Ξ Γ ( Ξ + 1 ) ∑ k = 0 p ( ( p − k + 1 ) Ξ − ( p − k ) Ξ ) ( b 2 Ξ E − b Ξ Q − c Ξ Q − d Ξ Q ) , I p + 1 c = e 0 + h Ξ Γ ( Ξ + 1 ) ∑ k = 0 p ( ( p − k + 1 ) Ξ − ( p − k ) Ξ ) ( α Ξ E + c Ξ Q − ( η Ξ + d Ξ + δ Ξ ) I ) , R p + 1 c = f 0 + h Ξ Γ ( Ξ + 1 ) ∑ k = 0 p ( ( p − k + 1 ) Ξ − ( p − k ) Ξ ) ( η Ξ I + σ Ξ E − d Ξ R + p Ξ S M Ξ ) .

We are now discussing the numerical results of the governing model reported with respect to the approximate solutions. To this end, under the Caputo fractional operator, we utilize the efficient Euler method to do the job. As described in the above table, the initial conditions and parameter values are used. As explained in Section 2.1, the initial conditions as well as the values of the parameters are used in performing the simulation results. In Figure 6, the physical perspective of the individual state variables of the model under the Caputo fractional operator is displayed. It can be noted in Figure 7 that the fractional order Ξ = is varied by 1 , 0.95 , 0.85 . The robust nature of the Caputo operator can easily be seen as an integer variant of the model. As shown in Figure 8, I ( t ) is also increasing for a decreasing varying β (disease transmission rate) values. Similarly, the effect of the mortality rate δ on I ( t ) is demonstrated in Figure 9. An increasing pattern in I ( t ) is observed for an increasing value of δ as in Figure 9(a). Similarly, a decreasing pattern in I ( t ) is noted for a decreasing δ value. In this case, an increasing decreasing pattern is portrayed.

$Figure 6 Profiles for behavior of (a) Susceptible (b) Exposed (c) Quarantined (d) Infected and (e) Recovered population for the Atangana–Baleanu version of the fractional model using the values of the parameters.$

Figure 6

Profiles for behavior of (a) Susceptible (b) Exposed (c) Quarantined (d) Infected and (e) Recovered population for the Atangana–Baleanu version of the fractional model using the values of the parameters.

$Figure 7 The dynamics of (a) Susceptible (b) Exposed (c) Quarantined (d) Infected and (e) Recovered population for different Ξ \Xi values.$

Figure 7

The dynamics of (a) Susceptible (b) Exposed (c) Quarantined (d) Infected and (e) Recovered population for different Ξ values.

$Figure 8 Behavior of the infectious class I ( t ) I\left(t) for (a) decreasing values of β \beta (disease transmission rate) and (b) increasing values of β \beta (disease).$

Figure 8

Behavior of the infectious class I ( t ) for (a) decreasing values of β (disease transmission rate) and (b) increasing values of β (disease).

$Figure 9 Behavior of the infectious class I ( t ) I\left(t) for (a) increasing values of δ \delta (mortality rate) and (b) decreasing values of δ \delta (mortality rate).$

Figure 9

Behavior of the infectious class I ( t ) for (a) increasing values of δ (mortality rate) and (b) decreasing values of δ (mortality rate).

7 Conclusion

COVID-19 model is fractionalized by considering the ABC operator. The harmonic mean type incidence rate is taken into consideration. The reproduction number of the model along with local stability is calculated. Fractional Hyers–Ulam stability approach with Mittag–Leffler Law has been analyzed. The existence and uniqueness of solutions are investigated using Banach fixed point theorem. A fractional numerical scheme for the coronavirus model in ABC sense was constructed. Sensitivity analysis has been used to identify the most highly sensitive parameter and it is noted that β and A , i.e., disease transmission rate and recruitment rate have a high impact on basic reproduction number. While the control policy parameter M has a less impact as compared to other parameters on the basic reproduction number.

Funding information: The authors gratefully acknowledge the support from the Petchra Pra Jom Klao PhD research scholarship funded by King Mongkut’s University of Technology Thonburi, Thailand.
Author contributions: All authors have accepted responsibility for the entire content of this manuscript and approved its submission.
Conflict of interest: The authors state no conflict of interest.

References

[1] Waris A, Atta UK, Ali M, Asmat A, Baset AJ. COVID-19 outbreak: current scenario of Pakistan. New Microbes New Infect. 2020 May 1;35:100681. 10.1016/j.nmni.2020.100681Search in Google Scholar

[2] Wang J, Pang J, Liu X. Modelling diseases with relapse and nonlinear incidence of infection: a multi-group epidemic model. J Biol Dyn. 2014 Jan 1;8(1):99–116. 10.1080/17513758.2014.912682Search in Google Scholar

[3] Wang J, Zhang R, Kuniya T. The stability analysis of an SVEIR model with continuous age-structure in the exposed and infectious classes. J Biol Dyn. 2015 Jan 1;9(1):73–101. 10.1080/17513758.2015.1006696Search in Google Scholar

[4] Castillo-Chavez C, Blower S, Van den Driessche P, Kirschner D, Yakubu AA, editors. Mathematical approaches for emerging and reemerging infectious diseases: models, methods, and theory. New York: Springer Science & Business Media; 2002 May 2. 10.1007/978-1-4757-3667-0Search in Google Scholar

[5] Zhao S, Xu Z, Lu Y. A mathematical model of hepatitis B virus transmission and its application for vaccination strategy in China. Int J Epidemiol. 2000 Aug 1;29(4):744–52. 10.1093/ije/29.4.744Search in Google Scholar

[6] Khan K, Zarin R, Khan A, Yusuf A, Al-Shomrani M, Ullah A. Stability analysis of five-grade Leishmania epidemic model with harmonic mean-type incidence rate. Adv Differ Equ. 2021 Dec;2021(1):1–27. 10.1186/s13662-021-03249-4Search in Google Scholar

[7] Lahrouz A, Omari L, Kiouach D, Belmaâti A. Complete global stability for an SIRS epidemic model with generalized non-linear incidence and vaccination. Appl Math Comput. 2012 Feb 5;218(11):6519–25. 10.1016/j.amc.2011.12.024Search in Google Scholar

[8] Li MY, Muldowney JS. Global stability for the SEIR model in epidemiology. Math Biosci. 1995 Feb 1;125(2):155–64. 10.1016/0025-5564(95)92756-5Search in Google Scholar

[9] Zaman G, Kang YH, Jung IH. Stability analysis and optimal vaccination of an SIR epidemic model. BioSyst. 2008 Sep 1;93(3):240–9. 10.1016/j.biosystems.2008.05.004Search in Google Scholar PubMed

[10] Zou L, Zhang W, Ruan S. Modeling the transmission dynamics and control of hepatitis B virus in China. J Theor Biol. 2010 Jan 21;262(2):330–8. 10.1016/j.jtbi.2009.09.035Search in Google Scholar PubMed

[11] Mwasa A, Tchuenche JM. Mathematical analysis of a cholera model with public health interventions. Biosyst. 2011 Sep 1;105(3):190–200. 10.1016/j.biosystems.2011.04.001Search in Google Scholar PubMed

[12] Khan A, Zarin R, Hussain G, Usman AH, Humphries UW, Gomez-Aguilar JF. Modeling and sensitivity analysis of HBV epidemic model with convex incidence rate. Results Phys. 2021 Mar 1 22:103836. 10.1016/j.rinp.2021.103836Search in Google Scholar

[13] London WP, Yorke JA. Recurrent outbreaks of measles, chickenpox and mumps: I. Seasonal variation in contact rates. Amer J Epidemiol. 1973 Dec 1;98(6):453–68. 10.1093/oxfordjournals.aje.a121575Search in Google Scholar PubMed

[14] Liu WM, Hethcote HW, Levin SA. Dynamical behavior of epidemiological models with nonlinear incidence rates. J Math Biol. 1987 Sep;25(4):359–80. 10.1007/BF00277162Search in Google Scholar PubMed

[15] Kamien MI, Schwartz NL. Dynamic optimization: the calculus of variations and optimal control in economics and management. Amsterdam: Courier corporation; 2012 Nov 21. Search in Google Scholar

[16] Shi Y, Wang Y, Shao C, Huang J, Gan J, Huang X, et al. COVID-19 infection: the perspectives on immune responses. Cell Death Differ. 2020 May;27(5):1451–4. 10.1038/s41418-020-0530-3Search in Google Scholar PubMed PubMed Central

[17] Liu X, Yang L. Stability analysis of an SEIQV epidemic model with saturated incidence rate. Nonlinear Anal Real World Appl. 2012 Dec 1;13(6):2671–9. 10.1016/j.nonrwa.2012.03.010Search in Google Scholar

[18] Van-den-Driessche P, Watmough J. Reproduction number and sub-threshold endemic equilibria for computational models of diseases transmission. Math Biosci. 2005;180:1–21. 10.1016/S0025-5564(02)00108-6Search in Google Scholar

[19] Kumar S. A new analytical modelling for fractional telegraph equation via Laplace transform. Appl Math Model. 2014 Jul 1;38(13):3154–63. 10.1016/j.apm.2013.11.035Search in Google Scholar

[20] Ghanbari B, Kumar S, Kumar R. A study of behaviour for immune and tumor cells in immunogenetic tumour model with non-singular fractional derivative. Chaos Soliton Fractal. 2020 Apr 1;133:109619. 10.1016/j.chaos.2020.109619Search in Google Scholar

[21] Goufo EF, Kumar S, Mugisha SB. Similarities in a fifth-order evolution equation with and with no singular kernel. Chaos Soliton Fractal. 2020 Jan 1;130:109467. 10.1016/j.chaos.2019.109467Search in Google Scholar

[22] Kumar S, Kumar R, Agarwal RP, Samet B. A study of fractional Lotka–Volterra population model using Haar wavelet and Adams-Bashforth-Moulton methods. Math Methods Appl Sci. 2020 May 30;43(8):5564–78. 10.1002/mma.6297Search in Google Scholar

[23] Kumar S, Ghosh S, Samet B, Goufo EF. An analysis for heat equations arises in diffusion process using new Yang-Abdel-Aty-Cattani fractional operator. Math Methods Appl Sci. 2020 Jun;43(9):6062–80. 10.1002/mma.6347Search in Google Scholar

[24] Veeresha P, Prakasha DG, Kumar S. A fractional model for propagation of classical optical solitons by using nonsingular derivative. Math Methods Appl Sci. 2020 Mar 10;13:65–72. 10.1002/mma.6335Search in Google Scholar

[25] Kumar S, Kumar R, Cattani C, Samet B. Chaotic behaviour of fractional predator-prey dynamical system. Chaos Soliton Fractal. 2020 Jun 1;135:109811. 10.1016/j.chaos.2020.109811Search in Google Scholar

[26] Kumar S, Ahmadian A, Kumar R, Kumar D, Singh J, Baleanu D, et al. An efficient numerical method for fractional SIR epidemic model of infectious disease by using Bernstein wavelets. Mathematics. 2020 Apr;8(4):558. 10.3390/math8040558Search in Google Scholar

[27] Sene N. Analysis of a fractional-order chaotic system in the context of the Caputo fractional derivative via bifurcation and Lyapunov exponents. J King Saud Univ Sci. 2021 Jan 1;33(1):101275. 10.1016/j.jksus.2020.101275Search in Google Scholar

[28] Sene N. Introduction to the fractional-order chaotic system under fractional operator in Caputo sense. Alex Eng J. 2021 Aug 1;60(4):3997–4014. 10.1016/j.aej.2021.02.056Search in Google Scholar

[29] Sene N. Mathematical views of the fractional Chua’s electrical circuit described by the Caputo-Liouville derivative. Rev Mex Fis. 2021 Jan 7;67(1 Jan-Feb):91–9. 10.31349/RevMexFis.67.91Search in Google Scholar

[30] Sene N. SIR epidemic model with Mittag-Leffler fractional derivative. Chaos Soliton Fractal. 2020 Aug 1;137:109833. 10.1016/j.chaos.2020.109833Search in Google Scholar

[31] Kumar S, Chauhan RP, Momani S, Hadid S. Numerical investigations on COVID-19 model through singular and non-singular fractional operators. Numer Method Partial Differ Equ. 2020 Jan 1;2020:29–48.10.1002/num.22707Search in Google Scholar

[32] Kumar S, Kumar R, Momani S, Hadid S. A study on fractional COVID-19 disease model by using Hermite wavelets. Math Methods Appl Sci. 2021 Feb 7;9(7):728–42.10.1002/mma.7065Search in Google Scholar PubMed PubMed Central

[33] Safare KM, Betageri VS, Prakasha DG, Veeresha P, Kumar S. A mathematical analysis of ongoing outbreak COVID-19 in India through nonsingular derivative. Numer Methods Partial Differ Equ. 2021 Mar;37(2):1282–98. 10.1002/num.22579Search in Google Scholar

[34] Khan MA, Ullah S, Kumar S. A robust study on 2019-nCOV outbreaks through non-singular derivative. Eur J Phys J Plus. 2021 Feb;136(2):1–20. 10.1140/epjp/s13360-021-01159-8Search in Google Scholar PubMed PubMed Central

[35] Bonyah E, Zarin R. Mathematical modeling of cancer and hepatitis co-dynamics with non-local and non-singular kernel. Commun Math Biol Neurosci. 2020 Nov 27;2020:1–19.Search in Google Scholar

[36] Zarin R, Khan A, Yusuf A, Abdel-Khalek S, Inc M. Analysis of fractional COVID-19 epidemic model under Caputo operator. Math Methods Appl Sci. 2021 Mar 25;21:123–39.10.1002/mma.7294Search in Google Scholar PubMed PubMed Central

[37] Mandal M, Jana S, Nandi SK, Khatua A, Adak S, Kar TK. A model based study on the dynamics of COVID-19: Prediction and control. Chaos Soliton Fractal. 2020 Jul 1;136:109889. 10.1016/j.chaos.2020.109889Search in Google Scholar PubMed PubMed Central

[38] ur Rahman G, Agarwal RP, Din Q. Mathematical analysis of giving up smoking model via harmonic mean type incidence rate. App Math Comput. 2019 Aug 1;354:128–48. 10.1016/j.amc.2019.01.053Search in Google Scholar

[39] Khan A, Zarin R, Inc M, Zaman G, Almohsen B. Stability analysis of leishmania epidemic model with harmonic mean type incidence rate. Eur Phys J Plus. 2020 Jun;135(6):1–20. 10.1140/epjp/s13360-020-00535-0Search in Google Scholar PubMed PubMed Central

[40] Rakočević MM. A harmonic structure of the genetic code. J Theor Biol. 2004 Jul 21;229(2):221–34. 10.1016/j.jtbi.2004.03.017Search in Google Scholar PubMed

[41] Ahlbom A. Biostatistics for epidemiologists. Florida: CRC Press; 2017 Nov 22. 10.1201/9781315138411Search in Google Scholar

[42] Birkhoff G, Rota C. Ordinary differential equations. New York: John Wiley & Sons; 1982.Search in Google Scholar

[43] Rao VS, Rao PR. Dynamic models and control of biological systems. London: Springer Science & Business Media; 2009 Jul 30. 10.1007/978-1-4419-0359-4Search in Google Scholar

[44] Atangana A, Baleanu D. New fractional derivatives with nonlocal and non-singular kernel: theory and application to heat transfer model. Thermal Sci. 2016 Jan 20;20(2):763–9. 10.2298/TSCI160111018ASearch in Google Scholar

[45] Khan A, Zarin R, Hussain G, Ahmad NA, Mohd MH, Yusuf A. Stability analysis and optimal control of COVID-19 with convex incidence rate in Khyber Pakhtunkhwa (Pakistan). Results Phys. 2021 Jan 1;20:103703. 10.1016/j.rinp.2020.103703Search in Google Scholar PubMed PubMed Central

[46] Taylor AE, Lay DC. Introduction to functional analysis. New York: John Wiley and Sons; 1980. INDICE DEL.; 120. Search in Google Scholar

[47] Ali Z, Shah K, Zada A, Kumam P. Mathematical analysis of coupled systems with fractional order boundary conditions. Fractals. 2020 Dec 10;28(8):2040012. 10.1142/S0218348X20400125Search in Google Scholar

[48] Li C, Zeng F. Numerical methods for fractional calculus. Shanghai: Chapman and Hall/CRC; 2019 Jan 23. 10.1201/b18503Search in Google Scholar

[49] Jajarmi A, Baleanu D. A new fractional analysis on the interaction of HIV with CD4+ T-cells. Chaos Soliton Fractal. 2018 Aug 1;113:221–9. 10.1016/j.chaos.2018.06.009Search in Google Scholar

[50] Baleanu D, Jajarmi A, Hajipour M. On the nonlinear dynamical systems within the generalized fractional derivatives with Mittag-Leffler kernel. Nonlinear Dyn. 2018 Oct;94(1):397–414. 10.1007/s11071-018-4367-ySearch in Google Scholar

Received: 2021-03-02

Revised: 2021-05-15

Accepted: 2021-07-27

Published Online: 2021-11-19

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/phys-2021-0062

Keywords for this article

epidemic model; stability analysis; numerical simulations; sensitivity analysis; COVID-19

Creative Commons

BY 4.0