Global analysis of a time fractional order spatio-temporal SIR model

We deal in this paper with a diffusive SIR epidemic model described by reaction–diffusion equations involving a fractional derivative. The existence and uniqueness of the solution are shown, next to the boundedness of the solution. Further, it has been shown that the global behavior of the solution is governed by the value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0$$\end{document}R0, which is known in epidemiology by the basic reproduction number. Indeed, using the Lyapunov direct method it has been proved that the disease will extinct for \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}R0<1 for any value of the diffusion constants. For \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}R0>1, the disease will persist and the unique positive equilibrium is globally stable. Some numerical illustrations have been used to confirm our theoretical results.

In this context, most SIR models have been traditionally investigated in an uniform distribution of populations which are generally formulated only by ordinary differential equations. This fact shows the possibility that the disease can outbreak over a spatial region. In reality the infected individuals have the greatest effect on spatially nearest susceptible persons. The outbreak of infectious diseases is influenced by the spatial movement of populations. The great development in transportation networks is among the main contributing factors in the growth of people's movement around the world. For these reasons, many recent researches have been devoted to the study of reaction-diffusion models (particularly the existence, uniqueness, positivity and stability of the equilibria). They have as goal predicting the evolution of diseases in relation to time and space simultaneously 4 .
Recently, fractional derivatives have several applications in many fields such as mechanics 5 , control theory 6 , bioengineering 7 and viscoelasticity 8 . We point out that derivative order (fractional) can be any positive real choice in order to best correspond to the available data 9 . Consequently, the systems of non-integer order differential equations or partial differential equations give a more realistic behavior 10 . Fractional-order-derivatives are used widely in epidemiology to describe disease evolution and, in most cases, are considered to be more precise than the classical derivative 11 . For example, the spread of the virus is generally discontinuous, so that they are not well described by systems of ordinary differential equations. Then fractional systems naturally deal with such a property of discontinuity 12 . In addition, different models have used fractional derivatives to better predict the outbreak of diseases with sufficient data, among these models we find SIR epidemic model with fractional derivative with Mittag-Leffler kernel 13 , hybrid variable-order fractional coronavirus (2019-nCov) 14 , a hybrid stochastic fractional order Coronavirus (2019-nCov) 15 . In 16 a survey for novel fractional biological models and the numerical methods used to study these models. In 11 the authors used the real data from the Florida Department of Health in the period from September 2011 to July 2014, they concluded that the absolute error between the solutions obtained statistically and that of the fractional model decreases more than those obtained by the model of integer derivative.
In the literature, several definitions of fractional derivatives have been used in different works 17 . Among the most popular non-integer derivatives is that of Riemann-Liouville. It is not often adequate for modeling physical systems because it does not keep the nullity of the derivative constant and the initial conditions of the Cauchy problem are given by fractional derivatives. Caputo presents another alternative preserving the derivative of the constant is null and the initial conditions remain expressed as in the classical case by derivatives of integer order 18 .
The novelties of this work can be summarized in the following two points: Firstly, we investigate the global behavior of more real extension's of a basic SIR model with memory effects measured by Caputo's fractional derivative in time of order 0 < α ≤ 1 , such that when α = 1 we obtain a classical model (without memory). More precisely, memory measured by the nonlocal operator of the fractional derivative highlights the other possibilities not included in the model formulation as fear from infection and the movement in the space, and closing stores, reducing the mobility of persons, and others, which makes fractional systems more realistic to describe the real life situations. Secondly, we have constructed Lyapunov functions to show the global stability of the equilibrium points in a more general framework where the proposed system takes into account the spatial behavior of populations and memory effect. Through this research, we show that our system is well posed in the sense that we prove the global existence, uniqueness and boundedness of the solution. By constructing suitable Lyapunov functions, the disease will be eradicated for R 0 ≤ 1 , and persist for R 0 > 1.
This research is structured as follows. In "Preliminaries" section, we remember some basic results for fractional calculus. The proposed model formulation is given in "Model formulation" section. Next, we show in "Existence and uniqueness of solutions" section the existence and uniqueness of a bounded solution. "Existence of equilibria and local stability" section is devoted for calculating all possible equilibrium states. The global behavior of the solution is the subject of interest in "Global stability" section. Numerical simulations of the considered model in agreement with theoretical results are illustrated in "Graphical representation" section.

Preliminaries
We now recall definitions of the Mittag-Leffler function and Caputo fractional time derivative. First, the Mittag-Leffler function, E α (z) , is defined as the family of entire functions of z given as whenever the series converges 19 , where Ŵ(·) is Gamma function .
(t − s) α−1 is a power law function.

Definition 2.2 ( 17 )
We consider α > 0 , and letting n ∈ N such that n − 1 < α ≤ n . The fractional derivative in sense of Caputo for α for a function f ∈ C n ([0, +∞), R) is with D = d dt for n = 1 . In particular, for 0 < α < 1 , we get For more details about the definition of fractional derivative in sense Caputo , we refer to 17 .

Model formulation
Let a bounded set in R with smooth boundary ∂� , and [0, T] is a finite interval. The classic SIR epidemic model 20 governed by reaction-diffusion equations takes the following form: , α > 0, z ∈ C, www.nature.com/scientificreports/ The total population N is divided into three compartments of the pathological state where S(t, x), I(t, x) and R(t, x) are respectively the densities of the susceptible population, infected population, removed population at time t and the spatial location x.
The positive constants , r and µ are respectively the entering flux into S-class, the recovery rate and the natural death rate. Persons in S-class acquire infection after a direct contact with person in I-class, with the rate βSI , where β is the transmission coefficient per unit of time. Positive constants 1 , 2 and 3 denote the diffusion coefficients for the susceptible, infected and removed individuals, respectively. represents the usual Laplacian operator.
In the actual world, and during an epidemic, there are numerous components that varies and can influence the outbreak of a disease and cannot be included in the model formulation, as the fear generated by the population from infection, weather (specific change in the weather help or reduce the spread of disease). This phenomena can be modeled by replacing the ordinary differential derivative by a fractional one, as it is used in understanding many real world phenomena, we cite for instance the research 21 . For the model (3.1), the state at any time t does not depend on the previous history. It is a markovian process (memory depends on time and corresponds to a Dirac function δ(t, s) ). However, the evolution and control of epidemic processes in human societies cannot be envisaged without any memory effect. When an idea is propagated within a human population, the experience or knowledge of individuals on this idea should influence their responses. To incorporate long-term memory effects in the classical SIR model (3.1), we convert it to an equivalent system of integral equations. We change the function δ(t, s) by the power law function M(t, s) which shows a slow decay so that the state of the system at early times also contributes to the evolution, afterward applying the fractional Caputo derivative 22 . One of our goals is to study the effect of vaccination on the basic reproduction number. We therefore introduced the term vaccination u into model (3.1). It is assumed that vaccination converts susceptible individuals into the removed class and confers immunity on them. Motivated by the above discussion, we introduce the time fractional derivative to the diffusive SIR model in the following manner: We consider that the model (3.2) is self-contained and there is a dynamic across the boundary but there is no emigration. Then the no-flux homogeneous Neumann boundary conditions are For epidemiological aspect, we consider that the following initial conditions of the three classes are positive The constant u refers to the vaccination rate. It is presumed that the vaccination transforms the persons in S-class to the removed class after acquiring immunity. We denote by ν the outward unit normal vector on the boundary ∂� and by ∂ ∂ν = ν.∇ the normal derivative.

Existence and uniqueness of solutions
Letting X = C(�, R) and X 3 be Banach spaces endowed by the uniform norms where J = (S, I, R) and J 0 = (S 0 , I 0 , R 0 ).  Noting that the two first equations does not depend on the class R(t, x) and then are uncoupled with the last equation of the system (3.2). Hence our attention will concentrated on the analysis of the following reduced system: for (t, x) ∈ [0, ∞) × �.

Existence of equilibria and local stability
The principal goal of this section is to determine the equilibria for (4.3). A crucial idea in epidemiology is the existence of significant threshold values quantifying and measuring an outbreak spread in a population. The given value is the basic reproduction number 26 . It is understood as the average number of newly cases of infection, generated by an person in I-class during infectious period, in a compartment entirely composed of susceptible individuals. From the definition of R 0 , we conclude the following results. . We then obtain S * = µ+r β and I * = µ+u β (R 0 − 1) . Hence, if R 0 > 1 , there exist a unique positive solution which is E * .
Next, we study the local stability of the disease-free equilibrium E f and the endemic equilibrium E * . The Jacobian matrix of system (4.3) at any equilibrium Ē = (S,Ī) is given by We recall that a sufficient condition for the local stability of Ē is where ξ i are the eigenvalues of JĒ (see 27 ). First, we establish the local stability of E f .

Theorem 5.2
If R 0 < 1 , then the disease-free equilibrium E f is locally asymptotically stable.
, www.nature.com/scientificreports/ We now establish the local stability of E * .

Theorem 5.3
If R 0 > 1 , then the endemic equilibrium E * is locally asymptotically stable.
Proof At equilibrium E * , the characteristic equation for the corresponding linearised system of model (4.3) is ξ 2 + a 1 ξ + a 2 = 0 , where and If R 0 > 1 then a 1 > 0 and a 2 > 0 . From 28 , we have the desired result.

Global stability
Our goal now is to study the global behavior for E f and E * using Lyapunov's functions. Firstly we give the definition of Lyapunov function. Consider the following fractional differential equation: with the initial condition: where D α t is the fractional derivative in Caputo sense of order α ∈ (0, 1] , the state variable is a positive vector of m elements u 1 , . . . , u m , and f : R m −→ R m is a function of class C 1 . Definition 6.1 ( 29 ) Let u * be an equilibrium point of system (6.1) such that �(u * ) a neighborhood of u * . Let V be a differentiable function defined on �(u * ) and with real value. We say that V is a Lyapunov function in u * , if it satisfies the following two properties: V (u * ) = 0 and V (u) > 0 in �(u * ) for all u = u * . Theorem 6.2 (LaSalle principle 30 ). Let u * be an equilibrium point of the system (6.1) and let V be a positive function of class C 1 defined in the neighborhood �(u * ) of u * . Then u * is asymptotically stable if: Moreover, if V (u) → ∞ , when �u� → ∞ , then u * is globally asymptotically stable.
Secondly, to prove the global stability of DFE, we need to use the following auxiliary lemma for the purpose of the application of Lyapunov function in the case of the fractional order systems: Lemma 6.3 ( 31 ). We put y(t) ∈ R * + be a continuous and derivable function. For all α ∈ (0, 1) and for t ≥ t 0 where is a positive function defined by �(y) = − ln(y) + y − 1 , y > 0.

Theorem 6.4 E f is globally asymptotically stable for
Proof Introducing the Lyapunov function: Calculating the fractional derivative of V in Caputo's sense, we have , www.nature.com/scientificreports/ Applying Green's formula, we get . Then the following two cases arise: • If R 0 < 1 , then I = 0.
• If R 0 = 1 , using the first eq. of (4.3) together with S = S f , we get then βS f I(t, x) = 0 . Thus, we obtain I = 0. Hence, the largest invariant set of (S, I) ∈ R 2 + : C 0 D α t V (t) = 0 is the singleton {E f } . Using LaSalle's invariance principle 30 , we conclude that E f is globally asymptotically stable.
Similarly, we shall show global stability of E * which is resumed in the following theorem Theorem 6.5 E * is globally asymptotically stable whenever exists.

Graphical representation
In this section, we present some graphical illustrations confirming our theoretical findings. The system (3.2)-(3.4) is numerically integrated by using the forward finite difference approximations to discretize the time-fractional derivative 32 and the centered finite difference schemes to approach the Laplacian's operator in one-dimensional space, then we can take = [0, L] . This method gives an accurate of order 2 − α in time and order 2 in space 32 . Let δ = T N and x = L n be the length of each time step and the space step respectively, for some large N and n, t l = lδ for l = 0, . . . , N and x i = i x for i = 0, . . . , n. We have and Then, we obtain the following scheme: with Note that unlike the usual derivative, the fractional derivatives are not local operators, i.e. for example to calculate S l i the number of susceptible at time l, we need all of its information up to the initial instant, and that comes from the summation term l j=1 (j which represents the memory effect. We also notice when α = 1 we obtain the discretization of classical model without memory. www.nature.com/scientificreports/ Next, we study the case without vaccination. Let L = 10 and T = 50 . We simulate system (3.2)-(3.4) with the following set of parameters: µ = 0.8 , = 0.9 , β = 0.1 , r = 0.02 , u = 0.0 , 1 = 2 = 3 = 0.2 , and the initial conditions S(0, x) = 1.0 , I(0, x) = 2.0 and R(0, x) = 3.0 . As a result we approximate the solutions of (3.2)-(3.4) for α = 1 , α = 0.8 and α = 0.6 that displayed respectively in Figs. 1, 2 and 3. We also calculate R 0 = 0.1372 . Hence, system (3.2)-(3.4) has a unique equilibrium E f = (1.12, 0, 0) . Using Theorem 6.4, E f is globally stable. In Fig. 4, we have fixed the space variable x to show the effect of the order α along the dynamics of the solution. We notice that all the solutions are globally asymptotically stable for different values of α not just for α = 1 . We also notice that the solution for α = 1 quickly converges to the equilibrium point E f . Since fractional derivatives describe reality well, we can say that the epidemic takes a longer duration to be stable. This is very important in terms of economics and the study of control strategies. Now, we consider µ = 0.2 and letting the same previous set of parameter. Then, R 0 = 2.0455 . From Theorem 6.5, E * is globally asymptotically stable. Figures 5, 6 and 7 illustrate this result for different values of α , which means biologically that the infection persists but it is under control. For easy comparison see Fig. 8.
Finally, let's study the effect of the vaccine on the basic reproduction rate R 0 . We are only interested in the endemic case for α = 0.8 . We note that the basic reproduction rate decreases when the vaccine rate increases (for example u = 0.8 we have R 0 = 0.4091 (see Fig. 9)). Consequently the number of infected individuals is  www.nature.com/scientificreports/ also decreasing. Besides, the number of removed individuals increases at the expense of susceptible people (see Fig. 10). It reflects the importance of the vaccine to eradicate the disease.

Conclusion
We dealt in this paper with the qualitative behavior of the solutions of a reaction-diffusion system under the influence of the fractional derivative α . Firstly, we investigated the global behavior of more real extension's of a basic SIR model with memory effects measured by Caputo's fractional derivative in time of order 0 < α ≤ 1 , such that when α = 1 we obtain the classical model (without memory). Secondly, we have constructed Lyapunov       www.nature.com/scientificreports/ functions to show the global stability of the equilibrium points in a more general framework where the proposed system takes into account the spatial behavior of populations and memory effect. Taking advantage of the Lyapunov function method we have shown that R 0 plays an important role in determining the global dynamics of the proposed model. We have established the global stability of the two equilibria: E f and E * for different values of α . From epidemiological point of view, this means that the infection will eradicated or persisted while respecting certain restrictions on the parameters. According to our theoretical analysis, we obtained the stability of the equilibria not only for the integer derivative ( α = 1 ) but also for all 0 < α ≤ 1 , which confirms the generality of our system. In addition, fractional derivatives have provided other means of predicting the progression of the disease and, in some cases, affecting the time required to reach stable states. Our future work is to control the vaccination term u to get a better optimal strategy with other fractional derivatives having a non singular kernel.

Methods
As an application of the fractional derivatives a diffusive SIR epidemic model is described. The fixed point theory is adopted for the results related existence and uniqueness of the solution and Lyapunov function theory is utilized for the stability analysis of proposed model. Numerical results are done for the verification of obtained results and it is surety that it will help the researcher in future related fractional order models.

Data availability
The database used and analyzed during the current study are available from the corresponding author on reasonable request.