Dynamics of Fractional-Order Epidemic Models with General Nonlinear Incidence Rate and Time-Delay

: In this paper, we study the dynamics of a fractional-order epidemic model with general nonlinear incidence rate functionals and time-delay. We investigate the local and global stability of the steady-states. We deduce the basic reproductive threshold parameter, so that if R 0 < 1, the disease-free steady-state is locally and globally asymptotically stable. However, for R 0 > 1, there exists a positive (endemic) steady-state which is locally and globally asymptotically stable. A Holling type III response function is considered in the numerical simulations to illustrate the effectiveness of the theoretical results.


Introduction
The study of the spread of disease is one of the main directions of mathematical modeling. There has been enormous effort to study the dynamics of infectious disease and various generalizations of the epidemic models that exist in the literature. More recently, there has been an increasing interest in epidemic models that involve time-delays and fractional-order derivatives [1][2][3][4][5]. The need for fractional-order differential equations (FODE) comes from real-world application problems. These equations are generalizations of ordinary differential equations and believed to take their origin from the question of L'Hopital to Leibniz in a report letter at the end of the 17th century [6]. Many famous mathematicians, such as Liouville, Riemann, Fourier, and Abel contributed to the Theory of Fractional Calculus. FODE models have advantages over classical ODE and/or delay differential equation models, since integer-order derivatives are used to gain information only about local features of a state, whereas fractional-order derivatives describe the whole space [7]. In other words, in the FODE models, the next certain position for a physical phenomenon does not only depend on the current state, but also on all historical states. Thus, fractional-order models not only give more realistic biological models involving memory, but also enlarge the stability region of the states. It is worth mentioning that FODEs have been successfully applied to various fields, such as finance [8], control theory [9], ecology [10], and medicine [11][12][13][14][15][16][17].
Incidence rate refers to the frequency with which a disease occurs over a specified period of time. It is numerically defined as the number of new cases of a disease within a time period, as a proportion of the number of people at risk for the disease. This rate provides the capacity to anticipate future incidents and plan accordingly. Thus, the incidence rate function type plays a key role in understanding the dynamics of infectious diseases. However, in most cases, it has been observed that the incidence rate starts to slow down after some time. In the last few decades, there have been studies considering nonlinear incidence rates given by a general function F (S, I) which satisfies certain mild conditions [18][19][20][21][22][23][24]. On the other hand, there are fractional-order epidemic models with a bilinear incidence rate and standard incidence rate, as well as specific nonlinear incidencerate functions [5,16,[25][26][27]. However, these models need to be modified by means of the general incidence-rate function. More recently, the authors in [26] considered a fractionalorder Susceptible-Exposed-Infected-Recovered (SEIR) model with a generalized incidence rate function of the type S f (I), where f satisfies some certain conditions. Furthermore, Lahrouz et al. studied a fractional-order Susceptible-Infected-Recovered (SIR) model with a general incidence rate function and carried out Mittag-Leffler stability and bifurcation analysis for the model with and without delays [24]. It should be noted that the authors considered the model with a delay on the recovered/removed population (R(t)) and investigated the global dynamics.
In the modeling of infectious diseases, it is assumed that there exists a latent period before an infected individual becomes infectious. In other words, it takes some time for an infected person to be able to transmit a disease to another susceptible host. This time-interval is called the latency period, and is modeled by delay differential equations. There are many interesting papers which study epidemic models with a constant delay in integer-order [28,29] and fractional-order differentiation [7,30].
Historically, the study of a so-called Susceptible-Infected-Recovered (SIR) model goes back to the work by Kermack and McKendrick [31]. Since then, this model has been widely studied to predict the dynamics of different infectious diseases. In the present paper, we introduce the fractional-order SIR model with a general incidence rate and time-delay. In Section 2, we provide the epidemic model and in Section 3, we study the local stability of the steady-states. In Section 4, we discuss the global stability of the endemic model. In Section 5, we display some numerical simulations, and then we conclude in Section 6.
Some preliminaries and definitions about fractional-order derivatives are given in the Appendix A.

Fractional Epidemic Model with a General Incidence Rate
Consider a system of delay differential equations for an SIR model Ψ is the recruitment rate, µ is the natural death rate, ρ is the recovery rate, δ is the diseaseinduced death rate, and τ > 0 is the delay due to the latency period. We assume the nonlinear incidence rate function F (S, I) satisfies the following conditions: is always positive, continuous, differentiable, and monotonically increasing: (A3) F (S, I) is concave with respect to I: One can easily check that the following incidence-rate functions satisfy the conditions (A1)-(A3) : βSI, βSI S + I + R , βSI 1 + α 1 S + α 2 I + α 3 SI , βS n I 1 + ηS n , n ≥ 2.
In this paper, we do not consider a specific incidence rate function, but a general one that satisfies mild conditions. Consequently, the results of this paper can be applied with any incidence rate function that satisfies the conditions (A1)-(A3). In Section 5, we implement the theoretical results with a Holling type III functional response, that is, F (S, I) = βS 2 I 1 + ηS 2 . After integrating from t 0 to t, the model (1) reads as where the constant kernel function K satisfies K(t − s) = 1. For some biological systems, fractional differential equations are more consistent with real phenomena than integerorder models. This is because fractional derivatives and integrals can be used to describe the memory and hereditary properties. In order to consider the historical effects in the above system, we modify the kernel function K as follows: where 0 < α ≤ 1 and Γ is the Gamma function, as specified in Appendix A. Introducing a power-law memory kernel function that allows near-past events to have a greater effect on the model than distant past events, by Definition A1, the modified model now reads as: By applying Lemma A1, see Appendix A, the last model is equivalent to the following SIR model: with initial conditions S(0) > 0, R(0) > 0, I(s) = χ(s) > 0, s ∈ [−τ, 0], and χ(s) is a smooth function. The fractional derivative α ∈ (0, 1] is defined by the Caputo sense, so introducing a convolution integral with a power-law memory kernel (3) benefits in describing memory effects in dynamical systems. The decaying rate of the memory kernel depends on α. A lower value of α corresponding to more slowly-decaying time-correlation functions leads a long memory. Therefore, as α → 1, the influence of memory decreases.
Under the condition (A1), it is easy to check that the model (5) has a disease-free steady-state solution E 0 = (S 0 , I 0 ) = (Ψ/µ, 0). In the analysis of epidemic models, it is of utmost importance to find so-called basic reproductive numbers, R 0 , which serves as a threshold value. We compute R 0 by using the next-generation matrix technique [32]: The existence of a unique endemic steady-state solution E * = (S * , I * ) usually depends on R 0 and satisfies the following algebraic equations: The following lemma [23] ensures the uniqueness of the positive steady-state solution holds for all S ∈ (0, S 0 ), then (5) has a unique positive endemic steady-state solution E * (S * , I * ). If R 0 ≤ 1, and then the disease-free steady-state solution E 0 (S 0 , I 0 ) is the only equilibrium state of the system (5).

Local Stability of the Epidemic Model
The corresponding linearized system of model (5) at any steady-state (S * , I * , R * ) is described as Taking the Laplace transform on both sides of (8) yields Here, S(s), I(s), R(s) are the Laplace transform of S(t), I(t), and R(t) with S(s) = L(S(t)), I(s) = L(I(t)) and R(s) = L(R(t)). The above Equation (9) can be written as and ∆(s) is the characteristic matrix for the system (8) at (S * , I * , R * ). The stability conditions are also related to the fractional-order α.
Proof. The characteristic equation at E 0 is described by • Case (i), when τ = 0.
Therefore, the disease-free equilibrium point is locally asymptotically stable if In this case, Equation (10) becomes Put Separate the real and imaginary parts of the above equation If we square and add both sides, one can have Obviously, w α > 0, cos πα 2 > 0 and the assumption ∂F ∂I (S 0 , 0) < ρ + µ + δ imply that the Equation (12) has no real solutions. Furthermore, we can infer that Equation (11) has no purely imaginary roots (eigenvalues) for any τ > 0. From Lemma 1 in [10] and Corollary 3 in [33], the equilibrium point E 0 is asymptotically stable for delay τ ≥ 0.
For τ = 0, it has been shown in [34] that the steady-state is locally asymptotically stable if all the eigenvalues satisfy the condition | arg λ i | > απ 2 . The fractional order enlarges the stability region of the solution of the model. Next, we discuss the local stability at an endemic equilibrium point E 1 , at which the characteristic equation is Therefore, where, The discriminant of (13) is given by Proof. Case (i) τ = 0, Equation (13) becomes Based on [35], the point E 1 is asymptotically stable if the following hold: Case (ii) τ = 0, let us assume that the equation (13) has pure imaginary roots s = ±iw, for τ > 0, then we will finally obtain that if A 2 Similarly to the proof of Theorem 1, it is hence omitted. Therefore, all the roots of Equation (13) are the negative real part for all τ > 0.

Holling Type III Functional Response
In this subsection, we consider the system (5) in the case with the Holling type III functional response. That is, we take F (S, I) = βS 2 I 1 + ηS 2 , as an incidence rate function. We derive the linearized system, and applying the Laplace transform, one can finally get the following characteristic matrix at (S * , I * , R * ), Proof. The characteristic equation at E 0 is described by • Case (i) when τ = 0.
Therefore, the disease-free equilibrium point is locally asymptotically stable if In this case, Equation (14) becomes Put s = iw, w > 0 Separate the real and imaginary parts of the above equation If we square and add both sides, one can have Obviously, w α > 0, cos πα 2 > 0 and the assumption β(S 0 ) 2 1+ηS 2 0 < ρ + µ + δ, imply that the Equation (16) has no real solutions. Furthermore, we can infer that Equation (15) has no purely imaginary roots for any τ > 0. From Lemma 1 in [10] and Corollary 3 in [33], the equilibrium point E 0 is asymptotically stable for delay τ ≥ 0.
The local stability at an endemic equilibrium point E 1 for this case, follows the same proof of Theorem 2.

Global Stability Analysis
Now, we study the global stability of the infection-free and endemic steady-states using a certain Lyapunov functional. Summing all equations in (5) yields One can easily show that the last inequality implies that where E α is the standard Mittag-Leffler function, see Appendix A for the definition. The last inequality implies that if S(0) In other words, the region is a positively invariant set of the system (5). Thus, we consider the model (5) in a biologically feasible region Ω.
Proof. Let us consider the following Lyapunov functional:

t] [I(t)] + I(t).
Note that V 1 (S 0 , 0) = 0 and (µ + ρ + δ)I α  1 (S(t), I(t)) in fractional-order along the trajectories of (5) yields: By theorem hypotheses, and from the monotonicity of F (S, I) with respect to S, we have The concavity of F (S, I) with respect to I leads to F (S, I) ≤ I ∂F (S, 0) ∂I , and hence we have Thus, 0 D α V 1 (S(t), I(t)) ≤ 0 for all S, I > 0. On the other hand, one can show that the largest invariant set of 0 D α V 1 (S(t), I(t)) = 0 is the singleton E 0 (S 0 , I 0 ), that is, the diseasefree steady-state solution. Thus, by the Lyapunov-Lasalle theorem for FODE (Lemma 4.6 in [11]), we conclude that E 0 (S 0 , I 0 ) is globally asymptotically stable. This finalizes the proof.
In what follows, the auxiliary lemma on the estimate of a Volterra-type Lyapunov function will be useful to prove the global stability of the endemic steady-state solution.
To carry out further stability analysis, we need the following condition: Theorem 5. Assume that (A1)−(A4) hold and R 0 > 1. Then, the endemic positive equilibrium state E * (S * , I * ) is globally asymptotically stable.
Consider the following Lyapunov function: I * ≥ 0 and F (S * ,I * ) x * is a Volterra-type Lyapunov function for any x, x * ≥ 0. Due to monotonicity of the function F (S, I), we have F (S(t), is a Lyapunov functional. Differentiating V 2 (S(t), I(t)) in fractional-order along the trajectories of (5) and implementing Lemma 2 together with (7) yields After factoring the last term of the above expression, we obtain Due to the monotonicity of F (S, I) with respect to S, we have Moreover, since x − 1 − ln x is a Volterra-type Lyapunov function with x * = 1, we have 1 − x + ln x ≤ 0. Thus, the second, third, and fourth lines in (17) are non-positive. Finally, the condition (A4) implies that Hence, 0 D α V 2 (S, I) ≤ 0 holds for S, I > 0, since µS * and F (S * , I * ) are positive numbers. Thus, by the similar arguments as in Theorem 4, we conclude that E * (S * , I * ) is globally asymptotically stable.

Numerical Simulations
To illustrate the theoretical results obtained in the previous sections, we provide a numerical scheme for finding the solution of the fractional-order ODEs. To this end, we follow the modified Adams-Bashforth-Moulton method which has proven to be an effective tool to numerically solve fractional-order delay differential equations [5,7,9]. Let us highlight the main features of the algorithm of the modified Adams-Bashforth-Moulton method. In this regard, consider the following fractional-order delay differential equation.
where x(t) = [x 1 (t), x 2 (t), · · · , x n (t)] T and f : J × R n × R n → R n is a Lipschitz-continuous function. This initial value problem is equivalent to the following Volterra-integral equation For given mesh points, T = {t −m , t −m+1, , · · · , t −1 , t 0 , t 1 , · · · t n } with step-size h = τ/m such that t −m = −τ, t 0 = 0 and t n = T let x F (t i ) = φ(t i ) and x F (t i − τ) = φ(t i h − mh) for i = −m, −m + 1, · · · , 0. Assume further that x F (t i ) ≈ x(t i ) are calculated for i = −m, −m + 1, · · · , n. Then the numerical scheme for the integral Equation (19) is given by Next, we consider the nodes t i , i = 0, · · · , n = 1, with the weight function (t n+1 − s) α−1 and apply the product trapezoidal quadrature formula to approximate the integral of (20) on the right-hand side. where Hence, we can formulate the numerical scheme for the IVP (18) as follows.
The term x F (t n+1 ) is approximated by the so-called predictor, x p h (t n+1 ), evaluated by the product rectangle rule by means of the Equation (20). where For the simulation part, we use the algorithm (22) and Holling type III function, F (S, I) = βS 2 I 1 + ηS 2 , as an incidence rate response for the fractional-order SIR model (5). One can easily check that the conditions (A1)-(A3) are satisfied. Thus, by using the next-generation matrix technique, one can compute the basic reproduction number as R 0 = βΨ 2 µ 2 + ηΨ 2 . Thus, let us write the numerical scheme (22) for the model (5) with Holling type III function: where q i,n+1 is defined as before.

Concluding Remarks
In the present paper, we studied local and global stability analyses of the delayed fractional-order SIR epidemic model with a general incidence rate function that satisfies some mild conditions. The fractional derivative is considered in the Caputo sense, which is suitable for initial-value problems. The disease-free steady-state solution of the underlying model is locally and globally asymptotically stable when R 0 ≤ 1, whereas when R 0 > 1, there exists a positive endemic steady-state solution which is locally and globally asymptotically stable. To illustrate our theoretical results, we considered the Hooling type III functional response as the incidence-rate function and used the modified Adams-Bashforth-Moulton method to simulate an example. Fractional-order differential equations are at least as stable as their integer-order counterparts. Further, the presence of a fractional differential order and time-delay in a differential equation can result in a significant increase in the complexity of the observed behavior, and the solution continuously depends on the previous states.