Asymptotic stability of an epidemiological fractional reaction-diffusion model

Lamia Djebara; Salem Abdelmalek; Samir Bendoukha

doi:10.1515/dema-2022-0224

Article Open Access

Asymptotic stability of an epidemiological fractional reaction-diffusion model

Lamia Djebara , Salem Abdelmalek and Samir Bendoukha

Published/Copyright: June 15, 2023

Published by

Become an author with De Gruyter Brill

Author Information Explore this Subject

From the journal Demonstratio Mathematica Volume 56 Issue 1

Abstract

The aim of this article is to study the known susceptible-infectious (SI) epidemic model using fractional order reaction-diffusion fractional partial differential equations [FPDEs] in order to better describe the dynamics of a reaction-diffusion SI with a nonlinear incidence rate describing the infection dynamics of the HIV/AIDS virus. We initially examined the nonnegativity, global existence, and boundedness for solutions of the proposed system. After determining that the proposed model has two steady states, we derived sufficient conditions for the global and local asymptotic stability of the equilibrium of the proposed system and their relationship to basic reproduction in the case of fractional ordinary differential equations and FPDEs by analyzing the eigenvalues and using the appropriately chosen Lyapunov function. Finally, we used numerical examples to illustrate our theoretical results.

Keywords: asymptotic stability; fractional calculus; basic reproduction number; fractional Lyapunov method; Lyapunov functional

MSC 2010: 35k45; 35k57

1 Introduction

In this article, we are interested in studying a generalized nonlinear fractional epidemic reaction-diffusion system model, where we present the following proposed fractional model describing the transmission of a communicable disease between individuals such as HIV and AIDS. The humans population is divided into two epidemiological classes denoted by susceptible-infective (SI), which is an extension of the work proposed in [1,2].

(1) D t α 0 C S ( t ) − d 1 Δ S = Λ − λ S ϕ ( I ) − μ S ≕ F 1 ( S , I ) , D t α 0 C I ( t ) − d 2 Δ I = λ S ϕ ( I ) − σ I ≕ F 2 ( S , I ) ,

defined over R + × Ω . In the context of this work, we denote by Ω an open bounded subset of R n with a piecewise smooth boundary ∂ Ω . In addition, Δ denotes the Laplacian operator over Ω , the parameters d 1 , d 2 ≥ 0 represent the diffusion coefficients, and D t α 0 C denotes the Caputo fractional derivative over ( 0 , + ∞ ) with the fractional differentiation order α being limited to the interval ( 0 , 1 ] . S ( x , t ) and I ( x , t ) represent, respectively, the relative density of the population of SI individuals at time t and location x . Λ > 0 denotes the birth rate, sometimes referred to as the recruitment rate of the population, and μ denotes the natural death rate. In addition, λ denotes the infection rate, and σ represents the recovery rate from the disease.

We assume bounded and nonnegative continuous initial data

(2) S 0 ( x ) = S ( x , 0 ) , I 0 ( x ) = I ( x , 0 ) in Ω ,

where S 0 ( x ) and I 0 ( x ) are in C 2 ( Ω ) ∩ C ( Ω ¯ ) , and the following Neumann boundary conditions on R + × ∂ Ω

(3) ∂ S ∂ ν = ∂ I ∂ ν = 0 .

Without going into too much detail regarding the physical interpretations of the parameters, let us assume that

μ > 0 , σ > 0 , λ > 0 , and Λ > 0 .

We would, however, like to emphasize that these conditions are realistic. The nonlinearity ϕ is assumed to be a continuously differentiable function on [ 0 , + ∞ ) satisfying the criteria

(4) ϕ ( 0 ) = 0 ,

and

(5) 0 < I ϕ ′ ( I ) ≤ ϕ ( I ) for all I > 0 .

It can be seen that when setting α = 1 in our present work, we obtain the model studied in [1,2] as a model on the spread of infectious diseases. Since it inception in 1927, this model has attracted the interest of many researchers. Some authors [3–5] studied a simple model proposed by nonlinear iterations of the formula S ϕ ( I ) with d 1 = d 2 = 0 and α = 1 . Similar work and results can be found in [6–8] with various formulas of the form S ϕ ( I ) . However, these studies only considered the simple diffusion-free case, while Djebara et al. [1] studied the asymptotic stability of the classical reaction-diffusion SI epidemic model, with its first-order time derivative, where d 1 ≠ d 2 , α = 1 and with a nonlinear incidence rate. Whereas Akdim et al. [9] studied the local and global stability of the disease-free and endemic equilibrium of our system in the case of d 1 = d 2 = 0 and α ∈ ( 0 , 1 ] according to the basic reproduction number with S ϕ ( I ) = S β 1 − β 2 I m + I I . Lu et al. [10] are the only ones who have studied and demonstrated the global existence, nonnegativity, and stability of solutions for class of fractional order derivatives susceptible-infective-recovered (SIR) epidemic models with a general incidence rate function, which satisfy stronger conditions than ours.

The fractional reaction-diffusion order systems on which we base our work are an extension of the classical integer order reaction-diffusion systems. It should be noted that it has advantages over the classical integer-order reaction-diffusion systems, in that it possesses memory and genetic properties that are not found in these latter systems, which are usually applied by a lot of biological systems. Furthermore, these systems have a more accurate description of population models than integer-order models.

Time-fractional systems of differential equations have been the focus of countless studies over the last few decades. This interest stems from the realization that such systems appear in a wide range of applications across various disciplines such as biological mathematics, physics, and biology. Most recently, time-fractional systems attracted the interest of researchers in the field of nonlinear dynamics, which has resulted in a vast amount of valuable research. Interested readers may refer to [11–15].

One of the most important issues to be studied for any type of model including system (1)–(3) is to determine its set of equilibrium and assess their local and global asymptotic stabilities. In this article, we rely on the basic stability theory of dynamical systems containing a fractional derivative along with the linearization and direct Lyapunov methods. In order to establish the system’s asymptotic stability, we identify its two positive equilibrium and assess their asymptotic stability in two different cases with respect to the reproduction number. Section 2 of this article gives some preliminary definitions and properties related to fractional calculus, while in Section 3, the positivity and boundedness of the solutions are studied. Then, Section 4 establishes sufficient conditions for the local asymptotic stability and instability of the disease-free equilibrium and the endemic equilibrium. In Section 5, we use an appropriate Lyapunov functional to prove that the two steady state solutions of the model are globally asymptotically stable. Finally, Section 6 presents some numerical simulations to illustrate our theoretical results of the stability of the solutions of the epidemiological fractional reaction-diffusion model (1)–(3).

2 Preliminaries in fractional calculus

Throughout this work, we adopt the following definitions, and we need the following lemmas and stability results to prove the main results, i.e., Propositions 7, 14, and 15 and Theorems 18 and 21.

Definition 1

[16,17] Define the Caputo fractional derivative of a function x ( t ) of order α > 0 for t > t 0 , which belongs to class C n by

(6) D t α t 0 C x ( t ) = 1 Γ ( n − α ) ∫ t 0 t x ( n ) ( τ ) ( t − τ ) α − n − 1 d τ = ∫ t 0 t ( t − τ ) n − α − 1 x ( n ) ( τ ) Γ ( n − α ) d τ = I n − α x ( n ) ( t ) ,

where Γ is the gamma function Γ ( β ) = ∫ 0 + ∞ e − t t β − 1 d t , n = min { i ∈ N / i > α } , and the fractional integral of x ( t ) , t > t 0 , of order γ ∈ R + is defined as follows:

I γ x ( t ) = ∫ t 0 t ( t − τ ) γ − 1 x ( τ ) Γ ( γ ) d τ .

Note that the Caputo derivative in time is a good choice to include long-term memory effects as well because the well-known Caputo derivative fulfills the following properties:

We have D t α t 0 C ( x ( t ) ) = 0 , for all x ( t ) be a constant function.
The Caputo derivative is linear, where
D t α t 0 C ( μ x ( t ) + ν y ( t ) ) = μ D t α t 0 C ( x ( t ) ) + ν D t α t 0 C ( y ( t ) )
for all x , y : [ a , b ] → R such that D t α t 0 C ( x ( t ) ) and D t α t 0 C ( y ( t ) ) exist in [ a , b ] .

Definition 2

Assuming two positives parameters γ , δ , and z ∈ C , the Mittag-Leffler function of γ and δ is defined as follows:

E γ , δ ( z ) = ∑ i = 0 ∞ z i Γ ( γ i + δ ) .

For δ = 1 , E γ , 1 ( z ) = E γ ( z ) = ∑ i = 0 ∞ z i Γ ( γ i + 1 ) and if γ = δ = 1 , then E 1 , 1 ( z ) = exp ( z ) = ∑ i = 0 ∞ z i i ! .

Lemma 3

[16,18] For two parameters, the Laplace transform L ( . ) of Mittag Leffler function is

(7) L ( t δ − 1 E γ , δ ( a t γ ) ) = s γ − δ s γ − a ,

where t ≥ 0 , s , a ∈ R , R ( s ) > ∣ a ∣ 1 γ , and

(8) L ( D t α 0 C V ( x ) ) = s α L ( V ( x ) ) − ∑ i = 0 n − 1 V ( i ) ( 0 ) s α − i − 1 .

Theorem 4

[19, Theorem 3.7, page 4] Consider the fractional Caputo autonomous dynamic system

(9) D t α t 0 C S ( t ) = F 1 ( S , I ) , D t α t 0 C I ( t ) = F 2 ( S , I ) in R + ,

with initial condition ( S 0 , I 0 ) , where α ∈ ( 0 , 1 ] , F = ( F 1 , F 2 ) : Ω ⊂ R 2 → R 2 . If F satisfies the local Lipschitz condition with respect to ( S , I ) , then there exists a unique solution of (9) on [ 0 , + ∞ ) × Ω . The constant ( S ∗ , I ∗ ) is said to be an equilibrium point of (9) if and only if

(10) F 1 ( S ∗ , I ∗ ) = F 2 ( S ∗ , I ∗ ) = 0 .

Lemma 5

The asymptotic stability of the point ( S ∗ , I ∗ ) established is subject to

(11) ∣ arg ( λ i ( J ( S ∗ , I ∗ ) ) ) ∣ > α π 2 , for a l l i = 1 , 2 ,

where arg ( ⋅ ) is the argument of a complex number, λ i are the eigenvalues, and J is the Jacobian matrix.

For more on the previous lemma, the reader may refer to [20], where the authors established general conditions for the asymptotic stability of fractional reaction-diffusion models.

Lemma 6

[11] If an equilibrium point ( S ∗ , I ∗ ) of (9) is locally asymptotically stable for the following system:

(12) d d t S ( t ) = F 1 ( S , I ) , in R + , d d t I ( t ) = F 2 ( S , I ) in R + ,

then, it is also locally asymptotically stable for (9).

Corollary 7

[11,20] In the diffusion case, if an equilibrium point ( S ∗ , I ∗ ) of (9) is locally asymptotically stable for the integer system

∂ S ∂ t − d 1 Δ S = F 1 ( S , I ) , ∂ I ∂ t − d 2 Δ I = F 2 ( S , I ) o n R + × Ω ,

then it is also locally asymptotically stable for

D t α 0 C S − d 1 Δ S = F 1 ( S , I ) , D t α 0 C I − d 2 Δ I = F 2 ( S , I ) ,

on R + × Ω .

Lemma 8

[21–24] Let D be a closed and bounded set. Each solution of t 0 C D t α U ( t ) = F ( U ) that starts from a point in D stays in D over time. If there exists a function V ( U ) : D → R with continuous first partial derivatives, then the following condition is satisfied:

D t α t 0 C V ( U ) ≤ 0 .

Let M be the largest invariant set of

E = { D t α t 0 C V ( U ) = 0 , U ∈ D } .

Then, if M = { E ∗ } , every solution in D tends to E ∗ as t goes to infinity.

Lemma 9

[25] Suppose that S ( t ) ∈ R + is a continuous and differentiable function. For all t ≥ t 0 ,

D t α t 0 C S 2 ( t ) ≤ 2 S ( t ) D t α t 0 C S ( t ) ,

with 0 < α ≤ 1 .

Lemma 10

[24,26] Suppose that S ( t ) ∈ R + is a continuous and differentiable function. Then, for all t ≥ t 0 , α ∈ ( 0 , 1 ] , and S ∗ ∈ R + , we have

D t α t 0 C S − S ∗ − S ∗ ln S S ∗ ≤ 1 − S ∗ S D t α t 0 C S .

3 Main results

3.1 Positivity and boundedness in the diffusion-free case

It is important to prove the existence and uniqueness of nonnegativity solutions and boundedness for our system before studying its stability. Let us assume that the initial conditions ( S 0 , I 0 ) ∈ R ≥ 0 2 . Note that for ( S , I ) ∈ R ≥ 0 2 , we have

D t α t 0 C S ( t ) ∣ S = 0 = F 1 ( 0 , I ) = Λ > 0 , D t α t 0 C I ( t ) ∣ I = 0 = F 2 ( S , 0 ) = S ϕ ( 0 ) = 0 .

Hence, the nonnegative quadrant R ≥ 0 2 is an invariant set. By dropping the diffusion terms, the proposed system reduces to the following system of ordinary differential equations:

(13) D t α 0 C S ( t ) = Λ − λ S ϕ ( I ) − μ S ≕ F 1 ( S , I ) , D t α 0 C I ( t ) = λ S ϕ ( I ) − σ I ≕ F 2 ( S , I ) in R + .

(14) S ( 0 ) = S 0 ≥ 0 , I ( 0 ) = I 0 ≥ 0 .

In the following subsections, we define an invariant region for the system, identify the system’s equilibrium with their relation to the basic reproduction number R 0 the most important arithmetic quantity in infectious diseases and biologically, R 0 can be interpreted as the average number of secondary infections produced during the period of infection, the value of which was previously calculated in [1] by means of the next generation matrix method formulated in [27], which is given by R 0 = ρ ( F V − 1 ) = λ σ Λ μ ϕ ′ ( 0 ) . We also demonstrate the global presence of solutions in time, and investigate the local stability of the system in ordinary differential equation (ODE) and partial differential equation (PDE) scenarios.

3.2 Invariant regions

Throughout this article, we use the total population size N ( t ) = S ( t ) + I ( t ) and let σ 0 = min ( σ , μ ) . We also define the region

D = { ( S , I ) : S , I ≥ 0 and S + I ≤ Λ σ 0 } .

The following proposition shows that D is an invariant region of the fractional-order system 13–14 that ensures the boundedness of the system’s solutions.

Proposition 11

The region D is a nonempty, attractive, and positively invariant and for the system 13–14 has a unique solution on R + .

Proof

We start by summing together all the equations of the fractional-order system 13–14, which yields the following fractional order derivative of N ( t ) :

D t α t 0 C N ( t ) = D t α t 0 C S + D t α t 0 C I ≤ Λ − σ 0 N .

Hence,

D t α t 0 C N ( t ) ≤ Λ − σ 0 N .

Applying the Laplace transform and using (8) imply that (see also [19, Lemma 3, page 4] by the comparison principle)

N ( t ) ≤ Λ σ 0 ( 1 − E α ( − σ 0 t α ) ) + N 0 E α ( − σ 0 t α ) ,

where E α is the Mittag-Leffler function. Substituting the value of N yields

( S + I ) ( t ) ≤ ( S + I ) ( 0 ) − Λ σ 0 E α ( − σ 0 t α ) + Λ σ 0

for t ≥ 0 . If the initial states satisfy ( S + I ) ( 0 ) ≤ Λ σ 0 , since 0 ≤ E α ( − σ 0 t α ) ≤ 1 , then ( S + I ) ( t ) ≤ Λ σ 0 , which implies that

limsup t → ∞ N ( t ) ≤ Λ σ 0 .

Hence,

limsup t → ∞ S ( t ) ≤ Λ σ 0

and

limsup t → ∞ I ( t ) ≤ Λ σ 0 , for t ≥ 0 .

Consequently, the solutions of system 13–14 are bounded for t ≥ 0 . This completes the proof of Proposition.□

As a result, region D is positively invariant and attractive within R ≥ 0 2 . Therefore, D is the biologically feasible region of the system as D is sufficient to consider the dynamics of the model within D , where the existence and uniqueness results hold for the system. The existence and uniqueness of solutions for system 13–14 in R + is deduced by using [28, Theorem 3.1, page 10] along with Lemma 4 (see also [1,2]), where we can conclude that F = ( F 1 , F 2 ) satisfies the local Lipschitz condition with respect to ( S , I ) ( t ) in R ≥ 0 2 .

3.3 Positivity and boundedness in the reaction-diffusion case

Proposition 12

Let α ∈ ( 0 , 1 ] , and for any nonnegative initial data ( S 0 , I 0 ) ∈ C ( Ω ¯ ) × C ( Ω ¯ ) , there exists a unique nonnegative global solution in time of the system (1)–(3). In addition, there exist two positive constants T and ℳ such that ∀ t > T ,

(15) ∥ S ( t , . ) ∥ L ∞ ( Ω ) + ∥ I ( t , . ) ∥ L ∞ ( Ω ) ≤ ℳ .

Proof

By condition (4), it follows that F 1 ( S , I ) = Λ − μ S − λ S ϕ ( I ) and F 2 ( S , I ) = λ S ϕ ( I ) − σ I , and satisfy the local Lipschitz conditions. Hence, using the contraction-mapping principle method (see [20, pp. 3], [29, Proposition 4, pp. 4], and [30, Theorem 3.1]), we can easily establish the existence and uniqueness of a mild solution to system (1)–(3) in [ 0 , T ) . It remains to prove that the solution is bounded. Integrating the equations of (1) over Ω , applying Green’s formula and adding the resulting two identities, we have

D t α 0 C Q = Λ ∣ Ω ∣ − μ ∫ Ω S ( t , x ) d x − σ ∫ Ω I ( t , x ) d x

and

(16) D t α 0 C Q ≤ Λ ∣ Ω ∣ − σ 0 Q ,

where σ 0 = min { σ , μ } , Q = ∫ Ω ( S ( t , x ) + I ( t , x ) ) d x . From the [19, Lemma 3, page 4], by using Laplace’s transform, we obtain

Q ≤ Λ σ 0 + ∫ Ω ( S 0 + I 0 ) d x − Λ σ 0 E α ( − σ 0 t α ) ≤ Λ σ 0 ( 1 − E α ( − σ 0 t α ) ) + Q 0 E α ( − σ 0 t α ) .

By taking into account 0 ≤ E α ( − σ 0 t α ) ≤ 1 , we obtain

(17) Q ( t , x ) ≤ Q 0 + Λ σ 0 .

Thus, we conclude that the solution ( S ( t , . ) , I ( t , . ) ) of the fractional order system (1)–(3) exists uniquely and globally in time. In addition, the solution is uniformly bounded by a positive constant for large time t (see [29, Remarque 4.2, page 5] and [20, Proposition 1 and Corollary 1, page 4]).□

4 Local asymptotic stability conditions

We start by identifying the equilibrium of the proposed system (1)– (3) and then study the stability of the system.

Definition 13

We call the point ( S ∗ , I ∗ ) a constant steady state solution of (1)–(3) under to the Neumann boundary condition (3) if

(18) F 1 ( S ∗ , I ∗ ) = 0 , F 2 ( S ∗ , I ∗ ) = 0 .

Note that by using Theorem 4 and [1], we end up with the following proposition, where E 0 = Λ μ , 0 and E ∗ = ( S ∗ , I ∗ ) are the equilibrium points with S ∗ = σ I ∗ λ ϕ ( I ∗ ) .

Proposition 14

Assuming that the incidence function ϕ ( I ) satisfies (4) and (5), system (1)–(3) has the one and only disease-free equilibrium E 0 = Λ μ , 0 . If R 0 > 1 , the system admits two equilibrium: E 0 and a positive endemic equilibrium E ∗ = ( S ∗ , I ∗ ) , where R 0 = λ σ Λ μ ϕ ′ ( 0 ) .

The following proposition establishes sufficient conditions for the local asymptotic stability of the two equilibrium points of the system 13–14.

Proposition 15

Assuming that the incidence function ϕ ( I ) satisfies (4) and (5), the following statements hold for the fractional-order system 13–14 .

Subject to R 0 < 1 , E 0 is the only locally asymptotically stable state of the system for α ∈ ( 0 , 1 ] .
Subject to R 0 > 1 , E 0 is unstable, and E ∗ is locally asymptotically stable for α ∈ ( 0 , 1 ] .

Proof

Applying Lemma 5 to the two cases R 0 < 1 and R 0 > 1 .

Evaluating the Jacobian matrix corresponding to the ODE system at E 0 = Λ μ , 0 , we find that
J ( E 0 ) = − μ − σ R 0 0 σ ( R 0 − 1 ) .
1. If R 0 < 1 ,
  ∣ arg ( λ 1 ( J ( E 0 ) ) ) ∣ = π > α π 2 as α ∈ ( 0 , 1 ]
  (because we have 0 < α π 2 ≤ . π 2 < π ).
2. On the other hand,
  ∣ arg ( λ 2 ( J ( E 0 ) ) ) ∣ = ∣ arg ( σ ( R 0 − 1 ) ) ∣ = π > α π 2
  as α ∈ ( 0 , 1 ] . It follows that the equilibrium point E 0 is locally asymptotically stable.
3. If R 0 > 1 , the eigenvalues λ 1 = − μ and λ 2 = σ ( R 0 − 1 ) > 0 and thus
  ∣ arg ( λ 2 ) ∣ = 0 < α π 2 ,
  because α is assumed positive. Then, the equilibrium E 0 clearly becomes asymptotically unstable for R 0 > 1 and α is assumed in the ( 0 , 1 ] .
Discuss now the asymptotic stability of the positive endemic equilibrium E ∗ = ( S ∗ , I ∗ ) , where S ∗ , I ∗ > 0 . The Jacobian matrix J ( E ∗ ) evaluated at the E ∗ = ( S ∗ , I ∗ ) equilibrium is given as follows:
J ( E ∗ ) = − λ ϕ ( I ∗ ) − μ − λ S ∗ ϕ ′ ( I ∗ ) λ ϕ ( I ∗ ) λ S ∗ ϕ ′ ( I ∗ ) − σ , = − F 0 − μ − G 0 F 0 G 0 − σ .
Because E ∗ satisfies (18), we find
tr ( J ( E ∗ ) ) = − Λ S ∗ − λ S ∗ ϕ ( I ∗ ) I ∗ − ϕ ′ ( I ∗ ) = − [ λ ϕ ( I ∗ ) + μ ] − σ + λ S ∗ ϕ ′ ( I ∗ ) = − Λ S ∗ − λ S ∗ ϕ ( I ∗ ) I ∗ − ϕ ′ ( I ∗ ) < 0 .
On the other hand,
tr ( J ( E ∗ ) ) = − ( μ + σ ) + G 0 − F 0 < 0 .
The determinant of the Jacobian may be given as follows:
det ( J ( E ∗ ) ) = λ σ ϕ ( I ∗ ) + μ σ − μ λ S ∗ ϕ ′ ( I ∗ ) ,
and using (18), we obtain
det ( J ( E ∗ ) ) = λ λ S ∗ ϕ ( I ∗ ) I ∗ ϕ ( I ∗ ) + μ λ S ∗ ϕ ( I ∗ ) I ∗ − ϕ ′ ( I ∗ ) ,
and from (5), we obtain
det ( J ( E ∗ ) ) > 0 .
On the other hand,
det ( J ( E ∗ ) ) = ( − F 0 − μ ) ( G 0 − σ ) + F 0 G 0 = σ ( F 0 + μ ) − μ G 0 > 0 .
The characteristic equation of the Jacobian matrix is
p ( λ ) = λ 2 − t r ( J ( E ∗ ) ) λ + det J ( E ∗ ) = 0 ,
and its discriminant is
ℛ = [ tr J ( E ∗ ) ] 2 − 4 det J ( E ∗ ) = [ − ( μ + σ ) + G 0 − F 0 ] 2 − 4 [ σ ( F 0 + μ ) − μ G 0 ] .
We study the different cases separately.
1. First, if ℛ > 0 , then the eigenvalues λ 1 and λ 2 are real and can be rewritten as follows:
  λ 1 , 2 = tr J ( E ∗ ) ± ℛ 2 .
  Note that the negativity of the eigenvalues rests on the sign of the trace tr ( J ( E ∗ ) ) . Because tr ( J ( E ∗ ) ) < 0 , then λ 1 = tr J ( E ∗ ) − ℛ 2 < 0 and, therefore, ∣ arg ( λ 1 ( J ( E ∗ ) ) ) ∣ = π . Since both eigenvalues are real, tr ( J ( E ∗ ) ) < 0 , and det ( J ( E ∗ ) ) > 0 , it is evident that
  ∣ arg ( λ 1 ( J ( E ∗ ) ) ) ∣ = ∣ arg ( λ 2 ( J ( E ∗ ) ) ) ∣ = π > α π 2
  as α ∈ ( 0 , 1 ] . It follows that the equilibrium E ∗ = ( S ∗ , I ∗ ) is asymptotically stable. Both cases tr ( J ( E ∗ ) ) > 0 and tr ( J ( E ∗ ) ) = 0 are not possible.
2. If ℛ = [ tr J ( E ∗ ) ] 2 − 4 det J ( E ∗ ) = 0 . Since the eigenvalues reduce to
  λ 1 , 2 = tr J ( E ∗ ) 2 < 0 .
  Then, because α ∈ ( 0 , 1 ] , we have
  ∣ arg ( λ 1 ( J ( E ∗ ) ) ) ∣ = ∣ arg ( λ 2 ( J ( E ∗ ) ) ) ∣ = π > α π 2 .
  Consequently, the equilibrium E ∗ = ( S ∗ , I ∗ ) is asymptotically stable. Both cases tr ( J ( E ∗ ) ) > 0 and tr ( J ( E ∗ ) ) = 0 are not possible.
3. Finally, if the discriminant ℛ < 0 , then the eigenvalues λ 1 and λ 2 are complex and can be rewritten as follows:
  λ 1 , 2 = tr J ( E ∗ ) ± i − ℛ 2 ,
  and hence tr ( J ( E ∗ ) ) < 0 , then by means of Corollary 7, E ∗ = ( S ∗ , I ∗ ) is asymptotically stable. Both cases tr ( J ( E ∗ ) ) > 0 and tr ( J ( E ∗ ) ) = 0 are not possible. The proof is complete.□

Now, let us move on to the complete system (1)–(5). For this, we may use the eigenfunction expansion method described by Casten and Holland [31].

Proposition 16

For the fractional-order system (1)–(3):

Subject to R 0 < 1 , E 0 is the only steady state of the system, and for all α ∈ ( 0 , 1 ] , E 0 is locally asymptotically stable. Note that the conditions of asymptotic stability are the same as those obtained in the case of free diffusion established earlier in Proposition 14.
Subject to R 0 > 1 , E 0 is unstable and for E ∗ , we have the following cases:
1. If d 1 = d 2 , then E ∗ is asymptotically stable and subject to the same conditions of the diffusion-free case in Proposition 14.
2. If d 1 ≠ d 2 , ℛ = [ tr J ( E ∗ ) ] 2 − 4 det J ( E ∗ ) > 0 , then
  1. E ∗ is an asymptotically stable constant and is steady state if d 1 < d 2 and ( F 0 + G 0 ) − ( σ − μ ) < 0 or d 1 > d 2 and ( F 0 + G 0 ) − ( σ − μ ) > 0 , where F 0 = λ ϕ ( I ∗ ) and G 0 = λ S ∗ ϕ ′ ( I ∗ ) .
  2. If d 1 < d 2 and ( F 0 + G 0 ) − ( σ − μ ) > 0 or d 1 > d 2 and ( F 0 + G 0 ) − ( σ − μ ) < 0 , E ∗ is asymptotically stable when the eigenvalues ζ 1 , 2 ( λ i ) = tr J i ± i 4 det J i − ( tr J i ) 2 2 satisfy ∣ arg ( ζ 1 ( λ i ) ) ∣ > α π 2 and ∣ arg ( ζ 2 ( λ i ) ) ∣ > α π 2 for all λ 01 < λ i < λ 02 .
  3. If ( F 0 + G 0 ) − ( σ − μ ) = 0 , E ∗ is asymptotically stable.
3. If d 1 ≠ d 2 and ℛ ≤ 0 , ℛ i = [ tr J i ( E ∗ ) ] 2 − 4 det J i ( E ∗ ) has two real roots λ 01 < λ 02 , then E ∗ is an asymptotically stable constant steady state for λ i > λ 02 or λ i < λ 01 . It follows that E ∗ is asymptotically stable when the eigenvalues ζ 1 , 2 ( λ i ) satisfy ∣ arg ( ζ 1 , 2 ( λ i ) ) ∣ > α π 2 and ∣ arg ( λ 2 ( J i ( E 0 ) ) ) ∣ > α π 2 for all λ i ∈ ( λ 01 , λ 02 ) .

Proof

In the presence of diffusion, the stability of E 0 and E ∗ reduces to applying Lemma 5 to the linearizing operator ℒ = D Δ + A . Note that A is the Jacobian matrix evaluated at the equilibrium point.

The equilibrium point E 0 = Λ μ , 0 satisfies
d 1 Δ S + Λ − λ S ∗ ϕ ( I ∗ ) − μ S ∗ = 0 , d 2 Δ I + λ S ∗ ϕ ( I ∗ ) − σ I ∗ = 0 ,
with Neumann boundaries
∂ S ∂ ν = ∂ I ∂ ν = 0 on R + × ∂ Ω .
The linearizing operator may be given as follows:
ℒ ( E 0 ) = d 1 Δ − μ − λ Λ μ ϕ ′ ( 0 ) 0 d 2 Δ + λ Λ μ ϕ ′ ( 0 ) − σ ,
where ( λ i ) i denotes the indefinite sequence of positive eigenvalues for the Laplacian operator Δ over Ω , with Neumann boundary conditions defined by 0 = λ 0 < λ 1 ≤ λ 2 ≤ λ 3 ≤ … ↗ + ∞ , see [31]. The stability of E 0 depends on the eigenvalues of the matrices
J i ( E 0 ) = − d 1 λ i − μ − λ Λ μ ϕ ′ ( 0 ) 0 − d 2 λ i + λ Λ μ ϕ ′ ( 0 ) − σ
for all i ≥ 0 , which are given for all i ≥ 0 by
r i 1 = λ 1 ( J i ( E 0 ) ) = − d 1 λ i − μ , r i 2 = λ 2 ( J i ( E 0 ) ) = − d 2 λ i + λ Λ μ ϕ ′ ( 0 ) − σ
for all i ≥ 0 , r i 1 < 0 , and r i 2 < 0 for R 0 < 1 leads to
arg ( λ 1 ( J i ( E 0 ) ) ) = arg ( λ 2 ( J i ( E 0 ) ) ) = π
and
∣ arg ( λ 1 ( J i ( E 0 ) ) ) ∣ > α π 2 , ∣ arg ( λ 2 ( J i ( E 0 ) ) ) ∣ > α π 2 ,
leading to the local stability of E 0 = Λ μ , 0 .
The case R 0 > 1
For E 0 = Λ μ , 0 . In this case, it is easy to see that r i 1 < 0 for all i ≥ 0 , but there exists i = 0 such that r i 2 = λ Λ μ φ ′ ( 0 ) − σ = σ ( R 0 − 1 ) > 0 for R 0 > 1 , where
∣ arg ( λ 1 ( J i ( E 0 ) ) ) ∣ > α π 2 , ∣ arg ( λ 2 ( J ( E 0 ) ) ) ∣ = 0 < α π 2
because α is assumed positive. Then, the equilibrium E 0 clearly becomes asymptotically unstable for R 0 > 1 and α is assumed in the ( 0 , 1 ] .

We can state that E ∗ = ( S ∗ , I ∗ ) is asymptotically stable if the eigenvalues of the linearized system satisfy the conditions of Lemma 5. Let us set
ℒ ( E ∗ ) = d 1 Δ − λ ϕ ( I ∗ ) − μ − λ S ∗ ϕ ′ ( I ∗ ) λ ϕ ( I ∗ ) d 2 Δ + λ S ∗ ϕ ′ ( I ∗ ) − σ .
The stability of E ∗ reduces to examining the eigenvalues
ξ j ( λ i ) , j = 1 , 2 , i = 1 , 2 , 3 …
of the following matrices
J i ( E ∗ ) = − d 1 λ i − F 0 − μ − G 0 F 0 − d 2 λ i + G 0 − σ
for all i ≥ 0 , where
F 0 = λ ϕ ( I ∗ ) > 0 , G 0 = λ S ∗ ϕ ′ ( I ∗ ) > 0 .
The characteristic equation corresponding to matrix J i ( E ∗ ) is as follows:
(19) det ( J i ( E ∗ ) − ζ ( λ i ) I ) = [ ζ ( λ i ) ] 2 − ( tr J i ( E ∗ ) ) ζ ( λ i ) + det J i ( E ∗ ) = 0 ,
where the trace of J i ( E ∗ ) is given as follows:
tr ( J i ( E ∗ ) ) = − λ i ( d 1 + d 2 ) − ( μ + σ ) + G 0 − F 0 = − λ i ( d 1 + d 2 ) + tr ( J ( S ∗ , I ∗ ) ) < 0 .
The determinant is
det ( J i ( E ∗ ) ) = d 1 d 2 λ i 2 + λ i H 0 + det ( J ( E ∗ ) ) > 0
for all i ≥ 0 , where
H 0 = d 1 ( σ − G 0 ) + d 2 ( F 0 + μ ) ≥ 0 .

By taking the discriminant of (19), we have

ℛ i = ( d 1 − d 2 ) 2 λ i 2 − 2 ( d 1 − d 2 ) [ ( σ − G 0 ) − ( F 0 + μ ) ] λ i + ℛ .

Hence,

(20) ℛ i = ( d 1 − d 2 ) 2 λ i 2 + 2 ( d 1 − d 2 ) H 1 λ i + ℛ ,

where H 1 = ( F 0 + G 0 ) − ( σ − μ ) and

ℛ = [ tr J ( E ∗ ) ] 2 − 4 det J ( E ∗ ) ℛ = [ − ( μ + σ ) + G 0 − F 0 ] 2 − 4 [ σ ( F 0 + μ ) − μ G 0 ] ,

because

det ( J ( E ∗ ) ) = λ σ ϕ ( I ∗ ) + μ σ − μ λ S ∗ ϕ ′ ( I ∗ ) = σ ( F 0 + μ ) − μ G 0

and

tr J ( E ∗ ) = − ( μ + σ ) + G 0 − F 0 .

The sign of ℛ i is important for the stability of E ∗ , and since ℛ i is a quadratic polynomial in λ i , its sign depends on its own discriminant; we study the different cases separately. The discriminant of ℛ i with respect to λ i is

ℛ λ = ( d 1 − d 2 ) 2 [ ( F 0 + μ ) − ( σ − G 0 ) ] 2 − ( d 1 − d 2 ) 2 ℛ ℛ λ = 4 ( d 1 − d 2 ) 2 F 0 G 0 .

We have two cases for ℛ λ :

The case d 1 = d 2 , if d 1 = d 2 , it follows that
ℛ i = ℛ .
Hence, the exact same conditions for fractional ordinary differential equations (FODEs) stability as described in Proposition 14 apply here.
The case d 1 ≠ d 2 , in this case, ( d 1 − d 2 ) 2 > 0 . Then, ℛ λ = 4 ( d 1 − d 2 ) 2 F 0 G 0 > 0 . Hence, ℛ i has two real roots λ 01 < λ 02 and we have two cases as follows:
1. If d 1 < d 2 and H 1 < 0 , we have 2 ( d 1 − d 2 ) H 1 = 2 ( d 1 − d 2 ) [ ( F 0 + G 0 ) − ( σ − μ ) ] > 0 . Since ℛ > 0 , then the solutions λ 01 and λ 02 of the equation ℛ i = 0 are both negative regardless of i , ( λ i ≥ 0 ). Hence, we have ℛ i > 0 for all i , and the roots of (19) are
  ζ 1 ( λ i ) = 1 2 [ tr ( J i ( E ∗ ) ) − ( tr J i E ∗ ) 2 − 4 det ( J i ( E ∗ ) ) ] ,
  and
  ζ 2 ( λ i ) = 1 2 [ tr ( J i ( E ∗ ) ) + ( tr J i E ∗ ) 2 − 4 det ( J i ( E ∗ ) ) ] .
  Since tr J i ( E ∗ ) < 0 , det J i ( E ∗ ) > 0 , which implies that ζ 1 and ζ 2 are negative ( ξ 1 , 2 ( λ i ) ∈ R − ) that satisfy ∣ arg ( ζ 1 ( λ i ) ) ∣ = π > α π 2 and ∣ arg ( ζ 2 ( λ i ) ) ∣ > α π 2 , which guarantees the asymptotic stability of E ∗ .
2. If d 1 < d 2 and H 1 > 0 , we have 2 ( d 1 − d 2 ) H 1 = 2 ( d 1 − d 2 ) [ ( F 0 + G 0 ) − ( σ − μ ) ] < 0 . Since ℛ > 0 , it is evident that ℛ i has two strictly positive roots, which will be denoted by λ 01 < λ 02 . The sign of ℛ i depends on λ i , which gives us two distinct cases:
3. Alternatively, for λ i > λ 02 or 0 ≤ λ i < λ 01 , we have ℛ i > 0 , and the roots of (19) are
  ζ 1 ( λ i ) = 1 2 [ tr ( J i ( E ∗ ) ) − ( tr J i E ∗ ) 2 − 4 det ( J i ( E ∗ ) ) ] ,
  and
  ζ 2 ( λ i ) = 1 2 [ tr ( J i ( E ∗ ) ) + ( tr J i E ∗ ) 2 − 4 det ( J i ( E ∗ ) ) ] ,
  and since tr J i ( E ∗ ) < 0 , det J i ( E ∗ ) > 0 , which implies that ζ 1 and ζ 2 are
  ξ 1 , 2 ( λ i ) ∈ R − ,
  and thus the arguments of the eigenvalues are ∣ arg ( ξ 1 , 2 ) ∣ = π .
4. Alternatively, if λ 01 < λ i < λ 02 , then ℛ i < 0 , and consequently, we have two complex eigenvalues given by
  ξ 1 , 2 ( λ i ) = tr J i ± i 4 det J i − ( tr J i ) 2 2 ,
  and the equilibrium E ∗ = ( S ∗ , I ∗ ) is asymptotically stable if the eigenvalues ζ 1 , 2 ( λ i ) satisfy (11).
5. If d 1 > d 2 and H 1 < 0 , we have 2 ( d 1 − d 2 ) H 1 = 2 ( d 1 − d 2 ) [ ( F 0 + G 0 ) − ( σ − μ ) ] < 0 . Since ℛ > 0 , it is evident that ℛ i has two strictly positive roots, which will be denoted by λ 01 < λ 02 . The sign of ℛ i depends on λ i , which gives us two distinct cases:
6. Alternatively, for λ i > λ 02 or 0 ≤ λ i < λ 01 , we have ℛ i > 0 , the roots of (19) are
  ζ 1 ( λ i ) = 1 2 [ tr ( J i ( E ∗ ) ) − ( tr J i E ∗ ) 2 − 4 det ( J i ( E ∗ ) ) ] ,
  and
  ζ 2 ( λ i ) = 1 2 [ tr ( J i ( E ∗ ) ) + ( tr J i E ∗ ) 2 − 4 det ( J i ( E ∗ ) ) ] ,
  and since tr J i ( E ∗ ) < 0 , det J i ( E ∗ ) > 0 , which implies that ζ 1 and ζ 2 are
  ξ 1 , 2 ( λ i ) ∈ R − ,
  and thus the arguments of the eigenvalues are
  ∣ arg ( ξ 1 , 2 ) ∣ = π .
7. Alternatively, if λ 01 < λ i < λ 02 , then ℛ i < 0 , and consequently, we have two complex eigenvalues given by
  ξ 1 , 2 ( λ i ) = tr J i ± i 4 det J i − ( tr J i ) 2 2 ,
  and the equilibrium E ∗ = ( S ∗ , I ∗ ) is asymptotically stable if the eigenvalues ζ 1 , 2 ( λ i ) satisfy (11).
8. If d 1 > d 2 and H 1 > 0 , we have 2 ( d 1 − d 2 ) H 1 = 2 ( d 1 − d 2 ) [ ( F 0 + G 0 ) − ( σ − μ ) ] > 0 . Since ℛ > 0 , the solutions λ 01 and λ 02 of the equation ℛ i = 0 are both negative regardless of i , ( λ i ≥ 0 ). Hence, we have ℛ i > 0 for all i , and the roots of (19) are
  ζ 1 ( λ i ) = 1 2 [ tr ( J i ( E ∗ ) ) − ( tr J i E ∗ ) 2 − 4 det ( J i ( E ∗ ) ) ] ,
  and
  ζ 2 ( λ i ) = 1 2 [ tr ( J i ( E ∗ ) ) + ( tr J i E ∗ ) 2 − 4 det ( J i ( E ∗ ) ) ] .
  Since tr J i ( E ∗ ) < 0 , det J i ( E ∗ ) > 0 , which implies that ζ 1 and ζ 2 are negative ( ξ 1 , 2 ( λ i ) ∈ R − ) that satisfy ∣ arg ( ζ 1 ( λ i ) ) ∣ = π > α π 2 and ∣ arg ( ζ 2 ( λ i ) ) ∣ > α π 2 , which guarantees the asymptotic stability of E ∗ .
9. If ( F 0 + G 0 ) − ( σ − μ ) = 0 , we have
  ℛ i = ( d 1 − d 2 ) 2 λ i 2 + ℛ > 0 ,
  for all i , and the roots of (19) are
  ζ 1 , 2 ( λ i ) = 1 2 [ tr ( J i ( E ∗ ) ) − ( tr J i E ∗ ) 2 − 4 det ( J i ( E ∗ ) ) ] ,
  and since tr J i ( E ∗ ) < 0 , det J i ( E ∗ ) > 0 , which implies that ζ 1 and ζ 2 are
  ξ 1 , 2 ( λ i ) ∈ R − ,
  and thus the arguments of the eigenvalues are
  □ ∣ arg ( ξ 1 , 2 ) ∣ = π .

5 Global asymptotic stability

Next, we study the global asymptotic stability of the two steady states E 0 and E ∗ . The global stability depends on the reproduction number R 0 , which is why we have decided to treat the scenarios R 0 < 1 and R 0 > 1 separately.

5.1 Global stability of E 0 when R 0 < 1

In the beginning, we state a lemma that was developed in [32], which will be useful later.

Lemma 17

Condition (5) implies that

(21) 0 < ϕ ( I ) I ≤ ϕ ′ ( 0 ) for a l l I > 0 .

Let us suggest the following candidate Lyapunov function (see [14,21]):

(22) ℒ ˜ ( t ) = ∫ Ω S − u ˜ − u ˜ ln S u ˜ + I ( t , x ) d x .

Theorem 18

For α ∈ ( 0 , 1 ] and R 0 < 1 , E 0 is globally asymptotically stable in D . Alternatively, for R 0 > 1 , E 0 is unstable.

Proof

Taking the time fractional order derivative in Caputo sense of (22) along the solution of the system, we obtain

D t α 0 C ℒ ˜ ( t ) = ∫ Ω D t α 0 C S − u ˜ − u ˜ ln S u ˜ + D t α 0 C ( I ) d x ,

where u ˜ = Λ μ . Using Lemma 10, we obtain

D t α 0 C ℒ ˜ ( t ) ≤ ∫ Ω 1 − u ˜ S D t α 0 C S d x + ∫ Ω D t α 0 C I ( t , x ) d x

and substituting D t α 0 C S , D t α 0 C I with its values, we obtain

D t α 0 C ℒ ˜ ( t ) ≤ d 1 ∫ Ω 1 − u ˜ S Δ S d x + ∫ Ω 1 − u ˜ S ( Λ − λ S ϕ ( I ) − μ S ) d x + d 2 ∫ Ω Δ I d x + ∫ Ω ( λ S ϕ ( I ) − σ I ) d x .

We apply Green’s formula with Neumann boundaries to expand the derivative to

D t α 0 C ℒ ˜ ( t ) ≤ H + H ˜ ,

where

H = − d 1 ∫ Ω u ˜ S 2 ∣ ∇ S ∣ 2 d x ≤ 0 ,

and

H ˜ = ∫ Ω 1 − u ˜ S ( Λ − μ S ) d x − λ ∫ Ω 1 − u ˜ S S ϕ ( I ) d x + ∫ Ω ( λ S ϕ ( I ) − σ I ) d x .

We can write

D t α 0 C ℒ ˜ ( t ) ≤ − μ ∫ Ω 1 − u ˜ S S − Λ μ d x − λ ∫ Ω ( S − u ˜ ) ϕ ( I ) d x + ∫ Ω λ S ϕ ( I ) d x − σ ∫ Ω I d x ≤ − μ ∫ Ω S − u ˜ S S − Λ μ d x − λ ∫ Ω ( S − u ˜ ) ϕ ( I ) d x + ∫ Ω λ S ϕ ( I ) d x − σ ∫ Ω I d x .

By defining u ˜ = Λ μ , we obtain

D t α 0 C ℒ ˜ ( t ) ≤ − μ ∫ Ω ( S − u ˜ ) 2 S − λ ∫ Ω ( S − u ˜ ) ϕ ( I ) d x + ∫ Ω λ S ϕ ( I ) d x − σ ∫ Ω I d x ≤ − μ ∫ Ω ( S − u ˜ ) 2 S + λ ∫ Ω u ˜ ϕ ( I ) d x − σ ∫ Ω I d x .

Then, using Lemma 17, we obtain

D t α 0 C ℒ ˜ ( t ) ≤ σ ( R 0 − 1 ) ∫ Ω I d x .

Clearly, if R 0 < 1 , then D t α 0 C ℒ ˜ ( t ) is negative. The proof is completed.□

5.2 Global asymptotic stability of E ∗ when R 0 > 1

First, let us state a necessary lemmas taken from [33], which will aid in what comes.

Lemma 19

[1] Assuming that ϕ satisfies criterion (5) and

(23) ℒ ( y ) = y − 1 − ln ( y ) , for all y > 0 ,

the inequality

(24) ℒ ϕ ( I ) ϕ ( I ∗ ) ≤ ℒ I I ∗

holds.

Lemma 20

[1] We have the following equivalents

1 − S ∗ S 1 − S S ∗ = − ℒ S S ∗ − ℒ S ∗ S , and

S ϕ ( I ) S ∗ ϕ ( I ∗ ) − I I ∗ 1 − I ∗ I + 1 − S ϕ ( I ) S ∗ ϕ ( I ∗ ) 1 − S ∗ S = − ℒ S ∗ S + ℒ ϕ ( I ) φ ( I ∗ ) − ℒ I I ∗ − ℒ S ϕ ( I ) I ∗ S ∗ ϕ ( I ∗ ) I .

Now we consider the candidate Lyapunov function

V ( t ) = ∫ Ω S ∗ ℒ S S ∗ + I ∗ ℒ I I ∗ d x = ∫ Ω S − S ∗ − S ∗ ln S S ∗ d x + ∫ Ω I − I ∗ − I ∗ ln I I ∗ d x .

Theorem 21

For α ∈ ( 0 , 1 ] and R 0 > 1 , E ∗ is globally asymptotically stable.

Proof

Taking the time fractional order derivative of (22), we obtain

D t α 0 C V ( t ) = ∫ Ω D t α 0 C S − S ∗ − S ∗ ln S S ∗ d x + ∫ Ω D t α 0 C I − I ∗ − I ∗ ln I I ∗ d x .

Using Lemma 10 and substituting D t α 0 C S , D t α 0 C I with its values, we obtain

D t α 0 C V ( t ) ≤ ∫ Ω 1 − S ∗ S D t α 0 C S d x + ∫ Ω 1 − I ∗ I D t α 0 C I d x ≤ d 1 ∫ Ω 1 − S ∗ S Δ S d x + d 2 ∫ Ω 1 − I ∗ I Δ I d x + ∫ Ω 1 − S ∗ S [ Λ − λ S ϕ ( I ) − μ S ] d x + ∫ Ω 1 − I ∗ I [ λ S ϕ ( I ) − σ I ] d x .

Similar to the previous scenario, we apply Green’s formula with Neumann boundaries to expand the derivative to

D t α 0 C V ( t ) ≤ H 1 + H 2 ,

where

H 1 = − d 1 ∫ Ω S ∗ S 2 ∣ ∇ S ∣ 2 d x − d 2 ∫ Ω I ∗ I 2 ∣ ∇ I ∣ 2 d x ≤ 0

and

H 2 = ∫ Ω 1 − S ∗ S [ Λ − λ S ϕ ( I ) − μ S ] d x + ∫ Ω 1 − I ∗ I [ λ S ϕ ( I ) − σ I ] d x .

Since E ∗ is the positive endemic equilibrium, it fulfills the following equations

λ S ∗ ϕ ( I ∗ ) + μ S ∗ = Λ , λ S ∗ ϕ ( I ∗ ) I ∗ = σ ,

and using the last relationship, we can rewrite H 2 as follows:

H 2 = ∫ Ω μ S ∗ 1 − S S ∗ 1 − S ∗ S d x + λ S ∗ ϕ ( I ∗ ) ∫ Ω S ϕ ( I ) S ∗ ϕ ( I ∗ ) − I ∗ I 1 − I ∗ I d x + λ S ∗ ϕ ( I ∗ ) ∫ Ω 1 − S ϕ ( I ) S ∗ ϕ ( I ∗ ) 1 − S ∗ S d x .

Applying Lemma 20 yields

H 2 = ∫ Ω μ S ∗ W 1 d x + λ ∫ Ω S ∗ ϕ ( I ∗ ) W 2 d x ,

where

W 1 = − ℒ S S ∗ − ℒ S ∗ S

and

W 2 = − ℒ S ∗ S + ℒ ϕ ( I ) ϕ ( I ∗ ) − ℒ I I ∗ − ℒ S ϕ ( I ) I ∗ S ∗ ϕ ( I ∗ ) I ,

which implies

H 2 = μ S ∗ ∫ Ω − ℒ S S ∗ − ℒ S ∗ S d x + λ S ∗ ϕ ( I ∗ ) ∫ Ω ℒ ϕ ( I ) ϕ ( I ∗ ) − ℒ S ∗ S d x − λ S ∗ ϕ ( I ∗ ) ∫ Ω ℒ I I ∗ + ℒ S ϕ ( I ) I ∗ S ∗ ϕ ( I ∗ ) I .

Then, by applying Lemma 19, we obtain D t α 0 C V ( t ) ≤ 0 . With thus, the proof is completed.□

6 Numerical simulations

The aim of this section is to present some numerical simulations to illustrate our theoretical results and also to investigate the nature of the solutions of the epidemiological fractional reaction-diffusion model (1)–(3) for distinct values of the fractional order α ∈ ( 0 , 1 ] .

6.1 Example 1

In the first example, we introduce the function

ϕ ( I ) = β 1 − β 2 I m + I I for all β 1 ≥ β 2 > 0 .

This results in the following problem:

(25) D t α 0 C S − d 1 Δ S = − β 1 − β 2 I m + I I S + Λ − μ S in ( 0 , ∞ ) × Ω , D t α 0 C I − d 2 Δ I = β 1 − β 2 I m + I I S − σ I in ( 0 , ∞ ) × Ω , S ( 0 , x ) = S 0 ( x ) , I ( 0 , x ) = I 0 ( x ) on Ω , ∂ S ∂ ν = ∂ I ∂ ν = 0 , on ( 0 , ∞ ) × ∂ Ω ,

which is a special case of the studied system (1)–(3). The system (25) is identical to the system proposed in [9] but with d 1 = d 2 = 0 , where the ODE scenario of this model was studied. The conditions (4) and (5) are clearly satisfied as follows:

ϕ ( 0 ) = 0 , ϕ ′ ( I ) = β 1 − β 2 I m + I + I − β 2 ( m + I ) + β 2 I ( m + I ) 2 = β 1 − 2 β 2 I m + I + β 2 I 2 ( m + I ) 2 = β 2 I ( m + I ) − 1 2 − 1 + β 1 β 2 = β 2 I ( m + I ) − 1 2 − β 2 + β 1 > 0 , ϕ ′ ( 0 ) = β 1 , β 1 − β 2 I m + I I − I 2 β 2 − β 2 I ( m + I ) = I ϕ ′ ( I ) ≤ ϕ ( I ) .

This is because

I ϕ ′ ( I ) = ϕ ( I ) − β 2 I 2 m + I + β 2 I 3 ( m + I ) 2 = ϕ ( I ) − β 2 m I 3 ( m + I ) 2 ≤ ϕ ( I ) .

System (25) possesses two constant steady states

E 0 = Λ μ , 0 and E ∗ = σ β 1 , μ σ ( R 0 − 1 ) .

Note that, for α ∈ ( 0 , 1 ] , the second steady state E ∗ exists only when the reproduction number R 0 = Λ μ σ ϕ ′ ( 0 ) = Λ μ σ β 1 > 1 and is globally asymptotically stable.

In addition, note that the first steady state E 0 is globally asymptotically stable if R 0 < 1 (Table 1).

Table 1

Simulation parameters for the first example: system (25) with different values for α = 1 and α = 0.9

	Set	Figure	m	S 0	I 0	d 1	d 2	B 1	B 2	σ	Λ	μ	R 0
FODEs	Set 1	Figure 1	20	320	20	—	—	0.0007	0.0005	0.25	16	0.05	0.896
Case	Set 2	Figure 2	20	320	20	—	—	0.002	0.001	0.25	16	0.05	2.56
FPDEs	Set 1	Figure 3	20	200 + cos ( x ) 300	250 + sin ( x ) 300	3	5 4	0.0007	0.0005	0.25	16	0.05	0.896
Case	Set 2	Figure 4	20	200 + cos ( x ) 300	250 + sin ( x ) 300	3	5 4	0.002	0.001	0.25	16	0.05	2.56
FODEs	Set 1	Figure 5	2	0.2 + cos ( x ) 7	0.5 + sin ( x ) 7	—	—	2 5	1 5	1	3 7	2 7	0.6
Case	Set 2	Figure 6	2	0.5 + cos ( x ) 8	2 + sin ( x ) 8	—	—	2 5	2 7	1 2	9 5	3 4	1.92

The following is the description of the results:

Figure 1 shows the solutions in the ODE case subject to set 1, with R 0 = 0.896 . R 0 ≤ 1 and E 0 is globally asymptotically stable for α = 1 and α = 0.9 .
Figure 2 depicts the solution in the ODE case subject to parameter set 2, where R 0 = 2.56 > 1 , which by Theorem 18 means that E ∗ is globally asymptotically stable for α = 1 and α = 0.9 .
Figure 3 depicts the solution in the PDE case subject to parameter set 3, where R 0 = 0.896 ≤ 1 . By Theorem 18, E 0 is globally asymptotically stable for α = 1 .
Figure 4 depicts the numerical solutions obtained by using the parameter of set 4 in the PDE case with R 0 = 2.56 > 1 for α = 0.9 .
Figure 5 shows the PDE solution obtained using parameter set 5. In this case, R 0 = 0.6 ≤ 1 and by Theorem 18, E 0 is globally asymptotically stable for α = 0.9 and α = 1 .
Figure 6 shows the PDE solution obtained using parameter set 6. Since R 0 = 1.92 > 1 , E ∗ is globally asymptotically stable.

$Figure 1 Numerical solutions of system (25) (ODE case) subject to the first set of parameters for α = 1 \alpha =1 and α = 0.9 \alpha =0.9 .$

Figure 1

Numerical solutions of system (25) (ODE case) subject to the first set of parameters for α = 1 and α = 0.9 .

$Figure 2 Numerical solutions of system (25) subject to set 2 of parameters for α = 1 \alpha =1 and α = 0.9 \alpha =0.9 .$

Figure 2

Numerical solutions of system (25) subject to set 2 of parameters for α = 1 and α = 0.9 .

$Figure 3 Numerical solutions of system (25) subject to set 3 of parameters for α = 1 \alpha =1 .$

Figure 3

Numerical solutions of system (25) subject to set 3 of parameters for α = 1 .

$Figure 4 Numerical solutions of system (25) (PDE case) subject to set 4 of parameters for α = 0.9 \alpha =0.9 .$

Figure 4

Numerical solutions of system (25) (PDE case) subject to set 4 of parameters for α = 0.9 .

$Figure 5 Numerical solutions of system (25) (PDE case) subject to set 5 of parameters for α = 1 \alpha =1 and α = 0.9 \alpha =0.9 .$

Figure 5

Numerical solutions of system (25) (PDE case) subject to set 5 of parameters for α = 1 and α = 0.9 .

$Figure 6 Numerical solutions of system (25) subject to set 6 of parameters for α = 1 \alpha =1 and α = 0.9 \alpha =0.9 .$

Figure 6

Numerical solutions of system (25) subject to set 6 of parameters for α = 1 and α = 0.9 .

6.2 Example 2

The second illustrative example that we are interested in is the fractional reaction-diffusion extension of the ODE SIR model studied in [34], which is a special case of (1)–(3) with ϕ ( I ) = β I 1 + k I , β > 0 , and k ≥ 0 . The resulting system is described as follows:

(26) D t α 0 C S − d 1 Δ S = − λ β S I 1 + k I + Λ − μ S in ( 0 , ∞ ) × Ω , D t α 0 C I − d 2 Δ I = λ β S I 1 + k I − σ I in ( 0 , ∞ ) × Ω , S ( 0 , x ) = S 0 ( x ) , I ( 0 , x ) = I 0 ( x ) on Ω , ∂ S ∂ ν = ∂ I ∂ ν = 0 , on ( 0 , ∞ ) × ∂ Ω .

The imposed conditions (4) and (5) can be easily checked. It is clear that

ϕ ( 0 ) = 0 and ϕ ( I ) > 0 for all I > 0 .

The derivative of ϕ ( I ) is also defined as follows:

ϕ ′ ( I ) = β I 1 + k I ′ = β ( 1 + k I ) − k β I ( 1 + k I ) 2 = β ( 1 + k I ) 2 > 0 and ϕ ′ ( 0 ) = β .

In addition, we have

I ϕ ′ ( I ) = I β ( 1 + k I ) 2 ≤ β I 1 + k I = ϕ ( I ) .

The constant steady states of system (26) are

E 0 = Λ μ , 0 and E ∗ = σ ∗ λ ϕ ( I ∗ ) , I ∗ = σ ( 1 + k I ∗ ) λ β , I ∗ ,

where I ∗ = μ ( R 0 − 1 ) ( λ β + k μ ) .

For α ∈ ( 0 , 1 ] , and provided that the reproduction number R 0 = λ Λ μ σ ϕ ′ ( 0 ) = λ Λ μ σ β > 1 , E ∗ exists and is globally asymptotically stable. On the other hand, for α ∈ ( 0 , 1 ] and R 0 < 1 , E 0 is globally asymptotically stable with no conditions.

6.3 Example 3

The third illustrative example that we are interested in is the fractional reaction-diffusion extension of the model studied in [7] and [35] in the case where d 1 = d 2 = 0 , α = 1 , and ϕ ( I ) = k I 1 + I β , which is a special case of (1)–(3). The resulting system is defined as follows:

(27) D t α 0 C S − d 1 Δ S = − λ k I 1 + I β S + Λ − μ S in ( 0 , ∞ ) × Ω , D t α 0 C I − d 2 Δ I = λ k I 1 + I β S − σ I in ( 0 , ∞ ) × Ω , S ( 0 , x ) = S 0 ( x ) , I ( 0 , x ) = I 0 ( x ) on Ω , ∂ S ∂ ν = ∂ I ∂ ν = 0 , on ( 0 , ∞ ) × ∂ Ω ,

for β > 0 and k > 0 . The conditions imposed can be verified as follows:

ϕ ( 0 ) = 0 , and ϕ ( I ) > 0 for all I > 0 , ϕ ′ ( I ) = k 1 + I β 2 > 0 , ϕ ′ ( 0 ) = k , I ϕ ′ ( I ) = v k 1 + I β 2 ≤ k I 1 + I β = ϕ ( I ) .

The steady states of system 27 are defined by E 0 = Λ μ , 0 and E ∗ = σ I ∗ λ ϕ ( I ∗ ) , I ∗ = σ ( β + I ∗ ) λ β k , I ∗ , where I ∗ = μ β ( R 0 − 1 ) ( λ k β + μ ) with the reproduction number R 0 = λ Λ μ σ ϕ ′ ( 0 ) = λ Λ μ σ k > 1 . Note that, for α ∈ ( 0 , 1 ] , if E ∗ exists, then it is globally asymptotically stable and that E 0 is globally asymptotically stable when R 0 < 1 (Table 2).

Table 2

Parameters for the numerical simulation of the third example: system (27) with different values for α = 1 and α = 0.9

	Set	S 0	I 0	d 1	d 2	λ	β	σ	Λ	μ	k	R 0
FODEs case	Set 1	1	1.2	−	−	1	2	1 2	5	2	2	10
	Set 2	4	6	−	−	2 3	1	1 3	6	3 7	1 3	9.33
	Set 3	4	6	−	−	1	2	5	7	2 5	2	7
	Set 4	0.5	3	−	−	1 2	1	1	3 5	2 7	2 5	0.6
FPDEs case	Set 1	0.2 + cos ( x ) 10	0.4 + sin ( x ) 10	3	5 4	1	2	1 2	5	2	2	10
	Set 2	6 + cos ( x ) 7	4 + sin ( x ) 8	5	2	1	2	3	6	1 3	2	12
	Set 3	2.6 + cos ( x ) 7	2.4 + sin ( x ) 8	2	1	2 3	1	1 3	6	3 7	1 3	9.33
	Set 4	4 + cos ( x ) 10	6 + sin ( x ) 10	3	5 4	1	2	5	7	2 5	2	7
	Set 5	0.7 + cos ( x ) 7	0.5 + sin ( x ) 8	13 4	2	1 2	1	1	3 5	2 7	2 5	0.6

Figure 7, 8, 9, and 10 are graphics for the third example.

$Figure 7 Numerical solutions of system (27) subject to four sets of parameters for α = 1 \alpha =1 .$

Figure 7

Numerical solutions of system (27) subject to four sets of parameters for α = 1 .

$Figure 8 Numerical solutions of system (27) subject to four sets of parameters for α = 0.9 \alpha =0.9 .$

Figure 8

Numerical solutions of system (27) subject to four sets of parameters for α = 0.9 .

$Figure 9 Numerical solutions of system (27) subject to five sets of parameters for α = 1 \alpha =1 .$

Figure 9

Numerical solutions of system (27) subject to five sets of parameters for α = 1 .

$Figure 10 Numerical solutions of system (27) subject to five sets of parameters for α = 0.9 \alpha =0.9 .$

Figure 10

Numerical solutions of system (27) subject to five sets of parameters for α = 0.9 .

7 Conclusion

In this work, we have generalized an integer-order model proposed earlier in Djebara et al. [1] to a fractional-order model of the reaction-diffusion type. We proved the local asymptotic stability of the two equilibrium by using Lemma 5 in the two cases, R 0 < 1 and R 0 > 1 . We also demonstrated the global asymptotic stability of the equilibrium, where we used some suitable Lyapunov functions and the fractional LaSalle’s theorem on the proposed model in the ODE and PDE cases. For future work, we believe that it will be very important to study the behavior of the system at the critical state R 0 = 1 .

Acknowledgment

The authors would like to thank the reviewers and editors for their careful reading and for their valuable critical feedback and suggestions, which helped us improve the clarity and quality of this article.

Conflict of interest: The authors state that there is no conflict of interest.

References

[1] L. Djebara, R. Douaifia, S. Abdelmalek, and S. Bendoukha, Global and local asymptotic stability of an epidemic reaction-diffusion model with a nonlinear incidence, Math. Meth. Appl. Sci. 45 (2022), 1–25, https://doi.org/10.1002/mma.8205. 10.1002/mma.8205Search in Google Scholar

[2] L. Djebara, S. Abdelmalek, and S. Bendoukha, Global existence and asymptotic behavior of solutions for some coupled systems via a Lyapunov functional, J. Acta Math. Sci. 39 (2019), no. 6, 1538–1550, https://doi.org/10.1007/s10473-019-0606-7. 10.1007/s10473-019-0606-7Search in Google Scholar

[3] H. W. Hethcote, Qualitative analysis for communicable disease models, Math Biosci. 28 (1976), 335–356, https://doi.org/10.1016/0025-5564(76)90132-2. 10.1016/0025-5564(76)90132-2Search in Google Scholar

[4] A. Korobeinikov and G. C. Wake, Lyapunov functions and global stability for SIR, SIRS and SIS epidemiological models, Appl. Math. Lett. 15 (2002), no. 8, 955–961, https://doi.org/10.1007/s11538-005-9037-9. 10.1016/S0893-9659(02)00069-1Search in Google Scholar

[5] F. Brauer and P. Van den Driessche, Models for transmission of disease with immigration of infectives, Math. Biosci. 171 (2001), 143–154, https://doi.org/10.1016/s0025-5564(01)00057-8. 10.1016/S0025-5564(01)00057-8Search in Google Scholar

[6] M. E. Alexander and S. M. Moghadas, Periodicity in an epidemic model with a generalized non-linear incidence, Math. Biosci. 189 (2004), 75–96, https://doi.org/10.1016/j.mbs.2004.01.003. 10.1016/j.mbs.2004.01.003Search in Google Scholar PubMed

[7] V. Capasso and G. Serio, A generalization of the Kermack-McKendrick deterministic epidemic model, Math. Biosci. 42 (1978), 43–61, https://doi.org/10.1016/0025-5564(78)90006-8. 10.1016/0025-5564(78)90006-8Search in Google Scholar

[8] J. Li, Y. Yang, Y. Xiao, and S. Liu, A class of Lyapunov functions and the global stability of some epidemic models with nonlinear incidence, Appl. Anal. Comput. 6 (2016), no. 1, 38–46, https://doi.org/10.11948/2016004. 10.11948/2016004Search in Google Scholar

[9] K. Akdim, A. Zetouni, and M. Zahid, The influence of awareness campaigns on the spread of an infectious disease: a qualitative analysis of a fractional epidemic model, Model. Earth Syst. Environ. 8 (2022), 1311–1319, https://doi.org/10.1007/s40808-021-01158-9. 10.1007/s40808-021-01158-9Search in Google Scholar PubMed PubMed Central

[10] Z. Lu, Y. Yu, G. Ren, C. Xu, and X. Meng, Global dynamics for a class of reaction-diffusion multigroup SIR epidemic models with time fractional-order derivatives, Nonlinear Anal. Modell. Control 27 (2022), no. 1, 142–162, https://doi.org/10.15388/namc.2022.27.25192. 10.15388/namc.2022.27.25192Search in Google Scholar

[11] D. Mansouri, S. Abdelmalek, and S. Bendoukh, On the asymptotic stability of the time-fractional Lengyel-Epstein system, Comput. Math. Appl. 78 (2019), 1415–1430, https://doi.org/10.48550/arXiv.1809.10544. 10.1016/j.camwa.2019.04.015Search in Google Scholar

[12] S. Bendoukha and S. Abdelmalek, The fractional Chua chaotic system: dynamics, synchronization, and application to secure communications, Int. J. Nonlinear Sci. Numer. Simulat. 20 (2019), no. 1, 77–88, https://doi.org/10.1515/ijnsns-2018-0195. 10.1515/ijnsns-2018-0195Search in Google Scholar

[13] R. Mezhoud, K. Saoudi, A. Zarai, and S. Abdelmalek, Conditions for the local and global asymptotic stability of the time-fractional Degn-Harrison system, IJNSNS 21 (2020), no. 7–8, 749–759, https://doi.org/10.1515/ijnsns-2019-0159. 10.1515/ijnsns-2019-0159Search in Google Scholar

[14] S. Ullah, M. A. Khan, and M. Farooq, A fractional model for the dynamics of TB virus, Chaos Solutions Fractals 116 (2018), 63–71, https://doi.org/10.1016/j.chaos.2018.09.001. 10.1016/j.chaos.2018.09.001Search in Google Scholar

[15] H. A. A. El-Saka, The fractional-order SIS epidemic model with variable population size, J. Egypt. Math. Soc. 22 (2014), 50–54, https://doi.org/10.1016/j.joems.2013.06.006. 10.1016/j.joems.2013.06.006Search in Google Scholar

[16] Y. Li, Y. Q. Chenb, and I. Podlubnyc, Stability of fractional-order nonlinear dynamic systems: Lyapunov direct method and generalized Mittag Leffler stability, Comput. Math. Appl. 59 (2010), 1810–1821, https://doi.org/j.camwa.2009.08.019. 10.1016/j.camwa.2009.08.019Search in Google Scholar

[17] A. Nabti and B. Ghanbari, Global stability analysis of a fractional SVEIR epidemic model, Math. Methods Appl. Sci. 44 (July 2021), 1–21, https://doi.org/10.1002/mma.7285. 10.1002/mma.7285Search in Google Scholar

[18] F. Jarad and T. Abdeljawad, Generalized fractional derivatives and Laplace transform, Discrete Contin. Dyn. Syst. Ser. S. 13 (March 2020), no. 3, 709–722, https://doi.org/10.3934/dcdss.2020039. 10.3934/dcdss.2020039Search in Google Scholar

[19] H. L. Li, L. Zhang, C. Hu, Y. L. Jiang, and Z. Teng, Dynamical analysis of a fractional-order predator-prey model incorporating aprey refuge, J. Appl. Math. Comput. 54 (2017), 435–449, https://doi.org/10.1007/s12190-016-1017-8. 10.1007/s12190-016-1017-8Search in Google Scholar

[20] R. Douaifia, S. Abdelmalek, and S. Bendoukha, Asymptotic stability conditions for autonomous time-fractional reaction-diffusion systems, Commun. Nonlinear Sci. Numer. Simulat. 80 (2020), 104982, https://doi.org/10.1016/j.cnsns.2019.104982. 10.1016/j.cnsns.2019.104982Search in Google Scholar

[21] H. A. A. El-Saka, A. A. M. Arafa, and M. I. Gouda, Dynamical analysis of a fractional SIRS model on homogeneous networks, Adv. Difference Equations 2019 (2019), 144, https://doi.org/10.1186/s13662-019-2079-3. 10.1186/s13662-019-2079-3Search in Google Scholar

[22] H. Miao, X. Abdurahman, Z. Teng, and C. Kang, Global dynamics of a fractional order HIV model with both virus-to-cell and cell-to-cell transmissions and therapy effect, IAENG Int. J. Appl. Math. 47 (2017), no. 1, 75–81. Search in Google Scholar

[23] M. A. Duarte-Mermoud, N. Aguila-Camacho, J. A. Gallegos, and R. Castro-Linares, Using general quadratic Lyapunov functions to prove Lyapunov uniform stability for fractional order systems, Commun. Nonlinear Sci. Numer. Simulat. 22 (2015), 650–659, https://doi.org/10.1016/j.cnsns.2014.10.008. 10.1016/j.cnsns.2014.10.008Search in Google Scholar

[24] A. Boukhouima, K. Hattaf, El. M. Lotfi, M. Mahrouf, D. F. M. Torres, and N. Yousfi, Lyapunov functions for fractional-order systems in biology: Methods and applications, Chaos Solutions Fractals 140 (2020), 110224, https://doi.org/10.1016/j.chaos.2020.110224. 10.1016/j.chaos.2020.110224Search in Google Scholar

[25] N. Aguila-Camacho, M. A. Duarte-Mermoud, and J. A. Gallegos, Lyapunov functions for fractional order systems, Commun. Nonlinear Sci. Number. Simulat. 19 (2014), 2951–2957, http://dx.doi.org/10.1016/j.cnsns.2014.01.022. 10.1016/j.cnsns.2014.01.022Search in Google Scholar

[26] C. Vargas-De-León, Volterra-type Lyapunov functions for fractional-order epidemic systems. Commun. Nonlinear Sci. Numer. Simulat. 24 (2015), 75–85, https://doi.org/10.1016/j.cnsns.2014.12.013. 10.1016/j.cnsns.2014.12.013Search in Google Scholar

[27] P. van den Driessche and J. Watmough, Reproduction numbers and sub-threshold endemic, Math. Biosci. 180 (2002), 29–48, https://doi.org/10.48550/arXiv.1809.10544. 10.1016/S0025-5564(02)00108-6Search in Google Scholar PubMed

[28] W. Lin, Global existence theory and chaos control of fractional differential equations, J. Math. Anal. Appl. 332 (2007), no. 1, 709–726, https://doi.org/10.1016/j.jmaa.2006.10.040. 10.1016/j.jmaa.2006.10.040Search in Google Scholar

[29] M. R. Ammi, M. Tahiri, M. Tilioua, A. Zeb, I. Khan, and M. Andualem, Global analysis of a time fractional order spatio temporal SIR model, Sci. Rep. 12 (2022), 5751, https://doi.org/10.1038/s41598-022-08992-6. 10.1038/s41598-022-08992-6Search in Google Scholar PubMed PubMed Central

[30] M. M. El-Borai, Some probability densities and fundamental solutions of fractional evolution equations, Chaos Solitons Fract. 14 (2002), no. 3, 433–440, https://doi.org/10.1016/S0960-0779(01)00208-9. 10.1016/S0960-0779(01)00208-9Search in Google Scholar

[31] R. G. Casten and C. J. Holland, Stability properties of solutions to systems of reaction-diffusion equations, SIAM J. Appl. Math. 33 (1977), 353–364, https://doi.org/10.1137/0133023. 10.1137/0133023Search in Google Scholar

[32] S. Abdelmalek, S. Bendoukha, and B. Rebiai, and On the stability and nonexistence of turing patterns for the generalized Lengyel-Epstein model, Math. Methods Appl. Sci. 40 (2017), no. 18, 6295–6305, https://doi.org/10.1002/mma.4457. 10.1002/mma.4457Search in Google Scholar

[33] R. P. Sigdel and C. Connell McCluskey, Global stability for an SEI model of infectious disease with immigration, Appl. Math. Comp. 243 (2014), 684–689, https://doi.org/10.1016/j.amc.2014.06.020. 10.1016/j.amc.2014.06.020Search in Google Scholar

[34] A. Lahrouz, L. Omari, and D. Kiouach, Global analysis of a deterministic and stochastic nonlinear SIRS epidemic model, Nonlinear Anal. Modell. Control 16 (2011), no. 1, 59–76, https://doi.org/10.15388/NA.16.1.14115. 10.15388/NA.16.1.14115Search in Google Scholar

[35] G. Li and W. Wang, Bifurcation analysis of an epidemic model with nonlinear incidence, Appl. Math. Comput. 214 (2009), 411–423, https://doi.org/10.1016/j.amc.2009.04.012. 10.1016/j.amc.2009.04.012Search in Google Scholar

Received: 2022-06-04

Revised: 2023-02-08

Accepted: 2023-03-29

Published Online: 2023-06-15

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/dema-2022-0224

Keywords for this article

asymptotic stability; fractional calculus; basic reproduction number; fractional Lyapunov method; Lyapunov functional

Creative Commons

BY 4.0