Chaotic dynamics in a novel COVID-19 pandemic model described by commensurate and incommensurate fractional-order derivatives.

Mathematical models based on fractional-order differential equations have recently gained interesting insights into epidemiological phenomena, by virtue of their memory effect and nonlocal nature. This paper investigates the nonlinear dynamic behavior of a novel COVID-19 pandemic model described by commensurate and incommensurate fractional-order derivatives. The model is based on the Caputo operator and takes into account the daily new cases, the daily additional severe cases, and the daily deaths. By analyzing the stability of the equilibrium points and by continuously varying the values of the fractional order, the paper shows that the conceived COVID-19 pandemic model exhibits chaotic behaviors. The system dynamics are investigated via bifurcation diagrams, Lyapunov exponents, time series, and phase portraits. A comparison between integer-order and fractional-order COVID-19 pandemic models highlights that the latter is more accurate in predicting the daily new cases. Simulation results, besides to confirming that the novel fractional model well fit the real pandemic data, also indicate that the numbers of new cases, severe cases, and deaths undertake chaotic behaviors without any useful attempt to control the disease. Supplementary Information: The online version contains supplementary material available at 10.1007/s11071-021-06867-5.

Introduction

From Wuhan in China to the whole world, COVID-19 epidemic was first broke out in December 2019, where its biological hazards as regards of lethality and propagation were not familiar [1]. The so-called coronavirus SARS-CoV-2 is liable for such disease, which mostly causes an inflammatory storm and a viral pneumonia with certain symptomatology [2-4]. In the interest of assessing the epidemic tendency, controlling the COVID-19 epidemic and minimizing its harms, it is necessary to have knowledge about the behavior of this disease in terms of the speed of virus infections, the duration of its symptoms prior to diagnosis, the time of its peak and the time of its low points. One of the main aspects that can definitely help decision-makers to face such crisis is proposing suitable mathematical models, which in their turns, can offer easy, fast, and effective access to the optimum forecasts and predictions [5].

In epidemiology field, a lot of schemes have been developed to model several infectious epidemics mathematically. The compartment models, which divide communities into certain major classes, are the most employed models. The interactions between those classes are mainly determined by certain beforehand mathematical formulas. In general, a complete model of an epidemic disease with its equations cannot be easily formulated within these formalisms due to the novelty of such disease and rapidly developing in its shape and behavior [2]. For instance, the so-called Susceptible–Exposed–Infectious–Removed model, or simply SEIR model, still cannot precisely estimate the broad range of infections resulting from COVID-19 [5, 6]. The reader may refer to the references [7-11] for gaining further details about such common model. In spite of all these research studies and many others, there are further complex models which have been, not long ago, proposed to describe the lightning-fast spread of several epidemic diseases around the world. More recently, Mangiarotti et al. have proposed in [2] a novel deterministic mathematical model of COVID-19 pandemic in light of the two official data sets given from two official institutions; the National Health Commission of the People’s Republic of China [12] and the Johns Hudson University [13]. These data have been recently collected from these two sources by taking into consideration the spread of this epidemic disease in each of China, South Korea, Japan and Italy, for the period from 21 January till 10-April 2020. The complete form of this model has been described by a nonlinear three-dimensional system consisting of the following states [2]:where C represents the number of daily new cases of COVID-19, S represents the number of daily additional severe cases (positive or negative), D is the number of daily deaths, and where , (), represent the system’s parameters which are mentioned below in Table 1.

Table 1

Parameter values of system (1)

Parameter	Value	Parameter	Value
α_1	-0.10530723	α_6	0.44040714
α_2	2.343×10^-5	α_7	0.16060376
α_3	0.15204	α_8	-0.00011493
α_4	-0.01451520	α_9	-1.215×10^-5
α_5	-0.20517824	α_10	0.2844499
α_11	2.38×10^-6

In fact, Mangiarotti et al. have employed the time-series plot together with the phase portraits to prove that their proposed COVID-19 model exhibiting chaos. In particular, through comparing their obtained results with the observed real data, they have found the existence of chaos in their model (see [2, 14]). Besides the integer-order model described in [2], some dynamical systems described by fractional-order derivatives have been recently introduced to model infectious epidemics, including the COVID-19 pandemic ([15-26]). For example, in [17] the analytical and computational aspects of a fractional-order COVID-19 model have been studied. The model consists of five compartments and includes a Mittag–Leffler kernel. In [18] a fractional-order COVID-19 model for examining the consequence of adaptive immune responses to the viral mutation is proposed. Specifically, in [18] three populations have been considered, i.e., the uninfected epithelial cells, the infected cells, and the SARS-CoV-2 virus. In [19] the COVID-19 pandemic is modeled via the fractional-order SIDARTHE mathematical model. In particular, in [19] the existence of a stable solution is proved and some fractional-order conditions for generating a control strategy are given. In [20] the dynamical behavior of a fractional COVID-19 model is investigated when applied to study the spread of the disease in some Brazilian cities. In [21] variable memory indexes are introduced in a SIRD epidemic model to estimate the COVID-19 pandemic. In particular, in [21] incommensurate fractional-order derivatives defined by a time-dependent function are taken into account when studying the dynamics of the model. In [22] a fractional-order SEIRD model for the spread of COVID-19 is presented. By using the real data of Italy (reported by the World Health Organization), the results in [22] show that the considered fractional model provides a better prediction than the corresponding integer model. In [23] an epidemic COVID-19 model based on the fractional Caputo–Fabrizio derivative is illustrated. The existence and uniqueness of the system solution is investigated by using the Picard–Lindelöf theorem. In [24] a generalized fractional-order SEIR model is proposed. Specifically, the paper shows that, according to the real data of the USA, the considered fractional model has a good prediction ability for the epidemic trend in the next 2 weeks. In [25] numerical simulations of fractional-order modeling of COVID-19 in the case of Wuhan (China) have been carried out. In particular, the Adams–Bashforth numerical scheme has been used in the simulations of the Caputo–Fabrizio fractional-order derivative. In [26] [N] a fractional-order SEIHDR model for COVID-19 with inter-city networked coupling effects is presented. Based on these considerations, in this work, we intend to operate the Caputo fractional-order operator on system (1) and then study and explore the resultant dynamic behaviors of its three states through using some analytical and numerical useful tools.

In this work, we intend to operate the Caputo fractional-order operator on system (1) and then study and explore the resultant dynamic behaviors of its three states through using some analytical and numerical useful tools. The proposed nonlinear fractional-order COVID-19 model will be considered here in view of the commensurate and incommensurate fractional-order cases. The stability of the equilibrium points of the proposed model will be completely analyzed by continuous varying the fractional-order derivative value, and its dynamical behaviors will be then compared with each other. For further details about the Caputo differential operator and how could be operated on a certain mathematical model consisting of nonlinear differential equations, the reader may refer to the references [27-29]. However, the remaining of this work is arranged in the following manner: The next section introduces some needed definitions and preliminary results that pave the way to understand the rest of the paper. In Sect. 3, the commensurate/incommensurate fractional-order COVID-19 model is established. In Sect. 4, the dynamic analysis of such proposed model is studied in view of its phase portraits, bifurcations, Lyapunov exponents and also in view of the stability of its equilibrium points. Section 5 introduces some discussions regarding to the time-series plot of the proposed fractional-order COVID-19 model. Finally, the conclusion of this work is set out in the last section.

Preliminaries

In this section, we intend to state certain key definitions in regard to the non-integer calculus reported in [30]. Such definitions will lay the foundation to fractionalize the integer-order COVID-19 model established in [2] according to the Caputo differential operator.

Definition 1

The integral operator of fractional-order q in the sense of Riemann–Liouville of the function is outlined as:where , , , and indicates the gamma function defined via the following convergent improper integral for any complex number z with strictly positive real part:

Definition 2

The differential operator of fractional-order q in the sense of Caputo of the function is outlined as:where , and .

We would point out that, even though there are several definitions of fractional derivatives in literature, we have selected the fractional Caputo derivative because it offers the advantage that its initial conditions take the same form as that for integer-order differential equations [30]. Namely, while other fractional operators need initial conditions containing limit values whose physical meanings are not very clear, the Caputo derivative requires standard initial conditions in terms of derivatives of integer order, which have clear physical interpretation [30]. However, in order to go forward in our presented work, we recall the following significant result that addresses the stability of the equilibrium points for commensurate fractional-order systems reported in [31].

Theorem 1

If the eigenvalues of the Jacobian matrix J at the equilibrium point E satisfy the following condition:then the following commensurate fractional-order system:will be asymptotically stable, where , and .

On the other hand, for addressing the stability of incommensurate fractional-order systems, a necessary condition has to be set up in a similar manner of commensurate order case. In order to attain this objective, let us consider the following incommensurate fractional-order dynamical system [31]:where , and , are positive integers such that and for . Define M as a least common multiple of ’s. Let be an equilibrium point of system (4) and define as a small perturbation from a fixed point for . Then, we have the following consecutive states:This, consequently, implies:One can easily observe that system (5) is equivalent to the following matrix form:where J is the Jacobian matrix evaluated at a given point p, and has the form:DefineThen the solution of the linear system (6) is asymptotically stable if all the roots of the equation satisfy the condition [32]. In particular, this condition is equivalent to the following inequality:Thus; an equilibrium point p of system (4) is asymptotically stable if condition (8) is satisfied. Actually, the term is typically called the instability measure for equilibrium points in fractional-order systems, or simply IMFOSs. Hence, the necessary condition for exhibiting a chaotic attractor for the fractional-order system (4) can be summarized as follows [31]:

Mathematical model

In light of the above preliminaries together with the COVID-19 model (1) proposed in [2], the following nonlinear fractional-order version is established:where C, S and D are the state-variables of system (1), and is the Caputo differential operator of order for which and . Observe that system (10) will be called commensurate fractional-order system if , otherwise it is called incommensurate one. Actually, the predictor–corrector method (PCM), which was introduced by Diethelm as an improved version of the Adams–Bashforth–Moulton algorithm (ABMA), can be employed to provide a numerical solution to a nonlinear system consisting of a number of fractional-order differential equations formulated via Caputo operator (see [33]). Furthermore, the PCM can be employed to find the Jacobian matrix of system (10), which would be here as follows:In view of the above result and the parameter values given in Table 1, one can obtain the eigenvalues of J at two equilibria of the system; the first one is the origin which yields the eigenvalues , whereas the second one is which yields, in its turn, the eigenvalues . More particularly, these two equilibrium points and can be classified according to the previous discussion as an unstable node and a saddle-focus node, respectively.

Dynamic analysis of model

In this part, we intend to examine the existence of a chaotic behavior of the commensurate and incommensurate fractional-order COVID-19 model which is previously proposed. Such examination will be carried out using some numerical tools such as constructing bifurcation diagrams, computing Lyapunov exponents and sketching phase portraits in 2D and 3D projections.

The commensurate fractional-order model

The investigation of the stability of equilibrium points for the fractional-order COVID-19 model is considered extremely necessary for better understanding and exploring different complicated behaviors of its dynamics. From this point of view, some numerical analysis are performed below in order to determine the behavior of the system’s trajectories through continuous varying the fractional-order derivative value. As we mentioned a little while ago, in view of selecting the system’s parameters , (), as previously given in Table 1, the eigenvalues , , of the Jacobian matrix J at the equilibrium point will be as . Thus, if one substitutes the two previous eigenvalues and into the condition given in Theorem 1, the following result will be gained:which immediately implies that system (10) will exhibit a chaotic behavior when . Such considerable result can be numerically confirmed by carrying out some numerical simulations. In particular, the bifurcation diagram of system (10) can be constructed to be as shown in Fig. 1a for according to the initial condition . From such diagram, we can see that system (10) exhibiting asymptotic stability when , whereas it starts losing its stability and begins its behavior from a periodic motion to chaos mode when . On the other hand, in light of the fact that a presence of any positive Lyapunov exponent of the fractional-order system shows a chaotic behavior for it dynamics, we observe that system (10) exhibits a chaotic behavior for as shown in Fig. 1b. To further clarification, the phase portraits in CS-plan is sketched and exhibited in Fig. 2 according to the same system’s parameters given in Table 1 and to the same initial condition . It can be seen from such figure that system (10) will be asymptotically stable according to the equilibrium point , when (see Fig. 2a). From the other side, Fig. 2b-2c shows that the system loses its stability and then begins to construct a scroll around the point when and . In addition, a chaotic attractor is exhibited for system (10) when (see Fig. 2d), whereas another complex chaotic attractor is appeared when . For completeness, Fig. (3) represents a 3D projection sketch of the chaotic attractor of system (10), while Fig. (4) represents different 2D projection sketches for the complex chaotic attractors generated from the dynamics of the system, confirming hence the previous reported results.

Fig. 1

a Bifurcation diagram of system (10) for , b Lyapunov exponents of commensurate system (10) for according to the system’s parameters given in Table 1 and according to the initial condition

Fig. 2

Phase portraits of system (10) in CS-plan according to the data given in Table 1 and according to the initial condition for: a , b , c , d and e

Chaotic attractor of system (10) in 3D-projection for according to the data given in Table 1 and according to the initial condition

Chaotic attractor of system (10) in different 2D projections for according to Table 1 and to the initial condition for: a CS-plan, b SD-plan and c CD-plan

The incommensurate fractional-order model

In a similar manner to the previous subsection, we intend here to study the existence of a chaotic behavior of the incommensurate fractional-order COVID-19 model by applying the same numerical tools used before such as constructing the bifurcation diagram and sketching the phase portraits of the trajectories of system (10) in 2D and 3D projections. More particulary, we explore the behavior of the dynamics of system (10) under continuous varying of the incommensurate fractional-order value. The system’s parameters are selected here to be as given in Table 1 and the initial condition is also assumed as . To study the stability of system (10) in its incommensurate order case, two bifurcation diagrams are plotted in Fig. 5a and Fig. 7a according to two corresponding cases; the first one takes the incommensurate fractional-order values as and , while the second one takes them as and . It can be seen from such figures that the equilibrium point is clearly asymptotically stable when and . On contrary, system (10) starts losing its stability and begins its behavior from a periodic motion to chaos mode when and . Furthermore, based on the plot of Lyapunov exponents shown in Fig. 5b and Fig. 7b, we observe that system (10) is chaotic when and . In connection with of the phase portraits of system (10) in its incommensurate orders, Fig. 6 and Fig. 8 represent the sketches of them on CS-plan according to different values of and . In particular, when the incommensurate fractional order , the incommensurate fractional-order system (10) will be asymptotically stable to the equilibrium point (see Fig. 6a. Such system will be asymptotically stable to a limit cycle when (see Fig. 6b), while a chaotic attractor will be appeared when (see Fig. 6c), and moreover the system under consideration will exhibit a complex chaotic attractor when (see Fig. 6d). In a similar manner to the previous discussion, system (10) in its incommensurate orders will be asymptotically stable to the equilibrium point when (see Fig. 8a). The same system will be asymptotically stable to a limit cycle when (Fig. 8b), whereas a chaotic attractor will be appeared when (see Fig. 8c) and it will be also exhibit a complex chaotic attractor when (see Fig. 8d). In addition to these simulations which are performed according to the data given in Table 1 and to the initial condition , a 3D projection sketch of the complex chaotic attractor of the incommensurate order system (10) is furthermore exhibited in Fig. 9 by considering two cases; the first one is performed when and , while the second one is performed when and .

Fig. 5

a Bifurcation diagram of incommensurate system (10) for and b Lyapunov exponents of incommensurate system (10) for , and according to the data given in Table 1 and to the initial condition

Fig. 7

a Bifurcation diagram of incommensurate system (10) for and b Lyapunov exponents of incommensurate system (10) for , and according to Table 1 and to the initial conditions

Fig. 6

Phase portraits of incommensurate system (10) according to Table 1 and to the initial conditions for: a , b , c and d

Fig. 8

Phase portraits of incommensurate system (10) according to Table 1 and to the initial conditions for: a , b , c and d

Fig. 9

Chaotic attractor of incommensurate system (10) in 3D projections according to Table 1 and to the initial condition for: a and and b and

For more clarification, we notice that some Lyapunov exponents are positive, particularly for the two cases: The first one occurs when , (see Fig. 5b), while the second one occurs when , (see Fig. 7b). From this standpoint and in view of the stability condition (9) reported in Sect. 2, we might take into account these two cases as follows:As a result of all previous numerical findings, we can conclude that the fractional-order version of the COVID-19 pandemic model exhibits chaotic behaviors in accordance with the two considered cases; the commensurate and the incommensurate fractional-order cases. It can be furthermore concluded that the chaotic ranges generated by system (10) in its commensurate and incommensurate orders are changed according to a change of these orders, and the minimum fractional-order value that can exhibit a chaos for the system under consideration is which is occurred when .

When and , then the least common multiple M is become as . This, consequently, implies , which yields . Hence, the IMFOS of system (10) would be as: This means that such system shows a chaotic behavior, which confirms the resultant numerical plot shown in Fig. 9a.

When and , we can obtain , which leads to deduce to be as . This, however, yields , and consequently implies the IMFOS of system (10), which would be in the form: Similarly, we notice that the system shows also a chaotic behavior, which confirms the numerical plot result shown in Fig. 9b.

Discussion

For further discussion and better understanding of the chaotic oscillations of the COVID-19 pandemic, we intend to use the time-series plot in our analysis. Here, the time used for performing some numerical comparisons is chosen from the day 21 till the day 90 because model (1), which was established for the outbreak of COVID-19 in China, was considered for the period that began on January 21, 2020 and ended on April 10, 2020 (see [2, 14]. The time series plots of the fractional-order COVID-19 system (10) are illustrated in Fig. 10 according to the following three cases: Daily new cases C (Fig. 10a), daily additional severe cases S (Fig. 10b), and daily deaths D (Fig. 10c). Within these numerical comparisons, the blue, green and red lines represent the time-series plots sketched when , and , respectively. In view of the above numerical results, we can continuing resume exploring the recent results obtained from the integer-order system and the observed real results given in Table 2 [2]. It can be observed that the maximum number of the daily new cases obtained from the fractional-order COVID-19 system is the same maximum number obtained from the real data which is . Such observation confirms definitely that the use of the fractional-order system in the expectation for new cases is better than the use of the integer-order system, which it previously expected that these cases will reach to . On the other hand, we observe that the two integer- and fractional-order systems expect approximately the same minimum and maximum numbers of the daily severe cases and the daily death cases, which are both close to the real data.

Fig. 10

Time series plot of the three considered cases of the fractional-order system (10) according to Table 1 and to the initial condition for: (blue lines), (green lines) and (red lines)

Table 2

The minimum and maximum numbers of different cases

Cases	Real data	Integer system	Fractional system
	(min,max)	(min,max)	(min,max)
Daily new cases C	(186,5000)	(186,8000)	(186,5000)
Daily additional severe cases S	(-1000,1500)	(-50,1200)	(-2,800)
Daily deaths (8)	(8,140)	(8,100)	(8,80)

In the same vein and due to the fact that confirms the fractional-order systems can include a memory effect unlike the integer-order systems, we try next to predict the dynamic of the fractional-order COVID-19 model through taking a wider range of time than that taken in [2]. In particular, by taking, e.g., the commensurate fractional-order values as , the time-series plot of the three state variables C, S, and D are calculated and shown in Fig. 11, for and . It can be seen from such figure that the chaotic oscillations still exist for the three considered cases. This means that the number of the new cases, severe cases and also death cases will continue taking a chaotic dynamical behavior and will not be decreased with the passage of time without any serious attempt to control the COVID-19 pandemic.

Fig. 11

Time series plot of fractional-order system (10) according to Table 1 and to the initial condition , when for: a , b , with noting that the blue line represents the daily new cases, the red lines represents the daily severe cases and the black lines represents the daily deaths

Conclusion

This paper has presented a novel COVID-19 pandemic model described by commensurate and incommensurate fractional-order derivatives. By analyzing the stability of the equilibrium points and by varying the values of the fractional-order derivative, the paper has shown that the proposed model exhibits chaotic behaviors. The system dynamics have been investigated via bifurcation diagrams, Lyapunov exponents, time series and phase portraits. A comparison between integer-order and fractional-order COVID-19 pandemic models has clearly shown that the latter is more accurate in predicting the daily new cases. Simulation results have also indicated that the numbers of new cases, severe cases and deaths undertake chaotic behaviors without any useful attempt to control the disease. This would help the decision makers to better understand the epidemiological behavior of the COVID-19 disease over time. Consequently, this would also help them to select measures for effectively monitoring and controlling such pandemic.

Below is the link to the electronic supplementary material.

Supplementary material 1 (rar 1947 KB)