Global Stability Condition for the Disease-Free Equilibrium Point of Fractional Epidemiological Models

: In this paper, we present a new result that allows for studying the global stability of the disease-free equilibrium point when the basic reproduction number is less than 1, in the fractional calculus context. The method only involves basic linear algebra and can be easily applied to study global asymptotic stability. After proving some auxiliary lemmas involving the Mittag–Lefﬂer function, we present the main result of the paper. Under some assumptions, we prove that the disease-free equilibrium point of a fractional differential system is globally asymptotically stable. We then exemplify the procedure with some epidemiological models: a fractional-order SEIR model with classical incidence function, a fractional-order SIRS model with a general incidence function, and a fractional-order model for HIV/AIDS.


Introduction
Fractional differential equations play an important role in modeling real-life phenomena. By replacing an integer-order derivative with a real-order fractional derivative, often we can fit the system of equations to the real data more efficiently because many dynamical systems cannot be completely described by ODEs. To mention a few of them, we refer applications to bioengineering [1,2], biology [3], Lévy motion [4], harmonic oscillators with damping [5], economy [6,7], and engineering [8,9]. In this work, we are particularly interested in applications of fractional calculus in epidemiological models. This topic has been intensively studied in the recent past, from the well formulation of the problem, the existence of equilibrium points, modeling, and forecasting of epidemiological systems. For example, Refs. [10][11][12] proposed fractional epidemiological models to study the spread of COVID-19 in different countries, Ref. [13,14] investigated the HIV infection, in Ref. [15], a varicella outbreak in China was considered, the spread of dengue fever outbreak in the Cape Verde islands was studied in [16], and in [17] a fractional measles model was proposed. Stability studies were given in e.g., [13,[18][19][20][21] and numerical methods in e.g., [22][23][24]. We also refer to [25] where a review of several fractional epidemiological models was carried out.
An important problem is the study of the global stability of the equilibrium points, in order to better understand the evolution of the disease over time. That is, the system will evolve to the equilibrium point, independently of the starting points. The study of local stability is a relatively simple matter, as it usually involves finding the eigenvalues of the Jacobian matrix and studying their sign. However, the question of global stability is not, in many cases, simple to answer as it usually involves constructing suitable Lyapunov-like functions and there is no routine on how to find them. We emphasize here the fact that the use of the Lyapunov stability theory to establish the global asymptotic stability for fractional differential equations is more complicated than for the ODEs (see, e.g., [26][27][28][29]). This was the motivation to investigate a new method to study global asymptotic stability for dynamical systems described by fractional differential equations. In [30], a novel method was presented. By writing the system in a matricial form and analyzing the matrices involved is such procedure, under some simple assumptions, we can ensure that the equilibrium point is globally stable if the basic reproduction number R 0 is less than 1. The aim of this work is to generalize the main result of [30] to the fractional setting. To our knowledge, this is the first work on global stability following this last approach to the problem.
The paper is organized as follows. In Section 2, we present some concepts and known results needed for this work. Section 3 presents the new contributions of this paper. After deducing some auxiliary lemmas, we prove the main result of this work, Theorem 2. Under some assumptions that can be easily verified for a wide range of epidemiological models given by fractional differential systems, we prove that, if the basic reproduction number is less than 1, then the equilibrium point is globally asymptotically stable. Lastly, in Section 4, we present three examples to show the utility of our research.

Preliminaries
We begin this section with some basic definitions and results of the fractional calculus needed in this work. For more details, we refer the reader to [31,32].
Definition 1. Let f : R + 0 → R be an integrable function. The (left-sided) Riemann-Liouville fractional integral of function f of order α is given by Next, we recall the definition of the generalized Mittag-Leffler function, which is a special function that generalizes the standard exponential function. The Mittag-Leffler function is of great importance in fractional calculus because it arises naturally in the solution of fractional-order differential and integral equations. Definition 3. The Mittag-Leffler function with two parameters is defined by , t ∈ C, with α, β ≥ 0. When β = 1, we define the one parameter Mittag-Leffler function E α (t) := E α,1 (t).
To understand the theory of fractional differential equations, one needs to know properties of these special functions. Its main properties and applications can be found, for example, in [33]. We emphasize here the fact that E α,β (t) can take negative values (cf. [34]).
Recently, we have observed an increasing interest in the Mittag-Leffler function for matrix arguments, since the solution of many systems of differential equations of noninteger order can be expressed using this matrix function. For theoretical properties and a survey on numerical approximation of the matrix Mittag-Leffler function, we recommend the recent paper [35] and the references cited therein.

Definition 4.
Given A ∈ C n×n , the matrix Mittag-Leffler function with two parameters is defined through the convergent series where α, β ≥ 0. If β = 1, we define the one parameter matrix Mittag-Leffler function E α (A) := E α,1 (A).

Remark 1.
For α = β = 1, the matrix Mittag-Leffler function is the matrix exponential, that is, A k k! . Unfortunately, as noticed in [36], there are several works where some properties of the matrix exponential were incorrectly extended to the matrix Mittag-Leffler function and then used to solve certain linear matrix fractional differential equations. One of the properties that cannot be extended to the matrix Mittag-Leffler function is the semigroup property: for given commutating matrices A and B, in general, we have E α (A + B) = E α (A) · E α (B). We note, however, that, if matrices A and B commute and α ≈ 1, We recall now two properties of the matrix Mittag-Leffler function that are useful in the present work (see [35]): if there exists a non-singular matrix P such that To finalize this section, we review some concepts on matrix theory.

Definition 5.
We say that a square matrix A is an M-matrix if the off-diagonal entries are nonpositive and the real parts of all eigenvalues are nonnegative.
Given a square matrix A, the set of eigenvalues of A is denoted by σ(A). The spectral bound of matrix A is defined as m(A) = max{Re(λ) : λ ∈ σ(A)}, where Re(λ) denotes the real part of λ, and the spectral radius of A is defined as ρ(A) = max{|λ| : λ ∈ σ(A)}.
The following result is a fundamental tool in the proof of Lemma 4. where Ω := (a − d) 2 + 4bc. Let e 1 := E α,β (λ 1 ) and e 2 := E α,β (λ 2 ). If Ω and c are not zero, the matrix Mittag-Leffler function of matrix A is We remark that, in Lemma 1, if Ω = 0 or c = 0, a simple formula for the Mittag-Leffler function of a diagonalizable matrix of order 2 can be easily obtained. Theorem 1. (cf. [37]) Let A ∈ R n×n . If the spectrum of A satisfies the relation then lim t→∞ E α (At α ) = 0.

Main Results
Suppose that the epidemiological model under study is described by the fractional differential system with nonnegative initial conditions X(0) = X 0 ∈ R m and I(0) = I 0 ∈ R n , where the components of the vector X denote the number of uninfected individuals (e.g., susceptible, recovered, vaccinated, etc.) and the components of I denote the number of the infected and infectious (the ones that can transmit the disease, such as the asymptomatic but infectious and active infected). In addition, we assume that function F is continuous, G is of class C 1 , and the fractional differential system (1) with initial conditions X(0) = X 0 and I(0) = I 0 admits a unique solution.
Let A := ∂G ∂I (X , 0) and assume that matrix A can be written in the form A = M − D, where M, D are two square matrices with M ≥ 0 (all entries are nonnegative), and D > 0 is a diagonal matrix. The following result was proven in [38]: The value R 0 := ρ(MD −1 ) plays an important role in epidemiological models, and it is known as the basic reproduction number. This number gives the average number of secondary cases produced by one infected individual in a population where all individuals are susceptible to the infection. The following result is well known in the literature. For the convenience of the reader, we present here one possible proof that follows from the fact that the scalar Mittag-Leffler function is completely monotonic ( [39,40]). Since This completes the proof.
The following result is also useful in this work.
from (2), we conclude that d dt E α,α (t) ≥ 0, proving the desired result. Now, we prove the following lemma that shows the applicability of our main result (Theorem 2).
Proof. With our assumptions and using the notations from Lemma 1, we have that b, c, Ω ∈ R + 0 , λ 1 , λ 2 ∈ R − 0 , and λ 2 ≥ λ 1 . First, suppose that Ω = 0 and c = 0. Hence, from Lemmas 2 and 3, we conclude that e 2 := E α,α (λ 2 ) ≥ e 1 := E α,α (λ 1 ) ≥ 0. It remains to be proved that Suppose that d ≥ a (the other case is similar). Then, we just need to prove the first inequality. Since both sides of the inequality are nonnegative, we have that proving the desired. Now, we suppose that c = 0 and a = d. If a < d, then we get and, if a > d, then If c = 0 and a = d, then The following result is a fundamental tool in the proof of Theorem 2.
Lemma 5. Let B ∈ R n×n be an invertible matrix and H : R m+n → R n be a continuous function. Suppose that the fractional differential equation with initial condition I(0) = I 0 ∈ R n , has a unique solution. Then, the solution of this initial value problem satisfies Proof. The proof follows the ideas from [41] (Theorem 7.2). First, observe that To compute Since we get the following: . Thus, Hence, we may conclude that proving the desired result.
We are now in conditions to present a new global stability condition for the DFE of system (1) when R 0 < 1. Knowing that an equilibrium point is globally asymptotically stable with respect to the system that describes the evolution of the uninfected individuals, and with some extra assumptions related to the matrices involved in the system associated with infected and infectious individuals, we can conclude that the equilibrium point is in fact globally asymptotically stable with respect to the complete system. Although the result imposes some restrictions in order to be applied, for many epidemiological models, it can be easily used, as we will illustrate in Section 4.
For C D α 0+ X(t) = F(X, 0), the vector X is globally asymptotically stable;

2.
Function G can be written as G(X, I) = A · I −Ĝ(X, I), whereĜ(X, I) ≥ 0 for all (X, I), and A = ∂G ∂I (X , 0) can be written as A = M − D, where M ≥ 0 and D > 0 is a diagonal matrix; 3.
Proof. First, observe that, since R 0 < 1, then m(A) < 0 (see [38]), and so matrix A is invertible. Since then, by Lemma 5, we get Since the real parts of the eigenvalues of matrix A are negative, from Theorem 1, we get lim t→∞ E α (At α ) = 0, and hence lim t→∞ I(t) = 0.
Since X is globally asymptotically stable with respect to C D α 0+ X(t) = F(X, 0), which in turn is the limiting system of C D α 0+ X(t) = F(X, I), then we get lim t→∞ X(t) = X (t), which completes the proof.

Remark 2.
Note that the assumption E α,α (A) ≥ 0 in Theorem 2 is trivially satisfied if matrix A has dimensions 1 or 2 (by Lemmas 2 and 4, respectively). In addition, if α ≈ 1, since E α,α (A) ≈ exp(A), the above condition also holds for any matrix A ∈ R n×n .

Examples
In this section, we illustrate our main result, Theorem 2, by considering three Caputo fractional-order compartmental models and show that the disease free equilibrium is globally stable in all the cases, whenever R 0 < 1. We stress that usually (see e.g., [21]) the global stability of the disease free equilibrium is proved by considering an appropriate Lyapunov function and LaSalle's invariance principle [42], which are often difficult to apply especially when the model has a considerable number of variables and parameters. Our main result allows us to prove the global stability of the disease free equilibrium in an easier and simpler way. For the numerical implementation of the fractional derivatives, we have used the Adams-Bashforth-Moulton scheme, which has been implemented in the Matlab code fde12 by Garrappa [43].

A Fractional SEIR Model with Traditional Incidence Rate
We start by considering a Caputo fractional-order version of the classical SEIR model that has been applied to describe the transmission dynamics of infectious diseases where there exists a significant latency period during which the individuals are infected but not yet infectious. During this period, the individual is in the so-called exposed compartment E, see e.g., [45]. The other compartments of the model are susceptible S, infected I, and recovered R, and each of them denotes a fraction of the total population. The following assumptions are considered: the birth and death rates are assumed to be equal, and denoted by µ; the incidence rate is the traditional one, given by βSI, where β represents the transmission rate; the latent period is denoted by ε; infected individuals recover at a rate γ and remain recovered with permanent immunity. All parameters are assumed to be positive. The model is given by the following system: The disease-free equilibrium of the model (3) is given by

Following the notation from Section 3, we have
The matrix A can be written as A = M − D with The point X = (1, 0) is globally asymptotically stable for the system of uninfected individuals: It is easy to verify that the function satisfies the second equation of (4). From [41] (Theorem 7.2), the solution of the first equation of (4) is the function Simple computations lead to Thus, (S(t), R(t)) → (1, 0) as t → ∞.
Through adequate numerical simulations, we illustrate the global stability of the disease free equilibrium, whenever R 0 < 1, considering different values of α and initial conditions. In Figure 1, we consider different values for α and the initial condition x 0 from (5).

A Fractional SIRS Model with General Incidence Rate
In the second example, we consider the Caputo fractional-order version of the classical SIRS model from [21], given by the following system: The model considers a homogeneous population divided into three subgroups: susceptible individuals S(t), infected and infectious individuals I(t), and recovered R(t), individuals at time t. The parameters Λ, β, µ, and r, represent the recruitment rate of the population, the infection rate, the natural death rate, and the recovery rate of the infected individuals, respectively. The rate that recovered individuals lose immunity and return to the susceptible class is represented by λ. While contacting with infected individuals, the susceptible become infected at the incidence rate where k 1 , k 2 , and k 3 are nonnegative constants [21]. We remark that system (6) admits a unique positive solution (see [21] (Theorem 7)). The disease-free equilibrium of the model (6) is given by In this case, following the notation from Section 3, Hence, We easily conclude that A < 0 whenever R 0 < 1.
In what follows, we prove that the first condition of Theorem 2 holds, that is, X = (Λ/µ, 0) is globally asymptotically stable for the system of uninfected individuals: The solution of the fractional differential equations (7) is the functions and Obviously R(t) → 0, as t goes to infinity. Now, we prove that S(t) → Λ/µ, as t → ∞. First, observe that For the other term inside the integral, we get To evaluate this integral, we use the known formula involving the Beta function: With the change of variable u = s/t, we get Thus, we prove that the solution S(·) is given by Observe that, as t goes to infinity, To sum up, it remains to prove that the double sum converges to zero. For that purpose, we recall the concept of the Mittag-Leffler function of two variables (cf. [46]) E α,β (x, y, ·). With such notation, we can write which converges to zero, as t goes to infinity, by [46] (Theorem 3.1). This proves the desired conclusion. We also remark that, by Lemma 2, E α,α (A) is nonnegative. Therefore, by Theorem 2, the disease-free equilibrium of the model (6) is globally asymptotically stable.
The stability of the disease free equilibrium for model (6) is illustrated in Figure 3.

A Modified Fractional SICA Model for HIV/AIDS
In this example, we consider a modified Caputo fractional-order model for HIV/AIDS, based on the model proposed in [14,47]. We show that this fractional model satisfies the conditions of Theorem 2 and, through some numerical simulations, we illustrate the global stability of the disease free equilibrium, when R 0 < 1.
In this model, the total population is assumed to be homogeneous and divided into four mutually-exclusive compartments: susceptible individuals (S); HIV-infected individuals with no clinical symptoms of AIDS but able to transmit HIV to other individuals (I); HIV-infected individuals under antiretroviral (ART) treatment (the so called chronic stage) with a viral load remaining low (C); and HIV-infected individuals with AIDS clinical symptoms (A). Analogously to the assumption made in [47], we consider that individuals in the chronic stage C have a very low viral load and do not transmit HIV infection [48], but, differently from [47], we assume that individuals with AIDS A, due to their higher viral load, may transmit HIV virus, at a rate η A β with η A > 1. Therefore, effective contact with people infected with HIV is at a rate λ, given by where β is the effective contact rate for HIV transmission. We assume that the recruitment rate is equal to the natural death rate and is denoted by µ. The following assumptions are the same as in [14]. HIV-infected individuals with no AIDS symptoms I progress to the class of individuals with HIV infection under ART treatment C, at a rate φ, and HIV-infected individuals with AIDS symptoms are treated for HIV at rate γ. Individuals in the class C leave for the class I, at a rate ω. HIV-infected individuals with AIDS symptoms A that start treatment move to the class of HIV-infected individuals I, moving to the chronic class C only if the treatment is maintained. HIV-infected individuals with no AIDS symptoms I that do not take ART treatment progress to the AIDS class A, at rate ρ. Only HIV-infected individuals with AIDS symptoms A suffer from an AIDS induced death, at a rate d. The total population at time t, denoted by N(t), is given by N(t) = S(t) + I(t) + C(t) + A(t). The Caputo fractional-order system that describes the previous assumptions is given by The disease free equilibrium of system (9) is given by Using the notation from Section 3, we have The matrix A can be written as In this case, we will prove that X = (1, 0) is globally asymptotically stable for the system of uninfected individuals: The solution of (11) is and Similarly to Example 4.1, we can prove that the disease-free equilibrium of the model (9) is globally asymptotically stable.
We now show, using numerical simulations, that for different values of α and initial conditions, the global stability of the disease free equilibrium holds, whenever R 0 < 1.
In Figure 4, we consider different values for α and the initial conditions x 0 from (12). The global stability of the disease free equilibrium (10) is illustrated in Figure 5, considering different initial conditions x i , i = 0, . . . , 7, given by (12).   (9), considering α = 0.9 and different initial conditions x i , i = 0, . . . , 7, from (12), on the x-axis S and on the y-axis I + C + A.

Conclusions
A new and simple result for the global stability for the disease-free equilibrium of fractional epidemiological models is presented. We highlight here the fact that the approach available in the literature so far involves the determination of an appropriate Lyapunov function, very laborious computations and, in the end, the application of LaSalle's invariance principle. Our new method uses only basic results from matrix theory and some well-known results from fractional-order differential equations.
We also remark that the applicability of our main result, Theorem 2, is only possible if the matrix A satisfies the condition E α,α (A) ≥ 0. We proved that, under the assumptions of Theorem 2, this condition holds if matrix A has dimensions 1 or 2. It would be interesting to check under what conditions we can guarantee that E α,α (A) ≥ 0 if the matrix A has dimensions greater than 2. Since many epidemiological models divide the population into subpopulations of epidemiological significance where the number of the infectious compartments are at most 2, our main result can be applied to a wide variety of epidemiological models.