Fractional Derivatives Applied to Epidemiology

ABSTRACT

We seek investigate the use of fractional derivatives, both analytically and through simulations. We present some models and perform investigations about them, starting with the classic model and the basic definitions to discuss difficulties in constructing a non-artificial fractional model. Also, we analyze the COVID-19 pandemic using a fractional epidemiological SIR model carefully constructed and present numerical results using MATLAB.

Keywords:
SIR model; Fractional Derivatives; COVID-19

1 INTRODUCTION

Throughout history, epidemics have devastated a large percentage of humanity. Important questions can be worked out through mathematical modeling: mathematicians study their evolution and try to elucidate essential questions, such as when the peak of the disease is expected and how many people will be infected in total.

The coronavirus disease 2019 (COVID-19) is a respiratory disease that spreads from person to person. As of July 28, 2020, the reported cases of COVID-19 on the planet exceeded 16 million ²⁵. Person-to-person spread is thought to occur mainly via respiratory droplets, as the spread of influenza ²². This suggests that the use of the epidemiological model SIR (Susceptible - Infected - Removed) is reasonable for understanding the spread of COVID-19.

As fractional calculus has proven to be a faithful tool in capturing the dynamics of the physical process of many scientific objects, being its most striking features the memory effect, here we revisit much of the mathematical and epidemiological theory in order to improve understanding of fractional models. The tool has already been used to model the dynamics of the new coronavirus in interesting works, such as ^{21), (31}. In this work, we are concerned with the precise construction of a fractional model following the seminal ideas of the SIR model creators. We analyze problems and difficulties that occur when simply replacing derivatives of whole order with derivatives of fractional orders and correcting the dimension artificially, such as, for example, the unexplained change in the total population. It was in Angstmann, Henry & McGann ⁵ that we found solid steps that could explain why and where fractional derivatives can arise, paying attention to the precise definition of each parameter and the dimensioning. We revisited the authors’ work and, in addition to starting the discussion of the reproduction number for their model, we applied the model to the data of the Brazilian and Italian COVID-19 pandemic.

We intend that the work will attract readers’ attention to the care with the use of fractional derivatives, presenting a non-artificial construction for an epidemiological model and its application. For this aim, in Chapter 2 we offer a re-presentation of the classic SIR model proposed by Kermack & McKendrick in 1927 ²⁰, continuing with special attention to the case of constant parameters and the definition of the number of reproduction of an infectious disease. In Chapter 3, we offer fractional calculus preliminaries and establish difficulties in defining the fractional model. In Chapter 4, we present a physical derivation of a fractional model, following the steps of Angstmann, Henry & McGann ⁵, where they use the probabilistic language of the Continuous-Time Random Walks (CTRW) and the Riemann-Liouville fractional derivative. Finally, in Chapter 5 we display several numerical results.

2 THE SIR MODEL

In 1927, the SIR (Susceptible-Infected-Removed) model was introduced in a remarkable way in a work proposed by Kermack & McKendrick ²⁰. The authors concluded that, although the causative agent does not lose its infectiousness and the population is not entirely infected, the end of an epidemic may result from a special relationship between population density, infectivity and rates of recovery and death. We present below a brief adaptation of this work.

2.1 General Theory

Initially, we uniformly discretize time considering $\Delta T=1$ and we assume that people are infected only when passing from one interval to another, not during the interval itself. We denote the number of infected at time t who have been infected for θ intervals by v _t,θ . The total number of infected at time t is therefore $I_t=\sum _{\theta =0}^t{v}_{t,\theta }$. The notation v _t is also used to indicate the rate of new infections in time t. In general, $v_t=\frac{v_{t,0}}{\Delta T}$, except at the origin, when a population of infectious I ₀, regardless of the model to be developed, is inserted into the total population. Thus, $v_{0,0}=v_0\Delta T+I_0$.

If $\psi \left(\theta \right)={\psi }_{\theta }$ denotes the rate of removal of the infectious compartment at the age of infection θ (that is, the sum of recovery and death rates), so the number removed from each θ − group is given by ${\psi }_{\theta }{v}_{t,\theta }=\frac{v_{t,\theta }-v_{t+1,\theta +1}}{\Delta T}$. Therefore, it follows that $v_{t,\theta }=v_{t-1,\theta -1}(1-\psi (\theta -1\left)\Delta T\right)$. Proceeding, we obtain

$v_{t,\theta }=v_{t-\theta ,0}{B}_{\theta },$ (2.1)

where $B_{\theta }=(1-\psi (\theta -1\left)\Delta T\right)(1-\psi (\theta -2\left)\Delta T\right)\cdots (1-\psi (0\left)\Delta T\right)$. It is important to note that, in order not to overload the notation, we write time $t+1,\theta -2$, etc., instead of $t+\Delta T,\theta -2\Delta T$ and so on.

Now, if $\varphi \left(\theta \right)={\varphi }_{\theta }$ denotes the infectivity rate of an infected person who has been infected for θ stages, the rate of new infected v _t should be equal to $S_t\sum _0^t{\varphi }_{\theta }{v}_{t,\theta }$, where S _t denotes the number of people not yet infected/ immunized at time t. Clearly,

$S_t=N-\sum _0^t{v}_t\Delta T-I_0,$ (2.2)

where N is the total population. If R _t denotes the number of removed (by recovery or death), then, disregarding vital dynamics processes, $S_t+I_t+R_t=N$. We have

$I_t=\sum _0^t{v}_{t,\theta }=\sum _0^t{B}_{\theta }{v}_{t-\theta }\Delta T+B_t{I}_0.$ (2.3)

Also note that, from Eq. 2.1, we have

$v_t=S_t\sum _1^t{\varphi }_t{v}_{t,\theta }=S_t\left(\sum _1^t{A}_{\theta }{v}_{t-\theta }\Delta T+A_t{I}_0\right),$ (2.4)

where we define $A_{\theta }={\varphi }_{\theta }{B}_{\theta }$ and assume ${\varphi }_0=0$, that is, a person is not infective at the moment of the infection. By other hand, for $t>0$, we have $v_t=\frac{S_t-S_{t+1}}{\Delta T}$. Therefore, it follows from Eq. 2.4 that

$\frac{S_t-S_{t+1}}{\Delta T}=S_t\left(\sum _1^t{A}_{\theta }{v}_{t-\theta }\Delta T+A_t{I}_0\right).$ (2.5)

Finally, we note that

$\frac{R_{t+1}-R_t}{\Delta T}=\left(\sum _1^t{C}_{\theta }{v}_{t-\theta }\Delta T+C_t{I}_0\right),$ (2.6)

where we define $C_{\theta }={\psi }_{\theta }{B}_{\theta }$.

Allowing $∆T\to 0$, we get the relationship $v_t=-d{S}_t/dt$ and the three equations that define the SIR model

$\frac{d{S}_t}{dt}=-S_t\left({\int }_0^t{A}_{\theta }{v}_{t-\theta }d\theta +A_t{I}_0\right);$ (2.7)

$I_t={\int }_0^t{B}_{\theta }{v}_{t-\theta }d\theta +B_t{I}_0;$ (2.8)

$\frac{d{R}_t}{dt}={\int }_0^t{C}_{\theta }{v}_{t-\theta }d\theta +C_t{I}_0,$ (2.9)

where we have, by the property of the product integral ³², $B_{\theta }=\exp \left(-{\int }_0^{\theta }\psi \left(a\right)da\right)$. It is observed that the upper limit t must be divided by the unit of time considered, becoming scalar. The reader will be able to verify that $\frac{d{S}_t+d{I}_t+d{R}_t}{dt}=0$, which means that the population is kept constant, regardless of the functions ψ, ϕ.

2.2 Constant rates

The analysis of the case in which the ϕ, ψ rates are constants β, γ, respectively, provides insights for understanding epidemics. In this case, the system 2.7 - 2.9 becomes, omitting the t index for simplicity

$\frac{dS}{dt}=-\beta SI;$ (2.10)

$\frac{dI}{dt}=\beta SI-\gamma I;$ (2.11)

$\frac{dR}{dt}=\gamma I.$ (2.12)

So $dR/dt=\gamma (N-S-R)$. In addition, $dS/dR=-(\beta /\gamma )S$, a separable equation that leads to the result $\log (S_0/S)=(\beta /\gamma )R$. Therefore, $S=S_0\exp (-(\beta /\gamma \left)R\right)$, following that

$\frac{dR}{dt}=\gamma (N-S_0\exp (-(\beta /\gamma )R)-R).$ (2.13)

From Eq. 2.13, we can obtain the total number of infected throughout the course of the disease by studying $d{R}_{\infty }/dt=0$, which is equivalent to $N-S_0\exp (-(\beta /\gamma \left)R_{\infty }\right)-R_{\infty }=0$. Considering $S_0\ne 0$, we get

$R_{\infty }=N+\frac{\gamma }{\beta }{W}_0\left(\frac{-S_0\beta e^{-N\beta /\gamma }}{\gamma }\right),$ (2.14)

where W ₀ represents the principal branch of the Lambert W function ⁹, not used by the classical authors ²⁰. We note that if the disease is early, $S_0\approx N$ and the epidemic does not occur if $N\le \frac{\gamma }{\beta }$, that is, if $N\frac{\beta }{\gamma }\le 1$. In fact, in this case $W_0\left(\frac{-S_0\beta e^{-N\beta /\gamma }}{\gamma }\right)\approx W_0\left(\frac{-N\beta }{\gamma }{e}^{-N\beta /\gamma }\right)=-N\beta /\gamma$ e $R_{\infty }\approx 0$. The constant $N\frac{\beta }{\gamma }$ is equivalent, in this model, to the so-called basic reproduction number, to be presented in the next section. It is possible to obtain some other results for small epidemics, which are not of interest to us here, but have been studied in the main reference ²⁰.

2.3 The Basic Reproduction Number

The reproduction number ℜ reflects how infectious a disease is in a given context. This constant is often confused with the basic reproduction number, ℜ₀ , which represents the reproduction number when there is no immunity or deliberate intervention in the transmission of the disease. On the other hand, ℜ is the effective reproductive number and is modified both by the decreasing of susceptible population and by the implementation of interventions and the infectivity changes that the virus may suffer. Under conditions of homogeneous population, ℜ is defined as the average number of infections that a single individual can generate during his infectious period ⁷. In traditional models, epidemics of an infection cannot occur when the ℜ₀ is less than 1 and established outbreaks will disappear if interventions or depletion of the susceptible part of the population are sufficient to keep ℜ below 1. Also according to ⁷, it is important to keep in mind that the ℜ₀ of an infection depends on the population analyzed, because the rates of contact between people may differ due to differences in population density and culture. The ℜ at a given time can also vary from place to place, as communities may differ in their levels of immunity and intervention.

3 APPLICATIONS OF FRACTIONAL DERIVATIVES TO SIR TYPE MODELS

Currently, fractional calculus has proven to be a faithful tool in capturing the dynamics of the physical process of many scientific objects, such as biological or ecological phenomenas and control systems ³⁴. One of its most striking features is the memory effect. Over the years, varied fractional operators have been defined with different focuses and applications. Although it is not used in this work, we highlight the recent introduction of the ψ-Hilfer fractional derivative, a fractional derivative with respect to a function ψ which incorporates a large class of fractional derivatives as particular cases ^{10), (11), (12}.

The fractional calculus has already been used to model the dynamics of the new coronavirus in interesting works, such as ^{21), (31}. Here, we related the memory effect to potential laws regarding the infectiousness (transmission) of the infected patient and the removal rate. That is, the longer infected, the less infectious and the more likely he will be removed from the infectious compartment. This effect is important and is not captured by classic models. We intend to study the orders of derivatives that best reproduce the data of the pandemic considering the variation in infectivity and removal rate. It is important to mention that the memory effect decreases when the fractional orders tend to 1 ³³.

3.1 Preliminaries of Fractional Calculus

The Mittag-Leffler functions are of important use in the theory of fractional calculus. We present the following definition ⁸:

Definition 1. (Two and Three-parameter Mittag-Leffler function) Let $z\in \mathrm{ℂ}$ , with $\alpha ,\beta ,\rho \in \mathrm{ℂ}$ three-parameters such that $Re\left(\alpha \right)>0,Re\left(\beta \right)>0,Re\left(\rho \right)>0$ . We define the three-parameter Mittag-Leffler function via the power series

$E_{\alpha ,\beta }^{\rho }\left(z\right)=\sum _{k=0}^{\infty }\frac{(\rho {)}_k}{\Gamma (\alpha k+\beta )}\frac{z^k}{k!},$ (3.1)

where (ρ) _k is the Pochhammer symbol. Particularly, when $\rho =1$ , we have the two-parameter Mittag-Leffler function, denoted simply by $E_{\alpha ,\beta }\left(z\right)$ . Note that this function generalizes the exponential function, being equal to the exponential when $\alpha =\beta =1$ .

The following results are worth:

Proposition 1. The Laplace transform of the function $t^{\beta -1}{E}_{\alpha ,\beta }^{\rho }\left(a{t}^{\alpha }\right)$ is given by:

$\mathcal{L}\left[t^{\beta -1}{E}_{\alpha ,\beta }^{\rho }\right(a{t}^{\alpha }\left)\right]\left(s\right)=s^{-\beta }(1-a{s}^{-\alpha }{)}^{-\rho },$ (3.2)

where$Re\left(s\right)>0$and$\left|a{s}^{-\alpha }\right|<1$.

Proposition 2. (Three-parameter Mittag-Leffler function derivative) The following identity is valid:

$\frac{d^k}{d{x}^k}{E}_{\alpha ,\beta }^{\rho }\left(z\right)=\left(\rho {)}_k{E}_{\alpha ,\beta +\alpha k}^{\rho +k}\right(z).$ (3.3)

The proofs of the propositions can be found at ⁸. Below we consider $[a,b]\subset \mathrm{ℝ},\alpha >0$, a function $f\in L^p[a,b],p\ge 1$, and $n-1<\alpha <n$, with $n\in \mathrm{\mathbb{N}}$. Also, Γ is the gamma function.

Definition 2. The Riemann-Liouville fractional integral of order α is defined for $t\in [a,b]$ by:

${}_{a}J_t^{\alpha }f\left(t\right)=\frac{{\displaystyle 1}}{{\displaystyle \Gamma \left(\alpha \right)}}{\int }_a^t(t-\theta {)}^{\alpha -1}f(\theta )d\theta .$ (3.4)

Proposition 3. (Fractional integral of the three-parameter Mittag-Leffler function ⁸ ) We have the following identity:

${}_{0}J_t^{\nu }\left[t^{\beta -1}{E}_{\alpha ,\beta }^{\rho }\right(a{t}^{\alpha }\left)\right]=t^{\nu +\beta -1}{E}_{\alpha ,\nu +\beta }^{\rho }\left(a{t}^{\alpha }\right).$ (3.5)

Definition 3. The Riemann-Liouville fractional derivative of order α is defined for $t\in [a,b]$ by:

${}_a^{R-L}{D}_t^{\alpha }f\left(t\right)=\frac{1}{\Gamma \left(n-\alpha \right)}\left(\frac{d^n}{d{t}^n}\right){\int }_a^t(t-\theta {)}^{n-\alpha -1}f(\theta )d\theta .$ (3.6)

In other words,

${}_a^{R-L}{D}_t^{\alpha }f\left(t\right)=D^n{[}_a{J}_t^{n-\alpha }f\left(t\right)],$ (3.7)

with D ⁿ representing the integer order derivative.

Definition 4. The Caputo fractional derivative of order α is defined for $t\in [a,b]$ by:

${}_a{}^{C}D_t^{\alpha }f\left(t\right)=\frac{1}{\Gamma (n-\alpha )}{\int }_a^t(t-\theta {)}^{n-\alpha -1}\frac{d^n}{d{\theta }^n}f(\theta )d\theta .$ (3.8)

In other words,

${}_a^C{D}_t^{\alpha }f\left(t\right)={}_a{J}_t^{n-\alpha }\left[D^{n}f\right(t\left)\right].$ (3.9)

The Riemann-Liouville fractional integral is the inverse to the right of both the Riemann-Liouville and the Caputo fractional derivatives ²⁶, that is,

${}_a{}^{C}D_t^{\alpha }{}_a{J}_t^{\alpha }f\left(t\right)={}_a^{R-L}{D}_t^{\alpha }{}_a{J}_t^{\alpha }f\left(t\right)=f\left(t\right).$ (3.10)

The Definition 2 can be seen as an application of the Laplace transform and the Convolution Theorem. In what follows, we use the notation J ^α for ${}_{0}J_t^{\alpha }$. Using Laplace transform properties, we can write, to $n\in \mathrm{\mathbb{N}},\mathcal{L}\left\{\right(J^{n}f\left)\right(t\left)\right\}=s^{-n}\mathcal{L}\left\{f\right(t\left)\right\}$. Therefore, it is reasonable to assume that, for fractional α, it also holds:

$\left(J^{\alpha }f\right)\left(t\right)={\mathcal{L}}^{-1}\left(s^{-\alpha }\mathcal{L}\left\{f\right(t\left)\right\}\right).$ (3.11)

Considering $p\left(t\right)=t^{\alpha -1}$, we have $\mathcal{L}\left\{p\right(t\left)\right\}=s^{-\alpha }\Gamma \left(\alpha \right)$. Therefore, by the convolution equation and by the linearity of the transform:

$\left(J^{\alpha }f\right)\left(t\right)=\frac{1}{\Gamma \left(\alpha \right)}{\mathcal{L}}^{-1}\{\mathcal{L}\{p\left(t\right)\}\mathcal{L}\{f\left(t\right)\left\}\right\}=\frac{1}{\Gamma \left(\alpha \right)}(p*f)=\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t(t-\theta {)}^{\alpha -1}f(\theta )d\theta ,$ (3.12)

as in Definition 2. Is important to note that when the fractional order tends to the integer order, the fractional derivative tends to the same result as the ordinary derivative. This is known by backward compatibility and this property for the Gru¨nwald-Letnikov, Riemann-Liouville and Caputo fractional derivatives can be seen on ²⁷.

3.2 Difficulties in Defining the Fractional Model

We studied some models in which authors generalize classic first order ODE’s by replacing the integer derivative on the left side with fractional derivatives, usually the Caputo one, once it allows conventional initial conditions and its derivative of a constant is bounded (specifically, equal to 0) ²⁸. We found that these models can produce very good estimates, as well as interesting equations from a mathematical point of view. However, we are interested in the following question: does the change in the order of derivatives automatically establish consistent models, with respect to the definition of parameters, units and balance? In 2010, an interesting article was published ¹³ proving that, in general, this cannot happen.

In the literature on the use of fractional derivatives in SIR-type models, we find some more common versions. Based on the model 2.10 - 2.12 or its various extensions, some authors (eg. ^{3), (30}) replace the derivative on the left side with the Caputo fractional derivative of order α, keeping the same parameters as the model with an integer derivative. For example, for the SIR model, it is considered

$D^{\alpha }S\left(t\right)=-\beta I\left(t\right)S\left(t\right)$ (3.13)

$D^{\alpha }I\left(t\right)=\beta I\left(t\right)S\left(t\right)-\gamma I\left(t\right)$ (3.14)

$D^{\alpha }R\left(t\right)=\gamma I\left(t\right),$ (3.15)

when β and γ are positive constants defined as in the system 2.10-2.12. The first difficulty noted with this definition concerns units. We note that D ^α f(t) has unit [time]^−α . Therefore, the right side of the equation must have this unit. Just abolishing the units does not solve the problem. According to ¹³, the units actually help to reveal the issue: a 1 order rate, for example, is a different type of that of order 1/2, in the same way that, in classical physics, the rate of order 1, velocity, is different from that corresponding to the order 2, the acceleration.

Other authors (eg. ^{2), (19}) correct the dimension by multiplying time constants or raising the parameters to α. In this last case, the fractional SIR model could be written as follows

$D^{\alpha }S\left(t\right)=-{\beta }^{\alpha }I\left(t\right)S\left(t\right)$ (3.16)

$D^{\alpha }I\left(t\right)={\beta }^{\alpha }I\left(t\right)S\left(t\right)-{\gamma }^{\alpha }I\left(t\right)$ (3.17)

$D^{\alpha }R\left(t\right)={\gamma }^{\alpha }I\left(t\right).$ (3.18)

Finally, we read authors (eg. ^{1), (15), (16), (17}) who, intending to extend the fractional models, seek a system with multiple orders, which, if on the one hand is of great mathematical interest, on the other hand makes the analysis even more complicated. More than the parameter definitions and the dimensional difficulty, the relationship between different orders can influence the total mass or population. The article ¹³ demonstrates, in a simple two-compartment system, that the use of different fractional orders without proper care leads to the violation of mass balance. The authors defend, in an article of the same year ¹⁴, that, both in the context of the pharmacokinetics then considered and in analog systems, an outgoing mass flow defined as a fractional order rate cannot appear as an inflow into another compartment, with a different fractional order, without violating the mass balance.

We will analyze the relationship between orders and the total population N of the fractional SIR model defined as

$D^{{\alpha }_1}S\left(t\right)=-{\beta }_{1}I\left(t\right)S\left(t\right)$ (3.19)

$D^{{\alpha }_2}I\left(t\right)={\beta }_{2}I\left(t\right)S\left(t\right)-{\gamma }_{2}I\left(t\right)$ (3.20)

$D^{{\alpha }_3}R\left(t\right)={\gamma }_{3}I\left(t\right).$ (3.21)

The fractional derivative is in Caputo sense and ${\alpha }_1,{\alpha }_2,{\alpha }_3\in (0,1]$, as most of the authors we had access considered. We define the parameters in some way that balances the units, for instance, ${\beta }_i={\beta }^{{\alpha }_i}$ and ${\gamma }_i={\gamma }^{{\alpha }_i}$ for $i\in \{1,2,3\}$. Applying Laplace transform and performing operations on the equations 3.19 - 3.21, we obtain

$\mathcal{L}\{S+I+R\}=\frac{S_0+I_0+R_0}{s}+\left[\frac{{\beta }_2}{s^{{\alpha }_2}}-\frac{{\beta }_1}{s^{{\alpha }_1}}\right]\mathcal{L}\left\{SI\right\}+\left[\frac{{\gamma }_3}{s^{{\alpha }_3}}-\frac{{\gamma }_2}{s^{{\alpha }_2}}\right]\mathcal{L}\left\{I\right\}.$ (3.22)

Applying the inverse transform, we finally write

$N=N_0+\left[\frac{{\beta }_2{t}^{{\alpha }_2-1}}{\Gamma \left({\alpha }_2\right)}-\frac{{\beta }_1{t}^{{\alpha }_1-1}}{\Gamma \left({\alpha }_1\right)}\right]*SI+\left[\frac{{\gamma }_3{t}^{{\alpha }_3-1}}{\Gamma \left({\alpha }_3\right)}-\frac{{\gamma }_2{t}^{{\alpha }_2-1}}{\Gamma \left({\alpha }_2\right)}\right]*I,$ (3.23)

That is,

$\begin{array}{c}N=N_0+{\int }_0^t\left[\frac{{\displaystyle {\beta }_2(t-\theta {)}^{{\alpha }_2-1}}}{{\displaystyle \Gamma \left({\alpha }_2\right)}}-\frac{{\displaystyle {\beta }_1(t-\theta {)}^{{\alpha }_1-1}}}{{\displaystyle \Gamma \left({\alpha }_1\right)}}\right]S\left(\theta \right)I\left(\theta \right)d\theta \\ +{\int }_0^t\left[\frac{{\displaystyle {\gamma }_3(t-\theta {)}^{{\alpha }_3-1}}}{{\displaystyle \Gamma \left({\alpha }_3\right)}}-\frac{{\displaystyle {\gamma }_2(t-\theta {)}^{{\alpha }_2-1}}}{{\displaystyle \Gamma \left({\alpha }_2\right)}}\right]I\left(\theta \right)d\theta \end{array}.$ (3.24)

Note that, if ${\alpha }_1={\alpha }_2={\alpha }_3$, the integrals in Eq. 3.24 become nulls, that is, the population remains constant, but the same may not be true in other cases. For example, if ${\alpha }_1<{\alpha }_2<{\alpha }_3$, the population goes up to large t. Conversely, if ${\alpha }_1>{\alpha }_2>{\alpha }_3$, the population decreases to large t. The big problem is not to change the population, which can happen due to vital dynamics, migration etc., but to justify the change. In this case, without migration or dynamics, people could not disappear or appear as illustrated in the Fig. 1-2.

Figure 1
Oscilation of N(t)

Figure 2:
The R compart. exceed N ₀.

Finally, we note that Eq. 3.19, 3.20 and 3.21 have different parameters. This indicates that the ℜ₀ of the SIR model, classically given by the constant $N\frac{\beta }{\gamma }$, is not trivial on the fractional case.

4 DERIVATION OF THE FRACTIONAL MODEL

In the previous section, we presented fractional models commonly used by authors in the field and made some considerations. Still other questions arise: can we replace the entire order with fractional orders on the left based on some mathematical theorem? If there are two different orders (say, coming from potential laws on infectiousness and recovery/removal), these should be considered in the equations of which compartments among the three or four used? In addition to the adequacy of the units, is there an epidemiological explanation for multiplying time constants or raising the parameters to the order used? We have not yet obtained precise justifications, although we have verified, through published works and simulations that we carry out, that these models are capable of representing reality with a good degree of precision.

As already mentioned, we believe that the emergence of fractional derivatives in models similar to the original comes from considering the functions of infectiousness and removal derived from potential laws. Next, we present a physical derivation of a fractional model, following the steps of Angstmann, Henry & McGann ⁵, where they use the probabilistic language of the Continuous- Time Random Walks (CTRW).

4.1 Model Derivation

Consider an individual infected since the time t'. If there are S(t) susceptible in time t, this infected person has a probability $S\left(t\right)/N$ that his contact is susceptible, considering the population homogeneous. Therefore, in the period of t to $t+∆T$ , the expected number of new infections per infected individual is given by $\sigma (t,t\text{'})\frac{S\left(t\right)}{N}\Delta T$ . The transmission rate per infectious individual $\sigma (t,t\text{'})$ depends both on the age of the infection, $t-t\text{'}$, and on the present time, t, which, among others, influences the amount of contacts of the infectious.

The probability that an individual infected at the moment t' is still infected at the moment t is given by the survival function $\Phi (t,t\text{'})$. Therefore, the flux of individuals to the I compartment at time t is recursively given by

$q^{+}(I,t)={\int }_{-\infty }^t\sigma (t,t\text{'})\frac{S\left(t\right)}{N}\Phi (t,t\text{'})q^{+}(I,t\text{'})dt\text{'}.$ (4.1)

The initial condition is obtained by the number of individuals infected at the time 0 and considering the time in which each individual has become infected. This is given by the function $i(-t\text{'},0)$ which represents the number of individuals who are still infectious at time 0 and who were originally infected at some point earlier $t\text{'}<0$. Then $q^{+}(I,t\text{'})=i(-t\text{'},0)/\Phi (0,t\text{'})$ for $t\text{'}<0$. For simplicity, we consider $i(-t,0)=i_0\delta (-t)$, where $\delta \left(t\right)$ is the Dirac delta function. So,

$q^{+}(I,t)={\int }_0^t\sigma (t,t\text{'})\frac{S\left(t\right)}{N}\Phi (t,t\text{'})q^{+}(I,t\text{'})dt\text{'}+i_0\sigma (t,0)\frac{S\left(t\right)}{N}\Phi (t,0).$ (4.2)

The infection rate $\sigma (t,t\text{'})$ is assumed to be a function of both the current time (due, for example, to containment measures), having an extrinsic infectivity ω, and the age of infection $t-t\text{'}$, having an intrinsic infectivity ρ. So we can write

$\sigma (t,t\text{'})=\omega \left(t\right)\rho (t-t\text{'}).$ (4.3)

Assuming that the natural death and the removal of an infected individual are independent processes, we can write the survival function as

$\Phi (t,t\text{'})=\varphi (t-t\text{'})\theta (t,t\text{'}),$ (4.4)

where $\varphi (t-t\text{'})$ is the probability that an individual infected since t ^' has not yet recovered or been killed by the disease at time t. Also, $\theta (t,t\text{'})$ is the probability that an infected individual since t' has not yet died of natural death (that is, independent of the disease) until time t. The θ function is constructed similarly to the B function in Section 2.1, given by

$\theta (t,t\text{'})=e^{-{\int }_{t\text{'}}^t\gamma \left(u\right)du},$ (4.5)

where γ is the death rate. Notice that

$\theta (t,t\text{'})=\theta (t,u)\theta (u,t\text{'}),\forall t\text{'}<u<t.$ (4.6)

Individuals in the I compartment at time t must have entered this compartment at some time before and remained in it until t. Therefore, we can express the number of infected individuals as follows:

$I\left(t\right)=I_0\left(t\right)+{\int }_0^t\Phi (t,t\text{'})q^{+}(I,t\text{'})dt\text{'}.$ (4.7)

The $I_0\left(t\right)$ function provides the number of individuals who were infected at 0 and remain infectious at t. This function can be written as follows:

$I_0\left(t\right)={\int }_{-\infty }^0\frac{\Phi (t,t\text{'})}{\Phi (0,t\text{'})}i(-t\text{'},0)dt\text{'}=\frac{\Phi (t,0)}{\Phi (0,0)}{i}_0=\Phi (t,0)i_0.$ (4.8)

With this information we can build the principal equations of the model. We start by deriving Eq. 4.7 through the Leibniz Rule, obtaining

$\begin{array}{r}\frac{dI\left(t\right)}{dt}=q^{+}(I,t)-{\int }_0^t\psi (t-t\text{'})\theta (t,t\text{'})q^{+}(I,t\text{'})dt\text{'}\\ -\gamma \left(t\right){\int }_0^t\varphi (t-t\text{'})\theta (t,t\text{'})q^{+}(I,t\text{'})dt\text{'}+\frac{d{I}_0\left(t\right)}{dt},\end{array}$ (4.9)

where $\psi \left(t\right)=-\frac{d\varphi \left(t\right)}{dt}$. This ψ has an important relationship with the continuous random variable X that provides the time of removal of the individual from the infectious compartment. The cumulative distribution of X, namely F defined by $F\left(t\right)=P(X\le t)$, is such that $F\left(t\right)=1-\varphi \left(t\right)$. Therefore, the probability density function of X is $\psi \left(t\right)=-\frac{d\varphi \left(t\right)}{dt}$.

Using Eq. 4.2 - 4.7, Eq. 4.9 can be rewritten as

$\begin{array}{r}\frac{dI\left(t\right)}{dt}=\omega \left(t\right)\left(\frac{S\left(t\right)}{N}{\int }_0^t\rho (t-t\text{'})\Phi (t,t\text{'})q^{+}(I,t\text{'})dt\text{'}+\rho \left(t\right)\Phi (t,0)i_0\right)\\ -{\int }_0^t\psi (t-t\text{'})\theta (t,t\text{'})q^{+}(I,t\text{'})dt\text{'}-i_0\psi \left(t\right)\theta (t,0)-\gamma \left(t\right)I\left(t\right)\end{array}.$ (4.10)

In what follows, we will remove the dependency of $q^{+}(I,t\text{'})$ from Eq. 4.10 by defining infectivity and recovery memory kernels. We will henceforth consider $i_0=1$ for simplicity. We can rewrite Eq. 4.7 as

$\frac{I\left(t\right)}{\theta (t,0)}-\varphi \left(t\right)={\int }_0^t\varphi (t-t\text{'})\frac{q^{+}(I,t\text{'})}{\theta (t\text{'},0)}dt\text{'}.$ (4.11)

The right side is now in convolution form. Taking the Laplace transform, we get

$\mathcal{L}\left\{\frac{I\left(t\right)}{\theta (t,0)}-\varphi \left(t\right)\right\}=\mathcal{L}\left\{\varphi \right(t\left)\right\}\mathcal{L}\left\{\frac{q^{+}(I,t)}{\theta (t,0)}\right\}.$ (4.12)

The first integral of Eq. 4.10 is the product of $\theta (t,0)$ with an integral that can be rewritten using the Laplace transform in the form $\mathcal{L}\left\{\rho \right(t\left)\varphi \right(t\left)\right\}\mathcal{L}\left\{\frac{q^{+}(I,t)}{\theta (t,0)}\right\}$. We can write, by Eq. 4.12,

$$\begin{gathered} \mathcal{L}\left\{\rho \right(t\left)\varphi \right(t\left)\right\}\mathcal{L}\left\{\frac{q^{+}(I,t)}{\theta (t,0)}\right\}=\frac{\mathcal{L}\left\{\rho \right(t\left)\varphi \right(t\left)\right\}}{\mathcal{L}\left\{\varphi \right(t\left)\right\}}\mathcal{L}\left\{\frac{I\left(t\right)}{\theta (t,0)}-\varphi \left(t\right)\right\} \\ =\mathcal{L}\left\{{\int }_0^t{K}_I(t-t\text{'})\left(\frac{I(t\text{'})}{\theta (t\text{'},0)}\right)dt\text{'}-\rho \left(t\right)\varphi \left(t\right)\right\}, \end{gathered}$$ (4.13)

where we define the infectivity memory kernel by:

$K_I\left(t\right)={\mathcal{L}}^{-1}\left\{\frac{\mathcal{L}\left\{\rho \right(t\left)\varphi \right(t\left)\right\}}{\mathcal{L}\left\{\varphi \right(t\left)\right\}}\right\}.$ (4.14)

Similarly, the second integral of Eq. 4.10 is the product of $\theta (t,0)$ with an integral that can be rewritten using Laplace transform in the form

$\mathcal{L}\left\{\psi \right(t\left)\right\}\mathcal{L}\left\{\frac{q^{+}(I,t)}{\theta (t,0)}\right\}=\frac{\mathcal{L}\left\{\psi \right(t\left)\right\}}{\mathcal{L}\left\{\varphi \right(t\left)\right\}}\mathcal{L}\left\{\frac{I\left(t\right)}{\theta (t,0)}-\varphi \left(t\right)\right\}\\ =\mathcal{L}\left\{{\int }_0^t{K}_R(t-t\text{'})\left(\frac{I(t\text{'})}{\theta (t\text{'},0)}\right)dt\text{'}-\psi \left(t\right)\right\},$ (4.15)

where we define the recovery memory kernel by

$K_R\left(t\right)={\mathcal{L}}^{-1}\left\{\frac{\mathcal{L}\left\{\psi \right(t\left)\right\}}{\mathcal{L}\left\{\varphi \right(t\left)\right\}}\right\}.$ (4.16)

Replacing Eq. 4.13 and 4.15 in Eq. 4.10, we get

$\begin{array}{r}\frac{dI\left(t\right)}{dt}=\omega \left(t\right)\frac{S\left(t\right)}{N}\left(\theta (t,0){\int }_0^t{K}_I(t-t\text{'})\left(\frac{I(t\text{'})}{\theta (t\text{'},0)}\right)dt\text{'}-\theta (t,0)\rho \left(t\right)\varphi \left(t\right)+\rho \left(t\right)\Phi (t,0)\right)\\ -\theta (t,0){\int }_0^t{K}_R(t-t\text{'})\left(\frac{I(t\text{'})}{\theta (t\text{'},0)}\right)dt\text{'}+\theta (t,0)\psi \left(t\right)-\psi \left(t\right)\theta (t,0)-\gamma \left(t\right)I\left(t\right)\end{array}.$ (4.17)

Finally, we simplify the set of equations for the SIR model:

$\frac{dS\left(t\right)}{dt}=\gamma \left(t\right)N-\omega \left(t\right)\frac{S\left(t\right)}{N}\theta (t,0){\int }_0^t{K}_I(t-t\text{'})\frac{I(t\text{'})}{\theta (t\text{'},0)}t\text{'}dt\text{'}-\gamma \left(t\right)S\left(t\right);$ (4.18)

$\begin{array}{r}\frac{dI\left(t\right)}{dt}=\omega \left(t\right)\frac{S\left(t\right)}{N}\theta (t,0){\int }_0^t{K}_I(t-t\text{'})\frac{I(t\text{'})}{\theta (t\text{'},0)}t\text{'}dt\text{'}\\ -\theta (t,0){\int }_0^t{K}_R(t-t\text{'})\frac{I(t\text{'})}{\theta (t\text{'},0)}dt\text{'}-\gamma \left(t\right)I\left(t\right);\end{array}$ (4.19)

$\frac{dR\left(t\right)}{dt}=\theta (t,0){\int }_0^t{K}_R(t-t\text{'})\frac{I(t\text{'})}{\theta (t\text{'},0)}dt\text{'}-\gamma \left(t\right)R\left(t\right),$ (4.20)

where we consider the same rate γ(t) of natural mortality in each compartment, with the birth rate equal to that. In this case the population remains constant.

We incorporate fractional derivatives into the model by choosing ψ(t) with potential law and ρ(t) related to the choice of ψ(t). In particular, we use the Mittag-Leffler function, according to Definition 1:

$\psi \left(t\right)=\frac{t^{\alpha -1}}{{\tau }^{\alpha }}{E}_{\alpha ,\alpha }\left(-{\left(\frac{t}{\tau }\right)}^{{}^{\alpha }}\right),$ (4.21)

for $0<\alpha \le 1$, where τ is a scale parameter. The asymptotic expansion of the Mittag-Leffler function ⁶ allows to affirm that $\psi \left(t\right)=O\left(t^{-1-\alpha }\right)$ for $t\to \infty$.

It is not difficult to prove that $E_{\alpha ,\alpha }\left(-{\left(\frac{t}{\tau }\right)}^{\alpha }\right)=\alpha E_{\alpha ,1+\alpha }^2\left(-{\left(\frac{t}{\tau }\right)}^{\alpha }\right)$. Thus, by Eq. 3.3, the correspondent survival function is

$\varphi \left(t\right)=E_{\alpha ,1}\left(-{\left(\frac{t}{\tau }\right)}^{\alpha }\right).$ (4.22)

Using Eq. 3.2, we can calculate the Laplace transform of the recovery memory kernel with Mittag-Leffler distribution, obtaining

$\mathcal{L}\left\{K_R\right(t\left)\right\}=\frac{\mathcal{L}\left\{\psi \right(t\left)\right\}}{\mathcal{L}\left\{\varphi \right(t\left)\right\}}=s^{1-\alpha }{\tau }^{-\alpha }.$ (4.23)

On the other hand, we have $\mathcal{L}{\{}_0^{R-L}{D}_t^{1-\alpha }f\left(t\right)\}=s^{1-\alpha }\mathcal{L}\{f\left(t\right)\left\}{-}_0{J}^{\alpha }f\right(t){\mid }_{t=0^{+}}$. We consider here ${}_{0}J^{\alpha }f\left(t\right){\mid }_{t=0^{+}}=0$, as in ⁴. Hereafter, for simplicity, we write D ^1−α in place of ${}_0{}^{R-L}D^{1-\alpha }$. From these considerations, a convolution with the recovery memory kernel can be written as

${\int }_0^t{K}_R(t-t\text{'})\frac{I(t\text{'})}{\theta (t\text{'},0)}dt\text{'}={\tau }^{-\alpha }{D}^{1-\alpha }\left(\frac{I\left(t\right)}{\theta (t,0)}\right).$ (4.24)

A fractional derivative can also be incorporated into the infectivity memory kernel. Note that Eq. 4.14 has a Laplace transform similar to Eq. 4.23 and it is interesting to consider

$\rho \left(t\right)=\frac{1}{\varphi \left(t\right)}\frac{t^{\beta -1}}{{\tau }^{\beta }}{E}_{\alpha ,\beta }\left(-{\left(\frac{t}{\tau }\right)}^{{}^{\alpha }}\right).$ (4.25)

As $\rho \left(t\right)\ge 0$ is required, we must have $0<\alpha \le \beta \le 1$. Using Eq. 4.25, we obtain the Laplace transform of the infectivity kernel:

$\mathcal{L}\left\{K_I\right(t\left)\right\}=s^{1-\beta }{\tau }^{-\beta },$ (4.26)

Following thar

${\int }_0^t{K}_I(t-t\text{'})\frac{I(t\text{'})}{\theta (t\text{'},0)}dt\text{'}={\tau }^{-\beta }{D}^{1-\beta }\left(\frac{I\left(t\right)}{\theta (t,0)}\right).$ (4.7)

Replacing Eq. 4.24 and 4.27 in the model 4.18 - 4.20, we obtain a fractional SIR model:

$\frac{dS\left(t\right)}{dt}=\gamma \left(t\right)N-\frac{\omega \left(t\right)S\left(t\right)\theta (t,0)}{N{\tau }^{\beta }}{D}^{1-\beta }\left(\frac{I\left(t\right)}{\theta (t,0)}\right)-\gamma \left(t\right)S\left(t\right);$ (4.28)

$\frac{dI\left(t\right)}{dt}=\frac{\omega \left(t\right)S\left(t\right)\theta (t,0)}{N{\tau }^{\beta }}{D}^{1-\beta }\left(\frac{I\left(t\right)}{\theta (t,0)}\right)-\frac{\theta (t,0)}{{\tau }^{\alpha }}{D}^{1-\alpha }\left(\frac{I\left(t\right)}{\theta (t,0)}\right)-\gamma \left(t\right)I\left(t\right);$ (4.29)

$\frac{dR\left(t\right)}{dt}=\frac{\theta (t,0)}{{\tau }^{\alpha }}{D}^{1-\alpha }\left(\frac{I\left(t\right)}{\theta (t,0)}\right)-\gamma \left(t\right)R\left(t\right).$ (4.30)

Note that if $\alpha =\beta =1$ and $\gamma \left(t\right)\equiv \gamma ,\omega \left(t\right)\equiv \omega$ are considered constant, we get the model 2.10 - 2.12, added only by the vital dynamics:

$\frac{dS\left(t\right)}{dt}=\gamma N-\frac{\omega S\left(t\right)I\left(t\right)}{N\tau }-\gamma S\left(t\right);$ (4.31)

$\frac{dI\left(t\right)}{dt}=\frac{\omega S\left(t\right)I\left(t\right)}{N\tau }-\frac{1}{\tau }I\left(t\right)-\gamma I\left(t\right);$ (4.32)

$\frac{dR\left(t\right)}{dt}=\frac{1}{\tau }I\left(t\right)-\gamma R\left(t\right).$ (4.33)

Finally, note that the probability distribution $F\left(t\right)=1-\varphi \left(t\right)$ is a Mittag-Leffler distribution $F(t;\alpha ,\tau )=1-E_{\alpha }\left(-{\left(\frac{t}{\tau }\right)}^{\alpha }\right)$. If $\alpha =\beta =1$, we have an exponential distribution and then (first moment) of the random variable X exists, with τ being exactly the average recovery time.

4.1.1 The Reproduction Number

The effective reproduction number ℜ(t') at time t' can be understood as the expected number of individuals that will be infected by an individual infected since t'. So, we can calculate

$\mathfrak{R}(t\text{'})={\int }_{t\text{'}}^{\infty }\sigma (t,t\text{'})\frac{S\left(t\right)}{N}\Phi (t,t\text{'})dt={\int }_{t\text{'}}^{\infty }\omega \left(t\right)\rho (t-t\text{'})\varphi (t-t\text{'})\theta (t,t\text{'})\frac{S\left(t\right)}{N}dt.$ (4.34)

If we consider that during the infeccious period of the individual we have $\frac{S\left(t\right)}{N}\approx \frac{S(t\text{'})}{N}$, we obtain from Eq. 4.25

$\mathfrak{R}(t\text{'})=\frac{S(t\text{'})}{N}{\int }_{t\text{'}}^{\infty }\omega \left(t\right)\theta (t,t\text{'})\frac{(t-t\text{'}{)}^{\beta -1}}{{\tau }^{\beta }}{E}_{\alpha ,\beta }\left(-{\left(\frac{(t-t\text{'})}{\tau }\right)}^{\alpha }\right)dt.$ (4.35)

As an example, we consider $\gamma \left(t\right)\equiv \gamma >0$ and an exponential decay of ω over time: $\omega \left(t\right)=\omega ·e^{-t/a}$. So,

$\omega \left(t\right)\theta (t,t\text{'})=\omega ·e^{-(\gamma +1/a)(t-t\text{'})}{e}^{-t\text{'}/a}.$ (4.36)

Therefore,

$\mathfrak{R}(t\text{'})=\frac{S(t\text{'})\omega }{N\tau }{e}^{-t\text{'}/a}{\int }_{t\text{'}}^{\infty }{e}^{-(\gamma +1/a)(t-t\text{'})}{\left(\frac{(t-t\text{'})}{\tau }\right)}^{\beta -1}{E}_{\alpha ,\beta }\left(-{\left(\frac{(t-t\text{'})}{\tau }\right)}^{\alpha }\right)dt.$ (4.37)

We simplify it by changing $t:=t-t\text{'}$. Verifying that the integral is in the form of a Laplace Transform, we use Eq. 3.2, finally obtaining

$\mathfrak{R}(t\text{'})=\frac{S(t\text{'})\omega }{N{\tau }^{\beta }}{e}^{-t\text{'}/a}\frac{(\gamma +1/a{)}^{\alpha -\beta }}{(\gamma +1/a{)}^{\alpha }+{\tau }^{-\alpha }}.$ (4.38)

Particularly, if $S_0\approx N$,

${\mathfrak{R}}_0=\frac{\omega (\gamma +1/a{)}^{\alpha -\beta }}{{\tau }^{\beta }(\gamma +1/a{)}^{\alpha }+{\tau }^{\beta -\alpha }}.$ (4.39)

We note that if $a\to \infty$, then $\omega \left(t\right)\approx \omega$ for every t and the expression of ℜ₀ becomes the expression given by ⁵ for the case $\omega \left(t\right)\equiv \omega$ constant,

${\mathfrak{R}}_0=\frac{\omega {\gamma }^{\alpha -\beta }}{{\tau }^{\beta }{\gamma }^{\alpha }+{\tau }^{\beta -\alpha }}.$ (4.40)

It is interesting to note that, in fractional order simulations with different orders or variable ω, the peak of infectious may no longer occurs when $\mathfrak{R}(t\text{'})=1$. Moreover, it is very important to states that the approximation $\frac{S\left(t\right)}{N}\approx \frac{S(t\text{'})}{N}$ is really crude and should be explored in our future works.

5 APPLICATIONS TO THE COVID-19 PANDEMIC

We constructed a slightly modification of the L1-scheme ²⁴ to discretize the fractional SIR model 4.28-4.30. Then, according to available data from the COVID-19 pandemic on Brazil and Italy, we apply through MATLAB the Least Squares Method using the FDIPA ¹⁸ algorithm in order to minimize the function $f(\omega ,\tau )=\left|\right|AI-I|{|}^2+||AR-R|{|}^2$, where AI is a vector formed by the data of the infected being monitored and, AR, by the data of those removed (recovered + deaths).

5.1 Discretization of the fractional order derivative

The fractional order $1-\alpha$ is assumed in the range $[0,1)$. The time interval $[0,t]$ is discretized as $0=t_0<t_1<\cdots <t_n=t$, where the time steps ∆T have the same size. Since we wanted to compare the model to the daily data, we choose each time step with size $∆T=1$. In particular, $t_i=i\Delta T=i$ (unit [time]) for all $i\in \{0,1,\cdots ,n\}$, so we performed the following discretization for the Riemann-Liouville derivative, with unit [time]^1−α

$$\begin{gathered} D^{1-\alpha }f\left(t_j\right)=\frac{1}{\Gamma (\alpha +1)}\left(\frac{d}{d{t}_j}\right){\int }_0^{t_j}(t_j-\theta {)}^{\alpha -1}f(\theta )d\theta \\ \simeq \frac{1}{\Gamma (\alpha +1)}\left[\sum _{k=0}^{j-1}f\left(k\right)\left[\right(j-k+1{)}^{\alpha }-2(j-k{)}^{\alpha }+(j-k-1{)}^{\alpha }]+f(j)\right] \\ =\frac{1}{\Gamma (\alpha +1)}\left[\sum _{k=0}^{j-1}\left[f\right(k+1)-f(k\left)\right]\left[\right(j-k{)}^{\alpha }-(j-k-1{)}^{\alpha }]+f\left(0\right)\left[\right(j+1{)}^{\alpha }-j^{\alpha }]\right]. \end{gathered}$$ (5.1)

5.2 Parameters

In the 4.28 - 4.30 system, the parameter ω(t) is related to contagion, being strongly affected by government measures, such as quarantine, as well as hygiene and readaptation measures in times of pandemic, both social and physical. Thus, we consider three situations in the simulations: 1) Constant value $\omega \left(t\right)\equiv \omega$, as considered in most models; 2) Exponential decay in value over time: $\omega \left(t\right)=\omega ·e^{-t/a}$; 3) Exponential decay with oscillations, representing relaxation and periodic hardening of the containment measures: $\omega \left(t\right)=\omega ·e^{-t/a}·{\cos }^2(t\pi /b)$.

The Brazilian population is given by $N=210147125$ inhabitants and the initial stage by $(S,I,R)=(N-1,1,0)$. As said, applying the Least Squares Method through the FDIPA algorithm we seek to minimize the function $f(\omega ,\tau )=\left|\right|AI-I|{|}^2+||AR-R|{|}^2$. It is noted that the absence of data about recovered in the first two months of the Brasilian pandemic prejudices the model.

5.3 Results

Initially we fixed $\alpha =\beta =1$ and $\omega \left(t\right)\equiv \omega$, considering a classic SIR model with constant parameters. The results are obtained in Fig. 3 (comparison between model and real data) and Fig. 4 (projection of the peak of infection). The data was updated until 07/25/2020 and the ℜ showed is the effective reproduction number on this date with the approximations 4.37-4.38.

Figure 3:
Model X Real data.

Figure 4:
Projection of the peak.

We observed that, even though ℜ₀ is small, the I and R curves do not tend to equilibrium as fast as the real curves. We believe that this is due to the fact that we disregard the laws of infection and recovery and, above all, because we take constant ω.

Next, in Fig.5-10, we consider ω with exponential decay, $a=70$ and different orders α, β . We emphasize that, as α decreases, the right tail of the I compartment becomes heavier. In Fig.11 and Fig.12, we take an oscilating exponential decay for ω, with $a=140$,$b=110$ and orders $\alpha =0.8$,$\beta =0.9$. In this test we approximate the ℜ₀ as in the last case, but we do not calculate the ℜ.

Figure 5:
Model X Real data.

Figure 6:
Projection of the peak.

Figure 7:
Model X Real data.

Figure 8:
Projection of the peak.

Figure 9:
Model X Real data.

Figure 10:
Projection of the peak.

Figure 11:
Model X Real data.

Figure 12:
Projection of the peak.

Finally, we used the orders $\alpha =0.9$,$\beta =0.95$ to carry out a test in relation to the Italian pandemic, which has the important characteristic of having its curve of infected individuals declining since the last week of April. The results of the comparisons between the models and the actual data ²³ are shown in Fig. 13 and 14. For Italy, we use $a=40$ and $\gamma =2.92·{10}^{-5}$ ³⁵. The available data has visible outliers.

Figure 13:
Model X Real data.

Figure 14:
Projection.

6 DISCUSSION AND CONCLUSIONS

We investigated the use of fractional derivatives, both analytically and through simulations. We present models and perform investigations about them, discussing difficulties and differences between classic and fractional models. Also, we analyzed the COVID-19 pandemic using a fractional epidemiological SIR model carefully built and performed a numerical solution using discretization and implementation in MATLAB. We observed that, among others effects, the use of fractional orders allows us accentuate the asymmetry of the curve (slow decline). This effect is being observed in the real data of many countries ²⁹.

We reiterate that the work is not in disagreement with any research cited. On the contrary, we respect the models obtained by the literature review, believing that they contribute with equations of great interest to Mathematics. Also, it is possible that some authors write their equations as mentioned in section 3.2 for reasons that we do not perceive or that are in literature that we do not access. Our goal is only to present questions and calculations that we believe to be of both pure and applied interest, as the definition of parameters and the units and population balances.

We emphasize that the choice of the Riemann-Liouville fractional derivative was due to the application that arose during the construction of the model, with attention to Eq.4.23. In the future we want to extend the mathematical discussion on the parameters effects (including the cases when α is close to 0, what corresponds to too heavy right tail on I compartment), more accurate discretizations and also on the use of different types of fractional operators. Moreover, we intend to discuss equilibrium points and better calculations for the reproduction number. A better understanding in this area can be useful not only for epidemiology, but also for the other applied areas in which fractional derivatives are used, such as thermoelasticity and diffusion.

Acknowledgments

We appreciate the support of the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Financing Code 001.

REFERENCES

• ¹
  A. Ahmad et al. Dynamical behavior of SIR epidemic model with non-integer time fractional derivatives: A mathematical analysis. Int. J. Adv. Appl. Sci, 5(1) (2018), 123-129.
• ²
  R. Almeida. Analysis of a fractional SEIR model with treatment. Applied Mathematics Letters, 84 (2018), 56-62.
• ³
  I. Ameen. “Fractional Calculus: Numerical Methods and SIR Models”. Ph.D. thesis, University of Padova, Padova, Italy (2017). PhD thesis.
• ⁴
  C.N. Angstmann, B.I. Henry & A.V. McGann. A fractional-order infectivity SIR model. Physica A: Statistical Mechanics and its Applications, 452 (2016), 86-93.
• ⁵
  C.N. Angstmann, B.I. Henry &A.V. McGann . A fractional-order infectivity and recovery SIR model. Fractal and Fractional, 1(1) (2017), 11.
• ⁶
  H. Bateman. “Higher transcendental functions”, volume 3. McGraw-Hill Book Company (1953).
• ⁷
  N.G. Becker et al. Using mathematical models to assess responses to an outbreak of an emerged viral respiratory disease. In “Final Report to the Australian Government Department of Health and Ageing. National Centre for Epidemiology and Population Health, Australian National University” (2006).
• ⁸
  E. Capelas de Oliveira. “Solved Exercises in Fractional Calculus”. Springer (2019).
• ⁹
  R.M. Corless et al. On the Lambert W function. Advances in Computational mathematics, 5(1) (1996), 329-359.
• ¹⁰
  J.V. da Costa Sousa & E. Capelas de Oliveira. On the ψ-Hilfer fractional derivative. Communications in Nonlinear Science and Numerical Simulation, 60 (2018), 72-91.
• ¹¹
  J.V. da Costa Sousa & E. Capelas de Oliveira. Leibniz type rule: ψ-Hilfer fractional operator. Communications in Nonlinear Science and Numerical Simulation , 77 (2019), 305-311.
• ¹²
  J.V. da Costa Sousa, G. Frederico & E. Capelas de Oliveira. ψ-Hilfer pseudo-fractional operator: new results about fractional calculus (2020).
• ¹³
  A. Dokoumetzidis, R. Magin & P. Macheras. A commentary on fractionalization of multicompartmental models. Journal of pharmacokinetics and pharmacodynamics, 37(2) (2010), 203-207.
• ¹⁴
  A. Dokoumetzidis R. Magin & P. Macheras. Fractional kinetics in multi-compartmental systems. Journal of pharmacokinetics and pharmacodynamics , 37(5) (2010), 507-524.
• ¹⁵
  M. Farman et al. Dynamical Behavior of Hepatitis B Fractional-order Model with Modeling and Simulation. Journal of Biochemical Technology, 10(3) (2019), 11.
• ¹⁶
  A.C.F.N. Gomes & A. de Cezaro. Modelo SIRC fracionário com múltiplas ordens para influenza. Scientia Plena, 15(4) (2019).
• ¹⁷
  S. Hasan et al. Solution of fractional SIR epidemic model using residual power series method. Applied Mathematics and Information Sciences, 13(2) (2019), 153-161.
• ¹⁸
  J. Herskovits. Feasible direction interior-point technique for nonlinear optimization. Journal of optimization theory and applications, 99(1) (1998), 121-146.
• ¹⁹
  M.R., Islam A. Peace D. Medina & T. Oraby. Integer versus fractional order seir deterministic and stochastic models of measles. International Journal of Environmental Research and Public Health, 17(6) (2020), 2014.
• ²⁰
  W.O. Kermack & A.G. McKendrick. Contributions to the mathematical theory of epidemics-I. 1927. (1991).
• ²¹
  M.A. Khan & A. Atangana. Modeling the dynamics of novel coronavirus (2019-nCov) with fractional derivative. Alexandria Engineering Journal, (2020).
• ²²
  K. McIntosh M.S., Hirsch & A. Bloom. Coronavirus disease 2019 (COVID-19). UpToDate Hirsch MS Bloom, 5 (2020).
• ²³
  Microsoft/Bing. COVID-19-Data (2020). URL URL https://raw.githubusercontent.com/microsoft/Bing-COVID-19-Data/master/data/Bing-COVID19-Data.csv [Online; accessed 29-July-2020].
  » https://raw.githubusercontent.com/microsoft/Bing-COVID-19-Data/master/data/Bing-COVID19-Data.csv
• ²⁴
  K. Oldham & J. Spanier. “The fractional calculus theory and applications of differentiation and integration to arbitrary order”. Elsevier (1974).
• ²⁵
  W.H. Organization & N. Authorities. “Coronavirus disease 2019 (COVID-19): situation report, 190”. World Health Organization (2020).
• ²⁶
  M. Ortigueira & J. Machado. Fractional definite integral. Fractal and Fractional , 1(1) (2017), 2.
• ²⁷
  M.D. Ortigueira & J.T. Machado. What is a fractional derivative? Journal of computational Physics, 293 (2015), 4-13.
• ²⁸
  I. Podlubny. “Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications”. Elsevier (1998).
• ²⁹
  , R. Ranjan. Estimating the final epidemic size for Covid-19 outbreak using improved epidemiological models. medRxiv, (2020).
• ³⁰
  D. Rostamy & E. Mottaghi. Numerical solution and stability analysis of a nonlinear vaccination model with historical effects. Hacettepe Journal of Mathematics and Statistics, 47(6) (2017), 1478-1494.
• ³¹
  A.S. Shaikh I.N. Shaikh & K.S. Nisar. “A mathematical model of covid-19 using fractional derivative: Outbreak in india with dynamics of transmission and control”. Preprints (2020).
• ³²
  A., Slavík. “Product integration, its history and applications”. Matfyzpress Prague (2007).
• ³³
  H. Sun, W. Chen, H. Wei & Y. Chen. A comparative study of constant-order and variable-order fractional models in characterizing memory property of systems. The European Physical Journal Special Topics, 193(1) (2011), 185.
• ³⁴
  K. Wang & J. Huang. High order fast algorithm for the Caputo fractional derivative. arXiv preprint arXiv:1705.06101, (2017).
• ³⁵
  Wikipedia contributors. Demographics of Italy - Wikipedia, The Free Encyclopedia (2020). URL URL ttps://en.wikipedia.org/w/index.php?title=Demographics_of_Italy&oldid=969810317 [Online; accessed 29-July-2020].
  » ttps://en.wikipedia.org/w/index.php?title=Demographics_of_Italy&oldid=969810317

Publication Dates

• Publication in this collection
  02 July 2021
• Date of issue
  May-Aug 2021

History

• Received
  31 July 2020
• Accepted
  17 Dec 2020