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

We provide a SEIR epidemic model for the spread of COVID-19 using the Caputo fractional derivative. The feasibility region of the system and equilibrium points are calculated and the stability of the equilibrium points is investigated. We prove the existence of a unique solution for the model by using fixed point theory. Using the fractional Euler method, we get an approximate solution to the model. To predict the transmission of COVID-19 in Iran and in the world, we provide a numerical simulation based on real data.


Introduction
Coronaviruses are crown viruses that can cause disease in humans and animals. In humans, several coronaviruses are known to cause respiratory illnesses such as common cold and more severe illnesses such as acute Middle East respiratory syndrome (SARS), severe acute respiratory syndrome (SARS), and a recently discovered disease COVID-19. A coronavirus (COVID-19) that was first identified in the Chinese city of Wuhan in 2019 is a new strain that has not been previously identified in humans. Snakes or bats have been suspected as a potential source for the outbreak, though other experts currently consider this to be unlikely. Fever, cough, shortness of breath, and breathing difficulties are initial symptoms of this infection. In the next steps, the infection can cause pneumonia, severe acute respiratory syndrome, kidney failure, and even death.
In this paper, we intend to investigate the spread of COVID-19 disease using the SEIR mathematical model with the Caputo fractional-order derivative. First, we analyze the model mathematically and then, in the numerical section, we present simulations for the release of COVID-19 in Iran and around the world. Also, to evaluate the advantage of using the fraction derivative, we make a comparison between the results of the model with the fractional-and integer-order derivative with the real data to determine which one provides a better approximation in this model. By calculating the model results for different orders of fractional derivative, we investigate the effect of derivative order on the behavior of the resulting functions and equilibrium points and resulting numerical values.
The structure of the paper is as follows. In Sect. 2 some basic definitions and concepts of fractional calculus are recalled. The SEIR model of fractional order for COVID-19 transmission is presented in Sect. 3, and the equilibrium points and the reproduction number are calculated. The stability of the equilibrium points is analyzed in Sect. 4. The existence and uniqueness of solution for the system is proved in Sect. 5. In Sect. 6, a numerical method for solving the model is described and a numerical simulation is presented.

Preliminary results and definitions
In this section, we recall some of the fundamental concepts of fractional differential calculus, which are found in many books and papers. Definition 1 ([38]) For an integrable function g, the Caputo derivative of fractional order ν ∈ (0, 1) is given by Also, the corresponding fractional integral of order ν with Re(ν) > 0 is given by Definition 2 ([24, 39]) For g ∈ H 1 (c, d) and d > c, the Caputo-Fabrizio derivative of fractional order ν ∈ (0, 1) for g is given by where t ≥ 0, M(ν) is a normalization function that depends on ν and M(0) = M(1) = 1. If g / ∈ H 1 (c, d) and 0 < ν < 1, this derivative for g ∈ L 1 (-∞, d) is given by Also, the corresponding CF fractional integral is presented by a, b), and 0 < ν < 1, then the fractional Atangana-Baleanu derivative in the Caputo sense is defined by where B(ν) denotes the normalization function satisfying B(0) = B(1) = 1 and E ν (·) is a one-parameter Mittag-Leffler function. Also, the Atangana-Baleanu fractional integral is given as follows: The Laplace transform is one of the important tools in solving differential equations that are defined below for two kinds of fractional derivative.

Definition 4 ([38])
The Laplace transform of the Caputo fractional differential operator of order ν is given by which can also be obtained in the form

Model formulation
In viral epidemic diseases, mathematical models are very important for predicting the transmission of the virus by considering its behavior in different regions for helping to manage the disease. Different mathematical models such as SIR, SEIR, SIRD, SEIRD, SIRS, SIRC, MSEIR, SEAIHRD, etc. are used to investigate the spread of diseases. According to the information published about COVID-19 by the World Health Organization, there are two types of people with the disease: one group has no symptoms and the other group has symptoms. Both groups transmit the disease to healthy people, and the sick people either recover or die. Of course, during this process, more groups can be considered, such as people who are hospitalized, but because we do not have accurate information about these groups for simulation, we decided to use the simple SEIR model. In this model, we divide people into four groups: susceptible people (S), exposed or asymptomatic infected people (E), symptomatic infected people (I), and recovered people (R) including improved people. The diagram for the dynamics of COVID-19 disease model is shown in Fig. 1. Based on Fig. 1, we consider the SEIR model for COVID-19 as follows: where ω = n × N , N is the total number of individuals and n is the birth rate μ: the death rate of people, β 1 : the transmission rate of infection from E to S, β 2 : the transmission rate of infection from I to S, λ: the transmission rate of people from E to I, δ: the mortality rate due to the disease, τ : the rate of recovery of infected people, with the initial conditions In this section, we moderate the system by substituting the time derivative with the Caputo fractional derivative. The ordinary derivative has an inverse second dimension s -1 and the fractional derivative D ν has a dimension of s -ν . To solve this problem, we use an auxiliary parameter θ that has a second dimension s and is called the cosmic time [41]. By the parameter, from a physical point of view, we will have [θ ν-1C D ν ] = [ d dt ] = s -1 . According to the explanation presented, the COVID-19 fractional model for t > 0 and ν ∈ (0, 1) is given as follows: where the initial conditions are

Nonnegative solution
we show that the closed set ϒ is the feasibility region of system (1).

Lemma 5
The closed set ϒ is positively invariant with respect to fractional system (1).
Proof To obtain the fractional derivative of total population, we add all the relations in system (1). So where Using the Laplace transform and Theorem 7.2 (and Remark 7.1) in [42], we obtain where N(0) is the initial population size. With some calculations, we get Consequently, the closed set ϒ is positively invariant with respect to fractional model (1).

Equilibrium points
To determine the equilibrium points of fractional-order system (1), we solve the following equations: By solving the algebraic equations, we obtain equilibrium points of system (1). The disease-free equilibrium point is obtained as E 0 = ( ω μ , 0, 0, 0, 0, 0). In addition, if R 0 > 1, then system (1) has a positive endemic equilibrium point E 1 = (S * , E * , I * , R * ), so that Also, R 0 is the basic reproduction number and is obtained using the next generation method [43]. To find R 0 , we first consider the system as follows: where .
At E 0 , the Jacobian matrix for F and V is obtained as follows: FV -1 is the next generation matrix for system (1) and the basic reproduction number is obtained from R 0 = ρ(FV -1 ), so we get This basic reproduction number R 0 is an epidemiologic metric used to describe the contagiousness or transmissibility of infectious agents.

Stability of equilibrium points
In this section we investigate the stability of equilibrium points. The Jacobian matrix of system (1) is obtained as follows: So, the Jacobian matrix of system at E 0 is

Theorem 6
The equilibrium point E 0 of system (1) is locally asymptotically stable if R 0 < 1 If R 0 < 1, since all of the parameters are positive, then Also, from R 0 < 1 we have Applying the Routh-Hurwitz criterion, E 0 is locally asymptotically stable. If R 0 > 1, then B < 0, and there is one positive real root for Eq.(ref2), then E 0 will be unstable.
The Jacobian matrix of system (1) at the endemic equilibrium point is The characteristic equation of matrix J(E 1 ) is obtained as follows: The eigenvalues of the characteristic equation are k 1 = -μ, k 2 = -(μ + δ + τ ) and the roots of the equation Since k 1 , k 2 are negative, so E 1 is locally asymptotically stable when two roots of Eq. (2) are negative, so it is enough to have B 1 > 0 and A 1 < 0.

Existence and uniqueness of solution
In this section we show that the system has a unique solution. First, we write system (1) as follows: By taking integral form on both sides of the above equations, we get We show that the kernels Q i , i = 1, 2, 3, 4, satisfy the Lipschitz condition and contraction.

Theorem 7
The kernel Q 1 satisfies the Lipschitz condition and contraction if the following inequality holds: Proof For S and S 1 we have Thus, for Q 1 the Lipschitz condition is obtained, and if 0 ≤ β 1 d 2 + β 2 d 3 + μ < 1, then Q 1 is a contraction.
In the same way, we can prove that Q j , j = 2, 3, 4, satisfy the Lipschitz condition as follows: where S(t) ≤ d 1 , and h 2 = β 1 d 1 + λ + μ, h 3 = τ + μ + δ, h 4 = μ are bounded functions. If for j = 2, 3, 4 we have 0 ≤ h j < 1, then Q j are contractions for j = 2, 3, 4. According to system (3), consider the following recursive forms: with the initial conditions S 0 (t) = S(0), E 0 (t) = E(0), I 0 (t) = I(0), and R 0 (t) = R(0). We take the norm of the first equation in the above system, then with Lipschitz condition (4), we have In a similar way, we obtain Thus, we can write that In the next theorem, we prove the existence of a solution.

Theorem 8 A system of solutions given by the fractional COVID-19 SEIR model (1) exists if there exists t 1 such that
Proof From the recursive technique and Eq. (5) and Eq. (6) we conclude that Thus, the system has a solution and also it is continuous. Now we show that the above functions construct solution for model (3). We assume that By repeating the method, we obtain Taking limit on the recent equation as n approaches ∞, we obtain B 1n (t) → 0. In the same way, we can show that B jn (t) → 0, j = 2, 3, 4. This completes the proof.
To show the uniqueness of the solution, we suppose that the system has another solution such as S 1 (t), E 1 (t), I 1 (t), and R 1 (t), then we have We take norm from this equation It follows from Lipschitz condition (4) that Thus Theorem 9 The solution of COVID-19 SEIR model (1) is unique if the following condition holds: Proof Suppose that condition (7) holds Then S(t) -S 1 (t) = 0. So, we obtain S(t) = S 1 (t). Similarly, we can show the same equality for E, I, R.

Numerical results
Using the fractional Euler method for Caputo derivative, we present approximate solutions for the fractional-order COVID-19 SEIR model [44]. We present simulations to predict the COVID-19 transmission in the world.

Numerical method
We consider system (1) in a compact form as follows: where w = (S, E, I, R) ∈ R 4 + , w 0 = (S 0 , E 0 , I 0 , R 0 ) is the initial vector, and g(t) ∈ R is a continuous vector function satisfying Lipschitz condition Applying a fractional integral operator corresponding to the Caputo derivative to equation (8), we obtain Set h = T-0 N and t n = nh, where t ∈ [0, T] and N is a natural number and n = 0, 1, 2, . . . , N . Let w n be the approximation of w(t) at t = t n . Using the fractional Euler method [44], we get where u n+1,j = (n + 1j) ν -(nj) ν , j = 0, 1, 2, . . . , n.
The stability analysis of the obtained scheme has been proved in Theorem (3.1) in [44].
Thus, the solution of system (1) is written as follows: u n+1,j f 2 (t j , w j ) ,

Case I: The world
To provide a numerical simulation, we must first determine the value of the parameters. The current birth rate for the world in 2020 is 18.077 births per 1000 people, and the death rate is 7.612 per 1000 people [45].  Fig. 2, so that every part is a week. Using this method, we obtain the parameters as follows: β 1 = 2 × 10 -11 , β 2 = 2.2 × 10 -9 , λ = 2.35 × 10 -5 , τ = 0.03. Table 1 compares the absolute and relative errors for I(t) concerning the fractional-and integer-order models with respect to the reported cases of infected people. From this table, you can observe that the Caputo model with ν = 0.97 provides more realistic results than the classic model with integer-order derivatives.
A comparison between the noninteger-order model with ν = 0.97 and the integer-order one with ν = 1 and the real data for the infected cases of the COVID-19 from 4 February to   12 May is also given in Fig. 3. The obtained results show that the answer of the fractionalorder model corresponds well with the real data and together with the results of Table 1 shows the advantage of using the fractional-order derivative instead of the correct order one. In Figs. 4 and 5, we have plotted the results of model (1) for different values of ν. In this simulation, the equilibrium point is E 1 = S * , E * , I * , R * = 2.5 × 10 6 , 7.52 × 10 9 , 2.76 × 10 6 , 2.9 × 10 9 .  Fig. 5 shows that approximately 50 days after 4 February 2020 the number of active cases ceases to increase and becomes stable in 2.76 × 10 6 .
Case II: Iran In the second case, we provide a numerical simulation using real data for the transmission model of COVID-19 in Iran. According to the WHO report, the total population of Iran on 18 February 2020 was N = 83392425, the birth rate for Iran in 2019 was 18.547 births per 1000 people, and the death rate was 4.866 per 1000 people. Thus, for every day, we have ω = n×N 365 = 4237.477 and μ = 0.004866 365 = 0.0000133315. Similarly, we assumed the mortality rate due to the disease in Iran is δ = 3.4 × 10 -2 and θ = 0.99. Since  Fig. 6, so that every part is three days. Using this method, we obtain the parameters as follows: β 1 = 1.1 × 10 -4 , β 2 = 3.3 × 10 -6 , λ = 1.02 × 10 -3 , τ = 0.03. In Figs. 7 and 8, we plotted the results of the system of COVID-19 transmission (1). As you can see in Figs. 7 and 8, the variables have different results in different amounts of ν but exhibit the same behavior. Figure 7 shows that two months after the virus is released almost the entire population is at risk for the disease. Figure 8 shows that the number of people with COVID-19 increases until 200 days and then stabilizes. Also, the forecast is that the number of infected people could rise to 280,000. Also, Fig. 8 shows that the number of people who have recovered or died also increases over time.

Conclusion
In this work, the SEIR epidemic model for the transmission of COVID-19 using the Caputo fractional derivative has been presented. The feasibility region of the system and equilibrium points have been calculated, and the stability of the equilibrium points has been investigated. The existence of a unique solution for the model by using fixed point theory has been proved. Using the fractional Euler method, an approximate answer to the model has been calculated. To predict the transmission of COVID-19 in the world and in Iran, the numerical simulations based on real data have been provided. Also in the numerical section, we have examined the advantage of using the fractional-order derivative instead of the integer-order one, and in Table 1 and Fig. 3, we have compared the results of the model with the fractional-and integer-order derivative and the real data. The results show that the fractional-order model has a better result in this modeling.