Dynamic Systems with Fractional Derivatives Applied to Interagent Populations Problems

ABSTRACT

The present work has as main objective to investigate the use of fractional calculus in the modeling of epidemic outbreaks of interacting populations. In particular, we propose a generalization of the SIR model with fractional derivatives to describe the dynamics of the COVID-19 epidemic outbreak in two cities with interacting populations. In special, we consider the dynamics of COVID-19 in the municipalities of Pelotas and Rio Grande, which are neighboring cities and are relatively geographically isolated from the rest of the state of Rio Grande do Sul.

Keywords:
population dynamics; fractional differential equations; SIRD model; coronavirus (COVID-19)

1 INTRODUCTION

The main motivation for the study of differential equations is that even simple equations describe important real-world problems, such as, for example, the growth of a population, the proliferation of diseases, the mass-spring systems, among others. ¹¹. In this context, fractional calculus, also known as non-integer calculus, plays an increasingly important role. Since the beginning of the theory of differential and integral calculus, mathematicians like Leibniz, Riemann, and Liouville developed their ideas for derivatives and integrals of non-integer orders. Although as old as usual calculus, it was only in the last decades that fractional calculus’s importance in modeling real-world phenomena, displaying non-local characteristics, was verified. In this context, several phenomena that can not be well explained by models based on traditional calculus, have been successfully described by models containing fractional derivatives and integrals ¹⁰.

Among the recent areas of application of fractional calculus, we can highlight the modeling of epidemics, in which models are proposed to describes the epidemic dynamic, to help to outline policies to control these diseases and to eradicate them as soon as possible ⁵. Based on this premise, due to the current scenario, where the world faces the pandemic caused by the coronavirus (COVID-19), an infectious disease caused by the severe acute respiratory syndrome virus 2 (SARS CoV-2) belonging to the SARS class. Although the SARS coronavirus (SARS-CoV) outbreak in 2003 was fatal for 9% of infected individuals, it spread to only 26 countries and resulted in about 8, 000 cases. However, the new coronavirus outbreak has become an unprecedented threat to worldwide health ⁹. Due to these indices, it was decided in this work to study a structured model of a compartmental epidemic of the SIRD type for interacting populations. The SIRD model is a generalization of the well-known SIR (susceptible-infected-recovered) compartmental model that includes a compartment for dead individuals, denoted by D.

In this context, one of the objectives of the present study is to analyze real data referent to the number of confirmed cases and deaths, related to COVID-19, in the municipalities of Pelotas and Rio Grande. Because they are neighboring cities, there are interactions between populations (flow of people moving between the two cities), we propose an interagent fractional SIRD model. The solution of the model will be obtained through Python, by employing a semi-analytical method described in Section 4. The data relating to the populations were extracted from the website https://www.riogrande.rs.gov.br/(Rio Grande City Hall) and http://pelotas.com.br/(Pelotas City Hall).

Due to the importance of understanding the proliferation of infectious diseases in population dynamics, there is a significant amount of work in the literature that addresses this subject. In particular, we point out some studies related to the dynamics of COVID-19 in Brazil and a proposed approach for different interacting populations.

The articles ^{3), (7} focus on the role of fractional calculus to model the dynamic of COVID-19. Both works consider a fractional SIR model. In ³ the role of the memory effect, introduced by fractional derivatives due to its non-local characteristics, is investigated in COVID-19 dynamics for several countries. A detailed study of the application of a fractional SIR model to the spread of COVID-19 in Brazil was carried out in ⁷. It is important to stress that one important difference between our present work and the references ^{3), (7} is that, while we consider two distinct populations in interaction, in ^{3), (7} all people in the country belongs to a single population.

In ², real data was analyzed regarding the number of confirmed cases and deaths caused by COVID-19 in the Brazilian states of Paraná (PR), Rio Grande do Sul (RS), and Santa Catarina (SC). The objective of the study was to understand some aspects of the spread of the disease in the southern region of Brazil by obtaining an effective Reproduction Number R _e for each state, a SIRD model was used for this. It was shown that this model, despite its simplicity, accurately describes data from the past and makes it possible to make reliable projections for the trends of the epidemic in each location. For the results, it was shown that, until June 6, 2020, SC was the only state in the South region with $R_e<1$, which indicates that the number of new cases tends to decrease in the hypothesis that this scenario is maintained. On the other hand, PR and RS presented $R_e>1$, so that the growth in the number of infected people should continue in the following weeks if new measures are not taken.

Finally, in the context of interacting populations, the work ⁶ proposes a structured SEIR (susceptible-exposed-infected-recovered) model for the interaction of n different populations to describe the spread of pandemic diseases such as COVID-19. For the authors, the proposed model has the flexibility to include geographically separate communities, as well as taking into account other groups (for example, age groups) and their interactions. Different assumptions about the dynamics of the proposed model were shown, which lead to a curve similar to a plateau of the infected total population, reflecting data collected in large countries such as Brazil. Such observations point to the following conjecture: ”The spread of COVID-19 disease from the capitals to the interior of Brazil may be responsible for the appearance of a plateau on the infected curve”.

2 FRACTIONAL CALCULUS

The non-integer order calculus, known as fractional calculus, began on September 30 of 1695, in a letter from Leibniz to L’Hôpital. In addition to Leibniz and L’Hôpital, other brilliant mathematicians, such as Euler, Lagrange, Laplace, Fourier, Abel, Heaviside, Liouville, among others, studied the subject leading to the first definitions of fractional derivatives and integrals ^{4), (8}. Despite being as old as conventional calculus, it was only in the last three decades that Fractional Calculation attracted more attention due to its applications in diverse areas. This is because the realistic modeling of several complex physical phenomena depends not only on instantaneous time but also on the history of the previous time. Lately, a large number of studies have been developed on the application of fractional differential equations, in various areas of applied sciences, such as fluid mechanics, viscoelasticity, biology, physics, engineering, etc. ¹.

In this work we will use the fractional derivative of Caputo, because among the several existing definitions of fractional derivatives, the Caputo definition is the more popular fractional derivative definition among physics and engineers because of its algebraic properties and, most important, because fractional differential equations with Caputo derivatives require usual initial (or boundary) conditions ⁴.

The Caputo derivative definition is directly related to the analytical continuation of the Cauchy formula for repeated integration, known as Riemann-Liouville fractional integral ⁴:

Definition 2.1 (Riemann-Liouville integral). For $\alpha >0(\alpha \in \mathrm{ℝ})$ , the operator ${}_a{J}_x^{\alpha }$ defined on $L_1\left(\right[a,b\left]\right)$ by

${}_a{J}_x^{\alpha }f\left(x\right)=\frac{1}{\Gamma \left(\alpha \right)}{\int }_a^x\frac{f\left(t\right)}{(x-t{)}^{1-\alpha }}dt,$ (2.1)

is called the left Riemann-Liouville fractional integral of order α.

Remark 2.1.It is important to notice that when α is a positive integer n, the fractional Riemann- Liouville integral (2.1) reduces to a usual integer order n-fold integration⁴.

The Caputo fractional derivative of a function f(x) is defined by taking a fractional Riemann- Liouville integral of order α of the $n=\left[\alpha \right]+1$-th integer-order derivative of f(x) ⁴:

Definition 2.2 (Caputo derivative). The left fractional derivative of Caputo of order $\alpha >0(\alpha \in \mathrm{ℝ})$ is defined by

${{}_a{}^{C}D}_x^{\alpha }f\left(x\right)=\frac{1}{\Gamma (n-\alpha )}{\int }_a^x\frac{f^{\left(n\right)}\left(t\right)}{(x-t{)}^{1+\alpha -n}}dt(n=[\alpha ]+1)$ (2.2)

where $f^{\left(n\right)}\left(t\right)=\frac{d^{n}f\left(t\right)}{d{n}^n}$ is a usual derivative of integer order n.

Remark 2.2. In particular, for $\alpha =1$ the Caputo derivative is reduced to a usual first-order derivative.

Finally, the Caputo fractional derivative (2.2) and the Riemann-Liouville fractional integral (2.1) satisfies the generalized Fundamental Theorem of Calculus given by:

Theorem 2.1 (Fundamental theorem of Caputo Calculus).For$0<\alpha \le 1$, and$f\in A{C}^1[a,b]$or$f\in C^1[a,b]$, the following equality holds:

${{}_{a}J}_x^{\alpha }{{}_a{}^{C}D}_x^{\alpha }f\left(x\right)=f\left(x\right)-f\left(a\right).$ (2.3)

The proof of Theorem 2.1 can be found, for example, in ⁴.

3 THE FRACTIONAL SIRD MODEL FOR TWO DISTINCT POPULATIONS

In 1927, Kermack and McKendrik published the article entitled ”A Contribution to the Mathematical Theory of Epidemics”, in which they introduced the epidemiological model known as SIR. This model describes the spread of infectious diseases in a population divided into subgroups of susceptible, infected, and recovered individuals, with each of these groups having their dynamics described by a Differential Equation ². Although the Kermack and McKendrick model describe the essence of the dynamics of an epidemic, it includes only the main elements of the process. More complex models include other categories of individuals in the population, such as the SIRD model (Susceptible-Infected-Recovered-Dead), which divides the removed people into Recovered and Dead. With that, other parameters are included, in this case, the recovery and mortality rates.

Let us consider the case of two distinct interacting populations (populations in neighboring cities). Let N _i be the number of individuals in the population $i(i=1,2)$ and $N_T=N_1+N_2$ the total integrated population. Let S _i , I _i , R _i and D _i be the fractions, in relation to the N _i , of the i population that are, susceptible, infected, recovered and deaths, respectively, in time t. In the model, the time evolution given by the following dynamic system to i = 1, 2 is considered.

$$\begin{gathered} {{}_0{}^{C}D}_t^{\alpha }{S}_i=-\sum _{j=1}^2{\beta }_{ij}{S}_i{I}_j+{\mu }_i(1-S_i) \\ {{}_0{}^{C}D}_t^{\alpha }{I}_i=\sum _{j=1}^2{\beta }_{ij}{S}_i{I}_j-({\gamma }_i+{\kappa }_i+{\mu }_i)I_i \\ {{}_0{}^{C}D}_t^{\alpha }{R}_i={\gamma }_i{I}_i-{\mu }_i{R}_i \\ {{}_0{}^{C}D}_t^{\alpha }{D}_i={\kappa }_i{I}_i+{\mu }_i{D}_i, \end{gathered}$$ (3.1)

where $0<\alpha \le 1$ is the order of the derivative, β _ij is the disease fractional transmission rate (proportional to the average contact rate within the population and within the population), κ _i is the fractional death rate of the population infected by the disease, γ _i is the fractional inverse of the average infectious period, to i = 1, 2. In addition, the µ _i fractional mortality rates are assumed to be equal to birth rates, so that the total N _i of the population is constant during illness. Therefore, one has to $S_i\left(t\right)+I_i\left(t\right)+R_i\left(t\right)+D_i\left(t\right)=1$. Furthermore, concerning the dimension of the parameters in the model it is important to notice that, since S _i , I _i , R _i and D _i are dimensionless, and the fractional derivative operator ${{}_0{}^{C}D}_t^{\alpha }$ has dimension of time^−α (as the integer order derivative $\frac{d}{dt}$ has dimension of time⁻¹), the fractional rates β _ij , κ _i , γ _i , and µ _i should have dimension of time^−α .

Regarding the existence of the equilibrium points (time independent solution), we have that when t → ∞ the model solution (3.1) tends to the equilibrium point. In the situation where this equilibrium point is not zero, the disease becomes endemic. To find the equilibrium point, we just need to analyze under what conditions the variables S _i and I _i become constant over time. We will call S _i ^∗ and I _i ^∗ these constant values. To find (S _i ^∗ , I _i ^∗), we make ${}_0{}^{C}S_i=0$ and ${}_0{}^{C}I_i=0$ in the first two equations of (3.1). Let us also consider, for simplicity, the symmetrical case. In this case we will have N ₁ = N ₂ = N, β ₁₁ = β ₂₂ = β , β ₁₂ = β ₂₁ = $\hat{\beta }$, γ ₁ = γ ₂ = γ, κ ₁ = κ ₂ = κ, µ ₁ = µ ₂ = µ, S ₁ = S ₂ = S ^∗ and I ₁ = I ₂ = I ^∗. So, from the equation (3.1) we have:

$$\begin{gathered} -\beta S^{*}{I}^{*}-\hat{\beta }{S}^{*}{I}^{*}+\mu (1-S^{*})=0 \\ \beta S^{*}{I}^{*}+\hat{\beta }{S}^{*}{I}^{*}-(\gamma +\kappa +\mu )I^{*}=0. \end{gathered}$$ (3.2)

Isolating I ^∗ in the second equation of (3.2), we get:

$I^{*}\left(\right(\beta +\hat{\beta })S^{*}-(\gamma +\kappa +\mu \left)\right)=0.$

Considering $I\ne 0$, we get the value of S ^∗:

$S^{*}=\frac{\gamma +\kappa +\mu }{\beta +\hat{\beta }}.$ (3.3)

Substituting S ^∗ in the first equation of (3.2), we obtain that I ^∗ is given by:

$I^{*}=\frac{\mu }{\gamma +\kappa +\mu }-\frac{\mu }{\beta +\hat{\beta }}.$ (3.4)

Therefore, considering the particular case with no deaths (κ = 0 e D _i = 0), the fixed point of the symmetric model is:

$P_1(S^{*},I^{*},R^{*},D^{*})=\left(\frac{\gamma +\mu }{\beta +\hat{\beta }},\frac{\mu }{\gamma +\mu }-\frac{\mu }{\beta +\hat{\beta }},1-\frac{\gamma }{\beta +\hat{\beta }}-\frac{\mu }{\gamma +\mu },0\right),$ (3.5)

that corresponds to an endemic equilibrium point, and where we use that $R^{*}=1-S^{*}-I^{*}$ when D ^∗ = 0.

4 RESULTS AND DISCUSSION

In this section, analyzes and comparisons are being carried out around our proposed model (3.1) for integer-order (α = 1) and fractional ($0<\alpha <1$) derivatives. First, we consider some properties of the fractional model for the particular case without the interaction between cities. Next, the fractional model with interaction between two populations is analyzed. For comparison with real data, we used data obtained through the website https://www.riogrande.rs.gov.br/ (Rio Grande City Hall) e http://pelotas.com.br/ (Pelotas City Hall).

In order to obtain a numerical solution for the generalized SIRD model (3.1), first we integrate both sides of (3.1) with a fractional integral ${{}_{0}J}_t^{\alpha }$ of order α. From the Fundamental theorem of Caputo Calculus, (2.3), we have:

$$\begin{gathered} S_i\left(t\right)=S_i\left(0\right)-\sum _{j=1}^2{\beta }_{ij}{{}_{0}J}_t^{\alpha }{S}_i{I}_j+{\mu }_i{{}_{0}J}_t^{\alpha }(1-S_i) \\ {I}_i\left(t\right)=I_i\left(0\right)+\sum _{j=1}^2{\beta }_{ij}{{}_{0}J}_t^{\alpha }{S}_i{I}_j-({\gamma }_i+{\kappa }_i+{\mu }_i){{}_{0}J}_t^{\alpha }{I}_i \\ {R}_i\left(t\right)=R_i\left(0\right)+{\gamma }_i{{}_{0}J}_t^{\alpha }{I}_i-{\mu }_i{{}_{0}J}_t^{\alpha }{R}_i \\ {D}_i\left(t\right)=D_i\left(0\right)+{\kappa }_i{{}_{0}J}_t^{\alpha }{I}_i+{\mu }_i{{}_{0}J}_t^{\alpha }{D}_i. \end{gathered}$$ (4.1)

The second and last step towards solving Problem 3.1 consists in defining discretized functions from S _i , I _i , R _i and D _i . For a positive integer L, let $t_n=t_0+nh(n=0,1,...,L)$, where $h=\frac{t_L-t_0}{L}$ and t ₀ = 0. Then, let $X_i^{\left(n\right)}=X_i\left(t_n\right)$ (where X states for S, I, R or D). The discretized functions are defined by:

$X_i^L\left(t\right)=X_i^{\left(n\right)}\text{if}{t}_n\le t<t_{n+1}.$ (4.2)

Since X _i (t) is a differentiable function, we have $\underset{L\to \infty }{\lim }{X}_i^L\left(t\right)=X_i\left(t\right)$. Consequently, the discretized functions $X_i^L\left(t\right)$ are an approximation for X _i (t) when L is large. We can now find the values $X_i^{\left(n\right)}$ from the following recurrence relations:

$$\begin{gathered} S_i^{\left(n\right)}=S_i^{\left(0\right)}-\sum _{j=1}^2{\beta }_{ij}\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^{t_n}\frac{S_i^L{I}_j^L}{(t_n-t{)}^{1-\alpha }}dt+{\mu }_i\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^{t_n}\frac{1-S_i^L}{(t_n-t{)}^{1-\alpha }}dt \\ {I}_i^{\left(n\right)}=I_i^{\left(0\right)}+\sum _{j=1}^2{\beta }_{ij}\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^{t_n}\frac{S_i^L{I}_j^L}{(t_n-t{)}^{1-\alpha }}dt-({\gamma }_i+{\kappa }_i+{\mu }_i)\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^{t_n}\frac{I_i^L}{(t_n-t{)}^{1-\alpha }}dt \\ {R}_i^{\left(n\right)}=R_i^{\left(0\right)}+{\gamma }_i\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^{t_n}\frac{I_i^L}{(t_n-t{)}^{1-\alpha }}dt-{\mu }_i\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^{t_n}\frac{R_i^L}{(t_n-t{)}^{1-\alpha }}dt \\ {D}_i^{\left(n\right)}=D_i^{\left(0\right)}+{\kappa }_i\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^{t_n}\frac{I_i^L}{(t_n-t{)}^{1-\alpha }}dt+{\mu }_i\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^{t_n}\frac{D_i^L}{(t_n-t{)}^{1-\alpha }}dt. \end{gathered}$$ (4.3)

Finally, by computing the integrals in (4.3) we obtain an approximated solution for (3.1):

$$\begin{gathered} S_i^{\left(n\right)}=S_i^{\left(0\right)}-\sum _{j=1}^2\frac{{\beta }_{ij}}{\Gamma (\alpha +1)}\sum _{m=1}^n{S}_i^{(n-m)}{I}_j^{(n-m)}{\Delta }_{n,m}^{\alpha }+\frac{{\mu }_i}{\Gamma (\alpha +1)}\sum _{m=1}^n(1-S_i^{(n-m)}){\Delta }_{n,m}^{\alpha } \\ {I}_i^{\left(n\right)}=I_i^{\left(0\right)}+\sum _{j=1}^2\frac{{\beta }_{ij}}{\Gamma (\alpha +1)}\sum _{m=1}^n{S}_i^{(n-m)}{I}_j^{(n-m)}{\Delta }_{n,m}^{\alpha }-\frac{{\gamma }_i+{\kappa }_i+{\mu }_i}{\Gamma (\alpha +1)}\sum _{m=1}^n{I}_i^{(n-m)}{\Delta }_{n,m}^{\alpha } \\ {R}_i^{\left(n\right)}=R_i^{\left(0\right)}+\frac{{\gamma }_i}{\Gamma (\alpha +1)}\sum _{m=1}^n{I}_i^{(n-m)}{\Delta }_{n,m}^{\alpha }-\frac{{\mu }_i}{\Gamma (\alpha +1)}\sum _{m=1}^n{R}_i^{(n-m)}{\Delta }_{n,m}^{\alpha } \\ {D}_i^{\left(n\right)}=D_i^{\left(0\right)}+\frac{{\kappa }_i}{\Gamma (\alpha +1)}\sum _{m=1}^n{I}_i^{(n-m)}{\Delta }_{n,m}^{\alpha }+\frac{{\mu }_i}{\Gamma (\alpha +1)}\sum _{m=1}^n{D}_i^{(n-m)}{\Delta }_{n,m}^{\alpha }. \end{gathered}$$ (4.4)

where

${\Delta }_{n,m}^{\alpha }=(t_n-t_{n-m}{)}^{\alpha }-(t_n-t_{n-m+1}{)}^{\alpha }.$ (4.5)

For all the numerical results presented in the work we consider h = 0.5, given a time step of half a day.

4.1 A single isolated population: model without the interaction between cities

Before considering the interaction between two cities, a preliminary analysis was made for only a single city (N ₂ = 0). In this case, we consider real data of the infected and recovered population, that is, the total number of cases registered in Rio Grande, in a period (t) of 284 days. For the study, we considered the initial time (t = 0), and the final time (t = 284). The choice of the initial time was made for the registration of the first case of infected people in the municipality, at 23 in March of 2020.

In Figure 1 are displayed the real data of the population of Rio Grande and the solution for the SIRD model with integer-order derivative (α = 1). The parameters used in the model in Figure 1 were γ ₁ = 2/3 (we assume an average time of 1.5 days (or 36 hours of contact with other people) for the infected person to be diagnosed and isolated from the rest of the population) and, for simplicity, κ ₁ = µ ₁ = 0.0. As an initial condition, we considered that there were 5 individuals initially infected (we assume that when the first case was diagnosed, this patient would have already infected other people). With these parameters, we find that the value β ₁₁ = 0.683 is the one that best describes the real cases.

Figure 1:
Comparison between the real data of the population of Rio Grande (Source: https://www.riogrande.rs.gov.br/ (Rio Grande City Hall)) and the numerical solution of the model without interactions and derivatives of order α = 1.

Finally, it can be seen that the curve does not describe well the cases reported by the official data, since the real dynamics show an approximately linear growth. According to ⁶ this linear growth may be an indication of the spread of the disease to neighboring cities and between neighborhoods of the same city. This fact leads to anomalous growth dynamics, which motivates the investigation of models of interacting populations and fractional models.

One of the main effects on the dynamics introduced by the fractional derivatives can be seen in Figure 2. It can be seen in this figure that for the integer-order derivative, there is a very narrow peak of infected people (I partition) for the epidemic. In the case of fractional derivatives, as the α order of the derivative decreases, there is an elongation in the period of the epidemic. This extension can describe, for example, the spread of the disease between neighborhoods and neighboring cities. Diffusion produces an increase in the period of the epidemic because, in each neighborhood and city, the peak of infection occurs at a different time ⁶, as we will see in our model when considering the interaction between two cities. For the results presented in Figure 2 we used N ₁ = 211965 (data from the population of Rio Grande obtained through the website https://www.riogrande.rs.gov.br/pagina/284-anos-riogrande-mantem-forte-potencial-economico-e-turistico/ (Rio Grande City Hall)), β ₁₁ = 1, γ ₁ = 0.5, κ ₁ = 0 (for simplicity, we do not consider death cases) and µ ₁ = 0.001. In addition, we considered that at t = 0 there were 2 infected individuals.

Figure 2:
Numerical solution for the infected population of the model without the interaction between cities, and derivatives of order α = 1, α = 0.75 and α = 0.5.

We present at Figure 3 the behavior of the dead (partition D) for the model without interaction. For this figure we use N ₁ = 211965, β ₁₁ = 1, γ ₁ = 0.5, κ ₁ = 0.01 and µ ₁ = 0.001. In addition, we consider again that at t = 0 there were 2 infected individuals. We observed again that the introduction of fractional derivatives in the system’s dynamics results in an extension of the time necessary to reach the peak in the number of deaths. In addition, the amplitude of this peak decreases with the order of the fractional derivative.

Figure 3:
Numerical solution for the dead population of the model without the interaction between cities, and derivatives of order α = 1, α = 0.75 and α = 0.5.

Finally, another important effect introduced by fractional derivatives in the dynamics is presented in Figure 4, where we consider for simplicity κ ₁ = 0. The equilibrium point (S _eq , I _eq , R _eq , 0) (where $S_{eq}=\frac{{\mu }_1+{\gamma }_1}{{\beta }_1},I_{eq}=\frac{{\mu }_1}{{\mu }_1+{\gamma }_1}-\frac{{\mu }_1}{{\beta }_{11}}$ and $R_{eq}=1-S_{eq}-I_{eq}$) do model is not affected by the order of the derivative, however, the oscillatory behavior of the model with integer derivative (α = 1) disappears in the model with fractional derivative (α = 0.75). For t → ∞ the model reaches equilibrium, and the disease becomes endemic to ${\mu }_1\ne 0$. For the results shown in Figure 4 we used N ₁ = 211965, β ₁₁ = 1, γ ₁ = 0.5, κ ₁ = 0 and µ ₁ = 0.02 (in this case, we use a large value for µ to facilitate graphical visualization of the balance for those infected). In addition, we considered that at t = 0 there were 2 infected individuals.

Figure 4:
Numerical solution of the model without interactions for derivatives order α = 1 and α = 0.75. In graphs (a), (b) and (c), S _eq , I _eq and R _eq represent the equilibrium values, respectively.

4.2 Interacting populations

Based on the geographical spread of the disease, it was decided to analyze the interaction between neighboring cities. For simplicity, let us consider the case where the two cities have the same total population, that is, N ₁ = N ₂. We consider that there were initially 2 infected in the city 1 and 0 infected in the city 2. Note that in Figure 5, for the graph (a), there is a difference in the time in which the epidemic occurs in each city since initially we had infected individuals only in the first city. The contamination of the population of the second city occurs initially due to the interaction with the first city. This difference in time leads to a phenomenon of widening of the peak and duration of the epidemic in the total population $N_1+N_2$, as pointed out in [6]. We can see this widening in the graph (b), where we present the total sum of infected $(I_1+I_2)$ of these interacting populations. The parameters used in this simulation were α = 1, γ ₁ = γ ₂ = 0.5, κ ₁ = κ ₂ = µ ₁ = µ ₂ = 0.0, β ₁₁ = β ₂₂ = 1 and β ₁₂ = β ₂₁ = 0.0005 (we assume that 1 for every 2000 inhabitant transits between the two cities daily).

Figure 5:
Graph of the fraction of the infected population between two interacting cities (a) for integer order derivative α = 1 (on the left); Graph of total infected population per day (b) (on the right).

The same widening phenomenon is also present in the model with fractional derivatives, as can be seen in Figure 6 for α = 0.75.

Figure 6:
Graph of the fraction of the infected population between two interacting cities (a) or fractional derivative of order α = 0.75 (on the left); Graph of total infected population per day (b) (on the right).

In Figures 7 and 8 we show the asymptotic behavior of models with integer derivatives (α = 1) and fractional derivatives (α = 0.75). For simplicity, we consider the symmetric case β ₁₁ = β ₂₂ = 1, β ₁₂ = β ₂₁ = 0.0005, γ ₁₂ = γ ₂₁ = 0.5, µ ₁ = µ ₂ = 0.2 and κ ₁ = κ ₂ = 0. In this symmetrical case, in equilibrium we have: lim_t→∞ S ₁ = lim_t→∞ S ₂ = S _eq = 0.51974013, lim_t→∞ I ₁ = lim_t→∞ I ₂ = I _eq = 0.01847153, lim_t→∞ R ₁ = lim_t→∞ R ₂ = R _eq = 0.46178833 and D ₁ = D ₂ = D _eq = 0. We observed again that one of the effects of the fractional derivative on the dynamics is the disappearance of the oscillatory behavior of the model variables before reaching the asymptotic equilibrium.

Figure 7:
Numerical solution of the model with interactions between two populations for derivatives of the order α = 1. In graphs (a), (b) and (c), S _eq , I _eq and R _eq represent the equilibrium values, respectively.

Figure 8:
Numerical solution of the model with interactions between two populations for derivatives of the order α = 0.75. In graphs (a), (b) and (c), S _eq , I _eq and R _eq represent the equilibrium values, respectively.

4.3 A comparison with real data

We now present a comparison with real data between the model of interacting populations with integer derivatives and the model with fractional derivatives. We used data obtained from the websites of https://www.riogrande.rs.gov.br/ (Rio Grande City Hall) and the http://pelotas.com.br/ (Pelotas City Hall). We consider the total number of cases recorded in Rio Grande and Pelotas over a period of 284 days. For the study, we considered the initial time at t = 0, and the final time at t = 284. The choice of the initial time coincides with the registration of the first case of infected people in the municipality of Rio Grande, being set at 23 in March of 2020. The last day corresponds to the 31 December 2020 record. Data for the first 2 months of 2021 were not included in the analysis as the municipality of Rio Grande did not release the data on a daily basis in January.

For a more realistic description of the real data of an epidemic, we should consider that the parameters defining the model (β _i , β _ij , γ _i , κ _i , µ _i ) are functions of time, since the value of these parameters depends, for example, on public policies and public awareness. But for the purpose of this work, which is to compare the results obtained by the model of interacting populations with fractional derivatives with the model with integer derivatives (α = 1), we will consider all these parameters as constants. In addition, we considered, as an initial condition, that at t = 0 there were 5 infected in Rio Grande and 0 in Pelotas. The beginning of the epidemic in Pelotas is due to the interaction with Rio Grande. We consider as population 1 the city of Rio Grande, where the first case was diagnosed, and as population 2 the city of Pelotas. In this case, we have N ₁ = 211695 and N ₂ = 343651 (data from the population of Pelotas obtained through the website https://www.pelotas.com.br/cidade/dados-gerais (Pelotas City Hall)).

In Figures 9 and 10 we present the comparison for the accumulated data of infected I _rg and $I_p(I+R+D)$ and of dead D _rg and D _p for the cities of Rio Grande and Pelotas, respectively. For simplicity, we set γ ₁ = γ ₂ = 2/3 and µ ₁ = µ ₂ = 0. In addition, we also consider that β ₁₂ = β ₁₁ /500 and β ₂₁ = β ₂₂ /500 (on average, 1 for every 500 inhabitants of the city has contact with the inhabitants of the neighboring city). In Figure 9 we present the comparison with the number of infected and killed for the model with integer derivatives (α = 1). The parameters that best describe the official data in this case are: β ₁ = 0.7, β ₂ = 0.706, κ ₁ = 0.017 and κ ₂ = 0.014. In Figure 10 we present the results of the interacting model with fractional derivatives of the order (α = 0.8). The parameters that best describe the actual data in this case are: β ₁ = 0.723, β ₂ = 0.738, κ ₁ = 0.016 and κ ₂ = 0.013. Although we believe that initially in Pelotas there was 0 infected, due to the interaction between the two cities, the disease spreads quickly to Pelotas.

Figure 9:
Infected population and number of deaths for the interacting cities of Rio Grande and Pelotas (Source: https://www.riogrande.rs.gov.br/ (Rio Grande City Hall) and http://pelotas.com.br/ (Pelotas City Hall)). In (a) the cumulative number of infected for the model with integer derivative (α = 1), and in (b) the number of deaths.

Figure 10:
Infected population and number of deaths for the interacting cities of Rio Grande and Pelotas (Source: https://www.riogrande.rs.gov.br/ (Rio Grande City Hall) and http://pelotas.com.br/ (Pelotas City Hall)). In (a) the cumulative number of infected for the model with fractional derivative (α = 0.8), and in (b) the number of deaths.

We observed that both cases (integer derivative and fractional derivative) describe the dynamics of the epidemic qualitatively well in the first 284 days in Rio Grande and Pelotas. Despite this, the fractional derivative model behaves quite differently from the integer derivative model when considering time intervals greater than 284 days. In Figure 11 we show the forecast of the total cumulative number of people who will be infected and killed. Figure 11 (a) shows the result for the number of infected for the integer (α = 1) and fractional (α = 0.8), derivative, and Figure 11 (b) corresponds to the model result for the number of deaths. The fractional derivative model predicts a higher number of infected and dead compared to the integer derivative model. As the epidemic in the cities of Rio Grande and Pelotas is still far from reaching equilibrium, we cannot determine which of the two models (fractional and integer) will best describe the entire dynamics of the epidemic.

Figure 11:
Numerical solution of the model with interactions between two populations for derivatives of order α = 1 (graph (a)) and α = 0, 8 (graph (b)).

5 CONCLUSION

The objective of the present work is to investigate the use of fractional calculus in the modeling of epidemics of interacting populations. We propose a SIRD model for two interacting cities with fractional derivatives. We made some preliminary analyzes, and as a continuation of the work, we will analyze the data, over a longer time period, regarding the cumulative number of confirmed cases and deaths caused by COVID-19 in the municipalities of Pelotas and Rio Grande. From the preliminary analysis, it was possible to observe that the classic SIRD model, without interaction between cities and without fractional derivatives, does not describe the official data of the municipality of Rio Grande as well, since the dynamics show a linear growth. This result motivates a more detailed investigation of the fractional interacting model that we propose.

Acknowledgments

C. A. Tavares was financed in part by Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), Brazil, Finance Code 001. M. J. Lazo thanks Fundação de Amparo à Pesquisa do Estado do Rio Grande do Sul (FAPERGS), Brazil, Grant No. 20/2551-00002721, and Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), Brazil, Grant No. 310386/2020-9 and Grant No. 425333/2018-3.

REFERENCES

• ¹
  O.A. Arqub & A. El-Ajou. Solution of the fractional epidemic model by homotopy analysis method. Journal of King Saud University Science (2013).
• ²
  R.M. da Silva. Using the SIRD model to characterize the COVID-19 spreading in the states of Paraná, Rio Grande do Sul, and Santa Catarina. Scientific Electronic Library Online, (2020).
• ³
  L.C. de Barros, M.M. Lopes, F.S. Pedro, E. Esmi, J.P.C. dos Santos & D.E. Sánchez. The memory effect on fractional calculus: an application in the spread of COVID-19. Computational and Applied Mathematics, 40 (2021), 72.
• ⁴
  K. Diethelm. “The Analysis of Fractional Differential Equations: An Application-Oriented Exposition Using Differential Operators of Caputo Type”. Springer-Verlag, Berlin Heidelberg (2010).
• ⁵
  A.V. Farias. Um estudo da modelagem epidemiológica SIR usando conceitos de derivadas de ordem inteira e fracionária. Dissertação de Graduação, Universidade Federal do Rio Grande (2017).
• ⁶
  M.J. Lazo & A. de Cezaro. Why can we observe a plateau even in an out of control epidemic outbreak? A SEIR model with the interaction of n distinct populations for Covid-19 in Brazil. Trends in Applied and Computational Mathematics, (2020).
• ⁷
  N.Z. Monteiro & S.R. Mazroche. Fractional Derivatives Applied to Epidemiology. Trends in Computational and Applied Mathematics , 22 (2021), 157-177.
• ⁸
  D.S. Oliveira. Introdução ao Cálculo Fracionário. Dissertação de Graduação, Universidade Federal do Rio Grande (2012).
• ⁹
  I. Owusu-Mensah, L. Akinyemi, B. Oduro & O.S. Iyiola. A fractional order approach to modeling and simulations of the novel COVID-19. Advances in Difference Equations, (2020).
• ¹⁰
  J. Sabatier, O.P. Agrawal & J.A.T. Machado. “Advances in Fractional Calculus: Theoretical Developments and Applications in Physics and Engineering”. Springer, Netherlands (2007).
• ¹¹
  N. Varalta. “Das Transformadas Integrais ao Cálculo Fracionário Aplicado á Equação Logística”. Master’s thesis, Universidade Estadual Paulista Júlio de Mesquita Filho (2014).

Publication Dates

• Publication in this collection
  27 June 2022
• Date of issue
  Apr-Jun 2022

History

• Received
  15 Mar 2021
• Accepted
  27 Oct 2021