Complex mathematical SIR model for spreading of COVID-19 virus with Mittag-Leffler kernel

This paper investigates a new model on coronavirus-19 disease (COVID-19), that is complex fractional SIR epidemic model with a nonstandard nonlinear incidence rate and a recovery, where derivative operator with Mittag-Leffler kernel in the Caputo sense (ABC). The model has two equilibrium points when the basic reproduction number R0>1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$R_{0} > 1$\end{document}; a disease-free equilibrium E0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$E_{0}$\end{document} and a disease endemic equilibrium E1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$E_{1}$\end{document}. The disease-free equilibrium stage is locally and globally asymptotically stable when the basic reproduction number R0<1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$R_{0} <1$\end{document}, we show that the endemic equilibrium state is locally asymptotically stable if R0>1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$R_{0} > 1$\end{document}. We also prove the existence and uniqueness of the solution for the Atangana–Baleanu SIR model by using a fixed-point method. Since the Atangana–Baleanu fractional derivative gives better precise results to the derivative with exponential kernel because of having fractional order, hence, it is a generalized form of the derivative with exponential kernel. The numerical simulations are explored for various values of the fractional order. Finally, the effect of the ABC fractional-order derivative on suspected and infected individuals carefully is examined and compared with the real data.


Introduction
The coronavirus is a new viral pneumonia, deadly and rapidly spreading infection that has put great panic around the globe since the break out late 2019. Within a short span of time, it affected every country all over the world, which urges every nation take severe action to control the wild spread of the virus as the severity of the disease will harm human life badly.
It is well known that infectious diseases are a massive threat for humans and also for national economies. Proper understanding of a disease dynamics is necessary, which plays a major role for controlling and eventually exterminate the infection in a community. Implementation of the correct applicable strategy against the disease transmission is another challenge. Mathematical modeling of the infectious diseases is one of the key tools in order to handle these challenges. A number of disease models have been established in the past, which enables us to investigate and control the spread of infectious diseases in a better way [1][2][3].
Regarding COVID-19, from the initial spread of COVID-19, many of the known mathematical models were used, or newly established by many researchers, For example Giordano et al. in [4] proposed the so called SIDARTHE model and concluded that restrictive social-distancing measures will need to be combined with widespread testing and contact tracing to end the ongoing COVID-19 pandemic. At a same time, Alshammari [5] proposed the SEYNHR compartmental model and applied real data for Saudi Arabia and obtained good predictions on a short term. The role of quarantining and isolation to control the spread of COVID-19 is studied by Memon et al. [6], where they used an extended SIR model. They show that quarantine and isolation are effective measures to control the outbreak.
From real data of the infectious disease, we know that the outbreak of the disease within the country or state for time is generally nonlinear, which tells us to design the system where we can study those dynamic nonlinear phenomena. By this system, we can be able to define the transmission of such a contagious disease, which helps us to interpret the remedial measures to stop or contain the spread of the contagious disease.
Most of the present day studies on mathematical modeling of infectious disease are based on integer-order differential equations (IDEs). However, more recently, many authors state that fractional-order models (FDEs) can describe natural phenomena more precisely than the integer-order differential equations because fractional-order derivatives and integrals describe the memory and hereditary properties of different phenomena [7]. Further, the classical IDEs cannot provide the result for non-integer values. In order to overcome the above limitation of integer-order derivatives, different types of fractionalorder operators were defined, some of the research can be found in the existing literature [8] and applications of these fractional-order operators are noted in [9], also Sene in [10,11] consider the Liu et al. [12] chaotic system with Caputo fractional derivative and show the Lyapunov exponent to characterize the nature of chaos and prove the dissipativity of the considered chaotic system. However, it is known that the classical form of the fractional derivatives with singular kernel may not be suitable to characterize the non-local dynamics in an appropriate manner. In order to overcome this problem, Atangana and Baleanu [13] establish a new model, the so called generalized ABC differential operator involving the Mittag-Leffler (ML) and non-local type kernel, which has an anti-derivative fractional integral operator [14]. This new definition of derivative is shown to be more efficient for the SIR model with generalized incidence rate compared to the other existing fractional models [15] and an exothermic reactions model having a constant heat source in porous media with power problem [16]. More recently, Atangana and Atangana [17] solved the system of ABC fractional partial differential equations to see the side effects of facemasks, which is necessary to protect from COVID-19 and concluded that appropriate facemasks that could help avoid re-inhaling enough CO 2 should be used every time one is in open air even when alone especially in a windy environment. Later, Khan et al. [18] explored the dynamics of COVID-19 with quarantine and isolations with real statistical cases reported in the mainland China, where the model is developed by using a fractal-fractional derivative in the Atangana-Baleanu sense. They indicate that their results are useful in the early erad-ication of the disease in the community. Similarly, Sane in [19] also consider the SIR epidemic model with Mittag-Leffler fractional derivative and delay, he uses a Lyapunov direct method and analyzes the global asymptotic stability of the disease-free equilibrium and the endemic equilibrium. Naik et al. [20] proposed an Atangana-Baleanu-Caputo operator for the transmission of COVID-19 epidemic, where they used real data from Pakistan. They show that fractional operators have many advantages over the existing non-integerorder types.
The studies from [13] to [20] provide us with enough motivation to study and analyze a new fractional SIR model involving the ABC operator with complete memory effects. To the best of our knowledge, this is the first time to study a non-local and non-singular derivative operator for the model of the SIR, which has nonlinear incidence and recovery rates, where we also examine the existence and uniqueness results with the Atangana-Baleanu derivative by using a fixed-point method.
The rest of this paper is organized as follows. Some preliminary results are given in Sect. 2. Section 3 is devoted to the new fractional SIR model and its main characteristics. We first derived the solution of the problem and then provided the existence and uniqueness result results in Sect. 4. In the last section, we demonstrate the numerical solutions for Eq. (10) by using the ABC-derivative in a graphical method. We also discussed the necessity of maintaining enough number of hospital beds for controlling of the infectious diseases.

Preliminary results
We note some background material for the Caputo and Atangana-Baleanu fractional derivatives and Laplace transform are presented [13,[21][22][23], which we use later.

Definition 2.2
Let f ∈ C n be a function, then the fractional Caputo derivative of order α is defined as [13] C where α ∈ (n -1, n) and n ∈ N. Obviously, as α → 1, C 0 D α t (f (t)) tends to the f (t).
where the term B(α) is the normalization function with B(0) = B(1) = 1. The associated fractional integral is given by If α = 0, the initial function is recovered in Eq. (3), we have the ordinary integral for α = 1.

Theorem 2.1 The Laplace transform of fractional Caputo derivative and Atangana-Beleanu-Caputo (ABC) derivatives is given by
and Property 2.1 The inverse Laplace transform of specific functions can be obtained:

Mathematical model and discussion
Let us consider the total population N(t) to consist of three sub-categories: susceptible S(t), infected I(t) and recovered R(t) individuals (SIR), where N(t) = S(t) + I(t) + R(t) as seen in Fig. 1. There are many factors involved in modeling the infectious diseases and these factors substantially affect the dynamical behavior of the models such as incident rate and recovery rate. In the classical SIR model, incident rate β S(t)I(t) N(t) and linear recovery rate μI(t) are used. However, as for COVID-19, many infectious diseases show multiple peaks or periodic oscillations during the outbreak or in progress. Therefore, several nonlinear incidence rates have been suggested since they can produce rich dynamics for the epidemic models [24,25]. Liu et al. [26] proposed the following form of nonlinear saturated incidence rate to include the effect of behavioral changes: where β * , l, h, k * > 0. This model is due to Ruan et al. [27,28]. More recently, a more general form was considered by Rao et al. [29].
On the other hand, in view of the diversity of public resources, it is not enough to study the effect of treatment capacity on a large scale at this time, if we want to further explore the influence of hospital bed number on disease transmission. Then we need to define the nonlinear recovery rate introduced by [30] ( where α 1 is the maximum per capita recovery rate, when health care resources suffice and there are few infectious individuals, α 0 is the minimum per capita recovery rate due to lack of clinical resources, b is the ratio between the number of hospital bed and the population, all the above parameters are non-negative. It has the following properties: Here, we assume individuals move from compartment S to I with transmission rates βIS k+I (l = h = 1 in Eq. (9)), and when an individual is infected, then the individual either recovers at a rate (I, b) = (α 0 + (α 1α 0 ) b b+I ) or dies at the rate of γ I. Then this complex SIR model can be represented by following nonlinear system of equations: where A = n × N, N refer to total number of people and n is the birth rate, μ is a positive constant representing the natural death of the population, γ is a positive constant representing the disease induced death rate, α 0 is the minimum per capita recovery rate due to the function of basic clinical resources and α 1 is the maximum per capita recovery rate due to the sufficient health care resource and few infectious individuals as well as the inherent property of a specific disease.
The model in (11) is not capable of representing the influence of memory effects with respect to time derivative. In order to account for the memory effect, in addition, we modify the fractional operator by an auxiliary parameter σ , having the dimension of second, to ensure that the right-and left-hand sides of the resultant equation possess the same dimension s -1 . Consequently, we can reformulate the system in (9) by using the ABC derivatives where 0 < α < 1. In the first two equations in (12), R(t) does not appear, it is therefore assumed the first two equations are sufficient to handle the full system. We modify the fractional operator by an auxiliary parameter σ , having the dimension of day, to ensure that the right-and left-hand sides of the resultant equation possess the same dimension s -1

Non-negative solution
In this section, we show that the feasibility region of the system (12) is given by the closed set = (S, I, T) : which is positively invariant. (13) is positively invariant with respect to the fractional model (12).

Lemma 3.1 The closed set defined in
Proof Adding all three equations in (12), we have Applying the Laplace transform on both sides, we obtain Using Theorem 2.1 and the properties of the Laplace transform, we got Hence now, applying the inverse Laplace transform and using Property 2.1, we find where m = αμ (B(α)+μ(1-α)) , since the Mittag-Leffler function has an asymptotic behavior, we obtain N(t) ≤Â μ = A μ as t → ∞. Hence, a solution of the fractional model (12) stays in for every t > 0. Consequently, the closed set is a positively invariant set regarding the fractional model (12).

Equilibrium points, basic reproduction number and local stability of equilibrium points
The equilibrium points of the system are the zeros of second side of Eq. (12), obviously, disease-free equilibrium point i.e., I = 0 and S = A/μ or E 0 = ( A μ , 0). Therefore, the model (12) has a threshold parameter R 0 , known as the basic reproduction number, this is defined as the number of secondary infections produced by a single infection in a completely susceptible population. It is not difficult to write the reproduction number for the diseasefree equilibrium point of Eq. (12) as [5] R 0 = βA kμ(α 1 + γ + μ) .

Theorem 3.1 The disease-free equilibrium point E 0 is locally asymptotically stable and if
Proof The jacobian matrix corresponding to system (12) is Then the solution of characteristic equation associated to J( A μ , 0) is given by Hence, if R 0 < 1, the disease-free equilibrium point E 0 is locally asymptotically stable, if R 0 > 1, it is unstable. Theorem 3.2 A sufficient condition for the endemic equilibrium E 1 = (S * , I * ) to be locally asymptotically stable is N 1 < 0 and N 2 > 0, where Proof See the book of Chou and Friedman [31], we here note that Therefore, local asymptotic stability of the endemic equilibrium point is a function of R 0 and the number of available hospital beds (b). (12), if R 0 < 1, i.e. βA < kμ(α 1 + γ + μ), the disease-free equilibrium E 0 = ( A μ , 0) is globally asymptotically stable.

Theorem 3.3 For the system
Proof Consider the Lyapunov function V (t) = I in R 2+ , Calculating the time fractional derivative of the V (t), we obtain Utilizing the system (12), we have ABC 0 The Lyapunov-Lasalle theorem implies that solutions in approach the largest positively in-variant subset of the set V = 0, i.e., the plane I = 0. In this plane, S → A d as t → ∞. Thus all solutions in the plane I = 0 go to the disease-free equilibrium E 0 . Therefore E 0 is globally asymptotically stable.

Existence and uniqueness of the solution
Applying the AB-fractional integral on both sides of (12), we obtain By using Proposition 2.1, we have Using Eq. (3) in the above, we find and where (34) The above system of Eqs. (32)-(33) is the solution of the systems of Eq. (12),

Theorem 4.2 Suppose that the following condition holds:
Then the fractional epidemic model (12) has a unique solution for t ∈ [0, T].
Proof Recursively, the expressions in (35) can be written as As usual taking the difference of successive approximations, we obtain S n = S n (t) -S n-1 (t) It is not difficult to see that S n (t) = n i S i and I n (t) = n i I i , taking the sup norm on both sides, and using the Lipschitz condition on kernel functions and triangle inequality, we obtain Starting from n = 1, we obtain by back substitution from (30) and (31), respectively, and I n ≤ It is clear from our hypothesis that S n → 0 and I n → 0 as n → ∞. We obtain a convergent series as follows: and I n (t) -I 1 (t) = I n (t) -I n-1 (t) + I n (t) -· · · -I 1 (t) ≤ I n (t) -I n-1 (t) + I n-1 (t) -I n-2 (t) + · · · + I 2 (t) - Substituting (55) into (53)-(54) and by integration, we have

Case I: The world
We need to use real data for COVID-19 to determine 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 [16]. The world's population on February was N = 7,610,105,452, so A = n×N 365 = 391,347, μ = 7.612 365×1000 = 2.08547 × 10 -5 , according the WHO report [32] the COVID-19 mortality rate γ = 0.034. The other parameters can be found by a curve fitting technique with real reported data for COVID-19. Figure 1(a) shows the development of susceptible individuals decrease with the time as expected for integer-order derivatives, b = 100 and R 0 = 2.618. Figure 1(b) shows the variations of infected individuals with the time, it increases with time and reaches the plateau region as seen in each country, where we used same parameters as Fig. 1(a). We compare the prediction fractional-order derivatives in Fig. 2. It is well known that an increase in the number of active cases and reaching the plateau region are different for each country and the number of suspected population decreases and reaches the steady case. This is also different for each country. On the other hand, from Fig. 2, it is clear that the active or suspected population not only is a function of time but also a function of the fractional order. This fact proves that the fractional-order derivative can be used the adjust the modeling to the different countries. The other point for fractional-order derivative can also be used as a free parameter to fix the modeling with the data. The effect of the hospital beds for controlling of the diseases are given in Fig. 3(a)-(b). Figure 3(a) shows the difference of the results between b = 100 and b = 10,000. In Fig. 3(a) we show the infected individuals' difference between b = 100 and b = 10,000 could reach 400,000 daily cases, if we assume that only 20% of infected individuals need medical treatment at hospital, then, if we multiply this one by COVID-19 mortality rate 0.034, we obtain 2720 people life can be saved if we provide sufficient hospital beds. Figure 3(b) shows the difference of the suspected individual results of b = 100 and b = 10,000, this is expected because active infected cases for b = 100 always are bigger than active infected cases for b = 10,000. We also note that our numerical results converge to the disease equilibrium point.

Case 2: Saudi Arabia
In the second case, we provide a numerical simulation using real data for the transmis-  Fig. 4(b), we explored the prediction of active cases for different values of fractional order, α = 1 (blue line), 0.9 (red line), 0.7 (black line) and 0.5 (green line), when we compare modeling prediction with the real data (WHO report active cases in Saudi Arabia [32]), we see that model prediction overestimate sometimes lower estimate, but again it seems that α = 0.9 produces the best result. At this point, we note that Comparison of between the results of integer order derivatives (blue line, α = 1), fractional order of derivatives (α = 0.9 (red line), α = 0.7 (black line) and α = 0.5 (green line) for b = 100, k = 5000 and R 0 = 2.618 even our suggested modeling seems simple, but we believe that prediction of this model can be used as a basis for more complicated model predictions. Furthermore, we have one more adjustable parameters in fractional-order derivatives, this gives us enormous flexibility. Further, in Fig. 5, we concentrate on the effect of half saturation constant k over the suspected and infected individual for several values of fractional order of the derivations. We can see as before the effect of the fractional order provides us with enough flexibility and suspected individuals decrease with time and eventually reach the equilibrium point as we expected again the time to reach the equilibrium point to depend on the fractional order of the derivative.

Conclusion
In this study, consideration is given to a complex mathematical fractional SIR model for spreading of the COVID-19 virus, we first show the feasibility region where the solution exists. The equilibrium points (disease-free and endemic equilibrium) of the system are found. The local and global asymptotic stability of the system are discussed. Using the fixed-point theorem, existence and uniqueness of solution is proven. A numerical solution of the system is performed for several values of the parameters involved in the modeling. We showed the importance of maintaining a sufficient number of hospital beds. We also discussed and compared the prediction of our modeling with real word data. Finally, we concluded that using the fractional derivative provides us with one more parameter, which can be used to simulate the real word data.