Fractional-Order SIR Epidemic Model for Transmission Prediction of COVID-19 Disease

Abstract

In this paper, the fractional-order generalization of the susceptible-infected-recovered (SIR) epidemic model for predicting the spread of the COVID-19 disease is presented. The time-domain model implementation is based on the fixed-step method using the nabla fractional-order difference defined by Grünwald-Letnikov formula. We study the influence of fractional order values on the dynamic properties of the proposed fractional-order SIR model. In modeling the COVID-19 transmission, the model’s parameters are estimated while using the genetic algorithm. The model prediction results for the spread of COVID-19 in Italy and Spain confirm the usefulness of the introduced methodology.

1. Introduction

Coronavirus disease 2019 (COVID-19) is an acute infectious disease of the respiratory system that is caused by the SARS-CoV-2 virus infection. The World Health Organization (WHO) issued the official name and explained it as CO for corona, VI for virus, D for disease, and 19 for the year of first identification. The infection symptoms are usually a fever, dry cough, fatigue, loss of appetite, and shortness of breath. COVID-19 was first recognized and described on December 8, 2019, in the Wuhan, Republic of China. As a result of rapidly spreading the virus into other countries and continents, on March 11, 2020, WHO classified COVID-19 as a pandemic. With both mortality and infectiousness rates, COVID-19 is recognized as the worst pandemic in modern times since the flu pandemic of the early 20th century.

Mathematical modeling theories of infectious disease transmission dynamics have prompted considerable research interest, as the worldwide issue of the COVID-19 outbreak has attracted scientists from different areas [1,2,3]. The Susceptible-Infected-Removed (SIR) epidemic model is a classical compartmental model that is widely used to predict the progress of infectious disease [4,5,6]. Its origins date to the fundamental work by Kermack and McKendrick in the early 20th century [7].

Although the SIR model may be useful in modeling several epidemic processes, it may not be sufficient to describe the COVID-19 spread dynamics. Therefore, there are some modifications to the SIR model used in the approximation of COVID-19 disease [8]. Kucharski et al. [9] used the SEIR (susceptible-exposed-infectious-removed) compartment model in order to study the transmission and control of COVID-19. Sarkar et al. [10] proposed an SEIR model extended by adding contact tracing and other interventions, such as quarantine, lockdown, and social distancing. The SIDARTHE (susceptible-infected-diagnosed-ailing-recognized-threatened-healed-extinct) model with differentiating between diagnosed and non-diagnosed individuals has been presented by Giordano et al. [11] in order to analyze the COVID-19 epidemic in Italy. In [12], Ndaïrou et al. introduced the SEIPAHRF (susceptible-exposed-symptomatic infectious-super spreaders-asymptomatic infectious-hospitalized-recovered-fatality) model that considers the class of super-spreaders in order to study the transmission dynamics for the outbreak that occurred in Wuhan, China.

On the other hand, it is well known that the implementation of fractional calculus in differential equations of various processes, which lead to fractional-order generalizations of the modeling approaches, can significantly improve the modeling performance in several processes. Hence, many scientists have paid attention to this mathematical analysis branch, as it may give a better understanding of some complex physical phenomena, where ordinary calculus fails to provide a precise explanation. Atangana and Gómez-Aguilar [13] applied fractional derivative to model chaos and statistics. Ghanbari and Atangana [14] used the Atangana—Baleanu definition for the fractional derivative to introduce the new fractional mask for image denoising. Morales-Delgado et al. [15] obtained solutions of a cancer chemotherapy effect model involving the Caputo-Fabrizio and Atanagana-Baleanu fractional derivatives. There are some papers considering the fractional-order derivate implemented in the compartment models. A fractional-order extended SEIR (susceptible-exposed-infectious-removed) model considering interaction among people and the reservoir of the infection (seafood market) has been presented and studied by Khan and Atangana [16]. Srivastava and G$\ddot{\mathrm{n}}$erhan [17] used the conformable fractional differential transform method to compute an approximate solution of the fractional-order SIR epidemic model of childhood disease. Abdo et al. [18] presented a nonlinear fractional-order model under the fractional Adam-BashForth method. In this work, we will implement the Grünwald-Letnikov fractional-order derivative into the classical SIR model in order to simulate and predict the COVID-19 outbreak.

The paper is organized, as follows. Having introduced the problem of modeling the COVID-19 transmission in Section 1, Section 2 presents the ’classical’ SIR model. The fractional-order generalization of the SIR model and its fixed-step implementation are introduced in Section 3. Section 4 presents the influence of fractional orders on model dynamical properties and presents virus transmission modeling in Spain and Italy. The conclusions of Section 5 complete the paper.

2. Preliminaries

The SIR model is one of the simplest compartment models, and many others are derivatives of its basic form. The total population N is divided into three epidemiological subpopulations: $\left(S\right)$ The number of susceptible individuals. When a susceptible and an infectious person come into contact, the susceptible individual contracts the disease and transitions to the infectious compartment, $\left(I\right)$ The number of infectious individuals. These individuals have been infected and they can infect susceptible individuals, $\left(R\right)$, the number of removed individuals. This subpopulation consists of individuals who are immune due to natural immunity, have recovered from the disease and gained immunity, or have died. The SIR model consists of three ordinary differential equations that are presented below.

$$\begin{array}{cc}\hfill \frac{dS}{dt}=& -\frac{\beta IS}{N}\hfill \\ \hfill \frac{dI}{dt}=& \frac{\beta IS}{N}-\gamma I\hfill \\ \hfill \frac{dR}{dt}=& \gamma I\hfill \end{array}$$

(1)

where $\beta$ and $\gamma$ are infection and recovery rate, respectively.

Birth and death dynamics are usually neglected in simple compartmental models, since the dynamics of an epidemic, such as the flu, are often much faster than the dynamics of birth and death, and the numerical simulations are often executed for short time intervals.

Note that the dynamics of infectious disease depends on the following ratio:

$$R_0=\frac{\beta }{\gamma }$$

(2)

where $R_0$ is the basic reproduction number. This ratio quantifies the expected number of cases that are directly generated by one individual in an infected population. When $R_0>1$, the infection will be able to start spreading, causing a potential outbreak. If $R_0<1$ the disease will cease to exists.

3. Fractional-Order Generalization of SIR Model

It is well known that the integer-order derivative can be generalized by the fractional-order counterpart. In this case, the derivative is presented in the form of ${\mathcal{D}}^{\alpha }$, where the order $\alpha \in \mathbb{R}$. Note that the derivative is usually defined in one of three definitions, e.g., Caputo, Rieman-Liouville, and Grünwald-Letnikov ones. In the paper, Grünwald-Letnikov [19] is exploited hereinafter, where

$${\mathcal{D}}^{\alpha }x\left(t\right)=h^{-\alpha }\underset{h\to 0_{+}}{lim}\sum _{j=0}^{\lfloor \frac{t}{h}\rfloor }{(-1)}^j\left(\genfrac{}{}{0pt}{}{\alpha }{j}\right)x\left(t-\frac{j}{h}\right)$$

(3)

with t is continuous-time, h is sampling interval and $\left(\genfrac{}{}{0pt}{}{.}{.}\right)$ denoting newton binomial, which (for $j>0$) can be calculated using the Pochhammer function as

$$\left(\genfrac{}{}{0pt}{}{\alpha }{j}\right)=\frac{{\alpha }_j}{j!}=\frac{{\prod }_{i=0}^{j-1}(\alpha -i)}{j!}=\frac{a_{j-1}}{(j-1)!}\frac{\alpha -j+1}{j}$$

(4)

In application, we use fixed step method, where we use nabla difference defined in discrete-time $\mathbb{T}=\{0,1,2,...\}$, such that $t=\mathbb{T}h$. The nabla difference is defined as

$$\begin{array}{ccc}\hfill {\nabla }^{\alpha }x\left(t\right)& =& h^{-\alpha }\sum _{j=0}^{\lfloor \frac{t}{h}\rfloor }{(-1)}^j\left(\genfrac{}{}{0pt}{}{\alpha }{j}\right)x\left(t-\frac{j}{h}\right)\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill {\nabla }^{\alpha }{x}_k& =& h^{-\alpha }\sum _{j=0}^k{(-1)}^j\left(\genfrac{}{}{0pt}{}{\alpha }{j}\right)x_{k-j}\hfill \end{array}$$

(6)

Note that, for satisfactory low h, we have ${\mathcal{D}}^{\alpha }x\left(t\right)\approx {\nabla }^{\alpha }x\left(t\right)$, so that we can effectively use the discrete-time difference to model continuous time processes.

It is well known that the incorporation of fractional-order derivative into model equations can lead to improving model parameters in several physical processes. In particular, various diffusion processes can be better modeled by use of fractional-order differential equations [20,21,22]. Taking into account that virus transmission is a similar class of process to the above, the fractional-order generalization of the SIR model can significantly improve the model performances. In this case, the integer-order differences are substituted by the fractional-order counterparts ${\alpha }_i\in \mathbb{R}$, $i\in \{1,2,3\}$, and the fractional-order SIR model can be presented in form

$$\begin{array}{cc}\hfill {\mathcal{D}}^{{\alpha }_1}S\left(t\right)=& -\frac{\beta I\left(t\right)S\left(t\right)}{N}\hfill \\ \hfill {\mathcal{D}}^{{\alpha }_2}I\left(t\right)=& \frac{\beta I\left(t\right)S\left(t\right)}{N}-\gamma I\left(t\right)\hfill \\ \hfill {\mathcal{D}}^{{\alpha }_3}R\left(t\right)=& \gamma I\left(t\right)\hfill \end{array}$$

(7)

where ${\alpha }_i\in (0,2)$, $i\in \{1,2,3\}$, $\beta$ and $\gamma$ are as in Equation (1). The fractional-order generalization of the SIR model dynamics is different from the ’classical’ counterpart. The fractional-order model’s several properties are different from in the integer-order one, e.g., the total population number N is time-various phenomena. Therefore, the SIR model’s fractional-order generalization should be considered to be a heuristics model with no direct dependency on physical phenomena. Moreover, it is worth emphasizing that, in this case, we have three additional parameters for estimation $\{{\alpha }_1,{\alpha }_2,{\alpha }_3\}$.

Taking into account that, in most countries, COVID-19 statistics are presented once per day, to the prediction of the virus transmission, we can use discrete-time nabla difference as in Equation (5) where sampling period $h=1$, so we have

$$\begin{array}{cc}\hfill {\nabla }^{{\alpha }_1}{S}_{k+1}=& -\frac{\beta I_k{S}_k}{N}\hfill \\ \hfill {\nabla }^{{\alpha }_2}{I}_{k+1}=& \frac{\beta I_k{S}_k}{N}-\gamma I_k\hfill \\ \hfill {\nabla }^{{\alpha }_3}{R}_{k+1}=& \gamma I_k\hfill \end{array}$$

(8)

Now, combining Equations (6) and (8), we can easily present model (8) in the form of

$$\begin{array}{cc}\hfill S_{k+1}=& \left(\frac{-\beta I_k}{N}+{\alpha }_1\right)S_k-\sum _{j=2}^{k+1}{(-1)}^j\left(\genfrac{}{}{0pt}{}{{\alpha }_1}{j}\right)S_{k-j}\hfill \\ \hfill I_{k+1}=& \left(\frac{\beta S_k}{N}-\gamma +{\alpha }_2\right)I_k-\sum _{j=2}^{k+1}{(-1)}^j\left(\genfrac{}{}{0pt}{}{{\alpha }_2}{j}\right)I_{k-j}\hfill \\ \hfill R_{k+1}=& \gamma I_k-\sum _{j=2}^{k+1}{(-1)}^j\left(\genfrac{}{}{0pt}{}{{\alpha }_3}{j}\right)R_{k-j}\hfill \end{array}$$

(9)

which can be easily used in the simulation of the process.

In order to model the COVID-19 transmission, we have model parameters $\{{\alpha }_1,{\alpha }_2,{\alpha }_3,\gamma ,\beta ,N\}$. The parameter N is an arbitrary given, but the rest parameters have to be selected in the optimization process, as described in the next part of the paper.

In the next section, we show the influence of the fractional-order ${\alpha }_i$, $i\in \{1,2,3\}$ on the model dynamical properties, and we present the model performance in the modeling of COVID-19 transmission in Italy and Spain.

4. Numerical Simulations

This section gives insight into the dynamics of the presented fractional-order SIR epidemic model (9) and shows how the fractional-order difference alters the numerical results. Firstly, we will present simulations with different fractional-orders to demonstrate how their values affect the compartment populations. Next, we will show numerical results based on the real data for Italy and Spain (see report hubs [23,24]).

4.1. Dynamical Properties of Fractional-Order SIR Model

Consider three sets of simulations with following parameters: $S_0=950,I_0=45,R_0=5$, $N=1000,\beta =0.4,\gamma =0.1,t=60$ and various ${\alpha }_{1,2,3}=\{0.2,0.4,0.6,0.8\}$. In each set, only one ${\alpha }_i$ changes while two remain are equal to 1. That way, we can see how the particular fractional order ${\alpha }_i$, $i\in \{1,2,3\}$ affects the model. Figure 1, Figure 2 and Figure 3 present the simulation results.

As we can see in Figure 1, Figure 2 and Figure 3, fractional orders ${\alpha }_1$ and ${\alpha }_2$ affect the whole model, while ${\alpha }_3$ only affects the removed compartment. It can be expected if we look at the model’s Equation (9). With more parameters, the presented model is more complicated but is also more versatile. The results presented above show that fractional-order SIR can be more adjustable in the modeling of real-life processes.

4.2. Transmission Modeling for Italy and Spain

In this section, we implement the fractional-order SIR model in order to predict COVID-19 transmission for Italy and Spain.

In order to model the spread of COVID-19, we have model parameters, as follows $\{{\alpha }_1,{\alpha }_2,{\alpha }_3,\gamma ,\beta ,N\}$. Accounting that the order ${\alpha }_3$ does not affect model output I, for simplicity, we assume that ${\alpha }_3=1$. Parameter N is arbitrary given for each country. Therefore, using square lost function, we have to estimate four parameters $\{{\alpha }_1,{\alpha }_2,\gamma ,\beta \}=argmin{\sum }_{k=0}^K{(\widehat{I_k}-I_k)}^2$, where $\widehat{I_k}$ and $I_k$ denotes modeled and actual number of infected individuals, respectively. In this paper, optimal model parameters are obtained the use of the Monte Carlo method based on the genetic algorithm.

The starting point of simulation regarding the Italy epidemic is 21 February, when there were 19 active cases and 0 removed, so $I_0=19$, $R_0=0$, $N=$ 60,461,828 and $S_0$ = 60,461,726. Using the above-presented optimization method for $K=150$, we obtain parameters, as follows $\beta =3.3108,\gamma =0.1333,{\alpha }_1=0.4091$ and ${\alpha }_2=0.9964$. Figure 4 presents the time plots of modeled and actual disease transmission.

As we can see in Figure 4, the fractional-order SIR model is effective in modeling the COVID-19 transmission in Italy. The Root Mean Square Error, defined as $RMS{E}_n=\sqrt{\frac{1}{n}{\sum }_{i=1}^n{({\widehat{I}}_k-I_k)}^2}$ with n denoting prediction range, $I_k$ and $\widehat{I_k}$ being the actual and modeled active cases, respectively, is as $RMS{E}_{150}=3.283\times {10}^3$. Note that, for the first 30 days, the modeled transmission is higher than the actual one. This can be a result of relatively low daily tests when compared to the rest of the period. In contrast, for $k>120$, the actual transmission is significantly higher than the modeled one. It may be the result of gradually easing the restrictions (from May 16, $k=85$).

The starting point of simulation regarding the Spain epidemic is 27 February, when there were 31 active cases and 0 removed, so $I_0=31$, $R_0=0$, $N=$ 46,757,240 and $S_0$ = 46,757,209. Using the optimization method for $K=80$, we obtain parameters, as follows $\beta =3.76,\gamma =0.2708,{\alpha }_1=0.5374$ and ${\alpha }_2=0.5983$. Figure 5 presents time plots of modeled and actual disease transmission.

The fractional-order SIR model is quite effective in modeling the COVID-19 transmission in Spain, as we can see in Figure 5. However, for $k>85$, the actual transmission is significantly higher than the modeled one. Again, it may result from easing the restrictions (from 2 May, $k=65$). Furthermore, a resurgence is occurring in Spain, and the model is not capable of modeling this process. It is confirmed by RMSE, which, for $n=80$, is equal to $RMS{E}_{80}=7.506\times {10}^3$, but for $n=150$ is as $RMS{E}_{150}=3.083\times {10}^4$. In this case, using time-varying model parameters may improve the fractional-order SIR model performances. However, using more complex models may be necessary in this case.

5. Conclusions

In this paper, the fractional-order generalization of the SIR epidemic model in order to predict the spread of the COVID-19 disease is presented. The model implementation is based on the Grünwald-Letnikov derivative and discretization while using the nabla fractional-order difference. We study the impact of fractional orders of the model derivatives on the dynamic properties of the proposed model. In modeling the COVID-19 transmission, the model’s parameters are estimated while using the genetic algorithm. The model prediction results for the spread of COVID-19 in Italy confirm the usefulness of the introduced methodology. However, modeling the results for Spain show that the model has some limitations. In the case of the second wave of the epidemic, the model’s prediction does not represent the real data. In this case, using time-varying parameters of the fractional-order SIR model or more complex pandemic models is necessary. Nevertheless, the study shows that the presented model is effective and it could be useful in future pandemics studies.