Fractional SEIRP model for COVID-19 dynamics incorporating social distancing and environment

We propose a Caputo-based fractional compartmental model for the dynamics of the novel COVID-19 epidemic taking into consideration social distancing and the inﬂuence of the environment. Using basic concepts such as continuity and Banach ﬁxed-point theorem, existence and uniqueness of solution to the proposed model was shown. Furthermore, we analyze the stability of the model in the context of Ulam-Hyers and generalized Ulam-Hyers stability criteria. The concept of next-generation matrices was used to compute the basic reproduction number R 0 , a number which determines the spread or otherwise of the disease into the general population. Numerical simulation of the disease dynamics was carried out using the fractional Adam-Bashforth-Moulton method to validate the obtained theoretical results.


Introduction
Epidemiological modeling of infectious disease using differential equations with integer-order to examine and investigate the transmission dynamics of epidemics is widely studied in literature [1,2,3,4,5]. Mathematical modeling of epidemics in literature reveals that nonlinear dynamical equations can give some important insight into the transmission dynamics or dynamical behaviors of disease spread. The recent COVID-19 outbreaks around the world have attracted a lot of interest in the mathematical modeling of this highly contagious disease by constructing realistic nonlinear compartmental mathematical models that are driven by data to better understand the transmission dynamics of the epidemics [6,7,8,9,10,11,12,13,14].
Fractional order differential equations are a very useful and powerful mathematical modeling tool to study biological and engineering systems. This is due to the fact that differential operators in these types of equations or models are related to systems with memory dynamics which naturally exists in most biological and engineering systems, see [15,16,17]. Differential equations that are characterized by Caputo fractional order derivative have been used to study and analysed several infectious diseases transmission dynamics such as HIV/AIDS [18,19,20,21,22,23,24,25,26,27,28,29,30], TB [31,32,33], Malaria [34,35,36], Dengue fever [37,38,39,40], Zika [41,42], Ebola [43,44,45] and Hepatitis B [46,47,48,49,50,51]. In [52], backward bifurcation dynamics is explored in a basic and realistic vaccination epidemic model. Therein, the author considered rigorous qualitative analysis of the formulated mathematical model and also presented numerical simulations to confirm the obtained theoretical results. [53] [53] extended and studied the cholera epidemic mathematical model formulated by the author in [54] to capture Caputo fractional order derivative. The author in [55] used Lyapunov functions that are of Volterra-type and investigated uniform asymptotic stability of some basic epidemic models (SIS, SIR, SIRS) and the well-known Ross vector-borne diseases in Caputo sense. A 3-compartmental deterministic modeling structure (normal weight-overweightobese) is constructed to explore obesity epidemic dynamics in a variable population using Caputo fractional-order nonlinear system [56]. The basic SEIR mathematical model driven by Caputo fractional derivative based on the assumption of variable population dynamics is studied in [57]. The authors gave a detailed qualitative stability analysis for their new and realistic deterministic model.
Due to usefulness and powerful nature of fractional order differential equation models in fractional calculus [58,59,60,61,62], some recent studies in the literature have considered the mathematical modeling of this deadly COVID-19 pandemic using Caputo based nonlinear equations, see [63,64,65]. A Caputo fractional-order deterministic epidemic model for COVID-19 infection is developed and studied in [66]. They used the well-known Banach contraction mapping principle to established the existence and uniqueness of the solution for the mathematical model. In a new mathematical modeling study [67], the authors proposed a SEIQRDP fractional-order deterministic model characterized by Caputo derivative to examine the novel coronavirus epidemic. A SEIPAHRF Caputo fractional-order compartmental model is proposed to analyze COVID-19 transmission dynamics in Wuhan [68]. They numerically solved their proposed model and compared it with the initial 67 days reported data on confirmed infected and death cases in Wuhan city. In recent studies, the authors in [69] proposed a Caputo-type nonlinear COVID-19 epidemiological model to explore the significance of lockdown dynamics in controlling the spread of the infectious disease. They further presented the uniqueness and existence of solutions for the mathematical model under lockdown using Banach and Schauder fixed theorems. In [70], a susceptible-exposed-infected-isolated-recovered compartmental modeling framework with constant total population dynamics is used to construct a new SEIQR Caputo fractional-order COVID-19 mathematical model. Owusu-Mensah and co-authors [71], have also studied a new COVID-19 epidemic model using Caputo derivative in fractional calculus. The well-known and powerful generalized form of Adams-Bashforth-Moulton iterative scheme was applied to numerically solve the formulated nonlinear fractional-order differential equation model.
The goal of the present paper is two folds: to establish both the mathematical and biological well-posedness of the integer-order model proposed in [72] and employ an approximate analytical technique to obtain long-term dynamics of the disease. Subsequently, we modify the model to a dimensionally consistent Caputo fractional-order model which has been touted in the literature to describe more efficiently, memory effect typically associated with the way human cells respond to diseases. The paper is organized as follows. Following the essential preliminaries on fractional calculus in Section 2, we examine in Section 3 both mathematical and epidemiological wellposedness of the integer-order model in [72], and obtain disease dynamics using approximate analytical technique proposed in [73], [74]. Section 4 of the paper is devoted to the formulation, well-posedness, local and global stability analysis of the fractional order model. Therein, using basic concepts such as continuity and Banach fixed-point theorem, existence and uniqueness of solution to the proposed model was shown. Furthermore, we discussed stability analysis of the model in the context of Ulam-Hyers and generalized Ulam-Hyers stability criteria. Also, the concept of next-generation matrices was used to compute the basic reproduction number R 0 , a number that determines the spread or otherwise of the disease into the general population. We conclude the paper with numerical simulation of the fractional order disease model using fractional Adam-Bashforth-Moulton; and discussion of results.

Preliminaries on fractional calculus
In this section, we introduce well-known definitions and Lemmata in fractional calculus that are relevant for the current article. The interested reader can see the monograph [75] and the article [76] for the proofs and further references.
Definition 2.1. [75] The fractional order integral of the function g ∈ L 1 ([0, b], R + ) of order α ∈ R + is defined the sense of Riemann-Liouville as [75] For a function g ∈ C n ([0, b]), the Caputo fractional-order derivative of order v of g is defined by where n = v + 1 and v denotes the smallest integer that is less or equal to v.

Model framework: Integer order model
It is well known that the ravaging Covid-19 virus is transmitted through human-to-human interactions and transmission through the environment. When an infected person sneezes or coughs, the virus is released into the immediate environment which, experts believe, stays viable for as long as five days. We, therefore, consider two interacting populations of humans and pathogens denoted by N (t) and P (t) respectively. At any time t, the total population N (t) is assumed to comprise the susceptible population S(t), the exposed E(t), the asymptomatic infected I A (t), the symptomatic infected I S (t) and the recovered population R(t). Thus, N (t) = S(t) + I A (t) + I S (t) + R(t). These interactions are depicted in Figure 1 below as reported in [72]. In the figure, b denotes the rate of at which human population is born into susceptible class S(t). The terms β 1 SP 1+α 1 P and β 2 S(I A +I S ) 1+α 2 (I A +I S ) represent rate at which the susceptible get infected by pathogens and through interactions with infectious asymptomatic I A (t) and infectious symptomatic I S (t). The denominators in the above terms factor in the adherence to recent experts advice on social distancing and the use of face masks which minimizes contact with infectious individuals, and transmission through the environment. Based on the above description, the authors in [72] proposed the following integer-order differential equation to model the dynamics of the transmission of covid-19 taking into account environment and social distancing, namely Therein, the authors analysed the equilibra and basic reproduction number of the model, and numerical simulation. However, the analysis in article [72] did not examine the essential well-posedness of the proposed model. Here in the present paper, for completeness, we begin our exposition by establishing qualitative properties of the integer order model which ensures that the model is both mathematically and epidemiological well-posed i.e. existence, uniqueness, positivity and boundedness of the solutions to the model. Furthermore, establishing the aforementioned properties allows us the extend the study in [72] to fractional order model, which accommodate the memory effect typically associated with epidemic disease models. In addition, we shall obtain approximate analytical solutions of the model via a multistage method proposed in [73], [74].

Well-posedness
For the fact that (2) models human populations, all the model parameters are assumed nonnegative, see Table 1 for values and description of the model parameters. It then remains to show that unique solution exists for the model and that the state variables remain bounded and non-negative for all time t > 0.
reveal that both the right hand side of (2) namely F and its Jacobian are continuous for t > 0. Thus, F satisfies a Lipschitz condition on R 6 + . The existence and uniqueness of solution for some time interval (0, t f ) follows from Picard-Lindelof Theorem. The state variables S(t), E(t), I A (t), I S (t), R(t), P (t) of (2) with non-negative initial data remains non-negative for all t > 0.
Proof. Let S(0) = S 0 > 0. It follows from the first equation of (2) that and upon integration, one obtains Similar arguments yield and as a consequence, non-negativity of the remaining state variables are obtained.
Thus considering the initial valued problem N (t) = b − µN, N (0) = N 0 and invoking comparison theorem [77], it follows that and consequently lim sup For the pathogen population, using the fact that I A + I S ≤ N ≤ b/µ, the last equation of (2) yields

Approximate analytical solution of the integer order model
Here, we shall adopt the multistage technique proposed in [73], [74] to compute approximate analytical solution to the integer order model (2). For that purpose, let us define U = 1 1+α 1 P and W = 1 1+α 2 (I A +I S ) so that (2) is transformed to an equivalent polynomial system which is now amenable to the proposed technique. By writing S(t) = ∞ n=0 S n t n , it follows from the above polynomial system that the coefficients S n are obtained recursively through In a similar manner, approximate analytical solutions of the remaining compartments are computed. Such series solution often have small interval of convergence. Therefore to improve convergence, we choose a safe step length h, (see [73] and [74]) and compute a piecewise continuous approximate solution which is convergent in the entire integration interval. The dynamics of the disease progression for each sub-population are displayed in Figures 2a-2d and Figures  3a-3b below, confirming the results in [72].

Fractional order model
Fractional derivative are generally believed to model disease epidemic more realistically because of their capability to capture memory effect often associated with human body response to diseases. Here, we propose a fractional order variant of (2) given by In the above, operator D α denotes Caputo fractional derivative of order 0 < α ≤ 1. Noting that all the model parameters except α 1 , α 2 and δ have dimensions 1/t α , we have raised these parameters to power of α for dimensional consistency emphasized by [38].

Qualitative properties of solution
In this section, we examine the mathematical and biological well-posedness of the fractional order model. In essence, we prove that solution of the fractional model is bounded and remains positive as long as a positive initial condition is given. Furthermore, we prove the existence and uniqueness of the solution to the modified model. Let X(t) = (S, E, I A , I S , R, P ) T and K(t, X(t)) = (φ i ) T , i = 1, 2, . . . 6 where . Then the dynamical system (5) can be written as In the above, the condition X(0) ≥ 0 is to be interpreted component-wise. Problem (6), which is equivalent to fractional differential equation (5), in turn has integral representation Next, we shall analyse model (5)

Positivity and boundedness of solution
For the fractional order model to be biologically well-posed, its solution is expected to be positive and bounded at all times. These properties are established in the sequel.
Finally for boundedness, proceeding as in the integer order case (see the proof of Theorem

Existence of unique and uniformly stable solution
Here, we establish existence, uniqueness and uniform stability of solutions to (4) through (5).
The following preliminary result is in order. Proof. From the first component of K, we observe that However, Similarly we obtain, Altogether, we have In a similar manner, one obtains For the remaining components of K, it holds For uniqueness, we shall use Banach contraction mapping principle on operator P defined in (8) above. For that purpose, we show that P is both a self map and a contraction. Firstly, by definition, sup t∈[0,b] ||K(t, 0)|| = b v . Let us now define κ > ||X 0 ||+Ωb v 1−ΩL K and a closed convex set B κ = {X ∈ E : ||X|| E ≤ κ}. Thus for self map property, it suffices to show that P B κ ⊆ B κ . So, let X ∈ B κ , then It then follows that operator P X ⊆ B κ and P is indeed a self-map. It remains at this point to prove that P is a contraction. Let X andX ∈ E satisfy (6). Then again using the result of Thus if ΩL K < 1 then P is a contraction mapping, and consequently by the Banach contraction mapping principle, P has a unique fixed point on [0, b] which is solution of (5). Uniformly Lyapunov stability of solution follows from [80, Theorem 3.2].

Equilibra and Basic reproduction number of the fractional order model
By setting the left hand side of (5) to zero, one obtains the equilbra points. Disease-free equilibrium points are those where I A = I S = 0 and P = 0. These immediately imply E = 0 and R = 0. Thus, we obtain b α −µ α S * = 0 or S * = b α /µ α . Consequently, the disease-free equilibrium point of the fractional-order model is (S * , E * , I * A , I * S , R * , P * ) = ( b α µ α , 0, 0, 0, 0, 0). The dynamic and stability of a disease model is usually governed by the basic reproduction number R 0 which simply put, measure the average number of secondary infections resulting from introducing one infected individual into an entirely susceptible population. In other words, it determines the spread or otherwise of the infection into the entire population. We follow the recipes described in [81] to compute the basic reproduction number for the fractional order model. The approach is based on next-generation matrices.
The infected subsystem of (5) comprises the second, third, fourth and sixth equations. Therefore, letting x = (E, I A , I S , P ) T , then the infected subsystem can be written in the form The matrices F and V represents the rate of production of new infections and transition of new infections respectively.
Linearizing the infected subsystem about the disease-free equilibrium point, we obtain the corresponding rate of infection and infection transition rate near the equilibrium, respectively, as The basic reproduction number is the spectral radius of matrix K = F V −1 , the total production of new infections over the course of outbreak. Using symbolic computation Maple17 software, As suggested by [81], the above matrix can be further reduced to a 2 × 2 matrix through a certain transformation matrix E, while still preserving the dominant eigenvalue. It is immediate Let us now define the reduced matrix K S by Consequently, it follows that ρ(K) = ρ(K S ) = 1 2 trace(K S ) + trace(K S ) 2 − 4 det(K S ) . Hence, The following theorem, which is immediate from [82,Theorem 2], elucidate the importance of the basic reproduction number R 0 derived above.
The above theorem 4.4 implies that the fractional order epidemic model (6) is controllable through appropriate control measures (via the parameters α 1 and α 2 ) such as minimizing contact with infectious individuals as well as shedding of the virus pathogen by infected individuals to the environment.

Global stability
We establish the global stability of the fractional model (6)  For clarity of the discussion that follows, let us introduce the inequality given by We say a functionX ∈ E is a solution of (10) if and only if there exists h ∈ E satisfying It is important to observe that by invoking (7) and property ii. above, simple simplification yields the fact that any functionX ∈ E satisfying (10) also satisfies the integral inequality Definition 4.5. The fractional order model (6) (and equivalently (5)) is Ulam-Hyers stable if there exists C K > 0 such that for every > 0 and for each solutionX ∈ E satisfying (10), there exists a solution X ∈ E of (6) with Definition 4.6. The fractional order model (6) (and equivalently (5)) is said to be generalized Ulam-Hyers stable if there exists a continuous function φ K : R + → R + with φ K (0) = 0, such that, for each solutionX ∈ E of (10), there exists a solution X ∈ E of (6) such that We now present our result on stability of the fractional order model.  (5)) is Ulam-Hyers stable and consequently generalized Ulam-Hyers stable.

Numerical simulation and discussion
This section provides some illustrative numerical simulations to explain the dynamical behavior of the Caputo fractional-order deterministic nonlinear COVID-19 mathematical model. The numerical solution of a nonlinear mathematical model using the appropriate iterative scheme is very important in mathematical modeling. For this purpose, we have used the Adams-type predictor-corrector iterative scheme constructed in [83,84] to solve fractional-order differential equations. The method relies on equivalent integral formulation of our fractional model (5) written in the form (6). Consider a uniform discretization of [0, b] given by t n = nh, n = 0, 1, 2, . . . , N where 0 < h = b/n denote the step size. Now, given any approximation X h (t i ) ≈ X(t i ), we obtain the next approximation X h (t i+1 ) using the Adams-type predictor-corrector iterative scheme as follows; Predictor : Corrector : For varying values of the derivative order α, the dynamics of the fractional order model are presented in Figures 4 below. We have used the model parameters in Table 1 as reported in [72] and initial data S(0) = 93000, E(0) = 1000, I A (0) = 50, I S (0) = 50, R(0) = 0, P (0) = 500. For α = 1, coinciding with the integer-order model, the obtained dynamics is consistent with those reported in Section 3 and [72]. In Figures 4, effect of the fractional order α on disease progression are demonstrated. Furthermore, the basic reproduction number of the diseasefree equilibrium point Σ = b α µ α , 0, 0, 0, 0, 0 = (1.9861, 0, 0, 0, 0, 0, ) for α = 0.5 was computed as R 0 = 0.87831 < 1, showing the fulfillment of the necessary and sufficient conditions for local asymptotic stability of the disease-free equilibrium of Theorem 4.4. We also found out that for fractional orders considered namely α ∈ {1, 0.9, 0.8, 0.7}, and with new infection rates α 1 = α 2 = 0.1, the corresponding computed R 0 ∈ [0.90850, 0.82590] showing that the COVID-19 pandemic is controllable and will effectively die out as long as there is compliance with social distancing/lockdown regulations, and if infectious and infected individuals are appropriately quarantined, thereby preventing contamination of the environment through virus shedding.
Recall that the constants α 1 and α 2 represent the proportion of interaction with the infectious environment and infectious individuals respectively. Consider a scenario of a contaminated environment α 1 = 0.05 where there is higher chances of contracting COVID-19 through the environment than through infectious individuals I A or I S , (α 1 = 0.1). In this situation, a rapid decline in the susceptible population is noticeable in figure 6a during the first ten days. This is expected as more people get exposed through virus-laden environment, explaining the rapid rise in the exposed population and the infected population within the same time period, figure 6b, 6c and figure 6d. In comparison with other scenarios, namely α 2 = 0.05, α 1 = 0.1 where there is higher risk of contracting COVID-19 virus through contacts with infectious individuals than through the environment; and α 2 = 0.1, α 1 = 0.1 depicting generally low infection rate from both sources, the number of exposed, asymptomatic infectious and symptomatic infectious populations peaked in the former.

Conclusion
We have proposed a fractional order epidemic model for COVID-19 infection dynamics based on the integer-order model of [72] incorporating social distancing measures and interaction with environment. Qualitative analysis of both integer order and fractional order models were carried out. In particular, we established the well-posedness of the fractional model via the Banach fixed point theorem. Using the concept of next-generation matrices for computing basic reproduction number R 0 , local asymptotic stability of the fractional model was proved. Furthermore, we showed the global stability of the model in the sense of Ulam-Hyer stability criteria. Numerical simulations confirmed the established properties of the proposed model, and more importantly elucidated the need for compliance with regulations on basic preventive measures such as social