Complex dynamics of a fractional-order SIR system in the context of COVID-19.

This paper proposes and analyses a new fractional-order SIR type epidemic model with a saturated treatment function. The detailed dynamics of the corresponding system, including the equilibrium points and their existence and uniqueness, uniform-boundedness, and stability of the solutions are studied. The threshold parameter, basic reproduction number of the system which determines the disease dynamics is derived, and the condition of occurrence of backward bifurcation is also determined. Some numerical works are conducted to validate our analytical results for the commensurate fractional-order system. Hopf bifurcations for the fractional-order system are studied by taking the order of the fractional differential as a bifurcation parameter. © Korean Society for Informatics and Computational Applied Mathematics 2021.

Introduction

Disease is an integral part of modern civilization. As long as human lives, illnesses will last in our society. But, the extent of these diseases can be reduced, and even we can eradicate the disease from the population. Among various types of diseases, infectious diseases are the most vulnerable diseases. Any contagious disease, caused by the pathogen of bacteria, viruses, fungi, or some parasites can transmit from an infected human to a susceptible human either through direct contact (viz. influenza, tuberculosis, rubella, measles, HIV-AIDS, etc.) or via some medium (viz. malaria, dengue, chikungunya, etc.). Among all infectious diseases, the recent pandemic COVID-19 caused by the virus SARS-CoV-2 (severe acute respiratory syndrome coronavirus 2), the first pandemic of the twenty-first century has affected the lifestyle of human population of every country throughout the globe. The disease COVID-19 was first noticed in the Wuhan city of the Republic of China in December of last year [38]. Till today except for very few countries, COVID-19 pandemic has worsened the situation all over the globe [37]). It should be noted that till May 2020, more than six million world populations have been affected by the disease COVID-19, where almost 50% closed cases, including 12% death (see [37]).

The rate of progression of infectious diseases including COVID-19 depends on the amount of pathogen living within a host, its rate of growth, and the interaction with the host’s immunity. A mathematical model can perfectly capture these relations. Mathematical models help to represent all the scenarios in a compact form. It also captures how disease propagates, its long-term behavior, and ultimately helps to determine possible interventions to restrain the infectious disease. The mathematical model of an infectious disease was first introduced by Kermack and McKendrick [19]. They have used a simple SIR model to formulate the mathematical model. After this theoretical epidemiology has come a long way. Several works on control strategies like vaccination, treatment, isolation were taken into consideration by many researchers [17, 20, 24]. But, in all the works done so far, researchers have considered the integer order system while constructing the mathematical model. Recently some researchers have also developed some model-based works to examine the detailed dynamical behavior of the COVID-19. The researchers including Ahmed et al. [3], Kucharski et al. [21], Prem et al. [32], etc. have also developed some mathematical models and presented rigorous analysis. In an article, Ribeiro et al. [33] have used some stochastic-based regression models to forecast the phenomena in ten most affected states of Brazil. Noting that isolation is a good measure to control the disease, Hellewell et al. [14] have used the operative techniques of COVID-19 diseases using separation. In the recent works, Mandal et al. [25, 27] have constructed and analyzed a mathematical model on the COVID-19 in the pandemic scenario of three states of India. In other recent work, a balance between lockdown and compliance to investigate the COVID-19 scenes has been proposed by Zegarra et al. [2]. Some more theoretical works on COVID-19 can be found in [25, 27, 34].

The shortcoming of the integer order system is that it does not rely on the previous history of the system. But, when a disease transmits, the susceptible population uses their memory to prevent infection. A dynamical system involving fractional-order derivatives bears the information regarding its present and also its past states [12]. Hence, fractional-order systems give more in-depth knowledge about the system than integer order systems. Also, as there are many developments in fractional calculus (See the articles [1, 4, 28, 29]), so, it is logical to apply a fractional-order differential equation to represent a biological phenomenon through mathematical modeling. According to Du et al. [10], if the parameter , the order of the fractional-order system tends to 0, the model has full memory, and when , the model system is known as memoryless. Hartley et al. [13] used fractional differentiation for the first time for a physical system in demonstrating the behavior of the semi-infinite lossy (RC) line. Torvik et al. [35] have studied that fractional derivatives appear naturally to describe some properties of motions of a Newtonian fluid. At present, there are several works involving the application of fractional-order systems in various fields including epidemic models [15, 26].

In human society, any epidemic can not be thought without the effect of memory. Whenever an infectious disease propagates in human society, the knowledge or experience of each people about that particular disease should affect their response [36]. People apply several precautions (e.g. vaccination, if it is available) for a disease if they have a past experience of the disease. In this regard, some important control measures can suppress the spreading of the disease. But, mere knowledge about a disease may not protect all the time. That is why people are willing to adopt some new paths to control these diseases. The memory of an infectious disease from the previous situation has more contribution on the recent situation as compared to no memory. It is anticipated that the long-time impact of memory decreases in future time more tardily than exponential decay, but generally, its behavior is a similar way of power-law damping function.

In this work, we have formulated a new SIR type epidemic model with the help of fractional differential equations. Although there are many ways to define a fractional-order differential equation in the sense of Riemann-Liouville, Grnwald-Letnikov, and the Caputo. In this work, we use the Caputo definition due to its similarity in using the initial conditions with the integer-order differential equation. From this perspective, we develop a mathematical model considering a system of fractional-order differential equations and analyze the thorough dynamics of the proposed model. Moreover, we apply this model to predict the COVID-19 situation of India after the second wave.

The whole manuscript is structured in the following way: An SIR type model with the fractional differential equation is formulated in Sect. 2. The equilibrium points, existence and uniqueness of solution, their dynamical behavior, stability of the system, and the occurrence of Hopf bifurcation are discussed in Sect. 3. Several numerical simulations are performed in Sect. 4, and the application of the model on some real-world data is discussed in Sect. 5, and lastly, in Sect. 6, we present a concise discussion and conclusion.

Model formulation

In this section, a new SIR type model is considered. We divide the whole population into susceptible S(t), infected I(t), and recovered R(t), the three mutually exclusive time-dependent classes. First, we consider that at any time t,  A be the newly recruitment rate, be the disease transmission rate, d be the natural death rate, be the rate of death due to the disease, be the natural recovery rate. We assume that some infected person will be recovered due to available treatment control. Let v denote the treatment control parameter for the infected population. The rate of recovery is considered as with representing the effectiveness of the treatment due to saturated type treatment control (see Jana et al. [16]). Every country, including the most developed countries in the world, has limitations in providing medical resources like proper medicines, availability of beds in hospitals, etc. Particularly, the contagious disease like COVID-19, which has appeared all of a sudden, it is not possible to provide treatment to every patient when the number of patients becomes very large. From this point of view, we have used the saturated type recovered rate due to treatment. A schematic diagram of the above assumptions is described in Fig. 1. Based on these assumptions, we define our model in the following manner:

Fig. 1

Flow chart of the disease transmission

with initial conditionsAccording to Du et al. [10] a memory process usually consists of two stages: the fresh stage and the working stage. The fresh memory is used to just remember things while the working memory helps to do cognitive works. In the fractional model, we capture the working memory. The critical point between the fresh stage and the working stage is usually not the origin. This observation is quite different from the traditional fractional models of one stage. For example, the fractional Maxwell’s model is a one-stage model. As the combination of two simple models, it has a more complicated expression than Eq. (2.1). We also find that the order of fractional. Also, in reality, no individual has the same memories, it changes from person to person. Once infected an individual gains knowledge about the disease which helps to obstruct the transmission of infectious disease. We consider this fact into the fractional form of the system (2.1). Here we assume that , and are such properties corresponding to the susceptible, infected and recovered population respectively, where . Whereas, the recovered will have better working memory than infected one but less working memory than susceptible. So, the fractional order form of the system (2.1) takes the formHere is the fractional-order derivative in Caputo sense [30, 31] for . In the system (2.2) we replace the parameters of the system (2.1) in such a way that the problem of dimension mismatch is solved.

For simplicity, we set .

Using these transformations into the system (2.2), we havewith initial conditions . Here is the fractional-order derivative in Caputo sense [30, 31] for . The system (2.3) is called to be an incommensurate fractional-order system. If , then the system is said to be a commensurate fractional-order system which is given bywith initial conditions . Here is the fractional-order derivative in Caputo sense [30, 31].

The analysis of fractional-order system

At this stage, we analyze the dynamical behavior of the fractional-order system (2.4).

Equilibrium points

Now we discuss the existence of nonnegative equilibrium points of the system (2.4). We observe that the system (2.4) has three equilibrium points. One is disease free equilibrium (DFE) and the other two equilibrium points are endemic equilibrium (EE) if they exist and where and are the roots of the following quadratic equationwhere and

Here the parameter is recognized as the basic reproduction number. Hence we can derive and .

Let be the discriminant of the Eq. (3.1).

If , (i.e. if ) then (3.1) has a unique positive root, . Thus we can say that the system (2.1) has a positive unique EE if .

If , then which implies . From this we can get a positive EE of the system (2.4) as if and only if . Now, gives the treatment control parameter to be .

If , (i.e. if ) and , then the system (2.1) possesses no EE. Also if and i.e. then there is no endemic equilibrium points.

We formulate the next theorem which is grounded on the above discussions.

Theorem 3.1

(i) The system (2.4) possesses two EE points if , , and no EE point if , , .

The system (2.4) possesses two equal EE points if , , and no EE points if and .

If the basic reproduction number, , then the system (2.4) possesses a unique EE point if i.e. , where v is the control parameter and no EE points if .

The system (2.4) possesses a unique EE point if .

From the above theorems we experience that the system (2.4) has two different equilibrium points even when . Hence we can conclude that the proposed model experiences a backward bifurcation at which shown in Fig. 2. Thus the condition is not sufficient to eradicate the disease from the population.

Fig. 2

Backward bifurcation of the system (2.4) with the parameter values , , , , , , and

Theorem 3.2

[5] The system (2.4) passes through a backward bifurcation at if and only if where is determined in subsequent steps.

Here the treatment control parameter v is selected as the backward bifurcation parameter. Observe that if the case (iii) of theorem (3.1) is true then the system (2.4) undergoes a backward bifurcation at . Hence it can be explained that the system possesses two EE point in some interval , where and . can be derived from i.e. . Now, gives where and Solving we have Now, this case occurs only when and . Thus we get,From the above conditions we have derived at the critical value of v as Thus the new threshold value at which the diseases eradicates from the population (say) is Therefore, depending on the above discussion we can state the next theorem.

Theorem 3.3

The system (2.4) holds

a unique EE point if and the disease will remain in the population.

two EE points if which implies that the system passes through a backward bifurcation.

a unique DFE point if which implies that the disease will die out.

Existence and uniqueness of the solutions

We now state and prove the following lemma which helps us to conclude that our constructed model has a unique solution.

Lemma 3.4

[22] Consider the fractional-order systemwhere If the Lipschitz condition with respect to x is satisfied by the function f(t, x), then the system (3.2) has a unique solution on the interval .

Theorem 3.5

The system (2.4) has a unique solution in the region where, , and .

Proof

We prove the existence and uniqueness criterion for the fractional-order system (2.4). For that purpose we consider the region where . Let and be arbitrary two points in the region and define the mapping by , whereFor any , we haveTherefore, the function H(X) fulfills the criteria of the Lipschitz’s condition with respect to . Then by using the Lemma 3.4, we can conclude that our system (2.4) has a unique solution with initial conditions .

Nonnegativity and boundedness of solutions

In mathematical biology, our main concern is on the nonnegative and bounded solutions. Hence, we try to find the solutions that possess the (i) non-negativity and (ii) bounded properties. For this purpose, we consider the following region .

Theorem 3.6

Every solutions of the fractional system (2.4) whose initial values start in i.e. are all nonnegative and bounded.

To prove the nonnegativity of the solution, we use Lemma 3.2 [15]. We have discussed the existence of solution of the system (2.1) in Theorem 3.1. As (2.4) is a homogeneous system of equations with initial conditions we can conclude that the solutions of (2.4) is non-negative. Furthermore, the solutions are nondecreasing in t. Now to prove the uniform boundedness of the solutions, we consider a mapping . Then we have,Using the Lemma 3.2 [15], we have . Therefore, all the solution of the fractional system (2.4) whose initial conditions start in the region is bounded and it lies in the region .

Dynamical behavior

In this portion of the research article, we study the stability conditions of each equilibrium point. Here, we assume that the fractional system (2.4) possesses unique EE point. The local stability investigation of the system (2.4) is done with the help of making linearization around each equilibrium point. The Jacobian matrix of the fractional system (2.4) at the DFE point isThe eigenvalues of the above Jacobian matrix are . Therefore, , where and if i.e. if . Again we have established that the system undergoes backward bifurcation iff Hence, by the Lemma 3.3 [15], we can conclude that the system is asymptotically stable around the DFE if . Thus, we can state the next theorem.

Theorem 3.7

The DFE of the fractional system (2.4) is locally asymptotically stable iff

Next, we discuss the local stability criteria of EE point . J, the Jacobian matrix at of the system (2.4) is given byOne eigenvalue of the matrix is and the other two eigenvalues are obtained from the matrixThe characteristic equation of the Jacobian matrix at the EE, is given bywhere , and

Let be the eigenvalues of the matrix . Now, will be negative ifThen the eigenvalues arePreviously, we have determined one eigenvalue whose argument i.e. , for . Now, we consider different cases considering different values of . From the above discussions we may state the next theorem about the stability of the EE point.

If and , then both the eigenvalues are negative. So, . Hence from Lemma 3.3 and Lemma 3.4 of the article by Jana et al. [15], the EE is asymptotically stable.

If and , then one of the eigenvalues will be non negative. Then, for i=2 or, 3, and . Hence, by Lemma 3.3 and 3.4 [15], the EE is unstable.

If the conditions and hold, then both the eigenvalues will be complex conjugates. , where . Therefore, the EE point will be asymptotically stable if . Hence, for the order of differentiation the EE is locally asymptotically stable.

If and , then . Therefore, EE will be asymptotically stable if . Hence, for the order of differentiation

If and , then . Then EE is locally asymptotically stable.

Theorem 3.8

The conditions for the local stability of EE of the system (2.4) are followed by:

The EE is locally asymptotically stable if and .

The EE is locally unstable if and .

The EE is locally asymptotically stable if and .

The EE is locally asymptotically stable if and .

The EE is locally asymptotically stable if and .

Global stability of DFE

We now demonstrate the asymptotic global stability of the DFE point.

Theorem 3.9

The DFE point of the fractional system (2.4) is globally asymptotically stable if

We prove the theorem using the Lemma 4.3 of the article [15]. Consider a function in the following way:   .

Then the solution of the fractional system (2.4) and the constructed function V together imply thatThe Lemma 4.3 [15] suggests that any solution of (2.1) starting in the region tends to largest invariant set . Hence, , the fractional system (2.4) changes to the equationsThe solution of (3.6) is

Then we have as . Further, we have already established that the DFE is locally asymptotic stable if Combining both results, we may conclude that the DFE point of the limit set (3.6) is globally stable asymptotically for . Hence the proof.

Analysis of Hopf bifurcation

The existence and analysis of Hopf bifurcation for a system is modeled by ordinary differential equations is studied by many researchers [17, 18]. In this portion, we analyze Hopf bifurcation of the fractional-order system (2.4) for different values of . Here we assume the effectiveness of the saturated treatment control parameter and the fractional-order as bifurcation parameters.

It is known that Hopf bifurcation may take place in spite of stability of a system at the critical value of bifurcation parameter for an integer order system. Therefore, all the roots of the Eq. (3.3) are all real. Hence, delivers us the necessary critical value of the bifurcation parameterIn this case, the Hopf bifurcation occurs at if the succeeding conditions are fitted i.e.Now we state some necessary conditions for the existence of Hopf bifurcation in our fractional system.

Here we take as the bifurcation parameter. Using the Lemma 3.4 [15], we can say that stability of the fractional system (2.4) depends on sign of . Based on this observation, we define the function byThe Lemma 3.4 [15] implies that the EE point will be stable asymptotically if and it will be unstable if .

Theorem 3.10

[23] If the bifurcation parameter crosses the critical value , then the fractional system (2.4) passes through a Hopf bifurcation at EE point , if the given singularity conditions [(i), (ii)] and also the transversality condition [(iii)] are satisfied.

the Jacobian J at the EE point of the system (2.4) induces complex conjugate eigenvalues with .

and

.

Incommensurate fractional-order model

Let us assume the incommensurate fractional system (2.4), and where and , and , are relatively prime for . Then by Lemma 3.5 [15], an equilibrium point E(S, I, R) is locally stable asymptotically iff for all the roots ’s of the characteristic equationwhere and J, Jacobian matrix evaluated at the same equilibrium point E(S, I, R).

Numerical experiments

In the numerical experiments, we apply predictor-corrector P(EC)E (Predict, multi-term(Evaluate, Correct), Evaluate) method which is the modified method of PECE (Predict, Evaluate, Correct, Evaluate) and Adams-Moulton algorithm [8, 9, 11] in the environment of MATLAB-16 software. For this purpose, we use the solver function as implicit fractional linear multi-step methods (FLMMs) for fractional-order systems. To perform numerical works, we choose the set of parameters asThese parameters yield the EE of the fractional system (2.4) as where . exists if which implies that the feasible region of is . The EE of the ordinary system will be stable if and gives . For , the eigenvalues contain negative real parts and hence . Therefore the EE point is asymptotically stable and converges to for all . If , for , the EE is unstable as the eigenvalues contain positive real parts. Now, we depict some figures for the solutions of the fractional system (2.4) in the presence of the initial conditions . If and , then it can be said that the EE point is stable by Theorem 3.8. If we choose then the eigenvalues are . When , EE point will be asymptotically stable and it will be unstable whenever . These behaviors are shown in Fig. 3 for , Fig. 4 for , and Fig. 5 for .

Fig. 3

Time series plot of S(t), I(t) and R(t) for of system (2.4) corresponding to the orders

Fig. 4

Time series plot of S(t), I(t) and R(t) for of system (2.4) corresponding to the orders

Fig. 5

Time series plot of S(t), I(t) and R(t) for of system (2.4) corresponding to the orders

Next we assume , the eigenvalues are . When the value of , EE point will be asymptotically stable and for the EE point will be unstable. Figure 6 depicts the behavior.

Fig. 6

Time series plot for of system (2.4) corresponding to the order respectively

From the above discussion it can be said that with changing values of , we can get different region of for which the EE points remain asymptotically stable. The fractional order system will be asymptotically stable if . Here, . So, is a decreasing function. The stability region in terms of are shown in Fig. 7.

Fig. 7

Stability region of the commensurate system (2.4)

Now we demonstrate Hopf bifurcation for the fractional system (2.4). To satisfy the singularity condition (i) of Theorem 3.10, we must have and of the equation (3.5). We see that Hopf bifurcation of the system (2.4) occurs for if the parameter representing the effectiveness of treatment control lies in the interval (0, 1.5). We choose then the critical value of as is determined from condition (ii) of Theorem 3.10 and we get the critical value . The third condition holds good since, at . Hence, we see that Hopf bifurcation undergoes at and with increasing value of from the orbits of the limit cycles also increase. Figure 8. depicts such behavior. The eigenvalues at EE of the incommensurate fractional system (2.3) can be determined by the help of the equation (3.7) which is, in general, a higher-order polynomial based on the values of . Now we choose . Then we obtain the EE point (10.89, 0.703, 2.769). We set the fractional orders as and obtain the eigenvalues from equation (3.7). We see that the absolute value of the minimum of arguments of the eigenvalues is . Hence, in this case, the eigenvalues are unstable. Again we choose and find the absolute value of the minimum of argument of the eigenvalues of the equation (3.7) as . Hence, in this case, the EE point is asymptotically stable. These behavior are presented in Fig. 9.

Fig. 8

Limit cycles for at

Fig. 9

Time series plot for of system (2.3) corresponding to the order and

Now we state the theorem for existence criteria of Hopf Bifurcation due to the parameter .

Here we define a function byFrom Theorem 2 in [7], we can say that the EE point will be stable asymptotically if and it will be unstable if . Now we state the theorem for existence criterion of Hopf bifurcation about the parameter .

Theorem 4.1

[23] If the bifurcation parameter crosses its critical value (say), then the system (2.4) passes through a Hopf bifurcation at the EE point , if the following singularity conditions (i) and (ii) and the transversality condition (iii) are satisfied.

the Jacobian J computed at the EE point of the system (2.4) determines a pair of complex conjugate eigenvalues with .

and

.

Numerical simulation of Hopf bifurcation

We now choose as a bifurcation parameter. Here also we obtain the eigenvalues from the equation (3.5). The first condition of Theorem 4.1 implies that the eigenvalues obtained from equation (3.5) must contain positive real parts. Then we must have and . Hence, a Hopf bifurcation of the system (2.4) occurs if

The critical value of is obtained from and which impliesLet us choose the order of fractional derivative then the equation yields the critical value of the bifurcation parameter as and at this point . Therefore, the transversality criteria [(iii)] is also satisfied. Thus, the existence of Hopf bifurcation confirms at and with increasing value of from the orbits of the limit cycles also increase. Figure 10. depicts such behavior.

Fig. 10

Limit cycles for at

Hopf bifurcation of the incommensurate system

Here we consider the incommensurate system (2.3) where the fractional orders are taken as bifurcation parameter. We study the stability criteria of the fractional system (2.3) by using Theorem 2 in [7]. It shows that the stability depends on the sign of . Here we see that the number of eigenvalues, i depends on the term which depends on the different fractional orders (generally, i is greater than the number of equations of the given system). Here we consider the orders of fractional derivative as the bifurcation parameters.

Now we define a functionFrom Theorem 2 in [7], we can see that EE point will be stable asymptotically if and it will be unstable if .

To examine the existence of the Hopf bifurcation for incommensurate fractional system, we state the pursuing theorem.

Theorem 4.2

[7] If the bifurcation parameters meets its critical value , then the system (2.3) passes through Hopf bifurcation at the EE point , if the given singularity criteria [(i)] and also the transversality criteria [(ii)] are satisfied.

, where is the directional derivative of the vector curve at the critical value of bifurcation parameter .

We calculate , the critical value of the parameter for fixed . Let us choose fixed and . Then the equation determine as (0.692, 0.94, 0.96). At the critical value , Fig. 11 establishes that the Hopf bifurcation occurs. Also if increases beyond the critical value , the limit cycles will be an attractor with larger radius.

Fig. 11

Limit cycles for at and

Numerical application for some real world data

In this section, we apply the above system (2.1) to predict the COVID-19 situation of India after the second wave. For this, we collect the data of daily (cumulative) active infected of India from 1st March 2021 to 15th August 2021 from the official websites of the Government of India [6] and taking the parameter set as and At this environment, to fit the above data with the above model, we use the MATLAB-16 software to drow the prediction graph up to February 2022 which is shown in Fig. 12 and we see that the number of daily active infected will be decreased and at 28th February 2022, the number of active infected will be approximately 19100. From this graph, we can predict the daily active infected of any day between this period of time.

Fig. 12

Future prediction of daily active infected of COVID-19 cases in India

Conclusions

In this research work, we have constructed and studied a new SIR epidemic model with the help of fractional differential equations where the disease transmission and disease treatment control get saturated after a certain time. Both commensurate and incommensurate fractional differential equations are used to investigate the complete nature and dynamics of infectious diseases. The uniqueness, positiveness, and uniform boundedness of the solutions of the system (2.4) have been established. The threshold parameter to discuss the infectious disease dynamics is commonly known as basic reproduction number is derived, and it has also been confirmed that the DFE is locally and also globally stable when the threshold condition holds. On the other hand, the system also passes through a backward bifurcation at the EE point for the threshold range Moreover, several computer simulations works enrich our theoretical results. A stability region for the corresponding fractional system is identified. Along with the Hopf bifurcation phenomenon considering the fractional derivative i.e., as bifurcation parameter is visualized graphically. We have further observed that with the increasing value of the length of the cycle increases, and this phenomenon established that by increasing the importance of , the system repels away from the equilibrium point. Moreover, the Hopf bifurcation scenario for the incommensurate fractional system is also derived with bifurcation parameters , , and . We also experienced the significance of with the help of the model fitting.

Also, we have drow a prediction curve of daily active infected of India from 1st March 2021 to 28 February 2022, i.e., of one year with help of the existing data from 1st March 2021 to 15th August 2021. From this prediction curve, we can conclude that if the present situation (i.e., the value of the parameter that we have assumed here) does not change, there will be no chance of a third wave before February 2022.