key: cord-0889467-zojglnbc
authors: Jahanshahi, Hadi; Munoz-Pacheco, Jesus M.; Bekiros, Stelios; Alotaibi, Naif D.
title: A fractional-order SIRD model with time-dependent memory indexes for encompassing the multi-fractional characteristics of the COVID-19
date: 2021-01-10
journal: Chaos Solitons Fractals
DOI: 10.1016/j.chaos.2020.110632
sha: a52b33874151c0aa5fb993d1c44524df2635f9ec
doc_id: 889467
cord_uid: zojglnbc

COVID-19 is a novel coronavirus affecting all the world since December last year. Up to date, the spread of the outbreak continues to complicate our lives, and therefore, several research efforts from many scientific areas are proposed. Among them, mathematical models are an excellent way to understand and predict the epidemic outbreaks evolution to some extent. Due to the COVID-19 may be modeled as a non-Markovian process that follows power-law scaling features, we present a fractional-order SIRD (Susceptible-Infected-Recovered-Dead) model based on the Caputo derivative for incorporating the memory effects (long and short) in the outbreak progress. Additionally, we analyze the experimental time series of 23 countries using fractal formalism. Like previous works, we identify that the COVID-19 evolution shows various power-law exponents (no just a single one) and share some universality among geographical regions. Hence, we incorporate numerous memory indexes in the proposed model, i.e., distinct fractional-orders defined by a time-dependent function that permits us to set specific memory contributions during the evolution. This allows controlling the memory effects of more early states, e.g., before and after a quarantine decree, which could be less relevant than the contribution of more recent ones on the current state of the SIRD system. We also prove our model with Italy’s real data from the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University.

COVID-19 was stated as a pandemic outbreak on January 22, 2020, by the World Health Organization. This new coronavirus affects the world population in several ways. On the one side, it is evident the health issues for older people and with comorbidities, mostly. Indeed, the COVID-19 is also breaking the economic systems, limiting human interaction, and above at all, testing our mental toughness. Therefore, researchers from several scientific areas have devoted to studying COVID-19 [1] [2] [3] [4] [5] . The purpose is to contribute to the understanding, analysis, identification, prediction, and so forth from many points of view. Any new idea is a step forward to cope with this health issue. This paper goes in the same direction but considering fractional-order dynamical systems and variable memory-indexes. To this end, we first found that a SIRD model is one of the main mathematical models for studying the COVID-19 [19] [20] [21] [22] [23] [24] [25] [26] [27] . Based on the fractional-order derivative with singular kernel we found that Rezapour et al. analyzed the transmission of COVID-19 in Iran and in the world using a SEIR model with constant fractional-order [22] ; Amhed et al. proposed a five-term dynamical system to understand the trade-off between the lockdown and the spread of the virus [23] ; Tuan et al. determined a reproduction number of COVID-19 from a fractional-order mathematical model; Shah et al. studied a constant fractional-order model with the Haar collocation method as numerical solver to reveal the transmission dynamics [24] ; Rajagopal et al. proposed a fractional-order SEIRD model for the spread of COVID-19 and compared it with the real data and integer-order cases [25] ; Abadias et al. analyzed the transmission of perturbations across the amino acids of a protein represented as an interaction network applying the Caputo derivative in a SI model [26] ; and Lu et al. presented a fractional SEI-HDR model based on the coupling effect of inter-city networks to cope with the mortality rates (exposure, infection and hospitalization) and the infectivity of individuals during the incubation period [27] .

Nevertheless, the mentioned works only considered constant fractional-orders q, i.e., the same memory indexes for the entire COVID-19 evolution. The heterogeneity of the countries population and the specific countermeasures imposed by the governments, such as mobility restriction, lockdown, experience or knowledge of individuals about that disease, quarantine, home office, homeschooling, and so forth, brings up the idea that memory indexes are not constant during the epidemic process. The remaining question is, what does implication has the time-dependent fractionalorders q (t) (variable memory indexes) in the prediction of COVID- 19 ? Herein, we give a positive answer.

Additionally, previous works showed that COVID-19 evolution follows scaling features similar to power-law [28] [29] [30] [31] [32] [33] . The fractalness, also know as scale-invariant and self-similar, is a universal property of nature that permits the interpretation of the growth of trees, snowflakes, etc. Fractals are formed by an identical pattern over different scales obtained from a simple recursive rule. In the real-world, the fractals have been observed both in the spatial and temporal domains. In the case of temporal signals, the analysis tools of scale invariance determine whether a frequency represents a transcendent role in a process. The scale-free signals present a power-law relationship with a spectral power at frequency f,

where β is related to Hurst exponent [34] . For the COVID-19 case, Altan and Karasu's seminal work focuses on diagnosing people infected with COVID-19 at an early stage using chest X-ray images and a hybrid model [35] . This model consists of 2D curvelet transform, chaotic salp swarm algorithm, and EfficientNet-B0 method. In this manner, the chest X-ray images were classified as normal, COVID-19 positive, and viral pneumonia with high accuracy compared to previous works. Because the approach in Altan and Karasu [35] is based on 2D curvelet transform, which is the best fit for the signals' multiscale behaviors, the observed results would imply a fractal nature for the COVID-19 epidemic. Additionally, Ref. [36] reported the spectral, nonlinear time series, fractal, and complexity analysis of vesicular (VB) and bronchial (BB) breath signals for COVID-19 auscultation. Hurst exponent of the signals revealed fractal signatures with two distinct powerlaw exponents between VB and BB signals, suggesting again a multifractality.

Materassi suggested some ideas about a geometric scenario where the generalized logistic equations could describe the COVID-19 infection [28] . Abbassi et al. found a power-law and scaling behaviors for the number of contamination cases as a function of the city population [31] . Gowrisankar et al. studied top 15 affected countries and detect the infection rates behave as a power-law growth. Also, they presented the generalized fractal dimensions of confirmed case and confirmed death of COVID-19, resulting with a different value for each country [33] . Ziff et al. argued that the current growth closely follows power-law kinetics, underlying fractal behaviors [29] . Additionally, they discovered the power-law exponents changes with the time, presumably by the containment efforts.

In fact, we analyzed the experimental time-series (from February 22 to August 26) of twenty-three countries using fractal formalism, and found that the time-series show the same scaling pattern and present a universality among countries to some extent. Besides, our results suggest that the time-series exhibit a spectrum of power-law exponents (no just a single one), which agrees with Refs. [29, 32, 33] where was discussed that the COVID-19 infection rate may cluster in three stages. In the former (growthstage), the infections obey a power-law behavior characterized by a scaling exponent at a certain speed. In the second stage (spread of the virus stage), the initial power-law range splits into a second one with a smaller scaling exponent. Then, the entire behavior presents distinct power-law exponents. Finally, the infections stop in the last-stage. It means that the COVID 19 could follow a memoryless evolution or with low memory (fractional-order very close to unity) at early stages. As time evolves, the memory implications will be more evident, becoming a pandemic process with diverse memory indexes.

Hence, we incorporate variable memory indexes in the proposed incommensurate fractional-order SIRD model. The memory contributions are updated by considering the fractional-orders as a time-dependent function q (t) . We choose the Caputo fractional derivative because it has a power-law kernel where its decaying rate depends directly on the fractional-orders. As the fractionalorder varies, the memory contributions, i.e., the power-law scaling features, control and emulate an evolution with multiple exponents. In this manner, we can switch between memory indexes with distinct strengths, i.e., long-and short-memory.

The main contributions of the proposed approach are 1. First work introducing variable memory indexes in a SIRD epidemic model. Therefore, the scenarios that alter the infection speed, such as containment measures, indirect protection, virus varieties, vaccination, and so forth, may be considered by our model to predict the COVID-19 propagation with more accuracy. It means that the rate of the infected, recovered, and dead populations not only can be established independently but also be adjusted during the COVID-19 progress. Those two characteristics are feasible since the memory strength is controlled by time-dependent fractional-orders defined as a piecewise linear constant function. 2. By adding incommensurate fractional-orders as a timedependent function plus a fractional derivative with a powerlaw kernel, the proposed fractional-order SIRD model considers various power-law exponents to estimate the COVID-19. This indicates that the diverse outbreak stages (growth, spread, and decrease), which are assumed to have multiple power-law exponents (multiscaling properties), can be analyzed. As a result, we present a model to elucidate the multifractal behavior of COVID-19. 3. Identification of universality in the COVID-19 progression among geographical regions by using the generalized Hurst exponents method. The results suggest that the governments' adopted measures to control the epidemic are similar among countries in the same areas. However, not all countries share physical borders. It may imply that cultural, social, and economic aspects are also crucial to understanding how the pandemic outbreak affects the populations.

The paper consists of five sections. Section 2 presents the experimental time-series analysis for 23 countries from different world regions using generalized Hurst exponents. Section 3 introduces the proposed fractional-order SIRD model with multiple memory indexes. It is also discussed its stability and nonnegative solutions along with a dedicated numerical solver. Section 4 presents a case of a study to compare the proposed model's predictions concerning Italy's real data. Section 5 concludes the work.

This section presents the analysis of the time-series sets for 23 countries. We consider the time-series of the daily confirmed cases (DCC) from January 22th (the date when WHO declared COVID19 as a pandemic outbreak) to August 26th. The data were retrieved from COVID-19 Dashboard by the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University (JHU) [37] .

It is well known that fractalness is a ubiquitous property of nature. This scale invariance appears not only in the spatial-domain but also in the temporal-domain. In this regard, several complex physical, human, and biological phenomena have been analyzed with this perspective [34, [38] [39] [40] [41] [42] . Long-memory dependence in time series or spatial data are associated with power-law correlations, which can be characterized by the Hurst ( H) coefficient [43] [44] [45] . When 0 ≤ H ≤ 0 . 5 , the time-series behaves very complex with low correlation ("short-memory"). For 0 . 5 < H ≤ 1 the timeseries has a long-memory dependence. Besides, to identify the multi-scaling properties of time series, the α-order Hurst exponent H α can be used instead [36, 46] .

A stochastic process X (t) is called multi-scaling if it has stationary increments and satisfies E( | X ( t ) | α ) = c (α) t τ (α)+1 for all t ∈ F, α ∈ L , with F and L intervals on the real line, τ (α) and c(α) functions with domain L , and E( . . . ) the expectation value [46, 47] . The function τ (α) is called the scaling function of the multi-scaling process. To analyze the DCC time-series, we have chosen the generalized Hurst exponent method [46] . This approach is based on the αth-order moments to detect the multiscaling complexity of the time-series, which avoid the sensitivity to any type of dependence in the data [4 8,4 9 ] .

Let us define the statistical evolution of time-series X (t) as

where the time interval τ can vary between ϑ and τ max . Then, the generalized Hurst exponent H(α) can be defined by

When H(α) is constant and independent of α, a uniscaling is described. It means that there is a linear behavior between αth-order moments, e.g., H(1) = H (2) . On the other hand, if H(α) is not constant, a multiscaling process is observed, and different exponents characterize the scaling of the α-moments. Additionally, the classical Hurst exponent arises when α = 1 , H(1) , whereas the generalized exponent at α = 2 , H(2) , is associated with the scaling of the autocorrelation function. Finally, the time-series may exhibit a power-law relationship, where the spectral density,

By applying Eq. (1) , we compute the generalized Hurst exponent for 23 countries, eight, seven, and eight from America, Europe, and Asia, respectively. For a sample of five countries, Fig. 1 (a) -(e), shows the (1st, 2nd, 3rd)-order moments in a loglog plot relating τ and K α (τ ) . The time interval (days) τ varies as follows. Turkey (ϑ, 13) , and Germany (ϑ, τ max ) = (1 , 11) . As can be noted, there is a scaling factor for K α (τ ) as τ varies, demonstrating the Eq. (1) . Besides, the timeseries sets reveal multiscaling behaviors since the results given in Fig. 1 (f) manifest a nonlinear trend for αH(α) with α ∈ (0 , 3) . We observe that the nonlinearity is more evident for certain countries than others, e.g., India, US, Saudi Arabia, and Turkey, contrary to Iran, Italy, UK and Germany. It is worthy to remark that we are not trying to predict the persistent or antipersistent evolution of COVID-19, but instead, we are arguing the hypothetical power-law behaviors, which have been reported in previous works [28] [29] [30] [31] [32] [33] . In this manner, our results suggest that the time-series exhibit a spectrum of power-law exponents (no just a single one), which agrees with Refs. [29, 32, 33] where was discussed that the COVID-19 infection rate evolve in three stages. Each one with different power-law exponent indicating that the power-law exponents changes with the time, presumably by the containment efforts.

Next, we focus on the H(1) and H(2) moments of Eq. (1) . Fig. 2 gives the estimated of H(1) and H(2) , respectively, for 23 countries strongly affected by the outbreak. We also compare the results with those obtained from Detrenead Function Analysis (DFA) [49] . Since the results are similar, they were omitted for plotting. Fig. 2 (b) depicts the generalized Hurst exponent H(2) at α = 2 which is related to the scaling of the autocorrelation function and to the power spectrum. We found that according to the geographical area, the countries present similar values for the generalized Hurts exponents H (2) , and therefore, we can classify them, roughly speaking, as follows: North America and Latin American countries (Chile, Ecuador, Argentina, Peru, Brazil, Colombia, Mexico, excluding US) are clustering at the leftmost edge of the graph due to presents the lowers generalized Hurst exponents H(2) < 0 . 10358 . Most of Asia countries such as Russia, Turkey, Saudi Arabia, Iran, India, excluding Bangladesh, Pakistan, and China, are clustering in the rightmost part of the graph, since they have the highers H(2) > 0 . 46673 . In the middle part of the graph, it can be found grouped Europe countries (Belgium, Netherlands, United Kingdom, Germany, Italy) with 0 . 27433 < H(2) < 0 . 46821 . As an exception, France and Spain are located with lower values.

In particular, we found that Latin American countries (Chile, Ecuador, Argentina, Peru, Brazil, Colombia, Mexico) presents similar H(2) ( Fig. 3 (a) ). It could suggest the adopted measures to face the pandemic, i.e. the government's decisions are analogous among those countries, and also may be associated that they share some cultural, social and economical issues, even although not all the countries share physical borders. Those behaviors repeat for other countries from different continents. For instance, Russia, Turkey, Saudi Arabia, Iran, India also have similar H(2) ( Fig. 3 (b) ), which could mean that they are addressing the pandemic outbreak in a similar fashion. We also detect that certain neighbor countries has likely H(2) in America and Asia, but this not necessarily happens in Europe. This could be associated to European countries faced with particular strategies, however, as previous mentioned, the European countries cluster with a quite similar H(2) ( Fig. 3 (c) ). Due to the results suggest that the outbreak presents different powerlaw exponents as the time evolves, we propose herein a mathematical model which can control the decaying rate, i.e. the memory effects.

One of the most straightforward mathematical models to study the progress of an epidemiological outbreak is the SIRD model [6, 7] . Fanelli and Piazza proposed a SIRD model to foresee the COVID-19 evolution that divide the population in four sub-types: susceptible ( S), infected ( I), recovered ( R ), and dead ( D ) individuals [8] . Its dynamics is governed by the following coupled ordinary 

in which r, a, and d are infection, recovery and death rates, respectively. The infection rate is a system's parameter that varies

The SIRD system (2) represents a Markovian process where the previous states do not affect the current conditions, i.e., a memoryless SIRD system. Nevertheless, the network-like human interactions (including the COVID-19) evolve to a pattern with memory properties were past events impact recent events [19] [20] [21] [22] [23] [24] [25] [26] [27] . To cope with that, we apply the fractional-order calculus as a mechanism to represent the memory of the biological system since the fractional-order derivative operators represent more accurately the real phenomena in nature [17, 18] .

In this manner, we present a fractional-order SIRD model to observe the memory effects in the outbreak. Let us rewrite the the system (2) as a function of time-dependent integrals:

where κ (t − t ) represents a time-dependent kernel which permit us to consider the long-term memory implications. As confirmed in Section 2 and supported by Refs. [28] [29] [30] [31] [32] [33] , the scaling behavior of COVID-19 should be pretty well described by an spectrum of power-law functions as the kernel. Hence, we choose κ (

with 0 < q ≤ 1 and (·) being the gamma function, to replace κ (t − t ) in Eq. (3) . Then, we obtain a fractionalorder SIRD model in the Caputo's sense as follows:

being 0 < q ≤ 1 . When dealing with practical problems, Caputo fractional-order formulations have the same initial conditions, which has a specific physical meaning. Besides, we adopt the Caputo fractional derivative because it has a power-law kernel where its decaying rate depends directly on the fractional-orders. As the fractional-order varies, the memory contributions, i.e., the powerlaw scaling features, control and emulate an evolution with multiple exponents. In this manner, we can switch between memory indexes with distinct strengths, i.e., long-and short-memory. As well known, for dynamical systems, the fractional-order q relates to a "memory index". It means that as the fractional-order q tends to zero the memory effects are more evident because of the chosen time-dependent kernel vanishes more slowly. On the other hand, when q is closer to 1, the implications of the memory are reduced. Finally, if q = 1 the system has no memory. In the previous studies, the memory property is almost all considered as a constant parameter and similar for the whole system [19] [20] [21] [22] [23] [24] [25] [26] [27] . This imply that the fractional-order q is the same for the entire analysis. Hence, the remaining question is how does different memory indexes influences the underlying dynamics of SIRD model to improve the description of the COVID-19 outbreak?

In this framework, we now consider that the fractional-order is not more a constant but rather it can vary continuously as a function of time q (t) . The variable-order operators not only have time-dependent memory but also memory order ( ϕ) when a step-function-type order variation is considered [50, 51] . Three scenarios may arise: the fractional-order q (t, ϕ) q (t) , which means there is time-dependent memory only; the fractional-order q (t, ϕ) q (ϕ) where the memory from past order is weak; and the fractional-order q (t, ϕ) q (t − ϕ) which implies the memory order is strong. Because the policies of countries to face the COVID-19 are heterogeneous, it would be beneficial to have more degrees of freedom in the model to adapt faster to radical changes. Therefore, we consider herein distinct memories indexes in two cases:

i) Incommensurate fractional-orders q (1) , . . . , = q (p) with p = 2 , . . . , d im, and d im being the number of system's dependent variables.

ii) Fractional variable-order q (t, ϕ) q (t) .

Let us to consider the following general short-memory system on the interval [ t 0 , T ] .

. . . with [ t 0 , x 0 ] being the initial condition for t ∈ [ t 0 , t 1 ] and first fractional-order q 0 . As demonstrated in Wu et al. [52] , the principle relies on taken into account only the memory effects from the initial times t i and fractional-orders q i with i = 1 , 2 , . . . for each time interval [ t i , t i +1 ] contrary to consider the initial condition t 0 . Hence, we can define a incommensurate variable-order fractional differential equation as follows

where n, ν, L are positive integers, t ∈ [ t 0 , T ] , L = nν, h = (T − t 0 ) /L, n is the number of segments for variable-order, m is the number of dependent variables, and the incommensurate fractional-order q (k ) j (t) is a piecewise constant function given by

being (k ) the fractional-order for each dependent variable in the selected time interval. By rewritten the SIRD model (4) in the form of Eq. (6) with m = 4 , we obtain the proposed variable-order fractional SIRD model:

with ˆ t = t jν , t ∈ [ t jν , t ( j+1) ν ] , and 0 < q (k ) j ≤ 1 .

Next, we prove that system (8) has non-negative solutions. Let R 4

recalling the Generalized Mean Value Theorem as state in Ref. [53] , we have the following Theorem.

given by (8) 

Proof. By applying Theorem 1, Remark 1 , and recalling Ref. [54] , it is straightforward to demonstrate that

, j = 0 , . . . , n − 1 , since we taken into account only the memory effects from the initial times t jν and fractional-orders q (k ) j contrary to consider the initial condition t 0 . Therefore, the nonnegative orthant is bounded and the vector field points into R 4 + . (r 0 , a, d) is unstable for all q (k ) j ∈ [0 , 1) when the initial population, S 0 , is higher than a threshold defined by the system parameters.

Proof. The characteristic equation of fractional-order model (8) at equilibrium point (S 0 , 0 , R 0 , D 0 ) evaluated in the Jacobian matrix J(E 1 ) is given as found that (r 0 , a, d) are always nonnegative because their biological meaning. Therefore, the coefficient of the reduced polynomial is negative when S 0 > (a + d) /r 0 .

By applying the Descartes rule of sings, we demonstrate that there is only one positive root, and therefore, the criterion | arg λ ii | > (r 0 , a, d) as free parameters, is unstable for

On the other hand, the numerical solution of system (8) is obtained by applying a dedicated predictorcorrector method given by [52] :

In general, virus infection needs close contact between the infected and the susceptible individuals. For the COVID-19, close contact is also required. Then, the people in large cities do not have the same contamination chance since the probability of meeting any other person in the city is very low. Regularly, the population follows a small-group clustering (family, friends, coworkers). In those scenarios, the basic integer-order SIRD model does not hold, as usual, in no-Markovian processes where the memory is not taken into account. However, the memory of past events has substantial implications for studying an outbreak evolution. It means that infected people follow a fractal relationship defined by a nonlinear generalized Hurst exponent as concluded in Section 2 . It was shown that these power-law exponents depend on the geography contagion and presumably also by the sociology of contagion. These conditions appear to be what gives the contagion network a different topology for each country. In this manner, we argued that the COVID-19 outbreak presents fractional kinetics. Additionally, in the early stages of a pandemic, the infection is more likely than exponential growth. Still, those behaviors become power-law relationships as individuals already infected might face increased immunity, and other individuals might have had a mild infection that confers them immunity. These various effects would inhibit the exponential growth of the virus. However, we also argued that containment measures strongly affect the evolution of contagion. Then, social distancing, quarantines, mobility restriction, lockdown, experience or knowledge of individuals about the disease, home office, homeschooling, and so forth brings up the idea that the COVID evolution presents multiple power-law exponents with time. Then, the proposed variable fractional-order SIRD model can capture those scenarios as demonstrated as follows.

We select Italy for demonstrating the incommensurate variableorder fractional SIRD model (8) . All experimental data were retrieved from CSSE at Johns Hopkins University [37] . First, we choose t 0 = 20 days after day one (22/01/2020) because it has been reported that presumably the delay in testing shifted the initial day [8] . Also, we set the initial condition for the susceptible population as S(t 0 ) = 2 . 66 × 10 5 due to they were the cumulative confirmed cases ( C C ) in (18/08/2020), which should be classified in Infected Regarding the infection rate r, it varies with time to follow the suggestion in Fanelli and Piazza [8] . In this manner, r(t ) = r 0 for t ≤ˆ t , and 0 . 66 r 0 e −(t−ˆ t ) / 2 + 0 . 34 r 0 for t > ˆ t , with r 0 = 1 . 2 × 10 −6 . Thereby, the parameters for the variable-order SIRD system (8) t 270 , t 540 , t 810 , t 1080 , t 1350 , t 1620 , t 1890 . To take into account the efforts (quarantine decree) to contain the COVID-19 imposed by Italian government 27 days (08/03/2020) after t 0 , we propose two sets of fractional-orders having as break-point t 270 . As a result, the following incommensurate fractional-orders are obtained,

and,

Figs. 4-6 present the numerical results for the proposed SIRD model (8) by applying the numerical method (10) . Fig. 4 shows the memory effects of the first set of fractional orders q (k ) 0 = (1 , 0 . 8 , 0 . 6 , 0 . 4 , 0 . 2) within t ∈ [ t 0 , t 270 ] on the estimation of the susceptible, infected, recovered, and dead populations in the SIRD system (8) . For each case in Fig. 4 , we varied the fractional-order of one population while fixed an integer-order for the remaining ones. For instance, in Fig. 4 

changes from 1 to 0.2. The time interval is from the initial day of the outbreak in Italy t 0 = (11 / 02 / 2020) to the breakpoint t 270 , which means 27 days (08/03/202) after t 0 when the Italian government imposed a first quarantine decree. As can be noted, the memory implications are low for the susceptible, recovered, and dead populations. This is not the case for the infected individuals where there is a time-delay for reaching the peak of infections as the fractional-order reduces. Those results agree well with the typical evolution of a outbreak at early stages where the number of infected individuals grows almost exponentially. This is expected since we consider only short-memory effects of 27 days after the initial day, and the containment measures has not been declared yet. After t 0 , the evolution is computed with integer-order similar to the standard memoryless SIRD model.

On the other hand, the memory of the second set of fractional orders q (k ) 1 , ... , 7 = (1 , 0 . 8 , 0 . 6 , 0 . 4 , 0 . 2) has strong effects after the break-point t 270 (08/03/202) as shown in Fig. 5 . Analogously, for each case in Fig. 5 , we change the fractional-order of one population while select an integer-order for the remaining ones. The rate of susceptible individuals decays faster as the fractional order minimizes ( Fig. 5 (a) ). For the infected population ( Fig. 5 (b) ), we found that the fractional-order controls the maximum peak along with the speed of contagions. Later, the number of recovered ( Fig. 5 (c) ) and dead ( Fig. 5 (d) ) people increases as a function of the fractional-order. The lower the memory index, the higher the prediction of recovery and death cases. As can be seen, different memory indexes for specific time intervals can change the estimate of the outbreak evolution. Those results are expected because we have incorporated long-term memory implications. The fractional-order was applied almost for the entire time window (162 days), i.e. from t 270 to t 1890 . From the day of the quarantine decree ( t 270 ), the infection rate follows a power-law behavior with a faster decaying rate than initial stage. Since infection rate decays, the recovered population increases whereas the susceptible people reduces. As a result, the proposed SIRD model capture the evolution of the outbreak with different power-law exponents using short-and long-memory indexes.

Next, we compare the Italy's real data and the prediction of the proposed fractional-order SIRD model based on variable fractionalorder for four scenarios: i) integer order, ii) a set of constant fractional-orders, iii) another group of constant fractional-orders, and iv) two sets of fractional-orders as a function of time. Fig. 6 (a) gives the simulation result when the system is considered in the integer-order domain, i.e., the system has no memory at all. We observe that there is no correlation between the experimental data and the prediction of the model, which agrees well with previous works [7, [19] [20] [21] [22] [23] [24] [25] [26] [27] . It indicates that the COVID-19 evolution and control cannot be considered without any memory effect. When the virus spreads within a human population, the experience of individuals and the containment efforts should influence its dynamics. Fig. 6 (b) illustrates the results using only the set of fractionalorder q (k ) 0 , ... , 7 for the susceptible, infected, recovered and death populations in the whole time interval t ∈ [ t 0 , t 1890 ] . As we can see, the approximation of the model to the data is not suitable using only one memory index. Fig. 6 (c) shows the evolution of SIRD model by considering the second set of fractional-orders q (k ) 0 , ... , 7 in t ∈ [ t 0 , t 1890 ] . We note a similar result from previous scenario since the fitting of the model to the experimental data remains not correct although the approximation improves considerably. Notwithstanding that the COVID-19 growth presents a power-law exponent, this may not be constant due to the governments' containment decisions. It was confirmed in Section 2 and supported by Refs. [28] [29] [30] [31] [32] [33] , where the COVID-19 time-series have various powerlaw exponents.The appearance of a fractional exponent and powerlaw behavior suggests an underlying fractal process. However, as stated in Section 2 , the generalized Hurst exponents have a non- (0 . 9 , 0 . 913 , 0 . 907 , 0 . 987) linear behavior inferring that the power-law progression of the COVID-19 in each country changes with time. For Italy, the day of the quarantine decree suggests a modification of the outbreak speed. The implication is that the COVID-19 grows momentarily as a power-law with a specific exponent but then slows down until another exponent. The averaged effect of this evolution yields a power-law behavior with multiple exponents. Therefore, a single one set of fractional-orders, i.e., a memory index with the same strength during all evolution, is not proper to predict the outbreak. Also, it suggests that if the pandemic spreads into new countries, another set of scaling exponents will determine the COVID-19 progression.

Finally, Fig. 6 (d) presents the evolution of the three populations: infected, recovered and deaths of the COVID-19 for Italy from (11/02/20) to (16/08/20) . It is worth noting that the best fit to COVID-19 real data is obtained when two memory indexes with distinct strength are considered. Therefore, time-dependent incommensurate fractional-orders represent more accurately the COVID-19 propagation.

In particular, the first set of fractional-orders (11) is used to compute the dynamics of the SIRD model (8) when t ∈ [ t 0 , t 270 ] . As was demonstrated in Saeedian et al. [7] , Ziff and Ziff [29] , Singer [32] , Gowrisankar et al. [33] at early stages of an epidemic, the dynamics of a SIRD model can be characterized with a short-memory evolution almost similar to a integer-order system. We find a similar result for the evolution of COVID-19 in Italy. For the first days, two fractional-order are about q ≥ 0 . 950 and the others are integer values indicating that the memory has lower implications. Therefore, the decaying rate of the power-law kernel is minimized. Then, the first set of fractional-orders in (11) allowed to represent the growth stage of the COVID in Italy where the total number of infections is in the power-law regime with a single scaling exponent. On the other hand, the second set of fractional-orders (12) is selected for t ∈ [ t 270 , t 1890 ] , the day after the decree of the Italian Government. From this date, the memory has strong effects in the evolution of the outbreak showing another power-law exponent. It is important to note that the second set of fractional-orders q (k ) 1 , ... , 7 is lower than the first set q (k ) 0 which points out that the outbreak is fading. This represents the slow down stage when the Italian government imposed strict control to reduce the spread of the COVID-19, and therefore, the initial power-law range split into a second one with a smaller scaling exponent. The total number is still increasing as a power-law but more slowly than before To compare the proposed approach, we compute the root-mean-square error 

in which rd and ap are the real data and its approximation with the model, respectively, and T is the time of accessible data. Table 1 gives the normalized RMSE for each model. The best fit is obtained for the incommensurate fractional variable-order model. As expected, the RMSE reduces because the proposed model con-siders two distinct power-law exponents to estimate the outbreak for Italy, which is closely related to the multiscaling characteristics of COVID-19 given in Section 2 .

A fractional-order SIRD (Susceptible-Infected-Recovered-Dead) model based on the Caputo derivative for incorporating the memory effects (long and short) as time-dependent functions in the outbreak progress has been introduced. Based on the experimental time series of daily confirmed cases for 23 countries, it was identified that the COVID-19 evolution presents various power-law exponents (no just a single one). Those multifractal signatures are shared among geographical regions. The countries tend to formed groups with a similar generalized Hurst exponent, H(2), suggesting that cultural, economic, political, and sociological aspects can be critical to understanding the outbreak. For all analyzed countries, the H(2) showed a nonlinear behavior proposing that the COVID-19 presents multiscaling properties. Our model has been tested on Italy's real data of the Center for Systems Science and Engineering (CSSE) at JohnsHopkins University. The estimate of the proposed SIRD model to susceptible, infected, recovered, and death populations were more suitable because it incorporates two scaling exponents, i.e., a variable decaying rate as the power-law kernel to consider different memory strengths. In this manner, the implications of the Italian government's quarantine decree were accurately approximated since the time-dependent fractional-orders fit well the evolution of the SIRD populations. Without loss of generality, many sets of fractional-orders may be used according to the countries' specific restrictions.

Since different governmental measures (quarantine decrees, mobility restrictions, lockdowns, home office, homeschooling, vaccination, etc.) and adherence of the people to those measures leads to strikingly different growth rates and scaling exponents, the proposed approach will allow studying those aspects that could change the propagation speed of COVID-19 by taking into account the memory effects at specific stages (growth, slow down, and saturation). In this manner, the proposed multi-fractional approach could help forecast the total number of infections more efficiently to prepare doctors and health care systems for different scenarios. In future studies, the trade-off between the control strategies and the number of infected people will be analyzed.

Authors declare that they have no conflict of interest. puter resources, technical advice and support provided by Laboratorio Nacional de Supercómputo del Sureste de México (LNS ), a member of the CONACYT national laboratories, with project no. 202001025C .

Modeling and prediction of the 2019 coronavirus disease spreading in China incorporating human migration data

A model for SARS-COV-2 infection with treatment

Optimal policies for control of the novel coronavirus (COVID-19)

Prediction of bifurcations by varying critical parameters of COVID-19

Modeling and prediction of COVID-19 in Mexico applying mathematical and computational models

Mathematical biology: I. an introduction

Memory effects on epidemic evolution: the susceptible-infected-recovered epidemic model

Analysis and forecast of COVID-19 spreading in China, Italy and France

SEIR modeling of the COVID-19 and its dynamics

Variable infectiousness in HFV transmission models

Non-Markov stochastic dynamics of real epidemic process of respiratory infections

Fractional order mathematical modeling of COVID-19 transmission

Non-Markovian infection spread dramatically alters the susceptible-infected-susceptible epidemic threshold in networks

Applicability of time fractional derivative models for simulating the dynamics and mitigation scenarios of COVID-19

Simulating non-Markovian stochastic processes

Dynamics of a fractional order mathematical model for COVID-19 epidemic

Clarify the physical process for fractional dynamical systems

The role of fractional calculus in modeling biological phenomena: a review

A mathematical model for COVID-19 transmission by using the Caputo fractional derivative

Novel fractional order SIDARTHE mathematical model of COVID-19 pandemic

Forecast analysis of the epidemics trend of COVID-19 in the USA by a generalized fractional-order SEIR model

SEIR epidemic model for COVID-19 transmission by Caputo derivative of fractional order

Analysis of Caputo fractional-order model for COVID-19 with lockdown

Haar wavelet collocation approach for the solution of fractional order COVID-19 model using Caputo derivative

A fractional-order model for the novel coronavirus (COVID-19) outbreak

Fractional-order susceptible-infected model: definition and applications to the study of COVID-19 main protease

A fractional-order SEIHDR model for COVID-19 with inter-city networked coupling effects

Some fractal thoughts about the COVID-19 infection outbreak

Fractal kinetics of COVID-19 pandemic

Scaling features in the spreading of COVID-19

Fractal signatures of the COVID-19 spread

Short-term predictions of country-specific COVID-19 infection rates based on power law scaling exponents

Can India develop herd immunity against COVID-19?

The suppression of scale-free fMRI brain dynamics across three different sources of effort: aging, t ask novelty and t ask difficulty

Recognition of COVID-19 disease from X-ray images by hybrid model consisting of 2D curvelet transform, chaotic salp swarm algorithm and deep learning technique

Nonlinear time series and principal component analyses: potential diagnostic tools for COVID-19 auscultation

An interactive web-based dashboard to track COVID-19 in real time

Short-time fractal analysis of biological autoluminescence

Low-dimensional chaos and fractal properties of long-term sunspot activity

Fractal mechanisms and heart rate dynamics: long-range correlations and their breakdown with disease

Coexistence of chaotic and complexity dynamics of fluctuations with long-range temporal correlations under typical condition for formation of multiple anodic double layers in dc glow discharge plasma

Characterizing chaos and multifractality in noise-assisted tumor-immune interplay

Vargovi ć E . Changes in the hurst exponent of heartbeat intervals during physical activity

Some comments on hurst exponent and the long memory processes on capital markets

Performance analysis of hurst exponent estimators using surrogate-data and fractional lognormal noise models: Application to breathing signals from preterm infants

Multi-scaling in finance

Multifractality in asset returns: theory and evidence

Assessment of long-range correlation in time series: how to avoid pitfalls

Estimating long-range dependence: finite sample properties and confidence intervals

Applications of variable-order fractional operators: a review

A review on variable-order fractional differential equations: mathematical foundations, physical models, numerical methods and applications

New variable-order fractional chaotic systems for fast image encryption

Generalized Taylor's formula

Global existence theory and chaos control of fractional differential equations

This work was funded by National Council of Science and Technology (CONACYT/Mexico), Proyecto Apoyado por el Fondo Sectorial de Investigación para la Educación under the project no. 258880 . Also, the Deanship of Scientific Research (DSR) at King Abdulaziz University , Jeddah, Saudi Arabia funded this project, under grant no. ( FP-205-42 ) . JMMP thankfully acknowledges com-