key: cord-0570278-i5wc8nxl
authors: Santos, Juan; Carcione, Jos'e; Savioli, Gabriela; Sciences, Patricia Gauzellino School of Earth; Engineering,; University, Hohai; Nanjing,; China,; Aires, Universidad de Buenos; Ingenier'ia, Facultad de; Petr'oleo, Instituto del Gas y del; Aires, Buenos; Argentina,; Mathematics, Departmet of; University, Purdue; States, United; Oceanography, National Institute of; Geophysics, Applied; OGS,; Trieste,; Italy,; Geof'isicas, Facultad de Ciencias Astron'omicas y; Plata, Universidad Nacional de La
title: An SEIR epidemic model of fractional order to analyze the evolution of the COVID-19 epidemic in Argentina
date: 2021-04-07
journal: nan
DOI: nan
sha: dba71db925e448b6cae0353713574a222b1d01f5
doc_id: 570278
cord_uid: i5wc8nxl

A pandemic caused by a new coronavirus (COVID-19) has spread worldwide, inducing an epidemic still active in Argentina. In this chapter, we present a case study using an SEIR (Susceptible-Exposed-Infected-Recovered) diffusion model of fractional order in time to analyze the evolution of the epidemic in Buenos Aires and neighboring areas (Regi'on Metropolitana de Buenos Aires, (RMBA)) comprising about 15 million inhabitants. In the SEIR model, individuals are divided into four classes, namely, susceptible (S), exposed (E), infected (I) and recovered (R). The SEIR model of fractional order allows for the incorporation of memory, with hereditary properties of the system, being a generalization of the classic SEIR first-order system, where such effects are ignored. Furthermore, the fractional model provides one additional parameter to obtain a better fit of the data. The parameters of the model are calibrated by using as data the number of casualties officially reported. Since infinite solutions honour the data, we show a set of cases with different values of the lockdown parameters, fatality rate, and incubation and infectious periods. The different reproduction ratios R0 and infection fatality rates (IFR) so obtained indicate the results may differ from recent reported values, constituting possible alternative solutions. A comparison with results obtained with the classic SEIR model is also included. The analysis allows us to study how isolation and social distancing measures affect the time evolution of the epidemic.

solutions. A comparison with results obtained with the classic SEIR model is also included. The analysis allows us to study how isolation and social distancing measures affect the time evolution of the epidemic.

We present an SEIR subdiffusion model of fractional order , with 0 < ≤ 1 to analyze the time evolution of the COVID-19 epidemic in Buenos Aires and neighboring areas (Region Metropolitana de Buenos Aires, (RMBA)) with a population of about 15 million inhabitants. RMBA consists of Ciudad Autónoma de Buenos Aires (CABA) plus forty municipalities covering an area of about thirteen thousand square kilometers, where some of these municipalities have rural areas. Thus, RMBA has an average population density of 1100 people/km 2 , but in CABA and many of its neighboring cities this number increases significantly. For example, CABA has a population density of about 14000 people/km 2 . In this work, we consider that RMBA has a uniform population distribution.

The epidemic started officially on March 9th with the number of cases and deaths still increasing at the day of writing (September 22th, 2020. The classical SEIR model ( = 1) has been used by Carcione et al. [1] and Santos et al. [2] to model the COVID-19 epidemic in Italy and Argentina, respectively.

Fractional calculus has been used to define diffusion and wave propagation models in biological and viscoelastic materials [3, 4, 5, 6, 7, 8, 9, 10] . One important property of the fractional-order SEIR model is that incorporates memory and hereditary properties, a behavior exhibited by most biological systems. The use of fractional order derivatives affects the duration of the epidemic, peaks of infected and dead individuals per day and number of number casualties.

Among other authors that have applied fractional calculus to obtain solutions of the SEIR model, we mention Scherer et al. [11] , that used a Grünwald-Letnikov time-discrete procedure, introduced by Ciesielski and Leszczynski [12] (CL method). Besides, Zeb et al. [13] presented an analysis of several numerical methods to solve the SEIR model of fractional order. For general works on fractional calculus including numerical methods, we refer to Podlubny [14] and Li and Zeng [15] .

We first formulate an initial-value problem (IVP) for the classical SEIR model ( = 1) and the SEIR subdiffusion equations of fractional order at the continuous level using the Caputo definition of the fractional derivative [6] . Existence and uniqueness of the solution of this IVP, with positive values, is demonstrated in [13] . The numerical solutions of the continuous IVP are computed by using the timeexplicit algorithm of Gorenflo-Mainardi-Moretti-Paradisi (GMMP method) [16, 17] . The conditional stability of the time-explicit GMMP method (and also of the CL method) was demonstrated by Murillo et al. [19] [see their equation (19) ]. The validation of the GMMP method is performed by comparison of its results against those of the classic SEIR model and those of the fractional Adams-Bashford-Moulton method (ABM method) as defined in [15] .

The parameters of the SEIR model are the birth and death rates, infection and incubation periods, probability of disease transmission per contact, fatality rate and initial number of exposed individuals. These parameters, together with the order of the fractional derivative, are obtained by fitting the number of fatalities officially reported. This is an inverse problem with an infinite number of solutions (local minima) honouring the data, which is solved by using a quasi-Newton technique for nonlinear least squares problem with the formula of Broyden-Fletcher-Goldfarb-Shanno [20] . The numerical simulations give an effective procedure to study the spread of the evolution of virus, analyze the effects of the lockdown measures and predict the peak of infected and dead individuals per day.

For 0 < ≤ 1, the time fractional Caputo derivative ( ( ) is defined as [3, 16, 17, 6] 

where Γ(·) denotes the Euler's Gamma function. Note that the Caputo derivatives of constant functions ( ) = 1 vanish and those of powers of , ( ) = are

The advantage of using the Caputo derivative in Caputo-type IVP's is that the initial conditions are the same as those of the classical ordinary differential equations. For details on the Caputo derivative and its relation with the Riemann-Liouville fractional derivative we refer to [6] . To approximate the time-fractional Caputo derivative, we use a backward Grünwald-Letnikov approximation at time = Δ , = 0, 1, , · · · , with = ( Δ ), Δ being the time step, as follows [16, 17] :

The coefficients

can be obtained in terms of Euler's Gamma function using the recurrence relation

The work by Abdullah et al. [18] presents an analysis of the fractional-order SEIR model formulated in terms of the Caputo derivative and its GMMP time discretization.

The IVP for the classic SEIR system of nonlinear ordinary differential equations is

with initial conditions (0), (0), (0) and (0). A dot above a variable indicates the time derivative, while ( ) is the number of live individuals at time , i.e., = + + + ≤ 0 , 0 being the total initial population. In (4) , is the number of individuals susceptible to be exposed while is the number of exposed individuals, in which the disease is latent; they are infected but not infectious. Individuals in the -class become infected ( ) with a rate and infected become recovered ( ) with a rate . People in the class do not move back to the class since lifelong immunity is assumed. Furthermore, 1/ and 1/ are the infection and incubation periods, respectively, Λ is the birth rate, is the natural per capita death rate, is the average fatality rate, and is the probability of disease transmission per contact. All of these coefficients have units of 1/time. Given the short period of the epidemic in Argentina (6 months at the time of writing), and that the average life expectancy is about 76 years, it is reasonable to assume that Λ = , so that the deaths balance the newborns.

Dead individuals ( ) are computed as ( ) = 0 − ( ), so that the dead people per unit time ( ), can be obtained as [21] :

Next, we reformulate the system (4) into a fractional-order system by using the Caputo derivative in (1):

The reproduction ratio, 0 , indicates the number of cases induced by a single infectious individual. When 0 < 1, the disease dies out; when 0 > 1, an epidemic occurs. Al-Sheikh [22] analyzes the behavior of the SEIR models in terms of 0 . For the SEIR model, 0 is given by [23] 

The infection fatality rate (IFR) is defined as

where this relation holds at all times, not only at the end of the epidemic.

An explicit conditionally stable GMMP algorithm for the fractional order system (6) is formulated as follows [16, 17] :

The results of the GMMP method (9)-(12) will be validated against the solution of the classical SEIR model ( = 1) and the Adams-Bashford-Moulton (ABM) time-explicit scheme as defined in [15] and included in the Appendix.

The results of the GMMP algorithm are cross-checked with those of the ABM solver for the classical SEIR model ( = 1 ) and SEIR models of fractional orders = 0.9 and 0.8. We use the following parameters, given in Chowel et al. [24] and used by Carcione et al. Figures 1-6 show the results of the four classes, S,E,I,R, and the dead and dead per day individuals computed by using the GMMP and ABM algorithms. First, an excellent agreement between the results of the two algorithms is observed for all values of the fractional order derivative . To quantify this agreement, we compute a mean squared relative error between the estimations of both methods. For example, in the computation of infected individuals, the following errors are obtained: 1.512 × 10 −5 for = 1, 9.880 × 10 −6 for = 0.9 and 1.053 × 10 −5 for = 0.8. In particular, the results for = 1 agree with those of Figures 1 and 2 in [1] . Figure 1 shows that decreasing the order of the fractional derivative causes a delay and an increase in the number of susceptible individuals. While for the classical model the number of infectious individuals vanish at long times, this is not the case for the orders = 0.8 and = 0.9 ( Figure 3 ). We run the simulator up to a very long time but the individuals do not vanish, so that the epidemic never ends (in theory). This happens because 0 ≥ 1. We run other examples with different parameters such that 0 < 1 and as expected the number of infectious individuals vanish and the epidemic dies out. For brevity these plots are not shown. The case 0 < 1 is analyzed in Subsection 4.2, when simulating the evolution of the epidemic in the RMBA using fractional derivatives. This value of 0 is associated with the strict lockdown imposed by the government, with a corresponding decrease in the number of infected individuals. Regarding the exposed infected classes (Figures 2-3) , a decrease in causes delays and reduces the amplitude of the peaks of these classes. Furthermore, as decreases the number of casualties increase as seen in Figure 4 while Figure 6 shows a delay and increase of the peak in the number of dead individuals per day. Also, note that Figure 5 shows a delay and decrease in the number of recovered individuals as the order of the fractional derivative decreases.

These simulations consider a single value of , the lockdown parameter. In a realistic case, is a function of time and the procedure is that every time changes, the algorithm has to be fully initialized from the beginning. Changing in the same time loop yields wrong results. This fact has been verified by cross-checking different algorithms and several fractional orders. 

We model the COVID-19 epidemic in the RMBA, with a population 0 = 14839026 individuals according to the 2010 Census (https://www.indec.gob.ar/indec/ web/Nivel4-Tema-2-41-135). The prediction of the time evolution of the epidemic is very difficult due to the uncertainty of the parameters defining the SEIR model. Virus properties such as the infectious and incubation periods ( −1 and −1 ) and life expectancy of an infected individual ( −1 ) lie in certain bounded intervals. Instead, the parameter is time dependent, due to changes according to the lockdown and social-distance measures imposed by the government. Most authors use the infectious individuals to calibrate the model, e.g., González-Parra et al. [25] , who model the AH1N1/09 influenza epidemic in Bogotá, Colombia and in the Nueva Esparta state in Venezuela.

Since the number of asymptomatic, undiagnosed infectious individuals in RMBA is unknown, we choose to calibrate the model with the number of officially reported casualties as the most reliable data, from day 1 (March 9, 2020) to day 198 (September 22th, 2020) (https://www.argentina.gob.ar/coronavirus/ informe-diario). Concerning the parameters, fractional order and initial conditions of the model, we assume = 3.6 × 10 −5 / day, corresponding to a life expectancy of 76 years. Changes in the parameter are associated with different measures of lockdown and social distance imposed by the goverment. Thus, we assume that is a piecewise constant function, where its variations are related to the inflection points observed in the curve of casualties. After the initial time 0 = 1 day, this curve shows two inflection points at times 1 = 31 day and 3 = 50 day. The fractional-order derivative , the values of , , , and the initial exposed individuals (0) are estimated by minimizing the 2 -norm between the simulated and actual casualties, which is an inverse problem with an infinite number of solutions due to the existence of local minima. The estimation is also performed for the classical case = 1. This inverse problem is solved by using a quasi Newton approximation technique for nonlinear least-squares problems, based on the formula of Broyden-Fletcher-Goldfarb-Shanno [20] . Application of this technique to solve inverse problems in reservoir engineering can be found in [26] . Table 1 shows ranges of the fractional derivative , of the parameters , , , and the initial exposed individuals (0) used in the inversion procedure. Table 2 displays the initial values and results of four outputs (Cases) of the fitting procedure. Table 1 Constraints and ranges of the estimation procedure Let us analyze four cases, resulting from the minimization algorithm. We obtained the SEIR parameters, the fractional order and the initial exposed humans values fitting the data. In all the cases, the initial number of infected individuals is assumed to be (0) = 100. Figures 7 and 8 show the dead individuals and dead individuals per day for Case 1. The inflection point at 1 = 30 day, related to a change of 0 from 3.178 to 0.688, shows a decay in the simulated curves, because of the effect of the lockdown. After 1 = 50 day, the curves exhibits a continuous increase in casualties due to the relaxation of the lockdown measures with 0 = 1.725. Figure 9 shows the behavior of all classes, with a a peak of 555 thousand infected individuals at day 188 (September 12th, 2020) while Figure 10 exhibits a death toll of 19000 people after 800 days (May 17th, 2022) and a peak of 234 casualties at day 188.

The parameters of Cases 2 and 3 in Table 2 also fit the data, with graphs similar to those in Figures 7 and 8 . Case 2 estimates peaks of 309 deaths and 285 thousand infected individuals at day 222 (October 16th, 2020). At day 800 (May 17, 2022), there are 34 thousand deaths and 7457 thousand recovered humans. This increase in the number of casualties is due to the higher infection fatality rate IFR and higher reproduction ratios 0 as compared with those of Case 1 (see Table 2 ).

Case 3, which corresponds to the classical SEIR model ( = 1), exhibits a peak of 171 casualties at day 184 (September 8th, 2020) and 607 thousand people infected. The end of the epidemic is consider the day at which the number of infected individuals is smaller than 1, which is day 594 (October 24th, 2021) for this case. At this day, the total number of recovered and dead individuals are 10157 thousand and 18 thousand, respectively, so that the total number of infected people at the end of the epidemic is 10175 thousand individuals. This is the case predicting the smallest number of casualties.

Finally, since the reported number of deceased people could possibly be underestimated due to undeclared cases and delays in the upload of official data, we also consider a case with 30 % more casualties to date (Case 4 in Table 2 ), giving IFR = 0.254 % and values of the parameters similar to those of Case 1. Besides, the peak occurs almost at the same day of Case 1 (day 187: September 11th, 2020) with 592 thousand infected individuals and 296 casualties. This peak of casualties and the death toll of 24400 individuals are approximately 30 % higher than those of Case 1. In the following, we compare the behavior of all classes for the different orders of the fractional derivative used in this analysis, i.e., = 1, 0.919 and 0.812. Figure  11 displays the number of infected individuals, where there is a delay and decrease of the peak values as the order of the fractional derivative decreases. This behavior is consistent with that observed in Figure 3 . Figure 12 shows an increase in the number of casualties by decreasing the order of the fractional derivative, with a 47 % increase between = 1 and = 0.812. Moreover, it can be seen that the curves stabilize at later times as the fractional order decreases. Finally, Figures 13 and 14 exhibit the estimated recovered and susceptible individuals for the three values of . Recovered individuals increase and, consequently, susceptible individuals decrease as the order of fractional derivative increases. The curves exhibit asymptotic values at later times as decreases, and the lower the value of the later individuals recover from the virus infection. Note that the general trends of Figures 11-14 are similar to those of the figures in Subsection 4.1, in spite of the fact that parameters obtained from the adjustment are different for the three cases.

In the four cases described above, we consider that the initial number of infected individuals is (0) = 100. Nevertheless, we tested other values: if (0) belongs to the interval [10, 150] a reasonable adjustment is obtained, with similar values to those shown in Table 2 and a slight delay on the infected individuals peak as (0) decreases. Outside this interval, the fit is poor and the results have no physical meaning. 

We use a fractional SEIR (Susceptible, Exposed, Infected, Recovered) diffusion model to analyze the evolution of the COVID-19 epidemic in Argentina, particularly in the Region Metropolitana de Buenos Aires (RMBA), where a significant number of the population is concentrated.

We solve the SEIR system of fractional order , 0 < < 1 and the classical ( = 1) SEIR model by using a time-explicit Gorenflo-Mainardi-Moretti-Paradisi (GMMP) method. To validate this method, the results were cross-checked with those of the time-explicit fractional Adams-Bashford-Moulton (ABM) method, obtaining an excellent agreement between the two schemes.

Assuming that the birth and death rates are balanced, the parameters that characterize the model are the infection and incubation periods, the probability of disease transmission per contact, the fatality rate and the initial number of exposed individuals. These parameters and the order of the fractional derivative are estimated by fitting the number of casualties officially reported. This inverse problem is solved by using a quasi-Newton technique for non-linear least-squares problem with the Broyden-Fletcher-Goldfarb-Shanno formula.

In all the simulations we used three lockdown parameters (denoted by ), associated with the different measures taken by the government during the evolution of the epidemic. One important conclusion related with this time-dependent parameter is that both the fractional GMMP and ABM algorithms need to be fully initialized from the beginning in order to obtain correct results.

Different cases have been analyzed since the inverse problem has an infinite number of solutions. We observe a similar behavior in all the cases, with a fatality rate IFR varying in the range, [0.175, 0.444]. After the 50th day of lockdown, it is observed a continuous increase in casualties due to the relaxation of the preventive social isolation and community circulation of the virus.

The numerical simulations in RMBA show that when the order of the fractional derivative decreases, i.e., higher subdiffusion of the virus, the duration of the epidemic is extended, and the peak of infected individuals and number of casualties increase. Furthermore, the classical SEIR model yield a smaller number of casualties and infected individuals with associated peaks located at earliest times as compared with those of the fractional-order cases.

The Adams-Bashford-Moulton explicit scheme for the fractional order SEIR equations is formulated as follows [15] Predictor +1 = (( + 1)Δ ) 0 + =0 , +1) 1 ( , , , ) Concerning the error of the numerical scheme ABM, Abdullah et al. [18] give a bound in terms of the time step size Δ . On the other hand, Li and Zeng [15] and Li et al. [27] show that the fractional forward Euler and ABM methods are stable and convergent of order one in Δ .

A simulation of a COVID-19 epidemic based on a deterministic SEIR model

A numerical simulation of the COVID-19 epidemic in Argentina using the SEIR model, submitted to Lat

Linear models of dissipation whose Q is almost frequency independent, Part II

The fundamental solutions for the fractional diffusion-wave equation

Time-domain seismic modeling of constant Q-wave propagation using fractional derivative

Fractional calculus and waves in linear viscoelasticity

Wave simulation in biological media based on the Kelvin-Voigt fractional-derivative stress-strain relation

Hysteresis cycles and fatigue criteria using anelastic models based on fractional derivatives

General fractional calculus evolution equations and renewal processes, Integral Equations of Fractional Order

On the relation between sources and initial conditions for the wave and diffusion equations

The Grünwald-Letnikov method for fractional differential equations

Proc. 15th Conf. on Computer Methods in Mechanics (Wisla, Polonia)

Comparison of Numerical Methods of the SEIR Epidemic Model of Fractional Order

Fractional differential equations

Numerical methods for fractional calculus

Time fractional diffusion: A discrete random walk approach

Convergence of the Grünwald-Letnikov scheme for time-fractional diffusion

Novel analytical and numerical techniques for fractional temporal SEIR measles model

On three explicit difference schemes for fractional diffusion and diffusion-wave equations

Practical optimization

On a new epidemic model with asymptomatic and dead-infective subpopulations with feedback controls useful for Ebola disease

Modeling and analysis of an SEIR epidemic model with a limited resource for treatment

Global dynamics of an SEIRS epidemic model with constant immigration and immunity

SARS outbreak in Ontario, Hong Kong and Singapore: the role of diagnosis and isolation as a control mechanism

A fractional order epidemic model for the simulation of outbreaks of influenza A(H1N1)

Comparison of optimization techniques for automatic history matching

Numerical approximation of nonlinear fractional differential equations with subdiffusion and superdiffusion