Fractional order SIR epidemic model with Beddington-De Angelis incidence and Holling type II treatment rate for COVID-19.

In this paper, an attempt has been made to study and investigate a non-linear, non-integer SIR epidemic model for COVID-19 by incorporating Beddington-De Angelis incidence rate and Holling type II saturated cure rate. Beddington-De Angelis incidence rate has been chosen to observe the effects of measure of inhibition taken by both: susceptible and infective. This includes measure of inhibition taken by susceptibles as wearing proper mask, personal hygiene and maintaining social distance and the measure of inhibition taken by infectives may be quarantine or any other available treatment facility. Holling type II treatment rate has been considered for the present model for its ability to capture the effects of available limited treatment facilities in case of Covid 19. To include the neglected effect of memory property in integer order system, Caputo form of non-integer derivative has been considered, which exists in most biological systems. It has been observed that the model is well posed i.e., the solution with a positive initial value is reviewed for non-negativity and boundedness. Basic reproduction number R 0 is determined by next generation matrix method. Routh Hurwitz criteria has been used to determine the presence and stability of equilibrium points and then stability analyses have been conducted. It has been observed that the disease-free equilibrium Q d is stable for R 0 < 1 i.e., there will be no infection in the population and the system tends towards the disease-free equilibrium Q d and for R 0 > 1 , it becomes unstable, and the system will tend towards endemic equilibrium Q e . Further, global stability analysis is carried out for both the equilibria using R 0 . Lastly numerical simulations to assess the effects of various parameters on the dynamics of disease has been performed. © Korean Society for Informatics and Computational Applied Mathematics 2021.

Introduction

The current outbreak of coronavirus disease COVID-19 got reported first from Wuhan, China on 31st Dec 2019. On March 11, 2020 WHO announced that the spread of coronavirus infection COVID-19 would be described as pandemic. On March 13,2020 WHO declared the Europe the center of the Pandemic outbreak. Since March 2020, while new cases appear to be evolved in China, the number of cases in rest of the world has been rising exponentially. Many governments have adopted stringent steps to avoid the spread of this new virus, such as travel bans, demanding social distance and shutting schools, bars, restaurants, and other businesses. The transmission rate is heterogeneous across countries and far exceeds the rate of recovery which allows for a rapid spread. Since the proper treatment facilities are not available for corona virus, therefore the only way to stabilize the outbreak is perhaps to reduce the rate of transmission by limiting exposure of public gathering or social distancing.

The complex behavior of infectious disorders has long been explored as good health is among the most significant concerns in the physical world. The epidemic dynamics are investigated by several mechanisms. A mathematical description of the dynamics helps us to explore the complexity of the disease critically, to analyze statistical and model data, to make predictions about newly emerging diseases like covid 19 and to determine the efficacy of the steps taken.

Kermack & McKendrick [1] developed the first compartmental model for the study of the disease dynamics, many researchers then presented numerous mathematical models such as SIR, SIRS, SEIR, SEIRS etc. (where S, I, E and R represents the susceptible, infected, exposed, and recovered individuals) [2-7].

In recent times, infectious simulation research has been diverted to fractional differential research due to its ability to offer a convincing analysis of certain nonlinear dynamics, the study of fractional order differential equations in the past decades has received a significant attention of researchers.

Since, 2019 novel coronavirus is discovered, many research articles in modeling of COVID-19 transmission have used ordinary differential equations to study its dynamics [8-10]. Though there are reasons that the non-integer models give better fit to the actual data, numerous models have been developed to study the dynamics of COVID-19 using fractional order derivatives [11-14].

Motivated by the above-mentioned literature, we formulated a fractional order SIR model for COVID-19 pandemic with nonlinear incidence rate known as Beddington–DeAngelis type incidence rate and Holling type II saturated treatment rate. Hattaf and Yousfi [15] introduced a Hattaf-Yousfi type of functional response and has the form where measures the force of infection of the disease are non-negative constants and β is a positive constant. This functional response generalizes many functional responses like bilinear incidence rate, saturated incidence rate, Beddington–DeAngelis functional response, Crowley- Martin functional response and Holling type II functional response. To incorporate the effects of measure of inhibition taken by susceptible such as wearing proper mask, personal hygiene and maintaining social distancing and the measure of inhibition taken by infectives can be quarantine, Beddington–De Angelis incidence rate has been chosen [16]. For any outbreak of new disease, treatment capacity is first very low and then grows slowly with improvement of hospital’s conditions, availability of effective drugs, vaccines etc. which are also the conditions arises in case of Covid 19. Since Holling type II function [17, 18] follow the same pattern, therefore it has been chosen as treatment rate for the present study. Hattaf [19] gave a new generalized definition of fractional derivative with non-singular kernel for both Caputo and Riemann–Liouville types. In the present work, fractional order derivative is taken in Caputo-sense. We have evaluated the model mathematically first, and then the simulations have been performed for the model.

The design of the paper shall be as follows: Sect. 2: Some preliminary on fractional order derivatives, Sect. 3: The SIR model of fractional order for COVID-19 transmission is presented, Sect. 4: Non-negativity of solutions, Sect. 5: Analyses of Mathematical Model, Sect. 6: Analysis of local stability behavior, Sect. 7: Analyses of Global stability behavior, Sect. 8: Numerical simulations are provided to support our analytical work, finally conclusion is presented in Sect. 9.

Some preliminaries on fractional order derivatives

We address in this section, some of the basic principles of fractional differential calculus.

Definition 1.

[20] “The Caputo fractional derivative of order of a function is defined as where is the Gamma Function and n is the positive integer such that . When , we have

We have the following definition for the Riemann–Liouville case:

Lemma 1

[21] “Consider the fractional-order system with initial condition where . The equilibrium points of system (1) are locally asymptotically stable if all eigenvalues of the Jacobian matrix evaluated at the equilibrium points satisfy ”.

Formulation of fractional order epidemic model

Several disease models were suggested and evaluated during the last few years. Here we are proposing a fractional order epidemic model for COVID-19 by using the Caputo fractional derivative with order , where . In epidemiological compartmental models, the fractional order is related to the memory or history of the most biological systems. It considers the present or past state information of the system. Also, it describes the memory effect in a cell biological structure. Whole population N(t) is categorized into three classes namely susceptible S(t), infectives I(t) and recovered R(t). We assume that recovered individuals cannot become infected again and, they cannot infect the susceptible population again.

Under these assumptions we proposed a Caputo- fractional order SIR model represented by the following non-linear differential equations: with non-negative initial conditions: and ,where the different model parameters are defined likewise (Table 1):

Table 1

Parameters with their descriptions

Parameter	Description	Unit
β	Coefficient of transmission between susceptible S and infected persons I	Per person per day
λ	Natural mortality rate of people in each class	Per day
μ	Constant rate of recruitment rate of susceptible due to immigration	Person per day
α	Infection mortality rate in infected people	Per day
γ	Rate of recovery in contaminated persons	Per day
β_1	Measure of inhibition (preventive measures taken by susceptibles such as wearing proper mask, maintaining personnel hygiene, maintaining social distance etc.)	Per person per day
β_2	Measure of inhibition (preventive measures taken by infectives such as treatment of infective individuals, quarantine)	Per person per day
σ	Treatment rate of disease	Per day
δ	Rate of limitation in treatment availability	Per person per day

Non-negative solutions: existence and uniqueness

Let and , we need the following lemmas to prove the Theorem 1 on non-negative solutions.

Lemma 2

“Generalized mean value theorem [22] Let and for . Then we have

where , ”.

Lemma 3

[20] “Let Define , where (z) denotes the two-parameter Mittag–Leffler function with parameters and . Then the Laplace transformation of y is given by ”.

Lemma 4

[20] “Let is an arbitrary real number. If then there is a constant such that, for all z in the complex plane, .

Remark 1

[22, 23] “If and for , it is clear from lemma 1 that if then the function is non-decreasing, and if then the function is non-increasing for all each ”.

Theorem 1

There exists a unique non-negative solution for fractional differential equation given by system (2) for and the solution will remain in region for

Proof

First, we check that all the solutions of the model (2) with its initial conditions are non-negative.

Suppose that for . Now we show that for all . We will prove this by contradiction, let for all . Then, there exists a such that for that , at , and for with sufficiently small . Now taking first equation of model (2), we obtain . From Lemma 2, for any , we have . Thus, we get , which is in contradiction to the fact that for . Hence, we have for all .

Next, we prove that for all . Once again, we do it by contradiction. Assume that is not true, then there exists a such that for that , at , and for with sufficiently small . Then we obtain.

.

From Lemma 2, for any , we have . Thus, we get , which is in contradiction to the fact that for . Hence, we have for all . Similarly, we can prove that for all

Thus, it is concluded that all the solutions of model (2) with initial conditions are non-negative. Now we prove boundedness of solutions.

Adding all the equations of the model (2), we obtain:

where . As , we have

Now consider the initial value problem , . Using comparison principle [24], we obtain the following inequality:

for all . Now applying Laplace transform to the initial value problem we obtain

Using Lemma 3, we obtain

Applying inverse Laplace transform in the above two equations, we get

, using we have

,

By Lemma 4, we obtain where is constant given in Lemma 4. Hence, as , we have with . Thus, the solutions are bounded and will remain in region M for Therefore, Theorem 1 is proved and solution remains in .

Mathematical model analyses

In this section, disease-free equilibrium and endemic equilibrium are first identified after which the stability results are observed for both the equilibriums. Thereafter the Basic Reproduction number is determined using next generation matrix which tells about persistence and removal of the disease. Finally, we have discussed the existence and stability of equilibrium points along with their local stability behavior.

Equilibria and their stability

As in the above non-linear model of three equations, for stability analysis it is logical to consider only the first two equations because these are independent of the third equation. The model is therefore reduced to: where .

To find the equilibrium points, set the right side of the system (3) to zero, we get

Thus, solving the above system, we obtain the two equilibria points called disease free equilibrium point (DFE) and endemic equilibrium point (EEP) specified as

DFE: showing no illness in the environment and all the people are susceptible only.

EEP: which is described in the later section.

Now, to obtain the behaviour of stability of these two equilibria, Basic Reproduction Number is required to be computed [25].

Determination of basic reproduction number

To compute the Basic Reproduction number using the next generation matrix method, we assume that where , represents the matrix of new infections coming and represents the matrix of transfer of individuals entering and leaving from the compartments.

The Jacobian matrices of and are evaluated at disease free equilibrium point and are given by

Now we need to compute the inverse of and then given by

This matrix is known as next generation matrix and the spectral radius of this matrix is the basic reproduction number (Ye & Xu, 2019). for our model is

The basic reproduction number is the average number of secondary cases caused by an infective in a fully susceptible population. The value of depends on both the disease and the host population.

Analysis of local stability behaviour

Here in this section the analysis of local stability behaviour of both the equilibrium points DFE and EEP are carried out.

Local stability analysis of DFE ()

Now we are concerned about the local stability of disease-free equilibrium point (DFE ). We claim the following theorem for this and prove this:

Theorem 2

= , the disease- free equilibrium of model (3) is locally asymptotically stable when the Basic Reproduction number  < 1 and unstable otherwise.

Proof

The local stability behaviour of the DFE point of the non-linear system (4) is computed by linearizing the system around the DFE (). Thus, we have obtained a linearized matrix of the system at as shown below

Now the characteristic equation of the Jacobian matrix is computed as where and .

First, we will prove that the roots () of the characteristic Eq. (6) which are given as are real

Therefore, both the roots and are real. We now show that when both the roots are negative.

Sum of the roots

Product of the roots .

Thus (when ).

and (when ).

Hence, we conclude that the roots of the characteristic Eq. (6) are negative when . Thus, alternately we can say that all the eigen values of the Jacobian matrix are of negative sign when Consequently all the roots of the Eq. (6) have negative real parts, thus, by the fractional Routh-Hurwitz criteria [21, 26] all roots follow:

|> .

The above Fig. 1 interprets Theorem 2 using the values given in Table 2.

Fig. 1

Basic Reproduction number vs Infected Population

Table 2

Model parameters with their values

S.No	Parameters	Value
1	β	0.003
2	λ	0.05
3	μ	5
4	α	0.06
5	γ	0.002
6	β_1	0.002
7	β_2	0.001
8	σ	0.02
9	δ	0.002

EEP () existence and study of local equilibrium

The conditions for the existence of the endemic equilibrium () are defined here and the endemic equilibrium’s local stability behavior is discussed in this section.

After solving the Eqs. (7) we obtain:

and is the root of the following equation: where the coefficients are given by

Theorem 3

For >1 there exist a unique endemic equilibrium of model (3).

Let  > 1

Considering Eq. (8) which is third degree polynomial given as

. We observe that the coefficient of , the leading coefficient, is negative. Therefore , also and for  > 1, is a continuous function of

Thus, the fundamental theorem of algebra implies that there exist at most three real positive roots of . Here, we examine only the case of a particular endemic equilibrium. Assuming  > 1 and noting that and , under the following signs of and a unique endemic equilibrium exists:

,

,

,

If any of the criteria (i)-(iii) is fulfilled, then there exist a unique , from which we can find the value of So, unique endemic equilibrium  =  has been demonstrated for  > 1

The local stability of endemic equilibrium is now examined here. To do so, we linearize the model (4) around endemic equilibrium and thus obtain the Jacobian matrix as shown below:

The characteristic equation of is given by

where

It can readily be seen that Jacobian matrix ’s eigen value have negative real parts if and only if and . Also, and if .Hence, by fractional Routh-Hurwitz criteria [21, 26] all the roots of Eq. (9) have negative real parts and satisfy the condition . This proves the following theorem.

Theorem 4

= , the endemic equilibrium of model (3) is locally asymptotically stable when the Basic Reproduction number  > 1 and unstable otherwise.

Analysis of global stability behaviour

Here in this section the analysis of global stability behaviour of both the equilibrium points DFE and EEP are carried out.

Global stability analysis of DFE () and Endemic equilibrium ()

Now we are concerned about DFE’s and EE’s global stability. We claim the following theorem for this and prove this:

Theorem 5

The disease- free equilibrium  =  of model (3) is globally asymptotically stable when .

We define a Lyapunov function by

Now we have

Using the value of in above equation, we have

This concludes that if , then we have Furthermore, we know that if and only if and . Thus the largest invariant set for is the singleton set {}, where and also all the solutions in converges to in accordance with the LaSalle’s invariance principle [25, 27–29]. So, is globally asymptotically stable when Hence Theorem 5 is verified.

Theorem 6

The endemic equilibrium  =  of model (3) is globally asymptotically stable when  > 1.

We define a Lyapunov function by

Now we have

Using the endemic conditions,

We get

This implies,

Therefore, . Also, the largest invariant set for is the singleton set { and all the solutions in converges to in accordance with the LaSalle’s invariance principle [25, 27–29]. Therefore, is globally asymptotically stable when  > 1. Hence Theorem 6 is verified.

Numerical simulation

In this section, to support our analytical work we performed numerical simulations using MATLAB 2012(b). The initial conditions that are used in the simulation are S (0) = 85, I (0) = 0 and R (0) = 0 [30] and the required parameters are given in Table 2. Parameters , , are estimated and rest are chosen from paper [30]. Using the predictor–corrector method [31], system (2) is simulated numerically by MATLAB for various values of the fractional order ranging from 0.4 to 1 to observe the dynamical behavior of susceptibles, infectives, and recovered population with changed order of derivative.

Figure 2 demonstrates the dynamic phenomena of susceptible class for various values of and it is found that susceptible are significantly influenced by decreasing the fractional order value. It has been seen that the increased value of fractional order is one of the reasons which held responsible for the increment in susceptible population in the community. Also, memory effect can be observed for different values of , for there is no memory effect on susceptibles and by lowering down the value of the memory affect becomes visible which is the fractional order derivatives property.

Fig. 2

Result of changing fractional order on susceptibles

Figure 3 indicates the pattern of infected persons for various fractional order values. From this figure it is noticed that there is a lower number of infected individulas for small order, as the order raises the number of infected individuals also increases and and then it starts decreaing as the fractional order decreases after hitting the peak value of I(t) = 68 for = 1. Also, we can observe the presence of memory effect because as we decrease the value of , the infected people are taking more time to recover.

Fig. 3

Result of changing fractional order on infectives

In addition, from both the figures (Figs. 2 and 3) it is seen that system converges to its steady state more rapidly as the order of fractional order derivatives decreases, which is in accordance with these papers (Kumar, 2020; Naik, 2020).

Figure 4 reveals the population of people recovered for various fractional orders and it is observed that the number of those recovered people is raised as we increase the order of fractional derivatives.

Fig. 4

Result of changing fractional order on recovered

Figures 5(i), (ii) and (iii) shows behavior of susceptibles on varying preventive measures taken by susceptibles for different values of fractional order and no measures of inhibition are adopted by infected individuals. It infers those preventive measures taken by susceptibles will be affecting the transmission of susceptibles into infected compartments. More the number of individuals following the preventive measures implies a greater number of susceptibles and lesser number of infectives. Despite contact occurring between susceptible and infective individuals, if preventive measures are taken by susceptibles, then the transmission from susceptible compartment to infective compartment is reduced largely.

Fig. 5

(i) Effect of varying measures taken by susceptibles on susceptibles with and no measures taken by infectives, (ii) Effect of varying measures taken by susceptibles on susceptibles with and no measures taken by infectives, (iii) Effect of varying measures taken by susceptibles on susceptibles with and no measures taken by infectives

In Fig. 5(iii) for fractional order = 1, an inverted peak in the number of susceptible population can be noticed, which progressively fades as the order of derivative reduces, which is due to the memory effect or the people’s encounter with the disease. For = 0.001, number of susceptible peoples were 76 for = 0.6 at time t = 50 days (Fig. 5(i)), while they were 58 and 43 respectively for = 0.8 (Fig. 5(ii)) and = 1 (Fig. 5(iii)). It is indicating that for the same as the fractional order decreases, the number of susceptibles increases ultimately number of infected populations decreases which is due to the people’s memory effect about the disease.

Figures 6(i), (ii) and (iii) shows the behavior of infected individuals by varying the level of preventive measures taken by susceptibles () when no inhibitory measure is taken by infected people for different values of fractional order and it is observed that on increasing the value of preventive measures by susceptibles number of infectives starts decreasing. Also, the memory effect can be seen by comparing these three graphs as the fractional order decreases the peak in the infectives gradually disappears.

Fig. 6

(i) Effect of varying preventive measures taken by susceptibles on infectives with and no measures taken by Infectives, (ii) Effect of varying preventive measures taken by susceptibles on infectives with and no measures taken by Infectives, (iii) Effect of varying preventive measures taken by susceptibles on infectives with and no measures taken by Infectives

Figures 7(a)(i,ii,iii) and 7(b)(i,ii,iii) shows the behavior of infectives by changing the measures of inhibition taken by infectives for various values of fractional order. In 7(a)(i,ii,iii) no preventive measures are taken by susceptibles while in 7(b)(i,ii,iii) susceptibles take preventive measures. It is evident that at any time if enough preventive measures are taken by susceptibles then number of infectives are less in comparison to when no preventative measures are taken by susceptibles. Also, it is depicted that as the preventive measures taken by infectives increases, number of infected populations starts decreasing. In Fig. 7a(ii) at time t = 50 number of infected individuals were 24 while in 7b(ii) they were 21 at same time for = 0.001 which shows that if susceptibles use enough preventive measures then their chances of entering the infected compartment reduces. Initially for lower order of fractional order no peak is there but as the fractional order increases, peak in the increasing infected population starts appearing which shows the effect of memory property.

Fig. 7

a(i) Effect of varying measures taken by infectives on infectives with and no measures taken by susceptibles, (ii) Effect of varying measures taken by infectives on infectives with and no measures taken by susceptibles, (iii) Effect of varying measures taken by infectives on infectives with and no measures taken by susceptibles, b(i) Effect of varying measures taken by infectives on infectives with , (ii) Effect of varying measures taken by infectives on infectives with , (iii) Effect of varying measures taken by infectives on infectives with

Figures 8(a)(i,ii,iii) shows the effect of change in preventive measures taken by infectives on susceptibles when no preventive measures were taken by susceptibles while in Figs. 8(b)(i,ii,iii) preventive measures were taken by susceptibles. First, if infected individuals would get proper and timely treatment then the rate of spreading the infection decreases due to which susceptible population also decreases. In Fig. 8(a)(iii) at time t = 50, number of susceptibles were 51 with = 0.001 and no safety measures were used by susceptibles while in Fig. 8(b)(iii) at the same time susceptible population is 45, here they have used the preventive measures with same . Second, when 8(a)(i,ii,iii) is compared with 8(b)(i,ii,iii) it is inferred that if susceptible population practice more preventive measures like wearing masks, keeping social distance, and maintaining sufficient hygiene, their risk of entering infected compartment decreases. When comparing graphs, the memory effect is immediately evident.

Fig. 8

a(i) Effect of varying measures taken by infectives on susceptibles with and no measures taken by susceptibles, (ii) Effect of varying measures taken by infectives on susceptibles with and no measures taken by susceptibles, (iii) Effect of varying measures taken by infectives on susceptibles with and no measures taken by susceptibles, b(i) Effect of varying measures taken by infectives on susceptibles with , (ii) Effect of varying measures taken by infectives on susceptibles with , (iii) Effect of varying measures taken by infectives on susceptibles with

Figures 9(a)(i,ii,iii) shows the effect of lockdown on infected population. Because of the lockdown the rate of transmission is reduced leading in a drop of infected population. Also, it is observed that for lower order of fractional order there is no peak but peak starts appearing as increases which shows the presence of memory effect. Similar effect of lockdown can be observed for susceptible population as shown in Figs. 9(b)(i,ii,iii), because of lockdown interaction will be less and the speed of transmission will be reduced, and a smaller number of people will be infected.

Fig. 9

a(i) Effect of lockdown on infectives with , (ii) Effect of lockdown on infectives with , (iii) Effect of lockdown on infectives with , b(i) Effect of lockdown on susceptibles with , (ii) Effect of lockdown on susceptibles with , (iii) Effect of lockdown on susceptibles with

Figures 10(i) shows the effect of immigration on infectives for different values of recruitment rate. As the pace of recruitment rate rises, so does the number of infectives; however, after reaching its peak, the number of infectives begin to decline and settles down into a steady state, likewise Fig. 10(ii) describes the effect of immigration on susceptibles and same impact can be observed on susceptibles i.e., the higher the number of immigrants, higher the number of susceptibles and after a period the number of susceptibles archives the steady state. Figure 11 demonstrates that when the transmission rate increases, so does the reproduction number, implying that the infected population will increase. Figure 12 illustrates that, the more the preventive actions like as vaccination, masking, ensuring social distance and personal hygiene are adopted by susceptible population, the lower the number of infected cases and thus reproduction number will be reduced.

Fig. 10

(i) Effect of immigration on infectives, (ii) Effect of immigration on susceptibles

Fig. 11

Effect of transmission of disease on basic reproduction number

Fig. 12

Effect of measures taken by susceptibles on basic reproduction number

Conclusion

This paper introduces and discusses a novel Caputo-fractional SIR model with Beddington–De Angelis incidence rate for COVID-19 and Holling type II type saturated treatment rate. The system (2) experience two equilibriums: disease-free equilibrium and endemic equilibrium. The , basic reproduction number studied the analysis of local stability for both the equilibria. The equilibrium of disease-free state has proven to be stable for i.e., disease diminishes and for , the disease-free equilibrium becomes unstable and the endemic equilibrium becomes stable under some conditions. Further, the global stability behavior of both the disease-free equilibria and the endemic equilibria is analysed and it is found that DFE is globally asymptotically stable when and the endemic equilibrium is globally asymptotically stable when.

Moreover, numerical results for the favour of the analytical work are obtained. It is noted that the variation in the values of non-integer order derivatives does not impact the equilibrium stability of the system but the time to achieve steady conditions or endemic equilibrium decreases i.e., better convergence can be accomplished by reducing the order of fractional derivatives or greater the values of slower is the convergence thus, the utility of the non-integer order derivative rather than the integer order model has been investigated. The rate of disease transmission will decrease if there is less contact between susceptibles and infectives, which could be achieved by taking adequate safety precautions both by susceptibles and infectives, as well as imposing lockdown. Moreover, if disease transmission is reduced, , the basic reproduction number declines. With this fractional order model and taking the authentic data of COVID-19, the dynamics of COVID-19 can be predicted.